calculating error while averaging multiple images

Hi, I have a simple problem and needs guidance. I have 15 image files. I want to calculate the average of these images and calculate error on these images. Any help is appreciated.

Answer 1

[Walter Roberson]

Do you want the average of the R, G, and B separately for each image?

Are you asking for averages to be taken over time for each pixel's R, G, B components?

Read the images into a 4D array, rows by columns by colorpanes by image_number . Then mean(Stack,4) or std(Stack, [], 4)

[Sumera Yamin]

Hi, i do not understand your point. i am not very much into image analysis and its terminologies, but i would say what i want to do is average out 10 images and calculate the position of peak intensity in the image with std as the error in the peak average. In terms of R,G and B, i guess i do not want average of R,G and B for each image seperately, rather average of peak intensity/pixel number.

[Walter Roberson]

projectdir = '.';   %directory images are in
dinfo = dir( fullfile(projectdir, '*.png'));
filenames = fullfile(projectdir, {dinfo.name});
nfiles = length(filenames);
image_contents = cell(nfiles,1);
for K = 1 : nfiles
    image_contents{K} = im2double(imread(filenames{K}));
end

fail = false;
try
    image_array = cat(4, image_contents{:});
catch ME
    fail = true;
    fprintf('Images are not all the same size, cannot proceed');
end

if ~fail
    mean_img = mean(image_array, 4);
    std_img = std(image_array, [], 4);
    gray_mean_img = rgb2gray(mean_img);
    [peak_gray_value, peak_location_idx] = max(gray_mean_img, 'all');
    [peak_row, peak_col] = ind2sub(size(gray_mean_img), peak_location_idx);
end

... but I have a suspicion that when you talk about peak intensity, that you mean you have images that represent graphs of some kind, and that you want to analyze the graphs represented in the images, extracting the data from them.

[Sumera Yamin]

Hi, many thanks for your answer. when i run your code it gives me following error.

"Error using rgb2gray>parse_inputs (line 80)

MAP must be a m x 3 array.

Error in rgb2gray (line 52)

isRGB = parse_inputs(X);

Error in Untitled (line 19)

gray_mean_img = rgb2gray(mean_img);"

you are right. i do have images with gaussian distribution and i want to extract data from these distributions meaning center of mass (or call it peak intensity), and the average of multiple images would give me error/standard deviation from this center of mass.

[Walter Roberson]

    gray_mean_img = mean_img;

Apparently you already have grayscale instead of RGB.

You should post a sample image.

[Sumera Yamin]

Sorry for my lack of understanding, this again gives me following error

Error using max

MAX with two matrices to compare and two output arguments is not supported.

Error in Untitled (line 20)

[peak_gray_value, peak_location_idx] = max(gray_mean_img, 'all');

[Walter Roberson]

    [peak_gray_value, peak_location_idx] = max(gray_mean_img, [], 'all');

[Sumera Yamin]

thanks, it works now

[Sumera Yamin]

what does 4 in the mean_img = mean(image_array, 4); stand for.

Also if i have a function like as follows. It is calculating the average of the images. will it be possible to add the calculation of standard deviation in ythe same function?

function [imageData] = meanImage(path, fileNamePrefix, fileNameSufix, n, mStart, mEnd)
    N = mEnd - (mStart - 1);
    firstImage = loadImage(path, fileNamePrefix, fileNameSufix, n, mStart);
    for mm = (mStart+1):mEnd
        currentImage = loadImage(path, fileNamePrefix, fileNameSufix, n, mm);
        firstImage = firstImage + currentImage;
    end
    imageData = firstImage ./ N;
end

[Walter Roberson]

The 4 means dimension #4.

The code I posted did

    image_array = cat(4, image_contents{:});

so it took all of the images and put them together along the 4th dimension. Then when you are taking the mean along the 4th dimension, you are taking the mean over all of the images.

I put the images together along the 4th dimension because I did not assume that they are grayscale images.

[Sumera Yamin]

i see, the code gives me mean image in terms of full matrix but i want to extract the mean and standard deviation of all the of image in terms of numbers of peak intensity.

[Walter Roberson]

function [imageData, imgstd] = meanImage(path, fileNamePrefix, fileNameSufix, n, mStart, mEnd)
    N = mEnd - (mStart - 1);
    currentImage = double(loadImage(path, fileNamePrefix, fileNameSufix, n, mStart));
    imgsum = currentImage;
    imgsqsum = imgsum.^2;
    
    for mm = (mStart+1):mEnd
        currentImage = double(loadImage(path, fileNamePrefix, fileNameSufix, n, mm));
        imgsum = imgsum + currentImage;
        imgsqsum = imgsqum + currentImage.^2;
    end
    mean_img = imgsum ./ N;
    imageData = cast(mean_img, class(currentImage));
    imgstd = sqrt((imgsqsum - 2 .* mean_img .* imgsum + mean_img.^2) ./ (N-1));
end

This code deliberately turns the mean image back into the datatype of the individual images, such as uint8. This is because it calculates the mean image in the same units as the original image, not in terms of the 0 to 1 range that could be used instead. Whether this is best for you depends on what you want to do with the mean image.

[Sumera Yamin]

so i took this part of code and did the furthur processing of my data. combining with the above function, the code is as shown below. The problem is that for i want to use both imageData and imgstd for calculating [maxzeile, maxspalte] and [maxzeile_err, maxspalte_err] respectively. My code however takes imageData for calculating both and hence giving me the same output. Could you please pointout how should i change the code to use both imagedata and imgstd.

