0

I am working with nifti images (CT scans) and when I normalize the Hounsfield Unit converted images, the output is just a full black image. My code to convert to HU scale is as follows:

def transform_to_hu(img_data, img_obj): 
    
    slope = img_obj.dataobj.slope
    intercept = img_obj.dataobj.inter
    
    img_data[img_data >= 1200] = 0
    images = slope * img_data.astype(np.float64)
    images += np.float64(intercept)
    
    return np.array(images, dtype=np.float64)

hu_scans = transform_to_hu(img_data, img_obj)

And then I normalize the HU converted images using following function:

def normalize(volume):
    level = 70 
    window = 200 
    
    max = level + window/2 
    min = level - window/2 
    volume = volume.clip(min,max)
    
    """Normalize the volume"""
    volume[volume < min] = min
    volume[volume > max] = max
    volume = (volume - min) / (max - min)
    volume = volume.astype("float32")
    return volume

normalized_image = normalized(hu_scans)

Why is output image black after normalization?

enter image description here

Dushi Fdz
  • 161
  • 3
  • 22
  • What is your pixel type? It looks like you're normalizing to the range of 0-1. If the pixel type is not floats, it'll all be zeroes. – Dave Chen May 25 '22 at 13:50
  • pixel data type is float64 – Dushi Fdz May 25 '22 at 14:23
  • It could be your display method. It could be assuming pixel intensities are on the 0-255 range, so 0-1 pixels will all look black. Trying multiplying everything by 255 just to see. – Dave Chen May 25 '22 at 15:37

1 Answers1

0

Assuming you're using nibabel, you can simply use the get_fdata() method to load the CT as a numpy array with float64 dtype. If your Nifti headers are set correctly, it should automatically apply slope and intercept scaling.

To window and normalize your array to [0, 1], you can use the following code:

import numpy as np
import skimage.exposure

def window(data: np.ndarray, lower: float = -125., upper: float = 225., dtype: str = 'float32') -> np.ndarray:
    """ Scales the data between 0..1 based on a window with lower and upper limits as specified. dtype must be a float type.

    Default is a soft tissue window ([-125, 225] ≙ W 350, L50).

    See https://radiopaedia.org/articles/windowing-ct for common width (WW) and center/level (WL) parameters.
    """
    assert 'float' in dtype, 'dtype must be a float type'
    clipped = np.clip(data, lower, upper).astype(dtype)
    # (do not use in_range='image', since this does not yield the desired result if the min/max values do not reach lower/upper)
    return skimage.exposure.rescale_intensity(clipped, in_range=(lower, upper), out_range=(0., 1.))

Using PIL, you can get an image from your data by using:

from PIL import Image

img = Image.fromarray((window(data[..., 55].T) * 255.).round().astype('uint8'))  # axial slice number 55
asdf
  • 1,006
  • 1
  • 9