How to do histogram equalization to DICOM images read with SimpleITK

Question

I'm trying to do histogram equalization to all images read from a *.nii.gz file.

I have tried this code:

import SimpleITK as sitk
flair_file = '/content/gdrive/My Drive/Colab Notebooks/.../FLAIR.nii.gz'

images = sitk.ReadImage(flair_file)
print("Width: ", images.GetWidth())
print("Height:", images.GetHeight())
print("Depth: ", images.GetDepth())

print("Dimension:", images.GetDimension())
print("Pixel ID: ", images.GetPixelIDValue())
print("Pixel ID Type:", images.GetPixelIDTypeAsString())

With this output:

Width:  240
Height: 240
Depth:  48
Dimension: 3
Pixel ID:  8
Pixel ID Type: 32-bit float

But when I try to do histogram equalization with OpenCV I get errors:

images_array = sitk.GetArrayFromImage(images)
gray = cv2.cvtColor(images_array[24], cv2.COLOR_BGR2GRAY)

Output:

error: OpenCV(4.1.2) /io/opencv/modules/imgproc/src/color.simd_helpers.hpp:92: error: (-2:Unspecified error) in function 'cv::impl::{anonymous}::CvtHelper<VScn, VDcn, VDepth, sizePolicy>::CvtHelper(cv::InputArray, cv::OutputArray, int) [with VScn = cv::impl::{anonymous}::Set<3, 4>; VDcn = cv::impl::{anonymous}::Set<1>; VDepth = cv::impl::{anonymous}::Set<0, 2, 5>; cv::impl::{anonymous}::SizePolicy sizePolicy = (cv::impl::<unnamed>::SizePolicy)2u; cv::InputArray = const cv::_InputArray&; cv::OutputArray = const cv::_OutputArray&]'
> Invalid number of channels in input image:
>     'VScn::contains(scn)'
> where
>     'scn' is 1

So, I have tried this other code:

images_array = sitk.GetArrayFromImage(images)
#gray = cv2.cvtColor(images_array[24], cv2.COLOR_BGR2GRAY)
output = cv2.equalizeHist(images_array[24])

But I get this error:

error: OpenCV(4.1.2) /io/opencv/modules/imgproc/src/histogram.cpp:3429: error: (-215:Assertion failed) _src.type() == CV_8UC1 in function 'equalizeHist'

How can I do histogram equalization to those DICOM images (maybe not using OpenCV, and using instead SimpleITK)?

UPDATE:
When I run this command:

print(images_array[24].shape, images_array[24].dtype)

I get this:

(240, 240) float32

What's print(images_array[24].shape, images_array[24].dtype)? — HansHirse
are you sure the images are BGR? my guess is that DICOM images are 16 bit greyscale, a lot are like that. If that is the case, you cannot do cv2.COLOR_BGR2GRAY, because it is already gray and does not have 3 channels. And you can't do cv2.equalizeHist because it is not 8 bit greyscale... Please print what HansHirse recommend to see what is happening. — api55
Hmmm, well I think you will have to implement it yourself :( At least is one channel, so it is not so complicated. Since it is a float mat, you have to create some bins (you decide the number of bins. Then get the cdf of the bins (ordered from less intensity to high intensity) and normalize it. When it is an 8 bit image is easier, since you only need to count the occurrence of each intensity, so in that case one has 256 bins and is easy to normalize.... are the values of your image discrete? what is the range? 0-255? 0-65536? — api55

Dave Chen Dave Chen · Accepted Answer · 2020-01-28T14:20:26

SimpleITK does have an AdaptiveHistogramEqualization function, and it does work on float32 images. Here's how you could use it:

new_images = sitk.AdaptiveHistogramEqualization(images)

Note that this would do equalization across the whole 3d image. If you wanted to do it on a slice-by-slice basis, it'd look something like this:

new_images = []
for z in range(images.GetDepth()):
    new_images.append(sitk.AdaptiveHistogramEqualization(images[:,:,z])

UPDATE: As pointed out by @blowekamp, AHE doesn't produce a global histogram equalization across the whole image, but a local equalization. Here is some example code that show's how to use the HistogramMatching function, as described by him, to achieve global histogram equalization.

import SimpleITK as sitk
import numpy as np

# Create a noise Gaussian blob test image
img = sitk.GaussianSource(sitk.sitkFloat32, size=[240,240,48], mean=[120,120,24])
img = img + sitk.AdditiveGaussianNoise(img,10)

# Create a ramp image of the same size
h = np.arange(0.0, 255,1.0666666666, dtype='f4')
h2 = np.reshape(np.repeat(h, 240*48), (48,240,240))
himg = sitk.GetImageFromArray(h2)
print(himg.GetSize())

# Match the histogram of the Gaussian image with the ramp
result=sitk.HistogramMatching(img, himg)

# Display the 3d image
import itkwidgets
itkwidgets.view(result)

Note that itkwidgets allows you to view a 3d image in a Jupyter notebook. You can see the histogram of the image there.

How to do histogram equalization to DICOM images read with SimpleITK

1 Answers