# Convert 3D contour point cloud to binary voxel map

**URL:** https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214
**Category:** Support
**Created:** [September 14, 2020, 4:11pm UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214 "2020-09-14T16:11:37Z")
**Posts on this page:** 9
**Page:** 1

<div class="post-metadata">

### Author: ![ferdymercury](https://discourse.vtk.org/user_avatar/discourse.vtk.org/ferdymercury/32/3870_2.png) [@ferdymercury](https://discourse.vtk.org/u/ferdymercury)
#### Post date: [September 14, 2020, 4:11pm UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/1 "2020-09-14T16:11:37Z")

</div>

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://github.com/Slicer/Slicer/blob/ca118c83f699e035e4fb712cdffb2888cb237a58/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L164>

[https://vtk.org/Wiki/VTK/Examples/Cxx/Utilities/PointInPolygon](https://vtk.org/Wiki/VTK/Examples/Cxx/Utilities/PointInPolygon)

> <https://stackoverflow.com/questions/52699495/check-if-a-point-is-inside-and-object-using-vtk>

> [@How to convert from 2D polygons (vtkPolydata) to a 3D volumen](https://discourse.vtk.org/t/how-to-convert-from-2d-polygons-vtkpolydata-to-a-3d-volumen/2175):
>
> How to convert from 2D polygons (vtkPolydata) to a 3D volumen? Check the following image: [polygons]

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.

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [September 14, 2020, 5:27pm UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/2 "2020-09-14T17:27:04Z")

</div>

RT structure sets are horrible. We worked on developing custom VTK filters for converting it to closed surface in [SlicerRT](http://slicerrt.github.io/) (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.

---

<div class="post-metadata">

### Author: ![ferdymercury](https://discourse.vtk.org/user_avatar/discourse.vtk.org/ferdymercury/32/3870_2.png) [@ferdymercury](https://discourse.vtk.org/u/ferdymercury)
#### Post date: [September 16, 2020, 7:38am UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/3 "2020-09-16T07:38:36Z")

</div>

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?

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [September 16, 2020, 12:29pm UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/4 "2020-09-16T12:29:56Z")

</div>

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:

> **[SlicerRt/SlicerRT](https://github.com/SlicerRt/SlicerRT/tree/master/BatchProcessing)**
>
> Open-source toolkit for radiation therapy research, an extension of 3D Slicer. Features include DICOM-RT import/export, dose volume histogram, dose accumulation, external beam planning (TPS), struc...

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.

---

<div class="post-metadata">

### Author: ![ferdymercury](https://discourse.vtk.org/user_avatar/discourse.vtk.org/ferdymercury/32/3870_2.png) [@ferdymercury](https://discourse.vtk.org/u/ferdymercury)
#### Post date: [September 20, 2021, 9:51pm UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/5 "2021-09-20T21:51:11Z")

</div>

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:

```auto
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:

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

```

which matches well the CT size.

Thanks in advance.

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [September 20, 2021, 11:10pm UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/6 "2021-09-20T23:10:00Z")

</div>

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

---

<div class="post-metadata">

### Author: ![ferdymercury](https://discourse.vtk.org/user_avatar/discourse.vtk.org/ferdymercury/32/3870_2.png) [@ferdymercury](https://discourse.vtk.org/u/ferdymercury)
#### Post date: [September 20, 2021, 11:31pm UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/7 "2021-09-20T23:31:05Z")

</div>

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](https://discourse.vtk.org/uploads/default/original/2X/d/d43d4bd6d64c2b13cf7dd54bde7c12f922645650.png)

The same nrrd file opened in an online viewer:  
 ![image](https://discourse.vtk.org/uploads/default/original/2X/e/e31145455ee0740b056aa2afbe17db3f7c90730d.png)

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

```auto
        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`.

---

<div class="post-metadata">

### Author: ![May922](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/m/f05b48/32.png) [@May922](https://discourse.vtk.org/u/May922)
#### Post date: [December 29, 2022, 8:18am UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/8 "2022-12-29T08:18:13Z")

</div>

Hello Andras, I want to use SlicerRT in js, how can i use it ?

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lassoan/32/50_2.png) [@lassoan](https://discourse.vtk.org/u/lassoan)
#### Post date: [December 31, 2022, 2:56pm UTC](https://discourse.vtk.org/t/convert-3d-contour-point-cloud-to-binary-voxel-map/4214/9 "2022-12-31T14:56:44Z")

</div>

I would recommend to start using SlicerRT using its GUI in 3D Slicer. Once you confirmed that it does what you need you can start using it via Python scripting or C++ (there is a complete example [here](https://github.com/PerkLab/SlicerSandbox/blob/master/ImportOsirixROI/ImportOsirixROI.py) for using SlicerRT for creating a solid segmentation object from a set of contours - for importing Osirix ROI files).
