How does Pyradiomics deal with several structures segmented in a nrrd file? [FEAT EXTRACTION]

[indranilmallick]

I am trying pyradiomics on a CT scan with DICOM data and a DICOM RTSS structure file. The RTSS file has several structures in it. I have successfully used plastimatch via command line to convert the dicom image files to a single 3d nrrd file and then converted the RTSS dicom file to another nrrd file.
When I perform a feature extraction with pydicom, I get some results, but it is a single set of numbers. I was expecting that I would get values for each structure in the nrrd structure set image. Does pydicom work for only one structure per mask?

Version (please complete the following information):

• OS: [WIndows 10]
• Python version: [e.g. 3.6]
• PyRadiomics version [e.g. 2.2.0]

[JoostJM]

For each case, PyRadiomics only computes the features for 1 ROI in your mask file (identified by the value of label).

See also this example of plastimatch. Multiple ROIs are assigned to the ith bit in the output volume, allowing for overlapping structures. Additionally, plastimatch can output a file showing your segments.

Be aware that PyRadiomics will not function correctly if you have overlapping segments! PyRadiomics only looks at the value as a whole, not a specific index. Overlapping structures in a mask are allowed by supplying a 4D volume as a mask. If there is no overlap it should work fine (i.e. possible labels are 1, 2, 4, 8, etc.)

If you know which labels are present in your mask, you can define them when using PyRadiomics in batch mode (column "Label").

[indranilmallick]

Thanks a lot @JoostJM
I understand the issue a lot better now.
What is still a bit confusing for me is finding out which labels represent which structure.
I had used a RTSS dicom file with 32 structures in it. The file contents are pasted below:

0|0 255 0|BODY
1|0 0 160|Vein_IVC
2|0 255 64|Vein_PV
3|255 0 0|HepVeinRight
4|255 0 0|HepVeinMid
5|255 0 0|HepVeinLeft
6|255 0 0|PortVeinLeft
7|255 0 0|PortVeinRight
8|10 10 250|LivSegI
9|128 255 255|LivSegVII
10|255 128 128|LivSegVIII
11|128 64 64|LivSegVI
12|0 150 0|LivSegV
13|0 150 0|LivSegII
14|0 0 160|LivSegIII
15|255 255 128|LivSegIVa
16|255 0 255|LivSegIVb
17|0 0 255|Vein_SMV
18|0 0 255|Vein_Splenic
19|255 255 0|Colon
20|255 0 255|BowelSmall
21|0 0 255|Stomach
22|187 255 187|Duodenum
23|128 64 64|Pancreas
24|0 255 255|GallBladder
25|0 127 255|Spleen
26|255 128 128|KidneyR
27|50 204 153|KidneyL
28|255 255 0|UreterR
29|255 255 0|UreterL
30|255 210 0|AdrenalR
31|255 210 0|AdrenalL
32|255 0 0|Art_Aorta
33|255 128 0|Art_Celiac
34|255 128 0|Art_SMA
35|255 128 0|Art_IMA
36|255 128 0|Art_Renal
37|0 0 255|Vein_Renal

However, the labels don't match the numbers here.
I am not sure how to figure out how to match the labels with the structures.

[JoostJM]

Can you share the converted mask file? Then I think I can figure out what's what.

[indranilmallick]

This is the link to the .nrrd file.
https://www.dropbox.com/s/mdavn9gxk4fpt5d/abdomen_rtss1.nrrd?dl=0
Will be great if this can be figured out.

[JoostJM]

The output file is indeed more or less like specified in the documentation of plastimatch. However, the datatype is uint8 (meaning 8 bits = 8 possible labels). To allow for more segments (38 in your case), the volume is transformed into a 4D volume, so 8 labels per channel.

To get something that Slicer and PyRadiomics understands, we need to do some restructuring:

import numpy as np
import SimpleITK as sitk

ma = sitk.ReadImage(r'abdomen_rtss1.nrrd')
ma_arr = sitk.GetArrayFromImage(ma)

labels = []
for channel in range(ma_arr.shape[3]):
    for bit in range(8):
        label_map = np.zeros(ma_arr.shape[:3], dtype='uint8')
        roi = np.bitwise_and(ma_arr[..., channel], 2 ** bit) > 0
        if np.sum(roi) > 0:
            label_map[roi] = 1
            labels.append(label_map)

mask_arr = np.array(labels).transpose(1, 2, 3, 0)
mask_arr.shape

mask = sitk.GetImageFromArray(mask_arr)
mask.CopyInformation(ma)

sitk.WriteImage(mask, r'segmentation.seg.nrrd', True)  # Specify True to enable compression

This piece of code iterates over all possible labels and if it finds voxels that match, stores them in a labelmap (be aware that this assumes a 3D volume, with uint8 datatype).

Then at the end, it makes 1 4D matrix and builds an image. As a final step, the geometric information of the mask is copied and the segmentation is stored as a segmentation. This you can load into slicer to see if it is correct.

The first value in the output file (0-37) should match the index in the segmentation. If you want to extract features for, say, "Art_IMA", run PyRadiomics with setting label_channel=35 (note that setting label=1 should not be changed).

Finally, you can also use the SlicerRT extension to load the RTSTRUCT directly into a slicer node and use it to make the segmentation files (which will then also automatically incorporate the name).

[indranilmallick]

Thanks a ton! @JoostJM
Works beautifully. Checked on Slicer and python output.
You have saved me hours of work, making new exports from the planning system. Really appreciate this.

[JoostJM]

Happy to help

[JE-Jimenez]

How would you incorporate that code into the batch file example "batchprocessing.py"?

[JoostJM]

@JE-Jimenez, the advised entrypoint for batchprocessing in PyRadiomics is via the command line (better error handling, multiprocessing support, etc). In that case you can specify a combination of image, mask and label in the csv file using the columns "Image", "Mask" and "Label"/"Label_channel", respectively ("Label" indicates a label value in a labelmap, "Label_channel" the 0-based index of the segment, see 3D Slicer on labelmaps and segmentations for the differences).

[JoostJM]

In batch processing.py, you can achieve something similar, though currently it is only provided for label. See L81-84
for details on parsing out the label value.