I am trying to solve the following problem. I have created a 3D segmentation with SimpleITK which I would like to split into two parts using a plane defined by the normal vector and point. My goal is to set the voxel-value of all voxels in the segmentation (binary mask) depending on whether the voxel is located above or below an arbitrary plane, i.e., a plane that is not perpendicular to the axes.
The only solution I can think of is to loop over all voxels, calculate its position in space, and check whether it is below or above the plane. However, I was hoping that there might be a more efficient approach. Does anyone have a hint for me which filters I could look into?
This should not be that slow. ITK already provides efficient and convenient method to calculate voxel’s position in space.
And to avoid it, you would need to have a significantly more complicated processing. E.g. for each scan-line find the voxel which intersects with the plane using binary search, then see which side of the scan-line is below plane, then mark all voxels of that part of scan-line accordingly, then mark the other part of scan-line the opposite way.
This was a fun exercise, see SimpleITK code below. Most of the code is data generation, and the rest is straightforward, the sign of (v-p)*n determines on which side of the plane (p,n) point v resides.
import SimpleITK as sitk
import numpy as np
# plane in 3D, line in 2D code works for both
segmentation = sitk.GetImageFromArray(np.random.random(*dimension)) > 0.8
point_on_plane = segmentation.TransformContinuousIndexToPhysicalPoint(np.array(segmentation.GetSize())/2.0)
plane_normal = np.array(segmentation.GetOrigin()) - np.array(point_on_plane)
plane_normal /= np.linalg.norm(plane_normal)
#get the coordinates of all the segmentation volume voxels
segmentation_coordinates = sitk.PhysicalPointSource(size=segmentation.GetSize(),
xyz = [(sitk.VectorIndexSelectionCast(segmentation_coordinates,i)-point_on_plane[i])*plane_normal[i] for i in range(segmentation_coordinates.GetDimension())]
projection_on_normal = xyz
for i in range(1,segmentation_coordinates.GetDimension()):
projection_on_normal += xyz[i]
segmentation_on_one_side = (projection_on_normal>0)*segmentation
thanks a lot for your help!! Your approach works like charm and is super fast.