Skip to content

Writing DICOM Series with Float Data - Rescale issue #2136

Open
@marcmorcos

Description

Hello, this seems very similar to the issue posted here: #1014

I notice an issue where I simply read a dicom series (CT image data, Float64) and then re-write it without changes to test the functionality.

If I explicitly set the rescale and intercept values to the original CT dicom values (-1000.035 and 0.09397), then the output series is correct - data ranges from ~-1000 to +2055.

However, if I set the rescale value to 0.01 and the intercept to 0 (per the documentation: https://simpleitk.readthedocs.io/en/master/link_DicomSeriesFromArray_docs.html) then the output series data ranges from zero to some irrelevant value depending on how many significant digits the rescale value is.

Am I missing something? I was under the impression that ITK would internally do the scaling but this seems not to be the case.
Thanks,
Marc
(Code below for reference)

def writeSeries(IMAGEobj):
    out_dir = "output_DICOM_"+time.strftime("%Y%m%d%H%M%S")
    os.mkdir(out_dir)

    #pixel_dtypes = {"int16": np.int16, "float64": np.float64}
    pixel_dtype = sitk.GetArrayFromImage(IMAGEobj).dtype # np.int16 #np.float64

    # Copy relevant tags from the original meta-data dictionary (private tags are
    # also accessible).
    tags_to_copy = [
        "0010|0010",  # Patient Name
        "0010|0020",  # Patient ID
        "0010|0030",  # Patient Birth Date
        "0020|000d",  # Study Instance UID, for machine consumption
        "0020|0010",  # Study ID, for human consumption
        "0008|0020",  # Study Date
        "0008|0030",  # Study Time
        "0008|0050",  # Accession Number
        "0008|0060",  # Modality
        "0018|5100",  # Patient Position (Marc)
    ]


    def writeSlices(series_tag_values, new_img, out_dir, i):
        image_slice = new_img[:, :, i]

        # Tags shared by the series.
        list(
            map(
                lambda tag_value: image_slice.SetMetaData(tag_value[0], tag_value[1]),
                series_tag_values,
            )
        )

        # Slice specific tags.
        #   Instance Creation Date
        image_slice.SetMetaData("0008|0012", time.strftime("%Y%m%d"))
        #   Instance Creation Time
        image_slice.SetMetaData("0008|0013", time.strftime("%H%M%S"))

        # Setting the type to CT so that the slice location is preserved and
        # the thickness is carried over.
        image_slice.SetMetaData("0008|0060", "CT")

        # (0020, 0032) image position patient determines the 3D spacing between
        # slices.
        #   Image Position (Patient)
        image_slice.SetMetaData(
            "0020|0032",
            "\\".join(map(str, new_img.TransformIndexToPhysicalPoint((0, 0, i)))),
        )
        #   Instance Number
        image_slice.SetMetaData("0020|0013", str(i))

        # Write to the output directory and add the extension dcm, to force
        # writing in DICOM format.
        writer.SetFileName(os.path.join(out_dir, str(i) + ".dcm"))
        writer.Execute(image_slice)

    ####

    writer = sitk.ImageFileWriter()
    # Use the study/series/frame of reference information given in the meta-data
    # dictionary and not the automatically generated information from the file IO
    writer.KeepOriginalImageUIDOn()

    modification_time = time.strftime("%H%M%S")
    modification_date = time.strftime("%Y%m%d")

    # Copy some of the tags and add the relevant tags indicating the change.
    # For the series instance UID (0020|000e), each of the components is a number,
    # cannot start with zero, and separated by a '.' We create a unique series ID
    # using the date and time. Tags of interest:
    direction = IMAGEobj.GetDirection()
    series_tag_values = [
        (k, moving_reader.GetMetaData(0, k))
        for k in tags_to_copy
        if moving_reader.HasMetaDataKey(0, k)
    ] + [
        ("0008|0031", modification_time),  # Series Time
        ("0008|0021", modification_date),  # Series Date
        ("0008|0008", "DERIVED\\SECONDARY"),  # Image Type
        (
            "0020|000e",
            "1.2.826.0.1.3680043.2.1125." + modification_date + ".1" + modification_time,
        ),  # Series Instance UID
        (
            "0020|0037",
            "\\".join(
                map(
                    str,
                    (
                        direction[0],
                        direction[3],
                        direction[6],
                        direction[1],
                        direction[4],
                        direction[7],
                    ),
                )
            ),
        ),  # Image Orientation
        # (Patient)
        ("0008|103e", "Created-SimpleITK"),  # Series Description
    ]

    if pixel_dtype == np.float64:
        print("Processing as float64")
        # If we want to write floating point values, we need to use the rescale
        # slope, "0028|1053", to select the number of digits we want to keep. We
        # also need to specify additional pixel storage and representation
        # information.
        rescale_slope = 0.1# 0.93977#0.001  # keep three digits after the decimal point
        series_tag_values = series_tag_values + [
            ("0028|1053", str(rescale_slope)),  # rescale slope
            ("0028|1052", "0"),#-1000.0358"),  # rescale intercept
            ("0028|0100", "16"),  # bits allocated
            ("0028|0101", "16"),  # bits stored
            ("0028|0102", "15"),  # high bit
            ("0028|0103", "0"),   # pixel representation
        ]  

    # Write slices to output directory
    list(
        map(
            lambda i: writeSlices(series_tag_values, IMAGEobj, out_dir, i),
            range(IMAGEobj.GetDepth()),
        )
    );
    print("...Done!")

def sitkDCM(dir, meta=False):
    print("Reading Dicom directory:", dir)
    reader = sitk.ImageSeriesReader()

    dicom_names = reader.GetGDCMSeriesFileNames(dir)
    
    reader.SetFileNames(dicom_names)
    
    reader.MetaDataDictionaryArrayUpdateOn() #MM
    reader.LoadPrivateTagsOn()               #MM

    image = reader.Execute()
    size = image.GetSize()
    print("Image size:", size[0], size[1], size[2])
    
    return image, reader



moving_image, moving_reader = sitkDCM(dir= "zzCase10_LFrontal")
print( sitk.GetArrayFromImage(moving_image).dtype ) ## PRINTS "INT16"
writeSeries(moving_image)

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions

      Writing DICOM Series with Float Data - Rescale issue · Issue #2136 · SimpleITK/SimpleITK