A Methodology of Three-Dimensional Medical Image Registration Based on Conformal Geometric Invariant

A three-dimensional multimodality medical image registration method using geometric invariant based on conformal geometric algebra (CGA) theory is put forward for responding to challenges resulting from many free degrees and computational burdens with 3Dmedical image registration problems.Themathematical model and calculationmethod of dual-vector projection invariant are established using the distribution characteristics of point cloud data and the point-to-plane distance-based measurement in CGA space. The translation operator and geometric rotation operator during registration operation are built in Clifford algebra (CA) space. The conformal geometrical algebra is used to realize the registration of 3D CT/MR-PD medical image data based on the dual vector geometric invariant. The registration experiment results indicate that the methodology proposed in this paper is of stronger commonality, less computation burden, shorter time consumption, and intuitive geometric meaning. Both subjective evaluation and objective indicators show that the methodology proposed here is of high registration accuracy and suitable for 3D medical image registration.


Introduction
Registration and fusion of medical images in different modalities are an important task in medical image processing, as it can be convenient for doctors to achieve treatment schedules, lesion locations, disease progress determination, and treatment assessment and provide more comprehensive information for understanding the process of medical image.In recent years, with improvements of medical devices performance, imaging information gradually tends to be a character of multiresolution color and multidimensional.Intracranial 3D medical image registration technique has important clinical value.For intracranial 3D medical images (such as CT, MR, SPECT, and PET), it has notable features: firstly, intracranial soft tissue is not easily deformed as it is protected by skull; as a result, whole image data can be viewed as a rigid body; secondly, intracranial image data is provided in distribution forms of sliced layers, and a three-dimensional can be rebuilt by stacking up these lamellar images in a fixed order.Because of these two rezones described above, outer contour of different types of brain images for a patient has high degree of similarity.The similarity of multimodal intracranial images' contour profile provides excellent conditions for high-precision registration.
Three-dimensional registration technique has characteristics of complex geometric transformation, difficult global optimization, slow convergent speed, and easiness to stuck at a local optimum.Some effective methodologies have been proposed in open literatures [1].The in vitro labeling brackets fixed in head are used to process registration in the literature [2].In the literature [3], Rangarajan et al. conducted registration by using mutual information method, and Staring et al. promote the idea of "multifeature mutual information" [4].On the other hand, self-similarity is introduced into mutual information calculation by Rivaz and Louis Collins; they call it self-similarity weighted mutual information [5].Iterative closest point (ICP) algorithm is primary in the area of 3D image registration, which is based 2 Mathematical Problems in Engineering on purely geometric model and proposed in the literature [6] by Besl and McKay in 1992.Chen and Medioni also proposed a similar algorithm in the literature [7] during the same period.After that, many modified algorithms based on ICP have been proposed, which can be seen in the literatures [8][9][10][11][12].In literature [8], a nonlinear optimization factor is utilized to increase the robustness of ICP, and the modified algorithm is entitled as LM-(Levenberg-Marquardt-) ICP.Contributions of literature [9,10] are that stability of the affine transformation in ICP is successfully solved by introducing Lie group exponential.In addition, ICP is extended to registration of multiple data in the literatures [11,12].In addition to methodologies described above, there are many other 3D registration algorithms, such as a methodology based on deformation invariant attribute vector (DIAV) proposed in the literature [13], a methodology based on voxel intensity information proposed in the literature [14], a mutual information image registration methodology based on torque spindle presented in the literature [15], a methodology based on genetic algorithms and segmentation rectangle described in the literature [16], a methodology based on stochastic differential equations proposed in the literature [17], a methodology based on angular-invariant feature of 3D point cloud reported in the literature [18], and algorithms based on combined mutual information and edge correlative deviation introduced by the literatures [19,20].Most of the current registration algorithms including those described above have not yet achieved full automation and need for manual interventions.Furthermore, the rapidity of search strategy and the effectiveness of feature extraction method also need to be improved.A significant increase in the amount of 3D data registration leads to a sharp rise in conversion complexity, and existing methods have larger computing cost.
Conformal geometric algebra (CGA) was first proposed in 1878 by Clifford; it is an algebra that can contain both vector and scalar operations.The perfect theoretic system of CGA is developed by Li Hong-bo based on the theory put forward by Clifford.Currently, CGA has been successfully utilized in many fields, such as computer graphics, computer vision, and robotics, which can be seen in the literatures [21][22][23].In the Literature [24], CGA has to be combined with ICP firstly and utilized for image registration.In this paper, a 3D medical image registration method based on conformal geometric invariant is firstly presented.The mathematical model and calculation method of dual-vector projection invariant are established in CGA space.The rotation operator in CA space is put forward and registration of 3D CT/mr-PD medical image data is realized.Experimental results show that the method presented here have characteristics of simple model, high execution efficiency, more intuitive geometric meaning, better registration performance, and not being easy to fall into the local extreme value.

