An ITK Implementation of Physics-based Non-rigid Registration Method

,

A clinically practical non-rigid registration method should take the following factors into account: speed, robustness and accuracy.The registration should be done within a short time to provide timely responses to the surgeons.The registration results should not be susceptible to image intensity inhomogeneity and artifacts.Moreover, the registration results should realistically reflect the physical deformation of the tissue.
In this paper, we present an ITK implementation of a physics-based non-rigid registration method [1].This method relies on: 1. a sparse displacement field and a Finite Element biomechanical model to estimate the entire deformation field.The unknown deformation field is represented by a parameterized piece-wise linear polynomials, whose parameters are estimated by approximating the scattered displacements, essentially a scattered data approximation approach.2. The sparse registration points are obtained by locating the centers of small image blocks with rich structural information, a typical feature point selection approach.3. The displacement associated with the registration point is obtained by comparing the block surrounding the registration point with the blocks located in a predefined window, i.e., block matching.
The sparse displacement field is characterized by sparsity and outliers, which compromises the accuracy of the estimation of the entire deformation field.To deal with sparsity of the deformation filed, the parameter estimation is regularized by a biomechanical model, which is capable of describing the physical deformation based on quite few data, i.e., the boundary condition.To make the estimation robust again outliers, the parameter are estimated as a Least Trimmed Squares (LTS) regression [5].More specifically, at each iteration estimate parameters first without any outliers, then identify the points with larger error as outliers, finally remove outliers from the data and re-estimate the parameters.To reduce the approximation error, the entire deformation field is estimated by an iterative method that gradually shifts from an approximation problem (minimizing the sum of a regularization term and a data error term) towards an interpolation problem (least square minimization of the data error term).
Among the three filters or steps in the PBNRR filter, block matching is the most computationally intensive, about 30% of the total execution time for an average case of non-rigid registration of preoperative and intra-operative brain MRIs (see Table 2).The block matching algorithm [1], for each block, loops over its search window to calculate the cross-correlation.In ITK implementation, we parallelized the Block Matching based on our previous work [2,3] within ITK multithreading/GPU framework to make full use of multi-core and GPU processors available to even average computing platforms like desktops and laptops.
In this paper, we first describe the principle of the physics-based non-rigid registration method.Users are referred to [1,2,3] for more details about the sequential and parallel algorithms.Then, we present its ITK implementation covering three independent filters and one main filter which puts together all three filters.
Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382Distributed under Creative Commons Attribution License Finally, we evaluate our implementation using 5 public datasets [12] on the registration of the preoperative MRI and the intra-operative MRI.
2 Physics-Based Non-Rigid Registration Method The physics-based non-rigid registration method is based on the following three steps: feature point selection, block matching, and robust finite element solver.Feature point selection identifies the small image blocks with rich structure information (see [1]), block matching estimates the displacements of these block to produce a sparse displacement filed, and the robust solver iteratively approximates the scattered data while rejecting outliers to produce an entire deformation field.

Feature point selection
The relevance of a displacement estimated with a block matching algorithm depends on the existence of highly discriminative structures within a block.The block variance is used to measure its relevance and only select a fraction of all potential blocks based on a predefined fraction.To avoid redundancy by the overlapping of blocks (i.e., eliminate blocks which are too "close to each other), a parameter of prohibited connectivity is used.Three connectivity patterns are supported in the ITK implementation: 6-connectivity, 18-connectivity, and 26-connectivity (see Section 3.1).
To address the aperture problem [6,7], the structural tensor of the block is calculated.The structural tensor reflects the distribution of the edge detections within the block, which will be incorporated into the Finite Element solver to make the estimated node displacement to favor the reduction of the deviation along the direction orthogonal to the edge direction.To avoid finding false correspondence (e.g., the tumor resection cavity), the block selection utilizes a mask image when necessary to exclude certain portions of the image while searching for the feature points (e.g., in the case of tumor resection).

Block matching
Block Matching is a well-known technique widely used in motion coding, image processing and compression [8,9,10].Block Matching is based on the assumption that a complex non-rigid transformation can be approximated by point-wise translations of small image regions.Considering a block B(O k ) in a floating image centered in O k and a predefined search window W k in a reference image, the Block Matching algorithm consists in finding the position O m in W k that maximizes a similarity measure M, which can be: mean square difference of intensity (MSD), mutual information (MI), and normalized cross correlation (NCC). ) We use NCC as the similarity measurement; however future implementation can include MSD or other similarity measures pertinent to the application of the PBNRR filter.Figure 1 illustrates the procedures of block matching.Note that the block can be specified in both floating and reference images depending on the application.For each block detected by feature point selection, we use Block Matching to calculate the displacement of the block.The location of the block that maximizes the similarity is obtained by exhaustive search.By assembling the individual displacement vectors one can create a sparse Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382Distributed under Creative Commons Attribution License displacement field D, which will be used by the solver to estimate the unknown displacement vector associated with the mesh nodes.
Figure 1.Block matching.For a small block in the floating image, find its corresponding block in a predefined search window in the reference or fixed image, then the displacement associated with the block can be calculated.The block can be specified in both floating and reference images depending on the application.

