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

Authors: Liang Hua; Kean Yu; Lijun Ding; Juping Gu; Xinsong Zhang; Tianyu Cheng; Hao Feng
Journal: Mathematical Problems in Engineering
 (2014)

---

## 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.


---

## Body

## 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 “e
1,” “e
2,” and “e
3,” respectively. However, five-dimensional conformal geometric algebra consists of traditional Euclidean geometric algebra and two newly added empty vectors, that is, “e
0” and “e
∞” that mean three-dimensional origin and point at infinity, respectively. Two empty vectors satisfy the following formula:
(1)
e
0
2
=
e
∞
2
=
0
.Inner product ofe
0 and e
∞ is
(2)
e
0
·
e
∞
=
-
1
.Geometric product ofe
0 and e
∞ is
(3)
e
0
e
∞
=
e
0
·
e
∞
+
e
0
∧
e
∞
=
-
1
+
e
0
∧
e
∞
.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 “X” can be expressed as
(4)
X
=
x
+
1
2
x
2
e
∞
+
e
0
,
where “x” is certain point in the three-dimensional Euclidean space. Supposed that three coordinates of point x are, respectively, “2”, “−3,” and “1,” then they can be expressed as (5) in the five-dimensional space of conformal geometric algebra:
(5)
X
=
2
e
1
-
3
e
2
+
1
e
3
+
14
2
e
∞
+
e
0
.Equation (5) is a space transformation from three-dimensional Euclidean space to five-dimensional conformal geometric space and can also be expressed as below:
(6)
X
=
e
∞
∧
x
.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 sphereS with center C and radius r can be expressed as below:
(7)
S
=
c
+
1
2
(
c
2
-
r
2
)
e
∞
+
e
0
.It can be easily found that pointX 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 planeP is
(8)
P
=
N
+
D
e
∞
,
where D is a distance from origin to plane P and N is unit normal vector of European geometry in this plane and can be denoted by N
=
a
e
1
+
b
e
2
+
c
e
3.


## 2.1. Conformal Geometric Algebra Expression of Ball and Point

In the five-dimensional conformal geometric space, a point denoted here by “X” can be expressed as
(4)
X
=
x
+
1
2
x
2
e
∞
+
e
0
,
where “x” is certain point in the three-dimensional Euclidean space. Supposed that three coordinates of point x are, respectively, “2”, “−3,” and “1,” then they can be expressed as (5) in the five-dimensional space of conformal geometric algebra:
(5)
X
=
2
e
1
-
3
e
2
+
1
e
3
+
14
2
e
∞
+
e
0
.Equation (5) is a space transformation from three-dimensional Euclidean space to five-dimensional conformal geometric space and can also be expressed as below:
(6)
X
=
e
∞
∧
x
.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 sphereS with center C and radius r can be expressed as below:
(7)
S
=
c
+
1
2
(
c
2
-
r
2
)
e
∞
+
e
0
.It can be easily found that pointX 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 planeP is
(8)
P
=
N
+
D
e
∞
,
where D is a distance from origin to plane P and N is unit normal vector of European geometry in this plane and can be denoted by N
=
a
e
1
+
b
e
2
+
c
e
3.


## 3. The Construction of Conformal Geometric Invariant

