Converting hand radiograph DICOM to PNG returns white/bright image

Question

I am converting hand X-rays in DICOM format to PNG format. The code below does this:

import os
import cv2
import pydicom
import numpy as np
from PIL import Image

inputdir = "P:/BoneDataset/DICOM-File/0-RefinedDICOM/"
outdir = 'P:/BoneDataset/DICOM-File/1-ConvertedPics/'

test_list = [f for f in os.listdir(inputdir)]

for f in test_list[:10]:
    ds = pydicom.read_file(inputdir + f) # read dicom image
    img = ds.pixel_array # get image array
    scaled_img = (np.maximum(img,0) / img.max()) * 255.0
    img = scaled_img.astype(np.uint8)
    cv2.imwrite(outdir + f.replace('.dcm','.png'),img)

The image below shows some of the result (Input (DICOM) --> Output (PNG)):

Input 1 --> Output 1

Input 2 --> Output 2

Input 3 --> Output 3

As you can see, I want the output images to look the same as the input X-ray, however, I get black and white output. Is this because of image threshold issue or something wrong with the file?

EDIT 1:

I tried the code suggested below, however it returns the same output as before, not as the same image as the input.

I've anonymized the DICOM files, hence you can find the DICOM dataset files HERE. Please use these DICOM files.

EDIT 2:

Trying the code edited update's code below works 100%. Also, using cv2.bitwise_not() also works after the line ds= img.pixel_array. Like this (Trying this works but the suggested answer works too):

test_list = [f for f in os.listdir(inputdir)]

for f in test_list[:10]:
    ds = pydicom.read_file(inputdir + f) # read dicom image
    img = ds.pixel_array # get image array
    img = cv2.bitwise_not(img)
    img = lin_stretch_img(img, 1, 99) # Apply "linear stretching"
    cv2.imwrite(outdir + f.replace('.dcm','.png'),img)

Hard to say, whether there is something wrong with the file, when you do not attach the file... Apart from that, you scale to 0..255 (i.e. uint8) and after that set top uint16. Why? You set negative values of the dicom image to 0. Why? — MSpiller, Oct 11 '22 at 12:15
The reason I cannot attach the file is due to ethical reasons for my research. So I edited it to uint8 however, it returns the same outputs. — Bathtub, Oct 11 '22 at 13:56
I think your issue might be that you are only scale the max. You need to also scale the min. Try using skimage.rescale_intensity. Use the 'input' option to get the input range automatically or adjust the input range to get the desired result. You should also remove the labels so that their brightness does not affect the max value. — fmw42, Oct 11 '22 at 16:12

Rotem · Accepted Answer · 2022-10-12T22:15:34.627

As commented, scaling the maximum may not be enough.
We may try "linear stretching", where some low percentile goes to 0, high percentile goes to 255, and levels in between are transformed linearly.

Additional option:
Since many DICOM images have black margins, and white labels, we may want to ignore the minimum and maximum values, when computing the percentiles.

There is no guarantee that the result is going to be the same as the "the input X-ray", but it is probably not going to be black and white.

For testing, I downloaded CASE 1 DICOM samples from here.

Code sample:

import cv2
import pydicom
import numpy as np

def lin_stretch_img(img, low_prc, high_prc, do_ignore_minmax=True):
    """ 
    Apply linear "stretch" - low_prc percentile goes to 0, 
    and high_prc percentile goes to 255.
    The result is clipped to [0, 255] and converted to np.uint8

    Additional feature:
    When computing high and low percentiles, ignore the minimum and maximum intensities (assumed to be outliers).
    """
    # For ignoring the outliers, replace them with the median value
    if do_ignore_minmax:
        tmp_img = img.copy()
        med = np.median(img)  # Compute median
        tmp_img[img == img.min()] = med
        tmp_img[img == img.max()] = med
    else:
        tmp_img = img

    lo, hi = np.percentile(tmp_img, (low_prc, high_prc))  # Example: 1% - Low percentile, 99% - High percentile

    if lo == hi:
        return np.full(img.shape, 128, np.uint8)  # Protection: return gray image if lo = hi.

    stretch_img = (img.astype(float) - lo) * (255/(hi-lo))  # Linear stretch: lo goes to 0, hi to 255.
    stretch_img = stretch_img.clip(0, 255).astype(np.uint8)  # Clip range to [0, 255] and convert to uint8
    return stretch_img


# https://www.visus.com/fileadmin/content/pictures/Downloads/JiveX_DICOME_Viewer/case1.zip
ds = pydicom.read_file('case1_008.dcm') # read dicom image
img = ds.pixel_array # get image array

img = lin_stretch_img(img, 1, 99)  # Apply "linear stretching" (lower percentile 1 goes to 0, and percentile 99 to 255).

cv2.imwrite('case1_008.png', img)

Output of your code:

Output of above sample code:

Output of img = lin_stretch_img(img, 0.01, 99.99) (may give better result):

Update

The polarity of the sample DICOM images is inverted.
The minimum value is intended to be displayed as white, and the maximum as black.

