Abstract
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 3D medical image registration problems. The mathematical model and calculation method 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.
1. 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 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–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–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.
2. Image Registration Conformal Position
In traditional Euclidean geometric algebra, there are three basis vectors, denoted here by “,” “,” and “,” respectively. However, five-dimensional conformal geometric algebra consists of traditional Euclidean geometric algebra and two newly added empty vectors, that is, “” and “” that mean three-dimensional origin and point at infinity, respectively. Two empty vectors satisfy the following formula:
Inner product of and is
Geometric product of and is
Five-dimensional conformal geometric algebra described above can represent geometric entities easily, that is, points, lines, surfaces, balls, and so forth.
2.1. 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 three-dimensional 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.
2.2. 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 .
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 “.” In the CGA, inner product can be used as measure of distance. Therefore, on the premise of centroid coincide with origin, the point “” in the rigid body “” and free plane “” are processed by inner product operation. Then measurement of distance between “” and “” is where 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 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, . 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 .
4. 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 “”. 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 “.” The vector “” can be coincident with the vector “” after rotating symmetrically around the plane of “.” According to the formula (25), the vector “” can be calculated:
Then a plane could be constructed perpendicular to vector “,” and the vector “” is the normal vector of the plane . It is easy to know that vector “” is also perpendicular to the plane . The vector “” and vector “” are symmetrically around the plane . 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 “” and is satisfied by the following formula:
5. Registration and Experimental Results
Detailed steps of registration methodology presented in this paper can be described as follows:(1)constructing two discrete point gatherings and according to 3D modality date of reference image and floating image;(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 and of a dual-vector projection invariant;(4) constructing a rotation operator that is suitable for the floating mode point gathering through geometric rotation algorithm construction methodology: where (5) through rotation operation, becomes a new point gathering . The mathematical expression of this new point gathering is (6)the rotation operation put forward above could not reach the registration precision in the parallel direction of “.” The third eigenvectors and , respectively, belonging to reference mode and floating mode gotten from formula (23) are used as rotation axes. The rotation operator “” can be gotten by using “” and “.” The construction of rotation operator is similar to the construction of rotation operator using the steps mentioned above. The point cloud of floating image is rotated twice by successively using rotation operations and , and the point cloud is finally gotten.
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, and . 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.
| (a) Cross-section of - |
| (b) Cross-section of - |
| (c) Cross-section of - |
| (a) Cross-section of - |
| (b) Cross-section of - |
| (c) Cross-section of - |
| (a) Cross-section of - |
| (b) Cross-section of - |
| (c) Cross-section of - |
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.
6. 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 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. 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (61273024, 61305031), Natural Science Fund of Jiangsu Province and “333” Fund (12KJB510023, KB2012227), and PAPD, “226” High-level Personnel Training Project of Nantong City.