Landmark-based spherical quasi-conformal mapping for hippocampal surface registration

Introduction

Alzheimer’s disease (AD) (1) is a common neurodegenerative disease which seriously endangers the health of older adults. Dementia, a major medical and social problem that diminishes individuals’ ability to function in daily life, afflicts more than 25 million people worldwide (2). To delay the progression of AD at its early stage, it is necessary to establish a biomarker which can accurately quantify the effect of AD. Research into structural magnetic resonance imaging (sMRI) (3) has identified potential noninvasive neurodegeneration biomarkers of AD through this modality (4). Among these, hippocampal atrophy measures from sMRI are widely used, as hippocampal morphometry changes are apparent in the early stages of memory decline and may anticipate progression to mild cognitive impairment (MCI) and AD (5,6). However, due to the differences in patient age, patient gender, and the software used to extract hippocampal surfaces, it is necessary to register the hippocampal surface differences between individuals (7) to evaluate the effects of AD. By extracting the landmarks that reflect the essential morphological properties of the hippocampus and by establishing a spherical quasi-conformal mapping mechanism, a one-to-one morphology correspondence between different hippocampal surfaces can be established, and the effects of AD can then be assessed at the group level.

Surface registration refers to aligning and matching surface models from two or more different data sources with landmarks to maintain consistency in the same coordinate system, thus providing a foundation for subsequent data analysis, visualization, and simulation. Bijective registration can identify meaningful one-to-one correspondences between different surfaces. The landmark-based hippocampal surface registration algorithm involves aligning essential anatomical features (landmarks) from different individuals. These anatomical features are usually landmarks that can be manually marked by experts or automatically marked by certain algorithms. For instance, Wong et al. (8) proposed a model to automatically extract two intrinsic landmark curves according to the hippocampal surface morphology based on the Laplace-Beltrami (LB) operator on Riemannian manifolds. However, the calculation cost of this model was high because it involved the gradient descents of two energy functions, which typically require thousands or more iterations to converge numerically. In addition, due to the lack of consideration of large variations between individuals, the extracted eigen-graphs between individuals lacked correspondences. Later, Chan et al. (9) applied the calibration of the eigen-graphs to improve the accuracy of the eigen-graph correspondences. However, this method excessively pursued the smoothness of the landmark curves and neglected the intrinsic morphological features of hippocampal surfaces. Moreover, the determination of landmarks in the head and the tail of hippocampal surface has not been solved. Therefore, greater attention should be paid to the computational efficiency and essential morphological feature correspondences for extracting the landmarks on the hippocampal surfaces.

After the landmarks on hippocampal surfaces are extracted, the common procedure is to select a high-quality registration method that can perform diffeomorphism and landmark matching on the parametric domain. Shen et al. (10) proposed a surface registration method that involves rotating the main axes of the first-order ellipsoid (FOE) obtained from each hippocampal surface to coincide with the two designated coordinate axes. However, this method only focuses on the overall orientation registration of the hippocampus and ignores the alignment of local anatomical feature. Zou et al. (11) proposed a landmark-based registration method that could achieve alignment between anatomical features by minimizing the thin plate bending energy. This method, however, leads to the loss of bijectivity because it is difficult to strictly guarantee conformal mapping under the premise of a complete matching of landmarks. Recently, various registration methods have begun to incorporate quasi-conformal mapping rather than conformal mapping. Lam and Lui (12) proposed a quasi-conformal registration method that used the landmark information to obtain diffeomorphic registration between two natural 2D images. This method is a viable landmark-based registration method for dealing with large deformation in the 2D parametric plane, but it has not been applied to the registration of three-dimensional (3D) hippocampal surfaces.

In this study, we aimed to extract the intrinsic landmarks through spectrum analysis and to establish an eigen-graph correspondence via the calibration method. We then further resolved to use the local stereographic projection and quasi-conformal mapping method to complete the accurate alignment of landmarks while minimizing the local geometric distortion. This paper provides three major contributions to this field:

• A method to extract the stable and essential landmarks of hippocampal surface based on the spectrum analysis, which may serve as platform for further landmark matching;
• A method based on local forward and inverse stereographic projection that can achieve 3D surface registration by means of 2D quasi-conformal mapping;
• A method that leverages the mechanism of quasi-conformal mapping to accurately align landmarks while maintaining the bijectivity of mapping under a condition of large deformation through minimizing the compound energy function.

Methods

This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013), and the data used in preparation of this article were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu) (13). ADNI is a publicly available dataset that researchers from around the world can share, and it is the product of the joint efforts of many researchers from a wide range of academic institutions and private companies (for the latest information, please refer to www.adniinfo.org).

Image acquisition and data preprocessing

The experimental data in our study were downloaded from the ADNI site. Its research participants were recruited from over 50 locations in the United States and Canada. For each participant, the T1-weighted magnetic resonance imaging (MRI) scan was obtained using a sagittal 3D magnetization-prepared rapid gradient echo (MPRAGE) sequence with the following parameters: repetition time (TR) =2,400 ms, inversion time (TI) =1,000 ms, image matrix =256×256, and slice thickness =1 mm. Based on T1-weighted MRI scanning, the hippocampal substructure was segmented using FMRIB’s Integrated Registration and Segmentation Tool (FIRST) (14), and the hippocampal surface was automatically reconstructed (15). In order to overcome the noise caused by image scanning and partial volume effects, progressive mesh (16) and loop subdivision (17) were applied to smooth the reconstructed hippocampal surface. Through the above process, we could obtain the triangular meshes {V_i, F_i} of the hippocampal surfaces, where i represents the i-th hippocampus.

Landmark extraction

