- Email: [email protected]
Automatic characteristic-calibrated registration (ACC-REG): Hippocampal surface registration using eigen-graphs
Automatic Characteristic-Calibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs
Journal Pre-proof
Automatic Characteristic-Calibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs Hei Long CHAN, Tsz Chun YAM, Lok Ming LUI PII: DOI: Reference:
S0031-3203(19)30443-1 https://doi.org/10.1016/j.patcog.2019.107142 PR 107142
To appear in:
Pattern Recognition
Received date: Revised date: Accepted date:
17 May 2019 5 November 2019 27 November 2019
Please cite this article as: Hei Long CHAN, Tsz Chun YAM, Lok Ming LUI, Automatic CharacteristicCalibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs, Pattern Recognition (2019), doi: https://doi.org/10.1016/j.patcog.2019.107142
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
HIGHLIGHT • Automatic intrinsic feature extraction on hippocampal mesh surface • Feature calibration by Eigen-graphs to improve consistency and accuracy of the feature correspondences on different subjects • Fast and accurate surface registration is available in just 5 seconds with minimal RMSE • Quantitative comparison with other existing methods demonstrates the effectiveness and efficiency of the proposed model on real data from the ADNI database
1
Automatic Characteristic-Calibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs Hei Long CHANa , Tsz Chun YAMa , Lok Ming LUIb,∗ a
222B, Lady Shaw Building, The Chinese University of Hong Kong, Shatin, Hong Kong 207, Lady Shaw Building, The Chinese University of Hong Kong, Shatin, Hong Kong
b
Abstract In this paper, we propose an efficient algorithm, the ACC-REG, to automatically extract intrinsic key characteristics on hippocampal mesh surfaces and hence compute an accurate registration mapping between them. Given a pair of hippocampal surface mesh, the proposed algorithm constructs the eigengraphs, an intrinsic feature on the surface, on each surface as its representative. The eigen-graphs are then calibrated along the longitudinal direction of the hippocampal surfaces. Accurately corresponded intrinsic characteristics on each hippocampus can thus be extracted. As a result, the two surfaces can be registered with improved accuracy and low computation cost. Experiments on ADNI data demonstrate the effectiveness of the proposed ACC-REG model over existing methods. Keywords: Feature Extraction, Feature Correspondence, Surface Registration, Surface deformation, Eigen-graph
∗
Corresponding author Email addresses: [email protected] (Hei Long CHAN), [email protected] (Tsz Chun YAM), [email protected] (Lok Ming LUI)
Preprint submitted to Pattern Recognition
January 4, 2020
1. Introduction
Surface based shape analysis has a wide range of applications. For instance, it is a key tool in medical science for disease classification and further diagnosis. With a reliable registration process, MR/CT images of different patients scanned at various time can be accurately corresponded and then compared, providing a promising environment for further analysis. One example is the registration between hippocampal surfaces which is very important to tracing the progression of the Alzheimer’s disease. While there are existing models particularly designed for hippocampal surface registration, accuracy is still hindered by different difficulties. On one hand, there are not much obvious and natural geometric features on the hippocampal surfaces that can be utilized to define landmark correspondences. On the other hand, scanning error due to hardware is still common and puts great challenges in extracting well-defined features on the hippocampal surfaces. In addition, large deformations of the hippocampus are expected in case of disease, which further ruins the stability and accuracy in the registration process. To address the above problems, a robust and automatic landmark based registration algorithm, ACC-REG, is proposed in this paper for accurate hippocampal surface registration. The intrinsic Laplace-Beltrami operator is utilized to generate feature correspondences on the hippocampus with high stability. Given two hippocampal surface meshes, the Helmholtz differential equation involving the Laplace-Beltrami operator is solved and hence an eigen-graph, an intrinsic feature on the hippocampal surface, is constructed on each of them. Inspired by the histogram matching technique from the 2
study of statistics, the eigen-graphs are then calibrated along the longitudinal direction on the hippocampus so that the accuracy of the correspondences between the two hippocampus are greatly improved. Then, well-corresponded intrinsic characteristic curves on the hippocampus are extracted. Finally, a fast and accurate surface registration algorithm is applied to compute the correspondence over the whole surface. The contribution of this paper mainly lies in three aspects. Firstly, the proposed characteristic curves are closely related to the unique geometry of the hippocampus. So the landmark correspondences are meaningful and robust. Besides, the extraction of the characteristic curves involves only a direct and linear process and hence efficiency is guaranteed. Finally, calibrating the eigen-graphs greatly reduces the influences of global and local deformations on the hippocampal surfaces, which in turn greatly improves the accuracy of the registration between the hippocampal surfaces, as being verified by various experiments. 2. Related Works
Registration between anatomical structures has a wide range of real applications, for example, disease classification [1, 2], morphometry analysis [3, 4], and surface reconstruction [5], etc. The research field is mainly divided into two main streams, namely, 2D and 3D registrations. For 2D registration, scanned images of tissues are one-to-one corresponded according to the depth of the scanning. In [6], Zhou et. al. proposed a novel technique for prealignment in both mono-modality and multi-modality image 3
registration based on statistical correlation of gradient information. In [7], Zhang et. al. proposed a multi-modality medical image registration model using support vector machines and incorporating with prior joint intensity distribution between the two imaging modality based on registered training images. In [8], Ma et. al. proposed a robust feature-guided Gaussian mixture model for image matching which can handle both rigid and non-rigid image transformation robustly. In [9], Ma et. al. proposed a feature mismatching removal scheme which removes mismatches from given putative image feature correspondences by maintaining the local neighborhood structures of those potential true matches. The proposed scheme can also be applied to medical imaging to seek for reliable correspondences between images of the same tissue. However, we note a general problem for utilizing 2D registration in medical analysis. Namely, the slicing of the tissue into 2D images is solely based on the scanning hardware and the gesture of the subject. Regarding to the fact that the position and the shape of the target tissue may have great variation among subjects, images of the same tissue at the same depth across subjects may not necessarily have a robust and meaningful point-wise correspondence. Comparing to 2D registration, 3D registration accounts more for geometric information to improve the robustness and accuracy of the registration. One commonly used registration model for spherical topology surface is the spherical harmonics (SPHARM) algorithm. For example, in [10], DV Vranic et. al. built a database of surface query and designed a SPHARM-based algorithm to retrieve similar objects using the surface database. In [11], Michael Kazhdan et. al. used the spherical harmonics to represent 3D shapes to solve
4
the rotation invariant shape matching problem. Later, other models are designed for surface registration. In [12], Gu et. al. proposed an algorithm for global conformal parametrizations based on the structure of the cohomology group of the holomorphic one-forms for surfaces. In particular, for genus zero surfaces, the proposed algorithm finds an unique mapping between any two manifolds by minimizing the harmonic energy of the mapping, which can then be applied to the cortical surface matching problem. In [13], Wang et. al. developed an algorithm which extends the mutual information method to automatically match general 3D surfaces to automate the matching of surface features. However, both [12, 13] are built in the absence of manual landmark correspondences on the surfaces, which does not suit our problem setting. In [14], Ma et. al. solves the problem of nonrigid point set registration by designing a robust transformation learning scheme. The principle is to iteratively establish point correspondences and learn the nonrigid transformation between two given sets of points. In [15], Ma et. al. presented a new point matching algorithm for robust nonrigid registration which iteratively recovers the point correspondence and estimates the transformation between two point sets. Their proposed model is general and has much applications, including medical analysis. In [16], Cheng-Tiao Hsieh presented a model to utilize the open source Point Cloud library (PCL) to develop a series of computational registration processes efficiently and robustly. The model is a key process of the Digital Face-Inspection (DFI) system and has lots of medical applications. In [17], Ilwoo Lyu et. al. designed a cortical surface registration model using then hierarchical spherical deformation which achieves group-wise surface correspondence with no template selection
5
bias. Recently, surface registration has also been incorporated with neural networks, for example, in [18]. It is also worth mentioning that high-dimensional registration is also under study. In [19], Lee et. al. proposed a new method to obtain registration between n-dimensional manifolds with very large deformations. Highdimensional registration also has medical applications to register between volumetric data, for example, in [20]. Finally, we would like to review some literature which lead to our proposed model. In [21], Wong et. al. proposed a novel approach for extracting two intrinsic feature curves on hippocampal surfaces. The model is used to study the geometric changes of hippocampus in case of the Alzheimer’s disease (AD). In [22, 23], Lui et. al. proposed a model to generate optimized conformal parametrizations of brain surfaces which match prior landmark curves exactly with shape-based correspondences between them. The idea is further extended to general surfaces in [23, 24]. These works provide a hint to our proposed model on defining and utilizing meaningful landmark correspondences on the hippocampal surfaces. In [25], Choi et. al. proposed a fast and efficient landmark-based surface registration model for genus zero surfaces. The proposed model computes the optimized spherical harmonic parametrization with consistent landmark alignment. It is not just accurate but also speedy as their model has a totally linear implementation. 3. The proposed model
The details of the proposed ACC-REG model is explained in this sec6
tion. Registration is a process to extract features on each subject, find the correspondence over the feature sets, and finally estimate and construct the surface correspondence based on the feature sets correspondence. In view of this, intrinsic features provide much stability and reliability over extrinsic features. Therefore, in this work, the eigen-graph is utilized as the (intrinsic) features for landmark correspondences. Given a hippocampal surface H, consider the Helmholtz equation on H: − ∆f (x) = λf (x),
(1)
where ∆ denotes the Laplace-Beltrami operator on H. The common geometry of hippocampus ensures the eigenvalues of H always form a non-negative increasing sequence 0 ≤ λ1 ≤ λ2 ≤ · · · → +∞. Therefore, one can define an eigen-graph on H as follows: Definition 1 (The Eigen Graph). The eigen graph of a hippocampal surface H is defined to be the eigenvector FH of the Helmholtz equation (1) such that the associated eigenvalue λ satisfies: λ > 0,
ˆ for all other eigenvalues λ ˆ λ≤λ
(2)
Figure (1) shows an example of the eigen-graph FH on H. The eigengraph FH stores much topological and geometric information of H. For instance, FH always has exactly 2 critical points on the two tips of H. Besides, the values of FH increases from one critical point to another critical point. This gives a parametrization of H as follows:
7
Figure 1: Example of the eigen-graph FH on a hippocampal surface H, function values is normalized from blue (0) to red (1)
Theorem 1. Let H be a hippocampal surface and FH : H → [0, 1] be the eigen-graph on H. Let p1 , p2 ∈ H be the two critical points of FH such that
FH (p1 ) < FH (p2 ). For any 0 < < 12 (FH (p2 ) − FH (p1 )), FH−1 ([FH (p1 ) + , FH (p2 ) − ])
(3)
is diffeomorphic to I × S 1 , where I = [FH (p1 ) + , FH (p2 ) − ] and S 1 is a circle. From the theorem, there exists a smooth mapping from I × S 1 to FH−1 (I), where I ⊂ FH (H \ {p1 , p2 }).
Hence, the hippocampal surfaces can be
parametrized by longitudinal and azimuthal coordinates. In view of this, Wong et. al. [21] proposed to extract two intrinsic feature curves γ1 , γ2 as landmark correspondences by minimizing the following energy: 1 {γ1 , γ2 } = argmax b−a
Z
a
b
|γ1 (t) − γ2 (t)|dt − c2
Z
a
b
|γ˙1 (t)|dt − c1
Z
a
b
|γ˙2 (t)|dt, (4)
where c1 , c2 ≥ 0 are the weighting parameters and [a, b] = f (H) is the range of the un-normalized first nontrivial eigenfunction on H. The motivation of the definition (4) is that the feature curves should be separated as much 8
(a) HP from normal people
(b) HP from AD patient
Figure 2: Examples of hippocampal surfaces: (left) From normal subject; (right) From patient
as possible while each of them should also be as smooth as possible. In Wong et. al.’s work, after extracting γ1 and γ2 , the hippocampal surfaces are parametrized onto the domain K = {(α, θ) ∈ [0, 1] × [0, 2π)} by minimizing the energy: Φ=
arg min Φ(Γ1 )={θ=0},Φ(Γ2 )={θ=π}
Z
H
||∇Φ||2 .
(5)
This is called the eigen-harmonic parametrization (EHP) of H. However, there are two major drawbacks of Wong et. al.’s model. Firstly, the implementation is inefficient as it involves gradient descents of two energy functional, with each of them usually takes over thousands of iterations to converge numerically. Secondly, it is realized that the progression of the eigen-graph FH on H from one tip to another varies significantly among different subjects. In fact, in real application, it is found that the geometry of the hippocampus from different subject can have large variations due to illness and defect (see figure (2)). In such situation, feature correspondences may be inaccurate and cause error in the registration. Therefore, we propose the following calibration of the eigen-graphs to improve the accuracy of the eigen-graph correspondences.
9
3.1. Calibration of the Eigen-Graph
Suppose FH1 and FH2 are eigen-graphs corresponding to two hippocampal surfaces H1 and H2 respectively. Using the longitudinal parametrization of H1 and H2 as described in Theorem (1) and normalizing the parametrization domain as [0, 1], one has: Theorem 2. For all t ∈ [0, 1], there exists an unique z ∈ R+ such that FHj (Hj (t)) = z on the slide Hj (t). Therefore, FHj can be re-defined as functions on [0, 1] for all j = 1, 2. Then, the cumulative distribution function FH∗ j of FHj can be defined by: FH∗ j (x)
=
Z
x
FHj (t)dt,
(6)
0
for all x ∈ [0, 1] and j = 1, 2.
Now, for each x1 ∈ [0, 1], there exists x2 ∈ [0, 1] such that FH∗ 1 (x1 ) =
FH∗ 2 (x2 ). By the monotonicity of FH∗ 1 and FH∗ 2 , there exists:
Definition 2 (Histogram Matching Function). Suppose that FH∗ 1 , FH∗ 2 : [0, 1] → R are two monotone increasing functions. The histogram matching
function G12 : [0, 1] → [0, 1] for FH∗ 1 , FH∗ 2 is defined by G12 (x1 ) = x2 ,
(7)
for all x2 ∈ [0, 1] such that FH∗ 1 (x1 ) = FH∗ 2 (x2 ), and for all x1 ∈ [0, 1]. Using the histogram matching function G12 , the eigen-graphs FH1 and FH2 can be calibrated as follows: 10
Definition 3 (Calibrated Eigen-Graph). Let H1 , H2 be two hippocampal surfaces, FH1 , FH2 be their eigen-graphs, FH∗ 1 , FH∗ 2 be their cumulative eigengraphs, and G12 : [0, 1] → [0, 1] be the histogram matching function matching FH∗ 1 to FH∗ 2 . The calibrated eigen-graph of H1 is defined to be FHH12 = FH1 ◦ G12 .
(a) FH1 on H1
(b) FH2 on H2
(8)
(c) FHH12 on H1
Figure 3: Example of calibrating the eigen-graphs on two hippocampal surfaces H1 and H2 .
The calibrated FHH12 can be viewed as a synchronized version of FH1 with respect to FH2 . Figure (3) demonstrates an example of the calibration process and figure (4) shows the histograms of FH1 and FH2 and the comparison between them. It can be seen that after the calibration, FHH12 has almost the same distribution as FH2 . Hence, any features correspondences generated using the calibrated eigen-graphs are more coherence and accurate. In particular, intrinsic characteristic curves can be extracted to provide wellmatched features for accurate registration. 3.2. Extraction of Intrinsic Characteristic Curves from Eigen-Graphs The extraction of feature characteristic curves mainly involves two parts. In the first part, two feature curves mutually opposite to each other are 11
(a) Histogram of FH1
(b) Histogram of FH2
(c) Comparing FH1 and FH2
(d) Histogram of FHH12
Figure 4: Demonstration and comparison of the histograms of FH1 and FH2 . Each graph is a histogram plot, with the y-axis being the number of vertices having the particular value in the eigen-graph indicated by the x-axis.
extracted as the main features on the hippocampus. Secondly, extra assisting feature curves are extracted based on the main feature curves to help improving the registration accuracy. To begin with, we extract two fundamental intrinsic characteristic curves: Definition 4 (Fundamental Intrinsic Characteristic Curves). Let H be a hippocampal surface and FH be its eigen-graph. The fundamental intrinsic characteristic curves on H are defined to be the curves γ1 , γ2 : [0, 1] → H such that {γ1 , γ2 } = argmax
Z
1
0
12
|γ1 (t) − γ2 (t)|2 dt.
(9)
By the geometry of hippocampal surfaces, one curve (γ1 ) must lie on the most curvy region along the longitudinal direction with the other one (γ2 ) on its opposite side. However, solving (9) directly requires high computation cost and may result in zigzag curves depending on the mesh triangulation. Therefore, the following notion of discrete curve on mesh is employed: Definition 5 (Discrete Curves on Surfaces). Let H be a hippocampal surface mesh or point cloud, that is, H = {Vi }N i=1 for some N ∈ N. A discrete curve γ on H is defined by γ = {vi }ni=1 ,
(10)
where {vi }ni=1 ⊂ {Vi }N i=1 and the connectivity of γ is given by v1 → v2 → v3 → · · · → vn .
(11)
Under this definition, after normalizing FH to be ranged [0, 1], γ1 , γ2 are extracted by γ1 = {vi }ni=1 , γ2 = {wi }ni=1 , where {vi , wi } = arg max |v − w|,
(12)
Ωi ()
where Ωi () = FH−1 ([ ni − , ni + ]) is the subset of H at the query level-set
1 FH−1 ( ni ) with bandwidth ∈ (0, 2n ). Figure (5)(a),(b) shows an example of
the two (discrete) intrinsic characteristic curves γ1 and γ2 on a hippocampal surface. 13
(a) γ1
(b) γ2
(c) γ1f it
(d) γ2f it
(e) γ1proj
(f) γ2proj
Figure 5: Examples of characteristic curves extracted on a hippocampal surface H. Vertices on discrete curves are connected sequentially for visualization.
Note that γ1 , γ2 are not smooth enough in general. To remedy the situation, the following polynomial fitting calibration scheme is used. Firstly, the polynomial fitting algorithm is applied on γj to get γjf it : [0, 1] → R3 , where γjf it (t) = (xj (t), yj (t), zj (t)) ∀t ∈ [0, 1], Z 1 f it γj = argmin |γj (t) − γjf it (t)|2 dt,
(13) (14)
0
for all j = 1, 2. From the examples shown in figure (14), γjf it are far smoother than γj .
14
However, there is a drawback that γjf it sits in R3 indeed and may not be lie on H. Therefore, a projection of γjf it onto H is applied to get γjproj , a smooth intrinsic characteristic curves lying on H, which is defined by: Definition 6 (Projected Smooth Intrinsic Characteristic Curves). Let H be a hippocampal surface, FH be its eigen-graph, and γ1f it , γ2f it be its two polynomial fitted intrinsic characteristic curves of degree k. The projected smooth intrinsic characteristic curves of degree k are defined to be the curves γjproj : [0, 1] → H such that γjproj
= argmin
Z
0
1
|γjproj (t) − γjf it (t)|2 dt,
(15)
for all j = 1, 2. Numerically, γ1f it and γ2f it partitions the curve domain [0, 1] into 0 < l1 < l2 < · · · < lm < 1,
(16)
where 1 lk = (FH (vik ) + FH (wik )), 2
k = 1, 2, . . . , m.
(17)
In order to capture the projected curve γjproj on H, we re-partition 1 1 2 n−1 [0, 1] = [0, ] ∪ [ , ] ∪ · · · ∪ [ , 1] n n n n
(18)
into equally spaced sub-intervals. Then, the projected curve γjproj is extracted by 0
0
{vi , wi } = arg min(||v − γ1f it ||0 + ||w − γ2f it ||0 ) v,w∈Ωi ()
15
(19)
Since Ωi () is a finite set, γ1proj and γ2proj can be extracted directly by solving (19) using the exhaustive search. Figure (5)(e),(f) demonstrates a pair of projected curves on H. Comparing figure (5)(e),(f) with figure (5)(c),(d), the projected curves are indeed less smooth but lie within H. Also, comparing figure (5)(e),(f) and figure (5)(a),(b), the projected curves are still smoother than the fundamental curves. The two projected curves γ1proj and γ2proj are selected as the (intrinsic) characteristic curves for accurate surface registration. In order to enhance the robustness of the registration, some more feature curves lying on the centre part of H are also extracted as characteristic curves, apart from γ1proj , γ2proj . Given the hippocampal surface H and the two projected smooth intrinsic characteristic curves γ1proj and γ2proj , an azimuthal coordinate function θ : H → [0, 2π) can be defined such that θ=0
on γ1proj ,
θ=π
γ2proj .
on
(20)
In view of this, more characteristic curves can be defined on H, as inspired by [21], by: Definition 7 (Assisting Characteristic Curves). Let H be a hippocampal surface, FH be its eigen-graph, γ1proj , γ2proj be its two projected smooth intrinsic characteristic curves of degree k, and θ be an azimuthal coordinate function on H. The assisting characteristic curves are defined to be the curves Γk : [0, 1] → H where: Γk = θ−1 (
k π), M
k = 1, 2, . . . , 2M − 1. 16
(21)
In discrete case, for each region Ωi (), by applying the principal component analysis, Ωi () can be projected onto the plane P (Ωi ()) spanned by its first two component vector. Note that Ωi () intersects with γ1proj , γ2proj at 0
0
vi , wi respectively. Let
1 0 0 di = (vi + wi ) 2
(22)
be the mid-point of the two vertices and ui be the unit vector at di , perpendicular to the P (Ωi ()), pointing to the direction for which FH is increasing. 0
0
For each vertex z ∈ Ωi (), let zproj , (vi )proj , (wi )proj , dproj and uproj be the i i 0
0
projection of z, vi , wi , di and ui onto the plane P (Ωi ()) respectively. Then, by computing −1 <
σ = cos
0
(zproj − dproj ), ((wi )proj − dproj )> i i , 0 proj proj |zproj − di | · |(wi )proj − di |
(23)
the azimuthal function θ : H → [0, 2π) can be assigned to z depending on 0
, zproj − dproj , uproj } by the orientation of Sk = {(wi )proj − dproj i i i
θ(z) =
σ
if Sk is positively-oriented,
2π − σ
(24)
otherwise.
Therefore, each Γk in (21) can be extracted by Γk = {zik }ni=1 where θ(zik )
k = arg min θ(z) − π M z∈Ωi ()
(25)
for each i = 1, 2, . . . , n. Note that by definition, Γ0 = γ1proj and ΓM = γ2proj .
17
Figure 6: Example of all the characteristic curves extracted on a hippocampal surface (black curves indicate γ1proj , γ2proj , green curves indicate Γk ’s)
Therefore, the set of all characteristic curves is: C = {Γ0 = γ1proj , Γ1 , Γ2 , . . . , ΓM = γ2proj , . . . , Γ2M −1 , }.
(26)
Figure (6) demonstrates all the characteristic curves on a hippocampal surface. 3.3. Surface Registration using Calibrated Characteristic Curves
Once the characteristic curves are computed, they can be utilized as landmark correspondences for the surface registration. To promote efficiency and accuracy, the FLASH algorithm [25] is employed in this step. The FLASH algorithm is a landmark-based registration algorithm designed particularly to register between surfaces having spherical topology. The FLASH algorithm mainly consists of two steps. In the first step, a spherical conformal parametrization function Ψi is computed for each Hi . For each of H1 , H2 , a parametrization ψj : Hi \ {Tj } → C is computed by
18
𝐻1
𝑃𝑁−1
𝜓1
𝑃𝑆
𝑃𝑆−1
𝑔
𝑃𝑁 1) 𝑅(𝐻
𝑃𝑆−1
𝑔𝐷
𝑃𝑆
𝑃𝑁−1
𝑓∘𝜑
(Landmark)
𝐻2
𝑃𝑁−1
𝜓2
𝑃𝑆
𝑔
(Landmark)
𝑃𝑆−1
Figure 7: Illustration of the pipeline of the FLASH algorithm
solving the linear system: P
[u,v]∈Hj
kuv (ψj (u) − ψj (v)) = 0
ψi (tj ) = bl
if u 6= tj1 , tj2 , tj3 ,
(27)
if l = 1, 2, 3;
l
where kuv is the coefficient of the discrete Laplace-Beltrami operator, Tj = [tj1 , tjj , tj3 ] is a triangle of the mesh Hj and bl is a preassigned coefficient. In order to minimize the angle distortion of the parametrization, ψj is further composed with stereographic projections PN , PS , the north pole projection and the south pole projection respectively, and also a conformality distortion minimization function g which can also be computed by solving a linear 19
system. Finally, the spherical conformal parametrization is given by Ψi = PS−1 ◦ g ◦ Ps ◦ PN−1 ◦ ψj
(28)
for each j = 1, 2. After obtaining the parametrization Ψj for each hippocampal surface Hj , let Bj = Ψj (Hj ) be the ball in R3 . By directly solving for the Mobius P n 0 0 2 transformation f (v) = av+b on B1 to minimize 2M i=1 |f (vi )−wi | , landmark f1 = f (B1 ), then B f1 , B2 are projected onto mismatching is eliminated. Write B
C by PN and then the landmark constrained harmonic map ϕ : C → C can be obtained by solving the linear system P 0 kuv (ϕ(u) − ϕ(v)) + λ(ϕ(u) − wi ) = 0 [u,v]∈H 1 P [u,v]∈H1 kuv (ϕ(u) − ϕ(v)) = 0 ϕ(tjl ) = tjl
0
if u = vi , i = 1, 2, . . . , 2M n, 0
0
if u ∈ / {v1 , . . . , vn , tj1 , tj2 , tj3 }, if l = 1, 2, 3, (29)
where λ > 0 is a weighting parameter. Afterwards, a few more stereographic projections and a quasi-conformal mapping is applied to minimize the angle ˆ : H1 → B2 is computed by distortion. Then, the pre-registration mapping R ˆ = P −1 ◦ gD ◦ PS ◦ P −1 ◦ f ◦ ϕ ◦ PN ◦ Ψ1 . R S N
(30)
Finally, the registration from H1 to H2 can be obtained by conventional interpolation methods, e.g. scatter interpolation function I : R3 → H2 . ˆ is the desired registration. Then, R := I ◦ R
Figure (7) demonstrates the pipeline of the FLASH algorithm. For the 20
detailed implementation, readers are referred to [25]. There are several advantage on using FLASH algorithm on registration between hippocampal surfaces. Firstly, the FLASH algorithm involves only solving sparse linear systems and hence efficiency of the algorithm is guaranteed, whereas in [21] it may take over thousands of iterations for the gradient descent for the harmonic energy minimization to converge numerically. Secondly, using the FLASH algorithm does not only minimizes the angle distortion as in [21] but also enforces exact landmark correspondences. Therefore, the power of the calibrated characteristic curves as landmark correspondences can be fully utilized. 3.4. Implementation of the ACC-REG algorithm
The pipeline of the proposed ACC-REG algorithm is concluded in the Algorithm 1 table for reference. 4. Experiments
In this section, experimental results are recorded to evaluate the proposed ACC-REG algorithm, in terms of improving the registration accuracy and also the computation cost. For this purpose, 45 hippocampal surfaces {Hj }45 j=1 from the ADNI database are extracted. Each surface mesh consists of approximately 9000 vertices, and is re-orientated to be centered at origin and resized so that max(||vj ||2 ) = 1. We perform several experiments. For every experiments, each of the 45 surfaces are subjected to a particular manual deformation which simulates abnormal deformations in case of disease. 21
Algorithm 1: ACC-REG Input: Hippocampal surfaces H1 , H2 , Parameters n, M, Output: Surface Registration Functions R : H1 → H2 1 Compute the Eigen-Graph FH1 , FH2 by solving (1) on H1 , H2 respectively; 2 Compute the histogram matching function G12 by (6,7); H1 3 Compose the Calibrated Eigen-Graph FH = FH1 ◦ G12 ; 2 4 for H1 , H2 do 5 Compute the Fundamental Intrinsic Characteristic Curves γ1 = {vi }ni=1 , γ2 = {wi }ni=1 by solving (12); 6 Apply polynomial fitting to obtain γ1f it , γ2f it on R3 by (14); 7 Compute the Projected Smooth Intrinsic Characteristic Curves 0 0 γ1proj = {vi }ni=1 , γ2proj = {wi }ni=1 by (19); 8 Compute the Assisting Characteristic Curves Γ1 , Γ2 , . . . , Γ2M −1 by (25); ˆ : H1 → R3 by the FLASH 9 Compute the pre-registration Function R algorithm using the two sets of {Γk }’s on H1 , H2 as landmark correspondences; ˆ 10 Compose the registration mapping R = I ◦ R. Different registration algorithms are applied to register the original surfaces to the deformed surfaces, including the spherical harmonics (SPHARM) algorithm, the MP-RPM algorithm [14], the EHP-based algorithm [21] and the proposed ACC-REG algorithm. Quantitative comparisons are done among the algorithms in terms of RMSE and also computation time. The eigengraphs generated from the proposed ACC-REG algorithm and from the original EHP-based algorithm [21] are also compared. Throughout all the experiments, the parameter setting for the proposed ACC-REG algorithm is: n = 10, M = 3, = 0.03. All the parameter setting for the FLASH algorithm follows that in [25]. All the parameter setting for
22
other algorithms follow that listed in the original publications. In particular, we use all the spherical harmonics up to the 45-th degree for the SPHARM algorithm. The registration results are visualized in the figures in this section. In each plot, the colors on the surface indicate the value of the eigen-graph, from blue (0) to red (1). The black rings on the surface is the level-set of the eigen-graph for reference. Experiment 1: In the first experiment, abnormal elongation of the hippocampus is simulated. For each hippocampus Hj , let Ω ⊂ Hj denotes the area for which FH (Ω) ≥ 0.7. By applying the PCA, the first component d is taken as the deformation direction. Each vertex vj ∈ Ω is then deformed
by vj∗ = vj + (FH (vj ) − 0.7)d. Figure (8)(a),(b) demonstrates an example of
such deformation. Figure (8)(c),(d),(e) shows the registration results of the MP-PRM algorithm [25], the EHP-based algorithm [21] and the proposed ACC-REG algorithm respectively. it can be seen that the distribution of the eigen-graph on the resultant surface is the most natural using the proposed algorithm. The statistics of the RMSE using different algorithms are recorded in table (1). Note that the mean RMSE between the original surfaces and the deformed surfaces is recorded in the row “Control” for reference to the scale of the deformation. Experiment 2: In the second experiment, a point diffusion deformation is applied to Ω = FH (Ω) ≥ 0.7 to simulate inflammation of the hippocampus. The simulation is as follows. The mean vertex v¯ = meanvj is computed over Ω. Then, the unit vector nj from v¯ to vj is computed for each vj ∈ Ω.
Finally, the deformed surface is defined by vj∗ = vj + 0.5nj . Same as the first
23
(a) Original surface
(b) Deformed surface
(c) Result of SPHARM algorithm
(d) Result of MP-PRM (e) Result of EHP-based (f) Result of proposed algorithm ACC-REG algorithm algorithm Figure 8: Sample registration results of the first experiment
experiment, different registration algorithms are used to register between the original surfaces and the deformed surfaces. One example of such deformation and the registration results by different algorithms is plotted in figure (9). The numerical results, i.e. the RMSE, for each algorithm is recorded in the table (1). Experiment 3: In the third experiment, the deformation is designed to simulate tumors on the hippocampal surfaces. For each surface Hj , on each region Ωi () as defined in (12), a random vertex is picked and a tumor is built centered at this vertex with height 0.1 along the outward normal direction. An example of such deformation and the visualization of the registration results due to different algorithms, and the statistical result of the RMSE of different algorithms are recorded in the figure (10), and the table (1) respectively. The above three experiments totally consists of 45 × 3 = 135 subjects 24
(a) Original surface
(b) Deformed surface
(c) Result of SPHARM algorithm
(d) Result of MP-PRM (e) Result of EHP-based (f) Result of proposed algorithm algorithm ACC-REG algorithm Figure 9: Sample registration results of the second experiment
which should be enough to draw a rigid conclusion. Referring to the table (1), it is evident that the proposed ACC-REG algorithm out-performs all other registration algorithms in terms of RMSE. This can be explained by the effectiveness of the calibrated eigen-graphs and hence the intrinsic feature correspondences extracted from the eigen-graphs. The high accuracy in the feature correspondences favors a lot to the registration process. Moreover, the parameters M, n in the ACC-REG algorithm can be altered to adjust the number of the feature curves used and the number of query vertices involved in each feature curve. This helps the ACC-REG algorithm to adapt well to different database with consistent accuracy. In view of the table (2), the proposed ACC-REG algorithm also has a merit of fast computation. To further investigate the computation efficiency of the proposed algorithm, it is recorded that the mean computation time for computing the eigen-graph is 0.27 seconds. The mean computation time 25
(a) Original surface
(b) Deformed surface
(c) Result of SPHARM algorithm
(d) Result of MP-PRM (e) Result of EHP-based (f) Result of proposed ACC-REG algorithm algorithm algorithm Figure 10: Sample registration results of the third experiment
for the landmark correspondences (i.e. Γk ’s) is 0.80 seconds and that for the FLASH algorithm is 3.12 seconds. In particular, the reason for the fast computation of the landmark correspondences is that, the dimension of the linear systems involved is greatly reduced by partitioning the surface H into Ωi ()’s and find the landmark vertices on each Ωi () one-by-one. Also, the FLASH algorithm itself is also known for the fast computation speed. Combining these together, the proposed ACC-REG algorithm has an advantage on the computation speed over other algorithms. Comparatively, the accuracy of the SPHARM model is limited by the maximum degree of spherical harmonics used. Higher degree of spherical harmonics store much more local information and hence gives more accurate registration result, but the computation cost significantly increases in using more spherical harmonics. As a result, the model cannot achieve a fast and accurate registration. 26
Algorithm Control SPHARM MR-PRM EHP-based ACC-REG
Experiment 1 0.1170 0.0527 0.0433 0.0474 0.0333
Experiment 2 0.0862 0.0754 0.0323 0.0524 0.0234
Experiment 3 0.1000 0.0523 0.0351 0.0482 0.0330
Table 1: Mean RMSE of the registration results in the experiments of different algorithms.
As for the MR-PRM algorithm, it is noted that its accuracy is highly comparable with our proposed ACC-REG model. However, the computation complexity of the MR-PRM algorithm is quite high. Although the authors has designed an alternative way to slightly reduce the computation complexity (see [14]), the computation time is still significantly higher than the proposed ACC-REG model. This puts a great challenge in dealing with batch-wise tasks. The EHP-based model lacks the calibration of the eigen-graph. This makes the landmark correspondences inaccurate. Besides, using arc length as the only regularization to the feature curve makes no reference to the geometry of the hippocampal surface. This further hinders the accuracy in the correspondences of the feature curves. The proposed ACC-REG model has a higher accuracy then the EHP-based model because it incorporates more geometric information in the model by calibrating the eigen-graphs and the feature curves by histogram matching and polynomial fitting respectively. Finally, the registration process of the ACC-REG model is faster than the EHP-based model mainly due to the application of the FLASH algorithm.
27
Algorithm Computation time
SPHARM 49.49s
MR-PRM 46.29s
EHP-based 8.15s
ACC-REG 4.19s
Table 2: Mean computation time of different algorithms.
5. Conclusion and Future Work
In this paper, an advanced surface registration algorithm called the ACCREG algorithm is proposed. The ACC-REG algorithm uses eigen-graphs to understand the geometry of the given hippocampal surfaces. Firstly, the eigen-graphs are mutually calibrated to synchronize their progressions along the corresponding surfaces. Secondly, on each surface, intrinsic characteristic curves are extracted and are calibrated using a projected polynomial fitting approach. Assisting feature curves are then extracted to further incorporate with more geometric information. Taking these characteristic curves as landmark correspondences, a fast and efficient registration model is utilized to register the given surfaces accurately. Through experiments, it is clear that ACC-REG helps to improve the accuracy of the registration process over other algorithms including the SPHARM algorithm, the MR-PRM algorithm and an EHP-based model. The ACC-REG algorithm also has a merit of fast computation, in which two surface meshes consisted of roughly 9, 000 vertices can be registered accurately within 5 seconds. In the future, further applications of the proposed ACC-REG algorithm including disease classification, registration of high genus surfaces, object simulation, etc., will be investigated.
28
Acknowledgement This work is supported by HKRGC GRF (Project ID: 14303414). Data used in preparation of this article were obtained from the Alzheimers Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found here: ADNI Acknowledgement List Reference References [1] Chan, Hei Long, Hangfan Li, and Lok Ming Lui, Quasi-conformal statistical shape analysis of hippocampal surfaces for Alzheimer s disease analysis, Neurocomputing 175 (2016): 177-187. [2] Gerardin, Emilie, et. al., Multidimensional classification of hippocampal shape features discriminates Alzheimer’s disease and mild cognitive impairment from normal aging, Neuroimage 47.4 (2009): 1476-1486. [3] Meintjes, E. M., et. al., A tensor-based morphometry analysis of regional differences in brain volume in relation to prenatal alcohol exposure, NeuroImage: Clinical 5 (2014): 152-160. [4] Liu, Ya-Shu, Han-Bing Yan, and Ralph R. Martin. As-rigid-as-possible surface morphing. Journal of computer science and technology 26.3 (2011): 548-557. 29
[5] Reuter, Martin, Franz-Erich Wolter, and Niklas Peinecke. LaplaceBeltrami spectra as “Shape-DNA” of surfaces and solids. Computer-Aided Design 38.4 (2006): 342-366. [6] Zhou, Wu, et. al., A novel technique for prealignment in multimodality medical image registration, BioMed research international 2014 (2014) [7] Zhang, Zhao, et al., Multi-modality medical image registration using support vector machines, 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference. IEEE, 2006: 6293-6296. [8] Ma, Jiayi, et al., Feature-guided Gaussian mixture model for image matching, Pattern Recognition 92 (2019): 231-245. [9] Ma, Jiayi, et al., Locality preserving matching, International Journal of Computer Vision 127.5 (2019): 512-531. [10] Vranic, Dejan V., Dietmar Saupe, and Jrg Richter. Tools for 3D-object retrieval: Karhunen-Loeve transform and spherical harmonics. 2001 IEEE Fourth Workshop on Multimedia Signal Processing (Cat. No. 01TH8564). IEEE, 2001: 293-298. [11] Kazhdan, Michael, Thomas Funkhouser, and Szymon Rusinkiewicz. Rotation invariant spherical harmonic representation of 3 d shape descriptors. Symposium on geometry processing. Vol. 6. 2003: 156-164. [12] Gu, Xianfeng, et. al., Genus zero surface conformal mapping and its application to brain surface mapping, IEEE Transactions on Medical Imaging 23.8 (2004): 949-958. 30
[13] Wang, Yalin, Ming-Chang Chiang, and Paul M. Thompson, Automated surface matching using mutual information applied to Riemann surface structures, International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, Berlin, Heidelberg, 2005: 666-674. [14] Ma, Jiayi, et al., Nonrigid point set registration with robust transformation learning under manifold regularization, IEEE transactions on neural networks and learning systems (2018). [15] Ma, Jiayi, et al., Robust estimation of nonrigid transformation for point set registration, Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. 2013: 2147-2154. [16] Hsieh, Cheng-Tiao., An efficient development of 3D surface registration by Point Cloud Library (PCL). 2012 International Symposium on Intelligent Signal Processing and Communications Systems. IEEE, 2012. [17] Lyu, Ilwoo, et al. Hierarchical spherical deformation for cortical surface registration. Medical image analysis 57 (2019): 72-88. [18] Liu, Heng, Jingqi Yan, and David Zhang. A neural network strategy for 3d surface registration. International Conference on Computational Science and Its Applications. Springer, Berlin, Heidelberg, 2006: 528536. [19] LEE, YINTAT, Ka Chun Lam, and Lok Ming Lui. Large deformation registration via n-dimensional quasi-conformal maps. arXiv preprint arXiv:1402.6908 (2014). 31
[20] Chan, Hei Long, and Lok Ming Lui. Detection of n-dimensional shape deformities using n-dimensional quasi-conformal maps. J. Geom. Imaging Comput 1.4 (2014): 395-415. [21] Wong, Tsz Wai, et. al., Intrinsic feature extraction on hippocampal surfaces and its applications, SIAM Journal on Imaging Sciences 5.2 (2012): 746-768. [22] Lui, Lok Ming, et. al., Landmark constrained genus zero surface conformal mapping and its application to brain mapping research, Applied Numerical Mathematics 57.5-7 (2007): 847-858. [23] Lui, Lok Ming, et. al., Optimized conformal surface registration with shape-based landmark matching, SIAM Journal on Imaging Sciences 3.1 (2010): 52-78. [24] Lui, Lok Ming, et. al., Shape-based diffeomorphic registration on hippocampal surfaces using beltrami holomorphic flow, International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, Berlin, Heidelberg, 2010: 323-330. [25] Choi, Pui Tung, Ka Chun Lam, and Lok Ming Lui, FLASH: Fast landmark aligned spherical harmonic parameterization for genus-0 closed brain surfaces, SIAM Journal on Imaging Sciences 8.1 (2015): 67-94.
32
Hei Long Chan received the M.Phil. and B.Sc. degree from the Department of Mathematics, the Chinese University of Hong Kong, Hong Kong in 2017 and 2014 respectively. He is currently a PhD student in the Chinese University of Hong Kong, Department of Mathematics. His research interests include computational quasi-conformal geometry, medical imaging, shape analysis, image processing and machine learning. Tsz Chun Yam received the M.Phil. and B.Sc. degree from the Department of Mathematics, the Chinese University of Hong Kong, Hong Kong in 2016 and 2014 respectively. His research interests include computational quasiconformal geometry, medical imaging and machine learning.
Lok Ming Lui received the Ph.D. degree in applied mathematics from the University of California at Los Angles, Los Angeles, CA, USA, in 2008. He is an Assistant Professor with the Department of Mathematics, Chinese University of Hong Kong (CUHK), Hong Kong. Before joining CUHK, he was a Post-Doctoral Scholar for 2 years with the Department of Mathematics, Harvard University, Cambridge, MA, USA. His current research interests include computational conformal and quasiconformal geometry, Teichmuller theory, surface registration, medical imaging, and shape analysis. 1
Preprint submitted to Pattern Recognition
April 11, 2019
Automatic Characteristic-Calibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs Hei Long CHANa , Tsz Chun YAMa , Lok Ming LUIb,∗ a
222B, Lady Shaw Building, The Chinese University of Hong Kong, Shatin, Hong Kong 207, Lady Shaw Building, The Chinese University of Hong Kong, Shatin, Hong Kong
b
This work is supported by HKRGC GRF (Project ID: 14303414). Data used in preparation of this article were obtained from the Alzheimers Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found here: ADNI Acknowledgement List
∗
Corresponding author Email addresses: [email protected] (Hei Long CHAN), [email protected] (Tsz Chun YAM), [email protected] (Lok Ming LUI) Preprint submitted to Pattern Recognition
November 4, 2019
Journal Pre-proof
Automatic Characteristic-Calibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs Hei Long CHAN, Tsz Chun YAM, Lok Ming LUI PII: DOI: Reference:
S0031-3203(19)30443-1 https://doi.org/10.1016/j.patcog.2019.107142 PR 107142
To appear in:
Pattern Recognition
Received date: Revised date: Accepted date:
17 May 2019 5 November 2019 27 November 2019
Please cite this article as: Hei Long CHAN, Tsz Chun YAM, Lok Ming LUI, Automatic CharacteristicCalibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs, Pattern Recognition (2019), doi: https://doi.org/10.1016/j.patcog.2019.107142
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
HIGHLIGHT • Automatic intrinsic feature extraction on hippocampal mesh surface • Feature calibration by Eigen-graphs to improve consistency and accuracy of the feature correspondences on different subjects • Fast and accurate surface registration is available in just 5 seconds with minimal RMSE • Quantitative comparison with other existing methods demonstrates the effectiveness and efficiency of the proposed model on real data from the ADNI database
1
Automatic Characteristic-Calibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs Hei Long CHANa , Tsz Chun YAMa , Lok Ming LUIb,∗ a
222B, Lady Shaw Building, The Chinese University of Hong Kong, Shatin, Hong Kong 207, Lady Shaw Building, The Chinese University of Hong Kong, Shatin, Hong Kong
b
Abstract In this paper, we propose an efficient algorithm, the ACC-REG, to automatically extract intrinsic key characteristics on hippocampal mesh surfaces and hence compute an accurate registration mapping between them. Given a pair of hippocampal surface mesh, the proposed algorithm constructs the eigengraphs, an intrinsic feature on the surface, on each surface as its representative. The eigen-graphs are then calibrated along the longitudinal direction of the hippocampal surfaces. Accurately corresponded intrinsic characteristics on each hippocampus can thus be extracted. As a result, the two surfaces can be registered with improved accuracy and low computation cost. Experiments on ADNI data demonstrate the effectiveness of the proposed ACC-REG model over existing methods. Keywords: Feature Extraction, Feature Correspondence, Surface Registration, Surface deformation, Eigen-graph
∗
Corresponding author Email addresses: [email protected] (Hei Long CHAN), [email protected] (Tsz Chun YAM), [email protected] (Lok Ming LUI)
Preprint submitted to Pattern Recognition
January 4, 2020
1. Introduction
Surface based shape analysis has a wide range of applications. For instance, it is a key tool in medical science for disease classification and further diagnosis. With a reliable registration process, MR/CT images of different patients scanned at various time can be accurately corresponded and then compared, providing a promising environment for further analysis. One example is the registration between hippocampal surfaces which is very important to tracing the progression of the Alzheimer’s disease. While there are existing models particularly designed for hippocampal surface registration, accuracy is still hindered by different difficulties. On one hand, there are not much obvious and natural geometric features on the hippocampal surfaces that can be utilized to define landmark correspondences. On the other hand, scanning error due to hardware is still common and puts great challenges in extracting well-defined features on the hippocampal surfaces. In addition, large deformations of the hippocampus are expected in case of disease, which further ruins the stability and accuracy in the registration process. To address the above problems, a robust and automatic landmark based registration algorithm, ACC-REG, is proposed in this paper for accurate hippocampal surface registration. The intrinsic Laplace-Beltrami operator is utilized to generate feature correspondences on the hippocampus with high stability. Given two hippocampal surface meshes, the Helmholtz differential equation involving the Laplace-Beltrami operator is solved and hence an eigen-graph, an intrinsic feature on the hippocampal surface, is constructed on each of them. Inspired by the histogram matching technique from the 2
study of statistics, the eigen-graphs are then calibrated along the longitudinal direction on the hippocampus so that the accuracy of the correspondences between the two hippocampus are greatly improved. Then, well-corresponded intrinsic characteristic curves on the hippocampus are extracted. Finally, a fast and accurate surface registration algorithm is applied to compute the correspondence over the whole surface. The contribution of this paper mainly lies in three aspects. Firstly, the proposed characteristic curves are closely related to the unique geometry of the hippocampus. So the landmark correspondences are meaningful and robust. Besides, the extraction of the characteristic curves involves only a direct and linear process and hence efficiency is guaranteed. Finally, calibrating the eigen-graphs greatly reduces the influences of global and local deformations on the hippocampal surfaces, which in turn greatly improves the accuracy of the registration between the hippocampal surfaces, as being verified by various experiments. 2. Related Works
Registration between anatomical structures has a wide range of real applications, for example, disease classification [1, 2], morphometry analysis [3, 4], and surface reconstruction [5], etc. The research field is mainly divided into two main streams, namely, 2D and 3D registrations. For 2D registration, scanned images of tissues are one-to-one corresponded according to the depth of the scanning. In [6], Zhou et. al. proposed a novel technique for prealignment in both mono-modality and multi-modality image 3
registration based on statistical correlation of gradient information. In [7], Zhang et. al. proposed a multi-modality medical image registration model using support vector machines and incorporating with prior joint intensity distribution between the two imaging modality based on registered training images. In [8], Ma et. al. proposed a robust feature-guided Gaussian mixture model for image matching which can handle both rigid and non-rigid image transformation robustly. In [9], Ma et. al. proposed a feature mismatching removal scheme which removes mismatches from given putative image feature correspondences by maintaining the local neighborhood structures of those potential true matches. The proposed scheme can also be applied to medical imaging to seek for reliable correspondences between images of the same tissue. However, we note a general problem for utilizing 2D registration in medical analysis. Namely, the slicing of the tissue into 2D images is solely based on the scanning hardware and the gesture of the subject. Regarding to the fact that the position and the shape of the target tissue may have great variation among subjects, images of the same tissue at the same depth across subjects may not necessarily have a robust and meaningful point-wise correspondence. Comparing to 2D registration, 3D registration accounts more for geometric information to improve the robustness and accuracy of the registration. One commonly used registration model for spherical topology surface is the spherical harmonics (SPHARM) algorithm. For example, in [10], DV Vranic et. al. built a database of surface query and designed a SPHARM-based algorithm to retrieve similar objects using the surface database. In [11], Michael Kazhdan et. al. used the spherical harmonics to represent 3D shapes to solve
4
the rotation invariant shape matching problem. Later, other models are designed for surface registration. In [12], Gu et. al. proposed an algorithm for global conformal parametrizations based on the structure of the cohomology group of the holomorphic one-forms for surfaces. In particular, for genus zero surfaces, the proposed algorithm finds an unique mapping between any two manifolds by minimizing the harmonic energy of the mapping, which can then be applied to the cortical surface matching problem. In [13], Wang et. al. developed an algorithm which extends the mutual information method to automatically match general 3D surfaces to automate the matching of surface features. However, both [12, 13] are built in the absence of manual landmark correspondences on the surfaces, which does not suit our problem setting. In [14], Ma et. al. solves the problem of nonrigid point set registration by designing a robust transformation learning scheme. The principle is to iteratively establish point correspondences and learn the nonrigid transformation between two given sets of points. In [15], Ma et. al. presented a new point matching algorithm for robust nonrigid registration which iteratively recovers the point correspondence and estimates the transformation between two point sets. Their proposed model is general and has much applications, including medical analysis. In [16], Cheng-Tiao Hsieh presented a model to utilize the open source Point Cloud library (PCL) to develop a series of computational registration processes efficiently and robustly. The model is a key process of the Digital Face-Inspection (DFI) system and has lots of medical applications. In [17], Ilwoo Lyu et. al. designed a cortical surface registration model using then hierarchical spherical deformation which achieves group-wise surface correspondence with no template selection
5
bias. Recently, surface registration has also been incorporated with neural networks, for example, in [18]. It is also worth mentioning that high-dimensional registration is also under study. In [19], Lee et. al. proposed a new method to obtain registration between n-dimensional manifolds with very large deformations. Highdimensional registration also has medical applications to register between volumetric data, for example, in [20]. Finally, we would like to review some literature which lead to our proposed model. In [21], Wong et. al. proposed a novel approach for extracting two intrinsic feature curves on hippocampal surfaces. The model is used to study the geometric changes of hippocampus in case of the Alzheimer’s disease (AD). In [22, 23], Lui et. al. proposed a model to generate optimized conformal parametrizations of brain surfaces which match prior landmark curves exactly with shape-based correspondences between them. The idea is further extended to general surfaces in [23, 24]. These works provide a hint to our proposed model on defining and utilizing meaningful landmark correspondences on the hippocampal surfaces. In [25], Choi et. al. proposed a fast and efficient landmark-based surface registration model for genus zero surfaces. The proposed model computes the optimized spherical harmonic parametrization with consistent landmark alignment. It is not just accurate but also speedy as their model has a totally linear implementation. 3. The proposed model
The details of the proposed ACC-REG model is explained in this sec6
tion. Registration is a process to extract features on each subject, find the correspondence over the feature sets, and finally estimate and construct the surface correspondence based on the feature sets correspondence. In view of this, intrinsic features provide much stability and reliability over extrinsic features. Therefore, in this work, the eigen-graph is utilized as the (intrinsic) features for landmark correspondences. Given a hippocampal surface H, consider the Helmholtz equation on H: − ∆f (x) = λf (x),
(1)
where ∆ denotes the Laplace-Beltrami operator on H. The common geometry of hippocampus ensures the eigenvalues of H always form a non-negative increasing sequence 0 ≤ λ1 ≤ λ2 ≤ · · · → +∞. Therefore, one can define an eigen-graph on H as follows: Definition 1 (The Eigen Graph). The eigen graph of a hippocampal surface H is defined to be the eigenvector FH of the Helmholtz equation (1) such that the associated eigenvalue λ satisfies: λ > 0,
ˆ for all other eigenvalues λ ˆ λ≤λ
(2)
Figure (1) shows an example of the eigen-graph FH on H. The eigengraph FH stores much topological and geometric information of H. For instance, FH always has exactly 2 critical points on the two tips of H. Besides, the values of FH increases from one critical point to another critical point. This gives a parametrization of H as follows:
7
Figure 1: Example of the eigen-graph FH on a hippocampal surface H, function values is normalized from blue (0) to red (1)
Theorem 1. Let H be a hippocampal surface and FH : H → [0, 1] be the eigen-graph on H. Let p1 , p2 ∈ H be the two critical points of FH such that
FH (p1 ) < FH (p2 ). For any 0 < < 12 (FH (p2 ) − FH (p1 )), FH−1 ([FH (p1 ) + , FH (p2 ) − ])
(3)
is diffeomorphic to I × S 1 , where I = [FH (p1 ) + , FH (p2 ) − ] and S 1 is a circle. From the theorem, there exists a smooth mapping from I × S 1 to FH−1 (I), where I ⊂ FH (H \ {p1 , p2 }).
Hence, the hippocampal surfaces can be
parametrized by longitudinal and azimuthal coordinates. In view of this, Wong et. al. [21] proposed to extract two intrinsic feature curves γ1 , γ2 as landmark correspondences by minimizing the following energy: 1 {γ1 , γ2 } = argmax b−a
Z
a
b
|γ1 (t) − γ2 (t)|dt − c2
Z
a
b
|γ˙1 (t)|dt − c1
Z
a
b
|γ˙2 (t)|dt, (4)
where c1 , c2 ≥ 0 are the weighting parameters and [a, b] = f (H) is the range of the un-normalized first nontrivial eigenfunction on H. The motivation of the definition (4) is that the feature curves should be separated as much 8
(a) HP from normal people
(b) HP from AD patient
Figure 2: Examples of hippocampal surfaces: (left) From normal subject; (right) From patient
as possible while each of them should also be as smooth as possible. In Wong et. al.’s work, after extracting γ1 and γ2 , the hippocampal surfaces are parametrized onto the domain K = {(α, θ) ∈ [0, 1] × [0, 2π)} by minimizing the energy: Φ=
arg min Φ(Γ1 )={θ=0},Φ(Γ2 )={θ=π}
Z
H
||∇Φ||2 .
(5)
This is called the eigen-harmonic parametrization (EHP) of H. However, there are two major drawbacks of Wong et. al.’s model. Firstly, the implementation is inefficient as it involves gradient descents of two energy functional, with each of them usually takes over thousands of iterations to converge numerically. Secondly, it is realized that the progression of the eigen-graph FH on H from one tip to another varies significantly among different subjects. In fact, in real application, it is found that the geometry of the hippocampus from different subject can have large variations due to illness and defect (see figure (2)). In such situation, feature correspondences may be inaccurate and cause error in the registration. Therefore, we propose the following calibration of the eigen-graphs to improve the accuracy of the eigen-graph correspondences.
9
3.1. Calibration of the Eigen-Graph
Suppose FH1 and FH2 are eigen-graphs corresponding to two hippocampal surfaces H1 and H2 respectively. Using the longitudinal parametrization of H1 and H2 as described in Theorem (1) and normalizing the parametrization domain as [0, 1], one has: Theorem 2. For all t ∈ [0, 1], there exists an unique z ∈ R+ such that FHj (Hj (t)) = z on the slide Hj (t). Therefore, FHj can be re-defined as functions on [0, 1] for all j = 1, 2. Then, the cumulative distribution function FH∗ j of FHj can be defined by: FH∗ j (x)
=
Z
x
FHj (t)dt,
(6)
0
for all x ∈ [0, 1] and j = 1, 2.
Now, for each x1 ∈ [0, 1], there exists x2 ∈ [0, 1] such that FH∗ 1 (x1 ) =
FH∗ 2 (x2 ). By the monotonicity of FH∗ 1 and FH∗ 2 , there exists:
Definition 2 (Histogram Matching Function). Suppose that FH∗ 1 , FH∗ 2 : [0, 1] → R are two monotone increasing functions. The histogram matching
function G12 : [0, 1] → [0, 1] for FH∗ 1 , FH∗ 2 is defined by G12 (x1 ) = x2 ,
(7)
for all x2 ∈ [0, 1] such that FH∗ 1 (x1 ) = FH∗ 2 (x2 ), and for all x1 ∈ [0, 1]. Using the histogram matching function G12 , the eigen-graphs FH1 and FH2 can be calibrated as follows: 10
Definition 3 (Calibrated Eigen-Graph). Let H1 , H2 be two hippocampal surfaces, FH1 , FH2 be their eigen-graphs, FH∗ 1 , FH∗ 2 be their cumulative eigengraphs, and G12 : [0, 1] → [0, 1] be the histogram matching function matching FH∗ 1 to FH∗ 2 . The calibrated eigen-graph of H1 is defined to be FHH12 = FH1 ◦ G12 .
(a) FH1 on H1
(b) FH2 on H2
(8)
(c) FHH12 on H1
Figure 3: Example of calibrating the eigen-graphs on two hippocampal surfaces H1 and H2 .
The calibrated FHH12 can be viewed as a synchronized version of FH1 with respect to FH2 . Figure (3) demonstrates an example of the calibration process and figure (4) shows the histograms of FH1 and FH2 and the comparison between them. It can be seen that after the calibration, FHH12 has almost the same distribution as FH2 . Hence, any features correspondences generated using the calibrated eigen-graphs are more coherence and accurate. In particular, intrinsic characteristic curves can be extracted to provide wellmatched features for accurate registration. 3.2. Extraction of Intrinsic Characteristic Curves from Eigen-Graphs The extraction of feature characteristic curves mainly involves two parts. In the first part, two feature curves mutually opposite to each other are 11
(a) Histogram of FH1
(b) Histogram of FH2
(c) Comparing FH1 and FH2
(d) Histogram of FHH12
Figure 4: Demonstration and comparison of the histograms of FH1 and FH2 . Each graph is a histogram plot, with the y-axis being the number of vertices having the particular value in the eigen-graph indicated by the x-axis.
extracted as the main features on the hippocampus. Secondly, extra assisting feature curves are extracted based on the main feature curves to help improving the registration accuracy. To begin with, we extract two fundamental intrinsic characteristic curves: Definition 4 (Fundamental Intrinsic Characteristic Curves). Let H be a hippocampal surface and FH be its eigen-graph. The fundamental intrinsic characteristic curves on H are defined to be the curves γ1 , γ2 : [0, 1] → H such that {γ1 , γ2 } = argmax
Z
1
0
12
|γ1 (t) − γ2 (t)|2 dt.
(9)
By the geometry of hippocampal surfaces, one curve (γ1 ) must lie on the most curvy region along the longitudinal direction with the other one (γ2 ) on its opposite side. However, solving (9) directly requires high computation cost and may result in zigzag curves depending on the mesh triangulation. Therefore, the following notion of discrete curve on mesh is employed: Definition 5 (Discrete Curves on Surfaces). Let H be a hippocampal surface mesh or point cloud, that is, H = {Vi }N i=1 for some N ∈ N. A discrete curve γ on H is defined by γ = {vi }ni=1 ,
(10)
where {vi }ni=1 ⊂ {Vi }N i=1 and the connectivity of γ is given by v1 → v2 → v3 → · · · → vn .
(11)
Under this definition, after normalizing FH to be ranged [0, 1], γ1 , γ2 are extracted by γ1 = {vi }ni=1 , γ2 = {wi }ni=1 , where {vi , wi } = arg max |v − w|,
(12)
Ωi ()
where Ωi () = FH−1 ([ ni − , ni + ]) is the subset of H at the query level-set
1 FH−1 ( ni ) with bandwidth ∈ (0, 2n ). Figure (5)(a),(b) shows an example of
the two (discrete) intrinsic characteristic curves γ1 and γ2 on a hippocampal surface. 13
(a) γ1
(b) γ2
(c) γ1f it
(d) γ2f it
(e) γ1proj
(f) γ2proj
Figure 5: Examples of characteristic curves extracted on a hippocampal surface H. Vertices on discrete curves are connected sequentially for visualization.
Note that γ1 , γ2 are not smooth enough in general. To remedy the situation, the following polynomial fitting calibration scheme is used. Firstly, the polynomial fitting algorithm is applied on γj to get γjf it : [0, 1] → R3 , where γjf it (t) = (xj (t), yj (t), zj (t)) ∀t ∈ [0, 1], Z 1 f it γj = argmin |γj (t) − γjf it (t)|2 dt,
(13) (14)
0
for all j = 1, 2. From the examples shown in figure (14), γjf it are far smoother than γj .
14
However, there is a drawback that γjf it sits in R3 indeed and may not be lie on H. Therefore, a projection of γjf it onto H is applied to get γjproj , a smooth intrinsic characteristic curves lying on H, which is defined by: Definition 6 (Projected Smooth Intrinsic Characteristic Curves). Let H be a hippocampal surface, FH be its eigen-graph, and γ1f it , γ2f it be its two polynomial fitted intrinsic characteristic curves of degree k. The projected smooth intrinsic characteristic curves of degree k are defined to be the curves γjproj : [0, 1] → H such that γjproj
= argmin
Z
0
1
|γjproj (t) − γjf it (t)|2 dt,
(15)
for all j = 1, 2. Numerically, γ1f it and γ2f it partitions the curve domain [0, 1] into 0 < l1 < l2 < · · · < lm < 1,
(16)
where 1 lk = (FH (vik ) + FH (wik )), 2
k = 1, 2, . . . , m.
(17)
In order to capture the projected curve γjproj on H, we re-partition 1 1 2 n−1 [0, 1] = [0, ] ∪ [ , ] ∪ · · · ∪ [ , 1] n n n n
(18)
into equally spaced sub-intervals. Then, the projected curve γjproj is extracted by 0
0
{vi , wi } = arg min(||v − γ1f it ||0 + ||w − γ2f it ||0 ) v,w∈Ωi ()
15
(19)
Since Ωi () is a finite set, γ1proj and γ2proj can be extracted directly by solving (19) using the exhaustive search. Figure (5)(e),(f) demonstrates a pair of projected curves on H. Comparing figure (5)(e),(f) with figure (5)(c),(d), the projected curves are indeed less smooth but lie within H. Also, comparing figure (5)(e),(f) and figure (5)(a),(b), the projected curves are still smoother than the fundamental curves. The two projected curves γ1proj and γ2proj are selected as the (intrinsic) characteristic curves for accurate surface registration. In order to enhance the robustness of the registration, some more feature curves lying on the centre part of H are also extracted as characteristic curves, apart from γ1proj , γ2proj . Given the hippocampal surface H and the two projected smooth intrinsic characteristic curves γ1proj and γ2proj , an azimuthal coordinate function θ : H → [0, 2π) can be defined such that θ=0
on γ1proj ,
θ=π
γ2proj .
on
(20)
In view of this, more characteristic curves can be defined on H, as inspired by [21], by: Definition 7 (Assisting Characteristic Curves). Let H be a hippocampal surface, FH be its eigen-graph, γ1proj , γ2proj be its two projected smooth intrinsic characteristic curves of degree k, and θ be an azimuthal coordinate function on H. The assisting characteristic curves are defined to be the curves Γk : [0, 1] → H where: Γk = θ−1 (
k π), M
k = 1, 2, . . . , 2M − 1. 16
(21)
In discrete case, for each region Ωi (), by applying the principal component analysis, Ωi () can be projected onto the plane P (Ωi ()) spanned by its first two component vector. Note that Ωi () intersects with γ1proj , γ2proj at 0
0
vi , wi respectively. Let
1 0 0 di = (vi + wi ) 2
(22)
be the mid-point of the two vertices and ui be the unit vector at di , perpendicular to the P (Ωi ()), pointing to the direction for which FH is increasing. 0
0
For each vertex z ∈ Ωi (), let zproj , (vi )proj , (wi )proj , dproj and uproj be the i i 0
0
projection of z, vi , wi , di and ui onto the plane P (Ωi ()) respectively. Then, by computing −1 <
σ = cos
0
(zproj − dproj ), ((wi )proj − dproj )> i i , 0 proj proj |zproj − di | · |(wi )proj − di |
(23)
the azimuthal function θ : H → [0, 2π) can be assigned to z depending on 0
, zproj − dproj , uproj } by the orientation of Sk = {(wi )proj − dproj i i i
θ(z) =
σ
if Sk is positively-oriented,
2π − σ
(24)
otherwise.
Therefore, each Γk in (21) can be extracted by Γk = {zik }ni=1 where θ(zik )
k = arg min θ(z) − π M z∈Ωi ()
(25)
for each i = 1, 2, . . . , n. Note that by definition, Γ0 = γ1proj and ΓM = γ2proj .
17
Figure 6: Example of all the characteristic curves extracted on a hippocampal surface (black curves indicate γ1proj , γ2proj , green curves indicate Γk ’s)
Therefore, the set of all characteristic curves is: C = {Γ0 = γ1proj , Γ1 , Γ2 , . . . , ΓM = γ2proj , . . . , Γ2M −1 , }.
(26)
Figure (6) demonstrates all the characteristic curves on a hippocampal surface. 3.3. Surface Registration using Calibrated Characteristic Curves
Once the characteristic curves are computed, they can be utilized as landmark correspondences for the surface registration. To promote efficiency and accuracy, the FLASH algorithm [25] is employed in this step. The FLASH algorithm is a landmark-based registration algorithm designed particularly to register between surfaces having spherical topology. The FLASH algorithm mainly consists of two steps. In the first step, a spherical conformal parametrization function Ψi is computed for each Hi . For each of H1 , H2 , a parametrization ψj : Hi \ {Tj } → C is computed by
18
𝐻1
𝑃𝑁−1
𝜓1
𝑃𝑆
𝑃𝑆−1
𝑔
𝑃𝑁 1) 𝑅(𝐻
𝑃𝑆−1
𝑔𝐷
𝑃𝑆
𝑃𝑁−1
𝑓∘𝜑
(Landmark)
𝐻2
𝑃𝑁−1
𝜓2
𝑃𝑆
𝑔
(Landmark)
𝑃𝑆−1
Figure 7: Illustration of the pipeline of the FLASH algorithm
solving the linear system: P
[u,v]∈Hj
kuv (ψj (u) − ψj (v)) = 0
ψi (tj ) = bl
if u 6= tj1 , tj2 , tj3 ,
(27)
if l = 1, 2, 3;
l
where kuv is the coefficient of the discrete Laplace-Beltrami operator, Tj = [tj1 , tjj , tj3 ] is a triangle of the mesh Hj and bl is a preassigned coefficient. In order to minimize the angle distortion of the parametrization, ψj is further composed with stereographic projections PN , PS , the north pole projection and the south pole projection respectively, and also a conformality distortion minimization function g which can also be computed by solving a linear 19
system. Finally, the spherical conformal parametrization is given by Ψi = PS−1 ◦ g ◦ Ps ◦ PN−1 ◦ ψj
(28)
for each j = 1, 2. After obtaining the parametrization Ψj for each hippocampal surface Hj , let Bj = Ψj (Hj ) be the ball in R3 . By directly solving for the Mobius P n 0 0 2 transformation f (v) = av+b on B1 to minimize 2M i=1 |f (vi )−wi | , landmark f1 = f (B1 ), then B f1 , B2 are projected onto mismatching is eliminated. Write B
C by PN and then the landmark constrained harmonic map ϕ : C → C can be obtained by solving the linear system P 0 kuv (ϕ(u) − ϕ(v)) + λ(ϕ(u) − wi ) = 0 [u,v]∈H 1 P [u,v]∈H1 kuv (ϕ(u) − ϕ(v)) = 0 ϕ(tjl ) = tjl
0
if u = vi , i = 1, 2, . . . , 2M n, 0
0
if u ∈ / {v1 , . . . , vn , tj1 , tj2 , tj3 }, if l = 1, 2, 3, (29)
where λ > 0 is a weighting parameter. Afterwards, a few more stereographic projections and a quasi-conformal mapping is applied to minimize the angle ˆ : H1 → B2 is computed by distortion. Then, the pre-registration mapping R ˆ = P −1 ◦ gD ◦ PS ◦ P −1 ◦ f ◦ ϕ ◦ PN ◦ Ψ1 . R S N
(30)
Finally, the registration from H1 to H2 can be obtained by conventional interpolation methods, e.g. scatter interpolation function I : R3 → H2 . ˆ is the desired registration. Then, R := I ◦ R
Figure (7) demonstrates the pipeline of the FLASH algorithm. For the 20
detailed implementation, readers are referred to [25]. There are several advantage on using FLASH algorithm on registration between hippocampal surfaces. Firstly, the FLASH algorithm involves only solving sparse linear systems and hence efficiency of the algorithm is guaranteed, whereas in [21] it may take over thousands of iterations for the gradient descent for the harmonic energy minimization to converge numerically. Secondly, using the FLASH algorithm does not only minimizes the angle distortion as in [21] but also enforces exact landmark correspondences. Therefore, the power of the calibrated characteristic curves as landmark correspondences can be fully utilized. 3.4. Implementation of the ACC-REG algorithm
The pipeline of the proposed ACC-REG algorithm is concluded in the Algorithm 1 table for reference. 4. Experiments
In this section, experimental results are recorded to evaluate the proposed ACC-REG algorithm, in terms of improving the registration accuracy and also the computation cost. For this purpose, 45 hippocampal surfaces {Hj }45 j=1 from the ADNI database are extracted. Each surface mesh consists of approximately 9000 vertices, and is re-orientated to be centered at origin and resized so that max(||vj ||2 ) = 1. We perform several experiments. For every experiments, each of the 45 surfaces are subjected to a particular manual deformation which simulates abnormal deformations in case of disease. 21
Algorithm 1: ACC-REG Input: Hippocampal surfaces H1 , H2 , Parameters n, M, Output: Surface Registration Functions R : H1 → H2 1 Compute the Eigen-Graph FH1 , FH2 by solving (1) on H1 , H2 respectively; 2 Compute the histogram matching function G12 by (6,7); H1 3 Compose the Calibrated Eigen-Graph FH = FH1 ◦ G12 ; 2 4 for H1 , H2 do 5 Compute the Fundamental Intrinsic Characteristic Curves γ1 = {vi }ni=1 , γ2 = {wi }ni=1 by solving (12); 6 Apply polynomial fitting to obtain γ1f it , γ2f it on R3 by (14); 7 Compute the Projected Smooth Intrinsic Characteristic Curves 0 0 γ1proj = {vi }ni=1 , γ2proj = {wi }ni=1 by (19); 8 Compute the Assisting Characteristic Curves Γ1 , Γ2 , . . . , Γ2M −1 by (25); ˆ : H1 → R3 by the FLASH 9 Compute the pre-registration Function R algorithm using the two sets of {Γk }’s on H1 , H2 as landmark correspondences; ˆ 10 Compose the registration mapping R = I ◦ R. Different registration algorithms are applied to register the original surfaces to the deformed surfaces, including the spherical harmonics (SPHARM) algorithm, the MP-RPM algorithm [14], the EHP-based algorithm [21] and the proposed ACC-REG algorithm. Quantitative comparisons are done among the algorithms in terms of RMSE and also computation time. The eigengraphs generated from the proposed ACC-REG algorithm and from the original EHP-based algorithm [21] are also compared. Throughout all the experiments, the parameter setting for the proposed ACC-REG algorithm is: n = 10, M = 3, = 0.03. All the parameter setting for the FLASH algorithm follows that in [25]. All the parameter setting for
22
other algorithms follow that listed in the original publications. In particular, we use all the spherical harmonics up to the 45-th degree for the SPHARM algorithm. The registration results are visualized in the figures in this section. In each plot, the colors on the surface indicate the value of the eigen-graph, from blue (0) to red (1). The black rings on the surface is the level-set of the eigen-graph for reference. Experiment 1: In the first experiment, abnormal elongation of the hippocampus is simulated. For each hippocampus Hj , let Ω ⊂ Hj denotes the area for which FH (Ω) ≥ 0.7. By applying the PCA, the first component d is taken as the deformation direction. Each vertex vj ∈ Ω is then deformed
by vj∗ = vj + (FH (vj ) − 0.7)d. Figure (8)(a),(b) demonstrates an example of
such deformation. Figure (8)(c),(d),(e) shows the registration results of the MP-PRM algorithm [25], the EHP-based algorithm [21] and the proposed ACC-REG algorithm respectively. it can be seen that the distribution of the eigen-graph on the resultant surface is the most natural using the proposed algorithm. The statistics of the RMSE using different algorithms are recorded in table (1). Note that the mean RMSE between the original surfaces and the deformed surfaces is recorded in the row “Control” for reference to the scale of the deformation. Experiment 2: In the second experiment, a point diffusion deformation is applied to Ω = FH (Ω) ≥ 0.7 to simulate inflammation of the hippocampus. The simulation is as follows. The mean vertex v¯ = meanvj is computed over Ω. Then, the unit vector nj from v¯ to vj is computed for each vj ∈ Ω.
Finally, the deformed surface is defined by vj∗ = vj + 0.5nj . Same as the first
23
(a) Original surface
(b) Deformed surface
(c) Result of SPHARM algorithm
(d) Result of MP-PRM (e) Result of EHP-based (f) Result of proposed algorithm ACC-REG algorithm algorithm Figure 8: Sample registration results of the first experiment
experiment, different registration algorithms are used to register between the original surfaces and the deformed surfaces. One example of such deformation and the registration results by different algorithms is plotted in figure (9). The numerical results, i.e. the RMSE, for each algorithm is recorded in the table (1). Experiment 3: In the third experiment, the deformation is designed to simulate tumors on the hippocampal surfaces. For each surface Hj , on each region Ωi () as defined in (12), a random vertex is picked and a tumor is built centered at this vertex with height 0.1 along the outward normal direction. An example of such deformation and the visualization of the registration results due to different algorithms, and the statistical result of the RMSE of different algorithms are recorded in the figure (10), and the table (1) respectively. The above three experiments totally consists of 45 × 3 = 135 subjects 24
(a) Original surface
(b) Deformed surface
(c) Result of SPHARM algorithm
(d) Result of MP-PRM (e) Result of EHP-based (f) Result of proposed algorithm algorithm ACC-REG algorithm Figure 9: Sample registration results of the second experiment
which should be enough to draw a rigid conclusion. Referring to the table (1), it is evident that the proposed ACC-REG algorithm out-performs all other registration algorithms in terms of RMSE. This can be explained by the effectiveness of the calibrated eigen-graphs and hence the intrinsic feature correspondences extracted from the eigen-graphs. The high accuracy in the feature correspondences favors a lot to the registration process. Moreover, the parameters M, n in the ACC-REG algorithm can be altered to adjust the number of the feature curves used and the number of query vertices involved in each feature curve. This helps the ACC-REG algorithm to adapt well to different database with consistent accuracy. In view of the table (2), the proposed ACC-REG algorithm also has a merit of fast computation. To further investigate the computation efficiency of the proposed algorithm, it is recorded that the mean computation time for computing the eigen-graph is 0.27 seconds. The mean computation time 25
(a) Original surface
(b) Deformed surface
(c) Result of SPHARM algorithm
(d) Result of MP-PRM (e) Result of EHP-based (f) Result of proposed ACC-REG algorithm algorithm algorithm Figure 10: Sample registration results of the third experiment
for the landmark correspondences (i.e. Γk ’s) is 0.80 seconds and that for the FLASH algorithm is 3.12 seconds. In particular, the reason for the fast computation of the landmark correspondences is that, the dimension of the linear systems involved is greatly reduced by partitioning the surface H into Ωi ()’s and find the landmark vertices on each Ωi () one-by-one. Also, the FLASH algorithm itself is also known for the fast computation speed. Combining these together, the proposed ACC-REG algorithm has an advantage on the computation speed over other algorithms. Comparatively, the accuracy of the SPHARM model is limited by the maximum degree of spherical harmonics used. Higher degree of spherical harmonics store much more local information and hence gives more accurate registration result, but the computation cost significantly increases in using more spherical harmonics. As a result, the model cannot achieve a fast and accurate registration. 26
Algorithm Control SPHARM MR-PRM EHP-based ACC-REG
Experiment 1 0.1170 0.0527 0.0433 0.0474 0.0333
Experiment 2 0.0862 0.0754 0.0323 0.0524 0.0234
Experiment 3 0.1000 0.0523 0.0351 0.0482 0.0330
Table 1: Mean RMSE of the registration results in the experiments of different algorithms.
As for the MR-PRM algorithm, it is noted that its accuracy is highly comparable with our proposed ACC-REG model. However, the computation complexity of the MR-PRM algorithm is quite high. Although the authors has designed an alternative way to slightly reduce the computation complexity (see [14]), the computation time is still significantly higher than the proposed ACC-REG model. This puts a great challenge in dealing with batch-wise tasks. The EHP-based model lacks the calibration of the eigen-graph. This makes the landmark correspondences inaccurate. Besides, using arc length as the only regularization to the feature curve makes no reference to the geometry of the hippocampal surface. This further hinders the accuracy in the correspondences of the feature curves. The proposed ACC-REG model has a higher accuracy then the EHP-based model because it incorporates more geometric information in the model by calibrating the eigen-graphs and the feature curves by histogram matching and polynomial fitting respectively. Finally, the registration process of the ACC-REG model is faster than the EHP-based model mainly due to the application of the FLASH algorithm.
27
Algorithm Computation time
SPHARM 49.49s
MR-PRM 46.29s
EHP-based 8.15s
ACC-REG 4.19s
Table 2: Mean computation time of different algorithms.
5. Conclusion and Future Work
In this paper, an advanced surface registration algorithm called the ACCREG algorithm is proposed. The ACC-REG algorithm uses eigen-graphs to understand the geometry of the given hippocampal surfaces. Firstly, the eigen-graphs are mutually calibrated to synchronize their progressions along the corresponding surfaces. Secondly, on each surface, intrinsic characteristic curves are extracted and are calibrated using a projected polynomial fitting approach. Assisting feature curves are then extracted to further incorporate with more geometric information. Taking these characteristic curves as landmark correspondences, a fast and efficient registration model is utilized to register the given surfaces accurately. Through experiments, it is clear that ACC-REG helps to improve the accuracy of the registration process over other algorithms including the SPHARM algorithm, the MR-PRM algorithm and an EHP-based model. The ACC-REG algorithm also has a merit of fast computation, in which two surface meshes consisted of roughly 9, 000 vertices can be registered accurately within 5 seconds. In the future, further applications of the proposed ACC-REG algorithm including disease classification, registration of high genus surfaces, object simulation, etc., will be investigated.
28
Acknowledgement This work is supported by HKRGC GRF (Project ID: 14303414). Data used in preparation of this article were obtained from the Alzheimers Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found here: ADNI Acknowledgement List Reference References [1] Chan, Hei Long, Hangfan Li, and Lok Ming Lui, Quasi-conformal statistical shape analysis of hippocampal surfaces for Alzheimer s disease analysis, Neurocomputing 175 (2016): 177-187. [2] Gerardin, Emilie, et. al., Multidimensional classification of hippocampal shape features discriminates Alzheimer’s disease and mild cognitive impairment from normal aging, Neuroimage 47.4 (2009): 1476-1486. [3] Meintjes, E. M., et. al., A tensor-based morphometry analysis of regional differences in brain volume in relation to prenatal alcohol exposure, NeuroImage: Clinical 5 (2014): 152-160. [4] Liu, Ya-Shu, Han-Bing Yan, and Ralph R. Martin. As-rigid-as-possible surface morphing. Journal of computer science and technology 26.3 (2011): 548-557. 29
[5] Reuter, Martin, Franz-Erich Wolter, and Niklas Peinecke. LaplaceBeltrami spectra as “Shape-DNA” of surfaces and solids. Computer-Aided Design 38.4 (2006): 342-366. [6] Zhou, Wu, et. al., A novel technique for prealignment in multimodality medical image registration, BioMed research international 2014 (2014) [7] Zhang, Zhao, et al., Multi-modality medical image registration using support vector machines, 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference. IEEE, 2006: 6293-6296. [8] Ma, Jiayi, et al., Feature-guided Gaussian mixture model for image matching, Pattern Recognition 92 (2019): 231-245. [9] Ma, Jiayi, et al., Locality preserving matching, International Journal of Computer Vision 127.5 (2019): 512-531. [10] Vranic, Dejan V., Dietmar Saupe, and Jrg Richter. Tools for 3D-object retrieval: Karhunen-Loeve transform and spherical harmonics. 2001 IEEE Fourth Workshop on Multimedia Signal Processing (Cat. No. 01TH8564). IEEE, 2001: 293-298. [11] Kazhdan, Michael, Thomas Funkhouser, and Szymon Rusinkiewicz. Rotation invariant spherical harmonic representation of 3 d shape descriptors. Symposium on geometry processing. Vol. 6. 2003: 156-164. [12] Gu, Xianfeng, et. al., Genus zero surface conformal mapping and its application to brain surface mapping, IEEE Transactions on Medical Imaging 23.8 (2004): 949-958. 30
[13] Wang, Yalin, Ming-Chang Chiang, and Paul M. Thompson, Automated surface matching using mutual information applied to Riemann surface structures, International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, Berlin, Heidelberg, 2005: 666-674. [14] Ma, Jiayi, et al., Nonrigid point set registration with robust transformation learning under manifold regularization, IEEE transactions on neural networks and learning systems (2018). [15] Ma, Jiayi, et al., Robust estimation of nonrigid transformation for point set registration, Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. 2013: 2147-2154. [16] Hsieh, Cheng-Tiao., An efficient development of 3D surface registration by Point Cloud Library (PCL). 2012 International Symposium on Intelligent Signal Processing and Communications Systems. IEEE, 2012. [17] Lyu, Ilwoo, et al. Hierarchical spherical deformation for cortical surface registration. Medical image analysis 57 (2019): 72-88. [18] Liu, Heng, Jingqi Yan, and David Zhang. A neural network strategy for 3d surface registration. International Conference on Computational Science and Its Applications. Springer, Berlin, Heidelberg, 2006: 528536. [19] LEE, YINTAT, Ka Chun Lam, and Lok Ming Lui. Large deformation registration via n-dimensional quasi-conformal maps. arXiv preprint arXiv:1402.6908 (2014). 31
[20] Chan, Hei Long, and Lok Ming Lui. Detection of n-dimensional shape deformities using n-dimensional quasi-conformal maps. J. Geom. Imaging Comput 1.4 (2014): 395-415. [21] Wong, Tsz Wai, et. al., Intrinsic feature extraction on hippocampal surfaces and its applications, SIAM Journal on Imaging Sciences 5.2 (2012): 746-768. [22] Lui, Lok Ming, et. al., Landmark constrained genus zero surface conformal mapping and its application to brain mapping research, Applied Numerical Mathematics 57.5-7 (2007): 847-858. [23] Lui, Lok Ming, et. al., Optimized conformal surface registration with shape-based landmark matching, SIAM Journal on Imaging Sciences 3.1 (2010): 52-78. [24] Lui, Lok Ming, et. al., Shape-based diffeomorphic registration on hippocampal surfaces using beltrami holomorphic flow, International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, Berlin, Heidelberg, 2010: 323-330. [25] Choi, Pui Tung, Ka Chun Lam, and Lok Ming Lui, FLASH: Fast landmark aligned spherical harmonic parameterization for genus-0 closed brain surfaces, SIAM Journal on Imaging Sciences 8.1 (2015): 67-94.
32
Hei Long Chan received the M.Phil. and B.Sc. degree from the Department of Mathematics, the Chinese University of Hong Kong, Hong Kong in 2017 and 2014 respectively. He is currently a PhD student in the Chinese University of Hong Kong, Department of Mathematics. His research interests include computational quasi-conformal geometry, medical imaging, shape analysis, image processing and machine learning. Tsz Chun Yam received the M.Phil. and B.Sc. degree from the Department of Mathematics, the Chinese University of Hong Kong, Hong Kong in 2016 and 2014 respectively. His research interests include computational quasiconformal geometry, medical imaging and machine learning.
Lok Ming Lui received the Ph.D. degree in applied mathematics from the University of California at Los Angles, Los Angeles, CA, USA, in 2008. He is an Assistant Professor with the Department of Mathematics, Chinese University of Hong Kong (CUHK), Hong Kong. Before joining CUHK, he was a Post-Doctoral Scholar for 2 years with the Department of Mathematics, Harvard University, Cambridge, MA, USA. His current research interests include computational conformal and quasiconformal geometry, Teichmuller theory, surface registration, medical imaging, and shape analysis. 1
Preprint submitted to Pattern Recognition
April 11, 2019
Automatic Characteristic-Calibrated Registration (ACC-REG): Hippocampal Surface Registration using Eigen-graphs Hei Long CHANa , Tsz Chun YAMa , Lok Ming LUIb,∗ a
222B, Lady Shaw Building, The Chinese University of Hong Kong, Shatin, Hong Kong 207, Lady Shaw Building, The Chinese University of Hong Kong, Shatin, Hong Kong
b
This work is supported by HKRGC GRF (Project ID: 14303414). Data used in preparation of this article were obtained from the Alzheimers Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found here: ADNI Acknowledgement List
∗
Corresponding author Email addresses: [email protected] (Hei Long CHAN), [email protected] (Tsz Chun YAM), [email protected] (Lok Ming LUI) Preprint submitted to Pattern Recognition
November 4, 2019