Elastic Image Registration Based on Domain Decomposition with Mesh Adaptation

Medical images are increasingly used within healthcare for diagnosis, planning treatment, and monitoring disease progression. The images acquired at different times, with different imaging modalities, from different subjects etc. often provide an additional clinical information that is not revealed in the separate images. The spatial relation between the images has to be found and this process is called image registration. In our contribution, we use elastic registration which assumes that the images are two different observations of an elastic body which is discretized by the finite element method. We are especially interested in the problems where the requirements on the registration prevent the application of standard FFT based solvers to the solution of auxiliary linear problems, which is the case when the part of the two observations can be related by a rigid body motion. Because the medical images usually contain a large area of background and a small area of changes, a regular discretization results in waste of computational resources due to the fine refinement of the space outside the region of interest (especially in 3D). To avoid this, we use coarser grid with local refinement that takes into account specific features of the images and their differences. The related elasticity problems are solved by TFETI, which is a variant of the Finite Element Tearing and Interconnecting (FETI) domain decomposition method for massively parallel numerical solution of elliptic Partial Differential Equations (PDE) with optimal complexity.


Introduction
Image registration is one of the challenging problems in image processing.Since the information gained from two (or more) images is usually of a complementary nature, proper integration of useful data obtained from the separate images is often desired.Thus there is a need to establish an exact point-to-point correspondence between the voxels in one image and those in the other image.Although manual alignment of images is possible, it is time-consuming (especially in more than two dimensions) and lacks reproducibility.In practice, the specific type of geometric transformation as well as the meaning of the correspondence depends on the specific application.
The medical image registration is a vital component of a large number of applications, see [1], [2] and [3].Such applications include the alignment of data sets from different modalities, comparison of followup scans to a base-line scan, alignment of pre-and post-contrast images, updating treatment plans for radiotherapy and surgery, monitoring of diseases, atlasbased segmentation, and creating models of anatomy.
The early applications of the medical image registration considered registering the brain images of the same subject acquired with different modalities (e.g.MRI and CT or PET).For these applications a rigid body approximation was sufficient as there is a relatively little change in brain shape or position within the skull over the relatively short periods between scans.Today rigid registration is often extended to include affine registration, which includes scale factors and shears, and can partially correct for calibration differences across scanners or gross differences in scale between subjects.Clearly most of the human body does not conform to a rigid or even an affine approximation and much of  the most interesting and challenging work in registration today involves the development of non-rigid registration techniques.Besides elastic registration we can mention TPS (thin-plate spline) [3], free form deformation [4], optical flow techniques [5], and fluid, diffusion or curvature registration [3].
One of the images is usually viewed as a reference (target, fixed, baseline) R and the other one as a deformable template (source, moving, floating) T .The optimal transformation ϕ is estimated by minimizing the cost function, called the distance measure (similarity measure or metric) D, which determines how much is the image R, in a certain sense, similar to the image T .So the problem can be formulated as A direct minimization of the distance measure has its drawbacks: a solution is not necessarily unique and it actually may not exist.Thus the problem (Eq.( 1)) is ill-posed.Moreover additional implicit assumptions on the transformation can arise, for example in medical images no additional cracks or folding of the tissue are allowed (the transformation should be diffeomorphic).Both situations can be solved by adding regularizer.Problem (Eq.( 1)) is a challenging optimization problem, in particular in 3D, whose solution requires effective solution of the associated elasticity problems.In our paper, we propose to use a variant of the FETI do-main decomposition method [6], which is widely used in the field of material sciences and structural mechanics.
The performance of the FETI method can be improved by effective mesh generation.In the case of images of human body, e.g.CT scans, there are many heterogeneous parts of the image domain in which it is not necessary to have fine mesh.Mesh adaptation is clearly recognized as an efficient and powerful method for improving the accuracy of the solution as well as for capturing the physical phenomena behavior.Reducing the number of Degrees of Freedoms (DoF) allows to substantially reduce the CPU time.For example, Haber at al. [7] used octrees as a basic structure for the underlying displacement field which requires a special and careful discretization of the variational form.
Using the adaptive mesh generation based on geometric attributes leads to the combination of coarse parts of the mesh with very fine parts identified by doctors or by some automatic prediction.The paper is focused on using an a priori estimate of the image based on geometric qualities to create an adaptive mesh for the image registration.Our adaptive mesh has much smaller number of DoF for similar accuracy of the solution.Further image registration is done on solid mesh.Though the reduction is essential in 2D (single image frame), it is even more relevant in 3D.To simplify the exposition, the algorithms are explained in 2D, but our approach can be applied to 3D registration.Detailed information on numerical experiments with human body scans can be seen in Sec. 5.