Finite element solver
Finite element solver is used to find the unknown displacement vector associated with the mesh nodes.The energy function is defined as, The first regularization term describes the stain energy of a linear elastic biomechanical model, and the second term describes the error between the simulated displacements and the real displacements, i.e., D.
 controls the balance of these two terms.U is the unknown node displacement vector with a size of n 3 .n is the number of nodes of the mesh.K is the mesh stiffness matrix of size n n 3 3  .The building of K has been well documented in [4].H is the linear interpolation matrix of size n p 3 3  , where p is the number of the registration points.S is the matching stiffness matrix of size p p 3 3  .S is an extension to the classical diagonal stiffness matrix, taking into account the matching confidence and the local structure distribution [1], whose 3 × 3 sub-matrix k S corresponding to registration point k is defined as: where k T is tensor structure of the block surrounding the registration point, which allows us to only consider the matching direction collinear to the orientation of the intensity gradient in the block.

The equation (1) can be solved by
Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382 Distributed under Creative Commons Attribution License leading to the linear system The above approximation formulation performs well in the presence of outliers but suffers from a systematic error.Alternatively, solving the exact interpolation problem based on noisy data is not adequate.The robust solver can take advantage of both approximation and interpolation to iteratively estimate the deformation from the approximation to the interpolation while rejecting outliers.The gradual convergence to the interpolation solution is achieved through the use of an external force F added to the approximation formulation of Equation ( 4), which balances the internal mesh stress:

This force
F is computed at each iteration i to balance the mesh internal force i KU , which leads to the iterative scheme:

ITK Implementation
The ITK implementation of the physics-based non-rigid registration method contains three filters: MaskFeaturePointSelection, BlockMatchingImageFilter, and FEMScatteredDataPointSetToImageFilter, which correspond to the above mentioned three components: feature point selection, block matching and robust finite element solver, respectively.These three filters can be used independently or connected together to perform non-rigid registration.Block Matching is parallelized using ITK v4 multithreading/GPU framework, for both multi-core and GPU, to accelerate the computation and gain some speedup.The robust solver is enhanced to allow the accommodation of different geometry elements in dealing with linear elastic problems by simply providing appropriate mesh.To implement non-rigid registration and achieve ease-of-use the three filters are combined into a single registration filter, PhysicsBasedNonRigidRegistrationMethod.

Feature point selection
MaskFeaturePointSelectionFilter (see Figure 2 for chart flow and Figure 3 for inheritance diagram) generates a list of feature points selected from a masked input image.It takes an Image and a mask Image as inputs and generates a PointSet of feature points as output.The feature points are physical centers of a small image blocks with higher variance.Optionally, a structure tensor is computed and stored as a pixel value for each feature point.The following optional parameters can be set:   The following is the normal usage of the filter: After Feature point selection and Block Matching, three point sets are available: feature point set with the structure tensor as the pixel value, block matching point set with the displacement as the pixel value, and the confidence point set with the similarity value as the pixel value.These three point sets will be used by the FEMScatteredDataPointSetToImageFilter to perform approximation and compute the displacement field everywhere.

Scattered data approximation
The class itk::fem::RobustSolver implements the solver presented in [1].FEMScatteredDataPointSetToImageFilter is a wrapper of RobustSolver.FEMScatteredDataPointSetToImageFilter is used to facilitate the use of RobustSolver by converting natural inputs such as mesh and feature points into specific FEMObject, providing built-in 2D and 3D rectilinear meshes, invoking RobustSolver to resolve the solution to produce a deformed FEMObject, and converting the deformed FEMObject into a deformation field.RobustSolver takes a FEMObject as input, then iteratively approximates the data (displacement) associated with the feature points while rejecting outliers, and finally outputs a deformed FEMObject.

