How to calculate mean in moving cubic volume?

2 次查看（过去 30 天）
SojM
SojM 2020-7-25
回答： Amy Van Wey Lovatt 2022-9-11
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?

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2020-7-25
编辑：Matt J 2020-7-25
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
  6 个评论
显示 4更早的评论
Matt J
Matt J 2020-8-2
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
SojM 2020-8-2
编辑：SojM 2020-8-2
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

请先登录，再进行评论。

更多回答（1 个）

Amy Van Wey Lovatt
Amy Van Wey Lovatt 2022-9-11
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

类别

在 Help Center 和 File Exchange 中查找有关 3-D Volumetric Image Processing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by