Elastic Image Registration
The elastic potential of deformation in the image registration has been introduced by Broit [8].His method can register images with local non-rigid geometric differences.Broit uses it to automatically find an optimal mapping between a CT image and an atlas of brain anatomy.The external forces have been derived by correlating the intensity-based properties in local regions in the source and target image.
Elastic registration is based on physical motivation that the images are two different observations of an elastic body with the template image corresponding the body in the reference configuration.The transformation ϕ : R 2 → R 2 of the image is split into the identity part and the displacement u : As the regularizer we use the linearized elastic potential where λ and µ are the Lame constants.This linear model assumes small deformations.For larger deformations it can be replaced by a viscous fluid model [9].
We gain regularized problem, which is a starting point for a stable numerical implementation (parameter α control the strength of the smoothness of the displacement versus the similarity of the images, in our case is included in Lame constant µ) where A necessary condition for a minimizer u of problem (Eq.( 4)) is that the Gateaux-derivative dJ [u; v] of J vanishes for all siutable perturbations v.
The partial differential operator associated with the Gateaux-derivative of the elastic potential is the wellknown Navier-Lame operator.The displacement of the elastic body is then obtained by the solution of partial differential equation The regularizer has the meaning of internal forces which implicitly enforce the displacement to be smooth.
The distance measure D represents external forces, so it pushes the deformable template into the direction of the reference.We choose the Sum of Squared Differences (SSD) where The images are represented by the compactly supported mappings R(x), T (x) : Ω → R, where Ω = (0, 1) 2 .Hence, T (x) and R(x) denotes the intensities of images at the spatial position x.We set R(x) = 0 and T (x) = 0 for all x / ∈ Ω.
If the images are monomodal, i.e., the intensities of corresponding pixels are supposed to be identical, the SSD is a reasonable measure for applications.But if the images are multimodal, then the SSD can fail.Alternative choices are the mutual information [10] or the normalized gradient field [11].
Modersitzki [3] solves this equation by the finite difference approximation with periodic boundary conditions, so that the resulting matrix is highly structured.This structure allows him to use the Fast Fourier transform to diagonalize and invert the matrix.Thus he gains very effective solution with complexity O(N log N ).But this approach doesn't work with the constraints, which destroy the structure of the matrix.
We discretize the problem using a finite element method with piece-wise affine basis functions on triangular elements.To approximate the image gradient, we use a convolution with an appropriate kernel of the Sobel operator.
Let us mention some examples for future use of elastic registration.Non-rigid image registration was already used for the automatic quantification of small changes in the volume of the anatomical structures in brain over time by means of the segmentation propagation for dementia progression [12], for monitoring of the response of the brain to drug therapy [13].Another example represents modeling organs motion during the respiratory cycle in radiotherapy [14] or the study thorax and lung deformation to develop tools to evaluate the reproducibility of ABC (Active Breath Control) [15] .

