nifti to dicom conversion introduces negative intensities

Hello @zivy,

I tried what you proposed, but when I try to set the slope to 0.001, I get an assertion failed Error from GDCM:

Assertion failed: slope == (int)slope, file C:\toolchain\libs\ITK\5.2.1\source\ITK\Modules\ThirdParty\GDCM\src\gdcm\Source\MediaStorageAndFileFormat\gdcmRescaler.cxx, line 67

Is there an alternative?

Thanks again!

Hi @Keyn34,

Not sure why that isn’t working as you are exercising the same functionality as the SimpleITK example (it’s just one wrapping layer over ITK/GDCM). One thing to clarify, the DICOM dictionary based approach should be applied to the original data, no need to apply the RescaleFilter.

Code snippet doesn’t show how the slope and intercept are set in the dictionary. When looking at the specific line reported in the error, it appears we shouldn’t have reached it as the parameter there is int and there is another method that deals with float. So possibly something in the way the values were set in the dictionary?

1 Like

Hey @zivy and thanks again,

That is what I am wondering about as well, it should be the same as in SimpleITK. And thanks, I was aware of that and just to confirm that I am not using the rescaling together with the dictionary approach.

I originally set the slope and intercept as

itk::EncapsulateMetaData<std::string>(sliceDictionary, DICOMTag::Slope, "0.001");
itk::EncapsulateMetaData<std::string>(sliceDictionary, DICOMTag::Intercept, "0.0");

and this caused the above mentioned Error.

When I set it as

itk::EncapsulateMetaData<float>(sliceDictionary, DICOMTag::Slope, 0.001);
itk::EncapsulateMetaData<float>(sliceDictionary, DICOMTag::Intercept, 0.0);

I don’t get errors, but the slope and intercept are not applied, it seems?

Am I misunderstanding something?
Thanks again!

Hi @Keyn34,

Not sure why the rescale slope and intercept are ignored when set as float.

Before we dig any deeper, can you try and slightly modify the SimpleITK script (read the nifti image instead of the array creation) and confirm that the script does what you want. Once that is done, we can try to work backwards and see where the settings differ.

1 Like

What happens if you set slope to 1000? The meaning might be inverse of what you seem to be trying.

If you do itk.WriteImage<itk::Image<unsigned short, 3>>(rescale->GetOutput(), "debug.nrrd"); or itk.WriteImage<itk::Image<double, 3>>(rescale->GetOutput(), "debug.nrrd");? It is better to visualize in Slicer, as ITK-SNAP converts intensities to short internally.

1 Like

Hello @zivy,

I modified the script to read the nifti 4D volume and removed the reading, because I just needed the output. In case it helps, below the script:

import SimpleITK as sitk
import sys
import time
import os

nifti_path = "path to nifti file"
output_path = "path to output folder"

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

    # 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", "PT")

    # (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, slice, volume)))))
    #   Instance Number
    image_slice.SetMetaData("0020,0013", str(slice))

    # Write to the output directory and add the extension dcm, to force
    # writing in DICOM format.
    image_name = str(volume) + '_' + str(slice) + '.dcm'
    print(image_name)
    writer.SetFileName(os.path.join(out_dir, image_name))

    writer.Execute(image_slice)


if __name__ == '__main__':

    # Read nifti file
    niftiReader = sitk.ImageFileReader()
    niftiReader.SetFileName(nifti_path)
    img_new = niftiReader.Execute()

    number_of_volumes = img_new.GetSize()[3]
    number_of_slices = img_new.GetSize()[2]

    print(str(number_of_slices)+"/"+str(number_of_volumes))

    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 = img_new.GetDirection()
    series_tag_values = [
        ("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
    ]

    rescale_slope = 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'),  # rescale intercept
        ('0028|0100', '16'),  # bits allocated
        ('0028|0101', '16'),  # bits stored
        ('0028|0102', '15'),  # high bit
        ('0028|0103', '1')]  # pixel representation

    # Write slices to output directory
    print(range(number_of_slices))
    print(range(number_of_volumes))

    for v in range(number_of_volumes):
        for s in range(number_of_slices):
            writeSlices(series_tag_values, img_new, output_path, s, v)

    sys.exit(0)

