# z-maps for comparing PET images

**URL:** <https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414>\
**Category:** Algorithms\
**Tags:** simpleitk, registration, python\
**Created:** [October 14, 2022, 12:47am UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414 "2022-10-14T00:47:04Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Keyn34](https://discourse.itk.org/letter_avatar_proxy/v4/letter/k/9fc29f/32.png) [@Keyn34](https://discourse.itk.org/u/Keyn34)\
**Post date:** [October 14, 2022, 12:47am UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414/1 "2022-10-14T00:47:04Z")

</div>

Hello all,

I am trying to compute [z-maps](https://rstudio-pubs-static.s3.amazonaws.com/52340_0e269434591d4de790f259a330d86c1a.html) 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

```auto
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?):

```auto
    # 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!

---

<div class="post-metadata">

**Author:** ![dchen](https://discourse.itk.org/user_avatar/discourse.itk.org/dchen/32/34_2.png) [@dchen](https://discourse.itk.org/u/dchen)\
**Post date:** [October 14, 2022, 4:36am UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414/2 "2022-10-14T04:36:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Keyn34](https://discourse.itk.org/letter_avatar_proxy/v4/letter/k/9fc29f/32.png) [@Keyn34](https://discourse.itk.org/u/Keyn34)\
**Post date:** [October 14, 2022, 11:09am UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414/3 "2022-10-14T11:09:51Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dchen](https://discourse.itk.org/user_avatar/discourse.itk.org/dchen/32/34_2.png) [@dchen](https://discourse.itk.org/u/dchen)\
**Post date:** [October 14, 2022, 1:31pm UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414/4 "2022-10-14T13:31:14Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lassoan](https://discourse.itk.org/user_avatar/discourse.itk.org/lassoan/32/27_2.png) [@lassoan](https://discourse.itk.org/u/lassoan)\
**Post date:** [October 14, 2022, 1:53pm UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414/5 "2022-10-14T13:53:53Z")

</div>

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.

 ![image](https://discourse.itk.org/uploads/default/original/2X/d/da0c86d389638cb854a5a9ad1bbb544ec2c25845.jpeg)

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):

```python
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)

```

---

<div class="post-metadata">

**Author:** ![Keyn34](https://discourse.itk.org/letter_avatar_proxy/v4/letter/k/9fc29f/32.png) [@Keyn34](https://discourse.itk.org/u/Keyn34)\
**Post date:** [February 10, 2023, 9:07am UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414/6 "2023-02-10T09:07:32Z")

</div>

Hello everyone,

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:

```auto
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

```

Thank you for all your suggestions, as always 🙂

---

<div class="post-metadata">

**Author:** ![zivy](https://discourse.itk.org/user_avatar/discourse.itk.org/zivy/32/1726_2.png) [@zivy](https://discourse.itk.org/u/zivy)\
**Post date:** [February 10, 2023, 5:09pm UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414/7 "2023-02-10T17:09:32Z")

</div>

Hello @Keyn34,

Thank you for sharing.

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:

```auto
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

```

To upgrade to the latest pre-release:

```auto
pip install --upgrade --pre SimpleITK --find-links https://github.com/SimpleITK/SimpleITK/releases/tag/latest

```

---

<div class="post-metadata">

**Author:** ![Keyn34](https://discourse.itk.org/letter_avatar_proxy/v4/letter/k/9fc29f/32.png) [@Keyn34](https://discourse.itk.org/u/Keyn34)\
**Post date:** [February 10, 2023, 9:51pm UTC](https://discourse.itk.org/t/z-maps-for-comparing-pet-images/5414/8 "2023-02-10T21:51:59Z")

</div>

Hi @zivy,

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 😄

Thank you so much!