PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Comput Graph. Author manuscript; available in PMC 2010 June 1.
Published in final edited form as:
Comput Graph. 2009 June 1; 33(3): 299–311.
doi:  10.1016/j.cag.2009.03.002
PMCID: PMC2802331
NIHMSID: NIHMS102478

Fourier method for large scale surface modeling and registration

Abstract

Spherical harmonic (SPHARM) description is a powerful Fourier shape modeling method for processing arbitrarily shaped but simply connected 3D objects. As a highly promising method, SPHARM has been widely used in several domains including medical imaging. However, its primary use has been focused on modeling small or moderately-sized surfaces that are relatively smooth, due to challenges related to its applicability, robustness and scalability. This paper presents an enhanced SPHARM framework that addresses these issues and show that the use of SPHARM can expand into broader areas. In particular, we present a simple and efficient Fourier expansion method on the sphere that enables large scale modeling, and propose a new SPHARM registration method that aims to preserve the important homological properties between 3D models. Although SPHARM is a global descriptor, our experimental results show that the proposed SPHARM framework can accurately describe complicated graphics models and highly convoluted 3D surfaces and the proposed registration method allows for effective alignment and registration of these 3D models for further processing or analysis. These methods greatly enable the potential of applying SPHARM to broader areas such as computer graphics, medical imaging, CAD/CAM, bioinformatics, and other related geometric modeling and processing fields.

Keywords: Spherical harmonics, spherical parameterization, SPHARM expansion, surface registration, spherical thin plate spline

1. Introduction

Spherical harmonics were first used as a type of parametric surface representation for radial or stellar surfaces r(θ, ϕ) by Ballard et al. [1, 2]. An extended method, called SPHARM, was proposed by Brechbühler et al. [3] to model more general shapes, where three functions of θ and ϕ were used to represent a surface. These spherical harmonic methods have recently received a lot attention, and have been studied and applied to applications in various fields including computer vision [1, 3], graphics [4, 5, 6, 7, 8], medical image analysis [9, 10, 11, 12], bioinformatics [13, 14, 15], and biology [16, 17, 18].

SPHARM is a highly promising Fourier method for modeling arbitrarily shaped but simply connected 3D objects, where protrusions and intrusions can be effectively handled. Thanks to its underlying spherical parameterization, SPHARM is suitable for many surface manipulation and analysis applications, including texture mapping, morphing, remeshing, compression, statistical modeling, and surface-based morphometry, because of the following reasons: (1) processing is much easier on the spherical domain than on an irregular mesh; (2) the spherical domain is continuous, allowing for more flexible processing on it than on a regular mesh; and (3) processing can be done not only in the spatial domain but also in the frequency domain.

A typical SPHARM processing pipeline includes three key steps: (1) spherical parameterization [3, 6, 7, 8, 19], (2) SPHARM expansion [3, 20], and (3) SPHARM registration [3, 21]. Spherical parameterization has been extensively studied and we will briefly review a few existing methods that are appropriate for different SPHARM modeling conditions in Section 2. However, the existing methods still have limitations in the SPHARM expansion and registration steps that prevent SPHARM from being applied to its full potential for large scale and flexible 3D modeling and analysis.

In this paper, we introduce a new SPHARM framework to overcome these limitations. In particular, we emphasize that a recently proposed, simple and efficient Fourier expansion method on the sphere is an essential component of our framework for modeling large scale 3D surfaces. In addition, we propose a new SPHARM registration method that aims to keep the important homological relationships between 3D models. Although SPHARM is a global descriptor, our experimental results show that the proposed SPHARM framework can accurately describe complicated graphics models and highly convoluted surfaces and the proposed registration method allows for effective alignment of these 3D models for further processing or analysis.

The rest of the paper is organized as follows. Section 2 briefly describes the SPHARM framework. Section 3 describes our SPHARM expansion method that can deal with large scale 3D models. Section 4 proposes a new SPHARM registration method that aims to align surface landmarks as well as minimize area and length distortions. Section 5 demonstrates our experimental results. Section 6 concludes the paper.

2. SPHARM description

The SPHARM method was proposed by Brechbüuhler et al. [3] to model arbitrarily shaped but simply connected 3D objects. It is essentially a Fourier transform technique that defines a 3D surface using three spherical functions and transforms them into three sets of Fourier coefficients in the frequency domain. Three steps are often involved in a typical SPHARM processing pipeline: (1) spherical parameterization, (2) SPHARM expansion, and (3) SPHARM registration.

Step 1: Spherical parameterization creates a continuous and uniform mapping from the object surface to the surface of a unit sphere, and its result is a bijective mapping between each point v on a surface and a pair of spherical coordinates θ and ϕ:

v(θ,ϕ)=(x(θ,ϕ),y(θ,ϕ),z(θ,ϕ))T.

In Fig. 1, Panel (a) shows a hippocampal surface extracted from a magnetic resonance imaging (MRI) scan and Panel (f) shows its spherical parameterization. This parameterization is an area preserving mapping computed using Brechbüuhler’s method [3].

Figure 1
Sample SPHARM reconstruction: (a) A hippocampal surface, (b-e) its SPHARM reconstructions using coefficients up to degrees 1, 5, 10 and 15, respectively, (f) its spherical parameterization, (g) a regular mesh grid on the sphere that is used for reconstructions ...

