Error: origin is outside image geometry

[yigitsoy]

Hi,

I have been trying to load DICOM SEG object to Slicer (Slicer-4.13), with the DCMQI extension. The GUI reports that it cannot be loaded. Then, when I try to convert the DICOM SEG using the following command

/Applications/Slicer.app/Contents/Extensions-29731/DCMQI/lib/Slicer-4.13/cli-modules/segimage2itkimage --inputDICOM 55ba63c8-67da807c-b085d75d-8d4d12ae-cdaa94ed_dicom_seg.dcm -t nifti --outputDirectory ./temp_dcm/

I get an error about the image origin, as the following

Row direction: 1 0 0
Col direction: 0 0.995396 -0.0958457
Z direction: -0 0.0958457 0.995396
Total frames: 10
Total frames with unique IPP: 10
Total overlapping frames: 0
Origin: [-118, -29.0505, -684.721]
Failed to get CIELab values - initializing to default 43803,26565,37722
Failed to get CIELab values - initializing to default 43803,26565,37722
Failed to get CIELab values - initializing to default 43803,26565,37722
Failed to get CIELab values - initializing to default 43803,26565,37722
ERROR: Frame 3 origin [-118, -29.0505, -682] is outside image geometry![0, -1, 3]
Image size: [512, 512, 10]
Fatal error encountered.
vtkDebugLeaks has found no leaks

Here is a basic script how I use highdicom to create the DICOM SEG object

import nibabel
import numpy as np
from pydicom.sr.codedict import codes
from pydicom.filereader import dcmread
from pydicom.uid import generate_uid

from highdicom.content import AlgorithmIdentificationSequence
from highdicom.seg.content import SegmentDescription
from highdicom.seg.enum import (
    SegmentAlgorithmTypeValues,
    SegmentationTypeValues
)
from highdicom.seg.sop import Segmentation

from pathlib import Path

# Path to directory containing single-frame legacy CT Image instances
series_dir = Path('./55ba63c8-67da807c-b085d75d-8d4d12ae-cdaa94ed_dcm/')
image_files = series_dir.glob('*')

# Read CT Image data sets from PS3.10 files on disk
image_datasets = [dcmread(str(f)) for f in image_files]


# # Create a binary segmentation mask
mask = np.zeros(
    shape=(
        len(image_datasets),
        image_datasets[0].Rows,
        image_datasets[0].Columns
    ),
    dtype=np.uint16
)
mask[10:20, 10:-10, 100:-100] = 1


# Describe the algorithm that created the segmentation
algorithm_identification = AlgorithmIdentificationSequence(
    name='test',
    version='v1.0',
    family=codes.cid7162.ArtificialIntelligence
)

# Describe the segment
description_segment_1 = SegmentDescription(
    segment_number=1,
    segment_label='first segment',
    segmented_property_category=codes.cid7150.Tissue,
    segmented_property_type=codes.cid7166.ConnectiveTissue,
    algorithm_type=SegmentAlgorithmTypeValues.AUTOMATIC,
    algorithm_identification=algorithm_identification,
    tracking_uid=generate_uid(),
    tracking_id='test segmentation of computed tomography image'
)

# Create the Segmentation instance
seg_dataset = Segmentation(
    source_images=image_datasets,
    pixel_array=mask,
    segmentation_type=SegmentationTypeValues.BINARY,
    segment_descriptions=[description_segment_1],
    series_instance_uid=generate_uid(),
    series_number=2,
    sop_instance_uid=generate_uid(),
    instance_number=1,
    manufacturer='Manufacturer',
    manufacturer_model_name='Model',
    software_versions='v1',
    device_serial_number='Device XYZ',
)

print(seg_dataset)

seg_dataset.save_as("./seg.dcm")

It would be great if you could provide some pointers as to where the problem could be. Thank you.

[yigitsoy]

Note that I do not get this error when I use a different series as the input.

[yigitsoy]

Not being 100% sure, from the second dimension value of the mapped index [0, -1, 3], I suspect this to be some sort of rounding issue due to the parsing of the frame origin from string here: https://github.com/QIICR/dcmqi/blob/master/libsrc/ImageSEGConverter.cpp#L735

[lorteddie]

Can you try setting the locale to C before calling the DCMQI code?

[yigitsoy]

Tried it, no success.

[fedorov]

@yigitsoy I tried the sample code you shared on a public dataset, and I am not able to reproduce the issue you reported. Can you either share a de-identified dataset, or (ideally) reproduce the issue on a public dataset?

[yigitsoy]

@fedorov thanks for trying it out. Unfortunately I cannot share the dataset. Furthermore, it is hard to reproduce on another dataset. I tried on different series we have and it worked just fine. I will try to reproduce this on a public dataset and share it.

[d-vogel]

I'm hitting the same issue and started from this issue.
I'm trying to import a DICOMDIR export from Brainlab Elements (they call it BrainDump IIRC).