Suppose that there is a rigid body “R” whose centroid is denoted by “c.” In the five-dimensional conformal geometric algebra space, it can be expressed as
(9)
C
=
x
c
e
1
+
y
c
e
2
+
z
c
e
3
+
1
2
(
x
c
2
+
y
c
2
+
z
c
2
)
e
∞
+
e
0
.The denotation “P” is supposed to be a free plane through centroid “c,” which can be expressed as
(10)
P
=
n
+
d
e
∞
.In formula (10), “d” is a distance from origin to the plane, and “n” is a unit normal vector of European geometry in this plane and can be denoted by “n
=
x
e
1
+
y
e
2
+
z
e
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
i” in the rigid body “R” and free plane “P” are processed by inner product operation. Then measurement of distance between “v
i” and “P” is
(11)
min
⁡
∑
i
n
(
v
i
·
p
)
2
,
where v
i
·
p can be expanded in detail as
(12)
v
i
·
p
=
v
i
·
n
+
1
2
(
x
i
+
y
i
+
z
i
)
2
·
d
,
where x
i, y
i, and z
i are triaxial values corresponding to the point i because the centroid of rigid body coincides with the origin; the free plane P that passes through the centroid also passes through the origin. At this condition, the distance from free-plane P to the origin is 0; that is, the value of d is 0. Equation (12) can be further simplified as below:
(13)
v
i
·
p
=
v
i
·
n
.The inner product operation between point cloud and normal vector “n” of plane “p” can be gotten:
(14)
v
i
·
n
=
[
x
i
e
1
+
y
i
e
2
+
z
i
e
3
+
1
2
(
x
i
2
+
y
i
2
+
z
i
2
)
e
∞
+
e
0
]
iiiiiiii
·
(
x
e
1
+
y
e
2
+
z
e
3
)
iiiiiiii
=
x
i
x
+
y
i
y
+
z
i
z
.The expansion of∥
v
i
·
n
∥
2 is satisfied by the following formula:
(15)
∥
v
i
·
n
∥
2
=
(
x
i
x
+
y
i
y
+
z
i
z
)
2
iiiiiiiiiiii
=
x
i
2
x
2
+
y
i
2
y
2
+
z
i
2
z
2
+
2
(
x
i
y
i
)
x
y
iiiiiiiiiiiiiiii
+
2
(
x
i
z
i
)
x
z
+
2
(
y
i
z
i
)
y
z
.Equation (11) transform into:
(16)
min
⁡
∑
i
=
1
n
(
x
i
x
+
y
i
y
+
z
i
z
)
2
.From (15), (11) can be expressed as below:
(17)
min
⁡
(
n
T
M
n
)
,
where
(18)
n
T
=
(
x
,
y
,
z
)
,
M
=
∑
i
=
1
n
(
x
i
2
x
i
y
i
x
i
z
i
x
i
y
i
y
i
2
y
i
z
i
x
i
z
i
y
i
z
i
z
i
2
)
.It can be seen thatM is a symmetric matrix, and n is a as unit vector; that is, n
T
n
=
1. Now the problem is transformed to find the minimum solution of three variables (x, y, and z) regarding n in objective function. Lagrangian function is utilized to resolve this problem, and Lagrangian function L can be designed as follows according to (17):
(19)
L
=
n
T
M
n
-
λ
(
n
T
n
-
1
)
.Then the partial derivatives of these four unknown parameters (x, y, z, and λ) are set to zero, respectively. The gradient function of L is satisfied by the following formula:
(20)
∇
L
=
2
(
M
n
-
λ
n
)
=
0
.Further it is simplified as(21)
M
n
=
λ
n
.Then, there are(22)
∂
L
∂
x
=
∑
i
=
1
n
[
x
i
2
x
+
(
x
i
y
i
)
y
+
(
x
i
z
i
)
z
]
-
λ
x
=
0
,
∂
L
∂
y
=
∑
i
=
1
n
[
y
i
2
y
+
(
x
i
y
i
)
x
+
(
y
i
z
i
)
z
]
-
λ
y
=
0
,
∂
L
∂
z
=
∑
i
=
1
n
[
z
i
2
z
+
(
x
i
z
i
)
x
+
(
y
i
z
i
)
y
]
-
λ
z
=
0
,
∂
L
∂
λ
=
-
(
x
2
+
y
2
+
z
2
-
1
)
=
0
.If we rewrite it into a matrix form:(23)
∑
i
=
1
n
(
x
i
2
x
i
y
i
x
i
z
i
x
i
y
i
y
i
2
y
i
z
i
x
i
z
i
y
i
z
i
z
i
2
)
(
x
y
z
)
=
λ
(
x
y
z
)
,
extracting the left matrix:
(24)
M
=
∑
i
=
1
n
(
x
i
2
x
i
y
i
x
i
z
i
x
i
y
i
y
i
2
y
i
z
i
x
i
z
i
y
i
z
i
z
i
2
)
,
it can be seen that vector n
=
(
x
,
y
,
z
)
T is satisfied by the definition formula of eigenvector of matrix M. The λ is the eigenvalue corresponding to M.And any eigenvector of matrixM is a unit vector. This is the mandatory condition that n is a unit vector. Therefore, if there is a nonzero eigenvalue and eigenvector (experimental data will prove this point later), then eigenvector of M 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 M that we have calculated. This column vector is the unit normal vector that corresponding the free-plane P.