Step 2: SPHARM expansion expands the object surface into a complete set of spherical harmonic basis functions Ylm, where Ylm denotes the spherical harmonic of degree l and order m and it is essentially a Fourier basis function defined on the sphere. The expansion takes the form:

v(θ,ϕ)=l=0m=llclmYlm(θ,ϕ),
(1)

where clm=(cxlm,cylm,czlm)T. The Fourier coefficients clm up to a user-desired degree can be estimated by solving a linear system. The object surface can be reconstructed using these coefficients, and using more coefficients leads to a more detailed reconstruction. A sample reconstruction case is shown in Fig. 1: The original model is shown in (a) and its SPHARM reconstructions using coefficients up to degrees 1, 5, 10 and 15 are shown in (b-e), respectively. Different spherical sampling schemes can be used to reconstruct the object. Reconstructions shown in (b) and (c) are created using a regular mesh grid on the sphere (g), and reconstructions shown in (d) and (e) are created using a level 3 icosahedral subdivision (h).

Step 3: SPHARM registration creates a shape descriptor (i.e., excluding translation, rotation, and/or scaling) from a normalized set of SPHARM coefficients, which are comparable across objects. The first two steps are necessary for modeling an individual shape, while the third step is optional and mainly required for pair-wise or group analysis of 3D models such as morphing or surface-based morphometry.

Although SPHARM has been successfully applied to a few medical imaging problems, numerous fundamental challenges remain before its potential can be fully utilized in general cases. In this paper, we address two major challenges: one for SPHARM expansion and the other for SPHARM registration. Section 3 and Section 4 discuss these two steps in more details, respectively. For completeness, here we briefly review the existing methods for spherical parameterization.

In spherical parameterization, the traditional method [3] aims to create an equal area mapping as well as minimize angle distortions by solving a constrained optimization problem. This is very effective in analyzing small and moderately-sized structures extracted from volumetric images (e.g., MRI, CT). However, its scalability is limited and, in addition, it is applicable only to voxel surfaces. CALD [7] is a newer method that extends the traditional method and can be applied to general triangular meshes. Conformal mapping has also been used for spherical parameterization in some SPHARM studies [6]. It has solid mathematical foundation but tends to introduce large area distortion that may not be ideal for establishing surface correspondence between models in statistical shape analysis. Two recent studies [8, 19] use progressive meshes and stretch metrics to minimize vector length distortion and their parameterization results look very promising for handling large scale graphics models. Shown in Fig. 2 is a sample result of [19]. Since spherical parameterization has been extensively studied, people can often find a method that fits their needs.

Figure 2
Sample spherical parameterization.

3. Large scale SPHARM expansion

Given a 3D model v(θ, ϕ) and a user-specified maximum degree Lmax, the task of SPHARM expansion is to extract coefficients clm in Eq. (1) for lLmax and |m | ≤ l. There are two types of approaches for computing clm: one uses numerical integration [14, 20]; the other formulates a linear system and solves it using least square fitting (LSF) [3].

The naive numerical integration method is inefficient and so not applicable to large models. Healy et al. [20] proposed a fast algorithm to accelerate the integration procedure using a divide-and-conquer strategy; and several studies [4, 6, 8] employed this method. However, to use Healy’s method, a preprocessing step is required to remesh the model using a regular spherical mesh grid (e.g., Fig. 1(g)). This step is inconvenient and also it tends to introduce additional remeshing errors.

In contrast, the LSF method [3] is very easy to implement and works directly on original object meshes even if they are irregular. The LSF method seems to be a more popular method in SPHARM studies, especially in brain imaging. The bottleneck of this method is its limited scalability, because most methods for solving large linear systems [22, 23] are either designed for sparse or symmetric matrices or not easy to implement. This might be the reason that SPHARM has been mostly used for modeling small or moderately sized 3D surfaces.

We describe an iterative residual fitting (IRF) method that overcomes this limitation and facilitates the opportunity of creating large scale SPHARM models. The basic idea behind the IRF method is simple and follows the properties of spherical harmonic transform. First, these harmonics form a coarse-to-fine hierarchy. If we just use a few low degree harmonics to expand a spherical function f (θ, ϕ), we get a low-pass filtered reconstruction. If we use more degrees, more details are included in the reconstruction. The IRF method takes advantage of this coarse-to-fine hierarchy. It starts from a low degree reconstruction and then iteratively adds more details into our model by involving higher degree harmonics. Second, spherical harmonics form an orthonormal basis and geometric information is stored in different frequency channels. Thus, if we first extract information from low frequency channels, the residual (i.e., f (θ, ϕ)–its reconstruction) will exactly contain information in high frequency channels. To add in more details to our model, we can simply use a few higher degree harmonics to fit the residual.

The IRF method breaks a large linear system into several small linear systems and thus SPHARM modeling can be easily done at a large scale on standard workstations with average configuration using a standard linear solver. We also believe that much larger SPHARM models can be created if one combines the IRF method with an enhanced large-scale linear solver. IRF was originally presented in [24, 9], where this method was tested on a few models with around 40, 000 vertices and it worked well. In this paper, we aim to demonstrate that the IRF method can be used to accurately model larger scale graphics models and to enable a broader range of surface processing applications. We also want to emphasize that IRF is an essential component in the proposed enhanced SPHARM framework for modeling large scale 3D surfaces.