Image Registration Conformal Position
In traditional Euclidean geometric algebra, there are three basis vectors, denoted here by " 1 , " " 2 , " and " 3 , " respectively.However, five-dimensional conformal geometric algebra consists of traditional Euclidean geometric algebra and two newly added empty vectors, that is, " 0 " and " ∞ " that mean three-dimensional origin and point at infinity, respectively.Two empty vectors satisfy the following formula: Inner product of  0 and  ∞ is Geometric product of  0 and  ∞ is Five-dimensional conformal geometric algebra described above can represent geometric entities easily, that is, points, lines, surfaces, balls, and so forth.

Conformal Geometric Algebra Expression of Ball and Point.
In the five-dimensional conformal geometric space, a point denoted here by "" can be expressed as where "" is certain point in the three-dimensional Euclidean space.Supposed that three coordinates of point  are, respectively, "2", "−3, " and "1, " then they can be expressed as ( 5) in the five-dimensional space of conformal geometric algebra: Equation ( 5) is a space transformation from threedimensional Euclidean space to five-dimensional conformal geometric space and can also be expressed as below: The space transformation described above can also be named affine point representation.
Sphere is most basic geometry in the conformal geometric algebra, through which other geometries can easily be extended.A sphere  with center  and radius  can be expressed as below: It can be easily found that point  in the five-dimensional conformal geometry space is actually a sphere whose radius is 0.

The Conformal Geometric Algebra Expression of Plane.
Mathematical expression of CGA for a plane  is where  is a distance from origin to plane  and  is unit normal vector of European geometry in this plane and can be denoted by  =  1 +  2 +  3 .

The Construction of Conformal Geometric Invariant
Suppose that there is a rigid body "" whose centroid is denoted by "." In the five-dimensional conformal geometric algebra space, it can be expressed as The denotation "" is supposed to be a free plane through centroid ", " which can be expressed as In formula (10), "" is a distance from origin to the plane, and "" is a unit normal vector of European geometry in this plane and can be denoted by " =  1 + 2 + 3 ." In the CGA, inner product can be used as measure of distance.Therefore, on the premise of centroid coincide with origin, the point "V  " in the rigid body "" and free plane "" are processed by inner product operation.Then measurement of distance between "V  " and "" is where V  ⋅  can be expanded in detail as where   ,   , and   are triaxial values corresponding to the point  because the centroid of rigid body coincides with the origin; the free plane  that passes through the centroid also passes through the origin.At this condition, the distance from free-plane  to the origin is 0; that is, the value of  is 0. Equation ( 12) can be further simplified as below: The inner product operation between point cloud and normal vector "" of plane "" can be gotten: The expansion of ‖V  ⋅ ‖ 2 is satisfied by the following formula: Equation ( 11) transform into: From ( 15), ( 11) can be expressed as below: where It can be seen that  is a symmetric matrix, and  is a as unit vector; that is,    = 1.Now the problem is transformed to find the minimum solution of three variables (, , and ) regarding  in objective function.Lagrangian function is utilized to resolve this problem, and Lagrangian function  can be designed as follows according to (17): Then the partial derivatives of these four unknown parameters (, , , and ) are set to zero, respectively.The gradient function of  is satisfied by the following formula: Further it is simplified as Then, there are If we rewrite it into a matrix form: extracting the left matrix: it can be seen that vector  = (, , )  is satisfied by the definition formula of eigenvector of matrix .The  is the eigenvalue corresponding to .
And any eigenvector of matrix  is a unit vector.This is the mandatory condition that  is a unit vector.Therefore, if there is a nonzero eigenvalue and eigenvector (experimental data will prove this point later), then eigenvector of  and the corresponding eigenvalue of − is the solution of (22).The eigenvector corresponding to the minimum value of objective function is usually the first column vector of eigenvector matrices in the Eigenvector matrix of matrix  that we have calculated.This column vector is the unit normal vector that corresponding the free-plane .