We used principal component analysis (PCA) (18) to generate the extreme poles of the hippocampal surface. Assuming that V(X,Y,Z) is a coordinate matrix of a hippocampal surface H, where n is the number of vertices of H, we could first obtain the first principal component A of V(X,Y,Z) according to eigenvalue decomposition, where A is the eigenvector corresponding to the maximum eigenvalue. We projected the original surface coordinates of H into the coordinate system determined by A and selected the two vertices with the farthest Euclidean distance as the two extreme vertices, which could be used as the North pole N and the South pole S, respectively, in the spherical parameter domain.

We then used the LB operator (19) to generate the eigen-graph of the hippocampal surface. According to the theory of differential geometry (20), the LB operator is defined as Δf=divgradf on a Riemannian manifold H, where f:H→ℝ. On compact manifolds, the LB operator’s first nontrivial eigenfunction satisfy the following:

−Δf=λf,λ>0[1]

The solutions of Eq. [1] represent the spatial part with an infinite number of eigenvalues λ and eigenfunction f pairs. Each f_i represents the vibration mode at a specific frequency λ_j where j represents the j-th frequency. Subsequently, the first eigenfunction f₁ according to the first nontrivial eigenvalue is defined as the eigen-graph of H (9). Because f₁ and λ₁ are dependent only on the Riemannian structure of the manifold, f₁ contains the intrinsic geometric information of H. It is worth noting that the value of f₁ monotonously increases from one tip to the other (8). The hippocampal surface can then be conceived of as a series of eigen-loops along the longitudinal direction which are constructed by the vertices with the similar value of f₁.

However, due to individual morphological differences of the obtained hippocampal surfaces, the correspondences between the individuals’ eigen-graphs might not be established. Therefore, it was necessary to calibrate the eigen-graphs of all the individuals according to the template surface. With the assumption that H₁ is a hippocampal surface, fH1 is its eigen-graph, H_T is the template hippocampal surface, and fHT is its eigen-graph, the calibration details of the eigen-graph of a participants according to the template’s eigen-graph could be obtained as follows. First, we normalized fH1 and fHT to f⌢H1 and f⌢HT on [0,1], respectively. The cumulative distribution functions of the two normalized eigen-graphs could then be calculated based on the following equation:

f⌢i*(x)=∫0xf⌢i(t)dt,i=H1,HT[2]

Owing to the monotonicity of f⌢i*(x), for all x₁ ∈ [0,1], there exists x₂ ∈ [0,1] such that f⌢H1*(x1)=f⌢HT*(x2). Based on the histogram matching theory (21), we could define the histogram matching function FH1HT:[0,1]→[0,1] by FH1HT(x1)=x2 for all x₂ ∈ [0,1] such that f⌢H1*(x1)=f⌢HT*(x2) and for all x₁ ∈ [0,1]. Finally, the calibrated eigen-graph of H₁ could be obtained as f⌢H1HT=f⌢H1⋅FH1HT. Here f⌢H1 and f⌢HT are the two normalized eigen-graphs of a participant and the template.

After the calibration procedure was determined, the eigen-graph correspondence between the individual and the template was greatly improved, providing a solid foundation for subsequent landmark matching. The results of eigen-graph calibration are shown in Figure 1. Figure 1A shows the uncalibrated individual eigen-graph (blue) and the template eigen-graph (red), while Figure 1B shows the calibrated individual eigen-graph (cyan) and the template eigen-graph (red). In Figure 1, the trend of the calibrated individual eigen-graph appears highly similar to that of the template, which is more conducive to subsequent registration.

Figure 1 The calibration of the eigen-graph. The horizontal axis represents the value of the eigen-graph, and the vertical axis represents the corresponding number of vertices. (A,B) The red histogram is the eigen-graph of template H_T; (A) the blue histogram is the eigen-graph of individual H₁, and (B) the cyan histogram is the calibrated eigen-graph of H₁.

Next, we aimed to locate two intrinsic feature curves which lay on the opposite sides of the most curved regions of H₁. We defined the intrinsic feature curves γ1,γ2 as follows:

{γ1,γ2}=argmax∫01|γ1(t)−γ2(t)|dt[3]

In discrete cases, we uniformly divided the calibrated eigen-graph f⌢H1HT∈[0,1] into m partitions, and then took the minimal interval ε on both sides of each segmentation value except for the two endpoints 0 and 1, which were considered to be the intersection points of two feature curves at both tips of the hippocampal surface. Each annular band with the width of 2ε was defined as an eigen-loop. Among each eigen-loop, the f⌢H1HT values on the corresponding triangular mesh vertices were approximately equal, and each eigen-loop intersected the two feature curves at two vertices, with the distance between these vertices being the farthest. We considered the two vertices to be a single landmark pair of each eigen-loop. First, we connected all the landmarks on each side of the hippocampal surface according to the Dijkstra algorithm (22) to obtain the two landmark curves {γ1,γ2}.

For exceptional cases, due to the large morphological noise interference on some hippocampal surfaces, some obtained landmark curves could not be stably located in the contralateral curvature area of the hippocampal surfaces, as shown in the representation of 3D space in Figure 2A. Therefore, the correction procedure for the landmark curves needed to be implemented. We corrected the landmark curves based on the angle constraints between the standard vector and the vector from each landmark pair. First, we selected the first eigen-loop LO₁ as the standard eigen-loop, and the standard vector s₁ pointed from v₁ to w₁ in LO₁, where v₁ and w₁ are the two landmarks of LO₁, respectively. To calibrate the landmarks on each eigen-loop, we used the vectors determined by the landmarks on the eigen-loop. Calculating the other vectors {si}i=2m constructed from each computed landmark pairs, with each vector pointing from γ₁ to γ₂, we finally determined whether each computed landmark pair needed to be adjusted according to Eq. [4]:

cos(s1,si)=s1⋅si|s1|⋅|si|[4]

Figure 2 Correction of the landmark pair on a RH in 3D space. (A) Two landmark curves that needed to be corrected. They were obviously rotated and crossed on the hippocampal surface. (B) The process of vector angle correction, where β₁ = β₂ and α = arc cos (s₁, s_i). (C) The new vector from v_new to w_new which were parallel to the standard vector s₁ after angle correction. (D) The corrected landmark curves. RH, right hippocampal surface; 3D, three-dimensional.

