Jovian
Sign In
Learn practical skills, build real-world projects, and advance your career
%matplotlib inline

import cv2
import pydicom
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.ndimage as ndimage

from pydicom.data import get_testdata_files
from skimage import measure, morphology, segmentation
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

INPUT_PATH = "D:\\Documents\\LCTSC Test Data"
patients = os.listdir(INPUT_PATH)

def load_scan(path):
    
    slices = [pydicom.dcmread(path + '\\' + s) for s in os.listdir(path)]
    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
    try:
        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
    except:
        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
        
    for s in slices:
        s.SliceThickness = slice_thickness
        
    
    return slices


    #patient = pydicom.dcmread("D:\\Documents\\LCTSC Test Data\\LCTSC-Test-S1-101\\03-03-2004-08186\\79262\\000010.dcm")
    #img = patient.pixel_array
    #plt.imshow(img, cmap=plt.cm.CMRmap)
def get_pixels_hu(slices):
    image = np.stack([s.pixel_array for s in slices])
    
    # Convert to int16 (from sometimes int16), 
    # should be possible as values should always be low enough (<32k)
    
    
    image = image.astype(np.int16)

    # Set outside-of-scan pixels to 0
    # The intercept is usually -1024, so air is approximately 0
    
    
    image[image == -2000] = 0
    
    # Convert to Hounsfield units (HU)
    for slice_number in range(len(slices)):
        
        intercept = slices[slice_number].RescaleIntercept
        slope = slices[slice_number].RescaleSlope
        
        if slope != 1:
            image[slice_number] = slope * image[slice_number].astype(np.float64)
            image[slice_number] = image[slice_number].astype(np.int16)
            
        image[slice_number] += np.int16(intercept)
    
    return np.array(image, dtype=np.int16)

The unit of measurement in CT scans is the Hounsfield Unit (HU), which is a measure of radiodensity. CT scanners are carefully calibrated to accurately measure this. From Wikipedia:

HU examples

hounsfield.png

By default however, the returned values are not in this unit. Let's fix this.

first_patient = load_scan(INPUT_PATH + '\\' + patients[10])
first_patient_pixels = get_pixels_hu(first_patient)
print (first_patient_pixels)
plt.hist(first_patient_pixels.flatten(), bins=80, color='c')
plt.xlabel("Hounsfield Units (HU)")
plt.ylabel("Frequency")
plt.show()

# Show some slice in the middle

new_pix = first_patient_pixels[60]

plt.imshow(new_pix, cmap='plasma_r')
plt.show()
[[[-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] ... [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024]] [[-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] ... [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024]] [[-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] ... [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024]] ... [[-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] ... [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024]] [[-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] ... [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024]] [[-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] ... [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024] [-1024 -1024 -1024 ... -1024 -1024 -1024]]]
Notebook Image
Notebook Image