# Overview

TubeTK provides algorithms for deformable registration of images depicting multiple organs in which the organs may have shifted, expanded, or compressed independently.

Traditional deformable registration imposes a uniform smoothness constraint on the deformation field. However, discontinuities in the deformation field are expected with sliding motion, and this constraint is not appropriate. This ultimately leads to registration inaccuracies.

TubeTK provides deformable image registration incorporating a deformation field regularization term that is based on anisotropic diffusion. A cost function ${\displaystyle C(u)}$ is a function of the current estimation of the deformation field ${\displaystyle u}$. ${\displaystyle C(u)}$ is iteratively optimized using finite differences and is the sum of of two terms:

• Intensity-based distance measure: captures intensity differences between the fixed image and the transformed moving image (sum of squared differences)
• Anisotropic diffusive regularization term: penalizes unrealistic deformation fields, while considering sliding motion

The anisotropic diffusive regularization is based on decomposing the deformation field into normal and tangential components, which are defined with respect to a given organ boundary along which sliding motion is expected to occur. These two components are handled differently:

• Motion normal to the organ boundary should be smooth both across organ boundaries and deep within organs. The motion normal to the organ boundary must be smooth in both the normal and tangential directions. The former condition enforces coupling between neighboring organs under the assumption that organs do not pull apart. The later forces smooth motion of individual organs.
• Motion tangential to the organ boundary should be smooth in the tangential direction within each individual organ. However, smoothness is not required across organ boundaries, therefore sliding transformations are not penalized.

These conditions are implemented by defining the anisotropic regularizer as:

${\displaystyle S_{\mathrm {a} }(u)={\frac {1}{2}}\sum _{l=x,y,z}\sum _{\mathbf {x} \in \Omega }\|P\nabla u_{l}(\mathbf {x} )\|^{2}+w\left(n^{T}\nabla u_{l}^{\perp }(\mathbf {x} )\right)^{2}}$

where

${\displaystyle P=I-wnn^{T}}$

and

• ${\displaystyle n}$ is the organ boundary in the vicinity of ${\displaystyle {\textbf {x}}}$,
• ${\displaystyle u(\mathbf {x} )}$ is the vector within the deformation field ${\displaystyle u}$ at location ${\displaystyle {\textbf {x}}}$
• ${\displaystyle \nabla u_{l}(\mathbf {x} )}$ is the gradient of the ${\displaystyle l}$-th component of ${\displaystyle u({\textbf {x}})}$
• ${\displaystyle u_{l}^{\perp }(\mathbf {x} )}$ is the component of ${\displaystyle u_{l}({\textbf {x}})}$ in the normal direction
• ${\displaystyle w}$ is a weighting term between the anisotropic diffusive and the diffusive regularizations, and decays exponentially from 1 to 0 as a function of distance to the organ boundary.

Close to organ boundaries, where ${\displaystyle w}$ is close to 1:

• ${\displaystyle \|P\nabla u_{l}(\mathbf {x} )\|^{2}}$ penalizes any discontinuities in the deformation field that are in the plane tangential to the organ boundary. This anisotropically smooths:
• Discontinuities in the deformation field's normal component that occur in the tangential plane
• Discontinuities in the deformation field's tangential component that occur in the tangential plane
• ${\displaystyle w\left(n^{T}\nabla u_{l}^{\perp }(\mathbf {x} )\right)^{2}}$ penalizes any discontinuities in the deformation field's normal component that occur in the normal direction
• Discontinuities in the deformation field's tangential component that occur in the normal direction are allowed: these are sliding motions!

The gradient is defined with respect to ${\displaystyle u}$ and is implemented in itkImageToImageDiffusiveDeformableRegistrationFilter and itkImageToImageDiffusiveDeformableRegistrationFunction using ITK's finite differences framework:

${\displaystyle c_{S_{\mathrm {a} }}\left(u(\mathbf {x} ,t)\right)=\sum _{l=x,y,z}{\textrm {div}}\left(P^{T}P\nabla u_{l}(\mathbf {x} )\right)(e_{l})+{\textrm {div}}\left(w\left(n^{T}\nabla u_{l}^{\perp }(\mathbf {x} )\right)n\right)n_{l}n}$