Construction of Rotation Operator Based on the Clifford Algebra
The geometric transformation of three-dimensional Euclidean space can be expressed in the Clifford space.The corresponding Clifford algebra space is expressed as " 3 ".The geometric transformation of -dimensional Euclidean space can also be expressed in the Clifford space.The corresponding Clifford algebra space is expressed as "  ." Two   based vector spaces and two vectors ,  ∈   are given.Two vectors of "" and "" could be multiple vectors.There is a rotation operator "" between vector "" and vector "." The vector "" could be coincident with vector "" by using rotation operator "." The construction of operator "" is introduced in the following.As shown in Figure 1, a given plane denoted by "" can be denoted by bivector "" in Clifford algebra space.The vector "" can be coincident with the vector "" after rotating symmetrically around the plane of ""; namely,  =  (−, , ) . ( The vector "" in formula ( 25) is a unit vector of plane , and  is a function which does geometrical product of vectors.
The vectors "" and "" shown in Figure 2 both belong to   .The vector "" is coincident with the vector "" by rotating "" degrees.The rotating process of vector "" could be operated by using the symmetry operation method of formula (25) twice successively.Then we can obtain the geometric rotation operator "." For the sake of simplicity, "" and "" have been normalized to a unit vector.The angle bisect vector between "" and "" is denoted by "V." The vector "" can be coincident with the vector "ℎ" after rotating symmetrically around the plane of " 1 ." According to the formula (25), the vector "ℎ" can be calculated: Then a plane  2 could be constructed perpendicular to vector ", " and the vector "" is the normal vector of the plane  2 .It is easy to know that vector "ℎ" is also perpendicular to the plane  2 .The vector "" and vector "ℎ" are symmetrically around the plane  2 .According to the principle of formula (25), we can obtain the vector "": The following formula can be deducted from formulas (26) and ( 27): In formula (28), "" is called rotation operator; namely, The operator " † " is the inversion of  which is equal to "(V, )" and is satisfied by the following formula: (30)

Registration and Experimental Results
Detailed steps of registration methodology presented in this paper can be described as follows: Problems in Engineering 5 (2) calculating centroids   and   of the two point gatherings constructed in the above step: where  is number of point clouds in the discrete point gatherings; (3) coinciding the origin and centroid of reference image and floating image respectively, the minimum distance value from point cloud data to dual-vector projection invariant can be calculated, and we can obtain the normal vectors   1 and   1 of a dual-vector projection invariant; (4) constructing a rotation operator  1 that is suitable for the floating mode point gathering =1 through geometric rotation algorithm construction methodology: where   To verify the methodology proposed in this paper, clinical data of patients obtained from a project that is executed by Vanderbilt University and named "retrospective evaluation of image registration" are employed as an example [25].In these data, CT image in the training data is selected as a reference image, and MR PD is selected as a floating image.The data information of CT and MR PD images is shown in Table 1.
The size and resolution of reference image are different from floating image in Table 1.In order to achieve registration precision and improve the registration rate, the resolution of reference image and floating image can be changed to the new value shown in Table 2.
The centroid of the two modalities should firstly be moved to the origin before rotation operation.After calculation, moving extents of point gatherings with respect to reference modality and floating modality are, respectively, Δ = 62 1 + 71 2 + 21 3 and Δ = 68 1 + 51 2 + 26 3 .Registration parameters and rotation operators of the sample calculated by methodology proposed here are illustrated in Table 3.
Personal computer with "E3 1320V2" CPU and 8 GB memory is used as hardware platform.The software named "MATLAB R2012a" is used to achieve registration.The calculation time of the whole registration is 10 minutes and 21 seconds.Figures 3 and 4 are, respectively, images of CT and MR PD datasets on the cross-section of three orthogonal after visualization. Figure 5 is three views of fusion image after registration.The effect of image fusion shown in Figure 5 indicates that the registration method proposed here is of high registration accuracy and suitable for 3D medical image registration.
The assessment standard named "retrospective registration assessment standard" established by Vanderbilt University is used as "gold standard" for registration.The registration error can be calculated by using "gold standard." The screws are used as mark points fixed in brain surface.The points marked by the screws are hidden before sharing the brain data with researchers.The registration methods used by researchers are unknown by evaluation experts, and the "gold standard" established by evaluation expert is also unknown by researchers.In order to use "gold standard, " the medical images are unified to the same resolution shown in Table 4. Comparing to the "gold standard, " the registration error using    the method put forward is shown in Table 5. Comparing to the method of "IM-ICP" [26], the registration error of this paper is lower and the registration precision is higher.

Conclusions
CT and MR data of internal soft tissue is not easily to be squeezed because it is protected by outer skull.This means that such figure has strong rigidity, and it is able to maintain a clear outline.The physics of rigid body is used as the basic of the new registration method proposed in this paper.The mathematical model and calculation method of dual-vector projection invariant are established using the distribution characteristics of point cloud data and the pointto-plane distance-based measurement in CGA space.The translation operator and geometric rotation operator during registration operation are built in Clifford algebra (CA) space.Experimental results show that the method presented here has characteristics of simple model, high execution efficiency, more intuitive geometric meaning, high precision, and not being easy to fall into the local extreme value.The geometric algebra used in image registration allows us to rotate the model directly, instead of calculating the three axes rotation parameters (, , and ).This means that we just need to rotate the model once instead of three times.So it has advantages of reducing complexity of calculation, achieving high precision, and processing with more efficiency.Further research is required to seek more representative geometric invariant.Moreover, improvement in the theory system of conformal geometric algebra and Clifford algebra is also required in further research.

( 6 )
the rotation operation put forward above could not reach the registration precision in the parallel direction of "  1 ." The third eigenvectors   3 and   3 , respectively, belonging to reference mode and floating mode gotten from formula (23) are used as rotation axes.The rotation operator " 3 " can be gotten by using "  3 " and "  3 ." The construction of rotation operator  3 is similar to the construction of rotation operator  1 using the steps mentioned above.The point cloud {V   }   =1 of floating image is rotated twice by successively using rotation operations  1 and  3 , and the point cloud {V #  }   =1 is finally gotten.

Figure 5 :
Figure 5: Three views of fusion image after registration.

Table 1 :
List of 3D medical image data information in Vanderbilt University.

Table 2 :
List of 3D medical image data information after changing scale.

Table 3 :
List of registration parameters and rotation operators during registration processing.

Table 4 :
List of 3D medical image data information unified to the same resolution.

Table 5 :
Contrasting of registration error.