If arccos(s1,si)>ζ, the landmark pair in LO_i needs to be adjusted according to the procedure depicted in Figure 2A. We first calculated the geometric center vertices m_i according to the vertices in LO_i. Starting from the geometric center point m_i, a straight line was created along the direction of vector s₁ that intersected the eigen-loop at two new vertices {vnew,wnew} which were regarded as the adjusted landmark pair of LO_i. This procedure is presented in Figure 2B,2C. Figure 2D shows the result of the corrected landmark curves. Using this correction procedure, we traversed the entire hippocampal surface and acquired the new landmark pairs. This preserved the landmarks determined by the original geometric features of the hippocampus and corrected a few irregular landmarks (the size of was adjustable). Finally, we connected the optimized landmark pairs as the N and S points in the order of N→vnew1→vnew2→…→vnewm→S and the two extreme stably located in the contralateral curvature area of the hippocampal surfaces.

Equal-area spherical parameterization

With the obtained landmarks, we selected the spherical parameterization method based on the barycentric coordinate theory (23) and the triangular mesh optimization method (24) to achieve equal-area mapping which could facilitate a surface comparison on the common spherical parameter space. The details of spherical parameterization are as follows.

For the triangular mesh of the hippocampal surface H₁, we first constructed the initial weight matrix D∈ℝn×n as follows. Neig(v_i) was defined as the neighbor list of the vertices vi∉γ˜1, with each element d_ij in D being the distance between v_i and v_j such that d_ij = ||v_i - v_j||₂. The weights of the nonadjacent vertices were set to 0. We then normalized d_ij→d_ij* within each 1-ring to satisfy∑j∈Neig(vi)dij*=1. To lock the azimuth angles of all the vertices belonging to γ˜1 to 0, where γ˜1 is the curve on the curved side of the hippocampal surface, we set all the weight elements in the row where these vertices were located at 0. Finally, we obtained the final weight matrix as D’ = D* + I, where I is an identity matrix with the same size as D*, with D* being the set of dij*.

We then used the barycentric coordinate theory in the Cartesian coordinate system to achieve the initial spherical mapping for the hippocampal surface on a unit sphere. The principal idea of this method is that the center vertex of the 1-ring neighborhood is collinear with the convex combination vertex of the 1-ring neighborhood and the center of the sphere (as shown in Figure 3). Thus, the following equations were used for the vertices with unknown positions on the sphere:

