MRI based brain morphometry analysis has gained extensive interest in the past decades . A lot of research works are focused on identifying very early signs of brain functional and structural changes for early identification and prevention of neurodegenerative diseases. Alzheimer’s disease (AD), which is the sixth-leading cause of death in the United States, and the fifth-leading cause of death among those age 65 and older as reported by Alzheimer’s Association in 2018 , has received much interest from researchers around the world. Early detection and prevention of AD can significantly impact treatment options, improve quality of life, and save considerable health care costs. As a non-invasive method, brain imaging study has great potential that will powerfully track disease progression and therapeutic efficacy in AD. For example, whole brain morphometry, hippocampal and entorhinal cortex volumes are among most promising candidate biomarkers in structural MRI analysis. However, missing at this time is a widely available, highly objective brain imaging biomarker capable of identifying abnormal degrees of cerebral atrophy and accelerated rate of atrophy progression in preclinical individuals at high risk for AD in who early intervention is most needed.
Computational geometric methods are widely used in medical imaging fields including virtual colonoscopy and brain morphometry analysis. Rooted in deep geometry analysis research, computational geometric methods may provide rigorous and accurate quantification of abnormal brain development and thus hold a potential to detect preclinical AD in presymptomatic subjects. Specifically, surface morphometry techniques, such as conformal mapping and area preserving mapping, have shown to be feasible and powerful tools in brain morphometry research.
To the best of our knowledge, this is the first work to propose the use of the surface foliation theory for brain morphometry analysis. We validate our method by classifying brain surfaces of patients with Alzheimer disease and healthy control subjects. Experimental results indicate the efficiency and efficacy of our proposed method. The main contributions are summarized as follows:
A novel brain surface morphometry analysis method is proposed based on surface foliation theory.
A set of new geometric features (i.e., heights and circumferences) computed by pants decomposition and conformal mapping of topological cylinders are also proposed for surface indexing and classification.
The proposed method is rigorous, geometric and automatic.
Ii Previous Works
Brain morphometry analysis plays a fundamental role in medical imaging . Many research works have investigated the brain morphometry analysis and shape classification. Winkler et al.  proposed that the surface area could serve as an important morphometry feature to study brain structural MRI images. Besides, numerous methods have been presented in order to describe shapes, including statistical methods , topology based methods , and geometry based methods . To solve real 3D shape problems, researchers have also proposed many shape analysis and classification methods. Zacharaki et al.  proposed the use of pattern classification methods for classifying different types of brain tumors. Recently, Su et al.  presented a shape classification method using Wasserstein distance. The method computed a unique optimal mass transport map between two measures, and used Wasserstein distance to intrinsically measure the dissimilarities between shapes.
is a generalization of vector field. In computer graphics field, Zhang et al. invented a vector field design system which could help users create various vector fields with control over vector field topology. The technique can be used in some applications such as example-based texture synthesis, painterly rendering of images, and pencil sketch illustrations of smooth surfaces. Recently, Campen et al.  proposed a method for bijective parametrization of 2D and 3D objects based on simplicial foliations. The method decomposed a mesh into one-dimensional submanifolds, reducing the mapping problem to parametrization of a lower-dimensional manifold. It was proved that the resulting maps are bijective and continuous. In isogeometric analysis field, Lei et al.  presented a novel quadrilateral and hexahedral mesh generation method using foliation theory. A colorable quad-mesh method was employed to generate the quadrilateral mesh based on Strebel differentials, which then leads to the structured hexahedral mesh of the enclosed volume for high genus surfaces.
Iii Theoretic Foundation
We briefly review the basic concepts and theorems in conformal geometry. Detailed treatments can be found in .
If a complex function satisfies the Cauchy-Riemann equation
then is called a holomorphic function. If is invertible, and is also holomorphic, then is a bi-holomorphic function. A surface with a complex atlas , such that all chart transition functions are bi-holomorphic, then it is called a Riemann surface, and the atlas is called a complex structure. All oriented metric surfaces are Riemann surfaces.
Definition 1 (Holomorphic Quadratic Differentials
Let be a Riemann surface, and a complex differential form, such that on each local chart with the local complex parameter , where is a holomorphic function. The is called a holomorphic quadratic differential.
According to Riemann-Roch Theorem, the linear space of all holomorphic quadratic differentials is complex dimensional, where the genus . A point is called a zero of , if vanishes. A holomorphic quadratic differential has zeros. For any point away from zero, we can define a local coordinates
which is the so-called natural coordinates induced by . The curves with constant real (imaginary) natural coordinates are called the vertical (horizontal) trajectories. The trajectories through the zeros are called the critical trajectories.
Definition 2 (Strebel)
Given a holomorphic quadratic differential on a Riemann surface , if all of its horizontal trajectories are finite, then is called a Strebel differential.
A holomorphic quadratic differential is Strebel, if and only if its critical horizontal trajectories form a finite graph . The Strebel differentials are dense in the space of all holomorphic quadratic differentials. Given a holomorphic quadratic differential , the natural coordinates in Eqn. 1 induce a flat metric with cone singularities (cone angles equal to ), which is denoted as . Hubbard and Masur proved the following existence of a Strebel differential with prescribed type and heights.
Theorem 1 (Hubbard and Masur)
Given non intersecting simple loops , and positive numbers , , there exists a unique holomorphic quadratic differential satisfying:
The critical graph of partitions the surface into cylinders, , such that is the generator of ,
The height of each cylinder equals to , .
The geometric interpretation of above theorem is as follows: each cylinder becomes a canonical flat cylinder under , whose height is . Therefore, Strebel’s theorem allows one to specify the type of and the height of each cylinder .
Harmonic Map Suppose is a graph, is an edge weight function. Let and be two points on the graph, and be the shortest distance between them. Let be a surface with a Riemannian metric . Consider a map . We say a point is a regular point, if its image is not any node of , otherwise a critical point. The set of all critical points is denoted as . For each regular point , we can find a neighborhood , and the restriction of the map on can be treated as a normal function . We choose an isothermal coordinates on , such that the metric has a special form . The harmonic energy is given by , where , and the area element is . The harmonic energy of the whole map is defined as
The critical point of the harmonic energy is called a harmonic map. Wolf  proved the existence and the uniqueness of the harmonic map.
Theorem 2 (Wolf)
In each homotopy class, the harmonic map exists and is unique. Furthermore, the Hopf differential induced by the harmonic map is a holomorphic quadratic differential, where is the complex isothermal coordinates of .
Conformal Module Suppose is a surface of genus . Given non-intersecting simple loops and positive numbers , the unique Strebel differential in Hubbard and Masur’s theorem induces a flat metric with cone singularities, and cylinders . The height and circumference of each cylinder are (, ). The set of all (, ) are the conformal modules.
Pants Decomposition Let be a closed surface of genus , represented by triangular mesh. Let be a set of admissible curves, which can be generated automatically or manually specified. User also specifies a height parameter for each admissible curve . These admissible curves decompose surface to a set of pants . The pants decomposition graph is then constructed in the following way:
each pants corresponds to a node in
each admissible curve connecting two pants corresponds to an edge in ; two pants may be the same, in that case, the edge becomes a loop
Fig. 1 illustrates pants decomposition and pants decomposition graph.
Discrete Harmonic Map to Graph We compute a unique harmonic map from surface to . The harmonic energy is defined as
where is vertex on , on , is edge, and is cotangent weight.
For each , by moving to the barycenter of its neighbors on graph , the energy will decrease monotonically, which is due to the following definition of barycenter. By iteratively doing so, the energy will attain its minimum value, at which point we obtain a harmonic map . Thm. 2 guarantees this harmonic map we obtained is the unique one. Fig. 2 (a) and (b) illustrate harmonic map from a human face surface (a) to its pants decomposition graph (b), and Fig. 2 (c) shows the surface foliation, where color indicates vertices’ target position on graph .
The initial map should be specified in the same homotopy class as the final harmonic map . Subgraph at a node consists of the node and all edges connecting to it. Then initial map can be obtained automatically in the following way: each pants is mapped to the subgraph at node of , then all pants maps are glued together to obtain .
Barycenter Calculation For each , we move to the barycenter of its neighbors. Calculating barycenter is done by minimizing energy
where the right are exactly the terms in that involve . can be calculated piecewisely. Then minimization of above energy boils down to minimum calculation of a set of quadratic functions.
Surfaces with Boundaries For surfaces with boundaries, we can either double cover those surfaces to obtain a closed surface, or we can add boundaries to the set of admissible curves, and such curves correspond to open edges on . Computation of harmonic map remains the same.
Extract Geometric Features A holomorphic quadratic differential can be induced from the harmonic map we obtained. Tracing the critical trajectories of and slicing surface along them, we obtain a set of topological cylinders, each corresponding to an input admissible curve. The set of heights and circumferences of those cylinders are topological invariants, which we propsose to use as geometric features for classification problems in next section.
Data Preparation The dataset used in our experiments includes images from 30 patients with Alzheimer disease (AD) and 30 healthy control subjects. The structural MRI images were from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) . The brain cortical surfaces were reconstructed from the MRI images by FreeSurfer. Then, a set of ‘Core 6’ landmark curves, including the Central Sulcus (CeS), Anterior Half of the Superior Temporal Gyrus (aSTG), Sylvian Fissure (SF), Calcarine Sulcus (CaS), Medial Wall Ventral Segment, and Medial Wall Dorsal Segment, are automatically traced on each cortical surface using the Caret package . In Caret software, the PALS-B12 atlas is used to delineate the “core 6” landmarks, which are well-defined and geographically consistent, when compared with other gyral and sulcal features on human cortex. The stability and consistency of the six landmarks were validated in . An illustration of the landmark curves on a left cortical surface is shown in Fig. 3 with two views. We show the landmarks with both the original and inflated cortical surfaces for clarity. A brain surface and its foliation are shown in Fig. 4 (a) and (b), respectively.
The experimental procedures involving human subjects described in this paper were approved by the Institutional Review Board.
Foliation Feature Visualization We illustrate the difference of feature values between a pair of subjects with AD and healthy control subject (CTL) using radar chart. Radar chart displays multi-variate data in a two-dimensional chart where multiple variables are represented on axes starting from the same point. As shown in Fig. 5, six pairs of heights(H) and circumferences(C) corresponding to core 6 landmarks, i.e., twelve features (labeled by ’H1’, ’C1’,…,’H6’, ’C6’) are associated with twelve corners on the radar chart. We find that the pair of the H4 height and C4 circumference features associated to landmark curve of medial wall dorsal segment have the largest difference between these two subjects’ radar charts represented by a blue color line and an orange color line respectively. Although more validations are warranted, our research results may help discover AD related brain atrophy patterns.
Classification We validated our method with brain surface classification on a dataset of brain cortical surfaces from 30 patients with Alzheimer disease and 30 healthy control subjects. The SVM method was employed as the classifier with 10-fold cross validation in our experiments. For each image, the input feature vector of the classifier includes 12 features. For comparison purpose, we also computed cortical surface area and cortical surface mean curvatures, two cortical surface features frequently adopted in prior structural MRI analyses . We also applied SVM as the classifiers for these two features. Experimental results are shown in Table I. Our proposed method achieved 78.33% correctness rate, which is better than the correctness rate 56.67% in the brain surface area based method and 55.00% in the brain surface mean curvature based method. Although multi-subject studies are clearly necessary, this experiment demonstrates that the foliation theory based geometric features may have the potential to quantify and measure AD related cortical surface changes.
|Classification method||Correctness rate|
|Brain Surface Area||56.67%|
|Brain Mean Curvature||55.00%|
In this paper, a novel set of brain surface features are proposed based on surface foliation theory. To validate our method, we applied our method on classifying brain cortical surfaces between AD and healthy control subjects, and the preliminary experimental results demonstrated the efficiency and efficacy of our method. In future, we will employ our method to explore brain morphometry related to mild cognitive impairment (MCI) and other medical imaging applications.
-  A. Winkler et al., “Measuring and comparing brain cortical surface area and other areal quantities,” Neuroimage, vol. 61, no. 4, pp. 1428–1443, 2012.
-  A. Association et al., “2018 alzheimer’s disease facts and figures,” Alzheimer’s & Dementia, vol. 14, no. 3, pp. 367–429, 2018.
W. Zeng, R. Shi, Y. Wang, S.-T. Yau, X. Gu, A. D. N. Initiative, et al.,
“Teichmüller shape descriptor and its application to alzheimer?s disease
International journal of computer vision, vol. 105, no. 2, pp. 155–170, 2013.
-  R. Osada, T. Funkhouser, B. Chazelle, and D. Dobkin, “Shape distributions,” ACM Transactions on Graphics (TOG), vol. 21, no. 4, pp. 807–832, 2002.
M. Hilaga et al., “Topology matching for fully automatic similarity estimation of 3d shapes,” inProceedings of the 28th annual conference on Computer graphics and interactive techniques. ACM, 2001, pp. 203–212.
-  M. Mahmoudi and G. Sapiro, “Three-dimensional point cloud recognition via distributions of geometric distances,” Graphical Models, vol. 71, no. 1, pp. 22–31, 2009.
E. Zacharaki et al., “Classification of brain tumor type and grade using mri texture and shape in a machine learning scheme,”Magnetic Resonance in Medicine, vol. 62, no. 6, pp. 1609–1618, 2009.
-  Z. Su, W. Zeng, Y. Wang, Z.-L. Lu, and X. Gu, “Shape classification using wasserstein distance for brain morphometry analysis,” in International Conference on Information Processing in Medical Imaging. Springer, 2015, pp. 411–423.
-  K. Strebel, Quadratic differentials. Springer, 1984.
-  E. Zhang, K. Mischaikow, and G. Turk, “Vector field design on surfaces,” ACM Transactions on Graphics (ToG), vol. 25, no. 4, pp. 1294–1326, 2006.
-  M. Campen, C. Silva, and D. Zorin, “Bijective maps from simplicial foliations,” ACM Transactions on Graphics, to appear, vol. 35, no. 4, p. 7, 2016.
-  N. Lei, X. Zheng, J. Jiang, Y. Lin, and X. Gu, “Quadrilateral and hexahedral mesh generation based on surface foliation theory,” Computer Methods in Applied Mechanics and Engineering, 2016.
-  X. Gu and S.-T. Yau, Computational conformal geometry. International Press Somerville, Mass, USA, 2008.
-  M. Wolf, “On realizing measured foliations via quadratic differentials of harmonic maps tor-trees,” Journal D?Analyse Mathematique, vol. 68, no. 1, pp. 107–120, 1996.
-  S. Mueller, M. Weiner, L. Thal, R. Petersen, C. Jack, W. Jagust, J. Trojanowski, A. Toga, and L. Beckett, “The alzheimer’s disease neuroimaging initiative,” Neuroimaging Clinics of North America, vol. 15, no. 4, pp. 869–877, 2005.
-  D. Van Essen et al., “An integrated software suite for surface-based analyses of cerebral cortex,” J Am Med Inform Assoc, vol. 8, no. 5, pp. 443–459, 2001.
-  D. C. Van Essen, “A Population-Average, Landmark- and Surface-based (PALS) atlas of human cerebral cortex,” Neuroimage, vol. 28, no. 3, pp. 635–662, Nov 2005.
-  S. Li, X. Yuan, F. Pu, D. Li, Y. Fan, L. Wu, W. Chao, N. Chen, Y. He, and Y. Han, “Abnormal changes of multidimensional surface features using multivariate pattern classification in amnestic mild cognitive impairment patients,” J. Neurosci., vol. 34, no. 32, pp. 10 541–10 553, Aug 2014.