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

# Sliding Organ Registration Use Case

Sliding organ registration allows registration of images that show sliding motion at surfaces that are locally planar.

Sliding Regularization Type: SlidingOrgan


First, we set the fixed and moving images.

Fixed Image Volume
Moving Image Volume


Then, we specify an output displacement field and/or an output volume.

Output Displacement Field
Output Volume


A surface model (in the space of the fixed image) will determine the surfaces at which sliding is allowed to occur.

Input Sliding Surface Boundary Model


We will output the normal vector image ${\displaystyle n(x)}$ and sliding weight image ${\displaystyle w(x)}$, which are derived from the input sliding surface boundary model. These images can be used instead of the input sliding surface boundary model on subsequent registrations, so that they will not have to be recomputed if you decide to redo the registration with different parameters.

Output Normal Vector Image
Output Sliding Weight Image


Optionally, specify an initial transform or initial displacement field

Initial Transform
Initial Displacement Field


Next, specify the time step (must be small enough to prevent numerical instability) and the weighting ${\displaystyle \alpha }$ between the intensity distance term and the regularization term

Time Step
Regularization Weighting


Minimize errors near image borders by specifying a background fill value equal to the background intensity in the moving image. This will be used if the moving image is resampled outside of its domain.

Background Fill Value


The sliding organ registration can be used in a multiresolution (coarse to fine) framework. The multiresolution shrink factors are halved in all three directions at each multiresolution level. The number of levels is specified by the comma-separated list of iterations to perform at each multiresolution level.

Iterations


You can also specify a stopping criterion that will cause the registration to halt before the specified number of resolutions is reached. The energy change stopping criterion is defined as the absolute energy change that occurs in a stopping criterion evaluation period (note energy decreases over time). Provide a stopping criterion mask image to calculate the energy change stopping criterion only on voxels whose mask value is non-zero.

Energy Change Stopping Criterion
Stopping Criterion Evaluation Period


The sliding weight image ${\displaystyle w(x)}$ is defined by a function that decays with increasing distance to the input sliding surface boundary model. Specify Lambda only to use exponential decay, or specify both Lambda and Gamma to use a dirac-shaped function.

Lambda
Gamma


# Sliding Geometries Registration Use Case

Sliding geometries registration allows registration of images that show sliding motion at locally planar surfaces or at tubular structures. You can also specify point geometries where sliding is not allowed (work in progress).

Sliding Regularization Type: SlidingGeometries


The workflow is similar to the sliding organ registration use case, except:

In addition to an input sliding surface boundary model, you can also provide an sliding tube spatial object (in .tre format, in the space of the fixed image) that will determine the tubes that can slide.

Input Sliding Tube Spatial Object


In addition to the output normal vector image and sliding weight image, output a sliding geometries image, which stores the geometry type (plane, tube, point) at each voxel of the fixed image. This can be input along with the normal vector image and sliding weight image on subsequent registrations, so that they will not have to be recomputed if you decide to redo the registration with different parameters.

Output Sliding Geometries Image


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

• None

# OLDTutorial

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