pythonimageopencvpydicommedical-imaging

Converting hand radiograph DICOM to PNG returns white/bright image

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)
Solution

  • 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:
    enter image description here

    Output of above sample code:
    enter image description here

    Output of img = lin_stretch_img(img, 0.01, 99.99) (may give better result):
    enter image description here

    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:
    enter image description here