The results are not really… useable I would say, but I can imagine that this is my fault due to potential errors I made when editing the script:

Although, when I look at one slice, I can see that the DICOM tags are written correctly:

Does this give a hint, or is there anything else I can try?

@dzenanz, thank you, but when I set the slope to 1000, nothing changes as it gets ignored again.

What exactly do you mean by that?

The same visual errors occur when I open the images with slicer, unfortunately.

Thank you again!

When you examine the images with Slicer, that eliminates ITK-SNAP’s reduction to 16-bit signed as a potential source of errors.

What I meant by itk.WriteImage is to check how the image looks like after intensity rescaling, but before converting/writing it to DICOM format.

1 Like

Thank you @dzenanz for clearing that up for me!

Good, so it is clear that the images themselves are cropped when it comes to their intensity.
I got it working now with the rescaler as well, which confirms that additionally. The result is shown below:

The speckles are gone, and the image looks good now. The issue I face now is that the intensity values are altered (obviously) and not the same anymore between the two datasets for the same cursor position. There is no way to retain the intensities and do the conversion without the speckle issue, I suppose?

Also, when I do the following:

std::cout << "Input image intensity range is " << rescaler->GetInputMinimum() << " to " << rescaler->GetInputMaximum() << std::endl;

I get as output: Input image intensity range is 3.40282e+38 to 0, which is kinda unexpected.

Thanks!

This will not work properly before rescaler->Update(), which might be called implicitly via writer->Update().

It is probably possible, but you must not use rescaling. Going via setting DICOM tags rescale slope (and rescale intercept) is your best bet. Whoever wrote this page in 2008 thought that writing float pixels with GDCM was impossible. I am not sure whether that has changed.

But looking at the DICOM standard, it should be possible. You might have to use DCMTK instead of GDCM.

1 Like

As far as I know, writing float with GDCM is still not possible. At least for me, it crashes every time I try to set the pixel type to float.

I tried to build ITK with DCMTK, but I get the same errors as in here. Is there any workaround, building DCMTK by itself, or applying the “patch” manually?
Edit: It was possible to first build DCMTK and afterwords set the flag in ITK to use system DCMTK.

Are there any examples on writing DICOM with DCMTK as ImageIO within and Series Writer? Because I tried replacing the GDCM IO by DCMTK IO and now no files get written.

Thanks!

@zivy, sorry to bump you again directly, but I also found this post about the same issue, where the intercept and slope won’t get updated. Maybe this issue wasn’t resolved till now?

Thank you again!

Hi @Keyn34,

SimpleITK/ITK can write float pixels to DICOM via GDCM using float values for slope, intercept. The DicomSeriesFromArray example illustrates this:

python DicomSeriesFromArray.py  output_dir float64

Can you share the original nifti file? If yes, I’ll try and see what settings enable saving it as DICOM with float values, what modifications need to be done to your current SimpleITK script. Once we figure that out you’ll be able to set your ITK code accordingly.

1 Like

@zivy,

Thank you so much. I send you a private message with a download link to access the nifti file.

Hello @Keyn34,

The rescale slope/intercept in the original script were not appropriate for your data. Please try the following script (does not retain the exact values, but nearly, you’ll have to play with this if you want an exact copy):

import SimpleITK as sitk

import sys
import time
import os

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", "PT")

    # (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)


if len(sys.argv) < 3:
    print("Usage: python " + __file__ + " <input_file> <output_directory>")
    sys.exit(1)

new_img = sitk.ReadImage(sys.argv[1])

