Aligning Segmentations and images using python

Support
1

Dear all
I have a problem to understand how to use Segmentations nrrd file with image files in Python.

My problem is:
I’m able to open a stack of images and its segmentations using pynrrd.
However I can’t understand how to align the segmentation with the stack of images.

Reading the header of the segmentation I can’t understand how to shift the segmentation to align with the original image.
The same operation when performed in 3DSlicer is straightforward, and the alignment perfect.
So, I suspect that there is something in the header of the file that I need to understand.

Moreover I’m sure that there is a scale factor between segmentation and images.
Instead, the segmentation is too big (shape (43, 27, 5)) respect to original image ((96, 96, 19)). You have to think that we are in a case of a prostate cancer.
It seems that the segmentation is about 2 times the right dimension.

Can you give some advices on how to solve this problem?
Please remember that I’m not an expert user.

Thank you very much for your help

Andrea

This is an example of the code used and the results obtained on real images:

# Code for images:

filename = ‘image.nrrd’
readdata, header = nrrd.read(filename)
print(readdata.shape)
print(header)

(96, 96, 19)
OrderedDict([(‘type’, ‘unsigned short’), (‘dimension’, 3), (‘space’, ‘left-posterior-superior’), (‘sizes’, array([96, 96, 19])), (‘space directions’, array([[ 2.08301068, 0.03477655, -0.01159939],
[-0.03635903, 2.044063 , -0.40095127],
[ 0.00810046, 0.69308656, 3.53264395]])), (‘kinds’, [‘domain’, ‘domain’, ‘domain’]), (‘endian’, ‘little’), (‘encoding’, ‘gzip’), (‘space origin’, array([-96.8901825 , -92.38396454, -23.10687065]))])

# Code for segmentation

filename_seg = ‘image.seg.nrrd’

readdata_seg, header_seg = nrrd.read(filename_seg)

print(readdata_seg.shape)
print(header_seg)

(43, 27, 5)
OrderedDict([(‘type’, ‘unsigned char’), (‘dimension’, 3), (‘space’, ‘right-anterior-superior’), (‘sizes’, array([43, 27, 5])), (‘space directions’, array([[-6.24903227e-01, -1.04329815e-02, -3.47982208e-03],
[ 1.09077264e-02, -6.13218923e-01, -1.20285387e-01],
[-8.10047081e-03, -6.93086206e-01, 3.53264208e+00]])), (‘kinds’, [‘domain’, ‘domain’, ‘domain’]), (‘encoding’, ‘gzip’), (‘space origin’, array([ 26.12298371, -18.2445476 , -4.61884032])), (‘measurement frame’, array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])), (‘Segment0_Color’, ‘0.501961 0.682353 0.501961’), (‘Segment0_ColorAutoGenerated’, ‘1’), (‘Segment0_Extent’, ‘0 42 0 26 0 4’), (‘Segment0_ID’, ‘Segment_1’), (‘Segment0_Name’, ‘Segment_T2_ADC II’), (‘Segment0_NameAutoGenerated’, ‘0’), (‘Segment0_Tags’, ‘TerminologyEntry:Segmentation category and type - 3D Slicer General Anatomy list~SRT^T-D0050^Tissue~SRT^T-D0050^Tissue~^^~Anatomic codes - DICOM master list~^^~^^|’), (‘Segmentation_ContainedRepresentationNames’, ‘Binary labelmap|’), (‘Segmentation_ConversionParameters’, ‘Compute surface normals|1|Compute surface normals. 1 (default) = surface normals are computed. 0 = surface normals are not computed (slightly faster but produces less smooth surface display).&Crop to reference image geometry|0|Crop the model to the extent of reference geometry. 0 (default) = created labelmap will contain the entire model. 1 = created labelmap extent will be within reference image extent.&Decimation factor|0.0|Desired reduction in the total number of polygons. Range: 0.0 (no decimation) to 1.0 (as much simplification as possible). Value of 0.8 typically reduces data set size by 80% without losing too much details.&Default slice thickness|0.0|Default thickness for contours if slice spacing cannot be calculated.&Fractional labelmap oversampling factor|1|Determines the oversampling of the reference image geometry. All segments are oversampled with the same value (value of 1 means no oversampling).&Oversampling factor|1|Determines the oversampling of the reference image geometry. If it’s a number, then all segments are oversampled with the same value (value of 1 means no oversampling). If it has the value “A”, then automatic oversampling is calculated.&Reference image geometry|-0.624903227496551;0.01090772636925;-0.00810047080951418;96.9144821166992;-0.0104329814853039;-0.613218922740892;-0.693086206443774;94.4632263183594;-0.00347982208181081;-0.1202853868845;3.53264207871727;-33.7047958374023;0;0;0;1;0;319;0;319;0;24;|Image geometry description string determining the geometry of the labelmap that is created in course of conversion. Can be copied from a volume, using the button.&Smoothing factor|0.5|Smoothing factor. Range: 0.0 (no smoothing) to 1.0 (strong smoothing).&Threshold fraction|0.5|Determines the threshold that the closed surface is created at as a fractional value between 0 and 1.&’), (‘Segmentation_MasterRepresentation’, ‘Binary labelmap’), (‘Segmentation_ReferenceImageExtentOffset’, ‘116 166 14’)])
2

