Writing DICOM Series with Float Data - Rescale issue

[marcmorcos]

Hello, this seems very similar to the issue posted here: https://github.com/SimpleITK/SimpleITK/issues/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 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)

[github-actions[bot]]

Thank you for contributing an issue! 🙏

Welcome to the ITK community! 🤗👋☀️

We are glad you are here and appreciate your contribution. Please keep in mind our community participation guidelines. 📜 Also, please check existing open issues and consider discussion on the ITK Discourse. 📖

This is an automatic message. Allow for time for the ITK community to be able to read the issue and comment on it.

[marcmorcos]

Sorry meant to post in SimpleITK.