Skip to main content

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.