How to apply Otsu method to a stacked images

Dear all,

I have 21 RGB microscopic images located in one folder. I want to apply the Otsu thresholding to each image in the folder (each image has different Otsu thresholding value).

How to make a "for loop" applied to the 21 images?

I tried this code but it finds the threshold value of the first image and applies it to all 21 images.

Any idea?

Meshoo

*****************************************************************

files = dir('*.tif');

for k = 1:(numel(files));

              thresh = [];

              image = imread(files(k).name);

              thresh = graythresh(image);

              newImage = im2bw(image, thresh);

              imwrite(newImage, ['NewImage\', files(k).name,]);

end

Answer 1

[Image Analyst]

Why are you applying monochrome thresholding to RGB images? You need to extract a color channel first, or use rgb2gray().

Assuming they're monochrome, I don't see any reason why this would not work. Take the semicolon off this line:

thresh = graythresh(image)

and see what values get written to the command window.

Also, don't use image as the name of your variable because you'll blow away the built-in image() function.

[Meshooo]

Thank you very much. It is working now after I removed the semicolon.

[Image Analyst]

Removing the semicolon didn't fix it. That only says that the value of thresh gets printed to the command line but does NOT change thresh at all. If it's not "fixed" then either you did something to fix it or it was never broken in the first place.

[Meshooo]

Humm..It is working very well now after removing the semicolon. My code is below.

I want to do the following steps:

Step1: read the whole tif images in the folder

Step2: get the green channel of each image

Step3: apply the Otsu thresholding to those green channels for binarization

Step4: find the unity of all these binary images in step3 together

I'm having problem how to do step4. In manual cases, if you have 3 images then the unity of these 3 images is

unityGreen = image_1 | image_2 | image_3;

Do you have any idea how do the unity for those binary images after we wrote them in the NewImage folder?

My current code is as follows:

files = dir('*.tif');  
for k = 1:(numel(files));
  imageI = imread(files(k).name);
  imageGreen = imageI(:,:,2);    
  threshGreen = [];
  threshGreen = graythresh(imageGreen)
  newImage = im2bw(imageGreen, threshGreen);
  imwrite(newImage, ['NewImage\', files(k).name,]);
end

Unity_Green = ??

Thank you in advance..

Sincerely, Sami

[Image Analyst]

Why can't you just say:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 14;

folder = 'C:\Users\Mustafa\Pictures';
files = dir(fullfile(folder, '*.tif'));
for k = 1:(numel(files));
  fullFileName = fullfile(folder, files(k).name);
  rgbImage = imread(fullFileName);
  [rows columns numberOfColorChannels] = size(rgbImage);
  subplot(2, 2, 1);
  imshow(rgbImage);
  title('RGB Image');
  % Get dimensions/size.
  if k == 1
    % Enlarge figure to full screen.
    set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
    [rows1 columns1 numberOfColorChannels1] = size(rgbImage);
  end
  % Get green channel or entire image if monochrome.
  if numberOfColorChannels > 1
    greenChannel = rgbImage(:,:,2);
  else
    greenChannel = rgbImage; % It's really monochrome
  end
  % Resize if necessary.
  if rows1~=rows || columns1 ~= columns
    greenChannel = imresize(greenChannel, [rows1 columns1]);
  end
  % Threshold via Otsu.
  greenThreshold = graythresh(greenChannel);
  binaryImage = im2bw(greenChannel, greenThreshold);
  % Accumulate into overall binary image.
  if k == 1
    Unity_Green = binaryImage;
    % Enlarge figure to full screen.
    set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
  else
    Unity_Green = Unity_Green | binaryImage;
  end
  subplot(2, 2, 2);
  imshow(binaryImage);
  title('Binary Image');
  subplot(2, 2, 3);
  imshow(Unity_Green);
  title('Unity_Green Image');
  message = sprintf('Do you want to continue');
  button = questdlg(message, 'Continue?', 'Yes', 'No', 'Yes');
  drawnow;  % Refresh screen to get rid of dialog box remnants.
  if strcmpi(button, 'No')
    break;
  end
end
% imwrite(binaryImage, 'NewImage\Unity_Green.png');

[Meshooo]

Thank you very much, but the Unity_Green should be one binary image that contains all the binarized green channels.

In my case I have 21 RGB images, so we separated the green channels and applied the Otsu method to each channel. So we got 21 binary images of those green channels let's call them BW_Green_1, BW_Green_2, BW_Green_3,......, BW_Green_21.

Now I want to show all these 21 binary images as a single binary image by finding their unity.

I can do that manually as

Unity_Green = BW_Green_1|BW_Green_2|BW_Green_3|BW_Green_4|BW_Green_5|..........BW_Green_21;

But I want to make some loop to do that because some times I will have less or more than 21 input images.

In case of only two images the imadd can be used. But we have 21 binary images that should be added together.

Thank you very much.

[Image Analyst]

Sorry - I carelessly relied to much on your code and programming just off the top of my head without actually testing it for robustness. I've corrected it above now so it should work.

[Meshooo]

Yes, it is working very well. Indeed thank you very much..

[Meshooo]

Hi again Image Analyst, I hope you still remember my old programing problem. I wanted to continue the same program but I am having a problem for how to use the for loop.

Last time you helped me to do the following 4 steps:

Step1: read the whole tif images in the folder

Step2: get the green channel of each image

Step3: apply the Otsu thresholding to those green channels for binarization

Step4: find the unity of all these binary images in step3 together

Step5: Apply a simple morphological filter to the image generated from step4 (thickening, bridging, opening, and closing)

the code for the above 5 steps can be written as follows:

%%STEP-1****************** Load Image ***************************
dname = which (uigetfile( {...
    '*.jpg;*.tif;*.gif;*.bmp;*.png', 'All image files(*.jpg,*.tif,*.gif,*.bmp,*.png)';
    '*.jpg;*.jpeg', 'JPEG files(*.jpg)';
    '*.gif', 'GIF files(*.gif)';
    '*.tif;*.tiff', 'TIFF files(*.tif)';
    '*.bmp', 'BMP files(*.bmp)';
    '*.png', 'PNG files(*.png)';
    '*.*', 'All Files (*.*)'}, 'Open an image'));
%

filelist = dir([fileparts(dname) filesep '*.tif']);
fileNames = {filelist.name}';

%%STEP2 to 5************* Loop Images **************************

for k = 1:(numel(filelist));
    imageI = imread(filelist(k).name);
    imageGreen = imageI(:,:, 2); %Find the green channel

      threshGreen = graythresh(imageGreen);

      if threshGreen < 0.05;  %to segment weak signal
          threshGreen = 0.1;
      end

      newImage = im2bw(imageGreen, threshGreen);  %Otsu Thresholding 

      if k ==1
          Unity_Green = newImage;
      else
          Unity_Green = Unity_Green | newImage; % Sum of all Green Channels
          se = strel('square',6); % Morphological Filter
          BW1_G = Unity_Green;
          BW2_G = Unity_Green; 
          BW1_G = bwmorph(BW1_G,'thicken');
          BW2_G = bwmorph(BW2_G,'bridge');
          BW3_G = imadd(BW1_G,BW2_G);
          fo_G = imopen(BW3_G, se);
          BW_G_F = imclose(fo_G, se);
          BW_G = im2bw(BW_G_F);
          [L, num] = bwlabel(BW_G, 8);  % Labeling using 8 connectivity
      end
  end

[Meshooo]

Now I want to add few more steps as follows:

Step6: find the centroid of all objects in the image we obtained from Step5 (BW_G). This can be written as:

s  = regionprops(BW_G, 'centroid');
centroids = cat(1, s.Centroid);

imtool(BW_G)
hold(imgca,'on')
plot(imgca,centroids(:,1), centroids(:,2), 'r.')
hold(imgca,'off')

Step7: multiply the binary image (BW_G) that we obtained from step 5 with each green channel of the input 21 images. Then apply a 3x3 median filter to every image. How to make a for loop to do that? I did it without for loop and it was very long as follows. G_1 is the green channel of the first image, and G_2 is the green channel of the second image and so forth.

space_1 = immultiply(BW_G, G_1);
Mean_Filter_1 = medfilt2(space_1,[3 3]); 

space_2 = immultiply(BW_G, G_2);
Mean_Filter_2 = medfilt2(space_2,[3 3]);

space_3 = immultiply(BW_G, G_3);
Mean_Filter_3 = medfilt2(space_3,[3 3]);  

space_4 = immultiply(BW_G, G_4);
Mean_Filter_4 = medfilt2(space_4,[3 3]);

space_5 = immultiply(BW_G, G_5);
Mean_Filter_5 = medfilt2(space_5,[3 3]);

space_6 = immultiply(BW_G, G_6);
Mean_Filter_6 = medfilt2(space_6,[3 3]);

space_7 = immultiply(BW_G, G_7);
Mean_Filter_7 = medfilt2(space_7,[3 3]);

space_8 = immultiply(BW_G, G_8);
Mean_Filter_8 = medfilt2(space_8,[3 3]);

space_9 = immultiply(BW_G, G_9);
Mean_Filter_9 = medfilt2(space_9,[3 3]);

space_10 = immultiply(BW_G, G_10);
Mean_Filter_10 = medfilt2(space_10,[3 3]);

space_11 = immultiply(BW_G, G_11);
Mean_Filter_11 = medfilt2(space_11,[3 3]);  

space_12 = immultiply(BW_G, G_12);
Mean_Filter_12 = medfilt2(space_12,[3 3]);

space_13 = immultiply(BW_G, G_13);
Mean_Filter_13 = medfilt2(space_13,[3 3]);  

space_14 = immultiply(BW_G, G_14);
Mean_Filter_14 = medfilt2(space_14,[3 3]);

space_15 = immultiply(BW_G, G_15);
Mean_Filter_15 = medfilt2(space_15,[3 3]);

space_16 = immultiply(BW_G, G_16);
Mean_Filter_16 = medfilt2(space_16,[3 3]);

space_17 = immultiply(BW_G, G_17);
Mean_Filter_17 = medfilt2(space_17,[3 3]);

space_18 = immultiply(BW_G, G_18);
Mean_Filter_18 = medfilt2(space_18,[3 3]);

space_19 = immultiply(BW_G, G_19);
Mean_Filter_19 = medfilt2(space_19,[3 3]);

space_20 = immultiply(BW_G, G_20);
Mean_Filter_20 = medfilt2(space_20,[3 3]);

[Meshooo]

Step8: To find the intensity value of the determined central point in each object of the green channels. In my case I have 21 images and therefore 21 green channel images. After segmentation we had 6 objects in the BW_G image that we found the centroid of each object in BW_G. Now I want to find the intensity value of this centroid in each green channel of the 20 input images. I did that without using the for loop and the code is as follows:

center = round(centroids(:,:,1)); %get the round value of all center pixles

Then applied to all the centroids of the 6 objects as follows:

%%******************** Centroid #1 *****************************

P_1 = impixel(Mean_Filter_1,(center(1,1)),(center(1,2)));   %the index intensity value of X,Y
PP_1 = P_1(:,1);

P_2 = impixel(Mean_Filter_2,(center(1,1)),(center(1,2)));   
PP_2 = P_2(:,1);

P_3 = impixel(Mean_Filter_3,(center(1,1)),(center(1,2)));   
PP_3 = P_3(:,1);

P_4 = impixel(Mean_Filter_4,(center(1,1)),(center(1,2)));   
PP_4 = P_4(:,1);

P_5 = impixel(Mean_Filter_5,(center(1,1)),(center(1,2)));   
PP_5 = P_5(:,1);

P_6 = impixel(Mean_Filter_6,(center(1,1)),(center(1,2)));   
PP_6 = P_6(:,1);

P_7 = impixel(Mean_Filter_7,(center(1,1)),(center(1,2)));   
PP_7 = P_7(:,1);

P_8 = impixel(Mean_Filter_8,(center(1,1)),(center(1,2)));   
PP_8 = P_8(:,1);

P_9 = impixel(Mean_Filter_9,(center(1,1)),(center(1,2)));   
PP_9 = P_9(:,1);

P_10 = impixel(Mean_Filter_10,(center(1,1)),(center(1,2)));   
PP_10 = P_10(:,1);

P_11 = impixel(Mean_Filter_11,(center(1,1)),(center(1,2)));   
PP_11 = P_11(:,1);

P_12 = impixel(Mean_Filter_12,(center(1,1)),(center(1,2)));   
PP_12 = P_12(:,1);

P_13 = impixel(Mean_Filter_13,(center(1,1)),(center(1,2)));   
PP_13 = P_13(:,1);

P_14 = impixel(Mean_Filter_14,(center(1,1)),(center(1,2)));   
PP_14 = P_14(:,1);

P_15 = impixel(Mean_Filter_15,(center(1,1)),(center(1,2)));   
PP_15 = P_15(:,1);

P_16 = impixel(Mean_Filter_16,(center(1,1)),(center(1,2)));   
PP_16 = P_16(:,1);

P_17 = impixel(Mean_Filter_17,(center(1,1)),(center(1,2)));   
PP_17 = P_17(:,1);

P_18 = impixel(Mean_Filter_18,(center(1,1)),(center(1,2)));   
PP_18 = P_18(:,1);

P_19 = impixel(Mean_Filter_19,(center(1,1)),(center(1,2)));   
PP_19 = P_19(:,1);

P_20 = impixel(Mean_Filter_20,(center(1,1)),(center(1,2)));   
PP_20 = P_20(:,1);

%%plot center#1

A = [PP_1, PP_2, PP_3, PP_4, PP_5, PP_6, PP_7, PP_8, PP_9, PP_10, PP_11,...
    PP_12, PP_13, PP_14, PP_15, PP_16, PP_17, PP_18, PP_19, PP_20];

%plot(A)

%%******************** Centroid #2 *****************************

P2_1 = impixel(Mean_Filter_1,(center(2,1)),(center(2,2)));   %the index intensity value of X,Y
PP2_1 = P2_1(:,1);

P2_2 = impixel(Mean_Filter_2,(center(2,1)),(center(2,2)));   
PP2_2 = P2_2(:,1);

P2_3 = impixel(Mean_Filter_3,(center(2,1)),(center(2,2)));   
PP2_3 = P2_3(:,1);

P2_4 = impixel(Mean_Filter_4,(center(2,1)),(center(2,2)));   
PP2_4 = P2_4(:,1);

P2_5 = impixel(Mean_Filter_5,(center(2,1)),(center(2,2)));   
PP2_5 = P2_5(:,1);

P2_6 = impixel(Mean_Filter_6,(center(2,1)),(center(2,2)));   
PP2_6 = P2_6(:,1);

P2_7 = impixel(Mean_Filter_7,(center(2,1)),(center(2,2)));   
PP2_7 = P2_7(:,1);

P2_8 = impixel(Mean_Filter_8,(center(2,1)),(center(2,2)));   
PP2_8 = P2_8(:,1);

P2_9 = impixel(Mean_Filter_9,(center(2,1)),(center(2,2)));   
PP2_9 = P2_9(:,1);

P2_10 = impixel(Mean_Filter_10,(center(2,1)),(center(2,2)));   
PP2_10 = P2_10(:,1);

P2_11 = impixel(Mean_Filter_11,(center(2,1)),(center(2,2)));   
PP2_11 = P2_11(:,1);

P2_12 = impixel(Mean_Filter_12,(center(2,1)),(center(2,2)));   
PP2_12 = P2_12(:,1);

P2_13 = impixel(Mean_Filter_13,(center(2,1)),(center(2,2)));   
PP2_13 = P2_13(:,1);

P2_14 = impixel(Mean_Filter_14,(center(2,1)),(center(2,2)));   
PP2_14 = P2_14(:,1);

P2_15 = impixel(Mean_Filter_15,(center(2,1)),(center(2,2)));   
PP2_15 = P2_15(:,1);

P2_16 = impixel(Mean_Filter_16,(center(2,1)),(center(2,2)));   
PP2_16 = P2_16(:,1);

P2_17 = impixel(Mean_Filter_17,(center(2,1)),(center(2,2)));   
PP2_17 = P2_17(:,1);

P2_18 = impixel(Mean_Filter_18,(center(2,1)),(center(2,2)));   
PP2_18 = P2_18(:,1);

P2_19 = impixel(Mean_Filter_19,(center(2,1)),(center(2,2)));   
PP2_19 = P2_19(:,1);

P2_20 = impixel(Mean_Filter_20,(center(2,1)),(center(2,2)));   
PP2_20 = P2_20(:,1);

%%plot center#2

A2 = [PP2_1, PP2_2, PP2_3, PP2_4, PP2_5, PP2_6, PP2_7, PP2_8, ...
    PP2_9, PP2_10, PP2_11, PP2_12, PP2_13, PP2_14, PP2_15, PP2_16, ...
    PP2_17, PP2_18, PP2_19, PP2_20];

%%******************** Centroid #3 *****************************

P3_1 = impixel(Mean_Filter_1,(center(3,1)),(center(3,2)));   %the index intensity value of X,Y
PP3_1 = P3_1(:,1);

P3_2 = impixel(Mean_Filter_2,(center(3,1)),(center(3,2)));   
PP3_2 = P3_2(:,1);

P3_3 = impixel(Mean_Filter_3,(center(3,1)),(center(3,2)));   
PP3_3 = P3_3(:,1);

P3_4 = impixel(Mean_Filter_4,(center(3,1)),(center(3,2)));   
PP3_4 = P3_4(:,1);

P3_5 = impixel(Mean_Filter_5,(center(3,1)),(center(3,2)));   
PP3_5 = P3_5(:,1);

P3_6 = impixel(Mean_Filter_6,(center(3,1)),(center(3,2)));   
PP3_6 = P3_6(:,1);

P3_7 = impixel(Mean_Filter_7,(center(3,1)),(center(3,2)));   
PP3_7 = P3_7(:,1);

P3_8 = impixel(Mean_Filter_8,(center(3,1)),(center(3,2)));   
PP3_8 = P3_8(:,1);

P3_9 = impixel(Mean_Filter_9,(center(3,1)),(center(3,2)));   
PP3_9 = P3_9(:,1);

P3_10 = impixel(Mean_Filter_10,(center(3,1)),(center(3,2)));   
PP3_10 = P3_10(:,1);

P3_11 = impixel(Mean_Filter_11,(center(3,1)),(center(3,2)));   
PP3_11 = P3_11(:,1);

P3_12 = impixel(Mean_Filter_12,(center(3,1)),(center(3,2)));   
PP3_12 = P3_12(:,1);

P3_13 = impixel(Mean_Filter_13,(center(3,1)),(center(3,2)));   
PP3_13 = P3_13(:,1);

P3_14 = impixel(Mean_Filter_14,(center(3,1)),(center(3,2)));   
PP3_14 = P3_14(:,1);

P3_15 = impixel(Mean_Filter_15,(center(3,1)),(center(3,2)));   
PP3_15 = P3_15(:,1);

P3_16 = impixel(Mean_Filter_16,(center(3,1)),(center(3,2)));   
PP3_16 = P3_16(:,1);

P3_17 = impixel(Mean_Filter_17,(center(3,1)),(center(3,2)));   
PP3_17 = P3_17(:,1);

P3_18 = impixel(Mean_Filter_18,(center(3,1)),(center(3,2)));   
PP3_18 = P3_18(:,1);

P3_19 = impixel(Mean_Filter_19,(center(3,1)),(center(3,2)));   
PP3_19 = P3_19(:,1);

P3_20 = impixel(Mean_Filter_20,(center(3,1)),(center(3,2)));   
PP3_20 = P3_20(:,1);

%%plot center#3

A3 = [PP3_1, PP3_2, PP3_3, PP3_4, PP3_5, PP3_6, PP3_7, PP3_8, ...
    PP3_9, PP3_10, PP3_11, PP3_12, PP3_13, PP3_14, PP3_15, PP3_16, ...
    PP3_17, PP3_18, PP3_19, PP3_20];

%%******************** Cetnroid #4 *****************************

P4_1 = impixel(Mean_Filter_1,(center(4,1)),(center(4,2)));   %the index intensity value of X,Y
PP4_1 = P4_1(:,1);

P4_2 = impixel(Mean_Filter_2,(center(4,1)),(center(4,2)));   
PP4_2 = P4_2(:,1);

P4_3 = impixel(Mean_Filter_3,(center(4,1)),(center(4,2)));   
PP4_3 = P4_3(:,1);

P4_4 = impixel(Mean_Filter_4,(center(4,1)),(center(4,2)));   
PP4_4 = P4_4(:,1);

P4_5 = impixel(Mean_Filter_5,(center(4,1)),(center(4,2)));   
PP4_5 = P4_5(:,1);

P4_6 = impixel(Mean_Filter_6,(center(4,1)),(center(4,2)));   
PP4_6 = P4_6(:,1);

P4_7 = impixel(Mean_Filter_7,(center(4,1)),(center(4,2)));   
PP4_7 = P4_7(:,1);

P4_8 = impixel(Mean_Filter_8,(center(4,1)),(center(4,2)));   
PP4_8 = P4_8(:,1);

P4_9 = impixel(Mean_Filter_9,(center(4,1)),(center(4,2)));   
PP4_9 = P4_9(:,1);

P4_10 = impixel(Mean_Filter_10,(center(4,1)),(center(4,2)));   
PP4_10 = P4_10(:,1);

P4_11 = impixel(Mean_Filter_11,(center(4,1)),(center(4,2)));   
PP4_11 = P4_11(:,1);

P4_12 = impixel(Mean_Filter_12,(center(4,1)),(center(4,2)));   
PP4_12 = P4_12(:,1);

P4_13 = impixel(Mean_Filter_13,(center(4,1)),(center(4,2)));   
PP4_13 = P4_13(:,1);

P4_14 = impixel(Mean_Filter_14,(center(4,1)),(center(4,2)));   
PP4_14 = P4_14(:,1);

P4_15 = impixel(Mean_Filter_15,(center(4,1)),(center(4,2)));   
PP4_15 = P4_15(:,1);

P4_16 = impixel(Mean_Filter_16,(center(4,1)),(center(4,2)));   
PP4_16 = P4_16(:,1);

P4_17 = impixel(Mean_Filter_17,(center(4,1)),(center(4,2)));   
PP4_17 = P4_17(:,1);

P4_18 = impixel(Mean_Filter_18,(center(4,1)),(center(4,2)));   
PP4_18 = P4_18(:,1);

P4_19 = impixel(Mean_Filter_19,(center(4,1)),(center(4,2)));   
PP4_19 = P4_19(:,1);

P4_20 = impixel(Mean_Filter_20,(center(4,1)),(center(4,2)));   
PP4_20 = P4_20(:,1);

%%plot center#4

A4 = [PP4_1, PP4_2, PP4_3, PP4_4, PP4_5, PP4_6, PP4_7, PP4_8, ...
    PP4_9, PP4_10, PP4_11, PP4_12, PP4_13, PP4_14, PP4_15, PP4_16, ...
    PP4_17, PP4_18, PP4_19, PP4_20];

%%******************** Centroid #5 *****************************

P5_1 = impixel(Mean_Filter_1,(center(5,1)),(center(5,2)));   %the index intensity value of X,Y
PP5_1 = P5_1(:,1);

P5_2 = impixel(Mean_Filter_2,(center(5,1)),(center(5,2)));   
PP5_2 = P5_2(:,1);

P5_3 = impixel(Mean_Filter_3,(center(5,1)),(center(5,2)));   
PP5_3 = P5_3(:,1);

P5_4 = impixel(Mean_Filter_4,(center(5,1)),(center(5,2)));   
PP5_4 = P5_4(:,1);

P5_5 = impixel(Mean_Filter_5,(center(5,1)),(center(5,2)));   
PP5_5 = P5_5(:,1);

P5_6 = impixel(Mean_Filter_6,(center(5,1)),(center(5,2)));   
PP5_6 = P5_6(:,1);

P5_7 = impixel(Mean_Filter_7,(center(5,1)),(center(5,2)));   
PP5_7 = P5_7(:,1);

P5_8 = impixel(Mean_Filter_8,(center(5,1)),(center(5,2)));   
PP5_8 = P5_8(:,1);

P5_9 = impixel(Mean_Filter_9,(center(5,1)),(center(5,2)));   
PP5_9 = P5_9(:,1);

P5_10 = impixel(Mean_Filter_10,(center(5,1)),(center(5,2)));   
PP5_10 = P5_10(:,1);

P5_11 = impixel(Mean_Filter_11,(center(5,1)),(center(5,2)));   
PP5_11 = P5_11(:,1);

P5_12 = impixel(Mean_Filter_12,(center(5,1)),(center(5,2)));   
PP5_12 = P5_12(:,1);

P5_13 = impixel(Mean_Filter_13,(center(5,1)),(center(5,2)));   
PP5_13 = P5_13(:,1);

P5_14 = impixel(Mean_Filter_14,(center(5,1)),(center(5,2)));   
PP5_14 = P5_14(:,1);

P5_15 = impixel(Mean_Filter_15,(center(5,1)),(center(5,2)));   
PP5_15 = P5_15(:,1);

P5_16 = impixel(Mean_Filter_16,(center(5,1)),(center(5,2)));   
PP5_16 = P5_16(:,1);

P5_17 = impixel(Mean_Filter_17,(center(5,1)),(center(5,2)));   
PP5_17 = P5_17(:,1);

P5_18 = impixel(Mean_Filter_18,(center(5,1)),(center(5,2)));   
PP5_18 = P5_18(:,1);

P5_19 = impixel(Mean_Filter_19,(center(5,1)),(center(5,2)));   
PP5_19 = P5_19(:,1);

P5_20 = impixel(Mean_Filter_20,(center(5,1)),(center(5,2)));   
PP5_20 = P5_20(:,1);

%%plot center#5

A5 = [PP5_1, PP5_2, PP5_3, PP5_4, PP5_5, PP5_6, PP5_7, PP5_8, ...
    PP5_9, PP5_10, PP5_11, PP5_12, PP5_13, PP5_14, PP5_15, PP5_16, ...
    PP5_17, PP5_18, PP5_19, PP5_20];

%%******************** Centroid #6 *****************************

P6_1 = impixel(Mean_Filter_1,(center(6,1)),(center(6,2)));   %the index intensity value of X,Y
PP6_1 = P6_1(:,1);

P6_2 = impixel(Mean_Filter_2,(center(6,1)),(center(6,2)));   
PP6_2 = P6_2(:,1);

P6_3 = impixel(Mean_Filter_3,(center(6,1)),(center(6,2)));   
PP6_3 = P6_3(:,1);

P6_4 = impixel(Mean_Filter_4,(center(6,1)),(center(6,2)));   
PP6_4 = P6_4(:,1);

P6_5 = impixel(Mean_Filter_5,(center(6,1)),(center(6,2)));   
PP6_5 = P6_5(:,1);

P6_6 = impixel(Mean_Filter_6,(center(6,1)),(center(6,2)));   
PP6_6 = P6_6(:,1);

P6_7 = impixel(Mean_Filter_7,(center(6,1)),(center(6,2)));   
PP6_7 = P6_7(:,1);

P6_8 = impixel(Mean_Filter_8,(center(6,1)),(center(6,2)));   
PP6_8 = P6_8(:,1);

P6_9 = impixel(Mean_Filter_9,(center(6,1)),(center(6,2)));   
PP6_9 = P6_9(:,1);

P6_10 = impixel(Mean_Filter_10,(center(6,1)),(center(6,2)));   
PP6_10 = P6_10(:,1);

P6_11 = impixel(Mean_Filter_11,(center(6,1)),(center(6,2)));   
PP6_11 = P6_11(:,1);

P6_12 = impixel(Mean_Filter_12,(center(6,1)),(center(6,2)));   
PP6_12 = P6_12(:,1);

P6_13 = impixel(Mean_Filter_13,(center(6,1)),(center(6,2)));   
PP6_13 = P6_13(:,1);

P6_14 = impixel(Mean_Filter_14,(center(6,1)),(center(6,2)));   
PP6_14 = P6_14(:,1);

P6_15 = impixel(Mean_Filter_15,(center(6,1)),(center(6,2)));   
PP6_15 = P6_15(:,1);

P6_16 = impixel(Mean_Filter_16,(center(6,1)),(center(6,2)));   
PP6_16 = P6_16(:,1);

P6_17 = impixel(Mean_Filter_17,(center(6,1)),(center(6,2)));   
PP6_17 = P6_17(:,1);

P6_18 = impixel(Mean_Filter_18,(center(6,1)),(center(6,2)));   
PP6_18 = P6_18(:,1);

P6_19 = impixel(Mean_Filter_19,(center(6,1)),(center(6,2)));   
PP6_19 = P6_19(:,1);

P6_20 = impixel(Mean_Filter_20,(center(6,1)),(center(6,2)));   
PP6_20 = P6_20(:,1);

%%plot center#6

A6 = [PP6_1, PP6_2, PP6_3, PP6_4, PP6_5, PP6_6, PP6_7, PP6_8, ...
    PP6_9, PP6_10, PP6_11, PP6_12, PP6_13, PP6_14, PP6_15, PP6_16, ...
    PP6_17, PP6_18, PP6_19, PP6_20];

[Meshooo]

Step9: Finally to plot the whole signals we obtained from A1 to A6 and I think this was the easiest part.

hold on
plot(A,'b-')
plot(A2,'c-')
plot(A3,'g-')
plot(A4,'y-')
plot(A5,'r-')
plot(A6,'m-')
hold off

Thank you very much for your kind old help. I hope you will have a time to check the code above.

[Meshooo]

I am afraid to post here the whole program without using for loop! I already wrote it but it is too long and scary..

[Image Analyst]

This would be difficult for me to reproduce. Why don't you just look at what's changing at each steps and make that an index k, and then put it in a loop over k? I can assure you your program will be about 21 times shorter.

[Meshooo]

Hi again and thank you for your reply. I am trying by myself and things are going well. However, I would like to ask you how to make two For-loops inside each other?

For example I want to load the whole 21 intensity images (Intensity) from the folder and find their intensity value in 6 points which I already know them.

I wrote the following program which works very well but I want to use a for loop because the number of points is not always 6.

Could you please check this code to make a for loop

   for k = 1:(numel(filelist)); 
        Intensity = imread(filelist(k).name);

        P_Centro_1 = []; %1st centroid
        P_Centro_1 = impixel(Intensity,(center(1,1)),(center(1,2)));
        P_Centro_1 = P_Centro_1(:,1); %Intensity value
        A1(:,k)= P_Centro_1;

        P_Centro_2 = []; %2nd centroid
        P_Centro_2 = impixel(Intensity,(center(2,1)),(center(2,2)));
        P_Centro_2 = P_Centro_2(:,1);
        A2(:,k)= P_Centro_2;

        P_Centro_3 = []; %3rd centroid
        P_Centro_3 = impixel(Intensity,(center(3,1)),(center(3,2)));
        P_Centro_3 = P_Centro_3(:,1);
        A3(:,k)= P_Centro_3;

        P_Centro_4 = []; %4th centroid
        P_Centro_4 = impixel(Intensity,(center(4,1)),(center(4,2)));
        P_Centro_4 = P_Centro_4(:,1);
        A4(:,k)= P_Centro_4;

        P_Centro_5 = []; %5th centroid
        P_Centro_5 = impixel(Intensity,(center(5,1)),(center(5,2)));
        P_Centro_5 = P_Centro_5(:,1);
        A5(:,k)= P_Centro_5;

        P_Centro_6 = []; %6th centroid
        P_Centro_6 = impixel(Intensity,(center(6,1)),(center(6,2)));
        P_Centro_6 = P_Centro_6(:,1);
        A6(:,k)= P_Centro_6;

    end

    hold on
    plot(A1,'b-')
    plot(A2,'c-')
    plot(A3,'g-')
    plot(A4,'y-')
    plot(A5,'r-')
    plot(A6,'m-')
    hold off

Thank you very much for your time.

[Meshooo]

The following program seems to be the answer, but how to store the results from every loop in A1?

for k = 1:21
    Intensity = imread(filelist(k).name);
    for c = 1:6
        P_Centro_1 = []; % list of centroids
        P_Centro_1 = impixel(Intensity,(center(c,1)),(center(c,2)));
        P_Centro_1 = P_Centro_1(:,1); %Intensity value
        A1(:,k)= P_Centro_1;
    end
end