To answer these questions we will render a basic 2D image with an interacive slicer to apply a rotation transform under different conditions. The answer is at the end.
In [1]:
## Boiler plate code common to many notebooks. See the TestFilesCommonCode.ipynb for details
from __future__ import print_function
%run TestFilesCommonCode.ipynb
In [2]:
img = sitk.ReadImage(fdata("cthead1.png"))
myshow(img, title="Input Image")
In [3]:
def callback(transform, angle, fixed, moving, title):
"""A call back method to resample an image with a active 2D Euler Transform"""
# Here we use the old SimpleITK 0.8 transform inferface which will work with composites too
p = list(transform.GetParameters())
p[0] = angle
transform.SetParameters(p)
resampler = sitk.ResampleImageFilter()
resampler.SetTransform(transform)
resampler.SetReferenceImage(fixed)
resampledMoving = resampler.Execute(moving)
myshow(sitk.Compose(resampledMoving, resampledMoving, sitk.Maximum(resampledMoving,fixed)), title=title)
In [4]:
tx = sitk.Euler2DTransform()
interact(lambda **kwargs: callback(transform=tx, fixed=img, moving=img, title="Origin Center Of Rotation", **kwargs),
angle=widgets.FloatSliderWidget(min=-pi, max=pi, step=.02, value=0.0)
)
The Above interaction is default state of a EulerTransform where the rotaion occours around the image's origin.
In [5]:
tx2 = sitk.Euler2DTransform(sitk.CenteredTransformInitializer(img,img,sitk.Euler2DTransform()) )
interact(lambda **kwargs: callback(tx2, fixed=img, moving=img, title="Mass Center Of Rotation", **kwargs),
angle=widgets.FloatSliderWidget(min=-pi, max=pi, step=.02, value=0.0) )
In [6]:
# here we create a fixed and moving image resampled onto a larger image size,
# and move the fixed image (shown in blue) to about the center of the new space
fixed = sitk.Image(img)
moving = sitk.Image(img)
f = sitk.ResampleImageFilter()
f.SetReferenceImage(fixed)
f.SetSize([512]*2)
moving = f.Execute(moving)
fixed.SetOrigin([40,40])
fixed = f.Execute(fixed)
title="ImageTitle"
myshow(sitk.Compose(moving, moving, sitk.Maximum(moving,fixed)), title=title)
In [7]:
# Case 1:
# Here we form a composite transform with the first transform solving the alignment
# The second transform is an identity
ctx1 = sitk.Euler2DTransform()
ctx1.SetCenter([45.8801,48.03]) # This parameter has no effect
ctx1.SetTranslation([-40,-40]) # This is the correct solution
ctx2 = sitk.Euler2DTransform()
#ctx2.SetCenter([45.8801+40, 48.03+40])
ctx = sitk.Transform(ctx1)
ctx.AddTransform(ctx2)
interact(lambda **kwargs: callback(ctx, fixed=fixed, moving=moving, title="TX2 Transform", **kwargs),
angle=widgets.FloatSliderWidget(min=-pi, max=pi, step=.01, value=0.0) )
In [8]:
# Case 2:
# Here the first transform is only part of the solution.
# Propagating the center, make the moving image rotate around that point in the fixed image domain.
ctx1=sitk.Euler2DTransform(sitk.CenteredTransformInitializer(fixed,moving,sitk.Euler2DTransform()) )
ctx1.SetTranslation([-20,-20])
invTx = sitk.Transform(ctx1)
invTx.SetInverse()
ctx2 = sitk.Euler2DTransform()
ctx2.SetTranslation([-20,-20])
ctx2.SetCenter(ctx1.GetCenter())
# This is not the correct center
#ctx2.SetCenter(invTx.TransformPoint(ctx1.GetCenter()))
ctx = sitk.Transform(ctx1)
ctx.AddTransform(ctx2)
interact(lambda **kwargs: callback(ctx, fixed=fixed, moving=moving, title="TX2 Transform", **kwargs),
angle=widgets.FloatSliderWidget(min=-pi, max=pi, step=.01, value=0.0) )
In [9]:
# Case 3:
# This demonstrate the behavior of Case 2 is consistent
# with resampling the first transform then applying the
# second transfrom the the resampled image.
ctx1=sitk.Euler2DTransform(sitk.CenteredTransformInitializer(fixed,moving,sitk.Euler2DTransform()) )
ctx1.SetTranslation([-20,-20])
f = sitk.ResampleImageFilter()
f.SetReferenceImage(fixed)
f.SetTransform(ctx1)
resampledMoving = f.Execute(moving)
invTx = sitk.Transform(ctx1)
invTx.SetInverse()
myshow(sitk.Compose(resampledMoving, resampledMoving, sitk.Maximum(resampledMoving,fixed)), title="Resampled")
ctx2 = sitk.Euler2DTransform()
ctx2.SetTranslation([-20,-20])
ctx2.SetCenter(ctx1.GetCenter())
#ctx2.SetCenter(invTx.TransformPoint(ctx1.GetCenter()))
interact(lambda **kwargs: callback(ctx2, fixed=fixed, moving=resampledMoving, title="Resampled Moving with Transform", **kwargs),
angle=widgets.FloatSliderWidget(min=-pi, max=pi, step=.01, value=0.0) )
For a Composite transformed defined as $\phi=T_0(T_1(\mathbf{x}))$ where $\phi:\Omega_{fixed}\rightarrow \Omega_{moving}$ then if $\mathbf{C}$ such that $\mathbf{C} \in \Omega_{fixed}$ is the correct center initialization for $T_0$, then the correct center initialization for $T_1$ is still $\mathbf{C}$, as it is set in the in the domain of the transforms.
In [ ]: