Differences after applying forward & inverse transform

Beginner Questions
1

Hi all,

I started using sitk a while ago and hope to use it for registration tasks in my project.
While playing around with it, I tried to apply a simple Euler3D transform to my image and revert it with the the inverse transform.

Expectation: I should get my original image back with minor errors.

Result: I get rather large absolute errors in the order of double digit percentages.

I assume I missed something and this is intended behavior?
Any tips and insights are appreciated.

My code runs in a Jupyter Notebook in VS Code:
python 3.12.2
sitk 2.3.1

import numpy as np
import matplotlib.pyplot as plt
import SimpleITK as sitk

fixed_image = sitk.ReadImage("../data/training_001_ct.mha", sitk.sitkFloat32)

p1_transform = sitk.Euler3DTransform()
p1_transform.SetCenter(fixed_image.TransformContinuousIndexToPhysicalPoint(np.array(fixed_image.GetSize()) / 2.0))

p1_transform.SetTranslation([0, 0, 10])
p1_transform.SetRotation(0, 0, np.deg2rad(5))

moving_image = sitk.Resample(fixed_image, p1_transform, interpolator=sitk.sitkLinear)
moving_resampled = sitk.Resample(moving_image, p1_transform.GetInverse(), interpolator=sitk.sitkLinear)

difference = sitk.GetArrayFromImage(fixed_image)-sitk.GetArrayFromImage(moving_resampled)

plt.imshow(np.abs(difference[15,:,:]), interpolation='none', vmax=700, cmap="Greys_r")
plt.title("|difference| ; fixed-moving")
plt.colorbar()

output

2

Hello @Mckay42,

Your expectation is not valid. What you are experiencing is due to interpolation and is not surprising when dealing with digital signals. Repeated interpolation is an issue. Also, because this specific image has a pretty large spacing in the z direction and you are translating in that direction the effects of using linear interpolation are more significant. Generally speaking repeated interpolation is discouraged (see this notebook section titled Repeated Resampling - Just Say No).

The choice of interpolation method also plays a role, available interpolators are listed here. Linear is usually a compromise between interpolation quality and computational complexity.
Code below illustrates this using a 2D image (easier to visualize):

import numpy as np
import SimpleITK as sitk

fixed_image = sitk.ReadImage("training_001_ct.mha", sitk.sitkFloat32)
mid_slice = fixed_image.GetDepth() // 2
fixed_image = fixed_image[:,:,mid_slice]

p1_transform = sitk.Euler2DTransform()
p1_transform.SetCenter(fixed_image.TransformContinuousIndexToPhysicalPoint(np.array(fixed_image.GetSize()) / 2.0))

p1_transform.SetTranslation([0, 10])
p1_transform.SetAngle(np.deg2rad(5))

interps = {"linear": sitk.sitkLinear,
           "BSpline3": sitk.sitkBSpline3}

for s, interp in interps.items(): 
    moving_image = sitk.Resample(fixed_image, p1_transform, interpolator=interp)
    moving_resampled = sitk.Resample(moving_image, p1_transform.GetInverse(), interpolator=interp)
    diff_abs = sitk.Abs(fixed_image - moving_resampled)[50:-50,50:-50] #cropping to ignore changes introduced at image edges by default value
    print(f"mean absolute difference using {s}: {np.sum(sitk.GetArrayViewFromImage(diff_abs))/diff_abs.GetNumberOfPixels()}")
2 Likes
3

Thanks for the quick response!
After playing around some more it was already dawning on me that his might not be the best idea.