Adaptive Mesh Generation
The local Adaptive Mesh Refinement (AMR) problem uses the local insertion of additional vertices in order to produce a mesh which provides higher accuracy of the solution.We place more nodes to the areas where we anticipate a large local error of the solution or require higher precision for the final solution.We have no difficulties with boundary condition issues because the area near the boundary is the image background which doesn't require adaptation.Different methods exist to perform the mesh adaptation, see [16] or [17].
Our approach locally refines the mesh size depending of the geometry of the image.
Because we are working with the finite element solvers based on the meshes of triangles, we have to partition the selected triangles and handle the transition between two zones with different level of refinement.The conformity of the mesh is preserved by using a particular templates of the elements in the transition area, see [18].By conformity we understand that the intersection of adjacent triangles is either a common vertex or common edge.Good quality of the mesh is preserved during the refinement process by using a method that does not refine a single angle more than once.
Figure 2 shows 2D image meshes generated by the adaptive refinement algorithm described further.The presence of the refined parts of the image corresponds to the differences between images, in our case to the tissue movement.A refinement algorithm is controlled by a local error indicator.There are many methods that could be used to define the refinement criterion, such as those based on the geometry of the domain or a priori or a posteriori error estimators for a particular solver or manual (expert) selection [7] and [19].In our case, we are using the error indicator resulting from very simple analysis.The error indicator, the binary image/logical matrix I, is assembled by thresholding the absolute difference between the images R and T using the given intensity threshold t, see Fig. 3.We use the image erosion and dilation to fill pixel-sized holes to improve compactness.The matrix I reads where τ ij denote the coordinates of vertex j of element τ i and values R(τ ij ) and T (τ ij ) are intensities (gray values) in given coordinates.The mesh is constructed in such way that the nodes are identical with pixels.We use linear interpolation whenever the values in elements are needed.The error function is defined as If the intersection of the element and the region of interest given by the binary image is not empty, then the error function is not zero and the element has to be refined.The AMR algorithm, see Alg. 1, begins with a given level of refinement L, threshold t, and with a coarse initial mesh.The finest possible level of the refinement corresponds to the image size.All the initial elements belong to the level l = 0.The first step of the algorithm is to assembly binary image based on the error indicator to define the criteria of refinement.In the next step, for each level of refinement, we mark all the elements identified by a refinement criteria as the error elements (ones to be refined).To avoid splitting of any angle in the initial mesh more than once, all elements sharing a node with error elements are marked as nodal elements.This step ensures that the resulting mesh will keep similar quality as the initial mesh and is illustrated in Fig. 4.Each element which shares an edge with a nodal element and is not an error or nodal one is marked as transition element.The marked elements are then divided into sub-elements that belong to level l +1 by particular pattern.Standard pattern is used for error and nodal elements.For the transition el-ements, there are used special transition patterns, see Fig. 6.Finally, all unrefined elements are copied to level l + 1.At the finest level, the rest of the elements are added.

Algorithm 1 AMR algorithm
Require: Coarse mesh, threshold t, level of refinement L. Ensure: Final mesh of level L.

4:
Mark elements with error E(τ i ) > 0 as error elements,

5:
Mark (nodal) elements sharing a node with an error element,

6:
Select elements sharing edge with nodal elements as transition elements,

8:
Refine error and nodal elements by standard pattern,

9:
Refine transition elements by particular transition pattern, 10: Compose mesh for level l + 1 as non-refined elements from level l and elements created by pattern subdivision.11: end for We are not using a template with two edges to the refinement.Such a refinement could rapidly decrease the quality of the mesh.During the marking step of the algorithm, we also mark all the edges of elements that would be refined as edges to refine.At this point we need to investigate if there is any element which has two edges to refine (f alse element).We mark such element as an element to refine (we are using the same tag as for nodal elements) by marking the third edge as one to refine and repeat this process until there are only elements with zero, one or all three edges to refine.This procedure is called edge propagation and is illustrated in Fig. 6.We use the mid-edge subdivision for refining the elements of the mesh.Each error and nodal element is split into four homothetic elements.The transition element in 2D could be only element with a single hanging node.In this case we connect the pending node to the opposite vertex, which splits the element into two ones as in Fig. 5.The algorithm works similarly in 3D.

