# Algorithm 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.

# References

• Sliding organ registration
• Sliding geometries registration
• D. F. Pace, M. Niethammer, and S. R. Aylward, “Sliding geometries in deformable image registration,” Workshop on Computational and Clinical Applications in Abdominal Imaging, Medical Image Computing and Computer-Assisted Intervention–MICCAI 2011, Lecture Notes in Computer Science, vol. 7029, 2012 – In Press – PMC in Process.

# Running Sliding Geometries 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.

TODO

# Parameters

• Input Images
• Fixed Image Volume: Fixed image for registration.
• Moving Image Volume: Moving image for registration.
• Output Settings
• Output Displacement Field: (optional) Filename to which to save the estimated displacement field, which maps from fixed image coordinates to moving image coordinates.
• Output Volume: (optional) Moving image resampled to fixed image coordinates using the estimated displacement field.
• Input Sliding Boundaries, Normals, Weights and Geometries
• Input Sliding Surface Boundary Model: (optional)Surface model defining surfaces where sliding discontinuities can occur, in the space of the fixed image. Either: (1) specify an input sliding surface boundary model and/or tube spatial object, or (2) provide a normal vector image and sliding weight image, and optionally a sliding geometries image. If both are provided, the input sliding surface boundary model and/or tube spatial object is used. Typically, do first registration that inputs a sliding surface boundary model and/or tube spatial object and outputs a normal vector image, sliding weight image and/or sliding geometries image, which prevent duplicate computation on subsequent registrations.
• Input Sliding Tube Spatial Object: (optional)Tube spatial object file (.tre) defining tube surfaces where sliding discontinuities can occur, in the space of the fixed image. Use with sliding geometries registration only. Either: (1) specify an input sliding surface boundary model and/or tube spatial object, or (2) provide a normal vector image, sliding weight image and sliding geometries image. If both are provided, the input sliding surface boundary model and/or tube spatial object is used. Typically, do first registration that inputs a sliding surface boundary model and/or tube spatial object and outputs a normal vector image, sliding weight image and sliding geometries image, which prevent duplicate computation on subsequent registrations.
• Input Normal Vector Image: (optional)Image specifying the normal vectors of geometries where sliding can occur, in the space of the fixed image. For sliding organ registration, this is an image of 3D vectors. For sliding geometries registration, this is an image of 3x3 matrices. Either: (1) specify an input sliding surface boundary model and/or tube spatial object, or (2) provide a normal vector image and sliding weight image, and optionally a sliding geometries image. If both are provided, the input sliding surface boundary model and/or tube spatial object is used. Typically, do first registration that inputs a sliding surface boundary model and/or tube spatial object and outputs a normal vector image, sliding weight image and/or sliding geometries image, which prevent duplicate computation on subsequent registrations.
• Input Sliding Weight Image: (optional)Image specifying the spatially-varing weight between the anisotropic diffusive (sliding) regularization and the diffusive (smoothing) regularization, in the space of the fixed image. Typically a function of distance to the sliding boundaries. Either: (1) specify an input sliding surface boundary model and/or tube spatial object, or (2) provide a normal vector image and sliding weight image, and optionally a sliding geometries image. If both are provided, the input sliding surface boundary model and/or tube spatial object is used. Typically, do first registration that inputs a sliding surface boundary model and/or tube spatial object and outputs a normal vector image, sliding weight image and/or sliding geometries image, which prevent duplicate computation on subsequent registrations.
• Input Sliding Geometries Image: (optional)Image specifying the spatially-varying sliding geometries, in the space of the fixed image. Use with sliding geometries registration only. Either: (1) specify an input sliding surface boundary model and/or tube spatial object, or (2) provide a normal vector image, sliding weight image and sliding geometries image. If both are provided, the input sliding surface boundary model and/or tube spatial object is used. Typically, do first registration that inputs a sliding surface boundary model and/or tube spatial object and outputs a normal vector image, sliding weight image and/or sliding geometries image, which prevent duplicate computation on subsequent registrations.
• Output Sliding Boundaries, Normals, Weights and Geometries
• Output Normal Vector Image: (optional)Image specifying the normal vectors of geometries where sliding can occur, in the space of the fixed image. For sliding organ registration, this is an image of 3D vectors. For sliding geometries registration, this is an image of 3x3 matrices. Typically, do first registration that inputs a sliding surface boundary model and/or tube spatial object and outputs a normal vector image, sliding weight image and/or sliding geometries image, which prevent duplicate computation on subsequent registrations.
• Output Sliding Weight Image: (optional)Image specifying the spatially-varing weight between the anisotropic diffusive (sliding) regularization and the diffusive (smoothing) regularization, in the space of the fixed image. Typically a function of distance to the sliding boundaries. Typically, do first registration that inputs a sliding surface boundary model and/or tube spatial object and outputs a normal vector image, sliding weight image and/or sliding geometries image, which prevent duplicate computation on subsequent registrations.
• Output Sliding Geometries Image: (optional)Image specifying the spatially-varying sliding geometries, in the space of the fixed image. Use with sliding geometries registration only. Typically, do first registration that inputs a sliding surface boundary model and/or tube spatial object and outputs a normal vector image, sliding weight image and or sliding geometries image, which prevent duplicate computation on subsequent registrations.
• Initialization of registration
• Initial Transform: (optional)Filename of transform used to initialize the registration. If both an initial transform and initial displacement field are given, the initial displacement field will be used.
• Initial Displacement Field: (optional)Filename of displacement field used to initialize the registration. If both an initial transform and initial displacement field are given, the initial displacement field will be used.
• Main Parameters
• Time Step: Duration of each step in the optimization.
• Regularization Weighting: Relative weighting of the regularization term compared to the intensity distance term.
• Background Fill Value: Background fill value for the output image.
• Stopping Condition Parameters
• Iterations (determines multiresolution levels): Comma-separated list of the maximum number of iterations to try at each multi-resolution level before failing to converge. The registration will stop when either the absolute energy change is reached in a stopping criterion evaluation period, or if the specified number of iterations have been executed.
• Energy Change Stopping Criterion: Stopping criterion is defined as when the absolute energy change over the stopping criterion evaluation period is less than the energy change stopping criterion. Specify a larger energy change stopping criterion to halt faster. If zero, the stopping criterion is defined only by the number of iterations.
• Stopping Criterion Evaluation Period: The number of iterations that elapse between evaluations of the energy change stopping criterion.
• Stopping Criterion Mask Image: (optional)Mesh specifying which voxels should be used when evaluating the stopping criterion. Only pixels with non-zero value will be included when calculating the energy determining the stopping criterion.
• Lambda: Controls the decay of the spatially-varying sliding weighting as a function of the distance to the closest sliding geometry surface. If gamma=-1, the decay is exponential as parameterized by lambda. Otherwise, the decay is a dirac-shaped function, in which case gamma should be positive.
• Gamma: Controls the decay of the spatially-varying sliding weighting as a function of the distance to the closest sliding geometry surface. If gamma=-1, the decay is exponential as parameterized by lambda. Otherwise, the decay is a dirac-shaped function, in which case gamma should be positive.
• Risky Expert-Only Parameters
• Do Not Perform Regularization: If on, do not perform displacement field regularization.
• Do Not Compute Intensity Distance Term: If on, regularize initial displacement field without computing the intensity distance term.
• Do Not Use Sliding Regularization: If on, use the diffusive regularization (i.e. Gaussian smoothing).
• Sliding Regularization Type: Formulation used for the sliding anisotropic diffusive regularization. SlidingOrgan registration handles sliding at surfaces. SlidingGeometries registration considers planar, tubular and point-like geometries.
• World Coordinate System: The world coordinate system. RAS if the fixed image and sliding geometries are aligned when visualized in 3D Slicer. LPS if the fixed image and sliding geometries are aligned when visualized in Paraview.

# 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