How do I plot the intensity profiles of a sequence of images?

I have a sequence of over 200 images of a light beam. I need to plot the 3-D intensity profiles of each image as a sequence so I see how the profile changes. From these profiles, I need to extract the maximum intensity coordinates and the full width half maxima in both the x and y axes. In the end, I'd like to have a table that lists all these values for each image. My code currently works properly for a single iteration so I'm assuming that the problem lies in the way my loop is setup. My code is down below. I've also attached an image of what comes up when I run my program. Any help would be greatly appreciated.

 %code 
fileFolder = fullfile(matlabroot, 'toolbox', 'images', 'imdata');
dirOutput = dir('/Users/Name/Documents/MATLAB/image_*.tiff');
fileNames = {dirOutput.name}'
numFrames = numel(fileNames)
I = imread(fileNames{1});
sequence =  zeros([size(I) numFrames], class(I));
sequence(:,:,1) = I;
% 
for p = 2:numFrames
    sequence(:,:,p) = imread(fileNames{p});  
end
sequenceNew = stdfilt(sequence, ones(3));
test = 1;
while test == 1
    test = 0;
    for d = 0001:numFrames
        %intensity plot info
        Int = imread(fileNames{d});
        [x,y] = size(Int);
        X = 1:x;
        Y = 1:y;
        [xx,yy] = meshgrid(Y,X);
        i = im2double(Int);

        %max coordinates and FWHM x and y from plot
        [mz,k] = max(i(:));
        [ix,jy] = ind2sub(size(i),k);
        [ix,jy,mz];

        %max intensity coordinates in microns
        %pixel size is 6.45x6.45 microns
        ix = 6.45*ix;
        jy = 6.45*jy;
        [ix,jy,mz];
        mhz = 0.5*mz;
        hk = 0.5*k;
        [ihx,jhy]=ind2sub(size(i),hk);
        [ihx,jhy,mhz];
        FWHMx = abs(ix-ihx)*2;
        FWHMy = abs(jy-jhy)*2;

        %table of data from plot
        f = figure;
        data = [ix,jy,mz,FWHMx,FWHMy];
        colnames = {'X(Zmax) [um]', 'Y(Zmax) [um]', ...
            'Zmax/Max Intensity', 'FWHMx [um]', 'FWHMy [um]'};
        t = uitable(f, 'Data', data, 'ColumnName', colnames, 'ColumnWidth', {120});
        t.Position(3) = t.Extent(3);
        t.Position(4) = t.Extent(4);

        %make subplot of table and intensity plot
        subplot(2,1,1) = mesh(xx,yy,i); colorbar;
        title('Intensity of Beam Image as a Function of Pixel Position');
        xlabel('X(pixels)');
        ylabel('Y (pixels)');
        zlabel('Intensity (double precision corrected) [a.u.]');

        pos = get(subplot(2,1,2, 'position'));
        delete(subplot(2,1,2))
        set(t,'units','normalized')
        set(t, 'position', pos)
    end
    test = 1
end

[Walter Roberson]

I do not understand why in your loop you have

Int = imread(fileNames{d});

You already read the information in: it is in sequence(:,:,d)

[Jessica Flores]

It is a bit redundant and not needed. Thanks for pointing that out. Even with that fixed, when I run the program it still only does one iteration. How can I fix this?

[Image Analyst]

What is the value of numFrames? Is it 1? It sounds like it.

[Jessica Flores]

When I run the program, it tells me the value of numFrames and it's 200.

Answer 1

[Walter Roberson]

Your line

hk = 0.5*k;

looks wrong to me. You are taking half of the linear index, which could land you in either half of the array "closer" to the origin. You can be sure that the column position will be halved (or half minus 1) but the row position could end up on either side. There is no reason I can think of to expect that the resulting position will have any relationship to the Full Width Half Maxima. You should be searching for a location whose value is closest to mhz (the half maximum) and then the distance of that point to the maxima gives you the full width at half maxima.

[Walter Roberson]

fileFolder = '/Users/Name/Documents/MATLAB';
tiffpattern = fullfile(fileFolder, 'image_*.tiff');
dirOutput = dir(tiffpattern);
fileNames = cellfun(@(s) fullfile(fileFolder, s), {dirOutput.name}, 'Uniform', 0) .';

Remember, dir() returns only the basic name in the .name element, and you need to prefix it with the path.

[Jessica Flores]

When I did this, the while loop continues to run indefinitely and now the figure doesn't come up. Also, numFrames is zero after I prefixed it with the path.

[Walter Roberson]

Please show your revised code.

