Convert 3D contour point cloud to binary voxel map

Support
#1

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.


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.

#2

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.

#3

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?

#4

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.

1 Like
#5

Hi Andras, I have a related question.

I managed to call the Python script you suggested from my C++ application, and then reading the NRRD back in C++. However, I am unsure what would be the best way to convert the NRRD file content into a raw boolean 3D matrix with the dimensions of the CT.

bool roi[61][61][46];

where roi[i][j][k] would tell you if voxel i,j,k is part of the ROI or not.

So far, I am doing:

vtkSmartPointer<vtkNrrdReader> reader = vtkSmartPointer<vtkNrrdReader>::New();
reader->SetFileName(fileName);
reader->Update();
vtkImageData* outp = reader->GetOutput();
int* x = (int*)outp->GetScalarPointer();
std::cout << outp->GetDimensions()[0] << " " << outp->GetDimensions()[1] << " " << outp->GetDimensions()[2] << std::endl;

for which I get:

"Target.nrrd"
61 61 46
"rectum.nrrd"
61 61 46
"urethra.nrrd"
61 61 46

which matches well the CT size.

Thanks in advance.

#6

You can extract a binary volume of a specific label value from a labelmap volume using simple thresholding.

#7

Thanks for the support! The thresholding works.

Now I am only getting a weird effect, as one slice of the ROI seems to contain 4 tiled slices

image
image608Ă—633 7.11 KB

The same nrrd file opened in an online viewer:
image

I guess I must be doing something wrong with the indices ?

        int* data = (int*)outp->GetScalarPointer();
        const int nX = outp->GetDimensions()[0];
        const int nY = outp->GetDimensions()[1];
        const int nZ = outp->GetDimensions()[2];
        std::cout << nX << " " << nY << " " << nZ << std::endl;
        const int threshold = 0;
        for(int k = 2; k < nZ ; ++k)
        {
            for(int j = 0; j < nY ; ++j)
            {
                for(int i = 0; i < nX ; ++i)
                {
                    std::cout << ( data[ i + j*nX + k*nX*nY ] > threshold ) << " ";
                }
                std::cout << std::endl;
            }
            std::cout << std::endl;
            break;
        }

EDIT: Ok, found the bug. It was that I was assuming int for “data”, but it should be shorts.