where ${\displaystyle e_{l}}$ is the ${\displaystyle l^{th}}$ canonical unit vector, i.e ${\displaystyle e_{x}=[1,0,0]^{T}}$

Further away from organ boundaries, ${\displaystyle w}$ approximates 0 and this tends to the diffusive regularization, which is equivalent of Gaussian smoothing. Therefore, uniformly smooth motion is required within each individual organ.

## Running Anisotropic Diffusive Deformable Registration

The anisotropic diffusive deformable registration algorithm is implemented as a command line module. It can be run from the command line, or graphically using 3D Slicer (see the TubeTK with Slicer page for setup instructions).

Its implementation can be found in TubeTK/Applications/CLI/tubeAnisotropicDiffusiveDeformableRegistration and TubeTK/Base/Registration.

## Parameters

• Iterations: The number of iterations to run the registration for. Up to a point, more iterations improve the registration results.
• Time Step: The duration of the time step used to solve the finite difference equations defined to solve the registration.
• Lambda: controls the exponential decay used to calculate the weighting value ${\displaystyle w}$, which decreases as a function of the distance to the organ boundary.
• Do not perform regularization: Whether or not to perform the regularization - should be False
• Do not use anisotropic regularization: Whether or not to use the anisotropic diffusive regularization. If True, the diffusive regularization (Gaussian smoothing) is used instead.
• Input Organ Boundary: The surface model that defines the organ boundary in the space of the fixed image.
• Input Normal Vector Image: The image of normal vectors derived from the organ boundary.
• Input Weight Image: The image of weighting values ${\displaystyle w}$,
• Initial Transform: An initial alignment of the fixed and moving images (NOT YET IMPLEMENTED),
• Input Fixed Image: The fixed image to be registered.
• Input Moving Image: The moving image to be registered (will be deformed to look more like the input fixed image).
• Output Deformation Field: The deformation field calculated to align the fixed and moving images - maps from the fixed coordinate frame to the moving coordinate frame, in order to deform the input moving image to look more like the input fixed image.
• Output Volume: The moving image after it is deformed into the fixed image frame.
• Output Normal Vector Image: The image of normal vectors calculated based on the input organ boundary.
• Output Weight Image: The image of weightings ${\displaystyle w}$ calculated based on the input organ boundary.

### Parameter usage notes

• The first time you run the anisotropic diffusive deformable registration, provide a surface model of the organ boundary, which may be created using Slicer's segmentation and Model Maker modules. Also, specify an output normal vector image and an output weight image. The next time you run the registration, you can set the input normal vector image and the input weight image to be the results of the first registration. Then, the normal vector image and the weight image will not have to be calculated, and the registration will be faster.

# Tutorial

## Tutorial Data

• In this example, we will register two artificial images of sliding tubes
• The fixed image, moving image and surface border model can be found in TubeTK/Data

## Tutorial Steps

