Create mask for 3D image based on plane

Algorithms
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