PixelSpacing incorrectly taken into account

VincentMorard commented 1 year ago

Hi all,

Thanks for all your work, rt-utils is a great tool !

I am working with a MR image and sometimes, depending of the aquisitiion we need to apply a rotation so that the images are "axialized". The result is that the pixelSpacing has a different value for x and y and the rtstruct generated is incorrect.

Here is a python script that allows to create a SEG and a RT Please set force=True in rt_utils when loading the dcm with pydicom, I had to create dcm file to have a self content script

import glob
import os
import pydicom_seg
import pydicom
import SimpleITK as sitk
import numpy as np
from rt_utils import RTStructBuilder
import json

# To update based on the path
dcm_dir = "C:\\Temp\\dicom_file_new"
seg_filename = "C:\\Temp\\segmentation.dcm"
rt_filename = "C:\\Temp\\segmentationRT.dcm"
template_filename = "C:\\Temp\\template.json"
isotrop = False

def process_image(image_data):
    out = np.where(image_data == 180, 1, 0).astype(np.uint16)
    return out

def create_input_dicom(dicom_dir, isotrop):
    os.makedirs(dicom_dir, exist_ok=True)

    header = '{"00080005": {"Value": ["ISO_IR 192"], "vr": "CS"}, "00080008": {"vr": "CS"}, "00080012": {"Value": ["20230721"], "vr": "DA"}, "00080013": {"Value": ["035751"], "vr": "TM"}, "00080014": {"Value": ["3.0.0"], "vr": "UI"}, "00080016": {"Value": ["1.2.840.10008.5.1.4.1.1.4"], "vr": "UI"}, "00080018": {"Value": ["2.25.115231125129161360422423700315339076734"], "vr": "UI"}, "00080020": {"Value": ["20150516"], "vr": "DA"}, "00080021": {"Value": ["20230721"], "vr": "DA"}, "00080022": {"Value": ["20150516"], "vr": "DA"}, "00080023": {"Value": ["20230721"], "vr": "DA"}, "00080030": {"vr": "TM"}, "00080031": {"Value": ["035751"], "vr": "TM"}, "00080032": {"vr": "TM"}, "00080033": {"Value": ["084601"], "vr": "TM"}, "00080050": {"Value": ["15060734050"], "vr": "SH"}, "00080060": {"Value": ["MR"], "vr": "CS"}, "00080070": {"Value": ["UNKNOWN"], "vr": "LO"}, "00080080": {"vr": "LO"}, "00080081": {"vr": "ST"}, "00080090": {"vr": "PN"}, "00081010": {"vr": "SH"}, "00081030": {"Value": ["PREPROCESS"], "vr": "LO"}, "0008103E": {"Value": ["PREPROCESS"], "vr": "LO"}, "00081040": {"vr": "LO"}, "00081048": {"vr": "PN"}, "00081050": {"vr": "PN"}, "00081060": {"vr": "PN"}, "00081070": {"vr": "PN"}, "00081090": {"Value": ["UNKNOWN"], "vr": "LO"}, "00081110": {"Value": [{"00081150": {"Value": ["1.2.840.10008.3.1.2.3.1"], "vr": "UI"}, "00081155": {"Value": ["1.2.124.113532.10.160.226.57.20150312.142717.6605503"], "vr": "UI"}}], "vr": "SQ"}, "00081111": {"Value": [{"00081150": {"Value": ["1.2.840.10008.3.1.2.3.3"], "vr": "UI"}, "00081155": {"Value": ["1.2.840.113619.6.322.240598304242724185143032673684777965063"], "vr": "UI"}}], "vr": "SQ"}, "00081120": {"Value": [{"00081150": {"Value": ["1.2.840.10008.3.1.2.1.1"], "vr": "UI"}, "00081155": {"Value": ["1.2.124.113532.10.160.160.59.20100828.175520.3010744"], "vr": "UI"}}], "vr": "SQ"}, "00100010": {"Value": [{"Alphabetic": "ANNON"}], "vr": "PN"}, "00100020": {"Value": ["ANNON"], "vr": "LO"}, "00100021": {"vr": "LO"}, "00100030": {"Value": ["20230101"], "vr": "DA"}, "00100032": {"vr": "TM"}, "00100040": {"vr": "CS"}, "00101010": {"vr": "AS"}, "00101030": {"vr": "DS"}, "00102110": {"vr": "LO"}, "00102160": {"vr": "SH"}, "00102180": {"vr": "SH"}, "001021B0": {"vr": "LT"}, "00104000": {"vr": "LT"}, "00120062": {"vr": "CS"}, "00180010": {"Value": ["15    g"], "vr": "LO"}, "00180020": {"Value": ["SE"], "vr": "CS"}, "00180021": {"Value": ["SK"], "vr": "CS"}, "00180022": {"vr": "CS"}, "00180023": {"Value": ["3D"], "vr": "CS"}, "00180025": {"Value": ["N"], "vr": "CS"}, "00180050": {"Value": [0.48828125], "vr": "DS"}, "00180080": {"Value": [600.0], "vr": "DS"}, "00180081": {"Value": [12.319], "vr": "DS"}, "00180083": {"Value": [1.0], "vr": "DS"}, "00180084": {"Value": [127.765418], "vr": "DS"}, "00180085": {"Value": ["1H"], "vr": "SH"}, "00180086": {"Value": [1], "vr": "IS"}, "00180087": {"Value": [3.0], "vr": "DS"}, "00180091": {"Value": [25], "vr": "IS"}, "00180093": {"Value": [100.0], "vr": "DS"}, "00180094": {"Value": [90.0], "vr": "DS"}, "00180095": {"Value": [244.141], "vr": "DS"}, "00181000": {"vr": "LO"}, "00181020": {"Value": ["3.0.0"], "vr": "LO"}, "00181030": {"vr": "LO"}, "00181040": {"Value": ["IV"], "vr": "LO"}, "00181088": {"Value": [0], "vr": "IS"}, "00181090": {"Value": [0], "vr": "IS"}, "00181094": {"Value": [0], "vr": "IS"}, "00181100": {"Value": [250.0], "vr": "DS"}, "00181250": {"vr": "SH"}, "00181310": {"Value": [0, 256, 256, 0], "vr": "US"}, "00181314": {"Value": [90.0], "vr": "DS"}, "00181315": {"Value": ["N"], "vr": "CS"}, "00181316": {"Value": [1.57316], "vr": "DS"}, "00185100": {"Value": ["HFS"], "vr": "CS"}, "0020000D": {"Value": ["2.25.255069905384861154924841019137563439343"], "vr": "UI"}, "0020000E": {"Value": ["2.25.205287721952789789408662186088058372087"], "vr": "UI"}, "00200010": {"vr": "SH"}, "00200011": {"Value": [350], "vr": "IS"}, "00200012": {"Value": [1], "vr": "IS"}, "00200013": {"Value": [1], "vr": "IS"}, "00200032": {"Value": [-92.620185852051, -157.00784301758, -132.19627380371], "vr": "DS"}, "00200037": {"Value": [1.0, 0.0, 0.0, 0.0, 1.0, 0.0], "vr": "DS"}, "00200052": {"Value": ["2.25.115231125129161360422423700315339076734"], "vr": "UI"}, "00200060": {"vr": "CS"}, "00201002": {"Value": [320], "vr": "IS"}, "00201040": {"vr": "LO"}, "00201041": {"Value": [-132.19627], "vr": "DS"}, "00280002": {"Value": [1], "vr": "US"}, "00280004": {"Value": ["MONOCHROME2"], "vr": "CS"}, "00280010": {"Value": [512], "vr": "US"}, "00280011": {"Value": [192], "vr": "US"}, "00280030": {"Value": [1.0, 2.0], "vr": "DS"}, "00280100": {"Value": [16], "vr": "US"}, "00280101": {"Value": [16], "vr": "US"}, "00280102": {"Value": [15], "vr": "US"}, "00280103": {"Value": [1], "vr": "US"}, "00280106": {"Value": [0], "vr": "SS"}, "00280107": {"Value": [739], "vr": "SS"}, "00281050": {"Value": [148.0], "vr": "DS"}, "00281051": {"Value": [11.0], "vr": "DS"}, "00281052": {"Value": [0.0], "vr": "DS"}, "00281053": {"Value": [1.0], "vr": "DS"}, "00281054": {"Value": ["US"], "vr": "LO"}, "00282110": {"Value": ["00"], "vr": "CS"}, "00321033": {"vr": "LO"}, "00380500": {"vr": "LO"}, "00400244": {"Value": ["20150516"], "vr": "DA"}, "00400245": {"Value": ["082225"], "vr": "TM"}, "00400253": {"Value": ["7246.1431760945"], "vr": "SH"}, "00400254": {"vr": "LO"}}'
    resZ = 0.48828125

    for i in range(512):
        ds = pydicom.dataset.Dataset.from_json(header)
        ds[0x0020, 0x0032].value[2] = ds[0x0020, 0x0032].value[2] + i * resZ
        ds.SOPInstanceUID = pydicom.uid.generate_uid(prefix=None)

        arr = np.zeros([512, 192]).astype(np.int16)
        arr[100:400, 30:150] = 150
        if i < 256:
            arr[50:200, 60:90] = 180
        else :
            arr[50:200, 60:90] = 130
        ds.PixelData = arr.tobytes()
        ds.is_little_endian = True
        ds.is_implicit_VR = True

        if isotrop:
            ds.PixelSpacing = [1, 1]
        ds.save_as(os.path.join(dicom_dir, f'MR{i:05}.dcm'))

