How does slicer process the DICOM RTSTRUCTURE for pyradiomics

[shuo0220]

@fedorov Hi Dr. Fedorov, I'm wondering what slicer does to process the DICOM RTStructure for pyradiomics. I'm asking because I was using dcmrtstruct2nii (https://pypi.org/project/dcmrtstruct2nii/) to convert dicom rtstructure to a nii file, and used dcm2nifti to convert the dicom CT image to a single nii file. When I calculated using pyradiomics directly on these nii files, i found some of the feature calculation were different from the ones calculated using slicer+pyradiomics on DICOMs. I wonder if there is a standard package that can convert dicom images + rtstructure to a sITK readable file and use pyradiomics directly? Thank you so much!

Originally posted by @fedorov in https://github.com/Radiomics/pyradiomics/issues/395#issuecomment-405599780:

@lisherlock no, you cannot use the DICOM files directly, unless you are processing a single DICOM file. SimpleITK does not sort the DICOM files in geometric order. If you are working with a volumetric acquisition, you should first preprocess the DICOM series to generate a volume saved in any of the formats readable by ITK, and then use that as the input image for pyradiomics.

You can use dcm2niix, or plastimatch convert to geometrically sort files in a DICOM series and save them as a volume.

If your input is multiple series, you can first sort the files into separate series using the dicomsort tool.

[fedorov]

When you are reading DICOM RTSS in Slicer, reconstructions is handled by the SlicerRT extension. If I understand correctly, SlicerRT is using plastimatch for reading DICOM RT-structs, but then also implements some additional heuristics for constructing close surface object. We had a discussion with Andras Lasso (@lassoan) and Csaba Pinter (@cpinter) over email (unfortunately, not in a public forum; email thread around Oct 7, 2017) about the differences, and here's the summary from Andras:

Conversion from planar contours to binary labelmap is a very complex operation, as in a general case you need to handle inter-plane interpolation/splitting/merging, intra-plane splitting, end-capping, keyhole resolution, etc. Plastimatch may be able to do some of these and there are also some naïve implementations, such as vtkRuledSurfaceFilter, but you won't get too far if you need robust processing of real-world contours.

Kyle (cc'd) (@Sunderlandkyl) worked on a conversion class for a few years, refining the algorithm whenever it failed on some new data sets that we received. The class is in vtkSegmentationCore (which only depends on VTK).

There is no side-by-side comparison of the capabilities of the current Plastimatch and SlcierRT. Here's the relevant point from Csaba from that same email thread:

We uploaded all the data that we used for testing initially and that came from users as special cases that failed to load but we fixed since then. Maybe it helps:

https://github.com/SlicerRt/SlicerRtData

Also, this was the response from Greg Sharp (@gregsharp), the man behind plastimatch:

I'm confident that Plastimatch works well for your task of converting rtss to image. It is mature, has lots of command line options, and is easy to use in scripts.

For my own purposes, and whenever pyradiomics users ask me for a recommendation, I always recommend plastimatch, just because I do not have hard evidence about advantages of SlicerRT, and running Slicer python scripts in batch mode is far from pleasant. But my experience with DICOM RT is overall rather limited.

If you have any evidence to recommend one tool over the other, or perhaps recommend dcmrtstruct2nii instead, it would be very interesting to understand how you came up to that conclusion.

Also cc'ing @pieper, since he was also involved in that email thread back in 2017. Hope I didn't forget anyone!

[cpinter]

Plastimatch only provides that data reading, which is the planar contours data. There is a planar contours to closed surface conversion in SlicerRT, which is indeed very complex (I wouldn't call it heuristics; it is a deterministic algorithm). It is needed to process the RTSTRUCT without data loss, because direct planar contours to labelmap conversion lacks basic considerations (i.e. that the real-world object does not contain staircases, or that end-capping is needed).

We created this script in SlicerRT that takes DICOM RTSS and converts it to nrrd https://github.com/SlicerRt/SlicerRT/tree/master/BatchProcessing You can run it from the command-line, so I don't quite understand the batch processing unpleasantness part.

If you are content with the rough conversion, you can choose not to use this and go with Plastimatch.

[pieper]

@shuo0220 I don't do much with RT but have known the people behind SlicerRT and plastimatch for a long time and they are world class experts who have put a lot of thought and effort into these tools. Of course it's best to validate that any of these tools are giving what you expect, and ideally you can report back what you find out.

[fedorov]

We created this script in SlicerRT that takes DICOM RTSS and converts it to nrrd https://github.com/SlicerRt/SlicerRT/tree/master/BatchProcessing You can run it from the command-line, so I don't quite understand the batch processing unpleasantness part.

@cpinter I am aware about the batch converter, and below is my response to essentially similar comment that either you or Andras made back in 2017, which I believe holds today:

In our usage scenario, the only need is conversion of RTSS in the batch mode. Running python scripts with "--run-python-script" always has a startup cost, so it will be one time per dataset, not one time per session; it requires importing everything into the DICOM database, and also requires setting up X. Since we don't need any of the Slicer GUI, I was just hoping that we could find something smaller. But if we know that Plastimatch is inferior to SlicerRT, perhaps these will be the costs to pay.

For developers, costs in using SlicerRT are much higher. I recently talked with someone from a group that is developing certain analysis platform, and they have a need to deal with RT. I mentioned SlicerRT to them, but they absolutely did not want to have Slicer as a dependency. I think overall it would be really nice to have more modularity in SlicerRT, like dealing with all IO aspects in Plastimatch, but I understand there are many priorities.

About staircase effect, that is not a problem for me, since for the use cases I deal with, I typically need to rasterize the contours to the image matrix.

I am not prescribing the specific approach (at least I hope it does not come across that way), I am just sharing my perspective, and trying to help user make an informed decision for their situation.

[lassoan]

I've had a look at https://pypi.org/project/dcmrtstruct2nii/, and it is one of those naive implementations. It does not seem to handle keyholes, cannot intepolate between slices (which may be necessary if you want to reduce loss of accuracy due to discretization of the contours), etc. It may work well for specific workflows, but not a generally usable, robust implementation.

I agree that it used to be more complicated to do batch processing that involved Slicer modules. However, many things have changed in the past few years. Processing using bash/batch/etc. scripts is getting less popular and instead people tend to preprocess data in Python, often using Jupyter notebooks. Since Slicer Preview versions support Python 3 and you can install any Python package in Slicer's Python environment and Slicer can be even used as a Jupyter notebook kernel, I think it is much easier to use Slicer in batch operations. There is no additional startup cost, as the whole Python script or Jupyter notebook kernel runs in the Slicer process. Probably we could further simplify our examples and maybe show how to iterate through all the rtstruct files in a folder in Python.

[cpinter]

This is a good point Andras about the jupyter notebooks. Maybe we should convert that batch processing command line script into a notebook.

[fedorov]

Great idea Andras! I think the batch use case is still important, at least for some users, but I can see how the notebook could address needs of many users. The example you mentioned would definitely be good to have.

[shuo0220]

Thank you @fedorov for cc'ing so many experts. Thank you, everyone, for all the great comments! A notebook for batch processing would be great. I'm working in radiation oncology. I feel like our treatment planning system is also moving toward python/jupyter notebook for scripting purposes. Thank you @lassoan for looking at the dcmrtstruct2nii, and pointing out the difference! Another general question is that, from a standardization standpoint, are there any guidelines/recommendations on how should handle the conversion between DICOM and NRRD/NIFTI. Thank you all again!

[shuo0220]

@fedorov @cpinter @lassoan Sorry to bother you all again. I have another question about the temporary nrrd files (for both label and image) during the feature extraction process of SlicerRadiomics in 3D slicer GUI. As shown in the attachment, I copied out the temporary nrrd files (before it got cleaned up) during during the feature extraction process of SlicerRadiomics in 3D slicer (green box). Then I created the binary label map manually in 3D slicer and saved them to nrrd in 3D slicer in an uncompressed way (oragne box). Then I used Pyradiomics in 3D slicer GUI to extract features on both sets (green and orange boxes), and some of the features have big difference. For example, the relative diff% of original-glcm-Autocorrelation is about 17%, the wavelet-HHL-firstorder-Median and wavelet-HHL-firstorder-Mean are different by 5 to 6 times with opposite sign (-0.004421166 Vs. 0.024718373) and (0.000377281 vs -0.001650913). Also, when I read the label nrrd files (highlighted by red) a package called pynrrd into numpy array, the two label nrrd gave me different matrix shape for this case. I'm wondering which nrrd images and labels should we trust? Also, i was also using this batch-process of RTSS conversion in command line (https://github.com/SlicerRt/SlicerRT/blob/master/BatchProcessing/BatchStructureSetConversion.py) and it also gave me different feature results compared to the temporary nrrd files in 3D slicer GUI. Thank you for your attention and time!

[JoostJM]

@shuo0220, can you share the nrrd files?

[shuo0220]

@JoostJM Here are 2 examples of the label nrrd files. the original dicom files and image nrrd files are bigger than 10MB. I can share via Box. Thank you! example.zip

[JoostJM]

@shuo0220 I ran the following snippet to compare your labelmaps (I had to resample them to compare them properly, as there are small differences in size etc.

import numpy as np
import SimpleITK as sitk

#case = ('1/slicer_saved_label.nrrd', '1/DIBAI_vtkMRMLLabelMapVolumeNodeB.nrrd')
case = ('2/slicer_save_label.nrrd', '2/BECHG_vtkMRMLLabelMapVolumeNodeB.nrrd')

ma1 = sitk.ReadImage(case[0])
ma2 = sitk.ReadImage(case[1])

rif = sitk.ResampleImageFilter()
rif.SetReferenceImage(ma1)
rif.SetInterpolator(sitk.sitkNearestNeighbor)

ma2 = rif.Execute(ma2)

ma1_arr = sitk.GetArrayFromImage(ma1)
ma2_arr = sitk.GetArrayFromImage(ma2)

print(np.array_equal(ma1_arr, ma2_arr))

For both cases, this returned an output of True, i.e. both labels encode exactly the same ROI.

Can you run something similar on your input volumes? Additionally, what happens when you run SlicerRadiomics using 1 input volume, but the 2 different labelmaps as input? or 1 labelmap, but 2 different input volumes?

[cpinter]

Regarding the batch conversion script: it makes a difference whether you tell it to use a reference or not, see https://github.com/SlicerRt/SlicerRT/blob/master/BatchProcessing/BatchStructureSetConversion.py#L280

It would help if you'd provide us some info about the geometries of the labelmaps that you say are different. What are the dimensions of the anatomical image? How are the labelmaps different?

[shuo0220]

Thank you so much for your replies. @JoostJM I have run the code on the input volumes, and it printed True as well. The same input volume with different labels will result different results (as attached), but different volume with same labels will have the same results (features). @cpinter I also applied "-u" function in the BatchProcessing, the result is also included in the attached excel file. But it also seems to be different. I'm new to this, and not sure if all the diagnostics included will provide you the info needed. Thanks again for your time. Example2.xlsx

[shuo0220]

I just did a few more examples. I attached the slicer saved label maps and the temporal built label maps for the 5 examples. If I don't do the resample as @JoostJM suggested, the size of the two label maps are a little different, and it seems like the temp label maps are always bigger than the slicer saved ones. The origin is different a little bit as well. I would assume when SlicerRadiomics extracts the features, it won't apply the resampling process?

Below is the print out of these cases: `Size: slicer_saved vs temp (53, 69, 109) vs (55, 71, 111) Origin: slicer_saved vs temp (207.96875, -247.50000000000003, -151.0) vs (207.10937500000003, -248.359375, -152.0)

Size: slicer_saved vs temp (62, 94, 85) vs (64, 96, 87) Origin: slicer_saved vs temp (194.5703125, -234.41406250000006, -134.0) vs (193.90625, -235.07812500000006, -135.0)

Size: slicer_saved vs temp (107, 113, 178) vs (105, 112, 176) Origin: slicer_saved vs temp (185.625, -343.125, -175.99999999999997) vs (186.5625, -342.1875, -174.99999999999997)

Size: slicer_saved vs temp (91, 125, 196) vs (89, 124, 195) Origin: slicer_saved vs temp (175.89000000000004, -294.294, -164.0) vs (175.89000000000004, -293.43600000000004, -162.99999999999997)

Size: slicer_saved vs temp (80, 63, 177) vs (82, 64, 177) Origin: slicer_saved vs temp (198.75, -337.5, -169.0) vs (198.75, -337.5, -170.0) `

Here are some additional info: i'm using nightly built 3D slicer version 4.11.0-2019-08-13. The temporary label maps I copied out were saved to here: C:\Users...\AppData\Local\Temp\Slicer When I save the nrrd manually, I did not do any more conversion after creating the binary label map from the RTStructure. I don't know if I should do that or not.

Thank you again for all your time! examples.zip

[JoostJM]

@shuo0220, I've taken a look at your Excel file and found the cause of your differences.

If you take a look at your settings, you'll see that you have enabled resampling to [1, 1, 1]. This means that PyRadiomics has to define a new grid of points and calculate values at those points. Although in all cases, your grid will have the same spacing, the alignment will be slightly different. In PyRadiomics the grid is always aligned with the corner of the origin voxel, which is located a different points. This slight difference in alignment, in turn, causes differences in the resampled values, and thus differences in the extracted values.

[shuo0220]

@JoostJM Thank you for pointing out this. The difference was gone when the resampling is disabled. Was the alignment difference because of the small difference in size between the two label maps, like you previously mentioned that you have to resample the label maps to compare them properly?

Also, I feel we have to use resampling in order to have a fair comparison among patients due to different in-plane resolution and slice thickness. If we want to use resampling ([1,1,1] or else), which one is more reliable or accurate in your opinion (SlicerRadiomics GUI w/ temporary nrrd Vs. BatchProcess Saved nrrd files? Thank you again very much! Really appreciate it.

[JoostJM]

Yes, because the origin was different, the corner was different too, causing a mismatch in alignment during resampling.

I fully agree that resampling is a necessary part to deal with voxel anisotropy. Both methods are reliable, though I'd advise to use one method throughout your dataset.

I am looking into the function itself to see if I can deal with the alignment issue.

[shuo0220]

@JoostJM Thank you so much for all your suggestions and help!