FEMScatteredDataPointSetToImageFilter
Figure 6 shows the flow chart and the inheritance diagram of this filter, respectively.FEMScatteredDataPointSetToImageFilter provides a built-in 2D quadrilateral and 3D hexahedron mesh if the input mesh is not available.Otherwise, just simply passes the input mesh to the converter.The natural inputs of the RobustSolver are mesh and point sets including mandatory feature points and optional confidence and tensor.itk FEM library requires a FEMObject as input.FEMScatteredDataPointSetToImageFilter converts the mesh and point sets into a FEMObject, which is undertaken by a member function:

RobustSolver
Given a 2-or 3-D scattered and noisy point set, in which each point is associated with a 2-D or 3-D displacement, RobustSolver is able to approximate the data while rejecting outliers, advance toward interpolation, and finally output a deformed FEMObject.The flow chart and inheritance diagram are described in Figure 7 and Figure 8, respectively.
RobustSolver also takes into account two optional point sets: the confidence and structural tensor.Confidence point set describes our confidence for each feature point using a value between 0 and 1 (0: not trustful, 1: completely trustful), which will make the solver behavior like a weighted Least Square.Tensor point set describes the distribution of the edge direction within a small block surrounding the feature point, which is used to avoid the aperture problem [6,7].The following codes show the typical usage of the RobustSolver.Outlier rejection as a LTS regression [5]: resolve U first, then detect outliers, remove outliers and resolve U again.The F is used to reset the strain energy to enable the mesh to be deformed further.The difference between the two parts is there is no outlier rejection in the approximation to interpolation part.
Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382Distributed under Creative Commons Attribution License Figure 8.The inheritance diagram of RobustSolver.RobustSolver supports both VNL solver and Itpack solver to resolve the linear system of equations.Compared to VNL solver, Itpacks runs faster, which is the default LS solver in RobustSolver.

Main registration filter
PhysicsBasedNonRigidRegistrationMethod is the filter, which connects the MaskFeaturePointSelection, BlockMatchingImageFilter, and FEMScatteredDataPointSetToImageFilter into one pipeline to perform non-rigid registration as shown in the following codes:

Experiments and Results
We conducted experiments on the registration between preoperative MRI and the intra-operative MRI (iMRI).The five datasets come from public cases from SPL of Harvard medical school [12].Table I  The MRI of the 5 public cases were acquired with a protocol: whole brain sagittal 3D-SPGR (slice thickness 1.3 mm, TE/TR=6/35 ms, FA=75 o , FOV=24 cm, matrix=256 x 256) [11].
In Table 2 we show the qualitative results of the PBNRR filter for the 5 cases which run in a workstation equipped with an Intel® Core™ i7 CPU 260 @ 2.80 GHz with 8 GiB of RAM.
As a measure of the registration accuracy we use the one directional Hausdorff Distance (HD) as it is implemented in the vtkHausdorffDistancePointSetFilter.The HD(1→2) before PBNRR corresponds to the error between canny points in pre-operative MRI and intra-operative MRI while the HD(1→2) after PBNRR corresponds to the error between canny points in warped pre-operative MRI and intra-operative MRI.The running time includes the time for the PBNRR filter and the time for creating and writing the warped pre-operative MRI, not including the time for generating the canny points and the calculation of the HD.In Figure 9 we present the quantitative results of the PBNRR filter for the same 5 cases we used throught this evaluation.To reproduce the results, users need to checkout the master branch and the patch with topic "A2D2PBNRR" if they use the version less than ITK4.3.

Conclusion
We present an ITK implementation of a physics-based non-rigid registration method.The three filters: MaskFeaturePointSelection, BlockMatchingImageFilter, and FEMScatteredDataPointSetToImageFilter can be used separately or combined together to do registration.MaskFeaturePointSelection identifies small blocks with rich structure information.BlockMatchingImageFilter is used to find the displacement associated with these blocks in order to produce a sparse deformation field, which is used by FEMScatteredDataPointSetToImageFilter to interpolate the entire deformation field.For each block MaskFeaturePointSelection stores the structure tensor, and BlockMatchingImageFilter stores the confidence, i.e., the cross-correlation value.Both structure tensor and confidence are incorporated into the FEMScatteredDataPointSetToImageFilter to deal with aperture problem and weight the least square formula.
To reduce the computational time of block matching, it is parallelized on multi-core and GPU, and the execution time is reduced from about 70 secs to 50 secs i. Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382Distributed under Creative Commons Attribution License

Future work
In the future we plan to provide a web-service for image-to-mesh conversion to generate over the WEB the mesh of the images.This service can maintain new functionality as we better understand the needs of the ITK community.However, for the near future we can help as many ITK users as we can with the finite element mesh generation of the images that might want to use in conjunction of this filter, since we strongly suggest users to provide an anatomically adapted mesh as input for the physics-based non-rigid registration although a default rectilinear mesh is provided inside.