## 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 “G
3”. The geometric transformation of n-dimensional Euclidean space can also be expressed in the Clifford space. The corresponding Clifford algebra space is expressed as “G
n.” Two G
n based vector spaces and two vectors a, b
∈
G
n are given. Two vectors of “a” and “b” could be multiple vectors. There is a rotation operator “T” between vector “a” and vector “b.” The vector “a” could be coincident with vector “b” by using rotation operator “T.” The construction of operator “T” is introduced in the following.As shown in Figure1, a given plane denoted by “π” can be denoted by bivector “B” in Clifford algebra space. The vector “a” can be coincident with the vector “c” after rotating symmetrically around the plane of “π”; namely,
(25)
c
=
G
(
-
u
,
a
,
u
)
.Figure 1
Symmetry operation.The vector “u” in formula (25) is a unit vector of plane π, and G is a function which does geometrical product of vectors.The vectors “a” and “b” shown in Figure 2 both belong to G
n. The vector “a” is coincident with the vector “b” by rotating “θ” degrees. The rotating process of vector “a” could be operated by using the symmetry operation method of formula (25) twice successively. Then we can obtain the geometric rotation operator “T.”Figure 2
Rotation operator method.For the sake of simplicity, “a” and “b” have been normalized to a unit vector. The angle bisect vector between “a” and “b” is denoted by “v.” The vector “a” can be coincident with the vector “h” after rotating symmetrically around the plane of “π
1.” According to the formula (25), the vector “h” can be calculated:
(26)
h
=
G
(
-
v
,
a
,
v
)
.Then a planeπ
2 could be constructed perpendicular to vector “b,” and the vector “b” is the normal vector of the plane π
2. It is easy to know that vector “h” is also perpendicular to the plane π
2. The vector “b” and vector “h” are symmetrically around the plane π
2. According to the principle of formula (25), we can obtain the vector “b”:
(27)
b
=
G
(
-
b
,
h
,
b
)
.The following formula can be deducted from formulas (26) and (27):
(28)
b
=
G
[
-
b
,
G
(
-
v
,
a
,
v
)
,
b
]
ii
=
G
[
G
(
b
,
v
)
,
a
,
G
(
v
,
b
)
]
ii
=
G
(
T
,
a
,
T
†
)
.In formula (28), “T” is called rotation operator; namely,
(29)
T
=
G
(
b
,
v
)
.The operator “T
†” is the inversion of T which is equal to “G
(
v
,
b
)” and is satisfied by the following formula:
(30)
G
(
T
†
,
T
)
=
1
.


