python-3.xnormalizationpixeldicommedical-imaging

Image note correctly displayed after transforming to Hounsfield Units

I am working on CT scans and specifically interested in the liver area. I am trying to convert pixel values to Hounsfield Units using the following function in python:

def transform_to_hu(slices): 
    images = np.stack([file.pixel_array for file in slices], axis=-1) #axis=-1 for depth in the last channel
    images = images.astype(np.int16)

    for n in range(len(slices)):
        
        intercept = slices[n].RescaleIntercept
        slope = slices[n].RescaleSlope
        
        if slope != 1:
            images[n] = slope * images[n].astype(np.float64)
            images[n] = images[n].astype(np.int16)
            
        images[n] += np.int16(intercept)
    
    return np.array(images, dtype=np.int16)

After transforming to HU, why does the image looks like it's separated into two regions?

enter image description here

Solution

  • Your n variable is the first index of the numpy array (which corresponds to coronal slices), while you iterate through the number of slices when you apply the rescale operation. Because the number of slices is less than the number of rows the rescale operation doesn't cover the entire volume.

    You should be iterating through the axial slices (i.e using the last index of the numpy array images[..., n]). Here's an example how to do it with pydicom's apply_rescale() function:

    from pydicom.data import get_testdata_file
from pydicom.pixel_data_handlers import apply_rescale
import numpy as np

ds = get_testdata_file("CT_small.dcm")
ds_stack = [ds, ds, ds]

images = np.stack([ds.pixel_array for ds in ds_stack], axis=-1).astype('float64')
for idx, ds in enumerate(ds_stack):
    images[..., idx] = apply_rescale(images[..., idx], ds)