Convert 3D contour point cloud to binary voxel map

ferdymercury · September 14, 2020, 4:11pm

I am writing a small C++ program, which contains:

a voxel grid (nX * nY * nZ)
a voxel (i,j,k) to 3D point (x,y,z) conversion matrix
an array of 3D points (x,y,z) defining a closed surface

As I am new to VTK, I was looking for help to see what would be the most efficient / shortest way to compute a boolean 3D matrix that tells me which points of the 3D voxel grid are inside the contour, and which ones are not. I would prefer not to use external batch-processes, just C++ API libraries like VTK or ITK.

I was taking a look at different examples (see below), and wanted to ask about what would be in your opinion the best way to go.

github.com

Slicer/Slicer/blob/ca118c83f699e035e4fb712cdffb2888cb237a58/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L164


// Compute output labelmap geometry based on poly data, an reference image
// geometry, and store the calculated geometry in output labelmap image data
if (!this->UseOutputImageDataGeometry)
  {
  if (!this->CalculateOutputGeometry(closedSurfacePolyData, binaryLabelmap))
    {
    vtkErrorMacro("Convert: Failed to calculate output image geometry!");
    return false;
    }
  }
else if (outputGeometryLabelmap)
  {
  std::string geometryString = vtkSegmentationConverter::SerializeImageGeometry(outputGeometryLabelmap);
  vtkSegmentationConverter::DeserializeImageGeometry(geometryString, binaryLabelmap, false);
  }

// Allocate output image data
binaryLabelmap->AllocateScalars(VTK_UNSIGNED_CHAR, 1);

void* binaryLabelmapVoxelsPointer = binaryLabelmap->GetScalarPointerForExtent(binaryLabelmap->GetExtent());
if (!binaryLabelmapVoxelsPointer)

https://vtk.org/Wiki/VTK/Examples/Cxx/Utilities/PointInPolygon

Note: the original XYZ contour data are coming from a DICOM RTstruct, so they are arranged as a stack of 2D closed contour slices. Maybe the PointInPolygon2D function might be the most straightforward solution?

Thanks in advance.

lassoan · September 14, 2020, 5:29pm

RT structure sets are horrible. We worked on developing custom VTK filters for converting it to closed surface in SlicerRT (with support for keyhole technique, branching, smooth end-capping, etc.). If you only need to process special cases (e.g., no branching, no keyholes) then you can find lots of naive implementations out there, but if you need a full solution, which is expected to work for all data sets, then the only open-source solution that I’m aware of is in SlicerRT.

You can link MRML library and SlicerRT logic classes into your own application, but it is much simpler to not dissect these libraries but use as is. You can run the conversion from command-line, or use it from 3D Slicer’s Python environment (just another virtual Python environment where you can install any other Python packages, etc.), with an application GUI, without a GUI, or using Slicer as a Jupyter notebook kernel.

ferdymercury · September 16, 2020, 7:38am

Thanks for the swift reply and detailed info.

What would be the CLI command to run the MRML conversion from RTstruct.dcm to binary label map? I guess I need to pass the voxel grid dimensions?

lassoan · September 16, 2020, 7:30pm

See a script for converting a DICOM folder structure containing of RT Structure sets to labelmap volumes in NRRD format (and optionally corresponding CTs as well) and usage instructions here:

It is recommended to provide input CT images as well, because RT structure set does not contain slice spacing information, so without CT images the algorithm just tries to guess the spacing from the closest contour distance, which is of course not always optimal.