Centroid calculation image stack (for loop)

I want to calculate the total centroid of an image stack. For a single image i could use the regionprops() to calculate it as shown here. For a single image my variable c1 is where the coordinates of all the centroids in the image are saved To calculate the total centroid I just summed all coordinates in the image up and divide it by amount of centroids in the picture. The first row should display the x and y coordinates of my total centroid for the first image, the second row for the second image, etc.

% single image 
Ilabel = bwlabel(Ibw);
stat1 = regionprops(Ilabel,'Area','centroid');
c1 = cat(1, stat1.Centroid); % c1 is 16x2
[m1, n1] = size(c1);  % m1= 16 amount of centroids in a single image
totalcentroid1 = sum(c1)/m1; 

For an image stack my idea was to write a for loop to do this process. My first problem is when i save all centroids for all imagestacks in a variable , that i don't know the exact amount of centroids in a single image to calculate the total centroid of each image. How can i do this ? In the following i tried to implement the calculation of a total centroid for an imagestack. There is error message in this code. How can i fix this issue

%% imagestack 

for i=0:options.ImageAmount-1
    Ilabel(:,:,i+1) = bwlabel(Ibw(:,:,i+1));
    stat1(:,:,i+1) =regionprops(Ilabel(:,:,i+1),'centroid');
    %c1(:,:,i+1) = cat(1,stat1.Centroid(:,:,i+1);   % unsure aboute this part 
    %[m1, n1] = size(c1(:,:,i+1)); 
    % totalcentroid1(:,:,i+1) = sum(c1(:,:,i+1)/m1(:,:,i+1); 
end

%%
Subscripted assignment dimension mismatch.

Error in delete (line 126)
    stat1(:,:,i+1) =regionprops(Ilabel(:,:,i+1),'centroid');
 

Answer 1

[Image Analyst]

Your formula for a single image is not even correct. The centroid of a binary image with 16 blobs is not the average of the 16 centroids of course.

Why don't you skip regionprops and just use your binary image itself? Since

mean x = (sum of x*graylevel at x) / (sum of graylevels at x)

or in other words, the mean x is simply the mean of all the x values. So

for sliceNumber = 1 : size(Ibw, 3)
        [rows, columns] = find(Ibw(:, :, sliceNumber));
        xCentroids(sliceNumber) = mean(columns);
        yCentroids(sliceNumber) = mean(rows);
end
meanXCentroid = mean(xCentroids)  % Mean column, x
meanYCentroid = mean(yCentroids)  % Mean row, y

[Lightisthenight]

Thank you sir for your fast reply! I still have a question about your third and forth line in the code :

xCentroids(sliceNumber) = mean(columns);
yCentroids(sliceNumber) = mean(rows);

A binary image means that each pixel is stored as a single digit 0 or 1. If I'm not mistaken the find() function returns all the nonzero elements a.k.a ones in my array Ibw. With xCentroids(1)= mean(columns) it should give me the average value of my array which is the x coordinate of the total centroid of my first image right? And meanXCentroid = mean(xCentroids) gives me the average x coordinate of the total centroid of all my images.Or does xCentroids(sliceNumber) = mean(columns) give me all the centroids in the first image ?

[Image Analyst]

find() gives you the row and columns of the non-zero pixels (the white pixels). So averaging those will give you the centroid. It does not matter if there is one blob, two blobs, thirty thousand blobs or whatever. In essence, it's like each pixel is its own blob. The mean x is this formula

xMean = sum(x*value) / sum(value)

That's in general, not just for images. You can see the values sort of "cancel out" and you're left with x. If you wanted the mean gray level, then you'd divide by the sum of the x instead of value:

meanGL = sum(x*GL)/sum(x)

See how the x "cancels out" and you're left with units of gray level. Anyway, that's not what you're asking, I just think it helps to visualize what gets "canceled out" to figure out what you need to divide by.

So, the mean is essentially doing the formula I gave first. It really doesn't matter if you have 1, 2, 30 thousand or a million blobs, or what their size it. It will give you the mean x and y. Just be sure that you don't make that all too common beginner mistake and think that (x,y) is (row, column). Arrays are references yourArray(y, x) and NOT yourArray(x, y) because row (y) is the first index to an array, not x (column).

So I keps track of the centroids for each slice, then averaged those to get the mean centroid over ALL slices.

[Lightisthenight]

Ok thank you for your detailed explanation.

for sliceNumber = 1 : size(Ibw, 3)
        [rows, columns] = find(Ibw(:, :, sliceNumber));
        xCentroids(sliceNumber) = mean(columns);
        yCentroids(sliceNumber) = mean(rows);
        meanXCentroid = mean(xCentroids)  % Mean column, x
        meanYCentroid = mean(yCentroids)  % Mean row, y
end

So if I want to get the mean centroid for each slice then i just need to include the two last lines in the for loop right?

[Image Analyst]

No, leave them outside the loop just as I had it. The mean centroids for each slice are stored in the xCentroids and yCentroids arrays.

[Lightisthenight]

Thank you very much sir!