Implement `NrrdReader` for `nrrd` and `seg.nrrd` files

[kbressem]

Is your feature request related to a problem? Please describe.
For one project we are working with seg.nrrd files to create pixelwise segmentations. However, when trying to load these files with LoadImaged I get the following error:

TemplateTypeError: itk.Image is not wrapped for input type `itk.Vector[itk.SS,3], int`.

To limit the size of the package, only a limited number of
types are available in ITK Python. To print the supported
types, run the following command in your python environment:

    itk.Image.GetTypes()

Possible solutions:
* If you are an application user:
** Convert your input image into a supported format (see below).
** Contact developer to report the issue.
* If you are an application developer, force input images to be
loaded in a supported pixel type.

    e.g.: instance = itk.Image[itk.RGBPixel[itk.UC], int].New(my_input)

* (Advanced) If you are an application developer, build ITK Python yourself and
turned to `ON` the corresponding CMake option to wrap the pixel type or image
dimension you need. When configuring ITK with CMake, you can set
`ITK_WRAP_${type}` (replace ${type} with appropriate pixel type such as
`double`). If you need to support images with 4 or 5 dimensions, you can add
these dimensions to the list of dimensions in the CMake variable
`ITK_WRAP_IMAGE_DIMS`.

Describe the solution you'd like
When reading seg.nrrd files with pynrrd the files are handled correctly. I would thus propose to implement a NrrdReader based on pynrrd.
I have already written the NrrdReader for private use and would be happy to make a PR.

Describe alternatives you've considered
I have iterated through all supported pixel types, but this either raises an error or leads to incorrect import of the images. For example, some images have format [4, H, W, 1]. When forcing pixel type to itk.Image[itk.RGBPixel[itk.UC], int], the image is read as [3, H, W, 1], thus omitting a layer. This is a problem, as the number of layers in the different seg.nrrd files varies.

[wyli]

could you please share one example segmentation here, so that we can replicate/learn from the use cases?

[kbressem]

Sure. Here it is: example_file.zip

[wyli]

Sure. Here it is: example_file.zip

Thanks, what would be the expected outcome? I'm using

import itk
import matplotlib.pyplot as plt
from monai.transforms import LoadImage

out = LoadImage(reader="itkreader", pixel_type=itk.Vector[itk.F,3])("./example.seg.nrrd")
print(out[0].shape)
plt.imshow(out[0][..., 0])
plt.show()

and the output shape (out[0].shape) is (2500, 2048, 1, 3) and in the meta data dictionary out[1]['original_channel_dim'] is set to -1 (channel last). In my understanding the outputs are correct for this case.

[kbressem]

Yes, this is correct for this case, as it has three channels. The problem with .seg.nrrd is, that the number of channels are inconsistent because, segmentations without overlap are merged to the same channel (or layer) of the file. So if we have four segmentations, and all overlap, the .seg.nrrd would have four channels. If only three segmentations overlap, the .seg.nrrd would have three channels.
I also tried forcing the pixel type to a specific format. But this leads to problems, as soon as the number of channels in the seg.nrrd change. For example,

import itk
import nrrd
from monai.transforms import LoadImage

out_itk = LoadImage(reader="itkreader", pixel_type=itk.Vector[itk.F, 3])("example_3_channel.seg.nrrd")
print(out_itk[0].shape)
>>> (2500, 2048, 1, 3)

out_nrrd = nrrd.read("example_3_channel.seg.nrrd")
print(out_nrrd[0].shape)
>>> (3, 2500, 2048, 1)

out_itk = LoadImage(reader="itkreader", pixel_type=itk.Vector[itk.F, 3])("example_4_channel.seg.nrrd")
print(out_itk [0].shape)
>>> (2500, 2048, 1, 3)

out_nrrd = nrrd.read("example_4_channel.seg.nrrd")
print(out_nrrd[0].shape)
>>> (4, 2500, 2048, 1)

WIth forcing ITK to read the image as itk.Vector[itk.F, 3] one channel in of example_4_channel.seg.nrrd is discarded and the segmentation(s) on in it are lost. I have iterated through all supported pixel types in ITK, none lead to a correct image format.
Here is the file: example_4_channel.seg.zip

[wyli]

I see, I think an NrrdReader might be useful here, please send a pull request...

CC @lassoan @thewtex in case they have better ideas about pixel_type in itk.imread in this case.

[lassoan]

We provide a seg.nrrd reader in slicerio Python package. It can be used in any Python environments. It parses all metadata, including terminology information (content of the segments using standard coded terms), so you can find a segment not by some label value but by standard codes. It can also extract voxels of selected segments as a numpy array, each with a desired label value, regardless of how many layers are in the segmentation and what label values are used internally.

[lassoan]

The main advantage of using standard coded entries for describing what is in a segment because then you don't need to hardcode across your data set that label value 2 means liver and label value 3 means tumor. These label values can only be used in a single data set because these arbitrary label values are chosen differently in each project. In contrast, if you specify that there is a (SCT, 10200004, 'Liver') and a (SCT, 4147007, 'Mass') segment in the segmentation then this information is remains valid forever, in every context. This is the same way how DICOM specifies what is in a segmentation.

DICOM recommends using SNOMED Clinical Terms (and you are granted free license to use all the terms that are part of the DICOM standard). You can find a list of codes that the DCMQI project assembled here:
https://github.com/QIICR/dcmqi/blob/ac7d0fea3a4458a2e14bae0bac886813051b6ad4/doc/segContexts/SegmentationCategoryTypeModifier-SlicerGeneralAnatomy.json#L903-L918

You can find more details in this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5675033/