How to calculate mean in moving cubic volume?

I have a stack of 2D images in 3D array format 100x100x500. Now I want to select a small cubic volume of 1x1x1 at center of the stack and increase this cubic volume step by step till I reach the array size, calculating mean at each stage. I know I need to use mat2cell and mean function. How can be done in corect for loop?

Answer 1

[Matt J]

I know I need to use mat2cell and mean function.

I don't think you do. Let's call your 3D volume, V:

[ci,cj,ck]=deal(51,51,251); %center coordinate ?

[I,J,K]=ndgrid(abs((1:100)-ci),abs((1:100)-cj),abs((1:500)-ck));

L=max(max(I,J),K)+1; %label matrix

stats = regionprops3(L,V,'MeanIntensity','Volume');

cumVolumes = cumsum(stats.Volume);
cumIntensities = cumsum(stats.Volume.*stats.MeanIntensity);

result = cumIntensities ./cumVolumes

[SojM]

Thanks for your quick reply. I do get a table with all meanvalues with your code. But when I checked the mean of entire volume by using mean(V(:)), it does not match with mean intensity of last value from the result table from your code.

[Matt J]

I've modified the code.

[SojM]

Perfect. This gives me exactly what I was looking for. Thanks.

[SojM]

Hi Matt, I am trying to use this code in for loop by changing the 'Z center coordinate (50 to 300)' keeping x and y center coordinate constant. I also want to change grid in Z such that it stays 100. How would I save all the results in struct array for all "i" in the loop? Could you please help me to fix following code?

I tried this, but this does not seem right

for i = 1:400
    
[ci,cj,ck]=deal(51,51,50+i); %center coordinate
[I,J,K]= ndgrid(abs((1:100)-ci),abs((1:100)-cj),abs((i:100+i)-ck)); % Z remain 100 but changes with center
L = max(max(I,J),K)+1; %label matrix
stats = regionprops3(L,V(:,:[i:100+i]),'MeanIntensity','Volume');
cumVolumes{i} = cumsum(stats.Volume);
cumIntensities{i} = cumsum(stats.Volume.*stats.MeanIntensity);
result{i} = cumIntensities ./cumVolumes

end 

[Matt J]

Something like this?

clear S
S(400)=struct('cumVolumes',[],'cumIntensities',[],'result',[]); %pre-allocate

for i = 1:400
    
    [ci,cj,ck]=deal(51,51,50+i); %center coordinate
    [I,J,K]= ndgrid(abs((1:100)-ci),abs((1:100)-cj),abs((i:100+i)-ck)); % Z remain 100 but changes with center
    L = max(max(I,J),K)+1; %label matrix
    stats = regionprops3(L,V(:,:[i:100+i]),'MeanIntensity','Volume');
    
    
    S(i).cumVolumes = cumsum(stats.Volume);
    S(i).cumIntensities = cumsum(stats.Volume.*stats.MeanIntensity);
    S(i).result = cumIntensities ./cumVolumes

end 

[SojM]

Thanks Matt. Yes, I was looking for somethign like that. That was very useful. Your code gave me error that "Undefined function or variable 'cumIntensities' ". I needed to define them before using as struct fields. It is working perfectly now. Following is the modified version of your code:

clear S
S(400)=struct('cumVolumes',[],'cumIntensities',[],'result',[]); %pre-allocate
for i = 1:400
    
    [ci,cj,ck]=deal(51,51,50+i); %center coordinate
    [I,J,K]= ndgrid(abs((1:100)-ci),abs((1:100)-cj),abs((i:100+i)-ck)); % Z remain 100 but changes with center
    L = max(max(I,J),K)+1; %label matrix
    stats = regionprops3(L,V(:,:,[i:100+i]),'MeanIntensity','Volume');
    
    cumVolumes = cumsum(stats.Volume);
    cumIntensities = cumsum(stats.Volume.*stats.MeanIntensity);
    result = cumIntensities ./cumVolumes;
    
    S(i).cumVolumes = cumVolumes;
    S(i).cumIntensities = cumIntensities;
    S(i).result = result;
end 

Answer 2

[Amy Van Wey Lovatt]

Here is another option, which I'm using to calculate contrast uptake in breast images. I'm not sure which runs faster.

function PE = PeakEnhancement(S0,S1)

% S0 is the pre-contrast phase

% S1 is phase 1 or phase 2 of the contrast, this is an NxMxP matrix.

% sub is the initial percentage of uptake

sub=(S1-S0)./S0;

PE=zeros(size(sub));

% PeakEhnancement (PE) is the average of the 9 closest voxels in 3 x 3 x 3 cube.

% PE is a place holder ensuring size(PE)=size(S1) and boundaryies are all zerp.

for i=2:length(sub(:,1,1))-1

for j=2:length(sub(1,:,1))-1

for k=2:length(sub(1,1,:))-1

PE(i,j,k)=mean(sub(i-1:i+1,j-1:j+1,k-1:k+1),'all');

end

end

end

end