Not being 100% sure, from the second dimension value of the mapped index [0, -1, 3], I suspect this to be some sort of rounding issue due to the parsing of the frame origin from string here: https://github.com/QIICR/dcmqi/blob/master/libsrc/ImageSEGConverter.cpp#L735

this is now here

Indeed there is something going on with the rounding or even more. I've applied the following changes:

--- a/libsrc/Dicom2ItkConverter.cpp
+++ b/libsrc/Dicom2ItkConverter.cpp
@@ -139,14 +139,17 @@ itk::SmartPointer<ShortImageType> Dicom2ItkConverter::nextResult()
                     m_groupIterator = m_segmentGroups.end();
                     return nullptr;
                 }
-                if (!itkImage->TransformPhysicalPointToIndex(frameOriginPoint, frameOriginIndex))
-                {
-                    cerr << "ERROR: Frame " << framesForSegment[frameIndex] << " origin " << frameOriginPoint
-                         << " is outside image geometry!" << frameOriginIndex << endl;
-                    cerr << "Image size: " << itkImage->GetBufferedRegion().GetSize() << endl;
-                    m_groupIterator = m_segmentGroups.end();
-                    return nullptr;
-                }
+      
+                slicer_itk::ContinuousIndex<double, 3> frameOrigin_cidx = itkImage->TransformPhysicalPointToContinuousIndex<double>(frameOriginPoint);
+                cerr << "Frame "       << framesForSegment[frameIndex] << 
+                        " Origin "     << frameOriginPoint << 
+                        " | Index: "   << frameOrigin_cidx << 
+                        " | Spacing: " << m_computedSliceSpacing << endl;
+                // if (!itkImage->TransformPhysicalPointToIndex(frameOriginPoint, frameOriginIndex))
+                // {
+                // m_groupIterator = m_segmentGroups.end();
+                    // return nullptr;
+                // }
                 // Handling differs depending on whether the segmentation is binary or fractional
                 // (we have to unpack binary frames before copying them into the ITK image)
                 const DcmIODTypes::Frame* rawFrame      = m_segDoc->getFrame(framesForSegment[frameIndex]);

I have a lot of DICOMSeg, but here is an example output:

######################################################################################
$EXPORTDATA/12111870/46863212/04716225
dcmqi repository URL: [email protected]:d-vogel/QuantitativeReporting.git revision: 6908e49 tag: v1.3.0-1-g6908e49
Row direction: 0.999999 0.000690647 -0.000792147
Col direction: -0.000793156 0.00146193 -0.999999
Z direction: -0.000689488 0.999999 0.00146248
Total frames: 9
Total frames with unique IPP: 9
Total overlapping frames: 0
Origin: [-136.178, -24.4365, 161.783]
Slice extent: 16.3628
Slice spacing: 2.16572
WARNING: SliceThickness is present and is 2. using it!
Will not merge segments: Splitting segments into 1 groups
Frame 0 Origin [-136.221, -8.07363, 161.729] | Index: [-0.119559, 0.297448, 8.18141] | Spacing: 2.16572
Frame 1 Origin [-135.969, -10.0313, 162.041] | Index: [0.831287, -0.899452, 7.2027] | Spacing: 2.16572
Frame 2 Origin [-136.155, -12.2658, 161.767] | Index: [0.118362, 0.128887, 6.08534] | Spacing: 2.16572
Frame 3 Origin [-135.997, -14.2671, 161.983] | Index: [0.714826, -0.699767, 5.08481] | Spacing: 2.16572
Frame 4 Origin [-136.18, -16.4469, 161.743] | Index: [0.0135769, 0.197681, 3.99478] | Spacing: 2.16572
Frame 5 Origin [-135.9, -18.2456, 162.016] | Index: [1.07144, -0.849658, 3.09553] | Spacing: 2.16572
Frame 6 Origin [-136.194, -20.4117, 161.792] | Index: [-0.049583, -0.00887286, 2.0124] | Spacing: 2.16572
Frame 7 Origin [-135.598, -22.271, 162.172] | Index: [2.20389, -1.46353, 1.08286] | Spacing: 2.16572
Frame 8 Origin [-136.178, -24.4365, 161.783] | Index: [0, 0, 0] | Spacing: 2.16572
Writing itk image to $OUTPUT/1.nrrd ... done

For all objects, the origin of the first or last frame is perfectly fine (index [0,0,0]) but then is diverges. It seems the direction matrix in the individual segmentation set (if there's even any) does not match with the info used to init the ITK image. However it's hard to make any sense of the missmatch, as the error is really not linear.

For now I will disable the check and compare the data obtained with DCMSeg to data exported with OpenIGTLink.

[fedorov]

It looks like the jitter is on the order of 0.5 mm, so maybe this is the issue of rounding? But if this is the case, it is a problem with the data. If individual frames have inconsistent origin, then I don't know what tool can make sense out of this to reconstruct a volume. @d-vogel any chance you can bring this up with the BrainLab folks and ask for an explanation?