Full specification of segmentation storage is available here: http://apidocs.slicer.org/master/classvtkMRMLSegmentationStorageNode.html

If you do not have overlapping segments then you may find it easier to export the segmentation to a single merged labelmap volume and process that volume.

4

Hi! I also have now this problem - labelmap volume is good for non-overlapping, but what one should do when the segments are overlapping?
My CT volume has [512, 512, 326] dimensions and .seg.nrrd has [(29, 342, 241, 83)] dimensions - how to match them?

5

In general, all image processing should happen in physical space and so it should not matter what voxel coordinate system is used. Mapping between physical (LPS) and voxel (IJK) coordinate systems is specified in the image header, by origin, spacing, axis directions (see some more explanation here).

If your software only performs trivial pixel-wise processing (in IJK space), such as masking, and you find that your input volumes have exact same geometry (you must always check this) then you may ignore the physical coordinate system. In this special case, you may find it useful that in recent Slicer preview releases, you can disable cropping segmentation to minimum necessary size in Save Data window, which will result in segmentations that have the same size as the input volume.

6

If I use this feature then my segmentations are same size as volume so they are 3d, not 4d, which means the overlapping segments will not be shown properly

7

If segments do not overlap then they are saved as a 3D volume. If there are overlapping segments then Slicer automatically splits them to non-overlapping layers and saves them as a 4D volume.

@Sunderlandkyl I’ve just noticed that the new “Layer” and “LabelValue” fields are not described in the segmentation file format documentation. Could you add them? (and maybe also add a note on how storage node collapses layers on save)

1 Like
8

Are there some specific settings in Slicer so that overlapping segments do not ‘eat’ each other? Im using the segment editor

9

You can choose overlap/overwrite in masking section after you select an effect.

10

Please, correct me if I’m wrong:
Slicer will by default in older versions crop segmentation to minimum necessary size. But can this information be retrieved from the header in seg.nrrd file? Could not find anywhere clear information about how to much this volumes.

Solution with re-saving could work if not the amount of data that needs to be re-saved.

Thanks

11

Yes, Slicer-4.10.2 always crops to minimum size. Cropping does not affect segmentation’s physical position: if we crop one voxel from we adjust the image origin as needed.

Fortunately, re-saving is really simple. Iterating through any number of scene files, and load and re-save each is probably shold not be more than 4-5 lines of Python code. You can google for code example to iterate through files in a folder hierarchy and find code for reading and writing scene file in the script repository (read, write).

1 Like
12

Thank you very much for your help and explanation, the trick with the script over all files worked like a magic! For anyone who will also encounter this, here are two lines that you need to loop over each file

loadedSegmentation = slicer.util.loadSegmentation(nrrd_file)
savedSegmantation = saveNode(loadedSegmentation, new_nrrd_file)
2 Likes