TubeTK/Documentation/Sliding Organ Registration
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 is a function of the current estimation of the deformation field . 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:
- is the organ boundary in the vicinity of ,
- is the vector within the deformation field at location
- is the gradient of the -th component of
- is the component of in the normal direction
- 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 is close to 1:
- 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
- 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 and is implemented in itkImageToImageDiffusiveDeformableRegistrationFilter and itkImageToImageDiffusiveDeformableRegistrationFunction using ITK's finite differences framework:
where is the canonical unit vector, i.e
Further away from organ boundaries, 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.
- 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 , 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 ,
- 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 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.
- 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
- Load the input data
- 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
- 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.
- 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)
- Perform the registration
- Click "Apply" - the registration will take some time
- 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'
- 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
- Models must be flipped along the z-axis before the registration will work correctly.
- Output deformation and normal vector images not correctly visualized in Slicer4
- 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
- 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
- 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
- Data assimilation using a gradient descent method for estimation of intraoperative brain deformation
- Songbai J Alex Hartov, David Roberts and Keith Paulsen
- Medical Image Analysis, 15(5), p 744-756, 2009
- Biomechanical models that simulate brain deformation are gaining attention as alternatives for brain shift compensation. One approach, known as the “forced-displacement method”, constrains the model to exactly match the measured data through boundary condition (BC) assignment. Although it improves model estimates and is computationally attractive, the method generates fictitious forces and may be ill-advised due to measurement uncertainty. Previously, we have shown that by assimilating intraoperatively acquired brain displacements in an inversion scheme, the Representer algorithm (REP) is able to maintain stress-free BCs and improve model estimates by 33% over those without data guidance in a controlled environment. However, REP is computationally efficient only when a few data points are used for model guidance because its costs scale linearly in the number of data points assimilated, thereby limiting its utility (and accuracy) in clinical settings. In this paper, we present a steepest gradient descent algorithm (SGD) whose computational complexity scales nearly invariantly with the number of measurements assimilated by iteratively adjusting the forcing conditions to minimize the difference between measured and model-estimated displacements (model-data misfit). Solutions of full linear systems of equations are achieved with a parallelized direct solver on a shared-memory, eight-processor Linux cluster. We summarize the error contributions from the entire process of model-updated image registration compensation and we show that SGD is able to attain model estimates comparable to or better than those obtained with REP, capturing about 74% to 82% of tumor displacement, but with a computational effort that is significantly less (a factor of 4-fold or more reduction relative to REP) and nearly invariant to the amount of sparse data involved when the number of points assimilated is large. Based on five patient cases, an average computational cost of approximately 2 minutes for estimating whole-brain deformation has been achieved with SGD using 100 sparse data points, suggesting the new algorithm is sufficiently fast with adequate accuracy for routine use in the operating room (OR).