(Guest Editors)
Intrinsic Parameterizations of Surface Meshes
Mathieu Desbrun USC
Mark Meyer Caltech
Pierre Alliez USC / INRIA
Abstract
Parameterization of discrete surfaces is a fundamental and widely-used operation in graphics, required, for in- stance, for texture mapping or remeshing. As 3D data becomes more and more detailed, there is an increased need for fast and robust techniques to automatically compute least-distorted parameterizations of large meshes. In this paper, we present new theoretical and practical results on the parameterization of triangulated surface patches.
Given a few desirable properties such as rotation and translation invariance, we show that the only admissible parameterizations form a two-dimensional set and each parameterization in this set can be computed using a sim- ple, sparse, linear system. Since these parameterizations minimize the distortion of different intrinsic measures of the original mesh, we call them Intrinsic Parameterizations. In addition to this partial theoretical analysis, we pro- pose robust, efficient and tunable tools to obtain least-distorted parameterizations automatically. In particular, we give details on a novel, fast technique to provide an optimal mapping without fixing the boundary positions, thus providing a unique Natural Intrinsic Parameterization. Other techniques based on this parameterization family, designed to ease the rapid design of parameterizations, are also proposed.
1. Introduction
Parameterization is a central issue in graphics. Parameter- izing a 3D mesh amounts to computing a correspondence between a discrete surface patch and an isomorphic planar mesh through a piecewise linear function or mapping. In practice, this piecewise linear mapping is simply defined by assigning each mesh node a pair of coordinates(u,v) refer- ring to its position on the planar region. Such a one-to-one mapping provides a flat parametric space, allowing one to perform any complex operation directly on the flat domain rather than on the curved surface. This facilitates most forms of mesh processing, such as surface fitting, remeshing, or texture mapping. This last application, for instance, is widely used in Graphics as it dramatically enhances the visual rich- ness of a 3D surface, both for overly simplified character meshes in game engines as well as for incredibly detailed complex surfaces in computer-generated feature films. Un- fortunately despite numerous existing parameterization tech- niques8,22and commercial applications (Maya, Softimage, Flesh), it usually takes several hours of tedious(u,v) adjust- ments for a talented user to map a texture correctly (i.e., with acceptable distortions) onto an arbitrary surface.
This failure can be partially blamed on the intrinsic diffi- culty of the problem at hand. Since we are basically trying to flatten a surface from 3D down to 2D, there is in general no perfect way to perform such a flattening without introducing some form of distortion. However, most existing techniques, supposed to minimize distortion, do not result in a visually
“smooth parameterization”, even more so for very irregular meshes issued from scanners as demonstrated in Figure 3.
The goal of this paper is to lay out a theoretical back- ground on patch parameterization, as well as to develop a
ψ
Figure 1: A piecewise linear mapping between a 3D mesh M and an isomorphic flat meshU, where a triangle on the mesh is mapped to a triangle in the parameterization.
practical implementation. We will first define the notion of distortion measures between two 1-ring neighborhoods, ex- hibit the unique family of admissible distortion measures, derive from it a novel family of intrinsic parameterizations of surface patches, and propose fast algorithms to implement these theoretical results. We will also exhibit the practical re- sults one can get at low cost using our novel technique, such as the boundary-free flattenings depicted in Figure 2.
1.1. Problem Statement and Conventions
In this paper, we will deal with the following problem: Given a piecewise linear mesh patchM, possibly with holes but non-closed, construct a piecewise linear mappingψ between M and an isomorphic planar triangulation U ∈ IR2 that best preserves the original, intrinsic characteristics ofM. Throughout the paper, we will denote by xithe 3D position of the ith node in the original meshM, and by uithe 2D position (parameter value) of the corresponding node in the 2D meshU. We will also use the self-explanatory notation:
xi= (xi,yi,zi)t, ui= (ui,vi)t. Parameterizing a mesh is there- fore providing the piecewise linear mappingψ (see Figure 1)
The Eurographics Association and Blackwell Publishers 2002. Published by Blackwellc Publishers, 108 Cowley Road, Oxford OX4 1JF, UK and 350 Main Street, Malden, MA 02148, USA.
Desbrun, Meyer, and Alliez / Intrinsic Parameterizations such as:
ψ : M →U xi→ ui
Texturing the meshM will then be as simple as pasting a picture onto the parameter domain, and mapping each tri- angle of the original meshM with the part of the picture present within the associated triangle in the parameter plane.
Figure 2: Left: A 3D surface (top) and its natural conformal parameterization (bottom). Right: Views of the 3D surface.
1.2. Background
Due to its primary importance for any subsequent mesh ma- nipulation, the subject of mesh parameterization has been researched for a number of years, and not only in Computer Graphics.
Computer Graphics A significant body of work on pa- rameterization has been published over the last ten years in Computer Graphics. Almost all techniques explicitly aim at producing least-distorted parameterizations, and vary only by the distortions considered and the minimization pro- cesses used. Early work used the notion of flattening to ob- tain an isomorphic planar triangulation1,25,39, often min- imizing discrete variables in the process, such as the ra- tio of angles between the 3D triangles and their associ- ated 2D versions3,38,4. Others considered spring-like ener- gies19,28,2,20,8,21,30,9,24that can be minimized quickly by a linear system solver when the boundary has been fixed to an arbitrary contour (with the noticeable exception of22where only a few internal points need to be fixed by the user).
The Discrete Conformal Parameterization (DCP) has been proposed independently by a number of authors29,7,15who derived the same linear condition for conformality either us- ing Differential Geometry, harmonic maps, or Finite Ele- ments. Here again, a boundary condition is needed to induce a conformal mapping.
Finally, one can use nonlinear formulations to define an
optimal parameterization 16,31. The MIPS method for in- stance finds a “natural boundary" that minimizes their highly non-linear energy 16. Unfortunately, this requires quite a computational effort (even if hierarchical solvers can be used17,34) for a result visually very close to the DCP. Sander and co-authors35proposed yet another nonlinear energy for the specific problem of texture stretch distortion.
Most of these techniques proposed to minimize a con- tinuous energy over a piecewise linear surface. However, the choice of the energy sometimes seems very arbitrary, and most of them may visually result in non-smooth param- eterizations and therefore non-smooth textured meshes, as demonstrated in Figure 3(b-c). Note that using35results in a smooth parameterization, but takes more than six minutes to converge since it requires a non-linear minimization.
Cartography Concurrently, cartographers have been deal- ing with the parameterization of non-flat surfaces for cen- turies, in order to represent our spherical earth as flat maps.
Their work has mainly focused on differential parameteriza- tion, and is therefore only marginally relevant in Computer Graphics in practice. It is however interesting to mention that it is well known in this field that a mapping of a curved surface can either be authalic (i.e., area-preserving), or con- formal (i.e., angle-preserving). No mapping of the earth can be isometric (i.e., distance-preserving): as it would have to be both authalic and conformal, and this is strictly impos- sible for non-developable surfaces like a sphere and most other geometries. This paper establishes similar results, but for discrete surfaces, extending the known continuous dif- ferential geometry results as well as providing insights for novel notions.
1.3. Overview
In this paper, we restrict ourselves to the parameterization of non-closed triangulated surfaces since many existing papers already describe different techniques to split a closed object into a series of patches (also called atlas of charts13,10,35,23).
We demonstrate that the set of desirable mappings for such patches form a simple low-dimensional space (Section 2).
Moreover, the two generative parameterizations of this space are the existing discrete conformal mapping and a novel dis- crete authalic mapping, and all other valid intrinsic param- eterizations can be found by simply solving a sparse linear system as detailed in Sections 3 and 4. We also demonstrate that they generate smooth texture mappings even for highly irregular meshes. We then show how easily one can find an optimal parameterization without fixing boundary points, providing a natural parameterization, by simply adding nat- ural boundary conditions. Finally, we quickly review the possible immediate extensions that one could do with this new parameterization family before concluding.
2. Distortion Measures for 1-Rings
We want to preserve as much of the intrinsic qualities of a surface as we possibly can during its parameterization, i.e., its flattening. This implies that we need to first define what these intrinsic qualities are for a discrete surface: minimal
(a) irregular mesh (b) 19 (c) 8 (d)35 (e) DAP (f) DCP Figure 3: Intrinsic Parameterizations: Previous parameterization techniques (b-c) are not robust to mesh irregularity, exhibit- ing large distortions for highly irregular, yet geometrically smooth meshes such as in (a). Non-linear techniques (d) can achieve much better results, but often require several minutes of computational time. In comparison, with the exact same boundary con- ditions, our technique quickly generates very smooth parameterizations, regardless of the mesh irregularity (sampling quality) as demonstrated by the two texture-mapped members (e-f) of the novel parameterization family (denoted Intrinsic Parameteri- zations) that we introduce in this paper.
distortion will then mean best preservation of these quali- ties. In this section, we restrict our investigation to the dis- tortion measures between simple 1-ring neighborhoods, and demonstrate that the appropriate measures actually form a low dimensional space.
2.1. Notion of Distortion Measure
As in the problem statement, letM be a simple mesh em- bedded in 3D consisting of a 1-ring neighborhood (i.e., a vertex and all its adjacent triangles), and letU be an iso- morphic mesh:U∼M (we use the symbol∼ to indicate isomorphy). Figure 4 shows the mapping between two sim- ple 1-ring neighborhoods. We define a distortion measure betweenM andUas a functional E taking two isomorphic triangulations as inputs, and returning a real value:
E : T ×T → IR (M,U)→ E(M,U).
This kind of functional is sometimes referred to as a mu- tual energy, as it can be seen as a measure of the energy required to distort one into the other. By the very definition of a distortion measure, E(M,·) must be minimum forM, as there is no mesh less distorted compared toM than it- self.† We therefore have the following inequality for every Usuch thatU∼M:
E(M,M) ≤ E(M,U) (1) For convenience, we will denoteφ the distortion of a 1- ring with itself: E(M,M) = φ(M). Thus, φ is a measure (sometimes called energy) of the triangulated surface. In order to further investigate what the appropriate distortion measures are for 1-rings, we now explore what the possible measures of a mesh are, since it will restrict the possible set of distortion measures.
† Note however that there generally exist other meshes, different fromM, that also achieve the same energy minimum.
2.2. Properties of Intrinsic Measures
A measure of a mesh is a functional φ which, given a piecewise-linear surface patch M, basically returns a
“score"φ(M). This functional must satisfy a few basic prop- erties, that we now go over.
Rotation and Translation Invariance Obviously, we want the functional to be invariant to any translation or rotation of the mesh. Since these affine transformations do not affect the geometry of the mesh, the measure should remain identical. This will consequently render the param- eterization independent of rotation and translation of the input mesh.
Continuity We also want the functional to be a dis- crete version of a continuous measure, consistent with the continuous, differential case. Thus, the functional needs to converge to a continuous measure as we get a finer and finer triangulation, under some possible additional conditions (such as bounded fatness, or more generally, non-degenerated triangulations). This is called conditional continuity, and is usually stated as:
φ(Mn) → φ(M) ifMn→M as n→ ∞.
Here again, this will induce a very natural property for our parameterizations.
Additivity A measure should also be additive, i.e.:
φ(M1∪M2) + φ(M1∩M2) = φ(M1) + φ(M2).
The measure with such a property has the desirable qual- ity of being intrinsic, that is to say, it only depends on the surface itself, not on its sampling. To illustrate this fact, consider the addition of one or several vertices onto the existing surface (along the boundary or inside a triangle for instance); it is easy to verify that the functional will still return the same measurement, since the real geome- try of the surface is not affected – only its discretization, hence the term intrinsic. This sampling-independent prop- erty will be particularly attractive when dealing with large
The Eurographics Association and Blackwell Publishers 2002.c
Desbrun, Meyer, and Alliez / Intrinsic Parameterizations meshes, since hierarchical solvers will prove particularly
efficient in solving for the parameterization.
2.3. Admissible Intrinsic Measures
Although the restrictions imposed on the notion of measure seem to be loose and natural, there are, surprisingly, only a small family of functionals that meet the requirements.
Minkowski Functionals of 2-Manifolds A set of well- known functionals satisfies all the previous conditions.
These are called the Minkowski functionals. For 2D sur- faces, there are three such functionals: the AreaφA, the Eu- ler characteristicφχ, and the PerimeterφP (length of the boundary) of a triangle mesh. It is straightforward to check that each of these functionals meet the three conditions we just listed. It is also interesting to notice that the two first ones (the perimeter being only a boundary measure) corre- spond respectively to the integrals of the determinants of the first— area element — and of the second— Gaussian curva- ture —fundamental forms12. These are well-known to be in- trinsic in the differential geometry sense, meaning that they could be computed by “inhabitants” of the surface having no knowledge of the actual embedding of the surface.
Admissible Functionals Since we are looking for measures over a triangulation, a result dating back to the previous cen- tury explicitly states the set of all admissible functionals. A triangulation (considered as a 2-manifold, and disregarding its embedding) belongs to the convex ring, since it is the union of a set of triangles, therefore a union of convex bodies (which does not mean that the triangulation itself has to be convex). On this convex ring, Hadwiger14has proven that the only functionals, defined over the convex ring, matching the three conditions we mentioned above are linear combi- nations of the Minkowski functionals‡. Therefore, the only admissible functionals fitting the three previous properties are linear combinations of Area, Euler characteristic, and Perimeter. The set of all admissible functionals is therefore a 3-dimensional space, and for any admissible measureφ, there exists a unique triplet of constants c1, c2, c3such that:
φ= c1φA+ c2φχ+ c3φP. (2) Valid Distortion Measures Between 1-rings
Let’s go back to our measures of distortion between two iso- morphic 1-rings. Since the distortion measures must match the intrinsic measures, this restrains the admissible set to a special subset of the general case proven in 33, because of the additional additivity and continuity conditions. We show in the next section that the simplest relevant distortion mea- sures form a two-dimensional space.
3. Optimal 1-Ring Flattening
We now introduce the only quadratic distortion measures that fit the requirements that we derived in the previous sec- tions. We show that their critical points, when a boundary
‡ Hadwiger’s book has never been translated in English. There are however several books36 and papers37,18 that clearly state the aforementioned theorem.
condition is imposed, can be found by solving a simple lin- ear equation. We start by developing the two most represen- tative optimal mappings, the discrete conformal parameter- ization (DCP) and the novel discrete authalic parameteriza- tion (DAP), and demonstrate how all the others can be de- duced from their formulation. We also point to some of the similarities between the differential and the discrete case.
3.1. Notion of Optimal Vertex Placement
We call an optimal 1-ring parameterization any mapping from a given 3D 1-ringMto an isomorphic 2D 1-ringUthat is the minimum of a distortion measure (as previously de- fined) for a fixed, given boundary mappingψ(∂M). There- fore, if a distortion measure E is known and if each boundary vertex xjhas a given parameter value uj, the condition for the 2D 1-ring to be minimally distorted (i.e., optimal) is sim- ply that E(M,U) is minimum over allU∼M, which yields this simple condition for the center node uiof the 1-ring in the parameter plane:
∂E
∂ui
= 0.
We shall now describe what the appropriate energies E are that define a distortion measure between two meshes. The first one is known under the name of Dirichlet energy.
3.2. Discrete Conformal Mapping
The first optimal mapping is actually already known. We now recall a bit of history and background to demonstrate the connection to our problem, and later build upon it.
Conformality on Differential Surfaces: While working on the area minimization problem introduced by the belgian physicist Plateau, Mrs. Rado, Douglas, and later Courant proposed the use of the Dirichlet energy of a mapping in- stead of the highly nonlinear area functional previously used (see29for a good overview). The simple idea behind this functional12 is that, in differential geometry, the area of a patchM is:
Area = 12RM| fu× fv| dudv ≤12RM | fu| | fv| dudv
≤ 14RM( fu2+ fv2) dudv = Dirichlet energy.
It is simple to verify that the first inequality becomes an equality iff fu· fv= 0 everywhere, while the second does iff| fu| = | fv| (deriving from the positivity of ( fu− fv)2). A further analysis12shows that the minimum of this energy (quadratic in the parameterization) is the area, and is only attained for conformal mappings, i.e., mappings where the two previously introduced conditions on f hold. Conformal- ity of the map equivalently means angle preservation since these conditions imply that any angle between two vectors on the parameter plane will be preserved through the map- ping. In other words, in the differential case, it is known that a conformal map will result from the minimization of the Dirichlet energy.
Dirichlet Energy on Triangulations: Pinkall and Polthier provided a formal derivation of the Dirichlet energy between two triangles in29 for piecewise linear parameterizations.
Summing the energies over the whole 1-ring, they found:
EA=
∑
oriented edges(i, j)
cotαi j|ui− uj|2, (3)
where|ui− uj| is the length of the edge (i, j) in the parame- ter domain, andαi jis the opposite left angle in 3D as shown in Figure 4. This nicely complements the differential case since this is also a quadratic energy in the parameterization, and that this discrete energy depends only on the angles of the original surface. Indeed, the only term depending on the original surface is the cotangent term. This energy is also equal, at its minimum, to the total surface areaφA(M) when applied on the identity map (i.e., when ui jis taken to be the actual 3D edge), and is therefore the exact equivalent of the Equation 1: EA represents a distortion measure that fits our requirements defined in Section 2.
Critical Point of EA: Mimicking the differential case, Pinkall and Polthier29proposed to define the discrete con- formal map to be the critical point (a.k.a. the minimum) of the Dirichlet energy. Since this energy is quadratic, the derivation results in a simple linear system, that has a prov- ably unique solution which is easy to compute once we fix the boundaries in the parameter domain. Notice that we can not formally claim this mapping to be angle-preserving, since there is generally no way to flatten a curved, discrete surface with a one-to-one correspondence of the 3D angles to the 2D angles. However, since the Dirichlet energy de- pends only on the 3D angles, and that in the differential case, the minimum of the Dirichlet energy is indeed conformal, this definition results in a visually satisfying parameteriza- tion as depicted in Figures 3(f) and 6(a). This explains the success (and the name) of this discrete conformal parameter- ization (DCP). Due to the simple formulation of this energy, deriving its critical point is rather simple, yielding the linear equation for the central node i:
∂EA
∂ui
=
∑
j∈N(i)
(cotαi j+ cotβi j) (ui− uj) = 0. (4) Again, we can note that the linear coefficients are func- tions of only the angles of the original surface. We will describe in detail the computations required to numerically solve for the parameterization in Section 4, as well as an ex- tension to natural boundary conditions.
2D Analogy Consider a mesh vertex in flatland (2D) and its immediate neighboring vertices. Any motion of this ver- tex uiin the plane will preserve the 1-ring area, as men- tioned in5. Therefore, computing the gradient of the 1-ring area with respect to uiwill provide a nontrivial equation (the gradient for each triangle being nonzero) that does always sum to zero (since the total area is constant). Not surpris- ingly, Appendix A shows that we find the same coefficients as in Equation 4. Indeed,φA is the area of mesh, and there- fore, the coefficients of∂EA(M,M)/∂ximust match those of∂EA(M,U)/∂ui— giving an alternate, simple derivation of the conformality condition.
3.3. Discrete Authalic Mapping
Similarly to EA(M, ·) matchingφA on the identity map of M, we now discuss the existence of a novel quadratic en- ergy Eχthat matches the Euler characteristicχon the iden-
x
ψ
x
ij ij ij ij
δ α u
i
j
γ
ui
j
β
Figure 4: A 3D 1-ring, and its associated flattened version.
tity map while being a valid distortion measure. Despite the relative simplicity of the optimality condition, we did not find any mention of it for the differential or the discrete case in the vast literature available. We will show, however, that this new condition has smoothness qualities similar to those of the conformality condition.
Hands-on Derivation Remember that the Euler charac- teristic is the integral of the Gaussian curvature. It is known12,26that the Gaussian curvature, hence the determi- nant of the second fundamental form, is equal to 2π−∑jθj
where theθj’s are the tip angles around ui. Therefore, a sim- ilar gradient computation can be done for the sum of the tip angles around ui. Indeed, for a flat triangulation, this sum also remains constant (and is equal to 2π) as uimoves within the plane. This time, we get new coefficients as proven in Appendix B. From this simple derivation, we now derive an appropriate energy Eχin the next paragraph.
Chi Energy on Triangulations Guided by the previous derivation, we introduce the following quadratic energy:
Eχ=
∑
j∈N(i)
(cotγi j+ cotδi j)
|xi− xj|2 (ui− uj)2 (5)
where the anglesγi jandδi jare defined in Figure 4. This energy is constant for a given 1-ring when evaluated on the identity map, and therefore can always be scaled and shifted to be equal toχ(1 for a closed 1-ring). Additionally, since it is quadratic, we can show that it is greater than (or equal to)χ(after the above scaling and shifting) for any other map with the same boundary. This energy therefore satisfies all the properties we required in Section 2.
Critical Point of Eχ Once again, the optimal parameteri- zation deriving from Eχis easily obtained when the center node uisatisfies:
∂Eχ
∂ui =
∑
j∈N(i)
(cotγi j+ cotδi j)
|xi− xj|2 (ui− uj) = 0. (6) Duality of EA and Eχ Now, the coefficients of both this optimality condition and of Eχ are shown (also in the ap- pendix) to be only functions of local areas of the 3D mesh.
This should not come as a complete surprise: remember that the Dirichlet energy, which derives from the determinant of the first fundamental form, a measure of the local area ex- tension 12, depends only on local angles and provides an angle-preserving mapping when minimized. Sinceχis (up to a constant) the integral of the determinant of the second fundamental, which is a measure of the local angle excess12, we have a dual situation here. The energy Eχis now depend- ing only on local areas, and we therefore denote it Discrete
c
The Eurographics Association and Blackwell Publishers 2002.
Desbrun, Meyer, and Alliez / Intrinsic Parameterizations
Authalic Parameterization (DAP) for the same reasons as the DCP. Solving for the optimality of Eχ, using numerical methods described in Section 4, results in smooth parame- terizations as shown on Figures 3(e) and 6(b). Just like the DCP does for the angles, the DAP tries to preserve the area structure of the original 1-ring.
3.4. General Discrete Parameterization
As mentioned in Section 3, the set of all admissible measures are linear combinations of the Minkowski functionals. For fixed boundary conditions, the only distortion measures pos- sible are linear combinations of the area and the angle distor- tion measures (note: the perimeter distortion does not induce a particular position for the center vertex being a lower-order (1D) distortion measure for the boundary only). Therefore, it results that the family of admissible, simple distortion mea- sures of a 1-ring is reduced to linear combinations of the two discrete distortion measures defined above. A general distor- tion measure E as we defined can thus always be written as:
E=λEA+ µEχ
whereλ and µ are two arbitrary real constants. The opti- mality condition will simply be a linear combination of the two optimality condition we have described above. We call this 2-dimensional space of optimal discretizations Intrin- sic Parameterizations, since they naturally derive from in- trinsic measures of the input mesh. As demonstrated in Fig- ure 3, they provide smooth parameterizations even on highly irregular meshes since they minimize intrinsic distortions.
Caveat: Although the DCP can easily be proven to be glob- ally optimal (and therefore, angle-preserving when the tri- angulation is fine enough), the DAP is, as far as we know, only locally optimal. This means that we should not expect the DAP to perfectly preserve the area distortion across the mesh, but only as best as possible between each 1-ring.
3.5. Connection to Barycentric Coordinates
There exists a direct connection between barycentric coor- dinates and parameterization. It was already noted in 8,11 that the coefficients of the usual linear systems used to pa- rameterize meshes can be interpreted as barycentric coor- dinates of each internal vertex within its 1-ring. If the lin- ear system for parameterization really represents barycen- tric coordinates, then any flat mesh will be its own parame- terization, since each vertex will not move from its original position within its 1-ring. Although this condition seems to be an obvious quality for a “good” parameterization, only a few previous techniques satisfy this simple criterion. On the other hand, any linear combinations of the coefficients in Equations 4 and 6 defines perfectly valid barycentric co- ordinates27. We additionally proved that there is no other possible barycentric coordinates with the same properties, due to Hadwiger’s theorem.
4. Parameterizing Meshes
We now discuss the practical implementation of the theoret- ical results presented above, along with convenient improve- ments to further aid in the design of good parameterizations.
We first give details on how to solve for the least-distorted parameterizations with a fixed boundary, then show how to interactively move the boundaries to further reduce the bur- den of designing nice texture mapped surfaces, and finally present a natural parameterization that automatically finds an optimal boundary.
4.1. Computing an Intrinsic Parameterization
Since the gradients of the energies introduced in Section 3.4 are linear, computing a parameterization reduces to solving a sparse linear system:
MU=
λMA+ µMχ
0 I
Uinternal Uboundary
=
0 Cboundary
= C
where U is the vector of 2D-coordinates to solve for (sep- arated for convenience into the internal vertices and the boundary vertices); C is a vector of boundary conditions that contains the positions where the boundary vertices are placed; and MA and Mχare sparse matrices whose coeffi- cients are given respectively by:
MAi j=
cot(αi j) + cot(βi j) if j∈N(i)
−∑k∈N(i)MAik if i= j
0 otherwise,
Mχi j=
cot(γi j) + cot(δi j) /|xi− xj|2 if j∈N(i)
−∑k∈N(i)Mχik if i= j
0 otherwise.
Note that this technique can handle an arbitrary number of boundary curves (they are simply additional boundary ver- tices) and therefore easily parameterize patches containing holes. Once the boundary points have been chosen (either automatically or by the user), the sparse system is efficiently solved using Conjugate Gradient with an appropriate pre- conditioning (we recommend SSOR or inverse diagonal pre- conditioning — see32).
Constraints The user may possibly want to constrain cer- tain points to given parameter values. This can be easily achieved using Lagrange multipliers. Each point constraint creates a linear equation relating the parametric values of the vertices of the enclosing triangle (using triangular barycen- tric coordinates) to the constrained position. We then add these additional constraints to the linear system using stan- dard Lagrange multipliers. The previous linear system is then augmented to the following system:
M (Mη)T
Mη 0
U η
=
C Cη
where Mηi jis 1 only if the jthconstraint constrains the ithver- tex and 0 otherwise, and Cηjis set to the jthconstrained posi- tion. Note that constraining a line is also possible by simply constraining the endpoints as well as the intersections of the line with the edges of the mesh.
Figure 5: Other Examples of Natural Conformal Maps: to demonstrate the conformality of the maps we obtain, we use an irregularly sampled mesh and observe that the symmetry is preserved despite the drastic change in sampling rate. The third natural parameterization uses the same mesh as in Figure 3. These four parameterizations were obtained in 0.8 s, 0.5 s, 1.8 s, and 0.3 s, respectively. See also color section.
4.2. Modifying Boundaries
In addition to simple fixed-boundary conditions, we can al- low the user to interactively modify the positions of bound- ary points while updating the parameterization in realtime.
This efficiency is achieved by taking advantage of the linear nature of our solution and precomputing how the parameter- ization responds to the movement of a boundary point.
Impulse response We first precompute the parameteriza- tions that result from placing one boundary point at(1,1) and all others at(0,0) (these correspond to the Green func- tions of our parameterization equation):
Mbi= ei, ∀i ∈ boundary
where eiis a 1D vector (i.e., a vector of scalars) containing 1 in ithposition and 0 elsewhere, while biis the unknown 1D vector.
Realtime Boundary Manipulation By solving this system once for every boundary point, we construct a set of “basis parameterizations” that describe how the parameterization is altered by a change in a single boundary position. Indeed, the parameterization can then be efficiently updated as the user manipulates the boundary by noting that:
C=
i∈boundary
∑
(uboundaryi )Tei
= M
i∈boundary
∑
(uboundaryi )Tbi
where uboundaryi is the position of ithboundary point. There- fore, the parameterization for a given set of boundary points can easily be reconstructed in realtime, allowing realtime boundary manipulation, as:
U=
∑
i∈boundary(uboundaryi )Tbi.
This novel feature provides an easy tool for a user to op- timize the design of a texture mapping on arbitrary surface patches.
Natural Boundaries / Natural Conformal Map
Interestingly, we can also solve a similar linear system while letting the computer pick the “best” boundaries. Earlier, we showed how to get a parameterization once a boundary was given, but we can also solve for an optimal conformal map- ping by imposing natural boundaries (also called Neumann boundaries). This requires only minor modifications to the prior algorithm, and due to the quadratic nature of the en- ergy, we will also obtain a unique solution. As demonstrated in Appendix A, we show that the derivative of the Dirich- let energy on a triangle with respect to one of its vertices is equal to the opposite edge rotated by 90 degrees (such that (x,y) → (−y,x)). A natural boundary condition is therefore to have the same property at the boundary. Summing over all adjacent triangles, the equation for the boundary point i (which we place into the matrix M) becomes:
∆
∑
i jkcotα(ui− uj) + cotβ(ui− uk) =
∑
∆i jk
R90(uk− uj) (7) Whereα and β are the angles at k and j, and R90 is a ro- tation by 90◦. Note that this property also holds for interior vertices as the terms on the left become the conformal con- dition (Equation 4) and the terms on the right sum to zero.
To complete the minimization, we need to constrain two vertices to fix the rotation and translation of the resulting minimum parameterization. In our implementation, we con- strain the two boundary vertices the farthest from each other to two arbitrary positions in the parameterization plane. This
The Eurographics Association and Blackwell Publishers 2002.c
Desbrun, Meyer, and Alliez / Intrinsic Parameterizations simple modification results in Natural Conformal Maps,
such as those depicted on Figures 2 and 5. Notice that these parameterizations take the same amount of time to compute as the fixed-boundary ones, offering a very nice tool for ini- tial flattening before minor editing of the boundaries if nec- essary. Note also that if the authalic coefficients are proven, in the future, to derive from a global energy, a similar treat- ment can be used to find natural boundaries.
5. Nonlinear Optimization of Maps
The theoretical and practical work introduced in this paper opens many avenues. Aside from the natural conformal map and the entire family of intrinsic parameterizations by vary- ing the parametersλand µ, we can also compute parame- terizations sufficiently close to these optimal ones to be vi- sually smooth, but potentially more appropriate for a given application. Therefore, we propose two simple algorithms to compute good parameterizations that minimize other types of functionals.
5.1. Near-Optimal Maps
We sometimes wish to minimize highly nonlinear energies while remaining within the space of the aforementioned in- trinsic parameterizations. In order to make this tractable, we can linearize the solution space by assuming that all solu- tions can be expressed as a linear combination of the two base intrinsic parameterizations:
U=λUA+ (1 −λ)Uχ. (8) Note that µ= 1 −λto ensure that the solution always inter- polates the same boundary as we varyλ. Since we restrained the vector space of solutions to only linear combinations of the intrinsic parameterizations, many nonlinear functionals can be minimized by a simple low-order polynomial mini- mization (often in real time). Below are two simple examples of such functionals.
Edge-Length Distortion Minimization is achieved by minimizing the nonlinear energy:
E=
∑
i j∈Edges
(|ui− uj|2
|xi− xj|2− 1)2.
By substituting the values for ui and uj from equation 8, the energy becomes a quartic polynomial inλ. This energy can then be minimized in real time using a 3rd order poly- nomial root finder to solve for: dEdλ =∑i j∈Edges4(|u|xi−uj|2
i−xj|2−
1)(ui−uj)·[(uAi−u
χ
i)−(uAj−uχj)]
|xi−xj|2 = 0.
Area Distortion Minimization can be achieved by mini- mizing: E=∑i jk∈Faces((A
param i jk
Ai jk3D )2−1)2. As in the edge length distortion minimization, this results in a quartic polynomial in lambda and can be efficiently solved using a simple root finder. An example resulting from this technique is depicted in Figure 6
5.2. Boundary Optimization
Similar to the way we minimized an energy by modifyingλ, we can, alternatively, minimize the energy by modifying the
Conformal Authalic Optimal lambda
Figure 6: Area Distortion Minimization can be achieved by optimizing the linear combinationλUA+ (1 −λ)Uχof the conformal and authalic parameterizations. The parameteri- zations (top) and the area distortion pseudo-coloring (mid- dle) demonstrate the quality of the optimization (see color section).
boundary of the parameterization. We first choose an appro- priate energy (edge length distortion, area distortion, etc.), and then take its derivative with respect to each of the bound- ary points. Note that the terms of the form∂ui/∂uboundaryp can be (pre)computed using the impulse response technique de- scribed in section 4.2. These derivatives are then used to per- form a gradient descent to find a local minimum of the speci- fied energy. Since the gradient descent is performed in terms of boundary points only (much fewer than the total number of points), this process is very efficient, and takes generally less than 10 seconds for several hundreds boundary vertices.
A sequence of boundary optimizations using this method is depicted in Figure 7.
6. Results and Conclusions
All the results we obtained were computed in less than 5 sec- onds for fixed-boundary and natural parameterizations, and less than 15 seconds for boundary-optimized maps. These techniques are thus very well suited for interactive parame- terization design. In summary, we have introduced the novel family of Intrinsic Parameterizations. We showed that they are the only parameterizations satisfying the proper con- ditions to make them easy to compute and robust to arbi- trary meshes with guaranteed smoothness quality. We have also proposed algorithms to exploit these parameterizations and automatically design optimal maps, with or without boundary conditions. However, many other extensions can be thought of. Using least-squares as in24, one could define other parameterizations very easily, based on the two main sets of coefficients introduced in this paper. Basically, any parameterization obtained by a linear system close to the co-
Figure 7: Boundary Optimization: after choosing a (non)linear functional to minimize over the parameterization, we can move the boundary points to perform a gradient descent and optimize the parameterization. Here, an initial irregular spherical strip is mapped to a circle, then evolves towards an optimized parameterization (1.5 s) minimizing edge-length distortion.
efficients of our Intrinsic Parameterizations will be visually smooth, providing us with a great deal of freedom to build other algorithms for parameterizations if more complex con- straints must be enforced.
Independently, Levy23has developed an alternate deriva- tion of the Natural Conformal Maps, in which he proves that the parameterization will not contain any folds (overlapping triangles). Therefore, future work may show that the Natu- ral Conformal Maps are a discrete analog of the "Riemann Mapping Theorem".
Additional future work will focus on clarifying the rela- tionship between our results and the existing body of work in Circle Packing, a technique which also provides the same kind of discrete mapping, but at the cost of a computationally expensive iterative process. Developing a good hierarchical solver as in6and in34could also speed up the process, mak- ing parameterization of extremely large meshes tractable.
Finding optimal charts on a closed surface to locally param- eterize a whole geometry is also of interest.
Acknowledgements The authors are extremely grateful to Mar- tin Rumpf for his insightful comments. Thanks to David Cohen- Steiner for pointing us to Hadwiger’s work, to Peter Schröder for support and instant German translation, to Jerrold Marsden, Anil Hirani and Andrei Khodakovsky for math support, to Tony DeRose for interesting discussions, to Isaac Cohen for pointers on differen- tial geometry, and to the reviewers for extremely insightful com- ments. We also greatly benefited from a discussion with the authors of35. The work reported here was supported in part by IMSC, an NSF Engineering Research Center (EEC-9529152), and a NSF CA- REER award (CCR-0133983).
References
1. BENNIS, C., VÉZIEN, J.-M., IGLÉSIAS, G., AND GAGA- LOWICZ, A. Piecewise Surface Flattening for Non-Distorted Texture Mapping. Computer Graphics (Proceedings of SIG- GRAPH 91) 25, 4 (July 1991), pp.237–246.
2. CAMPAGNA, S.,ANDSEIDEL, H.-P. Parameterizing Meshes with Arbitrary Topology. In Image and Multidimensional Dig- ital Signal Processing (1998), B. G. H. Niemann, H.-P. Seidel, Ed., pp.287–290.
3. C.L.WANG, C., CHEN, S.-F., FAN, J., AND M.F.YUEN, M. Two-dimensional Trimmed Surface Development using a Physics-based Model. Proceedings of the 25th Design Au- tomation Conference (Sept. 1999). Paper No. DETC99/DAC- 8634.
4. C.L.WANG, C., CHEN, S.-F.,ANDM.F.YUEN, M. Surface Flattening based on Energy Model. Computer Aided Design (to appear).
5. DESBRUN, M., MEYER, M., SCHRÖDER, P.,ANDBARR, A. H. Implicit Fairing of Arbitrary Meshes using Diffusion and Curvature Flow. In Proceedings of SIGGRAPH 1999 (1999), pp.317–324.
6. DUCHAMP, T., CERTIAN, A., DEROSE, T.,ANDSTUETZLE, W. Hierarchical Computation of PL Harmonic Embeddings".
Technical Report (July 1997).
7. ECK, M., DEROSE, T., DUCHAMP, T., HOPPE, H., LOUNS- BERY, M.,ANDSTUETZLE, W. Multiresolution Analysis of Arbitrary Meshes. Proceedings of SIGGRAPH 95 (August 1995), pp.173–182.
8. FLOATER, M. S. Parametrization and Smooth Approximation of Surface Triangulations. Computer Aided Geometric Design 14, 3 (1997), pp.231–250. ISSN 0167-8396.
9. FLOATER, M. S.,ANDREIMERS, M. Meshless parameteri- zation and surface reconstruction. Computer Aided Geometric Design 18, 2 (March 2001), pp.77–92.
10. GARLAND, M., WILLMOTT, A.,ANDHECKBERT, P. S. Hi- erarchical Face Clustering on Polygonal Meshes. 2001 ACM Symposium on Interactive 3D Graphics (March 2001), pp.49–
58.
11. GOTSMAN, C., AND SURAHHSKY, V. Guaranteed Intersection-free Polygon Morphing. Computer and Graph- ics 25, 1 (2001), pp.67–75.
12. GRAY, A., Ed. Modern Differential Geometry of Curves and Surfaces. Second edition. CRC Press, 1998.
13. GRIMM, C. M.,ANDHUGHES, J. F. Modeling Surfaces of Arbitrary Topology using Manifolds. Proceedings of SIG- GRAPH 95 (August 1995), pp.359–368.
14. HADWIGER, H. Vorlesungen Über Inhalt, Oberfläche und Isoperimetrie. Springer-Verlag, 1957.
15. HAKER, S., ANGENENT, S., TANNENBAUM, A., KIKINIS, R., SAPIRO, G.,ANDHALLE, M. Conformal Surface Pa- rameterization for Texture Mapping. IEEE Transactions on Visualization and Computer Graphics 6, 2 (April-June 2000), pp.181–189.
16. HORMANN, K., AND GREINER, G. MIPS: An Efficient Global Parametrization Method. In Curve and Surface De- sign: Saint-Malo 1999 (2000), P.-J. Laurent, P. Sablonnière, and L. L. Schumaker, Eds., Vanderbilt University Press, pp.153–162.
c
The Eurographics Association and Blackwell Publishers 2002.
Desbrun, Meyer, and Alliez / Intrinsic Parameterizations 17. HORMANN, K., GREINER, G.,ANDCAMPAGNA, S. Hierar-
chical Parametrization of Triangulated Surfaces. In Proceed- ings of Vision, Modeling and Visualization (1998), H.-P. S.
B. Girod, H. Niemann, Ed., pp.219–226.
18. K. MICHIELSEN, H. DERAEDT, J. F. Morphological Char- acterization of Spatial Patterns. Prog. Theor/ Phys. Suppl. 138 (2000), pp.453–548.
19. KENT, J. R., CARLSON, W. E.,ANDPARENT, R. E. Shape Transformation for Polyhedral Objects. Proceedings of SIG- GRAPH 92 (July 1992), pp.47–54.
20. KRISHNAMURTHY, V.,ANDLEVOY, M. Fitting Smooth Sur- faces to Dense Polygon Meshes. Proceedings of SIGGRAPH 96 (August 1996), pp.313–324.
21. LEE, A. W. F., SWELDENS, W., SCHRÖDER, P., COWSAR, L.,ANDDOBKIN, D. MAPS: Multiresolution Adaptive Pa- rameterization of Surfaces. Proceedings of SIGGRAPH 98 (July 1998), pp.95–104.
22. LÉVY, B. Constrained Texture Mapping for Polygonal Meshes. Proceedings of SIGGRAPH 2001 (August 2001), pp.417–424.
23. LÉVY, B.,ANDMAILLOT, J. Least Squares Conformal Maps for Automatic Texture Atlas Generation. ACM SIGGRAPH Proceedings (July 2002).
24. LÉVY, B.,ANDMALLET, J.-L. Non-Distorted Texture Map- ping for Sheared Triangulated Meshes. Proceedings of SIG- GRAPH 98 (July 1998), pp.343–352.
25. MAILLOT, J., YAHIA, H., ANDVERROUST, A. Interactive Texture Mapping. Proceedings of SIGGRAPH 93 (August 1993), pp.27–34.
26. MEYER, M., DESBRUN, M., SCHRÖDER, P.,ANDBARR, A. H. Discrete Differential-Geometry Operators for Triangulated 2-Manifolds, 2002. submitted, found at http://www.multires.caltech.edu/pubs/difGeoOps.pdf.
27. MEYER, M., LEE, H., BARR, A. H., ANDDESBRUN, M.
Generalyzed Barycentric Coordinates for Irregular N-gons.
Journal Of Graphic Tools (to appear) (2002).
28. PARIDA, L.,ANDMUDUR, S. Constraint-satisfying Planar Development of Complex Surfaces. Computer Aided Design 25, 4 (April 1993), pp.225–232.
29. PINKALL, U.,ANDPOLTHIER, K. Computing Discrete Min- imal Surfaces. Experimental Mathematics 2, 1 (1993), pp.15–
36.
30. PRAUN, E., FINKELSTEIN, A.,ANDHOPPE, H. Lapped Tex- tures. Proceedings of SIGGRAPH 2000 (July 2000), pp. 465–
470.
31. PRAUN, E., SWELDENS, W.,ANDSCHRÖDER, P. Consistent Mesh Parameterizations. Proceedings of SIGGRAPH 2001 (August 2001), pp.179–184.
32. PRESS, W., FLANNERY, B., TEUKOLSKY, S.,ANDVETTER- LING, W. Numerical Recipes in C: the Art of Scientific Com- puting, 2nd ed. Cambridge University Press, 1994.
33. RUMPF, M. A Variational Approach to Optimal Meshes. Nu- mer. Math., 72 (1996), pp.523–540.
34. SANDER, P., GORTLER, S., SNYDER, J.,ANDHOPPE, H.
Signal-Specialized Parameterization. MSR Technical Report MSR-TR-2002-27 (2002).
35. SANDER, P. V., SNYDER, J., GORTLER, S. J.,ANDHOPPE, H. Texture Mapping Progressive Meshes. Proceedings of SIG- GRAPH 2001 (August 2001), pp.409–416.
36. SANTALO, L. A. Integral Geometry and Geometric Probabil- ity. Addison-Wesley, 1976.
37. SCHMALZING, J.,ANDKERSCHER, M. Minkowsky Func- tionals in Cosmology. Generation of Large-Scale Structure in Cosmology (1997), pp.255–260.
38. SHEFFER, A.,AND DESTRULER, E. Surface Parameteriza- tion For Meshing by Triangulation Flattening. In Proceedings of the 9th International Meshing Roundtable, Sandia National Laboratories (Oct. 2000), pp.161–172.
39. SUN, M., AND FIUME, E. A Technique for Constructing Developable Surfaces. Graphics Interface ’96 (May 1996), pp.176–185.
Appendix A: Gradient of Area
A proof of the gradient of a triangle area with respect to a vertex, valid for an arbitrary embedding, can be found in 5. It was proven that for a triangle (P, A, B) we get:
α
β ε1ε
A P
B
H
2 ∇A=1
2((cotα) AP + (cotβ) BP) where∇denotes the gradient with re- spect to P. Additionally, since the area is equal to|AB| times the height |PH|, we have another simple ex- pression for the gradient (where “⊥” indicates a 90 degrees counter- clockwise rotation about the triangle’s normal):
∇A= |AB|∇(|PH|) = |AB| PH
|PH|= AB⊥. Appendix B: Gradient of Angle
Despite an extensive literature search, we have not found any pub- lished derivation for the gradient of one of a triangle’s angles with respect to its associated vertex. We therefore describe our derivation here.
LetT = (P, A, B) be a triangle, and let H be the orthogonal pro- jection of P onto the segment AB. We denote byε1the angle [APH, ε2the angle [HPB, andεthe angle ofT at P. Finally, we denote by αthe angle at A andβthe angle at B.
The gradient ofεwith respect to P can be decomposed into the sum of the gradients ofε1 andε2. Using the relation cos(ε1) =
|PH|/|PA|, the gradient can be computed as:
∇ε1 = arccos(|PH|/|PA|) = −|AH||PA|∇(|PH||PA|) (9)
= −|AH||PA|∇(|PH|) |PA|−∇(|PA|) |PH|
|PA|2 (10)
From the following identities: ∇|PA| = AP/|PA|, ∇|PH| = HP/|PH|, HP = HA + AP, and cotα= |AH|/|PH|, we obtain:
∇ε1 = −|AH||PA|h
HP
|PA| |PH|−|PA||PH|3AP i
=|PA|PA2 |PA|2
|AH| |PH|−|PH||AH|
+|PH| |AH|AH
=|PA|cotα2PA+|PH| |AB|AB (11) The gradient ofε2will cancel out the last term, leading to the simple formula:
∇ε= cotα
|PA|2PA+ cotβ
|PB|2 PB
Notice that the vector weights can be expressed only in terms of local areas: if K is the orthogonal projection of B onto PA, then cotα/|PA|2is equal to the area of the triangle(A, B, K) divided by twice the square of the total area of triangleT.
Summing the contribution due to each triangle of a 1-ring, we obtain, withθthe total angle around xi(in the notation of figure 4):
∇θ=
∑
j∈N1(i)
(cotγi j+ cotδi j)
||xi− xj||2 (xj− xi).
