z-maps for comparing PET images

Hello all,

I am trying to compute z-maps for PET images to test it as a comparison methodology/metric with SimpleITK on Python.
The formula is simple, but I always have a completely black image with an insanely broad (10^32) intensity range, according to ITK-SNAP. I tried loading the image in ITK-SNAP and AMIDE, but no difference.

Would anyone of you have a hint of what might cause that issue? Wrong datatype while reading the images, faulty mathematical operation on the images?

I compute the z-map simply by

z_map = img_comp- mean_ref / std_ref

For calculating the mean and the std of a group of images, I use the following code (maybe something is wrong here?):

    # Compute the mean
    # Start with the first image
    summed_image_mean = sitk.ReadImage(image_paths[0])

    # Sum all other images
    for image_index in range(1, len(image_paths)):
        current_image = sitk.ReadImage(image_paths[image_index])
        summed_image_mean += current_image

    number_of_images = len(image_paths)
    avg_val = 1 / number_of_images
    mean_image = summed_image_mean * avg_val

    # Compute the STD
    # Start with the first image
    diff_image = sitk.ReadImage(image_paths[0]) - mean_image
    summed_image_std = diff_image**2

    # Process all other images
    for image_index in range(1, len(image_paths)):
        diff_image = sitk.ReadImage(image_paths[image_index]) - mean_image
        summed_image_std += diff_image**2

    std_image = (summed_image_std / number_of_images)**0.5

I would be very grateful if you have any suggestions for me.
Thank you all in advance!

What pixel type are your images? If they’re 8 or 16 bit integers, doing all that math can have underflow or overflow problems. When calling ReadImage, you can tell it to convert to floating point pixels.

1 Like

Thank you for the suggestions, but unfortunately, no luck. I tried reading all images with sitk.sitkFloat64, but the result is still the same.

How about if you look at the intermediate results that you’re computing? Do they look ok? Does the summed_image_mean seem correct?

If you call sitk.Show(summed_image_mean) you can look at the image in Fiji. And in Fiji if you mouse over the pixels you can see the actual pixel values. That’s useful when you have a very large range in values.

An image might look all black, but that could be because of a few outlier values. Mapping a wide range of floats to unsigned chars for display can make most pixels map to 0.

2 Likes

I’ve tested the code snippet in 3D Slicer and it worked well for the MRHead sample data set with short scalar type. So, the issue is most likely with the input data or how you load the input data or display the output.

In Slicer all the image slices are already in a 3D array, so I’ve changed file reading to numpy array slicing, but otherwise left everything unchanged (although having everything in a 3D array would allow further simplification of the code):

volumeNode = getNode('MRHead')
voxelArray = arrayFromVolume(volumeNode)

###

# Compute the mean
# Start with the first image
summed_image_mean = voxelArray[0,:,:].copy()

# Sum all other images
for image_index in range(1, voxelArray.shape[0]):
    current_image = voxelArray[image_index,:,:]
    summed_image_mean += current_image

number_of_images = voxelArray.shape[0]
avg_val = 1 / number_of_images
mean_image = summed_image_mean * avg_val

# Compute the STD
# Start with the first image
diff_image = voxelArray[0,:,:] - mean_image
summed_image_std = diff_image**2

# Process all other images
for image_index in range(1, voxelArray.shape[0]):
    diff_image = voxelArray[image_index,:,:] - mean_image
    summed_image_std += diff_image**2

std_image = (summed_image_std / number_of_images)**0.5

###

processedVolume = slicer.modules.volumes.logic().CloneVolume(volumeNode, "Processed image")
import numpy as np
std_image_3d = np.zeros([1, std_image.shape[0], std_image.shape[1]])
std_image_3d[0] = std_image
updateVolumeFromArray(processedVolume, std_image_3d)
3 Likes