J2K Pixel Data is not decoded correctly when Pixel Representation doesn't match J2K sign

[amirbar]

Describe the bug
DCM pixel_data is not decoded correctly. The resulting pixels are clearly in the wrong range (see image)

Looking at few interesting dicom tags:

1. RescaleIntercept = 0
2. RescaleSlope = 1
3. PixelRepresentation = 1
4. TransferSytntax is JPEG2000Lossless.
5. BitsStored = 13
6. HighBit = 12
7. BitsAllocated = 16

It seems that like based on PixelRepresentation, the data should already be in int16 and the Rescale values imply there shouldn't be any change applied to the input data. However, the resulting values are clearly wrong.

Expected behavior
When loaded using Horos viewer, this is the result (showing default window):

Steps To Reproduce
mydcm.dcm.zip

Unzip the attached dicom and run the following code:

import pydicom
from matplotlib import pyplot as plt

dcm = pydicom.dcmread("mydcm.dcm")
arr = pydicom.pixel_data_handlers.util.apply_voi_lut(dcm.pixel_array, dcm)
plt.imshow(arr, cmap='gray')

or alternatively,

use apply_modality_lut to view the expected HU scale values.

Environment

module	version
platform	Linux-5.3.0-28-generic-x86_64-with-Ubuntu-16.04-xenial
Python	3.6.10 (default, Dec 19 2019, 23:04:32) [GCC 5.4.0 20160609]
pydicom	2.1.0.dev0
gdcm	3.0.4
jpeg_ls	module not found
numpy	1.16.3
PIL	4.3.0

[scaramallion]

Hmm, this is a strange one

• The J2K data is unsigned (Ssiz is 0x0C -> 0b00001100 which from Table A-11 in ISO/IEC 15444-1 means unsigned), but the Pixel Representation is 1. Unsigned CT data with an identity rescale operation is very unusual.
• The CT appears to be a Hitachi Supria, but looking at their conformance it seems unlikely that's where the compression originates
• The Bits Stored value is correct. The range of pixel values runs from 0 to 8191, which matches as well.
• Pillow, GDCM and pylibjpeg all produce the same (apparently incorrect) output, which implies the J2K image data is incorrect or openjpeg is at fault.

I'm a bit curious about the order of operations here, are you using Horos to view the attached dataset or is the dataset already in Horos and exported as the attached dataset? What application performed the compression?

[amirbar]

Very strange one, I agree.

are you using Horos to view the attached dataset or is the dataset already in Horos and exported as the attached dataset?

I imported it to Horos and use Horos only to view it. In addition, I was just able to view it through a cornerstone based web viewer (an internal tool we have). I tried to reverse engineer their open source code to see if there's any discrepancy, but I couldn't find any hint yet.

What application performed the compression?

I am not sure. However, (0002, 0013) Implementation Version Name value is 'OFFIS_DCMTK_354'. Could that be DCMTK? I also tried decompressing it manually using their tools but haven't succeeded.

[scaramallion]

DCMTK is likely the backend that was used to transfer the dataset.

[scaramallion]

Running GDCM's gdcmconv also produces the same (incorrect) output.

[scaramallion]

It seems like Horos is taking the unsigned J2K data then applying the 2s complement correction to make signed data in order to match the Pixel Representation. I guess that would make sense if the J2K Ssiz parameter doesn't match the Pixel Representation in the dataset (such as here), but I'm wondering how generalisable that would be. It's interesting, it didn't occur to me to do that under this situation (I have pylibjpeg-openjpeg printing a warning about the mismatch), but it makes sense.

I guess we could add config flags to the various handlers to do the conversion automatically.

@amirbar is it OK if we add your dataset to our tests?

[scaramallion]

In the meantime you can fix it by converting to 2s complement yourself:

import numpy as np
import pydicom
from pydicom.pixel_data_handlers.util import apply_voi_lut
from matplotlib import pyplot as plt

ds = pydicom.dcmread("mydcm.dcm")
ds.PixelRepresentation = 0 
bit_shift = ds.BitsAllocated - ds.BitsStored
arr = (ds.pixel_array << bit_shift).astype(np.int16) >>  bit_shift
ds.PixelRepresentation = 1
vis = apply_voi_lut(arr, ds)
plt.imshow(vis, cmap='gray')
plt.show()

[amirbar]

That's pretty neat! Sure, you can use it.

[amirbar]

BTW, one assumption in the code above is that ds.BitsAllocated == 16, is this always the case?
Also, @scaramallion, thank you for the quick replies and solutions.

[scaramallion]

BTW, one assumption in the code above is that ds.BitsAllocated == 16, is this always the case?

Oh, good point. This should be more general (for cases where there's a mismatch):

import numpy as np
import pydicom
from pydicom.pixel_data_handlers.util import apply_voi_lut, pixel_dtype
from matplotlib import pyplot as plt

ds = pydicom.dcmread("mydcm.dcm")
ds.PixelRepresentation = 0 
bit_shift = ds.BitsAllocated - ds.BitsStored
arr = ds.pixel_array
ds.PixelRepresentation = 1
dtype = pixel_dtype(ds)
arr = (arr << bit_shift).astype(dtype) >>  bit_shift
vis = apply_voi_lut(arr, ds)
plt.imshow(vis, cmap='gray')
plt.show()

Also, @scaramallion, thank you for the quick replies and solutions.

No problem!