text_plain_output_1.png image_output_1.png
import matplotlib.pyplot as plt

import numpy as np

import pydicom

def make_lut(pixels, width, center, p_i):
    slope = 1.0
    intercept = 0.0
    min_pixel = int(np.amin(pixels))
    max_pixel = int(np.amax(pixels))
    lut = [0] * (max_pixel + 1)
    invert = False
    if p_i == 'MONOCHROME1':
        invert = True
    else:
        center = max_pixel - min_pixel - center
    for pix_value in range(min_pixel, max_pixel):
        lut_value = pix_value * slope + intercept
        voi_value = ((lut_value - center) / width + 0.5) * 255.0
        clamped_value = min(max(voi_value, 0), 255)
        if invert:
            lut[pix_value] = round(255 - clamped_value)
        else:
            lut[pix_value] = round(clamped_value)
    return lut

def apply_lut(pixels_in, lut):
    pixels_in = pixels_in.flatten()
    pixels_out = [0] * len(pixels_in)
    for i in range(0, len(pixels_in)):
        pixel = pixels_in[i]
        pixels_out[i] = int(lut[pixel])
    return pixels_out

from pydicom.pixel_data_handlers.util import apply_voi_lut
image = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00012/T1w/Image-17.dcm')
pixels = image.pixel_array

# Plot a histogram of the raw pixel data
fig, axes = plt.subplots(nrows=1, ncols=1,sharex=False, sharey=False, figsize=(10,4))
plt.title('Pixel Range: ' + str(np.min(pixels)) + '-' + str(np.max(pixels)))
plt.hist(pixels.ravel(), np.max(pixels), (1, np.max(pixels)))
plt.tight_layout()
plt.show()

window_width_1 = np.max(image.pixel_array)
window_center_1 = window_width_1 / 2
lut = make_lut(image.pixel_array, window_width_1, window_center_1, image.PhotometricInterpretation)
image1 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0], pixels.shape[1]))
window_width_2 = 450
window_center_2 = 450
lut = make_lut(image.pixel_array, window_width_2, window_center_2, image.PhotometricInterpretation)
image2 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0], pixels.shape[1]))
window_width_3 = 900
window_center_3 = 90
lut = make_lut(image.pixel_array, window_width_3, window_center_3, image.PhotometricInterpretation)
image3 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0], pixels.shape[1]))

fig, axes = plt.subplots(nrows=2, ncols=2, sharex=True, sharey=True, figsize=(12, 12))
ax = axes.ravel()
ax[0].set_title('Default Image')
ax[0].imshow(image.pixel_array, cmap='gray')
ax[1].set_title(f'Width: {window_width_1} / Center: {window_center_1}')
ax[1].imshow(image1, cmap='gray')
ax[2].set_title(f'Width: {window_width_2} / Center: {window_center_2}')
ax[2].imshow(image2, cmap='gray')
ax[3].set_title(f'Width: {window_width_3} / Center: {window_center_3}')
ax[3].imshow(image3, cmap='gray')
plt.tight_layout()
plt.show()