Reading Dicom: Rescaling does not allow pixel format change

[PtiLuky]

Description

I have a dicom image stored as short, which has a rescaling into the "unsigned short" domain that cannot be loaded with ITK.

Bits allocated 16
Bits stored 16
High Bit 15
Pixel Representation 1

Rescale Intercept 32768
Rescale Slope 1
Rescale Type US

When itk::GDCMImageIO::ReadImageInformationis called, there is an exception raised "Pixel type larger than output type".

Steps to Reproduce

High level snippet:

    using ReaderType = itk::ImageFileReader<itk::Image<uint16_t, 2>>;
    typename ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName(a_path);
    reader->SetImageIO(itk::GDCMImageIO::New());
    reader->Update();

Minimal code snippet:

auto a_imageIO = itk::GDCMImageIO::New();
a_imageIO->SetFileName(dicomPath.c_str()); // this file has the tags listed above
a_imageIO->ReadImageInformation(); // This raises the exception and refuse to load

This is also reproducible with the code taken directly from void GDCMImageIO::InternalReadImageInformation() https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/IO/GDCM/src/itkGDCMImageIO.cxx#L538

       // snippet with code taken from "void GDCMImageIO::InternalReadImageInformation()"
        gdcm::ImageReader reader;
        reader.SetFileName(dicomPath.c_str());
        reader.Read();
        const gdcm::Image& image = reader.GetImage();
        const gdcm::PixelFormat& pixeltype = image.GetPixelFormat(); // This is "INT16" (value 5)

       // note this is slightly changed from GDCMImageIO::InternalReadImageInformation (where it's conditionally done)
        const float rescaleSlope = 1.;  // This part comes from the dicom tags reader, but is simplified here for ease of reproducibility
        const float rescaleIntercept = 32768.;  // This part comes from the dicom tags reader, but is simplified here for ease of reproducibility
        gdcm::Rescaler r;
        r.SetIntercept(rescaleIntercept );
        r.SetSlope(rescaleSlope );
        r.SetPixelFormat(pixeltype);
        auto outputpt = r.ComputeInterceptSlopePixelType(); // This is "UINT16" (value 4)
        // Note: end of the "changed" part, the following condition is in current ITK master branch

        if (pixeltype > outputpt) { // INT16 > UINT16, despite the fact that after rescaling result will fit in UINT16...
            itkAssertInDebugOrThrowInReleaseMacro("Pixel type larger than output type")
        }

Expected behavior

Load image correctly without any exception. Like in other readers such as VTK or MicroDicom.

Actual behavior

Exception raised "Pixel type larger than output type".
Printed logs:

"<redacted>\\InsightToolkit-5.2.1\\Modules\\IO\\GDCM\\src\\itkGDCMImageIO.cxx:524:\nITK ERROR: Pixel type larger than output type"

Reproducibility

All the time given the problematic input image.

Versions

Currently testing with 5.2.1, but I see the problematic check still done in master branch.

Environment

Windows 10, MSVC

Additional Information

I do not see any reason for this check to be done, it means if my stored data is int32 and I have a rescale/slope that project the result into a int8 range itk raises an error? I see nothing wrong with that...

   if (pixeltype > outputpt) {
        itkAssertInDebugOrThrowInReleaseMacro("Pixel type larger than output type")
    }

Would there be any workaround other than removing the rescaling from the tags, load the image, and then reapply it manually afterward?

[github-actions]

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.

[dzenanz]

What happens if you replace this: using ReaderType = itk::ImageFileReader<itk::Image<uint16_t, 2>>; by this: using ReaderType = itk::ImageFileReader<itk::Image<uint32_t, 2>>;? In other words, change your desired image type from 16-bit to 32-bit?

[PtiLuky]

@dzenanz unfortunately it changes nothing, I tried with a using ReaderType = itk::ImageFileReader<itk::Image<int64_t, 2>>; to be safe and had the same issue.

It seems in itkGDCMImageIO.cxx the output type is overwritten with scaling information, regardless what was provided as a template parameters

[dzenanz]

@malaterre and @issakomi might have some opinion about this.

[issakomi]

I have created a test file and could reproduce the issue exactly as described above.

Not sure what the lines should do

    if (pixeltype > outputpt)
    {
      itkAssertInDebugOrThrowInReleaseMacro("Pixel type larger than output type")
    }

