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.
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.
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)
Sorry for the late reply. There were a lot of conferences at that time.
I fixed the problem by using numpy and only performing a division where the values in the images are higher than 0.01 and lower than -0.01. There might be a better way to do it, but this works fine now.
Here is the code if anyone might need it in the future:
def compute_z_map(sample: SimpleITK.Image, mean: SimpleITK.Image, std: SimpleITK.Image) -> SimpleITK.Image:
# The diff image is computed as SimpleITK.Image
diff_image = sample - mean
# The needed images for division as numpy arrays:
diff_image_np = SimpleITK.GetArrayFromImage(diff_image)
std_image_np = SimpleITK.GetArrayFromImage(std)
# Division only at voxels higher than 0.01 and lower than -0.01
z_map_np = np.divide(diff_image_np, std_image_np, where=((std_image_np >= 0.01) | (std_image_np <= -0.01)))
# Create a SimpleITK image from the numpy z-map
z_map = SimpleITK.GetImageFromArray(z_map_np)
z_map.SetOrigin(diff_image.GetOrigin())
z_map.SetSpacing(diff_image.GetSpacing())
z_map.SetDirection(diff_image.GetDirection())
return z_map
There appears to be a logical bug in the code, the standard deviation cannot be negative, so the where statement is incorrect, or possibly there is a problem with the computation of the standard deviation and this is just hiding it.
Note that with the latest and greatest SimpleITK pre-release you can simplify the function considerably or just remove it as it becomes a three liner:
def compute_z_map(sample: SimpleITK.Image, mean: SimpleITK.Image, std: SimpleITK.Image, std_positive_threshold: int=0.01) -> SimpleITK.Image:
# The diff image is computed as SimpleITK.Image
diff_image = sample - mean
# very small standard deviation, replace by 1, z value is the original value
std[std<std_positive_threshold] = 1
return sample - mean / std
Oops, thank you for pointing that out. That is a residue of previous iterations of the functions where I tried to escape negative values in the diff_image. I never had the time to revisit the code once it worked. Highly appreciate that you pointed it out. Thanks! I even checked for sanity reasons that all std images I created only have positive intensities.
Wow, that makes the function super sleek! Actually, the diff_image = sample - mean can also be removed. Only two lines of code for something that took me weeks to figure out and to pin the issue down