#####################################
# Create the input dicom image
#####################################
create_input_dicom(dcm_dir, isotrop)

#####################################
# Read the original dicom image
#####################################
reader = sitk.ImageSeriesReader()
dcm_files = reader.GetGDCMSeriesFileNames(dcm_dir)
reader.SetFileNames(dcm_files)
image = reader.Execute()

#####################################
# Create dicom SEG object
#####################################
segmentation_data = process_image(sitk.GetArrayFromImage(image))
segmentation = sitk.GetImageFromArray(segmentation_data)
segmentation.CopyInformation(image)

# Write the template to file
template_json = '{"ContentCreatorName": "Reader1", "ClinicalTrialSeriesID": "Session1", "ClinicalTrialTimePointID": "1", "SeriesDescription": "Segmentation", "SeriesNumber": "300", "InstanceNumber": "1", "BodyPartExamined": "", "segmentAttributes": [[{"labelID": 1,    "SegmentDescription": "test",    "SegmentAlgorithmType": "SEMIAUTOMATIC",    "SegmentAlgorithmName": "Test",    "SegmentedPropertyCategoryCodeSequence": {"CodeValue": "123037004",     "CodingSchemeDesignator": "SCT",     "CodeMeaning": "Anatomical Structure"},    "SegmentedPropertyTypeCodeSequence": {"CodeValue": "91750005",     "CodingSchemeDesignator": "SCT",     "CodeMeaning": "1st Diagonal Coronary Artery"}}]], "ContentLabel": "SEGMENTATION", "ContentDescription": "Image segmentation", "ClinicalTrialCoordinatingCenterName": "dcmqi"}'
with open(template_filename, "w") as write_file:
    json.dump(json.loads(template_json), write_file)