Note that a degree l SPHARM model involves (l + 1)2 × 3 complex coefficients. Our experiments show that a degree 85 SPHARM model described by 22,188 complex coefficients can reasonably capture surface details of an original model with 100,002 vertices and 200,000 faces that is described by 900,006 variables in total. Additional quantifiable information about the reconstruction errors and the compression ratios is available in Table 1 of Section 5.1. For more details about IRF (e.g., projection errors, optimal frequency band selection), please refer to [9, 24].

Table 1
Effects of SPHARM degrees on the average error Ea and the compression ratio (CR). CR is defined as a ratio of (N × 3) versus M, where N is the total number of vertices (involving N ×3 variables) on the original model and M = (degree + ...

4. Landmark guided SPHARM registration

SPHARM registration is an important operation allowing for pairwise processing or group analysis across different SPHARM models and is a critical step in many applications (e.g., surface-based morphometry in medical imaging [9, 10, 11, 12, 15], evolutionary morphology in biology [16, 17, 18]). The traditional method [3] uses the first order ellipsoid for alignment, and works only if this ellipsoid is a real ellipsoid and may not work well in many other cases. SHREC [21] is a newer SPHARM registration method that minimizes the mean square distance between corresponding surface parts and works for general cases. Both methods are designed for relatively smooth surfaces without landmarks. However, landmarks often contain important prior knowledge about the objects. For example, they tend to have critical anatomical or biological meanings in medical and biological applications and should not be ignored in the registration procedure.

In this paper, we present a new SPHARM registration method that is guided by a set of pre-existing landmarks. Since SPHARM is a parametric model, registration should be done in both object and parameter spaces. In the object space, we first scale all the objects with landmarks to have a normalized size, and then apply the iterative closest points (ICP) method [25] to align landmarks together in a least square fashion. Now the key problem is how to align landmarks in the parameter space.

Given two SPHARM models, a template and an individual, we can distort the parameterization of the individual to match its landmarks with the corresponding landmarks on the parameterization of the template. This can be done by applying spherical thin plate spline (STPS) [26]. We observe that there are many different ways to apply STPS and the results are very different from one another. Note that the correspondence between SPHARM models is implied by the underlying parameterization: two points with the same parameter pair (θ, ϕ) on two surfaces are defined to be a corresponding pair. Thus, in order to create an ideal correspondence, our goal is to identify an underlying parameterization that is the least distorted. To achieve this goal, we present a few strategies that aim to find the “best” parameterization from these STPS results.

Zou et al. [27, 28] did a similar study, where they directly applied STPS to spherical conformal parameterizations but did not consider to reduce area and length distortions. In the field of the SPHARM surface modeling, we feel that an equal area mapping is more attractive than other mappings because we want to treat each area unit on object surface equally by assigning the same amount of parameter space to it (i.e., similar to arc-length parameterization for comparing 2D contours). Our method is designed to achieve this goal. Asirvatham et al. [29] did another similar study, where landmark constraints were used to guide the spherical parameterization procedure in a progressive mesh framework. In our case, we don’t want to re-parameterize the entire object, since the existing parameterization has already been optimized based on certain criteria. Our goal is simply to distort as little as possible the existing parameterization in order to match landmarks.

In the rest of this section, we describe our landmark-guided SPHARM registration method. We assume that objects are already aligned to one another in the object space and their SPHARM descriptions are known (pre-calculated). We first briefly describe STPS, then define the area and length distortions we try to minimize, and finally presents our methods.

4.1. Spherical thin plate spline

Spherical thin plate spline (STPS), defined in [26], is an extension of 2D thin plate spline to a spherical domain and has been used to deform spherical parametric domains [27, 28]. STPS on a spherical domain S2 minimizes a bending energy J(u), subject to u(Pi) = zi, i = 1, 2, … , n, where Pi [set membership] S2 and zi is the fixed displacement at Pi. This bending energy is formulated as

J(u)=02π0π(Δu(θ,ϕ))2sinϕdθdϕ,
(2)

where θ [set membership] [0, π] is latitude, ϕ [set membership] [0, 2π ] is longitude, and Δ is the Laplace-Beltrami operator.

The solution on the sphere is given by

un(P)=i=1nciK(P,Pi)+d
(3)

where c and d are determined by

cn+d=,c=0,

n is the n × n matrix with i, jth entry K(Pi, Pj), = (1, …, 1)T, and = (z1, …, zn)T. K(X, Y) between two arbitrary points X, Y [set membership] S2 is defined as follows:

K(X,Y)=14π01(logh)(11h)(112ha+h21)dh

where a = cos(γ(X, Y)) and γ(X, Y) is the angle between X and Y. However, this is not in a computable closed form expression. Therefore, thin plate pseudo-spline on the sphere R(X, Y), defined in [26], is used in this study, which is formulated as

R(X,Y)=12π[1(2m2)q2m2(a)1(2m1)]
(4)

and

qm(a)=01(1h)m(12ha+h2)1/2dh,m=0,1,.

In our case, we have m = 2 giving

q2(a)=12{In(1+1W)[12W24W]12W3/2+6W+1},
(5)

where W = (1 – a)/2. Therefore, the Eq. 3 becomes

un(P)=i=1nciR(P,Pi)+d,
(6)

and c and d are determined by

Rnc+d=,c=0,

where Rn is the n × n matrix with i, jth entry R(Pi, Pj), = (1, …, 1)T, and = (z1, …, zn)T.

Given n pairs of two corresponding points {Pi, Pinew}, i = 1, … , n on the sphere, the displacements (Δθi, Δϕi) at a set of points {Pi, i = 1, … , n} are calculated. With this displacements, Eqs. 4, 5, and 6 are used to compute displacement (Δθk, Δϕk) of any point Pk [set membership] S2 and the point Pk, located at p(θ, ϕ), is moved to a new point Pknew at p(θ + Δθk, ϕ + Δϕk). For convenience, we call this the STPS-based displacement scheme.

4.2. Mesh distortion measures

In this study, STPS is employed to match the landmark positions of an individual object in the spherical parametric domain to the corresponding landmark positions of the template object and subsequently transform the underlying parametric mesh of the individual object. This algorithm distorts the parametric mesh of the individual object and can introduce additional errors to some extent to the reconstructed shape of the individual object. Therefore, mesh distortion cost functions are employed and these distortion costs are used as selection criteria for finding best rotation angles in the proposed algorithm.

Distortion of the parametric mesh by STPS is measured by calculating the area distortion cost (ADC) and the length distortion cost (LDC). The overall and worst costs for the whole parametric mesh are measured to evaluate the performance of the proposed algorithms in Section 4.3.

4.2.1. Area distortion measures

The concept of area distortion cost introduced by [7] is employed as one of the performance measures in this study. Let M = {ti} be a triangle mesh in the parametric space and let ψ be a continuous STPS-based displacement scheme, which maps M to a distorted parametric mesh ψ(M) = {ψ(ti)}. A(·) is used to denote the area of a triangle or a mesh. The area distortion cost (ADC) Ca with respect to ψ is defined as follows:

For each triangle ti [set membership] M,

Ca(ti,ψ)=A(ψ(ti))A(ti).

This measures the local ADC of a single triangle.

For each mesh vertex v in M,

Ca(v,ψ)=tiMvA(ψ(ti))tiMvA(ti),

where Mv is the set of triangle incident upon v. This measures the local ADC around a single vertex.

For the whole parametric mesh M,

Ca(M,ψ)=tiMmax(Ca(ti,ψ),1Ca(ti,ψ))A(ψ(ti))A(ψ(M)).
(7)

This measures the overall ADC for the whole mesh. By taking

max(Ca(ti,ψ),1Ca(ti,ψ))

as the ADC contribution from each triangle, we treat contraction and expansion equally, and so always have Ca(M, ψ) ≥ 1.

The worst ADC is defined as follows

CaW(M,ψ)=max{max(Ca(ti,ψ),1Ca(ti,ψ))|tiM}.
(8)

4.2.2. Length distortion measures

To measure length distortion introduced by STPS, the stretch concept by Sander et al. [30] is adopted in this study. They considered the case of mapping from a planar domain to 3D surface and at any point in the planar domain, two singular values of the 3 ×2 Jacobian matrix were computed to represent the largest and smallest length distortions when a vector in the 2D domain was mapped to the 3D surface. In our case, the length distortion cost (LDC) Cs with respect to a given mesh mapping ψ from M to ψ(M), is defined as follows:

Given a mesh mapping ψ from M to ψ(M),

Cs(M,ψ)=tiM(Γ(ti)2+1γ(ti)2)A(ψ(ti))2A(ψ(M)),and
(9)
CsW(M,ψ)=max{max(Γ(ti),1γ(ti))|tiM},
(10)

where Γ(ti) and γ(ti) are the largest and smallest length distortions for a triangle ti and A(·) is used to denote the area of a triangle or a mesh.

In the above definitions, Cs(M, ψ) measures the average length distortion cost (LDC) for the whole mesh M and the worst LDC is defined by CsW (M, ψ). The largest and smallest length distortions are directly computed from the length of three corresponding sides between ti and ψ(ti). Again, contraction and expansion are equally treated in both definitions.

4.2.3. Calculation of minimum distortion cost

At each step of the proposed approach, after STPS algorithm is applied, we calculate the average and worst ADCs as well as the average and worst LDCs, defined by Eqs. 7, 8, 9, and 10, respectively. In current experiments, we try to minimize the worst ADC. If several STPS results have the same worst ADCs, then their average ADCs, the worst LDCs, and the average LDCs are compared in order. This rule is used to avoid the extreme contraction or expansion while maintaining the reasonable area distorting and length distortion. We plan to explore other criteria using these measures in future studies.

Algorithm 1 Basis STPS method for SPHARM registration
1:Map landmarks of the template and an individual onto a common parametric mesh, e.g. an icosahedral subdivision
2:Distort the individual parametric mesh using STPS to move its landmarks to the location of the template’s landmarks
3:Calculate new SPHARM expansion of the individual using the STPS-distorted parameterization

4.3. Alignment of parametric mesh using STPS

4.3.1. Basic STPS algorithm

The basic method is to directly apply the STPS algorithm to two sets of landmarks with known correspondence between the sets, defined in a parametric mesh. Alg. 1 describes this approach to the landmark-guided alignment. Zou et al. [27, 28] used this approach in their studies. However, this naive approach can severely distort the parametric mesh of an individual object and result in a distorted reconstruction without re-sampling the parametric mesh, especially when one or more landmarks are located near the north or south pole (see Fig. 3).

Figure 3
Application of basic STPS method. When one or more landmarks are closely located to the north or south pole, the parametric mesh of the individual can be severely distorted. The average cost a and the worst cost w are shown as (a, w) for ADC and LDC.

4.3.2. Hierarchical STPS algorithm

Algorithm 2 Hierarchical STPS method
1:i:=1; assume that an parameter net is given
2:repeat
3: Create icosahedral samples at level i for α’s and β’s
4: if i=1 then
5:  Rotate the parameter net using each (αβ0), apply STPS, and keep the top K candidates that minimize the distortion costs, defined in Section 4.2.3
6: else
7:  Keep only local icosahedral samples for α’s & β’s
8:  Rotate the parameter net of each top K candidate using each (αβ0), apply STPS, and select top K candidates
9: i:=i+1
10:until The best distortion costs do not improve
11:Return the best result in top K list

To avoid large distortion introduced by STPS, we employ a sampling-based strategy that rotates the individual’s landmarks on the sphere using Euler angles (αβγ) in order to find the best oriented landmarks for applying STPS. The rotation space can be sampled nearly uniformly using icosahedral subdivisions (see Fig. 1(h) for an icosahedral subdivision mesh at level 3). This assigns rotation angles to α and β. We always set γ = 0, since the rotation along the z axis does not affect the distortion level of an STPS result (See Fig. 4 for an example).

Figure 4
Effect of rotation along z-axis on the results of STPS. Rotation along z-axis does not affect the distortion level of the STPS results. The average cost a and the worst cost w are shown as (a, w) for ADC and LDC.

To reduce computation time, we employ a hierarchical sampling scheme proposed in [21] to sample the rotation space instead of using all the icosahedral samples. We tested both the hierarchical method and the simple icosahedral method and the former greatly outperformed the latter (see Fig. 5). For each sampling point, the hierarchical approach moves it to the north pole, and the entire parametric mesh is rotated accordingly. Then, STPS is applied to each rotated parametric mesh and the best K rotation angles are selected for minimizing the distortion costs, defined in Section 4.2. This process is repeated for a higher level of sampling mesh until a certain criterion is satisfied. Alg. 2 shows more details of this approach.

Figure 5
Comparison of the number of samples in rotation space (left) and running time (right) for hierarchical and non-hierarchical rotational approaches.

4.3.3. Landmark rotation on template

While running the hierarchical STPS approach, we observed that template’s landmarks, if located near the poles, could also result in large distortion of the individual’s parametric mesh. This observation suggests that both template’s and individual’s landmarks should be adjusted to be as far away from the poles as possible to avoid large STPS distortions. A simple approach is to rotate not only individual’s but also template’s parametric mesh and test STPS on all possible pairs. This is obviously a very time consuming procedure. To accelerate it, we develop the following method to adjust the orientation of the landmarks.

Let P1, P2, …, Pn be all the landmarks on a sphere. Given any point P on the sphere, we define the following objective function f(P):

f(P)=maxi=1n|P·Pi|

Our goal is to find a P so that the above objective function is minimized and then we rotate the parameter net so that P becomes the north pole. Let θi be the angle between P and Pi. The shorter distance between Pi and one of the poles can be measured by |cos(θ)| = |P·Pi |. The smaller |cos(θ)| is, the bigger the distance is. Thus, if we minimize the above objective function, we can adjust the landmarks to the best orientation where the distance between a pole and its nearest landmark is maximized.

Therefore, the proposed two methods in Alg. 1 and Alg. 2 are modified by rotating template’s landmarks and individual’s landmarks away from the poles using the above approach before these methods are applied to the parametric mesh. After the application of the proposed methods, the distorted mesh is rotated back by using the inverse of a rotation matrix which is applied to the template’s landmarks before the STPS application step. Thus, the individual result matches the original template instead of a rotated template.

For convenience, we use BSTPS and R-BSTPS to denote Alg. 1 and its modified version, respectively. Similarly, we use HSTPS and R-HSTPS to denote the Alg. 2 and its modification, respectively. Shown in Fig. 6 is sample performance comparison of applying these four methods to one case. In general, HSTPS and R-BSTPS outperform BSTPS, and R-HSTPS performs the best in terms of the mesh distortion cost.

Figure 6
Performance comparison of different STPS methods. The first row shows (a) the original template Torg and (b) the original individual Iorg, and (c-d) the results of BSTPS and HSTPS. The second row shows (e) the rotated template Trot and (f) the rotated ...

5. Experimental results

We demonstrate our methods using four highly convoluted cortical models and a few publicly available large scale graphics models. The graphics models and their spherical parameterizations [19] were downloaded from Hugues Hoppe’s website at Microsoft Research. The cortical models and their initial spherical parameterizations were generated by applying the FreeSurfer software (http://surfer.nmr.mgh.harvard.edu/) on real MRI scans. The vertex numbers of these models range from 35K to 140K, and accordingly the face numbers range from 70K to 280K. The homologous landmarks between two models used for SPHARM registration were manually defined.

In this section, we first show that SPHARM models can be used to accurately describe these large-scale 3D models, which allows for many useful applications. Then, we show that the proposed SPHARM registration method can effectively align two models while preserving their homologous landmark relationships as well as minimizing parametric distortions. The experiments were performed on normal laptop and desktop computers running Matlab 7.2. With pre-computed spherical harmonic bases, it took 5 to 35 minutes to compute a SPHARM expansion for any of these models. The performance has a potential to be improved if we implement these methods using a lower level language (e.g., C) than Matlab.

5.1. SPHARM representation of large-scale 3D models

As a Fourier method, a SPHARM representation is a global descriptor. Even if this is the case, the proposed IRF method enables large scale SPHARM expansion that can reasonably capture fine surface details for complicated 3D surfaces and accurately describe them. Following typical literature in computer graphics [31, 32], we define the maximum reconstruction error Em and the average reconstruction error Ea as follows: Em=max1kN{pkrk2},Ea=1N1kNpkrk22, where pk is the k-th vertex (out of N) on the original model and rk is the corresponding vertex on a reconstructed surface. Similar to [31], in our experiments, all the models are scaled to have a bounding box of [0, 1]3, which makes the reconstruction error measure scaling invariant. Fig. 7 shows four original models and their SPHARM reconstructions using coefficients up to degrees 10, 20 and 85, respectively. The Ea measures are also shown for each case and also summarized in Table 1. We can see that a degree 85 SPHARM reconstruction is a pretty accurate representation of the original model visually and numerically. Shown in Table 2 is the comparison of reconstruction errors on two similar models (with different number of vertices) between this work and Cheng et al. [31]. Cheng et al. aimed to fit surfaces to unorganized point data, while our goal was to fit a SPHARM model to a given surface. Despite different tasks, the description errors were achieved at a similar level.

Figure 7
Shown from left to right are original models and their SPHARM reconstructions (with reconstruction error Ea) of degrees 10, 20 and 85. SPHARM can accurately model large scale 3D surfaces as well as be used for interesting applications such as geometric ...
Table 2
Performance comparison in terms of maximum (Em) and average (Ea) errors between our method (using Degree 85 models) and Cheng et al. [31].

With the goal of fitting a parametric model to a discrete representation of a triangulated surface, SPHARM inherently provides a mechanism for surface interpolation. While some surface interpolation methods (e.g., [33]) can deal with surfaces of arbitrary topology, SPHARM is only applicable to genus zero surfaces. Since SPHARM is a mathematical model defined on a continuous spherical domain, it enables easier processing for many applications involving arbitrarily-shaped but simply-connected 3D objects.

First, it is a more compact representation than a triangulated surface in many cases. For example, the original gargoyle model has 100,002 vertices and 200,000 faces, and so it is described by 900,006 variables in total. However, a degree 85 SPHARM model of gargoyle is described only by 22,188 complex coefficients and can reasonably capture the overall shape together with many surface details of the original model. This compactness property can be used for geometric compression. Table 1 shows effects of SPHARM degrees on the description error Ea and the compression ratio (CR). Even though the CR estimation is very conservative (i.e., ignoring all the face information in the original model), decent CRs and Ea’s can still be achieved for degree 85 reconstructions of the tested models.

Second, one can operate not only in the spatial domain but also in the frequency domain. Taking a lower order SPHARM reconstruction can naturally achieve the goal of surface smoothing and filtering. See Fig. 7 for a few samples. This property has been used in a couple of prior studies [4, 8] for surface smoothing and filtering. It can also be used for level of details representation and transferring.

Third, the level of details applications can be done not only via the frequency domain but also via the spatial domain. A SPHARM reconstruction is essentially a remeshed original model. We can use different spherical sampling schemes with different sampling resolutions for SPHARM reconstruction. Fig. 8 shows several reconstruction cases using regular spherical mesh grids (Fig. 1(g)) and icosahedral subdivisions (Fig. 1(h)) at different sampling resolutions. Although not tested in our experiments, the adaptive sampling mesh scheme described in [8] seems like a promising method that can derive more accurate reconstructions with fewer vertices.

Figure 8
SPHARM reconstructions using different spherical sampling schemes: The top two rows use icosahedral subdivisions and the bottom two rows use regular spherical meshes. Potential applications include remeshing, multi-resolution modeling, and level of details. ...

5.2. SPHARM registration

Table 3 presents the performance comparison among four proposed approaches to apply STPS to five cases. In general, R-BSTPS and R-HSTPS outperform BSTPS and HSTPS respectively in terms of the mesh distortion costs, and R-HSTPS exhibits the best performance. It is noticed that, however, HSTPS shows the better performance than R-HSTPS in the first case. However, even in this case, the distorted mesh using R-HSTPS does not introduce any noticeable distortion to its reconstructed object. This observation suggests that the criteria for selecting the best orientation of the landmarks on the sphere require further investigation for both the template’s parameterization and the individual’s parameterization.

Table 3
Comparison of mesh distortion costs among the proposed methods: each entry (a, w) contains the average cost a and the worst cost b.

Fig. 9 and Fig. 10 show the R-HSTPS alignment results of two pairs of cortical models and three pairs of graphics models respectively and their distortion costs can be found in Table 3. Each two rows correspond to one alignment case, where the top row shows the objects and the bottom row shows their spherical parameterization. We use bump maps to visualize the correspondence between the object and its parameterization. Landmarks and a coarse mesh grid are also shown on each surface. The first column shows template models. The second and the third columns exhibit individual models before and after SPHARM registration respectively. Comparing an original individual with the corresponding template, you will notice that their landmarks are aligned only in the object space but not in the spherical parameter space. However, comparing a registered individual with the corresponding template, you can see that their landmarks are aligned not only in the object space but also in the parameter space.

Figure 9
Two alignment results between cortical models. The first two rows show the result of aligning an individual left hemisphere to its template, while the last two rows show the results for a right hemisphere case. The first column shows a template model. ...
Figure 10
Three alignment results between graphics models. The first two rows show the result of aligning armadillo to gargoyle, the second two rows show the result of aligning gargoyle to bunny, and the last two rows show the result of aligning cow to horse. The ...

This improvement of alignment can subsequently contribute to the feasibility of performing pair-wise processing or group analysis of these 3D models. For example, such a SPHARM registered result can help improve a morphing sequence between two 3D models. Results shown in Fig. 11 exhibit several morphing sequences between a template and an original individual (Rows 1, 3, and 5) or between a template and a registered individual (Rows 2, 4, and 6). All the templates and the original and registered individuals are reconstructed using their degree 85 SPHARM models. We can see that morphing the template to the registered individual achieves a much better effect than to the original individual, since the intermediate shape of morphing to the registered individual looks more natural than those of morphing to the original individual. Besides morphing, there are many other SPHARM applications requiring models being registered with landmark guidance, such as surface-based morphometry in biomedical imaging [9, 10, 11, 12, 15] and morphological analysis in evolutionary biology [16, 17, 18].

Figure 11
Comparing morphing using original and registered SPHARM models. The first, third, and fifth rows show the results of morphing templates to original individuals and the second, forth, and sixth rows show the results of morphing templates to registered ...

6. Conclusions

We have presented an enhanced spherical harmonic (SPHARM) surface modeling and processing framework, and demonstrated that it could be used to accurately model large scale 3D surfaces and it could effectively register these models with or without landmark constraints. The main contribution is twofold: (1) incorporation of our recently proposed large scale SPHARM expansion method into the framework for facilitating the possibility of using SPHARM to accurately model large scale 3D surfaces, (2) design of a new STPS-based SPHARM alignment method that can register SPHARM models together under pre-existing landmark constraints as well as minimize parametric distortions. Our experimental results show the effectiveness of the proposed methods. These methods greatly enable the potential of applying the highly promising SPHARM method to broader areas such as computer graphics, medical imaging, CAD/CAM, bioinformatics, and other related geometric modeling and processing fields. The proposed methods will be incorporated into SPHARM-MAT, a 3D shape modeling and analysis toolkit, and will be released at the Neuroimaging Informatics Tools and Resources Clearinghouse (NITRC) website http://www.nitrc.org/ in the near future.

While image registration is considered as one of the most important topic in image processing and analysis, surface registration plays a similarly important role in computer graphics and shape analysis. The reasonable choice for criterion used for quantifying registration quality is critical. Different applications may require different criteria for obtaining optimal results. For smooth surfaces without landmarks (e.g., hippocampus, ventricle), we can use traditional SPHARM registration methods [3, 21] to align them. In these cases, their underlying parameterizations are not distorted. However, for surfaces associated with homologous landmarks (e.g., many biological structures), the underlying parameterizations have to be distorted so that the corresponding landmarks can be registered together across different objects. While the proposed registration methods are designed for this purpose and aim to control the area and length distortions, it remains an interesting future topic to examine the effects of using alternative criteria (i.e., considering a different combination of the area, length and angle distortions) for developing other effective registration schemes.

Acknowledgments

This work was supported in part by NIBIB/NEI R03 EB008674-01, NIA R01 AG19771, NCI R01 CA101318 and U54 EB005149 from the NIH, Foundation for the NIH, and grant #87884 from the Indiana Economic Development Corporation (IEDC). Graphics models were downloaded from Hugues Hoppe’s website for [19] at Microsoft Research. The cortical models were generated by the FreeSurfer software http://surfer.nmr.mgh.harvard.edu/. We thank the anonymous reviewers for their insightful comments and suggestions that help improve the paper.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

1. Ballard DH, Brown CM. Computer Vision. Prentice-Hall; Englewood Cliffs, N.J.: 1982.
2. Schudy R, Ballard D. Towards an anatomical model of heart motion as seen in 4-D cardiac ultrasound data. 6th Conf on Comp App in Rad & Anal of Rad Im. 1979
3. Brechbühler C, Gerig G, Kubler O. Parametrization of closed surfaces for 3D shape description. Computer Vision and Image Understanding. 1995;61(2):154–170.
4. Bulow T. Spherical diffusion for 3D surface smoothing. IEEE Trans on Pattern Analysis and Machine Intelligence. 2004;26(12):1650–1654. [PubMed]
5. Funkhouser T, Min P, et al. A search engine for 3d models. ACM Transactions on Graphics. 2003;22(1):83–105.
6. Gu X, Wang Y, Chan TF, Thompson PM, Yau S-T. Genus zero surface conformal mapping and its application to brain surface mapping. IEEE Transaction on Medical Imaging. 2004;23(8):949–958. [PubMed]
7. Shen L, Makedon F. Spherical mapping for processing of 3d closed surfaces. Image and Vision Computing. 2006;24:743–61.
8. Zhou K, Bao H, Shi J. 3D surface filtering using spherical harmonics. CAD. 2004;36(4):363–375.
9. Chung MK, Dalton KM, Shen L, Robbins SM, Evans AC, Davidson RJ. Weighted fourier series representation and its application to quantifying amount of gray matter. IEEE Transactions on Medical Imaging. 2007;26(4):566–581. [PubMed]
10. Gerig G, Styner M. MICCAI’2001. Vol. 2208. LNCS; 2001. Shape versus size: Improved understanding of the morphology of brain structures; pp. 24–32.
11. Gerig G, Styner M, et al. Shape analysis of brain ventricles using SPHARM. IEEE MMBIA. 2001:171–178.
12. Shen L, Ford J, Makedon F, Saykin A. A surface-based approach for classification of 3D neuroanatomic structures. Intelligent Data Analysis. 2004;8(6):519–542.
13. Cai W, Shao X, Maigret B. Protein-ligand recognition using spherical harmonic molecular surfaces: towards a fast and efficient filter for large virtual throughput screening. J Mol Graph Model. 2002;20(4):313–328. [PubMed]
14. Ritchie DW, Kemp GJ. Fast computation, rotation, and comparison of low resolution spherical harmonic molecular surfaces. Journal of Comp Chem. 1999;20:383–395.
15. Shen L, Saykin AJ, Chung MK, Huang H. Morphometric analysis of hippocampal shape in mild cognitive impairment: An imaging genetics study. IEEE 7th International Symposium on Bioinformatics & Bioengineering. 2007:211–217.
16. McPeek MA, Shen L, Torrey JZ, Farid H. The tempo and mode of 3-dimensional morphological evolution in male reproductive structures. American Naturalist. 2008;171(5):E158–E178. [PubMed]
17. McPeek MA, Shen L, Farid H. The correlated evolution of 3-dimensional reproductive structures between male and female damselflies. Evolution. 2009;63(1):73–83. [PubMed]
18. Shen L, Farid H, McPeek MA. Modeling 3-dimensional morphological structures using spherical harmonics. Evolution. in press. [PMC free article] [PubMed]
19. Praun E, Hoppe H. Spherical parametrization and remeshing. ACM Transactions on Graphics. 2003;22(3):340–349.
20. Healy D, Rockmore D, Kostelec P, Moore S. Ffts for the 2-sphere - improvements and variations. The Journal of Fourier Analysis and Applications. 2003;9(4):341–385.
21. Shen L, Saykin AJ, Huang H, Makedon F. Efficient registration of 3d spharm surfaces. Fourth Canadian Conference on Computer and Robot Vision. 2007:81–88.
22. Barrett R, Berry M, et al. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. 2. SIAM; Philadelphia, PA: 1994.
23. Freund R, Golub G, Nachtigal N. Iterative solution of linear systems. In: Iserles A, editor. Acta Numerica 1992. Cambridge University Press; 1992. pp. 57–100.
24. Shen L, Chung MK. IEEE 3DPVT’06. UNC; Chapel Hill: 2006. Large-scale modeling of parametric surfaces using spherical harmonics.
25. Besl PJ, McKay ND. A method for registration of 3-D shapes. IEEE Trans on PAMI. 1992;14(2):239–256.
26. Wahba G. Spline interpolation and smoothing on the sphere. SIAM Journal of Scientific and Statistical Computing. 1981;2:5–16.
27. Zou G, Hua J, Dong M. Integrative information visualization of multimodality neuroimaging data. 15th Pacific Conference on Computer Graphics and Applications; 2007. pp. 473–476.
28. Zou G, Hua J, Muzik O. Non-rigid surface registration using spherical thin-plate splines. MICCAI International Conference on Medical Image Computing and Computer-Assisted Intervention; 2007. pp. 367–374. [PubMed]
29. Asirvatham A, Praun E, Hoppe H. Consistent spherical parameterization. Computer Graphics and Geometric Modeling (CGGM) 2005 Workshop; 2005.
30. Sander PV, Snyder J, Gortler SJ, Hoppe H. Texture mapping progressive mesh; SIGGRAPH 2001, Computer Graphics Proceedings; 2001. pp. 409–416.
31. Cheng K, Wang W, Qin H, Wong K, Yang H, Liu Y. Fitting subdivision surfaces to unorganized point data using sdm. Proc of the 12th Pacific Conf on Computer Graphics and Applications; 2004. pp. 16–24.
32. Marinov M, Kobbelt L. Optimization methods for scattered data approximation with subdivision surface. Graphical Models. 2005;67(5):452–73.
33. Zheng J, Cai YY. Interpolation over arbitrary topology meshes using a two-phase subdivision scheme. IEEE Transactions on Visualization and Computer Graphics. 2006;12(3):301–10. [PubMed]