Figure 4 .Figure 5 .
Figure 4. Flow chart of block matching.After the filter is created and inputs are set using SetFixedImage, SetMovingImage and SetFeaturePoints, the calculation is triggered by calling Update method.After update the method GetDisplacements returns a PointSet that contains coordinates of feature points as Point values and displacement vectors as Pixel values, GetSimilarities returns a PointSet that contains coordinates of feature points as Point values and similarity values as Pixel values.

Figure 6
Figure 6.The flow chart (left) and the inheritance diagram (right) of FEMScatteredDataPointSetToImageFilter.FEMScatteredDataPointSetToImageFilter takes mesh, feature points, confidence and structural tensor as inputs.Converter first converts these inputs into a FEMObject, and then invokes RobustSolver to produce a deformed FEMObject.This deformed Object is converted into the deformation filed by DeformationFildGenerator.

Figure 7 .
Figure 7.The flow chart of RobustSolver.RobustSolver includes two parts: outlier rejection and approximation to interpolation.Outlier rejection as a LTS regression [5]: resolve U first, then detect outliers, remove outliers and resolve U again.The F is used to reset the strain energy to enable the

Figure 9 .
Figure 9.The Qualitative results for the 5 cases of the PBNRR filter.Each row corresponds to a different case, and each column from left to right: the pre-operative MRI , the intra-operative MRI and the warped pre-operative MRI.
(see Figure4and Figure5(left) for chart flow and Figure5(right) for inheritance diagram) computes displacements of given points from one image to another.It takes fixed and moving Images as well as a PointSet of feature points as inputs.Pixel values of input point set, i.e., the structural tensors are not used in this filter.The feature points are expected to lie at least SearchRadius + BlockRadius voxels from the image boundary.This is usually achieved by using an appropriate mask during selection of feature points.The default output (0) is a PointSet with displacement vector stored as the pixel value.Additional output (1) is a PointSet containing similarities, i.e., the NCC value.The number of points in the output PointSet is equal to the number of points in the input PointSet.
NonConnectivity: defines connectivity pattern (VERTEX_CONNECTIVITY, EDGE_CONNECTIVITY or FACE_CONNECTIVITY) to a feature point.The default is VERTEX_CONNECTIVITY;  BlockRadius: radius measured in voxels over which the variance is computed, its default value is 1; Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382DistributedunderCreativeCommonsAttributionLicense SelectFractionfraction of points to select out of total eligible points, default is 0.05.Figure 2. Flow chart of feature point selection After the filter is created and inputs are set using SetInput and SetMaskImage, the calculation is triggered by calling Update method.After the Update, the method GetOutput returns a PointSet that contains coordinates of feature points as Point values and (optionally) structure tensors as Pixel values.Figure 3. Inheritance diagram of feature point selection Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382DistributedunderCreativeCommonsAttribution LicenseThe following is the normal usage of the filter:BlockMatchingImageFilterThe following optional parameters can be set:  BlockRadius: radius over which variance is computed, default is 1.  SearchRadius: radius of the search window, default is 3.Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382Distributed under Creative Commons Attribution License DoF) for the building of stiffness matrix K.After converting to FEMObject, RobustSolver is invoked to construct linear system of equations described by equation (6), resolve U of the linear system, and output a deformed FEMObject, which is used by DeformationFieldGenerator to produce a deformation field.The following codes show its typical usage.
}The material properties of the biomechanical model such as Young modulus and Poisson's ratio are specified in InitializeMaterials. InitializeNodes and InitializeElements are used to store the nodes and the elements of the mesh into containers of the FEMObject.The displacement associated with the feature points are stored as loads in the FEMObject by InitializeLoads, in which the correspondence and tensor will be stored, too, if they are provided by users.After initialization, FinalizeMesh should be invoked to Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382Distributed under Creative Commons Attribution License produce Degree of Freedoms (

Table 1 .
lists the patient information including the gender, tumor location, and histopathology.Patient information of five cases from SPL of Harvard medical school.

Table 2 .
(1,12,antitative results for the 5 cases are obtained by running PBNRR filter using 8 threads .In table3we show the running time for each component of the PBNRR filter for the case 4. The experiment has run for different number of threads(1,12, 24and GPU Quadro 6000 with 448 cores).
Latest version available at the Insight Journal link http://hdl.handle.net/10380/3382Distributed under Creative Commons Attribution License