{xi2+yi2+zi2=1αixi=Dij'∑j=1n−txj+Dik'∑k=1txkαiyi=Dij'∑j=1n−tyj+Dik'∑k=1tykαizi=Dij'∑j=1n−tzj+Dik'∑k=1tzk[5]

Figure 3 Geometric interpretation of Eq. [5]. The vector between the vertex v_o and the weighted average v_c of its neighbors (v₁, v₂, v₃, v₄, v₅, v₆) is collinear with the vector between the vertex v_c and the sphere’s center R_o.

where n is the number of the vertices in H₁, i is the i-th, t=2 is the number of the fixed vertices (N and S), and α_i is a coefficient which extends or compresses the center vertex to the convex combination center (see Figure 3).

By solving the linear equation system in Eq. [5], we obtained the coordinates of all vertices in the Cartesian coordinate system. Subsequently, through the transformation from the Cartesian coordinates to spherical coordinates, we could obtain the initial spherical parameterized result as follows:

{r=x2+y2+z2θ=arccoszx2+y2+z2φ=arctan(yx)[6]

where r =1.

After obtaining the initial spherical parameterization Ψ, we further optimized Ψ using the triangular mesh optimization algorithm to acquire an equal-area mapping Ψ→Ψ˜ on a unit sphere. Due to the sparsity of the triangular mesh near the N and S poles caused by the initial spherical mapping, it was necessary to optimize the triangular mesh to reduce the area distortion of the mesh unit before and after the mapping. The procedure of the optimization algorithm for the triangular meshes on spherical surfaces was as follows. First, we calculated the distance g_k = ||v_i − v_j||₂ between the adjacent vertices v_i and v_j belonging to γ˜1 in the order from the S point to the N point and accumulated them in sequence to obtain the accumulated distance Gkc=∑k=1ngk, where n is the number of the landmarks on γ˜1. Specifically, we obtained g₁ =0 for the S point. Second, we normalized Gkc to obtain Gkc′:[0,1] and used latk=π(1−Gkc') to calculate the new latitudes of each landmark on γ˜1 while keeping the longitude of each landmark unchanged. We then modified the original latitudes of the vertices on the 1-ring of the N point to the new latitudes as follows:

latinew=π−(π−latn−1)Gn−1c′gNi[7]

where g_Ni is the distance between the N point and each vertex v_i on the 1-ring of the N point.

Similarly, we modified the original latitudes of the vertices on the 1-ring of the S point to the new latitudes as follows:

latjnew=π−(π−lat2)G2c′gSj[8]

where g_Sj is the distance between the S point and each vertex v_j on the 1-ring of the S point. After we obtained the new positions of all the landmarks and the vertices on the 1-ring of the N and S points as the known condition, we used the displacement function (25) to move the remaining vertices to the new positions and thus achieve the goal of optimizing the triangular mesh on the unit sphere. Finally, the spherical parameterized result Ψ˜ became more uniform, and there was less area distortion between the adjusted triangular mesh of the spherical surface and the triangular mesh of the original hippocampal surface.

Landmark-based spherical registration

We used local stereographic projection and the quasi-conformal mapping principle (26) to obtain accurate landmark-based registration with bijectivity and smoothness by minimizing the energy function, including landmark mismatch energy term and Beltrami coefficient terms (27). Quasi-conformal mapping f:ℂ→ℂ could map an infinitesimal circle to an infinitesimal ellipse. The eccentricity K=ab=1+|μ(z)|1−|μ(z)| was bounded, where a is the major axis of the ellipse, and b is the minor axis of the ellipse. f satisfies the following Beltrami equation:

∂f(z)∂z¯=μ(z)∂f(z)∂z[9]

where µ is the Beltrami coefficient which can effectively control the bijectivity of the registration. With {Vi}i=1n and {Wi}i=1n assumed to be the landmarks of H₁ and H_T, respectively. Our goal was to find a mapping f: H₁→H_T satisfying f (V_i) = W_i (i =1, 2, …, n) while minimizing the topology distortion. We first projected the spherical parameterized result described in section of “Equal-area spherical parameterization” onto a 2D plane based on the local stereographic projection and then minimized the compound energy, including the landmark mismatch energy and topology distortion term, via quasi-conformal mapping.

We first used North pole stereographic projection F_N to project the local spherical parameterized result Ψ˜H1(z<ρ) and Ψ˜HT(z<ρ) to the 2D plane, respectively, where ρ controls the local size of the projection. We minimized a compound energy to obtain the optimized quasi-conformal mapping f as follows:

f=argminf:H1→HT{∫H1|∇μ(f)|2+α∫H1|μ(f)|2}[10]

subject to f (F_N (V_i)) = F_N (W_i) (i =1, 2, …, n) and ‖μ(f)‖<1, where α is a weight parameter which can adjust the weight between the two energy terms. The first term of Eq. [10] ensured the smoothness of the mapping f, while the second term controlled the conformal distortion of the mapping f. Subsequently, we took ν: µ (f) as the variable of the optimization of f. We could use penalty splitting method (28) to solve the above optimization problem:

E⌣(v,f)=∫H1|∇v|2+α∫H1|v|2+σ∫H1|v−μ(f)|2[11]

subject to the constraints that f (F_N (V_i)) = F_N (W_i) (i = 1, 2, …, n) and ‖v‖∞<1, where σ is a penalty parameter. For a large enough σ, ν was close to µ(f), and Eq. [11] approximated the solution of Eq. [10]. At the n-th iteration, we used the energy descent method (29) to obtain the optimal Beltrami coefficient ν_n+1 by minimizing E⌣(v,fn) while fixing f_n. Then, the ν_n+1 was fixed, and we used the linear Beltrami solver (LBS) (30) to obtain the landmark-matching quasi-conformal mapping f_n+1 whose Beltrami coefficient µ(f_n+1) was similar to ν_n+1. By continuing to iterate this process until ‖vn+1−vn‖∞<ξ where ξ is a small real number, we could obtain f_n such that it satisfied the requirement of maintaining bijectivity before and after mapping while aligning the landmarks. Finally, we used the inverse North pole stereographic projection to obtain the local registration result Φ₁ in the local spherical parameterization domain.

Following the above procedure, we used the South Pole stereographic projection F_S to project the local spherical parameterized result Φ1(z>−ρ) and Ψ˜HT(z>−ρ) to the 2D plane, respectively. Eq. [11] and LBS were applied to register another part of the spherical parameterization domain. After spherical landmark alignment, we used spherical interpolation (31) to reconstruct the hippocampus with 2,562 vertices, and thus all registered hippocampal surfaces had the same number of vertices and one-to-one correspondence to the morphological characteristics. The whole procedure is summarized in Algorithm 1.

Results

In this section, we described experimental results and compare our method with other registration algorithms.

Landmark extraction

In the generation of the eigen-loops, we set m as 50 and ε as 0.02 after many experiments. As shown in Figure 4, we used a color bar to display the eigen-graph changes in the left hippocampal surface (LH) and the right hippocampal surface (RH).

Figure 4 Eigen-loops display. The value of the eigen-graph increased from the tail to the head of the LH (A) and RH (B). LH, left hippocampal surface; RH, right hippocampal surface.

During the correction procedure for the landmark curves, we set ζ as 20°. To better display the extracted landmark results, we connected the extracted landmarks into the landmark curves, as shown in Figure 5. Each group was divided into upper and lower parts. The upper were the original surfaces, and the lower were hippocampal surfaces with the extracted landmark curves. For each individual hippocampal surfaces, the RH and the LH were displayed in the left and right columns, respectively. The extracted landmark curve γ˜1 was marked in red, and the extracted landmark curve γ˜2 was marked in blue. The number of eigen-loops for each group of individuals was 50. There were 102 landmarks including N and S. The two landmark curves could well establish the correspondence of the significant geometric features between the hippocampal surfaces of the different participants.

Figure 5 The extracted landmark curves on the hippocampal surfaces of the different groups in the 3D space. There were 50 eigen-loops in each hippocampus. (A-C) The upper parts of the original hippocampal surfaces of the NC, MCI, and AD groups, respectively, while the lower parts are the extracted landmark curves of the NC, MCI, and AD groups, respectively. 3D, three-dimensional; NC, normal cognition; MCI, mild cognitive impairment; AD, Alzheimer’s disease.

Spherical parameterization results

In this section, we demonstrated the results of spherical parameterization and its optimization, as shown in Figure 6. The coordinates of the adjacent points of the North and South poles appear marked in red in Figure 6A,6D. Figure 6E,6H show the coordinates (marked in red) after triangular mesh optimization. Figure 6B,6C show the initial spherical parameterization results, and Figure 6F,6G show the results after spherical optimization. It is obvious that the landmarks on landmark curve γ˜1 were located to the spherical coordinate with φ = 0, and the optimized spherical triangular mesh after spherical optimization was more uniform. This might also be explained by the fact that each triangle area on the object surface could be treated equally by assigning the same amount of spherical parameterization space (32) through the equal-area spherical mapping method. Thus, we achieved equal-area spherical mapping which could facilitate the surface comparison on the common parameter space by mapping the closed hippocampal surface onto a unit sphere as per barycentric coordinate theory and triangular mesh optimization.

Figure 6 Spherical parameterization and optimization results. (A) A detailed view of the North pole of the initial spherical parameterization. (B,C) The results of the initial spherical parameterization. (D) A detailed view of the South pole of the initial spherical parameterization. (E) A detailed view of the North pole of triangular mesh optimization. (F,G) The results of triangular mesh optimization. (H) A detailed view of the South pole of triangular mesh optimization.

To measure the area distortion induced by the spherical parameterization, we defined the area distortion index (ADI) as follows:

ADI=12n−4∑i=12n−4|Fi−F′i|[12]

where n is the number of vertices in the hippocampal spherical parameterization, F_i is the area of the i-th triangle before spherical parameterization, and F_i^’ is the corresponding area of F_i after spherical parameterization. The smaller the ADI was, the smaller the area of distortion. We also calculated the angle distortion between the hippocampal and spherical parameterization as follows:

Eangle=(Angle(Φ2)−Angle(H))∗π/180[13]

where Φ2 is the spherical parameterization, Angle(Φ2) is the angle of Φ2, and Angle(H) is the angle of hippocampus H.

We compared our spherical parameterization method with two other spherical parameterization methods through examining 145 with Aβ+ patients with AD and 249 Aβ– normal controls (NCs) (see Table 1). A positive Aβ reading was considered dementia induced by AD. The comparison results are shown in Table 2. Conformal spherical mapping (CSM) (33) is a proposed of method of iteratively eliminating folding using weighted Laplacian–Beltrami feature projection and improving the time efficiency of the algorithm of a genus-0 closed surface. The large-scale surface modeling (LSM) (34) is based on the Fourier transform and represents the 3D surface as a linear combination of a set of complex functions, thereby improving the accuracy of spherical parameterization. As shown in Table 2, our method had the smallest ADI for these individuals, which indicates that our spherical mapping method could better maintain equal-area mapping. Table 2 shows that our algorithm was not optimal in E_angle (rad), as it was based on quasi-conformal mapping which allowed for angular distortion.

Landmark-based spherical registration

During the landmark alignment process, we set α to 0.1, the penalty parameter σ was set to 10, and ξ was set to 0.01. In the local stereographic projection, we set ρ to 0.30. To clearly observe the registration process, we selected 10 landmark pairs for landmark alignment. The entire procedure of the landmark-based spherical registration is shown in Figure 7. In this figure, RH is used as an example to show the whole landmark-based spherical registration process. The red pentagrams inside Figure 7B-7H are the landmarks of H_T, and the blue dots are the landmarks of H₁ to be registered. After the landmark alignment, the blue dots moved to the red pentagram position. The registration process of the LH was the same as that of RH, except with the RH template being replaced by the LH template.

Figure 7 Illustration of the landmark-based RH spherical registration. (A) Landmark extraction. (B) Spherical parameterization. (C) North pole stereographic projection. (D) Landmark alignment. (E) Inverse North pole stereographic projection. (F) South pole stereographic projection. (G) Repeat of landmark alignment. (H) Inverse South pole stereographic projection. (K) Hippocampal surface reconstruction. The red stars are the landmarks of H_T, and the blue dots are the landmarks of H₁. RH, right hippocampal surface.

To compare the performance of our method with other registration methods, we defined the landmark matching errors as follows:

Eland=1k∑i=1k‖Φ2(Vi)−Wi‖2[14]

where Φ2 is the landmark-based spherical registration function described in Algorithm 1, k is the number of landmarks on the hippocampal surface, V_i is the landmark of the individual hippocampal surface on parametric sphere, and Wi is the landmark of the template hippocampal surface. The smaller E_land was, the better the landmark alignment was.

The original deformation degree of landmarks between the individual hippocampal surface and the template hippocampal surface were defined as follows:

Ds=1k∑i=1k|Vi−Wi|2[15]

We compared our method with other registration algorithms in terms of the E_land and the topology preservation performance during the registration process with mean |µ| and overlaps, where |µ| is the calculated result between the registered spherical surface and the original hippocampal surface. The results are summarized in Table 3, where D_s is 1.0069±0.3150 in the AD group and 0.8822±0.1792 in the NC group. The spherical harmonics by first-order ellipsoid (SPHARM-FOE) method (10) is not a traditional landmark-based registration algorithm, and its registers the individual surfaces by adjusting the long and short axis of the FOE corresponding to the hippocampal surface into a specified position. Although this method does not involve the process of landmark alignment, we could use the extracted landmarks using our method to investigate the performance of this method in landmark alignment.

The large-deformation diffeomorphic metric mapping (LDDMM) (35) optimizes the energy function on the manifold to achieve optimal matching between two images during deformation. We combined it with the thin-plate spline (TPS) algorithm (36) in registration to obtain a landmark guided LDDMM algorithm. The spherical TPS (STPS) method (11) is a landmark-based registration algorithm that achieves landmark alignment by minimizing the thin-plate bending energy. As seem in Table 2, the landmark matching error of our method in the AD and NC groups was smaller than that of other algorithms, and the mean |µ| conformed to the characteristics of quasi-conformal mapping in large deformation, indicating that our method had a more accurate correspondence with the landmarks in the presence of large deformation. The overlaps of the algorithms were 0, indicating that the registration algorithm could ensure bijectivity. These findings indicated that our method could better maintain hippocampal surface alignment and conformality compared to the other registration algorithms.

To visually evaluate the match between the participants and the template hippocampus, we used the checkerboard pattern to validate the registration algorithm, the results of which are shown in Figure 8. Figure 8 shows that after alignment of the landmarks, the spherical mesh had no overlap with the triangular mesh but did have angle distortion which was due to the quasi-conformal mapping that allowed for a small amount of angle distortion, indicating that our algorithm could maintain bijectivity.

Figure 8 Checkerboard pattern texture with landmark alignment. SPHARM-FOE, spherical harmonic by first-order ellipsoid; LDDMM, large deformation diffeomorphic metric mapping; STPS, spherical thin-plate splines.

AD-induced regions of interest (ROIs) of the hippocampus

To verify whether our surface registration algorithm could establish a one-to-one morphological correspondence between the individuals, we extracted the ROIs of the hippocampus to see whether the hippocampal CA1 subregion matched based on our registration method.

First, we uniformly divided f⌢H1HT∈[0,1] into n partitions with a sufficiently small width and then defined the geometric center line of the hippocampus as the center vertices of each partition connected in sequence. Among each partition, the thickness of each vertex on hippocampal surface was defined as the distance between each vertex and the center point (37). The feature extraction results are shown in Figure 9.

Figure 9 Geometric centerline and thickness measures of the LH and RH. (A) The geometric centerline of the LH. (B) The geometric centerline of the RH. (C) Color coding on each vertex of the LH. (D) Color coding on each vertex of the RH. The unit of the thickness is expressed in millimeters. LH, left hippocampal surface; RH, right hippocampal surface.

Second, we generated the ROIs (38) from the Aβ+ AD and Aβ− NC groups based on statistical group difference analysis, which included 145 Aβ+ patients with AD and 249 Aβ− NCs (see Table 1). Moreover, considering the influence of three factors of age, gender, and group, we used the linear model in the SurfStat software package (http://www.math.mcgill.ca/keith/surfstat) to obtain the intrinsic thickness features of each vertex. We used the chi-square test and one-way analysis of variance (ANOVA) (39) to compare the demographic and clinical data. It could be seen that the age and gender factors of the two groups were matched, and there was a significant difference in the Mini-Mental State Examination (MMSE) between the two groups.

Next, the ROIs of the LH and the RH could be generated from the vertices with the group thickness difference with a P value <0.0001 based on the permutation t-test with 5,000 random permutations. The results are shown in Figure 10, in which the ROIs (red regions) indicate that the vertices of the region had statistical differences (P<0.0001) according to a 5,000-permutation t-test uncorrected for multiple comparisons. The results also show that each extracted ROI of the LH was larger than that of the RH based on the three registration methods, which indicated that the LH was more sensitive to AD than was the RH. These findings are consistent with the previous studies (40,41). We also drew a heatmap of the P values of each vertex in the template hippocampus, as shown in Figure 11: the red areas (Figure 11A,11B) are the ROIs of the LH and RH, respectively, and the heatmaps (Figure 11C,11D) correspond to the P value of the LH and RH, respectively; the red area represents vertices with P values less than 0.0001, the orange area represents vertices with P values less than 0.001, the yellow area represents the vertices with P values less than 0.01, and the blue areas represent the remaining vertices with other P values.

Figure 10 The extracted ROIs of our method, SPHARM-FOE, STPS, and LDDMM. (A,E) The ROIs extracted with our method. (B,F) The ROIs extracted with SPHARM-FOE. (C,G) The ROIs extracted with STPS. (D,H) The ROIs extracted with LDDMM. (A-D) The LH. (E-H) The RH. The ROIs (red regions) indicated that the vertices of the region had statistical differences (P value <0.0001) after a 5,000-permutation t-test. ROI, region of interest; SPHARM-FOE, spherical harmonic by first-order ellipsoid; STPS, spherical thin-plate splines; LDDMM, large deformation diffeomorphic metric mapping; LH, left hippocampal surface; RH, right hippocampal surface.

Figure 11 The P value heatmap of ROIs. (A) The red areas are the LH ROIs. (B) The red areas are the RH ROIs. (C) The P value heatmap of the LH. (D) The P value heatmap of the RH. ROI, region of interest; LH, left hippocampal surface; RH, right hippocampal surface.

Additionally, we found that only the ROIs based on our registration method matched the CA1 subregion of the hippocampus that was considered as the essential part for the majority of hippocampus-dependent memories (42). This also verified that our registration method might accurately establish a one-to-one morphological correspondence between the individuals. The variables described in this section are available in the Appendix 1.

Discussion

From the framework of our registration algorithm, it could be seen that the extraction of essential landmarks on the hippocampal surface directly affected the accuracy of registration. In this section, we examine the impact of the two steps of landmark extraction—i.e., eigen-graph calibration and eigen-graph segmentation—on the subsequent registration results.

Effects of eigen-graph calibration

In section of “Landmark extraction”, we describe the calibration of the eigen-graph of individuals according to the eigen-graph of the template based on histogram matching theory. To verify whether eigen-graph calibration is beneficial for improving registration accuracy, we compared the classification results of the thickness features obtained after registration with and without the eigen-graph calibration step. We used support vector machine (SVM) (43) to perform binary classification on the thickness measurements of 100 patients with AD and 230 NCs, as shown in Table 4. The parameter settings for SVM were as follows: fivefold cross validation, a rough Gaussian kernel function, and a box constraint level of 1. The classification results were in Table 5, which shows that the classification results with the eigen-graph calibration step had better a classification performance than did those without the eigen-graph calibration step.

In addition, we used the uniform manifold approximation and projection (UMAP) (44) algorithm to visualize the categories of AD and NC participants in a 2D embedding. Figure 12A,12B show the embedding results without and with the calibration, respectively. Through visual inspection, we could observe that the spacing between different categories increased significantly when the eigen-graph calibration step was used. This indicates that the eigen-graph calibration can help to improve registration accuracy, thereby enhancing the effectiveness of morphology comparison between individuals.

Figure 12 The 2D UMAP embedding (A) without eigen-graph calibration and (B) with eigen-graph calibration. The blue dots represent the patients with AD, and the red triangles represent the NCs. UMAP, uniform manifold approximation and projection; AD, Alzheimer disease; NC, normal control; 2D, two-dimensional.

Effect of eigen-graph segmentation

In section of “Landmark extraction”, we describe the extraction of the landmark pairs in each eigen-loop with the width of 2ε. A small width ε might result in fewer vertices being included in the eigen-loop. Thus, the extracted landmark pair might not be including the two farthest points in the Euclidean distance. However, a large width ε might result in each eigen-loop containing more vertices. Thus, the extracted landmark pair did not have similar eigenvector value. To compare the impact of width ε on registration results, we used SVM to classify those with AD and NCs based on the different eigen-graph segmentations with different ε values. The classification results are shown in Table 6. According to the SVM classification results, a large or small width ε can affect the accuracy of our landmark-based surface registration. In this study, we selected a width ε equal to 0.02 based on practical experience.

Limitations

Our surface registration method can ensure the bijectivity and the smoothness in registration in the presence of large deformations and achieve the precise alignment of the landmarks. Nevertheless, our method still has a few limitations to consider. First, our registration algorithm is only based on the landmarks and does not incorporate other surface morphology measures, such as the mean curvature and Gaussian curvature. If we combine these surface morphological measures with the landmarks, we can theoretically effectively improve the accuracy of the surface registration. Second, the data for extracting ROI were insufficient, as we only used a sample of 145 Aβ+ patients with AD and 249 Aβ− NC individuals, which is not sufficient to fully characterize the general pattern of morphological changes induced by AD symptoms. Nevertheless, our current research results indicate that our registration method is more sensitive to detecting hippocampal morphological changes compared to other registration algorithms.

Conclusions

In this paper, a landmark-based spherical quasi-conformal mapping for hippocampal surface registration algorithm was proposed which could maintain precise alignment of the landmarks and the bijectivity in the presence of large deformation. The experiment results indicated that our registration algorithm can effectively achieve an accurate one-to-one mapping between the different hippocampal surfaces under large deformations. In the future, we plan to incorporate other surface morphology measures into our method and improve the accuracy and effectiveness of surface registration.

Acknowledgments

Data used in preparation of this article were obtained from the 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 the ADNI investigators can be found online (http://adni.loni.usc.edu/wpcontent/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf).

Funding: This work was supported by the National Natural Science Foundation of China (No. 62171209).

Footnote

Conflicts of Interest: The authors have completed the ICMJE uniform disclosure form (available at https://qims.amegroups.com/article/view/10.21037/qims-23-1297/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Scheltens P, Blennow K, Breteler MM, de Strooper B, Frisoni GB, Salloway S, Van der Flier WM. Alzheimer's disease. Lancet 2016;388:505-17. [Crossref] [PubMed]
2. Trejo-Lopez JA, Yachnis AT, Prokop S. Neuropathology of Alzheimer's Disease. Neurotherapeutics 2022;19:173-85. [Crossref] [PubMed]
3. Wáng YXJ, Ma FZ. A tri-phasic relationship between T2 relaxation time and magnetic resonance imaging (MRI)-derived apparent diffusion coefficient (ADC). Quant Imaging Med Surg 2023;13:8873-80. [Crossref] [PubMed]
4. Jack CR Jr, Slomkowski M, Gracon S, Hoover TM, Felmlee JP, Stewart K, Xu Y, Shiung M, O'Brien PC, Cha R, Knopman D, Petersen RC. MRI as a biomarker of disease progression in a therapeutic trial of milameline for AD. Neurology 2003;60:253-60. [Crossref] [PubMed]
5. Qiu C, Kivipelto M, von Strauss E. Epidemiology of Alzheimer's disease: occurrence, determinants, and strategies toward intervention. Dialogues Clin Neurosci 2009;11:111-28. [Crossref] [PubMed]
6. Andersen E, Casteigne B, Chapman W D, Creed A, Foster F, Lapins A, Shatz R, Sawyer R. Diagnostic biomarkers in Alzheimer’s disease. Biomark Neuropsych 2021;5:100041.
7. Gerardin E, Chételat G, Chupin M, Cuingnet R, Desgranges B, Kim HS, Niethammer M, Dubois B, Lehéricy S, Garnero L, Eustache F, Colliot OAlzheimer's Disease Neuroimaging Initiative. Multidimensional classification of hippocampal shape features discriminates Alzheimer's disease and mild cognitive impairment from normal aging. Neuroimage 2009;47:1476-86. [Crossref] [PubMed]
8. Wong TW, Lui LM, Thompson PM, Chan TF. Intrinsic feature extraction on hippocampal surfaces and its applications. SIAM J Imaging Sci 2012;5:746-68.
9. Chan HL, Yam TC, Lui LM. Automatic characteristic-calibrated registration (ACC-REG): Hippocampal surface registration using eigen-graphs. Pattern Recognit 2020;103:107142.
10. Shen L, Huang H, Makedon F, Saykin AJ. Efficient registration of 3D SPHARM surfaces. Fourth Canadian Conference on Computer and Robot Vision (CRV '07), Montreal, QC, Canada, 2007:81-8.
11. Zou G, Hua J, Muzik O. Non-rigid surface registration using spherical thin-plate splines. Med Image Comput Comput Assist Interv 2007;10:367-74.
12. Lam KC, Lui LM. Landmark-and intensity-based registration with large deformations via quasi-conformal maps. SIAM J Imaging Sci 2014;7:2364-92.
13. Jack CR Jr, Bernstein MA, Fox NC, Thompson P, Alexander G, Harvey D, et al. The Alzheimer's Disease Neuroimaging Initiative (ADNI): MRI methods. J Magn Reson Imaging 2008;27:685-91. [Crossref] [PubMed]
14. Patenaude B, Smith SM, Kennedy DN, Jenkinson M. A Bayesian model of shape and appearance for subcortical brain segmentation. Neuroimage 2011;56:907-22. [Crossref] [PubMed]
15. Katsanikas M, Wiggins S. 3D generating surfaces in Hamiltonian systems with three degrees of freedom–I. Int J Bifurcation Chaos 2024;34:2430004.
16. Pajarola R, Rossignac J. Compressed progressive meshes. IEEE Trans Vis Comput Graph 2000;6:79-93.
17. Cheng FH, Fan FT, Lai SH, Huang CL, Wang JX, Yong JH. Loop subdivision surface based progressive interpolation. J Comput Sci Technol 2009;24:39-46.
18. Greenacre M, Groenen PJF, Hastie T, d’Enza AI, Markos A, Tuzhilina E. Principal component analysis. at Rev Methods Primers 2022;2:100.
19. RouveyrolM.Spectral estimate for the Laplace-Beltrami operator on the hyperbolic half-plane. arXiv:2401.14977. 2024.
20. Guggenheimer HW. Differential geometry. Courier Corporation. 2012.
21. Bottenus N, Byram BC, Hyun D. Histogram Matching for Visual Ultrasound Image Comparison. IEEE Trans Ultrason Ferroelectr Freq Control 2021;68:1487-95. [Crossref] [PubMed]
22. Luo M, Hou X, Yang J. Surface optimal path planning using an extended Dijkstra algorithm. IEEE Access 2020;8:147827-38.
23. Gotsman C, Gu X, Sheffer A. Fundamentals of spherical parameterization for 3D meshes. ACM Transactions on Graphics 2003;22:358-63.
24. Shen L, Kim S, Saykin AJ. Fourier method for large-scale surface modeling and registration. Comput Graph 2009;33:299-311. [Crossref] [PubMed]
25. Miao X, Wang Z, Zhu Y, Jiang Z, Dong L, Wu M. Local modified mesh deformation based on radial basis functions for fluid solid interaction in reactor core. Nucl Eng Des 2023;401:112076.
26. Zeng W, Luo F, Yau S T, Gu X. Surface quasi-conformal mapping by solving Beltrami equations. In: Hancock ER, Martin RR, Sabin MA. editors. Mathematics of Surfaces XIII. Mathematics of Surfaces 2009. Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 2009:391-408.
27. Gutlyanskiĭ V, Ryazanov V, Nesmelova O, Yakubov E. On the Hilbert problem for semi-linear Beltrami equations. J Math Sci 2023;270:428-48.
28. Wang P, Wang Y, Huang B, Wang L, Zhang X, Leung H, Chanussot J. Poissonian blurred hyperspectral imagery denoising based on variable splitting and penalty technique. IEEE Transactions on Geoscience and Remote Sensing 2023;61:1-14.
29. Zheng JH, Wu CQ, Li Z, Wang LX, Wu QH. A gradient descent direction based-cumulants method for probabilistic energy flow analysis of individual-based integrated energy systems. Energy 2023;265:126290.
30. Lui LM, Lam KC, Wong TW, Gu X. Texture map and video compression using Beltrami representation. SIAM J Imaging Sci 2013;6:1880-902.
31. Kim KH, Shim PS, Shin S. An alternative bilinear interpolation method between spherical grids. Atmosphere 2019;10:123-33.
32. Shen L, Farid H, McPeek M A. Modeling three‐dimensional morphological structures using spherical harmonics. Evolution 2009;63:1003-16. [Crossref] [PubMed]
33. Xue J, Li Q, Liu C, Xu M. Parameterization of triangulated surface meshes based on constraints of distortion energy optimization. J Vis Lang Comput 2018;46:53-62.
34. Chung MK, Dalton KM, Shen L, Evans AC, Davidson RJ. Weighted fourier series representation and its application to quantifying the amount of gray matter. IEEE Trans Med Imaging 2007;26:566-81. [PubMed]
35. Wu J, Tang X. A Large Deformation Diffeomorphic Framework for Fast Brain Image Registration via Parallel Computing and Optimization. Neuroinformatics 2020;18:251-66. [Crossref] [PubMed]
36. Rohr K, Stiehl HS, Sprengel R, Buzug TM, Weese J, Kuhn MH. Landmark-based elastic registration using approximating thin-plate splines. IEEE Trans Med Imaging 2001;20:526-34. [Crossref] [PubMed]
37. Burggren AC, Zeineh MM, Ekstrom AD, Braskie MN, Thompson PM, Small GW, Bookheimer SY. Reduced cortical thickness in hippocampal subregions among cognitively normal apolipoprotein E e4 carriers. Neuroimage 2008;41:1177-83. [Crossref] [PubMed]
38. Wang X, Wang Y. Multiple medical image encryption algorithm based on scrambling of region of interest and diffusion of odd-even interleaved points. Expert Syst Appl 2023;213:118924.
39. Langenberg B, Janczyk M, Koob V, Kliegl R, Mayer A. A tutorial on using the paired t test for power calculations in repeated measures ANOVA with interactions. Behav Res Methods 2023;55:2467-84. [Crossref] [PubMed]
40. Camargo A, Azuaje F, Wang H, Zheng H. Permutation - based statistical tests for multiple hypotheses. Source Code Biol Med 2008;3:15. [Crossref] [PubMed]
41. Foley SF, Tansey KE, Caseras X, Lancaster T, Bracht T, Parker G, Hall J, Williams J, Linden DE. Multimodal Brain Imaging Reveals Structural Differences in Alzheimer's Disease Polygenic Risk Carriers: A Study in Healthy Young Adults. Biol Psychiatry 2017;81:154-61. [Crossref] [PubMed]
42. Bird CM, Burgess N. The hippocampus and memory: insights from spatial processing. Nat Rev Neurosci 2008;9:182-94. [Crossref] [PubMed]
43. Roy A, Chakraborty S. Support vector machine in structural reliability analysis: A review. Reliab Eng Syst Saf 2023;233:109126.
44. Vermeulen M, Smith K, Eremin K, Rayner G, Walton M. Application of Uniform Manifold Approximation and Projection (UMAP) in spectral imaging of artworks. Spectrochim Acta A Mol Biomol Spectrosc 2021;252:119547. [Crossref] [PubMed]