Calculate the mean of nonzero pixels

I have a 3D volume and a 3D binary mask. I want to calculate the mean of the pixels which correspond to the regions where the mask is equal to 1.

I thought to calculate the volume with the mask in order to have only the regions I am interested and then I calculated the mean.

I wrote the following piece of code but the mean value is extremely small and I think it's because I am taking into account the whole 3D volume (even the areas where the volume is zero) and not only the regions which I am interested.

r=volume; 
for i=1:size(r,1)
   for j=1:size(r,2)
      for k=1:size(r,3)
          masked_volume(i,j,k) = volume(i,j,k)*mask(i,j,k);
       end
   end
end
mean = mean(masked_volume(:))    

[Guillaume]

Note: there was an issue that prevented any sort of editing/commenting on this question and its answers. The issue has been fixed now.

Answer 1

[Guillaume]

The loop was pointless to start with. All you're doing is:

m = mean(volume .* mask, 'all');  %don't use mean as a variable name!

which indeed averages a lot of 0s. The proper way to do this:

m = mean(volume(mask));  %assuming that mask is of type logical

[Gina Carts]

the mask is not type logical. Can I do logical(mask) to convert it?

[Guillaume]

Indeed if the mask is double, you'll have to convert it to logical with logical(mask).

[Guillaume]

Following the question in the other answer, "But I would like to exclude also the pixels within the mask with values <0.001. Any help how to do that?"

Simple:

m = mean(volume(logical(mask) & volume >= 0.001));

[Gina Carts]

I would like to ask something more. How can I also add a threshold in the calculating mean?

i.e. I want to calculate the mean of the non-zero pixels as I did before using m = mean(volume(mask)). But I would like to exclude also the pixels within the mask with values <0.001. Any help how to do that?

[Gina Carts]

If the mask has regions with value 1, 2 and 3 how can I calculate the mean of the regions labeled with 1? And then mean for regions labeled with number 2 and 3 respectively?

When I load my mask I have 0 background, and then regions labeled with 1, 2 and 3. If I use the logical function to make my mask logical I only have 0 and 1.

[Guillaume]

This requires a different approach altogether. You need to use one of the aggregating functions splitapply, accumarray or groupsummary. Assuming your mask is not a mask but a label image with integer values from 0 to N:

objectsmean = accumarray(double(mask(:))+1, yourimage(:), [], @mean);
%or
objectsmean = splitapply(@mean, yourimage(:), double(mask(:))+1);
%or
objectsmean = groupsummary(yourimage(:), mask, 'mean');

accumarray is probably the fastest. The double(..) is here in case your mask is stored as an integer type and the number of objects is equal to intmax(class(mask)) which would cause an overflow when 1 is added.

In each case, the first value in the vector objectsmean will be the mean of the background.

Answer 2

[Cyrus Tirband]

Replace the last line by

meanval = mean(nonzeros(masked_volume));

That said, are mask and volume the same size? Those for loops are incredibly inefficient, and the whole snippet can be replaced by the one-liner if the dimensions of volume and mask are identical by using point-wise multiplication:

meanval = mean(nonzeros(volume.*mask));

[Gina Carts]

The mask and the volume have the same size. I have a question. I tried your suggested solution as well as the solution proposed by Guillaume. In order to use Guillaume suggestion I converted my mask to logical type using logical(mask). When I ran the two different solution you guys suggested I got different means and standard deviantions..

[Guillaume]

This solution multiplies the volume with the mask resulting in 0s everywhere that is mask (like your original code did). It then removes all the 0s, that is the 0s that are the result of masking but also all the 0s in the non-masked area. So yes, you'll get different results unless there's no 0s in the non-masked area.

My solution simply discards the masked pixels (with logical indexing) and average what's left so will give you the correct mean in all cases.

[Guillaume]

You're commenting on the wrong answer, making the conversation a bit more difficult.

Answer to that as a comment to my answer.