Problem Solution
For the solution of the discretized problem we use a variant of FETI method which is one of the most successful methods for the parallel solution of elliptic partial differential equations.This method is based on the decomposition of the spatial domain into nonoverlapping subdomains that are "glued" by Lagrange multipliers.Thus, after eliminating the primal variables, the original problem is reduced to a small, relatively well conditioned, typically equality constrained quadratic programming problem that is solved iteratively.Its drawback is difficulty to determine the kernels of stiffness matrices reliably in the presence of rounding errors.
TFETI method is easier to implement and preserves efficiency of the coarse grid of the classical FETI [20].The basic idea of TFETI is to simplify the inversion of stiffness matrices of subdomains by using Lagrange multipliers not only for gluing of the subdomains along the auxiliary interfaces, but also to enforce the Dirichlet boundary conditions.Thus, all the subdomains are floating and their stiffness matrices will have a priori known kernels-bases of rigid body motion.To apply the FETI based domain decomposition, we partition Ω into N subdomains Ω s as in Fig. 7 (now λ detones vector of Lagrange multipliers) and we denote by K s , f s and u s for s = 1, . . ., N , the subdomain stiffness matrix, the subdomain force and displacement vectors respectively.
We shall get the discretized problem where The matrix B and the vector c enforce the prescribed displacements on the part of the boundary with the imposed Dirichlet condition and the continuity of the displacements across the auxiliary interfaces.For easy implementation we choose the Dirichlet condition on whole boundary.We shall introduce the Lagrangian associated with problem (Eq.( 10)) by and the equivalent saddle-point problem By substituting where d := BK † f − c, e := −R T f , and K † denotes a generalized inverse matrix satisfying KK † K = K.As mentioned above the kernels R s of the local stiffness matrices K s are bases of the rigid body motion and can be formed directly where (x i , y i ), i = 1, ..., n s are the coordinates of the nodes of Ω s .Using R s , we can easily assemble the block diagonal basis R of the kernel of K as After denoting with relatively small and well conditioned block F as compared with K. We can obtain λ by solving where P = I − Q and Q = G T (GG T ) −1 G denote the orthogonal projectors on the image space of G T and on the kernel of G, respectively.This problem may be solved effectively by the conjugate gradient method.More details can be found in [20].

Numerical Experiments
Algorithm for creating the adaptive image-based mesh described in this paper and assembling matrices for TFETI was implemented in Matlab.We used the implementation of the TFETI method in FLLOP (FETI Light Layer On top of PETSc), the package for constrained quadratic programming and FETI domain decomposition, see [21].We used the Matlab Mesh Partitioning and Graph Separator Toolbox to divide the mesh into equal-sized components with a small number of connecting nodes as in Fig. 8.This toolbox contains Matlab code for several graph and mesh partitioning methods.Specifically, we used a simple and quite efficient geometric method [22].
The values of Lame's constants, λ and µ, are artificial and there is no general rule how to choose them.λ is usually set to zero to ensure that deformations are only effected in the directions of applied forces.Our choice µ = 1.5 was obtained experimentally.For testing, we used medical images from department of Oncology at the University Hospital of Ostrava and the numerical experiments were, for now, performed on personal computer with two processors.The resulting images T u can be seen in Fig. 9 and the resulting optimal transformation as displacement field in Fig. 10.In Tab. 1, we report the computational times for preprocessing image with creating adaptive mesh, time for assembling matrices for TFETI, and time for performing TFETI method for different resolutions of images in Fig. 1.We can observe that   the adaptive meshes cost much less time than the full meshes.

Conclusion
We have implemented the image registration algorithm based on the TFETI method with the local adaptive mesh generation in 2D.We use the SSD distance measure for deriving external forces from the intensities of monomodal image.The efficiency of the algorithm was demonstrated by numerical experiments.In the next stage we will implement the algorithm in 3D and exploit the massively parallel capability of TFETI on the supercomputer cluster Anselm.Using local adaptive mesh combining with the domain decomposition method could provide effective strategy for future computing with huge 3D meshes.

Fig. 1 :
Fig. 1: Registration example.Images in this paper are from Department of Oncology at the University Hospital of Ostrava.The images show the cuts in the chest inhale and exhale obtained by CT method.

Fig. 3 :
Fig. 3: Binary image based on image differences and intensity threshold.
(a) One level of local refinement.(b) Finer level of refinement.

Fig. 5 :
Fig. 5: Refinement patterns of the error and the transition elements in 2D and 3D.

Fig. 6 :
Fig. 6: Refinement patterns of the error and the transition elements in 2D and 3D.
Tab. 1: Type denotes type of mesh -full and regural or adaptive image-based mesh, Ns number of subdomains, N nod number of nodes (nodes on interfaces are doubled), T i computational time for preprocessing image with creating adaptive mesh, Tm time for assembling matrix for TFETI and Tt time for performing TFETI method.