• File -> Add Volume -> TubeTK/Data/SlidingOrganRegistrationWithRegularization_boxes_originalFixedImage.mhd
• File -> Add Volume -> TubeTK/Data/SlidingOrganRegistrationWithRegularization_boxes_originalMovingImage.mhd
• File -> Add Data -> Add File(s) -> TubeTK/Data/SlidingOrganRegistrationWithRegularization_boxes_normalSurfaceBorder.vtk -> OK
• Link the slice controls across all slice viewers by clicking the link button in the red slice viewer
• View the image in the 3D viewer by clicking the closed eye in the red slice viewer
• Set the foreground image to be the original fixed image, and the background image to be the original moving image
• You can toggle between the foreground and background images to see the difference between them
• Compared to the fixed image, the moving image translates the subset of the upper tube with increasing intensity to the right by four pixels, and the subset of the bottom tube with decreasing intensity to the left by four pixels.
• Original fixed image:
• Original moving image:
2. Open the Anisotropic Diffusive Deformable Registration module
3. Set the registration parameters
• Iterations: 500
• Time Step: 0.125
• Lambda: -0.1
• Do not perform regularization: OFF
• Do not use anisotropic regularization: OFF
• Input Organ Boundary: SlidingOrganRegistrationWithRegularization_boxes_normalSurfaceBorder.vtk
• Input Normal Vector Image: None
• Input Weight Image: None
• Initial Transform: None
• Input Fixed Image: SlidingOrganRegistrationWithRegularization_boxes_originalFixedImage
• Input Moving Image: SlidingOrganRegistrationWithRegularization_boxes_originalMovingImage
• Output Deformation Field: Create new volume (Volume)
• Output Volume: Create new volume (Volume1)
• Output Normal Vector Image: Create new volume (Volume2)
• Output Weight Image: Create new volume (Volume3)
4. Perform the registration
• Click "Apply" - the registration will take some time
5. View the registration results
• After the registration, the original fixed image and the deformed moving image (Volume1) should be similar
• Set the background image to the deformed moving image (Volume1). Toggling between the background and foreground images allows you to compare them.
• The window/level values for the original fixed image and the transformed moving image must be the same when you compare them visually. Open the 'Volumes' module, select the volume of interest, and change the two window/level parameters so that they match.
• If an image appears black when you know there should be data there, the window/level parameters are usually the culprit! Open the 'Volumes' module, select the volume of interest, and change the window/level setting to 'Automatic'
• Original fixed image:
• Transformed moving image:
• Subtraction between the original fixed image and the original moving image:
• Subtraction between the original fixed image and the transformed moving image:
• You can also view the output deformation field (an image of 3D vectors), the normal vector image (an image of 3D vectors) and the weight image (a scalar image) by viewing Volume, Volume1 and Volume3, respectively
• Compare to the diffusive regularization
• You can compare the results of the anisotropic diffusive regularization with those of the more traditional diffusive regularization (i.e. Gaussian smoothing) by setting the "Do not use anisotropic regularization" parameter to OFF.
• The results with our example are shown below:
• Transformed moving images are similar in both cases
• BUT the anisotropic diffusive regularization gives a deformation field that is much closer to the actual applied deformation

## Known issues

• There is a problem related to Slicer's coordinate transforms. If you use Slicer to create the surface models, you must rotate it 180 degrees in the x- and y- directions to create a new flipped model, and use the flipped model in the algorithm: The surface model and images should not be aligned in Slicer's 3D viewer. This flip is easily performed using Paraview's Transform filter. This issue will be addressed soon.
• Output deformation and normal vector images not correctly visualized in Slicer4

# Related Works

• Slipping objects in image registration: Improved motion field estimation with direction-dependent regularization
• Alexander Schmidt-Richberg, Jan Ehrhardt, Rene Werner and Heinz Handels
• MICCAI 2009, Lecture Notes in Computer Science, Volume 5761, pp.755-762, 2009
• Abstract:
• The computation of accurate motion fields is a crucial aspect in 4D medical imaging. It is usually done using a non-linear registration without further modeling of physiological motion properties. However, a globally homogeneous smoothing (regularization) of the motion field during the registration process can contradict the characteristics of motion dynamics. This is particularly the case when two organs slip along each other which leads to discontinuities in the motion field. In this paper, we present a diffusion-based model for incorporating physiological knowledge in image registration. By decoupling normal- and tangential-directed smoothing, we are able to estimate slipping motion at the organ borders while ensuring smooth motion fields in the inside and preventing gaps to arise in the field. We evaluate our model focusing on the estimation of respiratory lung motion. By accounting for the discontinuous motion of visceral and parietal pleurae, we are able to show a significant increase of registration accuracy with respect to the target registration error (TRE).
• An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences
• Hans-Hellmut Nagel and Wilfried Enkelmann
• IEEE Transactions on Pattern Analysis and Machine Intelligence, 8(5), pp 565-593, 1986
• http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4767833
• Abstract:
• A mapping between one frame from an image sequence and the preceding or following frame can be represented as a displacement vector field. In most situations, the mere gray value variations do not provide sufficient information in order to estimate such a displacement vector field. Supplementary constraints are necessary, for example the postulate that a displacement vector field varies smoothly as a function of the image position. Taken as a general requirement, this creates difficulties at gray value transitions which correspond to occluding contours. Nagel therefore introduced the oriented smoothness requirement which restricts variations of the displacement vector field only in directions with small or no variation of gray values. This contribution reports results of an investigation about how such an oriented smoothness constraint may be formulated and evaluated.
• A review of nonlinear diffusion filtering
• Joachim Weickert
• Scale-Space Theory in Computer Vision, Lecture Notes in Computer Science, Volume 1252, pp. 1-28, 1997