How to access RGB pixel arrays from DICOM files using pydicom?

Question

I try to access a DICOM file's RGB pixel array with unknown compression (maybe none). Extracting grayscale pixel arrays works completely fine.

However, using

import dicom
import numpy as np

data_set = dicom.read_file(path)
pixel_array = data_set.pixel_array
size_of_array = pixel_array.shape

if len(size_of_array ) == 3:     
    chanR = pixel_array[0][0:size_of_array[1], 0:size_of_array[2]]
    chanG = pixel_array[1][0:size_of_array[1], 0:size_of_array[2]]
    chanB = pixel_array[2][0:size_of_array[1], 0:size_of_array[2]]
    output_array = (0.299 ** chanR) + (0.587 ** chanG) + (0.114 ** chanB)

with the goal to convert it to an common grayscale array. Unfortunately the result array output_array is not containing correct pixel data. Contents are not false scaled, they are spatially disturbed. Where is the issue?

Rand

Maybe its BGR and not RGB? – Circumvent

Nope, the corresponding DICOM tag says "RGB", modality is "OT" for other, its an Patient's report converted to an image format. Resolution and image size of the result image does fit well. But each channel seems to consist of only partial the spatial information. – Rand

What is pixel_array.shape? – Kemp

I ask because other libraries that I've used to read RGB images (not dicom, but JPG, PNG, etc) generally return an array with shape (m, n, 3), not (3, m, n). If that is the case, then you would write chanR = pixel_array[:,:,0]. – Kemp

the corresponding DICOM tag says "RGB" If you are talking about Photometric Interpretation, that does not mean order. – Neurath

Refer dicom.nema.org/medical/dicom/2014c/output/chtml/part03/… and dicomiseasy.blogspot.in/2012/08/chapter-12-pixel-data.html – Neurath

Thank you, @A_J the Photometric Interpretation is in my case RGB, which does not mean order. I'll check your references and will then reply. – Rand

@WarrenWeckesser in my case, the shape is (3, m, n). – Rand

Answer 1

It is not RGB pixel array and the better way is converting to gray image.

The way to get CT Image is to get the attribute of pixel_array in CT dicom file. The type of elements in pixel_array of CT dicom file are all uint16.But a lot of tool in python, like OpenCV, Some AI stuff, cannot be compatible with the type.

After getting pixel_array (CT Image) from CT dicom file, you always need to convert the pixel_array into gray image, so that you can process this gray image by a lot of image processing tool in python.

The following code is a working example to convert pixel_array into gray image.

import matplotlib.pyplot as plt
import os
import pydicom
import numpy as np 
# Abvoe code is to import dependent libraries of this code

# Read some CT dicom file here by pydicom library
ct_filepath = r"<YOUR_CT_DICOM_FILEPATH>"
ct_dicom = pydicom.read_file(ct_filepath)
img = ct_dicom.pixel_array

# Now, img is pixel_array. it is input of our demo code

# Convert pixel_array (img) to -> gray image (img_2d_scaled)
## Step 1. Convert to float to avoid overflow or underflow losses.
img_2d = img.astype(float)

## Step 2. Rescaling grey scale between 0-255
img_2d_scaled = (np.maximum(img_2d,0) / img_2d.max()) * 255.0

## Step 3. Convert to uint
img_2d_scaled = np.uint8(img_2d_scaled)


# Show information of input and output in above code
## (1) Show information of original CT image 
print(img.dtype)
print(img.shape)
print(img)

## (2) Show information of gray image of it 
print(img_2d_scaled.dtype)
print(img_2d_scaled.shape)
print(img_2d_scaled)

## (3) Show the scaled gray image by matplotlib
plt.imshow(img_2d_scaled, cmap='gray', vmin=0, vmax=255)
plt.show()

And the following is result of what I print out.

Tow

Answer 2

You probably worked around this by now, but I think pydicom doesn't interpret planar configuration correctly.

You need to do this first:

img = data_set.pixel_array
img = img.reshape([img.shape[1], img.shape[2], 3])

From here on your image will have shape [rows cols 3], with the channels separated

Cachepot

Thank you very much! I'll give it a try soon. Actually my solution was writing of an ITK method. – Rand

Answer 3

As said by @Daniel since you have a PlanarConfiguration== 1 you have to rearrange your colors in columns through np.reshape and then converting to grayscale, for example using OpenCV:

import pydicom as dicom
import numpy as np
import cv2 as cv

data_set = dicom.read_file(path)
pixel_array = data_set.pixel_array
## converting to shape (m,n,3)
pixel_array_rgb = pixel_array.reshape((pixel_array.shape[1], pixel_array.shape[2], 3))
## converting to grayscale
pixel_array_gs = cv.cvtColor(pixel_array_rgb, cv.COLOR_RGB2GRAY)

Sucrose