Split Volume Along 2D Surface

3 次查看（过去 30 天）
Johannes Rebling
Johannes Rebling 2020-4-18
评论： darova 2020-4-19
Hi everyone,
I have what I feel is a simple problem but I just can't find a fast, vectorized way to split a 3d volume along the z indicies defined in a matrix, as described below.
I have a volume [X Y Z] in which in find an arbitrary surface. This surface has the dimension [X Y] and the value at each positon in said surface tells me where to split my volume along z. Here is the naive for-loop approach, but it's sloooow (I have many volumes and they are much larger than in this example).
volume = rand(500,600,300);
surface = round(rand(500,600).*150+1);
tic;
[topVol,botVol] = split_volume(volume,splitDepth);
done(toc);
function [topVol,botVol] = split_volume(volData,Z)
[nX,nY,~] = size(volData);
topVol = volData;
botVol = volData;
for iX = 1:nX
for iY = 1:nY
zIdx = Z(iX,iY);
topVol(iX,iY,zIdx:end) = 0;
botVol(iX,iY,1:zIdx) = 0;
end
end
end
  2 个评论
darova
darova 2020-4-18
Can you show the volume and surface?
Johannes Rebling
Johannes Rebling 2020-4-19
编辑：Johannes Rebling 2020-4-19
I am not sure what you mean. The volume can be any size, like in my working example aboe:
volume = rand(500,600,300);
Same for the "surface", it will always be the same size in x-y as the volume and can have any values between 1 and the number of elements in z in the volume (i.e. 1 to 300 in my example).
surface = round(rand(500,600).*150+1);

请先登录，再进行评论。

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

采纳的回答

Johannes Rebling
Johannes Rebling 2020-4-18
I have found a somewhat vectorized solution that is 15 times faster than the native for loops and still 5 times faster than the optimized for loops.
If anyone has an even faster way, let me know!
My code as it might be usefull to someone doing the same:
tempVol = rand(500,500,800);
[nX,nY,nZ] = size(tempVol);
splitPlane = uint16(randi([round(nZ/2)-10 round(nZ/2)+10],[nX, nY])); % line to split along...
tic;
botVolKeep1 = split_volume_naive(nX,nY,nZ,splitPlane);
toc;
tic;
botVolKeep2 = split_volume_for(nX,nY,nZ,splitPlane);
toc;
tic;
botVolKeep3 = split_volume_fast(nX,nY,nZ,splitPlane);
toc;
isequal(botVolKeep1,botVolKeep2)
isequal(botVolKeep1,botVolKeep3)
function [botVolKeep] = split_volume_naive(nX,nY,nZ,splitPlane)
botVolKeep = ones(nX,nY,nZ);
for iX = 1:nX
for iY = 1:nY
zIdx = splitPlane(iX,iY);
botVolKeep(iX,iY,1:zIdx) = 0;
end
end
end
function [botVolKeep] = split_volume_for(nX,nY,nZ,splitPlane)
botVolKeep = true(nX,nY,nZ);
for iY = 1:nY
for iX = 1:nX
zIdx = splitPlane(iX,iY);
botVolKeep(iX,iY,1:zIdx) = false;
end
end
end
function [botVolKeep] = split_volume_fast(nX,nY,nZ,splitPlane)
linZIdx = uint16(repmat(1:nZ,[nY,1])); % array from 1:nZ with nY rows
botVolKeep = false([nY,nZ,nX]);
for iX = 1:nX
splitLine = splitPlane(iX,:)';
% t = linZIdx > splitLine;
botVolKeep(:,:,iX) = linZIdx > splitLine;
end
botVolKeep = permute(botVolKeep,[3 1 2]);
end

更多回答（1 个）

darova
darova 2020-4-19
See this madness
clc,clear
[x1,y1,z1] = meshgrid(1:10); % volume boundaries
[x2,y2] = meshgrid(1:0.5:10); % surface boundaries
[T,R] = cart2pol(x2,y2);
z2 = 10*sin(R)./R+5; % create surface
z11 = interp2(x2,y2,z2,x1,y1); % find 'Z' value for each volume point
ix = z1 > z11; % half of a volume
surf(x2,y2,z2)
hold on
plot3(x1(ix),y1(ix),z1(ix),'.r')
plot3(x1(~ix),y1(~ix),z1(~ix),'.g')
hold off
axis vis3d
  2 个评论
Johannes Rebling
Johannes Rebling 2020-4-19
Thanks a lot for taking the time and providing a fully vectorized solution. It does indeed work and produces the same results as the naive double for-loop, but is only marginally faster and actually much slower than my other approaches.
For a800 x 800 x 400 volume (actually still smaller than the volumes I am processing), here is what I get on my fairly decent PC.
Naive for-loop implementation took 3.230 s.
Optimized for-loop implementation took 0.822 s (3.9 x faster).
Partially vectorized implementation took 0.292 s (11.1 x faster).
Fully vectorized implementation took 2.335 s (1.4 x faster).
Unfortunately, the meshgrid in 3d becomes quite slow and the fully vectorized approach also uses a lot more memory.
darova
darova 2020-4-19
Im out of ideas. It was my best

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by