Create mask for 3D image based on plane

Algorithms
1

Hi there,

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?

2

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.

3

Hello @jbi35,

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
dimension=3

segmentation = sitk.GetImageFromArray(np.random.random([512]*dimension)) > 0.8
segmentation.SetOrigin([-100]*dimension)
segmentation.SetSpacing([0.1]*dimension)

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(),
                                                   origin=segmentation.GetOrigin(),
                                                   spacing=segmentation.GetSpacing(),
                                                   direction=segmentation.GetDirection())
#(v-point_on_plane)*plane_normal
xyz = [(sitk.VectorIndexSelectionCast(segmentation_coordinates,i)-point_on_plane[i])*plane_normal[i] for i in range(segmentation_coordinates.GetDimension())]

projection_on_normal = xyz[0]
for i in range(1,segmentation_coordinates.GetDimension()):
  projection_on_normal += xyz[i]
segmentation_on_one_side = (projection_on_normal>0)*segmentation


sitk.Show(segmentation*255)
sitk.Show(segmentation_on_one_side*255)
1 Like
4

Dear @zivy,

thanks a lot for your help!! Your approach works like charm and is super fast.

Best

1 Like