Split Volume Along 2D Surface

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

[darova]

Can you show the volume and surface?

[Johannes Rebling]

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);

Answer 1

[Johannes Rebling]

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

Answer 2

[darova]

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

[Johannes Rebling]

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]

Im out of ideas. It was my best