For correcting the polarity, we may execute img = 255-img (after converting to uint8).

Checking if the polarity is inverted:
According to the documentation, if Photometric Interpretation equals 'MONOCHROME1', then the polarity is inverted ('MONOCHROME2' is not inverted).

MONOCHROME1
Pixel data represent a single monochrome image plane. The minimum sample
value is intended to be displayed as white after any VOI gray scale transformations have been performed. See PS3.4. This value may be used only when Samples per Pixel (0028,0002) has a value of 1. May be used for pixel data in a Native (uncompressed) or Encapsulated (compressed) format.

Inverting polarity if Photometric Interpretation is 'MONOCHROME1':

if ds[0x0028, 0x0004].value == 'MONOCHROME1':
    img = 255-img

The documentation also says we have to apply it "after VOI gray scale transformations".

Applying "VOI gray scale transformations" is described here:

img = apply_voi_lut(img, ds, index=0)

Updated code sample:

import cv2
import pydicom
from pydicom.pixel_data_handlers.util import apply_voi_lut
import numpy as np

def lin_stretch_img(img, low_prc, high_prc, do_ignore_minmax=True):
    """ 
    Apply linear "stretch" - low_prc percentile goes to 0, 
    and high_prc percentile goes to 255.
    The result is clipped to [0, 255] and converted to np.uint8

    Additional feature:
    When computing high and low percentiles, ignore the minimum and maximum intensities (assumed to be outliers).
    """
    # For ignoring the outliers, replace them with the median value
    if do_ignore_minmax:
        tmp_img = img.copy()
        med = np.median(img)  # Compute median
        tmp_img[img == img.min()] = med
        tmp_img[img == img.max()] = med
    else:
        tmp_img = img

    lo, hi = np.percentile(tmp_img, (low_prc, high_prc))  # Example: 1% - Low percentile, 99% - High percentile

    if lo == hi:
        return np.full(img.shape, 128, np.uint8)  # Protection: return gray image if lo = hi.

    stretch_img = (img.astype(float) - lo) * (255/(hi-lo))  # Linear stretch: lo goes to 0, hi to 255.
    stretch_img = stretch_img.clip(0, 255).astype(np.uint8)  # Clip range to [0, 255] and convert to uint8
    return stretch_img


# https://www.visus.com/fileadmin/content/pictures/Downloads/JiveX_DICOME_Viewer/case1.zip
ds = pydicom.read_file('1.2.392.200036.9125.9.0.152034855.3288075520.2287343482.dcm') # read dicom image
img = ds.pixel_array # get image array

# https://pydicom.github.io/pydicom/stable/old/working_with_pixel_data.html#voi-lut-or-windowing-operation
# Apply "VOI gray scale transformations":
img = apply_voi_lut(img, ds, index=0)

img = lin_stretch_img(img, 0.1, 99.9)  # Apply "linear stretching" (lower percentile 0.1 goes to 0, and percentile 99.9 to 255).

# https://dicom.innolitics.com/ciods/rt-dose/image-pixel/00280004
if ds[0x0028, 0x0004].value == 'MONOCHROME1':
    img = 255-img  # Invert polarity if Photometric Interpretation is 'MONOCHROME1'

cv2.imwrite('1.2.392.200036.9125.9.0.152034855.3288075520.2287343482.png', img)

Output sample:

Thank you so much for the code. I tried on the CASE 1 DICOM samples and it worked like a charm. However, when I tried on my hand radiograph, it still returns the same output (white and black) as mentioned in the post. I've anonymized the DICOM file, so please have a look at the edited post above! — Bathtub, Oct 12 '22 at 10:15
Updated the post with the link to my anonymized DICOM Dataset. Please try it if you can? — Bathtub, Oct 12 '22 at 11:37
It looks like the image polarity is inverted. Add `img = 65535 - img` after `img = ds.pixel_array`. Right now, I don't have time for researching the cause. — Rotem, Oct 12 '22 at 13:24
I updated the answer: We have to invert the polarity if Photometric Interpretation value is 'MONOCHROME1' (the update is below the 3 sample images). — Rotem, Oct 12 '22 at 22:17
I did some searching and found that `cv2.bitwise_not()` also works? So, putting that after `img = ds.pixel_array`. Is this acceptable or inverting the polarity more acceptable way? — Bathtub, Oct 12 '22 at 22:59
Also thank you so much for your help, the above works perfectly! — Bathtub, Oct 12 '22 at 22:59
`cv2.bitwise_not(x)` is equivalent to `255 - x` when `x` is `uint8`. `cv2.bitwise_not(x)` is equivalent to `65535 - x` when `x` is `uint16`. Using `cv2.bitwise_not(x)` is probably better solution, because we don't have to check the type (`uint8` or `uint16`...). Place it after `img = apply_voi_lut(img, ds, index=0)`. — Rotem, Oct 14 '22 at 08:05

Converting hand radiograph DICOM to PNG returns white/bright image

1 Answers1

Update

Linked