Pydicom: Wrong jpeg 2000 image from pillow handler
Description
jpeg 200 image does not look right when pydicom uses pillow as the handler as opposed to gdcm
Steps/Code to Reproduce
Couldn't upload dcm file, so I zipped it. Please unzip to continue
CT000008.dcm.zip
import numpy as np
from pydicom import dcmread
import matplotlib.pyplot as plt
plt.interactive(False)
ds = dcmread("CT000008.dcm")
print(ds.file_meta)
print(ds.pixel_array)
v = ds.pixel_array + -1024
v = np.clip(v, 0, 80)
plt.imshow(v, cmap='gray')
plt.colorbar()
plt.title('2000')
plt.savefig('with_pillow.jpeg')
plt.show()
Expected Results
Same code as above, just disable jpeg 2000 from pillow (when gdcm is installed)
import numpy as np
from pydicom import dcmread
from pydicom.config import image_handlers
print(image_handlers)
image_handlers[2].have_pillow_jpeg2000_plugin = False
import matplotlib.pyplot as plt
plt.interactive(False)
ds = dcmread("CT000008.dcm")
print(ds.file_meta)
print(ds.pixel_array)
v = ds.pixel_array + -1024
v = np.clip(v, 0, 80)
plt.imshow(v, cmap='gray')
plt.colorbar()
plt.title('2000')
plt.savefig('without_pillow.jpeg')
plt.show()
Actual Results
Versions
Darwin-17.2.0-x86_64-i386-64bit
('Python', '2.7.13 (default, Dec 18 2016, 07:03:34) \n[GCC 4.2.1 Compatible Apple LLVM 7.0.2 (clang-700.1.81)]')
('numpy', '1.14.0')
('pydicom', '1.0.2')
DavidBord
All 14 comments
In this comment on issue #539 we talk about some issues with jpeg2000/Pillow. The test suite only has jpeg 2000 images from MR modality.
When I uncompress your image with gdcmconv, I get the image you expect.
In pydicom's pillow_handler, there is a clause in the pixel decoder with the comment:
if (dicom_dataset.file_meta.TransferSyntaxUID in
PillowJPEG2000TransferSyntaxes and
dicom_dataset.BitsStored == 16):
# WHY IS THIS EVEN NECESSARY??
pixel_array &= 0x7FFF
This was necessary to get the MR_small_jp2k_lossless.dcm test image to match the MR_small.dcm source image. The comment means I have no idea what I'm doing...
When I run the following tests on your image and a gdcm decompressed version, CT_raw.dcm, I see the following behaviour:
import pydicom
import numpy.testing
gt = pydicom.read_file("CT_raw.dcm")
img = pydicom.read_file("CT000008.dcm")
gt_pixels = gt.pixel_array
img_pixels = img.pixel_array
numpy.testing.assert_array_equal(img_pixels, gt_pixels)
Arrays are not equal
(mismatch 99.8756408691%)
x: array([[24768, 24768, 24768, ..., 24768, 24768, 24768],
[24768, 24768, 24768, ..., 24768, 24768, 24768],
[24768, 24768, 24768, ..., 24768, 24768, 24768],...
y: array([[-2000, -2000, -2000, ..., -2000, -2000, -2000],
[-2000, -2000, -2000, ..., -2000, -2000, -2000],
[-2000, -2000, -2000, ..., -2000, -2000, -2000],...
For some reason:
In [22]: gt.BitsStored
Out[22]: 14
In [23]: img.BitsStored
Out[23]: 16
numpy.testing.assert_array_equal(img_pixels/4, gt_pixels)
Arrays are not equal
(mismatch 21.2753295898%)
x: array([[6192, 6192, 6192, ..., 6192, 6192, 6192],
[6192, 6192, 6192, ..., 6192, 6192, 6192],
[6192, 6192, 6192, ..., 6192, 6192, 6192],...
y: array([[-2000, -2000, -2000, ..., -2000, -2000, -2000],
[-2000, -2000, -2000, ..., -2000, -2000, -2000],
[-2000, -2000, -2000, ..., -2000, -2000, -2000],...
In [31]: max(gt_pixels[gt_pixels != img_pixels/4])
Out[31]: -2000
In [32]: min(gt_pixels[gt_pixels != img_pixels/4])
Out[32]: -2000
So, when the jpeg2000 compressed image is divided by 4, then the two arrays only disagree in the pixels with the PixelPaddingValue.
For some reason, something in the original image told gdcm to shift 2 bits and store only a 14 bit result. Nothing in the image tags indicates to me why that should be, but it is obviously there.
Maybe someone more familiar with CT / DICOM might know that answer?
rhaxton
on 30 Jul 2018
For reference the Transfer Syntax UID is 1.2.840.10008.1.2.4.90 - JPEG 2000 Image Compression (Lossless Only). One frame in one fragment, basic offset table has a value. Frame length is 105362 bytes compressed starting from offset 1690 to offset 107051.
Partial JPEG header (SOC + SIZ)
ff 4f ff 51 00 29 00 00
00 00 02 00 00 00 02 00
00 00 00 00 00 00 00 00
00 00 02 00 00 00 02 00
00 00 00 00 00 00 00 00
00 01 8d 01 01
0x8D is the SIZ marker's Ssiz component (1000 1101). From IEC 15444-1:2000:
Precision (depth) in bits and sign of the ith component samples. The precision is the precision of the
component samples before DC level shifting is performed (i.e., the precision of the original component
samples before any processing is performed). If the component sample values are signed, then the range of component sample values is -2^[(Ssiz AND 0x7F)-1] ≤ component sample value ≤ 2^[(Ssiz AND 0x7F)-1]- 1
And for binary values starting with 1:
1xxx xxxx Component sample values are signed values
(0x8D AND 0x7F) is 13, not sure where the other bit goes...
scaramallion
on 31 Jul 2018
Oh
x000 0000 — x010 0101 - Component sample bit depth = value + 1.
So bit depth is 14 then, which matches gdcm
scaramallion
on 31 Jul 2018
This might be a way to approach it
def jpeg2k_bit_depth(bytestream):
# Add appropriate bytestream checks to ensure jpeg2k
ssiz = bytestream[42:43]
if ssiz[0] & (1 << 7):
# Signed
bit_depth = (ssiz[0] & b'\x7f'[0]) + 1
else:
# Unsigned
bit_depth = ssiz[0] + 1
return bit_depth
import matplotlib.pyplot as plt
import numpy as np
from pydicom import dcmread
ds = dcmread('path/to/file')
arr = ds.pixel_array
# Shift 2 bits based difference 16 -> 14-bit as returned by jpeg2k_bit_depth
padded_pixels = np.where(arr & (1 << 14))
arr = np.right_shift(arr, 2)
arr = arr * ds.RescaleSlope + ds.RescaleIntercept
arr[padded_pixels] = ds.PixelPaddingValue
v = np.clip(arr, 0, 80)
plt.imshow(v, cmap='gray')
plt.colorbar()
plt.title('2000')
plt.show()
It looks reasonable, but I don't have access to GDCM at the moment so no idea how it compares. In particular the step where the overflow pixels are used to shorthand where the padding is required seems a bit dodgy.
scaramallion
on 31 Jul 2018
This is definitely related to signed images (PixelRepresentation=1). If I start with the emri_small.dcm test image and I change the PixelRepresentation to 1 and then gdcmconv to jp2k, then the image will decode improperly.
The weird thing is that when I gdcmconv --j2k -i emri_small_signed.dcm -o emri_small_j2k_signed.dcm, the jpeg header indicates signed 16 bit values, but the gdcmconv --raw -i emri_small_j2k_signed.dcm -o emri_small_gdcm_raw.dcm does not change the BitsStored from 12 to 16 like it did with the CT image you posted....
With a signed image, masking out the high bit results in a correct image for any positive pixel values, but negative values are wrong.
rhaxton
on 2 Aug 2018
I don't think this bug is specific to 14 bit images. As long as the bits allocated is 16 and the bits stored is < 16, PIL handler will return the values such that the whitest white possible at the stored bit depth is now 2^16 -1 when decompressing the image from jpeg2k, essentially upscaling whatever was stored into bit stored 16. Without changing the tags. Which makes your viewer display it as burned white, as from its view the entire image has been bit shifted.
MinkBlaxx
on 9 Aug 2018
@DavidBord is it OK if we include the dataset you supplied as one of our testing datasets?
scaramallion
on 12 Oct 2018
@scaramallion yes
DavidBord
on 14 May 2019
The 14-bit 512x512 JPEG2000 file from the sample dataset: 693.jpg
scaramallion
on 20 Dec 2019
Interestingly I get the correct value if I flip the MSB of each SV (and then shift >> 2 to account for the different bit depth) rather than zeroing it. The unit tests are passing (but maybe that says more about our pillow unit tests than anything else).
if (transfer_syntax in PillowJPEG2000TransferSyntaxes and
ds.BitsStored == 16):
# WHY IS THIS EVEN NECESSARY??
# Zero MSb: b01111111 11111111
#arr &= 0x7FFF
# Flip MSb: b10000000 00000000
arr ^= 0x8000
And then doing
arr = ds.pixel_array
arr = numpy.right_shift(arr, 2)
Flipping the MSB fixes #539 as well (with 16-bit signed JPEG2000.dcm).
I wish we had way more test files for JPEG2k...
~Still have no idea why zeroing is required, let alone flipping~
OK, by XORing and shifting >> we are just undoing the conversion Pillow does when converting from N-bit signed data to 16-bit unsigned (as @MinkBlaxx has said).
scaramallion
on 21 Dec 2019
Have you checked the test images at ftp://medical.nema.org/MEDICAL/Dicom/DataSets/WG04?
Though these are quite old, and not many...
Edit: well, not so bad: 40 j2k-encoded images with different values for BitsStored, PixelRepresentation and PhotometricInterpretation.
mrbean-bremen
on 21 Dec 2019
Bit shifted arrays use:
arr = np.right_shift(arr, ds.BitsAllocated - ds.BitsStored)
Reversible J2K (1.2.840.10008.1.2.4.90) equal to uncompressed?
Filename Bits Equal w/ flip Equal w/ zero
emri_small_jpeg_2k_lossless.dcm u16 Yes* Yes
MR_small_jp2klossless.dcm i16 Yes No
CT2_J2KR i16 Yes No
MR2_J2KR u12 Yes, with shift Yes, with shift
RG1_J2KR u15 Yes, with shift Yes, with shift
RG3_J2KR u10 Yes, with shift Yes, with shift
US1_J2KR u8 Yes Yes
* After setting ds.BitsStored to 16 to match data (was 12)
Irreversible J2K (1.2.840.10008.1.2.4.91) equal to uncompressed?
Filename Bits Equal w/ flip Equal w/ zero
CT000008.dcm i14 Yes, with shift** No
JPEG2000.dcm i16 Yes No
SC_rgb_gdcm_KY.dcm u8 Yes Yes
CT2_J2KI i16 Yes No
MR2_J2KI u12 Yes, with shift Yes, with shift
RG1_J2KI u15 Yes, with shift Yes, with shift
RG3_J2KI u10 Yes, with shift Yes, with shift
US1_J2KI* u8 Yes Yes
* unable to decompress, viewed image looks same as GDCM
** After setting ds.BitsStored to 14 to match data (was 16)
I think that's conclusive enough to justify a PR with the changes.
scaramallion
on 22 Dec 2019
Related issues
tenggyut
·
3Comments
razorx89
·
8Comments
oliverzzm
·
12Comments
amandagb
·
8Comments
vsoch
·
9Comments
Most helpful comment
I don't think this bug is specific to 14 bit images. As long as the bits allocated is 16 and the bits stored is < 16, PIL handler will return the values such that the whitest white possible at the stored bit depth is now 2^16 -1 when decompressing the image from jpeg2k, essentially upscaling whatever was stored into bit stored 16. Without changing the tags. Which makes your viewer display it as burned white, as from its view the entire image has been bit shifted.