# Full text of "Optimization of Surface Registrations using Beltrami Holomorphic Flow"

## See other formats

optimization of Surface Registrations using Beltrami Holomorphic Flow tL.M. Lui, T.W. Wong, W. Zeng, X.F. Gu, P.M. Thompson, T.F. Chan, S.T. Yau Abstract ^^ In shape analysis, finding an optimal 1-1 correspondence between surfaces within a large class of admissible ^^ bijective mappings is of great importance. Such process is called surface registration. The difficulty lies in the ^^ fact that the space of all surface diffeomorphisms is a complicated functional space, making exhaustive search JL* for the best mapping challenging. To tackle this problem, we propose a simple representation of bijective surface ^! maps using Beltrami coefficients (BCs), which are complex-valued functions defined on surfaces with supreme I norm less than 1. Fixing any 3 points on a pair of surfaces, there is a 1-1 correspondence between the set of ^— I surface diffeomorphisms between them and the set of BCs. Hence, every bijective surface map can be represented ,_, by a unique BC. Conversely, given a BC, we can reconstruct the unique surface map associated to it using the ^^ Beltrami Holomorphic flow (BHF) method. Using BCs to represent surface maps is advantageous because it is a much simpler functional space, which captures many essential features of a surface map. By adjusting BCs, we O r^ equivalently adjust surface diffeomorphisms to obtain the optimal map with desired properties. More specifically, a BHF gives us the variation of the associated map under the variation of BC. Using this, a variational problem over the space of surface diffeomorphisms can be easily reformulated into a variational problem over the space of BCs. This makes the minimization procedure much easier. More importantly, the diffeomorphic property is ^ always preserved. We test our method on synthetic examples and real medical applications. Experimental results ^1^ demonstrate the effectiveness of our proposed algorithm for surface registration. <N CO Index Terms in ^^ Beltrami coefficient, Beltrami holomorphic flow, surface diffeomorphism, surface registration, shape analysis. o optimization 1. Introduction d Surface registration is a process of finding an optimal 1-1 correspondence between surfaces satisfying certain constraints. It is of great importance in different research areas, such as computer graphics and medical imaging. For example, in medical imaging, surface registration is always needed for statistical shape analysis, morphometry and the processing of signals on brain surfaces (e.g., denoising and filtering). In many cases, a surface must be non-rigidly aligned with another surface, while matching various features lying on both surfaces. Finding an optimal surface registration that best matches the required constraints is difficult, especially on convoluted surfaces such as the human brain. It is therefore necessary to develop an effective algorithm to compute the best surface registration. In order to obtain the best 1-1 correspondence between two surfaces, an optimized surface registrations is often required. Optimization of surface registrations is the process of selecting an optimal surface ^ Lok Ming Lui, Department of Mathematics, Harvard University and UCLA{malmlui@ math.harvard.edu} 1 diffeomorphism within a large class of admissible smooth mappings to best satisfy certain properties. It can usually be formulated as a variational problem in the form: min Eoif) (1) where FDifr = {/ : S'l ^ 5^2 : / is a diffeomorphism} is the space of all diffeomorphisms from surface Si to surface 5^2. Solving this type of variational problem is generally difficult, since the space of all surface diffeomor- phisms Foiff is a complicated functional space. For instance, Foiff is inherently infinite dimensional and has no natural linear structure. Constructing an efficient optimization scheme in such space that guarantees to obtain a minimizer is a big challenge, and a loss of bijectivity of the surface maps (overlapping) is often observed during the optimization process. To solve this problem, it is necessary to develop a simple representation of surface diffeomorphisms which helps to simplify the optimization procedure. In this paper, we propose a simple representation of surface diffeomorphisms using Beltrami coefficients (BCs). The BCs are any complex- valued functions defined on surfaces with L^-norm strictly less than 1. Fixing any 3 points on a pair of surfaces, there is a one-to-one correspondence between the set of surface diffeomorphisms and the set of BCs. Hence, every bijective surface map can be represented by a unique BC. Conversely, given a BC, we propose to reconstruct the unique surface map associated to it using the Beltrami Holomorphic flow (BHF) method introduced in this paper. The BHF formulates the variation of the surface maps under the variation of BCs. Hence, variational problems of surface diffeomorphisms can be easily reformulated into variational problems of BCs in the form: min E(a) (2) where Fbc = {/x: S'l ^ D: ||/i||oo < 1} is the set of BCs. The space of BCs is a much simpler functional space that captures the essential features of surface maps. There are no restrictions that BCs have to be 1-1, surjective or satisfy some constraints in their Jacobians. By adjusting BCs, we can adjust surface registrations accordingly using BHF to obtain surface maps with the desired properties. This greatly simplifies the minimization procedure. More importantly, the surface maps obtained is guaranteed to be diffeomorphic (bijective and smooth) during the optimization process. We have applied our proposed algorithm on synthetic examples and real medical applications for surface registration, which demonstrate the effectiveness of our proposed method. In summary, our work contributes to the following three aspects: • We propose a simple representation of surface diffeomorphisms to facilitate the optimization of surface registrations. • We develop a reconstruction algorithm of the surface diffeomorphism from a given BC, using BHF. This completes the representation scheme and allows us to move back and forth between BCs and surface diffeomorphisms. • With BHF, we formulate variational problems of surface maps into variational problems of BCs. This greatly simplifies the optimization procedure. A flow chart summarizing the framework proposed in this paper is shown in Figure [T} Surface maps between surfaces S1 and S2 with the same topology Beltrami representation > Section 4.2 Space of BCs representating surface maps between S1 and S2 BHF optimization ^ Section 4.4 Optimal BC satisfying the desired properties for the surface registration BHF reconstruction > Section 4.3 Reconstructed surface registration satisfying the desired properties Fig. 1 A FLOW CHART SUMMARIZING THE FRAMEWORK PROPOSED IN THIS PAPER. 2. Previous Work Surface registration has been extensively studied by different research groups. Most methods compute the optimal surface registration by minimizing certain kinds of energy functional . In this section, we briefly describe some related methods commonly used. Conformal surface registration has been widely studied to obtain smooth 1-1 correspondence between surfaces and minimize angular distortions Oll7lll8lll9l lfT0lllT3l l ll20lll2^ Conformal maps are usually com- puted using variational approaches to minimize some energy functional, such as the harmonic energy [|7]| and the least square energy based on the Cauchy-Riemann equation [.20,1 . A 1-1 correspondence between surfaces can be obtained in the optimal state. However, the above registration cannot map anatomical features, such as sulcal landmarks, consistently from subject to subject. To obtain a surface registration that matches important landmark features, landmark-based diffeo- morphisms are often used. Optimization of surface diffeomorphisms by landmark matching has been extensively studied. Gu et al. \J2 improves a conformal parameterization by composing an optimal Mobius transformation so that it minimizes a landmark mismatch energy. The resulting parameterization remains conformal, although features cannot be perfectly matched. Wang et al. lfT8lll29l proposed a variational framework to compute an optimized conformal registration which aligns landmarks as well as possible. However, landmarks are not matched exactly and diffeomorphisms cannot be guaranteed when there is a a large amount of landmark features. Durrleman et al. f2\fi\ developed a framework using currents, a concept from differential geometry, to match landmarks within surfaces across subjects, for the purpose of inferring the variability of brain structure in an image database. Landmark curves are not perfectly matched. Tosun et al. Il26ll proposed a more automated mapping technique that attempts to align cortical sulci across subjects by combining parametric relaxation, iterative closest point registration, and inverse stereographic projection. Glaunes et. al [l6l [rT2l proposed to generate large deformation diffeomorphisms of a sphere onto itself, given the displacements of a finite set of template landmarks. The diffeomorphism obtained can better match landmark features. Lui et al. 11191 proposed to compute shape-based landmark matching registrations between brain surfaces using the integral flow method. The one parameter subgroup within the set of all diffeomorphisms is considered and represented by smooth vector fields. Landmarks can be perfectly matched and the correspondence between landmark curves is based on shape information. Leow et al. IITSll proposed a level- set-based approach for matching different types of features, including points, 2D and 3D curves represented as implicit functions. These matching features in the parameter domain were then pulled back onto surfaces to compute correspondence fields. Later, Shi et al. [l24l computed a direct harmonic mapping between two surfaces by embedding both surfaces as the level- set of an implicit function, and representing the mapping energy as a Dirichlet functional in 3D volume domains. Although such an approach can incorporate landmark constraints, it has not been proven to yield diffeomorphic mappings. Since there may not be well-defined landmarks on surfaces, some authors proposed driving features into correspondence based on shape information. Lyttelton et al. [21 J computed surface parameterizations that match surface curvature. Fischl et al. JU improved the alignment of cortical folding patterns by minimizing the mean squared difference between the average convexity across a set of subjects and that of the individual. Wang et al. [27J computed surface registrations that maximize the mutual information between mean curvature and conformal factor maps across subjects. Lord et al. lfT6l matched surfaces by minimizing the deviation from isometry. In most situations, one has to pay extra attention to ensure the optimal map computed is diffeomorphic. Hence, developing an effective optimization algorithm that guarantees to give diffeomorphic surface registrations is necessary. This motivates us to look for a simple representation of surface diffeomorphisms which helps to simplify the optimization procedure. 3. Theoretical Background In this section, we describe some basic mathematical concepts related to our algorithms. For details, we refer readers to [5 J and lfT4ll . A surface S with a conformal structure is called a Riemann surface. Given two Riemann surfaces M and A^, a map f : M ^ N h conformal if it preserves the surface metric up to a multiplicative factor called the conformal factor. An immediate consequence is that every conformal map preserves angles. With the angle-preserving property, a conformal map effectively preserves the local geometry of the surface. A generalization of conformal maps is quasi-conformal maps, which are orientation- preserving diffeomorphisms between Riemann surfaces with bounded conformality distortion, in the sense that their first order approximations takes small circles to small ellipses of bounded eccentricity [5J. Thus, a conformal homeomorphism that maps a small circle to a small circle can also be regarded as quasi- conformal. Figure [2] illustrates the idea of conformal and quasiconformal maps. Mathematically, / : C ^ C is quasi-conformal provided that it satisfies the Beltrami equation: df , .df Hippocampus Circle pattern on the parameter domain Conformal map (Circles to circles) (A) (B) Qiiasi-coyjformal map (Circles io ellipses) (C) (D) Fig. 2 Illustration of a conformal map and a quasiconformal map. (A) shows a hippocampal surface. A circle packing PATTERN IS PLOTTED ON THE PARAMETER DOMAIN AS IN (B). (C) SHOWS A CONFORMAL PARAMETERIZATION, WHICH MAPS CIRCLES ON THE PARAMETER DOMAIN TO CIRCLES ON THE SURFACE. (D) SHOWS A QUASICONFORMAL PARAMETERIZATION, WHICH MAPS CIRCLES ON THE PARAMETER DOMAIN TO ELLIPSES ON THE SURFACE. for some complex valued functions /x satisfying ||/x||oo < 1. In terms of the metric tensor, consider the effect of the puUback under / of the usual Euclidean metric ds\\ the resulting metric is given by: Qf ridsi) which, relative to the background Euclidean metric dz and dz, has eigenvalues (1 + |/x|)^ '^ dz\ (4) and (1 — |/i|)^||^| ./iis called the Beltrami coefficient, which is a measure of non-conformality. In particular, the map / is conformal around a small neighborhood ofp when /x(p) = 0. Infinitesimally, around a point p, f may be expressed with respect to its local parameter as follows: f{z) = f{p) + Mp)^ + Mp)z = f{p)+ fz{p){z + fl{p)z). (5) Obviously, / is not conformal if and only if /x(p) 7^ at p. Inside the local parameter domain, / may be considered as a map composed of a translation to f{p) together with a stretch map S{z) = z + iJ^{p)z, which is postcomposed by a multiplication of fz{p), which is conformal. All the conformal distortion of S{z) is caused by /x(p). S{z) is the map that causes / to map a small circle to a small ellipse. From /i(p), we can determine the angles of the directions of maximal magnification and shrinking and the amount of them as well. Specifically, the angle of maximal magnification is arg(/x(p))/2 with magnifying factor 1 + |/i(p)|; The angle of maximal shrinking is the orthogonal angle (arg(/x(p)) — 7r)/2 with shrinking factor 1 — |/x(p)|. The distortion or dilation is given by: K = ii + Hp)\)/ii-Hp)\). (6) Thus, the Beltrami coefficient fx gives us important information about the properties of the map (See argfi Fig. 3 Illustration of how the Beltrami coefficient /i measures the distortion of a quasi-conformal mapping that maps a SMALL circle TO AN ELLIPSE WITH DILATION K. Figure [3]). Now, suppose /x and a are the Beltrami coefficients of the quasiconformal maps /^ and /"^ respectively. Then the Beltrami coefficient r of the composition map f^ = f'^o (/^)"^ can be computed as: T ^-^1,0 1 — fia 9 {n-\ (7) where 6 = ^ and p = ^f^^{z). In particular, if f^ is the identity, that is, if a = 0, then 4. Main Algorithm (8) In this section, we discuss in detail the main algorithms in this paper. Our goal is to look for a simple representation scheme for the space of surface diffeomorphisms, with the least number of constraints possible, to simplify the optimization process. A. The Beltrami Holomorphic Flow In this part, we describe two theorems about the Beltrami Holomorphic Flow(BHF) on the sphere §^ and the unit disk D. All the algorithms developed in this paper are mainly based on these theorems. Theorem 4.1 (Beltrami holomorphic flow on §^J.* There is a one-to-one correspondence between the set of quasiconformal diffeomorphisms of S^ that fix the points 0, 1, and oc and the set of smooth complex-valued functions /i on §^ with \\ii\\^ = k < 1. Here, we have identified S^ with the extended complex plane C. Furthermore, the solution /^ to the Beltrami equation depends holomorphically on /x. Let {/x(t)} be a family of Beltrami coefficients depending on a real or complex parameter t. Suppose also that /x(t) can be written in the form IJi{t){z) = ii{z) + tv{z) + te{t){z) (9) for z G C, with suitable /i in the unit ball of ^^(C), v, e{t) G L^(C) such that || e{t) ||oo^ as t ^ 0. Then for all t(; G C, /mW(^) = ff^^yj) + tV{r, u){w) + o{\t\) (10) locally uniformly on C as t ^ 0, where ^ ^^ ^ ^ yc/^(^)(/^(^)-i)(/^(^)-/^H) Proof: This theorem is due to Bojarski. For detailed proof, please refer to [2]. ■ Theorem 4.1 states that any diffeomorphism of §^ that fixes 0, 1 and oc can be represented uniquely by a Beltrami coefficient. In fact, the 3-point correspondence can be arbitrarily set, instead of fixing 0, 1 and oc only. This can be done easily by composing Mobius transformations to the diffeomorphism. Let / : S^ ^ S^ be any diffeomorphism of §^. Picking any 3-point coresspondence {a, b,c ^ §^} ^ {/(a), /(6), /(c) G S^}, we can look for unique Mobius transformations (j)i and 02 that map {a, 6, c} and {/(a), /(6), /(c)} to 0, 1, oc respectively. Then, the composition map / := (j)2ofo(j)~^ is a diffeomorphism of §^ that fixes 0, 1 and oc and can be represented by a unique Beltrami coefficient. In other words, given a diffeomorphism / of §^ and any 3-point correspondence, we can represent / uniquely by a Beltrami coefficient. The theorem also gives the variation of the diffeomorphism under the variation of the Beltrami coeffi- cient. In order to adjust the diffeomorphism, we can simply adjust the Beltrami coefficient by using the variational formula. Theorem |4.1| can be further extended to diffeomorphisms of the unit disk D. Theorem 4,2 (Beltrami holomorphic flow on Dj.* There is a one-to-one correspondence between the set of quasiconformal diffeomorphisms of D that fix the points and 1 and the set of smooth complex- valued functions /x on D for which \\fi\\oo = k < 1. Furthermore, the solution /^ depends holomorphically on fi. Let {/x(t)} be a family of Beltrami coefficients depending on a real or complex parameter t. Suppose also that /x(t) can be written in the form //(t)(z) = /i(z) + tiy{z) + te{t){z) (12) for z G D, with suitable /x in the unit ball of ^^(D), z/, e{t) G L^(D) such that || e{t) ||oo^ as t ^ 0. Then for all ^ G D r^'\w) = r{w) + ty(/^ v){w) + o{\t\) (i3) locally uniformly on D as t ^ 0, where rH(rH-i) l/(/^z.)H \2 ^^^W')^^^)f ,.,,+ r ^ ^^^^_ (14) f^{z){f^^{z) - l){f^{z) - f^^iw)) Jo fi^(z){l - Mz)){l - f^{z)ft^{w)) Proof: The proof of this theorem can be found in the Appendix Theorem |4.2| states that any diffeomorphism of D that fixes 2 points (i.e. and 1) can be represented uniquely by a Beltrami coefficient. Again, the 2-point correspondence can be arbitrary. Let g: B) ^B) he a diffeomorphism of D. Given any 2-point correspondence {a, 6 G D} f^ {g{ci),g{h) G D}, we can find two unique Mobius tranformations ^i and 02 of D that map {a, h} and {g{o),g{h)} to {0, 1} respectively. Then, the composition map ^ := (p2 ^ 9 ^ 0r^ i^ ^ diffeomorphism of D that fixes and 1 and can be represented by a unique Beltrami coefficient. Theorem |4.2| also gives the variation of the diffeomorphism of D under the variation of the Beltrami coefficient. Therefore, we can again adjust the diffeomorphism of D by adjusting the Beltrami coefficient, which is a much simpler functional space. Theorem |4.1| and Theorem |4.2| can be extended to genus closed surfaces and open surfaces with disk topology. Therefore, they can be applied to represent general surface diffeomorphisms. This will be discussed in Section 14. 2[ B. Representation of Surface Diffeomorphisms using BCs As mentioned earlier, it is crucial to look for a simple representation for the space of all surface diffeomorphisms so that the optimization procedure can be simplified. Surface registration is commonly represented by 3D coordinate functions in M^. This representation requires lots of storage space and is difficult to manipulate. For example, 3D coordinate functions have to satisfy a constraint in the Jacobian J (namely, J > 0) in order to preserve the 1-1 correspondence of surface maps. The Jacobian constraint is a complicated partial differential inequality. Enforcing this constraint adds extra difficulty in manipulating and adjusting surface maps. It is therefore important to have a simpler representation with as few constraints as possible. Theorem 4J_ and 4^ allow us to represent surface diffeomorphisms of §^ and D by Beltrami coefficients. The theorems can be further extended to genus closed surfaces and open surfaces with disk topology. Let Si and 5^2 be two genus closed surfaces with a 3-point correspondence between them: {pi,P2,P3 ^ Si} ^ {?i, ?2, ?3 ^ 5^2}. By Riemann mapping theorem. Si and 5^2 can both be uniquely parameterized by conformal maps (j)i: Si ^ S^ and (p2' S2 ^ S^ respectively, such that (j)i{pi) = 0, (pi{p2) = I, 0i(p3) = oc and (p2{Qi) = 0,02(^2) = 1,02(^3) = 00. The conformal parameterizations can be computed using the discrete Ricci flow method [fTT1|. Given any surface diffeomorphism f : Si ^ S2. The composition map / '= 02o/o0r^ : §^ ^ S^ is a diffeomorphism from S^ to itself fixing 0, 1 and oc. By Theorem 4.1 , / can be uniquely represented by a Beltrami coefficient Jl defined on S^. Hence, / can be uniquely represented by a Beltrami coefficient /x := /I o (/)~^ defined on Si. In other words, we have proven the following: Corollary 4.3: Let Si and S2 be two genus closed surfaces. Suppose /: S'l ^ 6^2 is a surface diffeomorphism. Given 3-point correspondence {pi,P2,Ps ^ Si} ^ {f(Pi)i f(P2)i /(Ps) ^ 5^2}, / can be represented by a unique Beltrami coefficient fi : Si ^ C. Similarly, Theorem |4.2| can be extended to open surfaces with disk topology. Let Mi and M2 be two genus open surfaces. Given two points correspondence {pi,P2 ^ ^1} ^ {?i, ?2 ^ ^2} between them. We can again uniquely parameterize Mi and M2 conformally to map the corresponding points to and 1. Denote them by (f)i : Mi -^ D and 02 : M2 -^ D. The composition map / := 02 o / o 0r^ : D ^ D is a diffeomorphism of D fixing and 1. Again, / can be uniquely represented by a Beltrami coefficient Jl defined on D. Hence, / can be uniquely represented by a Beltrami coefficient ji \= ]1 o (p~^ defined on Ml. So, we have the following Corollary: Corollary 4.4: Let Mi and M2 be two genus open surfaces with disk topology. Suppose / : Mi ^ M2 is a surface diffeomorphism. Given 2-point correspondence {pi,P2 G Mi} ^ {f{Pi)^ f(P2) ^ M2}, f can be represented by a unique Beltrami coefficient /x: Mi ^ C. Corollary |4.3| and |4.4| allows us to represent diffeomorphisms of genus closed surfaces and open surfaces with disk topology using Beltrami coefficients. Thus, we can use the Beltrami coefficient jHf associated uniquely to such diffeomorphism / to represent /. First of all, we need to compute the Beltrami coefficient Jlj of the composition map / = 02 o / o 0~^ : D ^ D, where D is the common conformal parameter domain of the surfaces.. Mathematically, Jlj is given by the following formula: Then, the Beltrami coefficient /x/ can be computed by /x/ := /Iro0~^ : S'l ^ C. /i/ is a complex-valued functions defined on Si with ||/x/||oo < 1. There are no restrictions on jif that it has to be 1-1, surjective or satisfy some constraints on the Jacobian. With this representation, we can easily manipulate and adjust surface maps without worrying about destroying their diffeomorphic property. In practice, surfaces are commonly approximated by discrete meshes comprising of triangular or rectangular faces. They are parameterized onto the mesh D in C. Then the partial derivatives (or gradient) of the map can be discretely approximated on each face of D. By taking average, the partial derivatives and hence the Beltrami coefficient can be computed on each vertex. The detailed numerical implementation can be found in the Appendix. Besides adjusting surface maps preserving diffeomorphism, another advantage of Beltrami coefficients is that they consist of two real functions only, namely the real and imaginary parts. Compared to the representation of surface maps using 3D coordinate functions, this representation reduces 1/3 of the original storage space. The computational algorithm can be summarized as follows: Algorithm 1. Beltrami Representation of Surface Diffeomorphisms Input: Surface diffeomorphism f: Si ^ S2; points correspondence {pi} ^ {qi = /{Pi)}- Output: Beltrami representation jif'. Si ^ C of f : Si ^ S2' 1) Compute the conformal parameterizations of Si and S2 that map {pi} and {qi} to consistent locations on the parameter domain D. Denote them by (j)i'. Si ^ D and (j)2'- S2 ^ D 2) Set / = 02o/o0r^- D ^ D and compute the Beltrami coefficient Jlj by Equation 3) Compute the Beltrami coefficient /Hf: Si ^ C using jif := Jljo 0~^. 15 10 Brain 1 Iteration 10 Brain 2 Iteration 15 Beltrami Coefficient Iteration 20 Fig. 4 Beltrami representation and reconstruction of a surface diffeomorphism / on the brain surface. The top-left figure shows a surface diffeomorphism between two different brain surfaces. the top-right figure shows the Beltrami representation /h of /. Here, the colormap of |/i| is shown. The bottom row shows the reconstructed map AFTER DIFFERENT NUMBER OF ITERATIONS N USING BHF RECONSTRUCTION. WHEN N = 20, THE MAP CLOSELY RESEMBLES THE ORIGINAL MAP (THE BLACK DOTS SHOW THE EXACT POSITIONS UNDER THE ORIGINAL MAP.) C. Reconstruction of Surface Diffeomorphisms from BCs Given the Beltrami coefficient jjl defined on Si. It is important to have a reconstruction scheme to compute the associated quasi-conformal diffeomorphism /^. This allows us to move back and forth between BCs and surface diffeomorphisms. We propose the Beltrami holomorphic flow (BHF) method to reconstruct the surface diffeomorphism f^\ Si ^ S2 associated with a given /x. BHF iteratively flows the identity map to /^. In this part, we describe the BHF reconstruction method in detail. The variation of /^ under the variation of /x can be expressed explicitly. Suppose Jl{z) = fi{z)+tu{z) + o{\t\) where z = x + ly e C. Then, f^^^'^w) = f^{w) + tV{f^, u){w) + o(|t|), where V{f^,u){w) = / K{z,w)dxdy, Jd (16) and K(z, w) f^(w)(f^(w)-l) ( u(z){{f^^),{z))^ \f^'{z){f^'{z)-l){f^^{z)-f^{w)) f^{w){f^{w)-l) f u{z){{fnz{z))^ + Hz)((f^)4z)r J^{z){f^'{z)-l){f^{z)-f^{w)) ' ff-(z)(l-fi-(z))(l-ff-(z)ff-(w)) We can also write V{f^, ^){^) ^s: GiUi + 021^2 D = S\ D (17) 'D \ G3Z/1 + G4Z/2 (18) where u = ui + 1^2 and Gi, G2, G3, G4 are real- valued functions defined on D. Here, we identify A + iB 11 HIGH low Irtiiial map Reltr^mi CoelTifient 5 lteratjnit$ 10 Iterations 15 Iterations 20 IleriitLt>aii Fig. 5 Beltrami representation and reconstruction of a surface diffeomorphism / on hippocampal surfaces. The top-left figure shows a surface diffeomorphism between two different hippocampal surfaces. the top-right figure shows THE Beltrami representation /h of /. The colormap of |/i| is shown. The bottom row shows the reconstructed map AFTER DIFFERENT NUMBER OF ITERATIONS N USING BHF RECONSTRUCTION. WHEN N = 20, THE MAP CLOSELY RESEMBLES THE ORIGINAL MAP (THE BLACK DOTS SHOWS THE EXACT POSITIONS UNDER THE ORIGINAL MAP.) Using this fact, we propose to use BHF to iteratively flow the identity map to /^. Given the parame- terizations (/)i: Si ^ D and 02 : 5^2 ^ D, we look for the map /^ = 02 o /^ o (pi^ : D ^ D associated uniquely with /I = /x o (j)~^ : D ^ C f^ can then be obtained by /^ = (/)2^ o f^ o (j)^. We start with the identity map Id of which the Beltrami coefficient is identically equal to 0. Let A^ be the number of iterations. Define Jlk = kJl/N, fc = {0, 1, 2, . . . A^}. Let f^^ be the map associated with Jlk. Note that /^^ = Id and /^^ = /^. Equation 16 allows us to iteratively compute f^^ and thus obtain a sequence of maps flowing from Id to /^. The iterative scheme is given by: (19) The computational algorithm of the reconstruction scheme can be summarized in Algorithm 2. The detailed numerical implementation can be found in Appendix. Algorithm 2. Reconstruction of Surface Diffeomorphisms from BCs Input: Beltrami Coefficient ji on Si; conformal parameterizations of Si and S2: 0i and 02,* Number of iterations N 12 Output: Surface dijfeomorphism f^: Si ^ S2 associated to fi. 1) Set k = 0; 7^° = Id. 2) Set Jlk := kJl/N; Compute J^^+i = p^ + l/(7^^ ^); A; = A; + 1. 3) Repeat Step 2 until k = N; Set f^ := c/)^^ oJ^ocjy^-.Si^ S2. Figure |4] and [5] illustrate the idea of reconstructing surface diffeomorphisms from BCs on human brain surfaces and hippocampal surfaces respectively. BHF computes a sequence of surface maps {/^^} converging to /^. The approximation of /^^ is more accurate with a smaller time step, or equivalently, a larger number of iterations A^. Figure |6] shows the error of the reconstructed map /^^ versus different number of iterations A^ used in the BHF process. The error is defined as Error = sup ||/^^ — fW where / is the original map. As expected, the error decreases as A^ increases. In practice, the approximations are very accurate when A^ > 15. In our experiments, we set A^ = 20. o Error =svii,\\f''-f\\ ( f ^ original map; f^^ = reconstructed map ) -• — • — •— • — • — • — •- Number of iterations in Beltrami Holomorphic flow Fig. 6 The error of the reconstructed map f^^ versus the number of iterations used in the BHF PROCESS. D. BHF Optimization of Surface Registrations We have described a simple representation scheme for surface diffeomorphisms using BCs. The space of BCs is a simple functional space with the least amount of constraints. There are no restrictions requiring BCs to be 1-1, surjective or satisfy some constraints on its Jacobian. With BCs, we can easily manipulate and adjust surface maps, while ensuring the diffeomorphic(l-l, onto and smooth) property of the surface registration. Theorem 4.1 and 4.2 give us the variation of surface maps under the variation of their BCs (Equation T6| and [17]). This allows us to perform optimization on the space of BCs, instead of working directly on the space of surface diffeomorphisms. The diffeomorphic property of the optimal surface registration can also be easily ensured during the optimization process. 13 Given an energy functional E defined on the space of surface diffeomorphisms, we can easily refor- mulate E and redefine it on the space of BCs. With the BHF variation, we can derive the Euler-Lagrange equation on E to optimize BCs iteratively. To demonstrate the idea, we consider a simple example of optimizing surface maps between two human brain surfaces. (A) Brain 1 (B) Brain 2 I te ratio 11$ (C) Energy in each iteration Fig. 7 Illustration of BHF optimization scheme on brain surfaces. This example shows the optimization result of MATCHING TWO FEATURE FUNCTIONS Fi AND F2 ON THE TWO BRAIN SURFACES. THE BLUE GRID REPRESENTS THE INITIAL MAP, WHILE THE BLACK GRID REPRESENTS THE OPTIMIZED MAP. Example 4,1: Consider two different human brain surfaces 5'i and 5^2 as shown in Figure [7| Denote the conformal parameterizations of them hy (j)i\ Si ^ D and (j)2'. S2 ^ D, \n surface registration, it is often important to find an optimal 1-1 correspondence that matches some intensity feature functions defined on each surfaces. Let Fi\ Si ^ M and F2 : 6^2 ^ M be two intensities (functions) defined on 5^1 and S2 respectively. As an illustration, we define Fi and F2 as Fi := (/)^^(5.2x^ + 3.3?/^) and F2 := 0^^(6.8x2 + 2.82/). We propose to find f : Si ^ S2 minimizing E{f) = j sS^i{w) - F2{f {w))f + \li{w)\^ dw. The optimized map / is a quasi-conformal map that best matches Fi and F2 while preserving the conformality as good as possible. We can formulate the energy functional to be defined on the space of BCs over the conformal parameter domain D. That is. E{fi) = / (FiH - F^in? + H^Wdw (20) D The Euler-Lagrange equation can be derived as follow: = - j 2(Fi - F2(/^))VF2(r)^|,=o/'^+*'^ -2ii-vdw JdJd\B ) ' \ G3U1 + G^U2 ^_ [ ( f ( AG^ + BGs ^^ Jd \Jd V AG2 + BG^ (21) 14 where ( | = 2(Fi - F2(/^))VF2; /x = /xi + z/X2 and z/ = z/i + zz/2. So, the descent direction for /x = /xi + i/X2 is ^ = / (AGi + SG3) dw - 2/xi and ^ = / (AG2 + SG4) ^^ - 2/X2.; (22) (it J^ (it Jjj We can iteratively optimize the energy E as follow: .-^ = .^ + rft f ^-(^^^^ + ^^^^^ '" - '^^ ^ (23) Figure [7] shows the experimental result for this example. (A) shows the standard grid on Brain 1. The standard grid is mapped by the initial map to Brain 2, which is shown as the blue grid. We optimize the map such that it minimizes the energy functional. The resulting map is plotted as the black grid. (C) shows the energy at each iteration. It decreases as the number of iterations increases. This shows that our BHF optimization algorithm can iteratively optimize the energy functional. ■ Therefore, with BHF, we can perform optimizations over the space of BCs, which is a much simpler functional space with least amount of contraints, and simplify the optimization procedure significantly. 5. Applications In this section, we outline some applications of our proposed optimization algorithm to surface regis- tration. These applications are motivated from practical problems encountered in medical imaging. A. Optimized Conformal Parameterization with Landmark Matching With BHF, we first develop an algorithm to effectively compute landmark-matching optimized conformal maps between surfaces. A landmark-matching optimized conformal map refers to a map that matches corresponding landmarks across surfaces, while preserving conformality as much as possible. It is very important for research applications in computational anatomy. For example, in human brain mapping, neuroscientists are often interested in finding a 1-1 correspondence between brain surfaces that matches sulcal/gyral landmark curves, which are important anatomical features ifTTl . Besides matching these brain features, they also want the maps to preserve local geometry as much as possible. Conformal maps are best known to preserve local geometry and hence are commonly used. However, landmark matching cannot be guaranteed under conformal maps. Therefore, it is of interest to look for maps which are as conformal as possible and match landmarks well. Most existing algorithms for computing landmark-matching optimized conformal maps cannot ensure exact landmark matching. Some existing algorithms can align landmarks consistently, but bijectivity is usually not guaranteed especially when a large number of landmark constraints are imposed [18J. Here, we introduce a variational approach to compute an optimized conformal map iteratively by minimizing the L^ norm of a Beltrami coefficient /x (the Beltrami energy). Since /x is a measure of local distortion in conformality, our proposed algorithm is in fact looking for the best landmark-matching map, which is as 15 conformal as possible. Also, a map is bijective as long as ||/x||oo < 1- This can be easily controlled and guaranteed in each iteration by minimizing the Beltrami energy in our algorithm. c^i Conformvil map (bo Surface 2 [^ a D D Fig. 8 This figure shows the framework of the landmark-matching optimized conformal parameterization algorithm. Given two surfaces Si and 5^2 with the same topology. Denote the corresponding landmark curves on Si and S2 by {Cl} and {Cl} respectively. We first parameterize Si and S2 conformally onto a common parameter domain L) (= D or S^ = C). Let (j)i\ Si ^ D and 02 : 5^2 ^ i^ be the parameterizations. We proposed to look for two maps Lpi\ D ^ D and Lp2'. D ^ D such that Lp^^ ( i = 1, 2) maps landmarks {(t)i{Cl)} onto the consistent landmarks {Ck} on D (see Figure |8J), and that it minimizes the following energy functional: E{^i) \^^^^ (24) D Equation |24] ensures that each landmark-matching parametrization Lpi has the least conformality dis- tortion. Hence, the local geometric distortion under ipi is minimized. Starting from the conformal map with // = 0, the energy also ensures the property that |/x|oo < 1 and so the diffeomorphic property of the minimizer is guaranteed. A landmark-matching map / between Si and S2 can then be obtained by the composition map: / := 0^^ 0LP20 Lp^^ o 0^. We can compute the Euler-Lagrange equation of Equation 24 with respect to /i^. as follow: (25) 7 r /~i ^\t=oE{fi^^ +tv)= —\t = 0|/i^^ + tv\^ Re(/i^jRe(^)) + Im(/x^JIm(^))] -2fi^.. Hence, we can iteratively minimize E{fi^.) The derivative in Equation 25 is negative when v by the following scheme: M. n+l K -2/i",dt. (26) The detailed computational algorithm can be described as follow: 16 Iterations (D) Iterations (E) Fig. 9 Illustration of landmark-matching optimized conformal parameterizations of synthetic surfaces with 1 landmark. the blue curves in (a) and (b) represent the landmarks on the two surfaces. the landmark on surface a cannot be mapped onto the one on surface b under the conformal map(the black curve in (b)). with optimized conformal parameterization, the corresponding landmarks on each surface can be exactly matched, as shown in (C). The PERCENTAGE CHANGE IN ENERGY FUNCTIONALS OF THE OPTIMIZED CONFORMAL PARAMETERIZATIONS FOR SURFACE A AND B ARE SHOWN IN (D) AND (E). Algorithm 3. Optimized Conformal Parameterization with Landmark Matching Input: Surfaces Si and S2; Landmark curves Cl on Si, C^ on 5^2. Output: Optimized conformal parameterization ^pi and 9^2 of Si and S2 with landmark matching, 1) Compute the initial map 99^ that aligns landmark curves {0i(C^)} to {C/e} on D. Set n = 0. 2) Compute the Beltrami coefficient /x^. of Lp^. Let ji^^ = /x^. — 2/x^.rft. 3) Compute Vn = V{lp^^ — 2/i^.) using the BHF formula. 4) Let (f^^^{p) = (f^{p) + 5{p)Vn{p)dt , where 5 is a smooth delta function on D that is equal to zero around {Cl} and one elsewhere. This ensures landmarks are matched in each iteration. Set n = n+L 5) Repeat Step 2 to Step 5. If \E{iii^-l^) - ^(/x^JJI < e, stop. We have tested our proposed method on synthetic data as v^ell as real medical data. Figure [9] show^s the result of matching tw^o synthetic surfaces vv^ith one landmark on each surface. The blue curves on (A) and (B) represent the landmarks on the tw^o surfaces. Under a conformal map, the landmark on surface A cannot be mapped exactly onto the one on surface B (the black curve in (B)). Using our proposed 17 method, the corresponding landmarks on each surface can be exactly matched, as shown in (C). (D) and (E) show the percentage change in energy functional of the optimized conformal parameterizations for surface A and B. The energies decrease as the number of iterations increases. This indicates a decrease in the conformalilty distortion. Figure [TO] shows the Beltrami coefficient of each optimized conformal parameterization. The colormap shows the norm of the Beltrami coefficient. Note that the norm of the Beltrami coefficient is very small except near the landmark curve. It means the conformality distortion is accumulated around the landmarks as expected. (A) (B) Fig. 10 The Beltrami coefficients of the optimized conformal parameterizations of 2 synthetic surfaces fixing 1 landmark. the norms of the beltrami coefficients are plotted as colormap, which are very small except near the landmark curves. it means that the conformality distortion is accumulated around the landmarks. We have also tested our algorithm on synthetic surfaces with five landmarks as shown in Figure 11 Again, the landmarks cannot be exactly matched under a conformal map (see black curves in (B)). However, they are exactly matched using our proposed algorithm. As shown in (D) and (E), the per- centage change in energies decreases as the number of iterations increases, meaning that conformality distortion is progressively reduced. Figure [T2| shows the Beltrami coefficients of the optimized conformal parameterizations fixing landmarks. Again, the norm of the Beltrami coefficient is very small except near the landmark curves. Finally, we have tested our algorithm on real cortical hemispheric surfaces extracted from brain MRI scans, acquired from normal subjects at 1.5 T(on a GE Signa scanner). Figures [13]^ A) and (B) show 2 different brain surfaces with 3 major sulcal curves labeled on each of them (see the blue curves). Under a conformal map, landmarks on Brain 1 and Brain 2 are not exactly matched (see the black curves in (B)). They are, however, exactly matched using our proposed algorithm as shown in (C). (D) and (E) show the percentage change in the energies of the optimized conformal parameterizations of the surfaces. The energies decrease as the number of iterations increases. This shows that the conformality distortion is gradually reduced. Figure [T4| shows the Beltrami coefficients of the optimized conformal parameterizations of the 2 brain surfaces. Again, the norms of the Beltrami coefficients are very small except near the sulci curves. 18 Iterations (D) a 15C 200 Iterations (E) Fig. 11 Landmark-matching optimized conformal parameterizations of 2 synthetic surfaces fixing 5 landmarks. The blue curves on (a) and (b) represent the landmarks on the surfaces. the landmark on surface a cannot be mapped to the landmark on surface b under the conformal map(black curves in (b)). with optimized conformal parameterization, the corresponding landmarks on each surface can be exactly matched (shown in (c)). the percentage change in energies of the optimized conformal parameterizations for surface a and surface b are plotted in (d) and (e). B. Hippocampal Registration with Geometric Matching In medical imaging, there are cases w^here anatomical landmark features cannot be easily defined on some brain structures. In such cases, landmark-matching constraints cannot be used as a criterion to establish good correspondence betvv^een surfaces. Finding the best 1-1 correspondence betw^een these structures becomes challenging. One typical example is the hippocampus(HP), w^hich is an important structure in the human brains. It belongs to the limbic system and plays important roles in long-term memory and spatial navigation. Surface-based shape analysis is commonly used to study local changes of HP surfaces due to pathologies such as Alzheimer disease (AD), schizophrenia and epilepsy [|25ll . On HP surfaces, there are no w^ell-defined anatomical landmark features. High-field structural or functional imaging, vv^here discrete cellular fields are evident Il30ll . is still not routinely used. Finding meaningful registrations betv^een HP surfaces becomes challenging. It is thus important to develop methods to look for good registrations betw^een different HP surfaces vv^ithout landmarks. To achieve this, w^e develop an algorithm to automatically register HP surfaces v^ith complete geometric matching and avoid the need to 19 Fig. 12 The Beltrami coefficient of the optimized conformal parameterization of 2 synthetic surfaces fixing 5 landmarks. The norm of the Beltrami coefficient is very small except near the landmark curve. It means that the conformality distortion is accumulated around the landmarks. manually label landmark features. This is done by optimizing a compounded energy, v^hich minimizes the L^ norm of the Beltrami coefficient and matches curvatures defined on each surface. Given tv^o hippocampal surfaces Si and 5^2. The compounded energy E shape is defined mathematically as E. shape {^A a f \^\^+p f(H,- H^inf +7/(^1-^2 Jd Jd Jd inr (27) w^here Hi, H2 are the mean curvatures on Si, S2 respectively, defined on the common parameter domain D, and Ki, K2 are the Gaussian curvatures. The first integral minimizes the conformality distortion of the surface registration, and ensures the diffeomorphic property of the minimizer by controlling /x. The second and third integrals ensure the optimized registration matches the curvatures as much as possible. It turns out that E'shape is a complete shape index vv^hich measures the dissimilarity betvv^een tw^o surfaces. Specifically, E'shape = if and only if Si and 5^2 are geometrically equal up to a rigid motion. Therefore, surface map minimizing E'shape is the best registration that matches the geometric information as much as possible. We can minimize E'shape in Equation 27 iteratively, using the proposed BHF optimization algorithm. The Euler-Lagrange equation of Equation [27] can be computed as foUow^s: — |f=0^shape(M) + 7 d L dt «=o H2ir^n) (Ki - K^ir^n) ^2a ^i-v-2j3 {Hi JD JD dt (28) 2 j {afi{w) - JiiPH + jK)-l^'], i/3H + jK) ]}-v{w) 20 (D) Iterations (E) Fig. 13 Landmark-matching optimized conformal parameterizations of cortical hemispheric surfaces with 3 major sulcal landmarks. the blue curves on (a) and (b) represent the landmarks on the two surfaces. under a conformal map, the landmarks on brain a cannot be correctly mapped onto landmarks on brain b (black curves in (b)). with landmark-matching optimized conformal parameterization, the corresponding landmarks on each surface can be exactly matched as shown in (c). the percentage change in energies of the optimized conformal parameterizations FOR Brain A and Brain B are shown in (D) and (E). vv^here J^ • := J^ • dw and Sz^ '^ Sd^ ^^ ^^^ defined as the integral over variables w and z respectively, H 18 [Hi - H2{f^))VH2{f^). K := {Ki - X2(/^))VX2(/^), G^ is as defined in Equation -2(/x(^) - j^[{H + K)-G, det{H + K, G)] ). The derivative in Equation 28 is negative vv^hen v Hence, wt can iteratively minimize E{fi) by the foUov^ing iterative scheme: /x^ = -2(a/x^ - f[{f3H'' + 7ic^) • G^, det{f3H'' + 7^^^, G^)] ) dt II n+l (29) The detailed computational algorithm can be described as foUovy^s: Algorithm 4. BHF Registration with Geometric Matching Input: Hippocampal surfaces Si and S2, step length dt, threshold e Output: Geometric matching registration f^ and the shape index E{f^) 1) Compute the conformal parameterizations of Si and 5^2. Denote them by (j)i'. Si ^ D and 02 2) Set ip^ := Id: D ^ D and n = 0. 21 Brain 1 - -0.; - -O.' - -0.- [ I Brain 2 Fig. 14 The Beltrami coefficient of the optimized conformal parameterizations of 2 cortical hemispheric surfaces with 3 major sulcal landmarks. the norms of the beltrami coefficients are shown as colormap, which are very small except near the landmark curves. as expected, the conformality distortion is accumulated around the landmarks. 3) Compute the Beltrami coefficient jjJ^ of Lp^ (e.g. jjl = 0). Update jjJ^'^ by Equation 4) Compute: Vn = V{(f^^ fi^'^^ — jjJ\) using Equation 16 Let (p^^^ = 99^ + 14- Set n = n-\-l. 29 5) Repeat Step 3 to Step 5. If \E{fi^-^^) - E{fi^)\ < e, Stop. We have tested our algorithm on 212 HP surfaces automatically extracted from 3D brain MRI scans v^ith a validated algorithm [23 J. Scans w^ere acquired from normal and diseased (AD) elderly subjects at 1.5 T (on a GE Signa scanner). In our experiments, w^e set a = 1 and /3 = 7 = 2. Experimental results show that our proposed algorithm is effective in registering HP surfaces with geometric matching. The proposed shape energy can also be used to measure local shape difference between HPs. Figure [TSj^A) shows two different HP surfaces. They are registered using our proposed BHF algorithm with geometric matching. The registration is visualized using a grid map and a texture map, which shows a smooth 1-1 correspondence. The optimal shape index E'shape is plotted as a colormap in (B). E'shape effectively captures the local shape difference between the surfaces. (C) shows the shape energy in each iteration. With the BHF algorithm, the shape energy decreases as the number of iterations increases. (D) shows the curvature mismatch energy {E = //3(i^i — i^2(/))^ + 7(^1 — ^2(/))^)- It decreases as the number of iterations increases, meaning that the geometric matching improves. (E) shows the Beltrami coefficient of the map in each iteration, which shows the conformality distortion of the map. Some conformality is intentionally lost to allow better geometric matching. Figure 16 shows the BHF registration between two normal HPs. The complete shape index E'shape is plotted as colormap on the right. Again, E'shape can accurately capture local shape differences between the normal HP surfaces. Figure [T7] shows the BHF hippocampal registrations between normal elderly subjects and subjects with Alzheimer's disease. The BHF registrations give smooth 1-1 correspondence between the HP surfaces. We can use the complete shape index E'shape to detect local shape differences between healthy and unhealthy 22 HIGH LOW m \ Energy' = -HMnf + -<iJ<i -lUif) B.29 0.23 Energ>' = -Ji'^tr.).)" «U n.265 ^ 026 ^ t2i \ 0,2SS D2fi 023 ^ 120 WO 160 180 Iteration (C) 30 100 120 no 160 139 200 Iteration (D) Fig. 15 Shape registration with geometric matching using Beltrami Holomorphic Flow (BHF). The registration is visualized as grid map and texture map as shown in (a). the optimal shape energy is shown in (b). the percentage changes of the shape energy, curvature mismatch energy and beltrami energy after different number of iteration are shown in (d), (e) and (f) respectively. subjects. We also study the temporal shape changes of normal and AD HP surfaces, as shov^n in Figure [TSJ We compute the deformation pattern of its HP surfaces for each subject, measured at time = and time =12 months (see [l22l for longitudinal scanning details). The left tw^o panels show^ the temporal deformation patterns for tw^o normal subjects. The middle tw^o panels shovv^ the temporal deformation patterns for tv^o AD subjects. The last column shovv^s the statistical significance p-map measuring the difference in the deformation pattern betw^een the normal (n=47) and AD (n=53) groups, plotted on a control HP. The deep red color highlights regions of significant statistical difference. This method can potentially be used to study factors that influence brain changes in AD. 6. Conclusion In this paper, v^e propose a simple representation of bijective surface maps using Beltrami coeffi- cients(BCs), w^hich helps the optimization process of surface registrations. To complete the representation scheme, wt develop a reconstruction algorithm of the surface diffeomorphism from a given EC using the Beltrami holomorphic flow method. This allows us to move back and forth between BCs and surface diffeomorphisms. By formulating the variation of the associated surface map under the variation of BC, we reformulate variational problems over the space of surface diffeomorphisms into variational problems over 23 HIGH LOW Normal 1 Normal 2 Shape index Fig. 16 BHF REGISTRATION BETWEEN TWO NORMAL SUBJECTS. THE SHAPE INDEX ^shape IS PLOTTED ON THE RIGHT, WHICH CAPTURES LOCAL SHAPE DIFFERENCES. Normal Alzheimer 1 Shape index Normal Alzheimer 2 Shape index Fig. 17 BHF REGISTRATION BETWEEN 2 NORMAL SUBJECTS AND 2 SUBJECTS WITH ALZHEIMERS DISEASE. THE LOCAL SHAPE DIFFERENCES CAPTURED BY ^shape ARE PLOTTED ON THE SURFACES. the space of BCs. It greatly simplifies the optimization procedure. More importantly, a bijective surface map is always guaranteed during the optimization process. Experimental results on synthetic examples and real medical applications show the effectiveness of our proposed algorithms for surface registration. 7. Appendix I. Numerical Implementation In this part, we give detailed numerical implementation on how the proposed algorithms can be com- puted. In practice, all surfaces are represented by meshes which consist of vertices, edges, and triangu- lar/rectangular faces. In our iterative scheme, the functions and their partial derivatives are defined on each vertex and then linearly interpolated to define their values inside each triangular/rectangular face. 7. Computation of the Beltrami Coefficient 24 Normal 1 Normal 2 Alzheimer 1 Alzheimer 2 p^U.Ol p-valuc Shape energy for deformation pattern (from t=0 to t= 12 Months) Stiitistk'iil sijfnificance p-niap Fig. 18 Temporal HC shape changes of normal and subjects with Alzheimers disease. Let / = (/i, /2) be the diffeomorphism defined on the parameter domain D. The Beltrami coefficient /x/ associated uniquely to / can be computed as follows (see Equation 15): r.dfi a/2 df2 , dfi dfi a/2 df2 dfi ^^ = [^ a^ - a^^ + '^^ + a^^]/[( a^ + a^^ + '^^ " a^^^' (30) In order to compute /i/, we simply need to approximate the partial derivatives at each vertex: Dxfi{v) ^{v) and Dyfi{v) ^ |^(^). We first approximate the gradient Vrfi on each face T by solving: Vi -vo V2 -Vo ^Th ft{vi)-fi{vo) liTi-iTol fi{v2)-fi{vo) |'i^2-'i^o| (31) where [iTq, iTi] and [w^, wl] are two edges on T. After the gradient Vr/i have been computed for each face T, Dxfi{v) and Dyfi{v) can be computed by taking average as follows: (32) DJi{v) where A'^ is the set of all faces around the vertex v. Hence, the Beltrami coefficient fJ^fiv) can be computed by: .^^ (DJliv) - Dyf2{v)) + i{Dj2{v) + Dyfljv)) ,.,. 5] VT/i/|iVd (D./i(i;) + DyMv)) + liDMv) - DyMv)) 2. Computation of the BHF Reconstruction For the BHF reconstruction algorithm, the most important step is the computation of the variation V{f^, v) of /^ under the variation of /x. We will discuss the computation of V{f^, v) for D = D. The computation for Z) = §^ = C is similar. From Equation 16 and 17 V{f^,v){w) = / K{z,w)dxdy D 25 where r{w){r{w)-i) K{z^ w) = 71 u{z){{rUz)f ^ u{z){{f^Uz)f f>^{z){f'^{z) - l){f'^{z) - //^H) /A.(^)(l - /A'(^))(l - fl^{z)f^^{w)) Now, f^ and v are both defined on each vertex v. Also, {f^)z{v) can be approximated as: For each pair of vertices (iT, t(;), K{y,w) can be easily approximated. In case K{y,w) is singular, we set K{v, w) = 0. Now, for each vertex v, we define A^ as A^= Y^ Area{T)/NT, (35) where A^ is the number of vertices on T. That is, Nt = 3 if T is a triangle and A^^ = 4 if T is a rectangle. Then, V{f^,u) can be approximated by: V{r,u){w) = Y,Ki^^^)Av (36) V II. Proof of Theorem 1X2) : To prove the theorem, we need the following lemma: Lemma 7.1: Let /: D — )■ D be a diffeomorphism of the unit disk fixing and 1 and satisfying the Beltrami equation /^ = fxfz with fj, defined on D. Let / be the extension of / to C defined as (f{z), if|^|<l, fiz)=r\'' '- (37) Then / satisfies the Beltrami equation h = i^fz (38) on C, where the Beltrami coefficient /i is defined as . u(z), if \z\ < 1, /i(^) = { ,_l^ ' ' - (39) j/x(l/z), if \z\ > 1. Proof: We need to prove / satisfies the Beltrami equation: h = A/. (40) 26 Clearly, / satisfies equation (38) on D. Outside D, we consider / and / as functions in z and z. Note that: d_ dz fiz,z) d_ dz f{z,z) We have: df{z, z) d 1 -2 d dz dzf{i/z,l/z) '^ ' ' ' ' dz'^ ' ' ' ' -2 d -f{l/z, 11 z) ^/(l/^, lA) = -/(l/^, lA) (-l/^^)/.(lM 11 z) Also, = z-^j{\l-zMz)- j^{\l-zMz). df{z,z) d 1 -2 d ^ -f{l/z,l/z) —f{l/z,l/z) dz dzf(i/z,l/z) ■>^ ' ^ ' ' dz'^ ' ' ' ' -2 d -f{l/z, 1/z) Q-J{l/z, 1/z) = -f{l/z, 1/z) {-l/z')Ml/z, 1/z) z-'fil/z,l/z) Ml/z,l/z) = z-'f{l/z,l/z) fiil/z)Ml/z,l/z) Now from Equation |42 Thus, we have. Ml/z,l/z) = z'fil/z,l/z) ■2df{z^ dz df{z,z) dz z-'fil/z,l/z) ^iil/z)Ml/z,l/z) = z-'f{l/z, 1/z) fi{l/z)z'f{l/z, 1/z) z^ , .. df{z,z) .df{z,z) = =2/«(lA) ^^ = Kz)- ■^d f{z,z) dz dz dz (41) (42) (43) (44) (45) a) Proof of Theorem \4.2\ According to Quasiconformal TeichmuUer Theory, there is a one-to-one correspondence between the set of quasiconformal homeomorphisms of C fixing 3 points and the set of smooth complex- valued functions /x on D for which sup |/x| = /c < 1. If a diffeomorphism / on C satisfies equations (38)(39), then l//(l/z) also satisfies the same equation. By the uniqueness of the solution according to Theorem 4.1, we must have f{z) = l//(l/z). On 5D, z = 1/z. This implies f{z) = l//(z), and hence \f{z)\ = 1 on 5D. Therefore, by restricting the solution of equation (38) on C fixing 0, 1 and oc to D, we can get a diffeomorphism of D fixing and 1. Equation \\\) can then be applied to get the variational formula V{f^, v) of /^ under the variation v of /x. To get V{f^, u), we 27 evaluate the integral in equation ( [TT] ). Jc V HzKHz) - i)iMz) - Mw)) dxdy ^''^^^f'^^^'^^' d.dyU ^^^^mmt^^.d.dy (46) /P(^)(/M(^) _ l){f,^z) - f,(yj)) ' ^ y^VD HZ){MZ) - 1){HZ) - Mw)) Now, outside the disk D, K^) = |JRiM and ^ = z-'f{l/z, 1/z)-'m1/z, 1/z) (47) We have: -mrwy d.dy Hz){Mz) - i){nz) - Mw)) y ^ icw Hz){f^{z) - i){Mz) - Hw)) y HzKMz) - i)iMz) - Mw)) ' ^ yc\D Mz){Hz) - i){.Hz) - Hw)) f,(z){Mz)-l){Mz^Mw)) y f^iii^y {f^{ii^)~ -i){f^{i/^)~ -Mw)) z\ ^(^W')^(^)y ,.,,+ /^ ,^,^ Jj, f''iz){f'^{z) - l){ft^{z) - fi^iw)) Lf^'{z){l-f^^{z)){l-f^^{z)f^^{w)) Substituting Equation 22 into Equation [TTJ we get an integral flow equation on D, which is given by f^{z){f^{z) - l){Hz) - f^{w)) h Hz){l - f^{z)){l - Hz)f^{w)) References [1] S. Angenent, S. Haker, A. Tannenbaum, and R. Kikinis. On the Laplace-Beltrami operator and brain surface flattening. IEEE Transaction of Medical Imaging, 18(8):700-711, 1999. [2] S. Durrleman, X. Pennec, A. Trouve, R Thompson, and N. Ay ache. Measuring brain variabiUty via sulcal Hnes registration: A diffeomorphic approach. Medical Image Computing and Computer-Assisted Intervention (MICCAI 2007) Lecture Notes in Comput. Sci. 4791, Springer-Verlag, Berlin, Heidelberg, page 675682, 2007. [3] S. Durrleman, X. Pennec, A. Trouve, P. Thompson, and N. Ayache. Inferring brain variabihty from diffeomorphic deformations of currents: An integrative approach. Medical Image Analysis, 12:626637, 2008. [4] B. Fischl, M. I. Sereno, R. B. Tootell, and A. M. Dale. High-resolution intersubject averaging and a coordinate system for the cortical surface. Human Brain Mapping, 8:272-284, 1999. [5] F. Gardiner. Quasiconformal Teichmilller Theory. American Mathematics Society, 2000. [6] J. Glaunes, M. Vaillant, and M. I. Miller. Landmark matching via large deformation diffeomorphisms on the sphere. Journal of Mathematical Imaging and Vision, 20:179-200, 2004. 28 [7] X. Gu, Y. Wang, T. Chan, P. Thompson, and S.-T. Yau. Genus zero surface conformal mapping and its appUcation to brain surface mapping. IEEE Transactions on Medical Imaging, 23(8):949-958, 2004. [8] S. Haker, S. Angenent, A. Tannenbaum, R. Kikinis, G. Sapiro, and M. Halle. Conformal surface parameterization for texture mapping. IEEE Transaction of Vision and Computer Graphics, 6(8): 181189, 2000. [9] M. Hurdal and K. Stephenson. Cortical cartography using the discrete conformal approach of circle packings. Neurolmage, 23:S1 19S128, 2004. [10] M. K. Hurdal and K. Stephenson. Discrete conformal methods for cortical brain flattening. Neuroimage, 45:86-98, 2009. [11] M. Jin, J. Kim, F. Luo, and X. Gu. Discrete surface Ricci flow. IEEE Transactions on Visualization and Computer Graphics, 14(5): 1030-1043, 2008. [12] S. Joshi and M. Miller. Landmark matching via large deformation diffeomorphisms. IEEE Transactions on Image Processing, 9(8): 13571370, 2000. [13] L. Ju, J. Stern, K. Rehm, K. Schaper, M. Hurdal, and D. Rottenberg. Cortical surface flattening using least square conformal mapping with minimal metric distortion. IEEE International Symposium on Biomedical Imaging, pages 11-SO, 2004. [14] O. Lehto and K. Virtanen. Quasiconformal Conformal Mappings in the Plane. Springer- Verlag New York, 1973. [15] A. Leow, C. Yu, S. Lee, S. Huang, H. Protas, R. Nicolson, K. Hayashi, A. Toga, and P. Thompson. Brain structural mapping using a novel hybrid implicit/explicit framework based on the level-set method. Neurolmage, 24(3):910-927, 2005. [16] N. A. Lord, J. Ho, B. Vemuri, and S. Eisenschenk. Simultaneous registration and parcellation of bilateral hippocampal surface pairs for local asymmetry quantification. IEEE Transactions on Medical Imaging, 26(4):471-478, 2007. [17] L. Lui, Y. Wang, T. Chan, and P. Thompson. Brain anatomical feature detection by solving partial differential equations on general manifolds. Discrete and Continuous Dynamical Systems B, 7(3):605-618, 2007. [18] L. Lui, Y Wang, T. Chan, and P. Thompson. Landmark constrained genus zero surface conformal mapping and its application to brain mapping research. Applied Numerical Mathematics, 5:847-858, 2007. [19] L. M. Lui, S. Thiruvenkadam, Y Wang, T. Chan, and P. Thompson. Optimized conformal parameterization of cortical surfaces using shape based matching of landmark curves. SIAM Journal of Imaging Science, 3(l):52-78, 2010. [20] B. Lvy, S. Petitjean, N. Ray, and J. Maillot. Least squares conformal maps for automatic texture atlas generation. ACM Transactions on Graphics (TOG), 21(3):362 - 371, 2002. [21] O. Lyttelton, M. Bouchera, S. Robbinsa, and A. Evans. An unbiased iterative group registration template for cortical surface analysis. Neurolmage, 34:1535-1544, 2007. [22] J. Morra and et al. Automated mapping of hippocampal atrophy in 1-year repeat mri data in 490 subjects with alzheimer's disease, mild cognitive impairment, and elderly controls. Neuroimage, 45(1):S3-15, 2009. [23] J. H. Morra, Z. Tu, L. G. Apostolova, A. E. Green, C. Avedissian, S. K. Madsen, N. Parikshak, X. Hua, A. W Toga, J. CUfford R. Jack, M. W Weiner, and P. M. Thompson. Validation of a fully automated 3d hippocampal segmentation method using subjects with alzheimer's disease, mild cognitive impairment, and elderly controls. Neuroimage, 43(l):59-68, 2008. [24] Y Shi, P. M. Thompson, I. Dinov, S. Osher, and A. W Toga. Direct cortical mapping via solving partial differential equations on implicit surfaces. Medical Image Analysis, 11 (3): 207-223, 2007. [25] P. M. Thompson, K. M. Hayashi, G. L de Zubicaray, A. L. Janke, S. E. Rose, J. Semple, M. S. Hong, D. H. Herman, D. Gravano, D. M. Doddrell, and A. W Toga. Mapping hippocampal and ventricular change in alzheimer's disease. Neurolmage, 22(4): 1754-66, 2004. [26] D. Tosun, M. Rettmann, and J. Prince. Mapping techniques for aligning sulci across multiple brains. Medical Image Analysis, 8:295309, 2004. [27] Y Wang, M.-C. Chiang, and P. M. Thompson. Automated surface matching using mutual information applied to Riemann surface structures. Proceeding in Medical Image Computing and Computer- Assisted Intervention MICCAI 2005, 2:666-61 A, 2005. [28] Y Wang, L. Lui, X. Gu, K. Hayashi, T. Chan, A. Toga, P. Thompson, and S. Yau. Brain surface conformal parameterization using riemann surface structure. IEEE Transactions on Medical Imaging, 26(6): 853-865, 2007. [29] Y Wang, L. M. Lui, T. F. Chan, and P. M. Thompson. Optimization of brain conformal mapping with landmarks. Proceeding in Medical Image Computing and Computer-Assisted Intervention MICCAI 2005, pages 675-683, 2005. [30] M. M. Zeineh, S. A. Engel, P. M. Thompson, and S. Y Bookheimer. Dynamics of the hippocampus during encoding and retrieval of face-name pairs. Neurolmage, 299(5606):577-580, 2003.