## 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{
v
i
r
}
i
=
1
n
r and {
v
i
f
}
i
=
1
n
f according to 3D modality date of reference image and floating image;
(2)
calculating centroidso
r and o
f of the two point gatherings constructed in the above step:
(31)
o
=
1
n
∑
i
=
1
n
v
i
∈
G
3
,
wheren 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 vectorsn
1
r and n
1
f of a dual-vector projection invariant;
(4)
constructing a rotation operatorT
1 that is suitable for the floating mode point gathering {
v
i
f
}
i
=
1
n
f through geometric rotation algorithm construction methodology:
(32)
T
1
=
(
n
1
r
)
*
η
1
,
where(33)
η
1
=
(
n
1
r
)
*
+
(
n
1
f
)
*
∥
(
n
1
r
)
*
+
(
n
1
f
)
*
∥
;
(5)
through rotation operation,{
v
i
f
}
i
=
1
n
f becomes a new point gathering {
v
i
f
′
}
i
=
1
n
f. The mathematical expression of this new point gathering is
(34)
v
i
f
′
=
T
1
v
i
f
T
1
*
(
i
=
1,2
,
…
,
n
f
)
;
(6)
the rotation operation put forward above could not reach the registration precision in the parallel direction of “n
1
f.” The third eigenvectors n
3
r and n
3
f, respectively, belonging to reference mode and floating mode gotten from formula (23) are used as rotation axes. The rotation operator “T
3” can be gotten by using “n
3
r” and “n
3
f.” The construction of rotation operator T
3 is similar to the construction of rotation operator T
1 using the steps mentioned above. The point cloud {
v
i
f
}
i
=
1
n
f of floating image is rotated twice by successively using rotation operations T
1 and T
3, and the point cloud {
v
i
f
#
}
i
=
1
n
f 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.Table 1
List of 3D medical image data information in Vanderbilt University.
Mode
The size of 3D image data
The spatial resolution (mm)
The total number of cloud points
x
y
z
x
y
z
CT
512
512
29
0.6535
0.6535
4.0000
164280
MR_PD
256
256
26
1.2500
1.2500
4.0000
177767The size and resolution of reference image are different from floating image in Table1. 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.Table 2
List of 3D medical image data information after changing scale.
Mode
The scale of 3D image data
The spatial resolution (mm)
x
y
z
x
y
z
CT
128
128
45
2.614
2.614
2.614
MR_PD
128
128
45
2.5
2.5
2.5The 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,Δ
r
=
62
e
1
+
71
e
2
+
21
e
3 and Δ
f
=
68
e
1
+
51
e
2
+
26
e
3. Registration parameters and rotation operators of the sample calculated by methodology proposed here are illustrated in Table 3.Table 3
List of registration parameters and rotation operators during registration processing.
Mode
Parameter
Parameter values
Rotation operators
Reference mode
n
1
r
0.0060
e
1
-
0.0103
e
2
+
0.9999
e
3
T
1
=
-
5.1608
×
1
0
-
5
e
12
+
0.00149
e
23
-
0.00413
e
31
(
n
1
r
)
*
0.0060
e
1
+
0.0103
e
2
-
0.9999
e
3
n
3
r
0.0097
e
1
-
0.9999
e
2
-
0.0104
e
3
(
n
3
r
)
*
0.0097
e
1
+
0.9999
e
2
+
0.0104
e
3
Floating mode
n
1
f
-
0.0023
e
1
-
0.0133
e
2
+
0.9999
e
3
T
3
=
-
0.06538
e
12
+
0.00159
e
23
-
0.00070
e
31
(
n
1
f
)
*
0.0023
e
1
+
0.0133
e
2
-
0.9999
e
3
n
3
f
-
0.1209
e
1
-
0.9926
e
2
-
0.0135
e
3
(
n
3
f
)
*
0.1209
e
1
+
0.9926
e
2
+
0.0135
e
3Personal 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. Figures3 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.Three views of CT image.
(a)
Cross-section ofx-y
(b)
Cross-section ofx-z
(c)
Cross-section ofy-zThree views of MR_PD image.
(a)
Cross-section ofx-y
(b)
Cross-section ofx-z
(c)
Cross-section ofy-zThree views of fusion image after registration.
(a)
Cross-section ofX-Y
(b)
Cross-section ofX-Z
(c)
Cross-section ofY-ZThe 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 Table4. 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.Table 4
List of 3D medical image data information unified to the same resolution.
Mode
The scale of 3D image data
The spatial resolution (mm)
x
y
z
x
y
z
CT
334
334
116
1
1
1
MR_PD
320
320
104
1
1
1Table 5
Contrasting of registration error.
Registration error (mm)
Minimum error
Maximum error
Average error
Method in this paper
1.74
5.46
3.00
IM-ICP method
1.29
6.58
4.22


## 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 (x, y, and z). 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.



---

*Source: 573702-2014-08-26.xml*