function [maxzeile, maxspalte, maxzeile_err, maxspalte_err, xRMS, yRMS] = findBeam(beamIn, windowsize, threshold, maxIterations)

%rough beam position
BEAMFOUND = 0;
WINDOWSIZE = windowsize/2;
THRESHOLD = round(threshold*WINDOWSIZE^2);
n = 0;
beamImageSize = size(beamIn);

while (BEAMFOUND == 0 && n <= maxIterations)

	[maxwerte, maxzeilen] = max(beamIn); 
	[maxwert, maxspalte] = max(maxwerte);
	maxzeile = maxzeilen(maxspalte);

   if (maxzeile <= WINDOWSIZE)
      yLow = 1;
   else
      yLow = maxzeile-WINDOWSIZE;
   end

   if (maxzeile + WINDOWSIZE >= beamImageSize(1))
      yHigh = beamImageSize(1);
   else
      yHigh = maxzeile+WINDOWSIZE;
   end

   if (maxspalte <= WINDOWSIZE)
      xLow = 1;
   else
      xLow = maxspalte-WINDOWSIZE;
   end

   if (maxspalte + WINDOWSIZE >= beamImageSize(2))
      xHigh = beamImageSize(2);
   else
      xHigh = maxspalte+WINDOWSIZE;
   end

   yLow = abs(round(yLow)); %disp(yLow)
   yHigh = abs(round(yHigh)); %disp(yHigh)
   xLow = abs(round(xLow)); %disp(xLow)
   xHigh = abs(round(xHigh)); %disp(xHigh)
   %disp('------')
   
   window = beamIn(yLow:yHigh,xLow:xHigh); %figure(n+5); imagesc(window);
    
   % Intensity cut
   thre = 0.75 * max(max(window));
   beamArea = sum(sum(window>thre));

   if (beamArea >= THRESHOLD)
		BEAMFOUND = 1;
   else
	beamIn(yLow:yHigh,xLow:xHigh) = 0; %beamIn(yLow:yHigh,xLow:xHigh)>10000;
   end
   n = n + 1;
end

[xRMS, yRMS] = beamRMS(window);

   beamImageSize = size(beamIn(yLow:yHigh,xLow:xHigh));
   beamX = 1:beamImageSize(2);
   beamY = 1:beamImageSize(1);
   totalMass = sum(sum(beamIn(yLow:yHigh,xLow:xHigh)));
   centerOfMass = zeros(1,2);
   for i=1:length(beamX)
      for  j=1:length(beamY)
      	centerOfMass = centerOfMass + double(window(j,i)).*double([beamY(j),beamX(i)]);
      end
   end
   centerOfMass = centerOfMass / totalMass;
   centerOfMass = round(centerOfMass);
   centerOfMass = centerOfMass + [yLow,xLow];
     maxzeile = centerOfMass(1);
   maxspalte = centerOfMass(2);
%error in center of mass:
   beamImageSize_err = size(beamIn(yLow:yHigh,xLow:xHigh));
   beamX_err = 1:beamImageSize_err(2);
   beamY_err = 1:beamImageSize_err(1);
   totalMass_err = sum(sum(beamIn(yLow:yHigh,xLow:xHigh)));
   centerOfMass_err = zeros(1,2);
   for i=1:length(beamX_err)
      for  j=1:length(beamY_err)
      	centerOfMass_err = centerOfMass_err + double(window(j,i)).*double([beamY_err(j),beamX_err(i)]);
      end
   end
   centerOfMass_err = centerOfMass_err / totalMass_err;
   centerOfMass_err = round(centerOfMass_err);
   centerOfMass_err = centerOfMass_err + [yLow,xLow];
   maxzeile_err = centerOfMass_err(1);
   maxspalte_err = centerOfMass_err(2);
end

Answer 2

[Image Analyst]

See my attached demo. It averages RGB and gray scale images and gives statistics about them.. Adapt as needed.

[Sumera Yamin]

Many thanks for your demo script. It indeed does the average, but it does not give the error/standard deviation on average, which is my original prob

[Image Analyst]

You can use my camera noise analyzer program to do that. I've attached it. Browse to your folder of images then select some and analyze them. Let me know if anything on it does not work.

[Image Analyst]

Sumera, regarding your "Answer", your attached screenshot does not show an error. Does it say something in the command window? Usually I also print the error to the command window as well as popping it up in a message box. Clear Mask lets you remove any mask you'd drawn over a subimage with the Step 3 button so that you are again analyzing the entire window instead of just analyzing inside the drawn mask area.

[Sumera Yamin]

It says it says error in beam analyze call back. error reported by matlab is index exceeds matrix dimensions. i have attached the image of sccreenshot with error.

[Image Analyst]

Are all your images RGB and the same size? Please attach two of them so I can see what went wrong, because for me, with a folder of RGB images all of the same size, it went fine.

[Sumera Yamin]

yes all my images are of smae sizes but grey scale. please see the attached images

[Image Analyst]

I've adapted it to handle gray scale images. See attatchment.

[Sumera Yamin]

hi, unfortunately, still gives me following error

[Walter Roberson]

Replace isfolder() with isdir() on older MATLAB.

[Sumera Yamin]

hi, thanks for your answer, this error gone now, but gui give me another (attached) message, however if i ignore the message, the gui does analyze the images.

[Image Analyst]

See attached.

Answer 3

[Sumera Yamin]

many thanks for this analyzer, it gives me attached error. Also what is mask and what does clear mask does?