MATLAB: Split Volume Along 2D Surface

segmentsplitvectorizevolume

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

Best Answer

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
Related Question