template = pydicom_seg.template.from_dcmqi_metainfo(template_filename)
writer = pydicom_seg.MultiClassWriter(
    template=template,
    skip_empty_slices=False,  # Don't encode slices with only zeros
    skip_missing_segment=False,  # If a segment definition is missing in the
                                 # template, then raise an error instead of
                                 # skipping it.
)

dcm_files = glob.glob(dcm_dir + '\\*')
source_images = [pydicom.dcmread(x, stop_before_pixels=True, force=True) for x in dcm_files]
dcm = writer.write(segmentation, source_images)
dcm.save_as(seg_filename)

#####################################
# Compute the RT struct
#####################################
# Swap axis from Z, Y, X to X, Y, Z
segmentation_data = np.swapaxes(segmentation_data, 0, 2)
segmentation_data = np.swapaxes(segmentation_data, 0, 1)

roi_mask = np.where(segmentation_data == 1, True, False)
rtstruct = RTStructBuilder.create_new(dicom_series_path=dcm_dir)
rtstruct.add_roi(mask=roi_mask, name='seg test', color='cc0000')
rtstruct.save(rt_filename)

When I set isotrop=True -> the x, y resolution is the same and the RT struct creation is fine. with MITK:

When I set isotrop=False -> the x, y resolution is not the same and the RT is not at the correct place with MITK:

VincentMorard commented 11 months ago

Changing the line 154 and 184 of the image_helper to : column_spacing, row_spacing = first_slice.PixelSpacing fixed my problem

plesqui commented 11 months ago

Hello @VincentMorard,

Thanks for opening this issue and for the detailed description.

I think your issue could be fixed if the 'synthetic' DICOM that you create would store correctly the PixelSpacing of the swapped axes. That should remove the need to have to modify the source-code within the image_helper module, wouldn't it?

Please let me know if doing something like that would also work.

Regarding the force=True of pydicom:

On one hand, I think it is safer have it as False in the tool, as we ensure dicom files are compliant and therefore we know what to expect within the tool, maximizing the chance of it working as intended.

On the other hand, I acknowledge that often in research we had to do some hacks to get newer stuff to work. So I also understand the need for more flexibility.

This is what I propose: I'll work on a way to implement an optional flag to allow the user to import non-compliant dicom files, with the warning that the tool might not work as intended if used like that. Then I'll bring that into the next release.

If you would like to work on that solution and create a PR, that would be fantastic as well and greatly appreciated!

Thank you for using our tool,

Pedro

qurit / rt-utils

PixelSpacing incorrectly taken into account #92