# Write the 3D image as a series
# IMPORTANT: There are many DICOM tags that need to be updated when you modify
#            an original image. This is a delicate opration and requires
#            knowledge of the DICOM standard. This example only modifies some.
#            For a more complete list of tags that need to be modified see:
#                  http://gdcm.sourceforge.net/wiki/index.php/Writing_DICOM
#            If it is critical for your work to generate valid DICOM files,
#            It is recommended to use David Clunie's Dicom3tools to validate
#            the files:
#                  http://www.dclunie.com/dicom3tools.html

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 = new_img.GetDirection()
series_tag_values = [
    ("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 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.
arr_view = sitk.GetArrayViewFromImage(new_img)
i_min = arr_view.min()
i_max = arr_view.max()

# the bulk intensity data range is [-2^15,+2^15], stored using 2's complement
# we rescale the actual data to fit in this discrete range.
rescale_slope = (i_min-i_max)/65536
rescale_intercept = i_min
series_tag_values = series_tag_values + [
        ('0028|1053', str(rescale_slope)),  # rescale slope
        ('0028|1052', str(rescale_intercept)),  # rescale intercept
        ('0028|0100', '16'),  # bits allocated
        ('0028|0101', '16'),  # bits stored, 
        ('0028|0102', '15'),  # high bit
        ('0028|0103', '1')]  # pixel representation

# Write slices to output directory
list(map(lambda i: writeSlices(series_tag_values, new_img, sys.argv[2], i),
         range(new_img.GetDepth())))
sys.exit(0)
2 Likes

Thank you @zivy, I just tested your script with the slight modification to just take the first volume for a quick run and the results look very good already. Really thank you very much for that - very helpful.

The remaining question would be how to apply the rescale and intercept to the written series within ITK for C++? Because here I still face the issue that the slope and intercept do not get applied.

Hello @Keyn34,

Happy to hear that you’re closer to solving the problem.

So, I don’t think there is much you need to do beyond a correct configuration of the meta-data information. Going over the C++ snippets above I already identified one difference. The C++ code has PixelRepresentation as 0 and the SimpleITK code has it as 1.

My advice is that you run the SimpleITK code and print out the contents of the meta-data dictionary and do the same for the C++ code. Compare the two and you’ll see what you need to set to get the C++ code to work.

2 Likes

Hello @zivy,

Thank you so much again for all your help so far.

When I set PixelRepresentation to 1, which is also carried out to the written .*dcm files, nothing changes. Everything I set in the C++ code is carried over, except the rescale slope and intercept. Below is an image for comparing the header information:

Do the other tags matter in that regard?

Also, what I noticed is the following: I have set the PixelType for my DICOM images like this using OutputImagePixelType = float;. This works fine, but the images are ALWAYS interpreted as short or int when read by ITK-SNAP. Are these expected or are we facing a possible error here, where I just set my datatypes wrongly?

An additional indicator for that (I think) is that when I take another data set, where one time set the datatype to unsigned short (gets displayed not correctly), and one time to signed short ( gets displayed correctly) the Pixel Representation tag is only set rightly for the correctly displayed one:


Left: correctly displayed and correctly set Pixel Representation
Right: not correctly displayed due to being unsigned with intensities starting at -1024 and not correctly set PixelRepresentation tag

Again, thank you a lot. Feels like we are getting close.

Hello @Keyn34,

using OutputImagePixelType = float; is the right setting. What you are experiencing is due to the use of ITK-SNAP. I believe it converts the data during reading. If you open the image in Slicer you’ll see the original float values.

Not sure what is causing the writer to ignore your setting of the slope/intercept values, that is very strange. Will require some more debugging on your side (see that the GDCMImageIO actually accesses these values during writing walk through this portion of the code with your debugger).

1 Like

Thank you @zivy,

Okay, then I will switch to Slicer altogether for checking the results. I just found it weird that I can read the nifti file, but the dicom images are wrapped regarding their intensity - interesting.

Thank you again, I will try walking through the GDCMImageIO part during the writing and report back if I find something.

Anything else there is that I did wrong potentially?

Hello @zivy,

I got it working. I am not quite sure what the tipping point was, but my best bet is that the pixel type of the output image must be float, and setting the metadata must be of type std::string. After I changed that back, it was working flawlessly.

itk::EncapsulateMetaData<std::string>(sliceDictionary, DICOMTag::Slope, std::to_string(rescaleSlope));
itk::EncapsulateMetaData<std::string>(sliceDictionary, DICOMTag::Intercept, std::to_string(rescaleIntercept));

Thank you so much for your support @zivy and @dzenanz as well of course! Great to have this forum with you guys.

2 Likes