I do not understand why you have coded a while loop instead of just using the "for d" ?

[Jessica Flores]

I didn't originally have the while loop. Here's the revised code.

fileFolder = '/Users/JessFlores/Documents/MATLAB/XrayEye_*.tiff';
tiffpattern = fullfile(fileFolder, 'XrayEye_*.tiff');
dirOutput = dir(tiffpattern);
fileNames = cellfun(@(s) fullfile(fileFolder, s), {dirOutput.name}, 'Uniform', 0) .';
numFrames = numel(fileNames)
for p = 2:numFrames
    sequence(:,:,p) = imread(fileNames{p}); 
end
sequenceNew = stdfilt(sequence, ones(3));

      for d = 0001:numFrames
          %intensity plot info
          sequence(:,:,p);
          [x,y] = size(sequence(:,:,p));
          X = 1:x;
          Y = 1:y;
          [xx,yy] = meshgrid(Y,X);
          i = im2double(sequence(:,:,p));

          %maz and FWHM x and y from plot
          [mz,k] = max(i(:));
          [ix,jy] = ind2sub(size(i),k);
          [ix,jy,mz];

          %max intensity coordinates in microns
          %pixel size is 6.45x6.45 microns
          ix = 6.45*ix;
          jy = 6.45*jy;
          [ix,jy,mz];
          mhz = 0.5*mz;
          [ihx,jhy]=ind2sub(size(i),k);
          [ihx,jhy,mhz];
          FWHMx = abs(ix-ihx)*2;
          FWHMy = abs(jy-jhy)*2;

          %table of data from plot
          f = figure;
          data = [ix,jy,mz,FWHMx,FWHMy];
          colnames = {'X(Zmax) [um]', 'Y(Zmax) [um]', ...
              'Zmax/Max Intensity', 'FWHMx [um]', 'FWHMy [um]'};
          t = uitable(f, 'Data', data, 'ColumnName', colnames, 'ColumnWidth', {120});
          t.Position(3) = t.Extent(3);
          t.Position(4) = t.Extent(4);

          %make subplot of table and intensity plot
          subplot(2,1,1) = mesh(xx,yy,i); colorbar;
          title('Intensity of Beam Image as a Function of Pixel Position');
          xlabel('X(pixels)');
          ylabel('Y (pixels)');
          zlabel('Intensity (double precision corrected) [a.u.]');

          pos = get(subplot(2,1,2, 'position'));
          delete(subplot(2,1,2))
          set(t,'units','normalized')
          set(t, 'position', pos)
      end

[Walter Roberson]

You start reading files "for p = 2", skipping the first file, but you use the data in that spot because you start the frame count "for d = 1"

There is no point in calculating the following values in your "for d" loop:

          [x,y] = size(sequence(:,:,p));
          X = 1:x;
          Y = 1:y;
          [xx,yy] = meshgrid(Y,X);

Those must be the same for all images.

You calculate

[ix,jy] = ind2sub(size(i),k);

and later

   [ihx,jhy]=ind2sub(size(i),k);

at a time when you have not changed "k" or "i". Your ihx and jhy are therefore going to be the same as what you calculated at the time of the first call. But in meantime you have multiplied ix and jy by 6.45. Then you do

FWHMx = abs(ix-ihx)*2;

which reflects the ix that has been multiplied by 6.45 and subjects the one that has not. The difference is going to be (6.45-1)*ihx, which you then multiply by 2, giving 10.90 * ihx. And you call that your FWHMx. But at no point have you exactly searched the data to find the place that has half magnitude for z.

As I wrote above,

You should be searching for a location whose value is closest to mhz (the half maximum) and then the distance of that point to the maxima gives you the full width at half maxima. That would look like,

[mz,k] = max(i(:));
mhz = 0.5*mz;
[ix,iy] = ind2sub(size(i),k);
[minhalfdiff, poshalfdiff] = min(abs(i - mhz));
[ihx, ihy] = ind2sub(size(i), poshalfdiff);
halfwidth = sqrt((ihx - ix).^2 + (ihy - iy).^2);

Now convert ix, iy, ihx, ihy, halfwidth to microns by multiplying by the 6.45 for presentation purposes.

Answer 2

[Thorsten]

As far as I understood is the problem not the algorithm but that the loop seems to run only once.

First I would simply add a

 d

within the for-loop to check whether the loop really runs only once. I tested a similiar program and it worked well.

Next you could add a

 pause

after the subplot to see whether there are some subplots that are simply not displayed because the loop keeps running faster than the time needed for display.

You may then check Walter's comments on your algorithm.