There is no > operator defined in gdcmPixelFormat, there are == and != operators for PixelFormat and PixelFormat::ScalarType. pixeltype is PixelFormat, outputpt is PixelFormat::ScalarType, not sure what is >-compared and how.

The rest seems to be just fine, Rescaler calculated the correct format.

[dzenanz]

Last touched by #1882, and introduced by f00506b.

[dzenanz]

Maybe explicitly handle discrepancy between bits stored and bits allocated, and remove this assert?

[issakomi]

Last touched by #1882, and introduced by f00506b.

In the first version f00506b there was the assert(pixeltype <= outputpt);. There are many asserts in GDCM that look like a kind of reminder and are for debug (#ifndef NDEBUG) only. So probably replacement with itkAssertInDebugOrThrowInReleaseMacro was too strict.
The comparison is probably related to the enum ScalarType ? Not sure about it.

Edit:

Maybe explicitly handle discrepancy between bits stored and bits allocated, and remove this assert?

I also saw the comment in the #1882:

When pixel type is larger than the output type, it is likely
something wrong in the dataset currently being read.

An example is if the following two headers are "out-of-sync"
(0028,0100) BitsAllocated
(0028,0101) BitsStored

I don't understand what "out-of-sync" means and how this comparison can prevent it. BTW, in this particular example the input type is signed short, the output type is unsigned short, this is unusual (but not necessarily a mistake) and not "larger".

[PtiLuky]

@issakomi I know the data set is unusual but it is not a mistake, we had a manufacturer that produced this kind of files for X-rays (that are strictly positive values) (yes, they could have written with a Pixel Representation of 0 (unsigned), but they did not).

From the compiler pov, yes I am fairly certain the > operator ends up comparing the enum implicitly casted into a int (notice how it's an enum and not an enum class, thus allowing implicit conversion to/from int).

If I understand correctly ideally this check shall be replaced with a another check that still verifies the sanity of the data, without blocking this situation?

I can remove this assert locally and change for a test such as following in my local clone of ITK, but I would not feel confident enough to submit a PR of this.

bool ptLarger = false;
switch (outputpt) 
{
  case gdcm::PixelFormat::UINT8:
   ptLarger = pixeltype > gdcm::PixelFormat::INT8;
   break;
  case gdcm::PixelFormat::UINT12:
   ptLarger = pixeltype > gdcm::PixelFormat::INT12;
   break;
  case gdcm::PixelFormat::UINT16:
   ptLarger = pixeltype > gdcm::PixelFormat::INT16;
   break;
  case gdcm::PixelFormat::UINT32:
   ptLarger = pixeltype > gdcm::PixelFormat::INT32;
   break;
  case gdcm::PixelFormat::UINT64:
   ptLarger = pixeltype > gdcm::PixelFormat::INT64;
   break;
  default:
   ptLarger = pixeltype > outputpt;
    break;
}

if (ptLarger )
{
  itkAssertInDebugOrThrowInReleaseMacro("Pixel type larger than output type")
}

Also do you think there can later be issues with buffer size while reading the data later because of ignoring this check?

[issakomi]

IMHO, this exception and assert could have been removed. It was just an assert in debug mode only 15 years ago. From the PR that turned it into an exception in release mode, it is not clear what the problem was and how to reproduce it. Even if there are some datasets that crash GDCMIO at this point (I have never seen such datasets), this is not a proper fix for the bug. The assumption that pixel sizes are being compared is wrong, I don't know what is compared and why the input can not be larger, honestly.
I hope Mathieu Malaterre, the owner of GDCM, will look into this issue, he approved the PR 5 years ago.

[malaterre]

I don't know what is compared and why the input can not be larger, honestly.

The original assert() was a naive check to verify we would never turned 16bits stored image into 8bits after rescale operation. I guess a new check should be implemented using simply bitsAllocated. One should pay attention to the FLOAT32/FLOAT64 case.

[dzenanz]

@PtiLuky do you want to give it a try?

[PtiLuky]

Uhm, my best guess would be the snippet I proposed above, using just the "bitsAllocated" would allow to cast f32 into i32, so testing explicitly the values for UINT8, UINT12, UINT16, UINT32 and UINT64 seems the best if we want to keep.

I don't specially want to be the author of this contribution since I'm not very familiar with the whole contribution workflow of this repo.

But if nobody wants to take it, I can have a look.