- Email: [email protected]
A unified approach for fairing arbitrary polygonal meshes
Graphical Models 66 (2004) 160–179 www.elsevier.com/locate/gmod
A unified approach for fairing arbitrary polygonal meshes Guiqing Li,a,b,c,* Hujun Bao,a and Weiyin Mab a
The State Key Lab of CAD&CG, Zhejiang University, Hangzhou, China Department of Manufacturing Engineering and Engineering Management, City University of Hong Kong, Kowloon, Hong Kong, China c School of Computer Science and Engineering, South China University of Technology, Guangzhou, China
b
Received 31 May 2002; received in revised form 17 January 2004; accepted 5 March 2004 Available online 6 May 2004
Abstract Arbitrary polygonal meshes play an important role in 3D model representation and surface modeling. Most of the previous smoothing approaches are, however, mainly designed for triangular meshes. This paper proposes a unified approach to address fairing problems of arbitrary polygonal meshes. The approach is built upon an assumption that all face centroids of a smoothed mesh locate on its noise-free surface and remain unchanged during the smoothing process in order to prevent the resulted mesh from shrinkage and preserve surface features. Boundary smoothing is considered in a similar way as well but the process is independent to mesh processing such that two noisy meshes sharing a boundary still adjoin continuously after smoothing. Finally, a Taubin-like inflation operator is employed to further improve the ability of anti-shrinkage. Although very simple iterative equation is employed in the smoothing procedure, experiments demonstrate that the approach is able to not only effectively address the problem of shrinkage and distortion but also preserve features very well. 2004 Elsevier Inc. All rights reserved. Keywords: Arbitrary mesh smoothing; Laplacian operator; Mesh modeling; Mesh split
*
Corresponding author. E-mail addresses: [email protected] (G. Li), [email protected] (H. Bao), [email protected] (W. Ma). 1524-0703/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.gmod.2004.03.001
G. Li et al. / Graphical Models 66 (2004) 160–179
161
1. Introduction With advancement of 3D data acquisition techniques such as laser scanners, 3D photography, and CT-scanning, very large meshes are commonly used as surface representation for graphics applications today for their simplicity and flexibility. However, digitization processes often introduce noise into mesh attributes such as vertex position due to some physical factors associated with measurement, thus removing noise is required before performing any other operations on the model. On the other hand, to yield free-form surfaces for fairness, it is also necessary to smooth control nets automatically. The fairness of a surface can usually be measured by certain energy functional of geometric attributes such as curvature and curvature variation. The essence of surface smoothing is therefore directly or indirectly to minimize these energies [10,11,20,35,36]. Noticing geometric curvature quantities can be approximated by combination of second derivatives under the assumption that a parametric surface is rather flat [13], some authors employed all kinds of lower or higher derivatives to replace geometric intrinsics for fairing surfaces. Minimizing the membrane or thin plate energy, for example, can lead to surfaces of excellent quality [13,16,26,37]. However, approaches directly involving functional minimization often exhibit low efficiency due to expensive computational costs. Fortunately, as a variational problem can usually be transformed into PDEs, relaxation iterative techniques can often be employed to build an iterative process for finding numerical solutions of the PDEs in which each iterative step moves vertex positions by a small increment. Approaches based on the Laplacian operator fall into this category. Taubin [30] pioneered such a kind of algorithms by reducing the problem of mesh smoothing to low-pass filtering. Actually, most of the previous smoothing techniques fall into this framework [4,9,14,18,19,21,34,38,39]. Some of them will be further discussed in Section 2. To produce surfaces of high aesthetic quality, Schneider and co-worker [27,28] proposed a completely intrinsic approach based on PDEs that no longer comes from functional minimization. The approach is able to reproduce surfaces with constant mean curvature such as spheres, cylinders, and minimal surfaces. Based on the linear heat equation, Clarenz et al. [3] formulated an anisotropic diffusion approach for enhancing geometric features such as creases and corners of meshes smoothed. Bajaj and Xu generalized the above method to a two-dimensional manifold embedded in Rk so as to put vertex positions of the meshes and functions defined on the meshes into a unified smoothing scheme in which the local analytic form of the Loop subdivision was employed to represent the mesh surfaces and associated vector functions. Statistics-based algorithms were pioneered by Peng et al. [24]. Their approach is limited to smooth semi-regular meshes and pre-processing is needed to convert the meshes to multiresolution representations before fairing. The Weiner filter is then applied to smooth the directional detail coefficients. Denoised results are finally reconstructed from the smoothed details. Recent statistics-based approaches were motivated by the successful application of bilateral filtering to image smoothing [1,6,33]. Based on robust statistics and local first-order predictors of surfaces, Jones
162
G. Li et al. / Graphical Models 66 (2004) 160–179
et al. [12] proposed a non-iterative feature-preserving approach. With this method, face normals should be evaluated first since prediction is made via projecting the estimated vertex onto mesh faces. Bilateral mesh denoising proposed by Fleishman et al. [7] also employed the same technique but provided a simpler predictor. They evaluated the offset from estimated vertex to its underlying position, i.e., the predictor, by casting its neighbors onto its normal line. For large-scale noise, it is not easy to obtain good estimation of the underlying normals even if the surface normals are first smoothed as in [12]. In addition, Ohtake et al. [23] smoothed face normals of the noisy mesh using Gaussian filter. The result was then obtained by fitting the faired normals. By introducing anisotropic neighborhoods, their approach was adapted to smooth models with sharp feature. Though the method was argued fully automatic, the standard deviation r of Gaussian-like functions and the positive constant k of the modified distance have to be determined through experiment for smoothing face normals. The method of Tasdizen et al. [29] bore the same framework, but the anisotropic diffusion of normals was realized by minimizing the total curvature of normal maps instead of statistics-based. However, all the above-mentioned approaches aim at addressing smoothness of triangular mesh models while in practice arbitrary polygonal meshes are often encountered in some applications [17] such as surface design based on control nets. One may convert a polygonal mesh into a triangular mesh and then apply smoothing algorithms developed for triangle meshes. It is, however, not appropriate to do so and it is better to directly remove the noise before any processing similar to the case for processing triangle meshes. This paper addresses the smoothing of arbitrary polygonal meshes. To prevent shrinkage and preserve detailed features, the smoothing algorithm is designed such that all face centers pffiffiffi remain unchanged after smoothing. The noisy model is first subdivided with a 3 split operator by inserting face centers as new vertices. Then the refined mesh is smoothed under the constraint that the newly inserted vertices remain unchanged. The denoised mesh is finally obtained by replacing vertex positions of the noisy model with those corresponding vertices of the smoothed refined mesh. As the approach can deal with arbitrary polygonal meshes, triangular meshes can also be processed under this unified framework. Main contributions of this paper can be summarized as follows: • Arbitrary polygonal mesh smoothing: An approach is proposed to smooth arbitrary polygonal meshes. In the algorithm, a very simple linear combination based on Catmull–Clark subdivision scheme is explored to move vertices of the noisy mesh in the iterative process. • Boundary smoothing: A constrained curve smoothing approach is developed by interpolating midpoints of boundary edges such that the boundary can be smoothed simultaneously. • Convergence treatment: A formal analysis on convergence of the algorithm is discussed. • Inflation operation: A Taubin-like operator is designed to further enhance the anti-shrinkage effect.
G. Li et al. / Graphical Models 66 (2004) 160–179
163
The rest of the paper is organized as follows. In Section 2, some closely related work is reviewed. Our smoothing approach is subsequently described in Section 3 in further details. Some mathematical treatment is discussed in Section 4. To further improve the effect of anti-shrinkage, we introduce an inflation operator similar to that of [30] in Section 5. Comparison between the presented approach and two other algorithms is discussed in Section 6. Finally, conclusions are drawn in the last section.
2. Notations and related work In this section, we first formulate some notations and review several related approaches on mesh smoothing. 2.1. Formal conventions An arbitrary polygonal mesh with vertex set V ¼ f1; 2; . . . ; ng is represented by M ¼ ðK; P Þ, where K 2V (2V expresses the set consisting of all subset of V ) is a simplicial complex consisting of subsets of V and P is a mapping from V to coordinates of its vertices, i.e., P : V ! R3 . Elements of K containing only one integer, two integers are called vertices and edges of M, respectively. All other elements are defined as faces of M. For instance, fig indicates vertex with index i while fi; jg refers to edge with endpoints fig and fjg, and fi; j; kg represents a triangle with vertices fig, fjg, and fkg. Generally, small letters such as e or f are used to represent edges or faces. We will also use vi to indicate P ðiÞ, the coordinates of fig in the rest of this paper. If fi; jg 2 E, then fig and fjg are said adjacent or neighboring. Set consisting of all adjacent vertices of fig is called its 1-neighborhood in associated with mesh M denoted by NM ðiÞ while the number of its neighbors is called valence and read as jNM ðiÞj. Accordingly, BNM ðiÞ is used to express the set of the boundary adjacent vertices of a boundary vertex fig. The set of all faces sharing vertex fig is defined as the face 1-neighborhood of fig and denoted as FM ðiÞ. In addition, the degree of a face is defined as its vertex number. We also assume that meshes are 2D manifold. No more than two faces share an edge, each vertex belongs to at least an edge and each edge belongs to at least a face. 2.2. Related work Let f ðx; yÞ be a differentiable function and D be the Laplacian operator Df ¼ fxx þ fyy . A variational problem of membrane energy can be described as follows Z min ðfx2 þ fy2 Þ dA; f
where dA is an area element on the surface z ¼ f ðx; yÞ. From variational theory we know that the differential form of the above variational problem is the Euler– Lagrange equation
164
G. Li et al. / Graphical Models 66 (2004) 160–179
Df ¼ 0: For any vertex vi of an irregular mesh, the simplest discrete form of the Laplacian operator can be expressed as X Dvi ¼ vi þ wij vj ; ð1Þ j2NM ðiÞ
where wij ¼
1 : jNM ðiÞj
ð2Þ
One can then use the following relaxation formula to iteratively find the approximate solution of the discrete Euler–Lagrange equation Dvi ¼ 0: ¼ vold þ kDvi ; vnew i i
ð3Þ
where k is a positive relaxation factor. We refer the smoothing method based on Eq. (3) as Laplacian smoothing. Notice that Eq. (2) defines uniform weights for all neighboring vertices. In literature, one may find an approach that uses the reverse of the edge length to replace the weights of Eq. (2) [4,8]: X 1 1 wij ¼ : ð4Þ jvi vj j j2N ðiÞ jvi vj j M
One may also find a sophisticated set of weights that is related to the surface curvature, an intrinsic attribute of the mesh geometry [4,5,32]: wij ¼ P
ðcot aij þ cot bij Þ ; j2NM ðiÞ ðcot aij þ cot bij Þ
ð5Þ
where aij and bij are the two angles opposite to edge ði; jÞ. Many variants of Laplacian smoothing were proposed due to the fact that the Laplacian smoothing may introduce serious shrinkage and distortion and result in over-smoothing. Some of the approaches involve Laplacian operator of higher order [30] and some others use constraint-based approaches [4,18,19,34]. Based on techniques for signal processing and eigenvalue analysis of the matrix derived from the Laplacian smoothing formula, Taubin found that the shrinkage problem could be eliminated by iteratively computing Eq. (3) following an inflation step vnew ¼ vold þ lDvi ; i i
ð6Þ
where l is a negative relaxation factor satisfying l > k. However, when vi is located on the boundary, boundary shrinkage cannot be avoided due to a very strong tangential component of Dvi . Though successful in preventing shrinkage, TaubinÕs scheme runs into difficulty in preserving features [21,22]. Taubin [32] addressed boundary shrinkage by replacing Dvi with its normal component in which the normal can be either estimated by averaging the normals of adjacent faces or specified as constraints by users. Though it is able to maintain the same convergent speed for
G. Li et al. / Graphical Models 66 (2004) 160–179
165
internal and boundary vertices using the same iterative framework, the approach may produce cracks on boundaries shared by two meshes. Another variant was proposed in [34]. The method also includes two steps—a variation of the Laplacian and a so-called HC-modification. The first step smoothes the mesh by moving all vertices of the old mesh based on the Laplacian smoothing criteria and the second step pushes the vertices of the resulting mesh of the first step back towards their original positions by certain distance. The approach is able to produce good results but the inflation step lacks obvious geometric meanings. When using explicit iterative formulae to smooth complex models, one has to perform the iteration in fine steps in order to obtain satisfactory results. Desburn et al. [4] established an implicit Euler scheme based on the Euler–Lagrange equation to allow stepwise refinements and consequently improve stability. A method that preserves the mesh volume was developed to prevent the smoothed object from shrinkage in the same literature. Guskov et al. [9] also presented a mesh smoothing approach that used a non-uniform iterative procedure. As the weights of the iterative procedure depend on both geometry and connectivity, the method of Guskov et al. can produce smooth results without affecting too much the triangle shapes. However, further enhancement is needed for preserving detailed features. Considering that shrinkage and distortion of the Laplacian operator is mainly introduced by over-relaxation on vertices, Liu et al. [18] introduced constraints as penalty by keeping centroids of the mesh triangles fixed when fairing the mesh using the Laplacian operator. The approach is able to effectively keep the model from shrinkage but may lead to over-preservation of feature details. Recently, Liu et al. [19] proposed another constrained fairing algorithm by keep the local volume unchanged. The approach, however, results in an optimization problem with non-linear constraints [19]. In addition, the local volume of one region may be transferred to another area. This generally results in feature translation and serious distortion. Moreover, the method may suffer from over-smoothing. These drawbacks can be observed in Fig. 13. Ohtake et al. [21,22] also modified vertices of a noisy mesh in a way similar to Laplacian smoothing. However, a curvature-based function must be properly chosen to control the moving speed. To overcome the drawback of the approach, a userdefined threshold has to be introduced to locally control the geometric features. It is usually a difficult task for a user to determine such a threshold. First, it is difficult to recognize proper features visually. Second, even proper features have been identified, it is still inconvenient to set up such a relationship between the threshold and the identified features. As mentioned in Section 1, all the above approaches are designed for smoothing triangular meshes. The method based dual mesh resampling proposed in [31] admits arbitrary polygonal meshes as input and performs Laplacian smoothing on primal and dual spaces alternately. Anti-shrinkage is achieved by designing dual resampling matrices such that the invariant subspace associated with the largest eigenvalue of the multiplication of the matrices is maximized. The method is, however, only applied to closed models.
166
G. Li et al. / Graphical Models 66 (2004) 160–179
A unified smoothing approach, which addresses the smoothing of arbitrary polygonal meshes of two-dimensional manifold with triangular meshes as a special case, is described in the following sections.
3. A unified scheme for arbitrary mesh smoothing In this section we develop a unified scheme for arbitrary polygonal mesh smoothing. Our approach pffiffiffiis subdivided into three phases. The initial noisy mesh is first subdivided using a 3 split operator. An iterative smoothing procedure is then applied to the refined mesh via interpolation to the new vertices by employing a large stencil to compute the new position of old vertices. Finally, a Taubin-like inflation operator is designed to improve the anti-shrinkage ability of the approach. The first two phases are described in this section and the last phase will be discussed in Section 5. 3.1. Linear refinement Let M be the initial noisy mesh to be smoothed and V be its vertex set. To use the hypothesis that the face centroids of the noisy mesh arepexactly located on the anticffiffiffi ipated geometry, we subdivide the noise mesh using the 3 subdivision split operator described pffiffiffiin [15] and denote the resulted mesh as M ¼ ðV ; E; F ; P Þ. The 3 subdivision is originally designed to get richer resolution levels under a fixed maximal mesh size. The face number only increases by about two times in each refinement step, while pffiffiffi other schemes often raise the face number by three times. The split operation of 3 subdivision can be described as follows [15]. For an arbitrary polygonal mesh M, a 1-to-n split is performed for every face of degree n by inserting a new vertex at its center and then connecting the center to each vertex of the face. Next, edge flipping is applied to each pair of triangles shared by an old edge to form a new triangulation (see Figs. 1B–D). In the following, the vertices of M are categorized into two subsets: one consists of old vertices (belonging to V ) and the other new vertices (belonging to V V ). Obviously, the valence of each old vertex remains unchanged. We perform smoothing operations on M with the constraint of interpolation to the new vertices. It is expected that M tends to be smooth while M is faired repeatedly. pffiffiffi Note that it is not necessary to recomputed old vertices like the 3 subdivision for our purpose due to the linearity of refinement. Moreover, M is not induced by simply
Fig. 1. From left to right: (A) Original mesh; (B) insert face centers and connections to old vertices; (C) swap connection of center pairs of faces sharing an old edge; and (D) The final mesh after removing the old edges. Note that we did not treat boundary case in this example.
G. Li et al. / Graphical Models 66 (2004) 160–179
167
triangulating the face polygons of M as all edges of M have been removed in split operation. In fact, we will obtain a dual mesh of M if dropping all old vertices and associated edges and triangles in M. Generally, the dual one of a mesh can retain its rough shape. This is the reason why our approach can preserve the shape and detail features of the model to be smoothed. Here, the meaning of a dual mesh of a mesh is similar to that of the dual graph in graph theory. For a given mesh M with a simplicial complex K ¼ ðV ; E; F Þ, the simplicial complex KD ¼ ðVD ; ED ; FD Þ of its dual mesh MD ¼ ðKD ; PD Þ is defined as VD ¼ F , V ¼ FD and ðf1 ; f2 Þ 2 ED iff f1 and f2 shares an edge in M while for f 2 VD ; PD ðf Þ is defined as the center point of the face f in M. 3.2. A new interpolatory smoothing operator Because all adjacent vertices of an old vertex are new ones, ordinary smoothing algorithms (For example, the Laplacian smoothing), which compute the new position of a vertex by only using its adjacent vertices on M will fail to approximate the underlying geometry due to the fact that for each vertex in V , interpolating its all adjacent vertices in M (i.e., those new ones around it) will decouple the relation between adjacent vertices in M. Therefore in order to obtain a desirable result one should have to explore a larger stencil in M to compute the new position of old vertices. We choose the 1-neighborhood and the flaps of a vertex as such a stencil, as shown in Fig. 2. Then for each old vertex vi of M, we define the following variation associated with vi : rn X sn X vk þ vk vi n k2N ðiÞ n k2N ðiÞ M M rn X sn X ¼ ðrn sn Þvi þ cf þ vk ; n f 2F ðiÞ n k2N ðiÞ
Dvi ¼ ð1 rn sn Þvi þ
M
M
where n ¼ jNM ðiÞj is the valence of vi , cf is the centroid of face f , and rn , sn are two parameters to be determined later.
Fig. 2. The 1-neighborhood and the flaps of a vertex (marked by a solid circle).
168
G. Li et al. / Graphical Models 66 (2004) 160–179
Considering removing all edges between adjacent vertices of vi in Fig. 2, we derive a quadrilateral 1-neighborhood of vi as shown in Fig. 3. This is just a configuration used as the weight mask of the Catmull–Clark scheme [2]. Here, a simplest group of weights of this scheme [25] is used, i.e., we choose 3 rn ¼ ; 8
sn ¼
1 : 16
Noting that all face centers are fixed while smoothing, namely, for each face f of M we always have (denote its center after the lth iteration as clf , by the way, we will use superscript l to indicate the lth iteration): 1 X 0 clf ¼ v; jf j k2f k where v0k represents the original position of vertex fkg of the noisy mesh M, from the lth smoothing iteration formula vlþ1 ¼ vli þ kDvli , we have i "
vlþ1 i
# X 1 X X 7 3 1 ¼ vli þ k vli þ v0 þ vl : 16 8n f 2F ðiÞ jf j k2f k 16n k2N ðiÞ k M
ð7Þ
M
If M is a triangular mesh, especially, due to jf j ¼ 3 for all f Eq. (7) can be further simplified to # " k 1 X lþ1 l l 0 0 l vi ¼ vi þ ð4v þ vk Þ : ð8Þ 7vi þ 2vi þ 16 n k2N ðiÞ k M
Analogously, for quadrilateral meshes, Eq. (7) can be reduced to " # X X k 1 1 vlþ1 ¼ vli þ ð6v0 þ 2vlk Þ þ 3v0 ; 14vli þ 3v0i þ i 32 n k2N ðiÞ k n k2O ðiÞ k M
M
where OM ðiÞ denotes the set of all opposite vertices of fig.
Fig. 3. The quadrilateral face 1-neighborhood of the vertex extracted from Fig. 2.
ð9Þ
G. Li et al. / Graphical Models 66 (2004) 160–179
169
In fact, for a noisy mesh its face centers are usually not completely noise-free thus we run our approach multiple rounds with several iterative steps each time. The input of next round is the output of current round. The process can be executed repeatedly until terminative conditions are satisfied. Figs. 4–6 give examples to smooth
Fig. 4. A triangular Venus model. The first two images illustrate the original mesh. The second one is the noisy mesh. The rest are, respectively, smoothing results using our approach in 20, 40, and 60 iterative steps with k ¼ 0:5.
Fig. 5. A quadrilateral torus. The first two images are, respectively, original and noisy meshes. The third and the fourth are smoothing results of our method with 30 and 60 iteration steps, respectively, with a relaxation factor of 0.5.
Fig. 6. An arbitrary polygonal mesh model. The first two images are, respectively, the original and the noisy meshes. The rest images are smoothing results with 40, 80, and 120 iteration steps, respectively, using our method with a relaxation factor 0.5.
170
G. Li et al. / Graphical Models 66 (2004) 160–179
triangular, quadrilateral, and arbitrary polygonal models, respectively, by our method. 3.3. Boundary smoothing Few smoothing algorithms have carefully treated boundary noise although often claimed for being able to do so. We extend the technique for meshes discussed in Section 3.3 to smooth boundary curves. To guarantee that two patches can still stitch seamlessly, we employ a similar strategy used in subdivision schemes, namely, computing the new position of a boundary vertex with only its boundary proximate vertices. Finally, boundary polygons may be independently faired as noisy curves. To preserve features and prevent shrinkage while fairing a noisy polygon P , we also insert the midpoint of each edge of P into P to yield a new polygon P and perform a constrained smoothing operation on P , i.e., smoothing P with interpolation to the inserted vertices. Here, we also call the inserted vertices as new vertices of P and others old vertices. For a boundary vertex vi of P , by the same reason as selection of the computation stencil of Dv for an interior vertex v, we choose a five-point scheme (2-neighborhood) to define its variation as follows: a b Dvi ¼ ð1 a bÞvi þ ðvi;1 þ vi;1 Þ þ ðvi;2 þ vi;2 Þ vi ; 2 2 where vi;1 , vi;1 are two neighbors of vi on the boundary and vi;2 , vi;2 are, respectively, neighbors of vi;1 , vi;1 not equal to vi , as illustrated in Fig. 7. a and b are two parameters required to be determined. Notice that vi;1 , vi;1 are new vertices of P thus remain fixed while smoothing, namely, we always have vli;1 ¼ ðv0i þ v0i;2 Þ=2 and vli;1 ¼ ðv0i þ v0i;2 Þ=2. Therefore, we can derive b l a 0 lþ1 l l l l l 0 0 vi ¼ vi þ kDvi ¼ vi þ k ð a bÞvi þ ðvi;2 þ vi;2 Þ þ ðvi;2 þ 2vi þ vi;2 Þ : 2 4 Noting that for a boundary vertex vi , vi;2 and vi;2 are its two boundary adjacent vertices in the noisy mesh M, we can finally conclude a fairing iteration formula for vi " !# X b X l a lþ1 l l 0 0 2vi þ vi ¼ vi þ k ð a bÞvi þ v þ vk : ð10Þ 2 k2BN ðiÞ k 4 k2BN ðiÞ M
M
To determine a and b, we consider the subdivision approach for generating B-splines of degree 7 (B-splines of lower degree have no five-vertex stencil). Noting that weights used to update control vertices of B-splines of degree 7 are {1,28,70,28,1}/ 128, we put a ¼ 2 28=128 ¼ 7=16 and b ¼ 2 1=128 ¼ 1=64. Fig. 8 demonstrates an example for smoothing an open mesh with Eq. (10).
Fig. 7. A boundary vertex and its 2-neighborhood.
G. Li et al. / Graphical Models 66 (2004) 160–179
171
Fig. 8. Open mesh smoothing. The first two images are, respectively, original and noisy mashes. The rest are smoothing results using our method produced in 12 and 24 iteration steps, respectively.
4. Convergence discussion As our method has a very simple linear iteration form, it is not difficult to perform continuity analysis. Let us first consider Eq. (7). Without loss of generality, assume jMj ¼ m and denote vl ¼ ðvl1 ; vl2 ; . . . ; vlm ÞT , then from Eq. (7) we have vlþ1 ¼ Svl þ r;
ð11Þ
where S ¼ ðsij Þi;j¼1;...;m with sij defined by 8 i ¼ j; < 1 7k 16 k sij ¼ 16n j 2 NM ðiÞ; : 0 others T
and r ¼ ðr1 ; r2 ; . . . ; rm Þ with ri defined by ! 3k X 1 X 0 ri ¼ v 8n f 2F ðiÞ jf j k2f k M
is a constant vector. Therefore matrix S will be responsible for the convergence of Eq. (7). Note that all elements of S are positive for 0 < k < 16=7 and its row sum X 7k k 3k ¼1 : jsij j ¼ 1 þ 16 16 8 j¼1;...;m P So, for 0 < k < 1, we always have j¼1;...;m jsij j < 1. Hence, we can deduce the following conclusion. Property 1. For parameter 0 < k < 1, Eq. (7) is convergent. As special cases Eqs. (8) and (9) are also convergent.
5. Taubins inflation Though there is no shrinkage in iterative process with Eq. (7), regarding the face centers as points of the underlying surface generally introduces small shrinkage due to the fact of linear approximation. It is easy to understand that the coarser the noisy mesh, the more the shrinkage as the approximate error becomes large. This can be illustrated by the example in Fig. 9. In the figure, (A) is the original mesh while
172
G. Li et al. / Graphical Models 66 (2004) 160–179
Fig. 9. Example for inflation operation. From left to right: (A) The original mesh; (B) the noisy mesh; and (C and D) results created separately without and with inflation operation.
Fig. 10. Efficiency comparison between TaubinÕs algorithm and the proposed approach for smoothing the cow and the apple models. In the figures, Taubin Eq. (i) indicates that the weights defined by Eq. (1) are used for the kjl algorithm. T, the cost (seconds) for the iterations; N, the number of iterations.
(C) is the result using our method. To correct this defect, we perform TaubinÕs inflation operator one time to enlarge the mesh M after finishing the iterative process of Eq. (1) via looking at M as a triangle mesh. Namely, for all vertices of M, the following equation is executed k 1 X 1 X 0 vnew ¼ vold þ v; ð12Þ i i jPB k 1 n f 2F ðiÞ jf j k2f k M
where k=ðjPB k 1Þ is the relaxation factor of the inflation operation. Fig. 9D displays the smoothed mesh generated by performing our smoothing algorithm and the above inflation operation alternately with jPB ¼ 0:3. Figs. 11 and 12 show other two examples generated by combining Eqs. (7) and (12). The first example is an arbitrary polygonal mesh and the second a rectangular mesh. These examples also demonstrate the excellent capability of our approach to exactly preserve the model features. For instance, the ripples of the body and head of the original ant model are completely preserved after fairing as shown in Fig. 12.
6. Results and discussions In contrast to other smoothing algorithms, our approach has not only the advantage of being capable of fairing arbitrary polygonal meshes, but also strong ability in preserving features when used for smoothing triangular meshes. Here, we give an
G. Li et al. / Graphical Models 66 (2004) 160–179
173
Fig. 11. The Venus model. The first two images of the top row are the original model and the third one is the noisy mesh. The rest images are smoothing results with each derived by performing our method on the last smoothed mesh with 10 iteration steps. The relaxation factor used is 0.5.
example to compare our approach with Laplacian smoothing, Taubin smoothing, and volume-preserving smoothing of [19]. A cow model of 2094 vertices and 5804 faces is tested in the example. From Figs. 13E and F we can observe that the Laplacian smoothing suffers from serious shrinkage. With increasing of iterative times, cowÕs legs become more and more slender and almost disappear finally. On the other hand, LiuÕs method leads to over-smoothing with substantial distortion introduced by local volume transfer. Contrary to the Laplacian smoothing, clumsy cowÕs legs are generated using this approach as shown in Figs. 13G and H. For TaubinÕs kjl smoothing, three kinds of weights, respectively, defined by Eqs. (2), (4), and (5) are examined. It seems that the results of our approach have little difference from those of TaubinÕs algorithm using the first kind of weights (see Figs. 13I–N). However, when zooming to the model closely, we can clearly observe their difference as shown in Fig. 14. The results produced with our method more faithfully preserve the triangle shape of the original model than that with the kjl algorithm associated with weights of Eq. (2) that leads to over-smoothing (see Figs. 14A–C). The kjl algorithm of the last two types of weights could faithfully regenerate triangle shapes (see Figs. 14D and E), but results in instability in this example. This exhibits
174
G. Li et al. / Graphical Models 66 (2004) 160–179
Fig. 12. An Ant model. The first two images of the top row are, respectively, the original and the noisy meshes. The rest images are smoothing results with each derived by performing our method on the last smoothed mesh with 10 iteration steps.
Fig. 13. Comparison with other methods. For all approaches, 10 and 50 iteration steps are performed. (A and B) The original and the noisy meshes. (C and D) Smoothing results with the presented method; (E and F) Laplacian smoothing; (G and H) Liu et al.Õs method; (I and J), (K and L), and (M and N) are produced using TaubinÕs algorithm with weights of Eqs. (2), (4), and (5), respectively.
G. Li et al. / Graphical Models 66 (2004) 160–179
175
Fig. 14. Local zooming. (A) Original smooth model; (B–D) Locally zoomed views corresponding to Figs. 13D, J, L, and N. It demonstrates that the proposed method causes the least distortion compared to TaubinÕs algorithm.
Fig. 15. Smoothing of a scanned model. (A) The scanned data; (B) result generated using our method in 50 smoothing steps and 5 inflating steps (every 10 smoothing steps followed by 1-step inflation); (C) TaubinÕs algorithm in 100 iteration steps with weights of Eq. (2); and (D) the method of Liu et al. in 20 steps of iteration.
burrs in some regions for large numbers of iteration step. It is probably caused by the geometry-dependent weights that may behave ill-conditioned due to bad triangle shapes. Fig. 15 depicts an example for smoothing an apple model with 36,330 vertices and 72,664 faces. The mesh was extracted from data scanned using a CT-scanner. As a comparison, results produced using the kjl algorithm [30] and the volume-preserving approach [19] are also presented. Fig. 15B is generated after 50 steps of iteration and 5 steps of inflation using our method (10 steps of iteration followed by 1-step inflation). Figs. 15C and D illustrate the smoothing results using the approaches of [30] and [19], respectively. The volume-preserving approach of Liu et al. reaches the desired result with the least iteration steps (20 iterations). The proposed algorithm achieves the same visual effect in 50 steps. Though 100 steps have been performed, TaubinÕs algorithm still exhibits undulation on the apple surface. Fig. 16 illustrates another example for comparing the result of proposed approach with that of Goskov et al. [9]. Fig. 16C is generated using their software BEDA with maximal blur while Fig. 16B is the same as Fig. 15B. Employing geometry-independent weights, the presented approach is more efficient for each individual smoothing step than TaubinÕs kjl algorithm using the weights of Eqs. (4) and (5). In addition, the proposed approach uses one inflation for a group of iteration steps, but one inflation step should be added for each round of smoothing iteration for the latter approach in order to prevent shrinkage. It can therefore be conclude that TaubinÕs method is computationally more expensive when
176
G. Li et al. / Graphical Models 66 (2004) 160–179
Fig. 16. Comparison between the proposed method and Goskov et al.Õs algorithm: (A) the original mesh; (B) result generated with the proposed method; and (C) result generated with Goskov et al.Õs algorithm.
performing the same number of smoothing steps. Fig. 10 uses two examples to demonstrate this conclusion. On the other hand, among the three methods, the volume-preserving approach seems most computing-expensive. The computation cost is about 15 times more than that of the proposed approach for performing the same number of iterations. However, its rapidly convergent speed remedies the disadvantage to some extent. To evaluate the smoothing quality of the presented approach, we first define approximation errors. Suppose Mo be a smooth model and M be the mesh generated by adding white noise onto Mo while Ms be the smoothed mesh by denoising M. Considering that both the distance error and the shape error of triangles are desired between Mo and Ms , we linearly refine Mo and Ms one time, respectively, to obtain Mo1 ¼ ðVo1 ; Eo1 ; Fo1 ; Po1 Þ and Ms1 ¼ ðVs1 ; Es1 ; Fs1 ; Ps1 Þ. Noticing that both Mo1 and Ms1 have the same topology and their corresponding vertices have the same index, it is reasonable to define the following three types of errors between Mo and Ms : 1
Eave ðMo ; Ms Þ ¼
jVo j 1 1 X P ðiÞ P 1 ðiÞ; o s 1 jVo j i¼1
Emax ðMo ; Ms Þ ¼ maxPo1 ðiÞ Ps1 ðiÞ; i2Vo1
Emin ðMo ; Ms Þ ¼ minPo1 ðiÞ Ps1 ðiÞ: i2Vo1
Table 1 shows the results using these errors to test TaubinÕs kjl algorithm with three kinds of weights and the presented approach for the Cow model depicted in Fig. 13. In the table, TaubinÕs algorithm with different weights defined in Eqs. (2), (4), and (5), respectively, is examined while the proposed smoothing without and with two steps of TaubinÕs inflation is conducted. From the table we can observe that the presented approach with inflations achieves the best approximation errors for almost all cases in this example. For all the aforementioned examples, a PC with a Pentium 1.6 GHz CPU and 256 MB RAM running Windows XP operating system is used for testing and evaluation.
G. Li et al. / Graphical Models 66 (2004) 160–179
177
Table 1 Comparison of the approximation errors between the results of TaubinÕs algorithm and those of the present approach for the Cow model Iteration steps
25
Error types
ave.
max
min
ave.
max
min
ave.
max
min
Taubin Eq. (2) Taubin Eq. (4) Taubin Eq. (5) Our approach Our approach + 2 inflations
0.00633 0.00891 0.00893 0.00658 0.00640
0.04048 0.03597 0.02545 0.03214 0.02758
0.00027 0.00030 0.00049 0.00030 0.00022
0.00741 0.00926 0.00930 0.00658 0.00640
0.04409 0.04032 0.02847 0.03220 0.02765
0.00020 0.00062 0.00062 0.00030 0.00023
0.00877 0.00988 0.01030 0.00761 0.00714
0.04736 0.10305 0.04883 0.04111 0.03859
0.00030 0.00055 0.00046 0.00053 0.00044
50
100
7. Conclusions In literature, one can find many papers on mesh fairing. However, most of them are developed for smoothing triangular meshes. This paper presents a unified approach for fairing arbitrary two-manifold polygonal meshes. As the dual mesh of a mesh can roughly determined its geometric shape, our smoothing scheme is able to prevent shrinkage and preserve features by interpolating the vertices of the dual mesh. The scheme explores very simple weights, which are only dependent on mesh connectivity, for linear iteration and thus has a concise form. Moreover, we have carefully treated the boundary case in this paper. Mathematical discussion shows that our scheme is convergent for a large range of relaxation factor. As assuming that the face centers of the noisy mesh locate at the corresponding anticipated surface only leads to linear approximation to the original geometry and therefore may introduce shrinkage. To further eliminate this kind of shrinkage, a Taubin-like inflation operator is explored to enlarge the resulted mesh. Experiments have demonstrated the effectiveness and robustness of our approach. In TaubinÕs smoothing algorithm, the inflation operation is built as the dual operator of the Laplacian. We only used it to enlarge the intermediate meshes here and there is no meaning of signal processing. Therefore it is more difficult to determine the relaxation factor l so as to produce an appropriate scale rate. To find a more measurable inflation operator is our future work. In addition, it is also interesting to incorporate crease enhancement techniques to further preserve sharp features. At last, approaches based on signal processing and resampling are probably valuable for arbitrary polygonal meshes in further investigation.
Acknowledgments The work described in this paper is supported by Postdoctoral Foundation of China, Strategic Research Grants Project #7001074 from City University of Hong Kong, National Natural Science Foundations of China (Grant Nos. 60373034, 60021201, 60133020, and 69925204), and the 973 program of China (Grant No. 2002CB312102). We are grateful to Dr. Xinguo Liu for letting us use his program
178
G. Li et al. / Graphical Models 66 (2004) 160–179
of the volume-preserving algorithm. Models of the Venus, mannequin, and Beethoven bust came from http://www.mrl.nyu.edu/~biermann/subdivision/ and the ant mesh was downloaded from http://www.3dcafe.com/asp/ meshes.asp. References [1] D. Barash, A fundamental relationship between bilateral filtering, adaptive smoothing, and the nonlinear diffusion equation, IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (6) (2002) 844–847. [2] E. Catmull, J. Clark, Recursively generated B-spline surfaces on arbitrary topological meshes, Computer-Aided Design 10 (6) (1978) 350–355. [3] U. Clarenz, U. Diewald, M. Rumpf, Anisotropic geometric diffusion in surface processing, in: Proceedings of IEEE Visualization Õ00, 2000, pp. 397–405. [4] M. Desbrun, M. Meyer, P. Schr€ oder, A.H. Barr, Implicit fairing of irregular meshes using diffusion and curvature flow, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1999, pp. 317–324. [5] M. Desbrun, M. Meyer, P. Schr€ oder, A.H. Barr, Discrete differential-geometry operators in nD. Multi-res modeling group TR, Caltech, September 2000. Available from http://www.multires. caltech.edu/pubs/pubs.htm. [6] F. Durand, J. Dorsey, Fast bilateral filtering for the display of high-dynamic-range images, ACM Transactions on Graphics (SIGGRAPH 2001) 21 (3) (2002) 257–266. [7] S. Fleishman, I. Drori, D. Cohen-Or, Bilateral mesh denoising, ACM Transactions on Graphics (SIGGRAPH 2003) 22 (3) (2003). [8] K. Fujiwara, Eigenvalues of Laplacians on a closed Riemannian manifold and its nets, Proceedings of the AMS 123 (8) (1995) 2585–2594. [9] I. Guskov, W. Sweldens, P. Schr€ oder, Multiresolution signal processing, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1999, pp. 325–334. [10] M. Halstead, M. Kass, T. DeRose, Efficient, fair interpolation using Catmull–Clark surfaces, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1993, pp. 35–44. [11] H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, J. Schweitzer, W. Stuetzle, Mesh optimization, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1993, pp. 19–26. [12] T.R. Jones, F. Durand, M. Desbrun, Non-iterative, feature-preserving mesh smoothing, ACM Transactions on Graphics (SIGGRAPH 2003) 22 (3) (2003) 943–949. [13] L. Kobbelt, Discrete fairing, in: Proceedings of the Seventh IMA Conference on the Mathematics of Surface, 1997, pp. 101–131. [14] L. Kobbelt, J. Vorsatz, U. Labsik, H.-P. Seidel, A shrink wrapping approach to remeshing polygonal surfaces, Computer Graphics Forum (EUROGRAPHICS 1999) 18 (3) (1999) 119–130. pffiffiffi [15] L. Kobbelt, 3-Subdivision, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 2000, pp. 103–112. [16] L. Kobbelt, G. Taubin, Digital mesh processing. (SIGGRAPH Course Notes#17, electronic version). Publications Department, ACM Inc., Los Angeles, CA, 2001. [17] A. Khodakovsky, P. Alliez, M. Desbrun, P. Schr€ oder, Near-optimal connectivity encoding of 2manifold polygon, Graphical Models 64 (3/4) (2002) 147–168. [18] X. Liu, H. Bao, P.-A. Heng, T.-T. Wong, Q. Peng, Constrained mesh fairing, Computer Graphics Forum 20 (2) (2001) 115–124. [19] X. Liu, H. Bao, H.-Y. Shum, Q. Peng, A novel volume constrained smoothing method for meshes, Graphical Models 64 (3/4) (2002) 169–182. [20] H.P. Moreton, C.H. Sequin, Functional optimization for fair surface design, in: Computer Graphics Proceedings, ACM SIGGRAPH, 1992, pp. 167–175. [21] Y. Ohtake, A.G. Belyaev, I.A. Bogaevski, Polyhedral surface smoothing with simultaneous mesh regularization, in: Proceedings of the Geometric Modeling and Processing, Hong Kong, 2000, pp. 229–237.
G. Li et al. / Graphical Models 66 (2004) 160–179
179
[22] Y. Ohtake, A.G. Belyaev, I.A. Bogaevski, Mesh regularization and adaptive smoothing, ComputerAided Design 33 (11) (2001) 789–800. [23] Y. Ohtake, A.G. Belyaev, H.-P. Seidel, Mesh smoothing by adaptive and anisotropic Gaussian filter applied to mesh normals, in: G. Greiner, H. Niemann, T. Ertl, B. Girod, H.-P. Seidel (Eds.), Proceedings of Vision, Modeling, and Visualization 2002, Erlangen, Germany, 2002, pp. 203–210. [24] J. Peng, V. Strela, D. Zorin, A simple algorithm for surface denoising, in: IEEE Visualization Õ01, 2001, pp. 107–112. [25] H. Prautzsch, G. Umlauf, A G2 -subdivision algorithm, in: G. Farin, H. Bieri, G. Brunnet, T. DeRose (Eds.), Geometric Modeling, Computing Suppl., vol. 13, Springer-Verlag, Berlin, 1998, pp. 217–224. [26] R. Schneider, L. Kobbelt, Generating fair meshes with G1 boundary conditions, in: Proceedings of Geometric Modeling and Processing, Hong Kong, 2000, pp. 251–261. [27] R. Schneider, L. Kobbelt, Mesh fairing based on an intrinsic PDE approach, CAD 33 (11) (2001) 767– 777. [28] R. Schneider, L. Kobbelt, Geometric fairing of irregular meshes for free-form surface design, Computer Aided Geometry Design 18 (4) (2001) 359–379. [29] T. Tasdizen, R. Whitaker, P. Burchard, S. Osher, Geometric surface smoothing via anisotropic diffusion of normals, in: Procceedings of IEEE Visualization Õ02, 2002, pp. 125–132. [30] G. Taubin, A signal processing approach to fair surface design, in: Computer Graphics Proceedings, ACM SIGGRAPH, 1995, pp. 351–358. [31] G. Taubin, Dual Mesh Resampling, Pacific Graphics, Tokyo, Japan, 2001, pp. 180–189. [32] G. Taubin, Linear anisotropic mesh filtering. IBM Research Report RC-22213, 2001. [33] C. Tomasi, R. Manduchi, Bilateral filtering for gray and color image, in: Proceedings of the 1998 IEEE International Conference on Computer Vision, 1998. [34] J. Vollmer, R. Mencl, H. Muller, Improved Laplacian smoothing of noisy surface meshes, Computer Graphics Forum (EUROGRAPHICS) 18 (3) (1999) 131–138. [35] W. Welch, A. Witkin, Variational surface modeling, in: Computer Graphics Proceedings, ACM SIGGRAPH, 1992, pp. 157–166. [36] W. Welch, A. Witkin, Free-form shape design using triangulated surfaces, in: Computer Graphics Proceedings, ACM SIGGRAPH, 1994, pp. 247–256. [37] G. Westgaard, H. Nowacki, A process for surface fairing in irregular meshes, Computer Aided Geometry Design 18 (7) (2001) 619–638. [38] Z. Wood, M. Desbrun, P. Schr€ oder, D. Breen, Semi-regular mesh extraction from volumes, in: Proceedings of IEEE Visualization, 2000, pp. 275–282. [39] H. Zhang, E. Fiume, Mesh smoothing with shape or feature preservation, in: J. Vince, R. Earnshaw (Eds.), Advances in Modeling, Animation, and Rendering, 2002, pp. 167–182.
A unified approach for fairing arbitrary polygonal meshes Guiqing Li,a,b,c,* Hujun Bao,a and Weiyin Mab a
The State Key Lab of CAD&CG, Zhejiang University, Hangzhou, China Department of Manufacturing Engineering and Engineering Management, City University of Hong Kong, Kowloon, Hong Kong, China c School of Computer Science and Engineering, South China University of Technology, Guangzhou, China
b
Received 31 May 2002; received in revised form 17 January 2004; accepted 5 March 2004 Available online 6 May 2004
Abstract Arbitrary polygonal meshes play an important role in 3D model representation and surface modeling. Most of the previous smoothing approaches are, however, mainly designed for triangular meshes. This paper proposes a unified approach to address fairing problems of arbitrary polygonal meshes. The approach is built upon an assumption that all face centroids of a smoothed mesh locate on its noise-free surface and remain unchanged during the smoothing process in order to prevent the resulted mesh from shrinkage and preserve surface features. Boundary smoothing is considered in a similar way as well but the process is independent to mesh processing such that two noisy meshes sharing a boundary still adjoin continuously after smoothing. Finally, a Taubin-like inflation operator is employed to further improve the ability of anti-shrinkage. Although very simple iterative equation is employed in the smoothing procedure, experiments demonstrate that the approach is able to not only effectively address the problem of shrinkage and distortion but also preserve features very well. 2004 Elsevier Inc. All rights reserved. Keywords: Arbitrary mesh smoothing; Laplacian operator; Mesh modeling; Mesh split
*
Corresponding author. E-mail addresses: [email protected] (G. Li), [email protected] (H. Bao), [email protected] (W. Ma). 1524-0703/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.gmod.2004.03.001
G. Li et al. / Graphical Models 66 (2004) 160–179
161
1. Introduction With advancement of 3D data acquisition techniques such as laser scanners, 3D photography, and CT-scanning, very large meshes are commonly used as surface representation for graphics applications today for their simplicity and flexibility. However, digitization processes often introduce noise into mesh attributes such as vertex position due to some physical factors associated with measurement, thus removing noise is required before performing any other operations on the model. On the other hand, to yield free-form surfaces for fairness, it is also necessary to smooth control nets automatically. The fairness of a surface can usually be measured by certain energy functional of geometric attributes such as curvature and curvature variation. The essence of surface smoothing is therefore directly or indirectly to minimize these energies [10,11,20,35,36]. Noticing geometric curvature quantities can be approximated by combination of second derivatives under the assumption that a parametric surface is rather flat [13], some authors employed all kinds of lower or higher derivatives to replace geometric intrinsics for fairing surfaces. Minimizing the membrane or thin plate energy, for example, can lead to surfaces of excellent quality [13,16,26,37]. However, approaches directly involving functional minimization often exhibit low efficiency due to expensive computational costs. Fortunately, as a variational problem can usually be transformed into PDEs, relaxation iterative techniques can often be employed to build an iterative process for finding numerical solutions of the PDEs in which each iterative step moves vertex positions by a small increment. Approaches based on the Laplacian operator fall into this category. Taubin [30] pioneered such a kind of algorithms by reducing the problem of mesh smoothing to low-pass filtering. Actually, most of the previous smoothing techniques fall into this framework [4,9,14,18,19,21,34,38,39]. Some of them will be further discussed in Section 2. To produce surfaces of high aesthetic quality, Schneider and co-worker [27,28] proposed a completely intrinsic approach based on PDEs that no longer comes from functional minimization. The approach is able to reproduce surfaces with constant mean curvature such as spheres, cylinders, and minimal surfaces. Based on the linear heat equation, Clarenz et al. [3] formulated an anisotropic diffusion approach for enhancing geometric features such as creases and corners of meshes smoothed. Bajaj and Xu generalized the above method to a two-dimensional manifold embedded in Rk so as to put vertex positions of the meshes and functions defined on the meshes into a unified smoothing scheme in which the local analytic form of the Loop subdivision was employed to represent the mesh surfaces and associated vector functions. Statistics-based algorithms were pioneered by Peng et al. [24]. Their approach is limited to smooth semi-regular meshes and pre-processing is needed to convert the meshes to multiresolution representations before fairing. The Weiner filter is then applied to smooth the directional detail coefficients. Denoised results are finally reconstructed from the smoothed details. Recent statistics-based approaches were motivated by the successful application of bilateral filtering to image smoothing [1,6,33]. Based on robust statistics and local first-order predictors of surfaces, Jones
162
G. Li et al. / Graphical Models 66 (2004) 160–179
et al. [12] proposed a non-iterative feature-preserving approach. With this method, face normals should be evaluated first since prediction is made via projecting the estimated vertex onto mesh faces. Bilateral mesh denoising proposed by Fleishman et al. [7] also employed the same technique but provided a simpler predictor. They evaluated the offset from estimated vertex to its underlying position, i.e., the predictor, by casting its neighbors onto its normal line. For large-scale noise, it is not easy to obtain good estimation of the underlying normals even if the surface normals are first smoothed as in [12]. In addition, Ohtake et al. [23] smoothed face normals of the noisy mesh using Gaussian filter. The result was then obtained by fitting the faired normals. By introducing anisotropic neighborhoods, their approach was adapted to smooth models with sharp feature. Though the method was argued fully automatic, the standard deviation r of Gaussian-like functions and the positive constant k of the modified distance have to be determined through experiment for smoothing face normals. The method of Tasdizen et al. [29] bore the same framework, but the anisotropic diffusion of normals was realized by minimizing the total curvature of normal maps instead of statistics-based. However, all the above-mentioned approaches aim at addressing smoothness of triangular mesh models while in practice arbitrary polygonal meshes are often encountered in some applications [17] such as surface design based on control nets. One may convert a polygonal mesh into a triangular mesh and then apply smoothing algorithms developed for triangle meshes. It is, however, not appropriate to do so and it is better to directly remove the noise before any processing similar to the case for processing triangle meshes. This paper addresses the smoothing of arbitrary polygonal meshes. To prevent shrinkage and preserve detailed features, the smoothing algorithm is designed such that all face centers pffiffiffi remain unchanged after smoothing. The noisy model is first subdivided with a 3 split operator by inserting face centers as new vertices. Then the refined mesh is smoothed under the constraint that the newly inserted vertices remain unchanged. The denoised mesh is finally obtained by replacing vertex positions of the noisy model with those corresponding vertices of the smoothed refined mesh. As the approach can deal with arbitrary polygonal meshes, triangular meshes can also be processed under this unified framework. Main contributions of this paper can be summarized as follows: • Arbitrary polygonal mesh smoothing: An approach is proposed to smooth arbitrary polygonal meshes. In the algorithm, a very simple linear combination based on Catmull–Clark subdivision scheme is explored to move vertices of the noisy mesh in the iterative process. • Boundary smoothing: A constrained curve smoothing approach is developed by interpolating midpoints of boundary edges such that the boundary can be smoothed simultaneously. • Convergence treatment: A formal analysis on convergence of the algorithm is discussed. • Inflation operation: A Taubin-like operator is designed to further enhance the anti-shrinkage effect.
G. Li et al. / Graphical Models 66 (2004) 160–179
163
The rest of the paper is organized as follows. In Section 2, some closely related work is reviewed. Our smoothing approach is subsequently described in Section 3 in further details. Some mathematical treatment is discussed in Section 4. To further improve the effect of anti-shrinkage, we introduce an inflation operator similar to that of [30] in Section 5. Comparison between the presented approach and two other algorithms is discussed in Section 6. Finally, conclusions are drawn in the last section.
2. Notations and related work In this section, we first formulate some notations and review several related approaches on mesh smoothing. 2.1. Formal conventions An arbitrary polygonal mesh with vertex set V ¼ f1; 2; . . . ; ng is represented by M ¼ ðK; P Þ, where K 2V (2V expresses the set consisting of all subset of V ) is a simplicial complex consisting of subsets of V and P is a mapping from V to coordinates of its vertices, i.e., P : V ! R3 . Elements of K containing only one integer, two integers are called vertices and edges of M, respectively. All other elements are defined as faces of M. For instance, fig indicates vertex with index i while fi; jg refers to edge with endpoints fig and fjg, and fi; j; kg represents a triangle with vertices fig, fjg, and fkg. Generally, small letters such as e or f are used to represent edges or faces. We will also use vi to indicate P ðiÞ, the coordinates of fig in the rest of this paper. If fi; jg 2 E, then fig and fjg are said adjacent or neighboring. Set consisting of all adjacent vertices of fig is called its 1-neighborhood in associated with mesh M denoted by NM ðiÞ while the number of its neighbors is called valence and read as jNM ðiÞj. Accordingly, BNM ðiÞ is used to express the set of the boundary adjacent vertices of a boundary vertex fig. The set of all faces sharing vertex fig is defined as the face 1-neighborhood of fig and denoted as FM ðiÞ. In addition, the degree of a face is defined as its vertex number. We also assume that meshes are 2D manifold. No more than two faces share an edge, each vertex belongs to at least an edge and each edge belongs to at least a face. 2.2. Related work Let f ðx; yÞ be a differentiable function and D be the Laplacian operator Df ¼ fxx þ fyy . A variational problem of membrane energy can be described as follows Z min ðfx2 þ fy2 Þ dA; f
where dA is an area element on the surface z ¼ f ðx; yÞ. From variational theory we know that the differential form of the above variational problem is the Euler– Lagrange equation
164
G. Li et al. / Graphical Models 66 (2004) 160–179
Df ¼ 0: For any vertex vi of an irregular mesh, the simplest discrete form of the Laplacian operator can be expressed as X Dvi ¼ vi þ wij vj ; ð1Þ j2NM ðiÞ
where wij ¼
1 : jNM ðiÞj
ð2Þ
One can then use the following relaxation formula to iteratively find the approximate solution of the discrete Euler–Lagrange equation Dvi ¼ 0: ¼ vold þ kDvi ; vnew i i
ð3Þ
where k is a positive relaxation factor. We refer the smoothing method based on Eq. (3) as Laplacian smoothing. Notice that Eq. (2) defines uniform weights for all neighboring vertices. In literature, one may find an approach that uses the reverse of the edge length to replace the weights of Eq. (2) [4,8]: X 1 1 wij ¼ : ð4Þ jvi vj j j2N ðiÞ jvi vj j M
One may also find a sophisticated set of weights that is related to the surface curvature, an intrinsic attribute of the mesh geometry [4,5,32]: wij ¼ P
ðcot aij þ cot bij Þ ; j2NM ðiÞ ðcot aij þ cot bij Þ
ð5Þ
where aij and bij are the two angles opposite to edge ði; jÞ. Many variants of Laplacian smoothing were proposed due to the fact that the Laplacian smoothing may introduce serious shrinkage and distortion and result in over-smoothing. Some of the approaches involve Laplacian operator of higher order [30] and some others use constraint-based approaches [4,18,19,34]. Based on techniques for signal processing and eigenvalue analysis of the matrix derived from the Laplacian smoothing formula, Taubin found that the shrinkage problem could be eliminated by iteratively computing Eq. (3) following an inflation step vnew ¼ vold þ lDvi ; i i
ð6Þ
where l is a negative relaxation factor satisfying l > k. However, when vi is located on the boundary, boundary shrinkage cannot be avoided due to a very strong tangential component of Dvi . Though successful in preventing shrinkage, TaubinÕs scheme runs into difficulty in preserving features [21,22]. Taubin [32] addressed boundary shrinkage by replacing Dvi with its normal component in which the normal can be either estimated by averaging the normals of adjacent faces or specified as constraints by users. Though it is able to maintain the same convergent speed for
G. Li et al. / Graphical Models 66 (2004) 160–179
165
internal and boundary vertices using the same iterative framework, the approach may produce cracks on boundaries shared by two meshes. Another variant was proposed in [34]. The method also includes two steps—a variation of the Laplacian and a so-called HC-modification. The first step smoothes the mesh by moving all vertices of the old mesh based on the Laplacian smoothing criteria and the second step pushes the vertices of the resulting mesh of the first step back towards their original positions by certain distance. The approach is able to produce good results but the inflation step lacks obvious geometric meanings. When using explicit iterative formulae to smooth complex models, one has to perform the iteration in fine steps in order to obtain satisfactory results. Desburn et al. [4] established an implicit Euler scheme based on the Euler–Lagrange equation to allow stepwise refinements and consequently improve stability. A method that preserves the mesh volume was developed to prevent the smoothed object from shrinkage in the same literature. Guskov et al. [9] also presented a mesh smoothing approach that used a non-uniform iterative procedure. As the weights of the iterative procedure depend on both geometry and connectivity, the method of Guskov et al. can produce smooth results without affecting too much the triangle shapes. However, further enhancement is needed for preserving detailed features. Considering that shrinkage and distortion of the Laplacian operator is mainly introduced by over-relaxation on vertices, Liu et al. [18] introduced constraints as penalty by keeping centroids of the mesh triangles fixed when fairing the mesh using the Laplacian operator. The approach is able to effectively keep the model from shrinkage but may lead to over-preservation of feature details. Recently, Liu et al. [19] proposed another constrained fairing algorithm by keep the local volume unchanged. The approach, however, results in an optimization problem with non-linear constraints [19]. In addition, the local volume of one region may be transferred to another area. This generally results in feature translation and serious distortion. Moreover, the method may suffer from over-smoothing. These drawbacks can be observed in Fig. 13. Ohtake et al. [21,22] also modified vertices of a noisy mesh in a way similar to Laplacian smoothing. However, a curvature-based function must be properly chosen to control the moving speed. To overcome the drawback of the approach, a userdefined threshold has to be introduced to locally control the geometric features. It is usually a difficult task for a user to determine such a threshold. First, it is difficult to recognize proper features visually. Second, even proper features have been identified, it is still inconvenient to set up such a relationship between the threshold and the identified features. As mentioned in Section 1, all the above approaches are designed for smoothing triangular meshes. The method based dual mesh resampling proposed in [31] admits arbitrary polygonal meshes as input and performs Laplacian smoothing on primal and dual spaces alternately. Anti-shrinkage is achieved by designing dual resampling matrices such that the invariant subspace associated with the largest eigenvalue of the multiplication of the matrices is maximized. The method is, however, only applied to closed models.
166
G. Li et al. / Graphical Models 66 (2004) 160–179
A unified smoothing approach, which addresses the smoothing of arbitrary polygonal meshes of two-dimensional manifold with triangular meshes as a special case, is described in the following sections.
3. A unified scheme for arbitrary mesh smoothing In this section we develop a unified scheme for arbitrary polygonal mesh smoothing. Our approach pffiffiffiis subdivided into three phases. The initial noisy mesh is first subdivided using a 3 split operator. An iterative smoothing procedure is then applied to the refined mesh via interpolation to the new vertices by employing a large stencil to compute the new position of old vertices. Finally, a Taubin-like inflation operator is designed to improve the anti-shrinkage ability of the approach. The first two phases are described in this section and the last phase will be discussed in Section 5. 3.1. Linear refinement Let M be the initial noisy mesh to be smoothed and V be its vertex set. To use the hypothesis that the face centroids of the noisy mesh arepexactly located on the anticffiffiffi ipated geometry, we subdivide the noise mesh using the 3 subdivision split operator described pffiffiffiin [15] and denote the resulted mesh as M ¼ ðV ; E; F ; P Þ. The 3 subdivision is originally designed to get richer resolution levels under a fixed maximal mesh size. The face number only increases by about two times in each refinement step, while pffiffiffi other schemes often raise the face number by three times. The split operation of 3 subdivision can be described as follows [15]. For an arbitrary polygonal mesh M, a 1-to-n split is performed for every face of degree n by inserting a new vertex at its center and then connecting the center to each vertex of the face. Next, edge flipping is applied to each pair of triangles shared by an old edge to form a new triangulation (see Figs. 1B–D). In the following, the vertices of M are categorized into two subsets: one consists of old vertices (belonging to V ) and the other new vertices (belonging to V V ). Obviously, the valence of each old vertex remains unchanged. We perform smoothing operations on M with the constraint of interpolation to the new vertices. It is expected that M tends to be smooth while M is faired repeatedly. pffiffiffi Note that it is not necessary to recomputed old vertices like the 3 subdivision for our purpose due to the linearity of refinement. Moreover, M is not induced by simply
Fig. 1. From left to right: (A) Original mesh; (B) insert face centers and connections to old vertices; (C) swap connection of center pairs of faces sharing an old edge; and (D) The final mesh after removing the old edges. Note that we did not treat boundary case in this example.
G. Li et al. / Graphical Models 66 (2004) 160–179
167
triangulating the face polygons of M as all edges of M have been removed in split operation. In fact, we will obtain a dual mesh of M if dropping all old vertices and associated edges and triangles in M. Generally, the dual one of a mesh can retain its rough shape. This is the reason why our approach can preserve the shape and detail features of the model to be smoothed. Here, the meaning of a dual mesh of a mesh is similar to that of the dual graph in graph theory. For a given mesh M with a simplicial complex K ¼ ðV ; E; F Þ, the simplicial complex KD ¼ ðVD ; ED ; FD Þ of its dual mesh MD ¼ ðKD ; PD Þ is defined as VD ¼ F , V ¼ FD and ðf1 ; f2 Þ 2 ED iff f1 and f2 shares an edge in M while for f 2 VD ; PD ðf Þ is defined as the center point of the face f in M. 3.2. A new interpolatory smoothing operator Because all adjacent vertices of an old vertex are new ones, ordinary smoothing algorithms (For example, the Laplacian smoothing), which compute the new position of a vertex by only using its adjacent vertices on M will fail to approximate the underlying geometry due to the fact that for each vertex in V , interpolating its all adjacent vertices in M (i.e., those new ones around it) will decouple the relation between adjacent vertices in M. Therefore in order to obtain a desirable result one should have to explore a larger stencil in M to compute the new position of old vertices. We choose the 1-neighborhood and the flaps of a vertex as such a stencil, as shown in Fig. 2. Then for each old vertex vi of M, we define the following variation associated with vi : rn X sn X vk þ vk vi n k2N ðiÞ n k2N ðiÞ M M rn X sn X ¼ ðrn sn Þvi þ cf þ vk ; n f 2F ðiÞ n k2N ðiÞ
Dvi ¼ ð1 rn sn Þvi þ
M
M
where n ¼ jNM ðiÞj is the valence of vi , cf is the centroid of face f , and rn , sn are two parameters to be determined later.
Fig. 2. The 1-neighborhood and the flaps of a vertex (marked by a solid circle).
168
G. Li et al. / Graphical Models 66 (2004) 160–179
Considering removing all edges between adjacent vertices of vi in Fig. 2, we derive a quadrilateral 1-neighborhood of vi as shown in Fig. 3. This is just a configuration used as the weight mask of the Catmull–Clark scheme [2]. Here, a simplest group of weights of this scheme [25] is used, i.e., we choose 3 rn ¼ ; 8
sn ¼
1 : 16
Noting that all face centers are fixed while smoothing, namely, for each face f of M we always have (denote its center after the lth iteration as clf , by the way, we will use superscript l to indicate the lth iteration): 1 X 0 clf ¼ v; jf j k2f k where v0k represents the original position of vertex fkg of the noisy mesh M, from the lth smoothing iteration formula vlþ1 ¼ vli þ kDvli , we have i "
vlþ1 i
# X 1 X X 7 3 1 ¼ vli þ k vli þ v0 þ vl : 16 8n f 2F ðiÞ jf j k2f k 16n k2N ðiÞ k M
ð7Þ
M
If M is a triangular mesh, especially, due to jf j ¼ 3 for all f Eq. (7) can be further simplified to # " k 1 X lþ1 l l 0 0 l vi ¼ vi þ ð4v þ vk Þ : ð8Þ 7vi þ 2vi þ 16 n k2N ðiÞ k M
Analogously, for quadrilateral meshes, Eq. (7) can be reduced to " # X X k 1 1 vlþ1 ¼ vli þ ð6v0 þ 2vlk Þ þ 3v0 ; 14vli þ 3v0i þ i 32 n k2N ðiÞ k n k2O ðiÞ k M
M
where OM ðiÞ denotes the set of all opposite vertices of fig.
Fig. 3. The quadrilateral face 1-neighborhood of the vertex extracted from Fig. 2.
ð9Þ
G. Li et al. / Graphical Models 66 (2004) 160–179
169
In fact, for a noisy mesh its face centers are usually not completely noise-free thus we run our approach multiple rounds with several iterative steps each time. The input of next round is the output of current round. The process can be executed repeatedly until terminative conditions are satisfied. Figs. 4–6 give examples to smooth
Fig. 4. A triangular Venus model. The first two images illustrate the original mesh. The second one is the noisy mesh. The rest are, respectively, smoothing results using our approach in 20, 40, and 60 iterative steps with k ¼ 0:5.
Fig. 5. A quadrilateral torus. The first two images are, respectively, original and noisy meshes. The third and the fourth are smoothing results of our method with 30 and 60 iteration steps, respectively, with a relaxation factor of 0.5.
Fig. 6. An arbitrary polygonal mesh model. The first two images are, respectively, the original and the noisy meshes. The rest images are smoothing results with 40, 80, and 120 iteration steps, respectively, using our method with a relaxation factor 0.5.
170
G. Li et al. / Graphical Models 66 (2004) 160–179
triangular, quadrilateral, and arbitrary polygonal models, respectively, by our method. 3.3. Boundary smoothing Few smoothing algorithms have carefully treated boundary noise although often claimed for being able to do so. We extend the technique for meshes discussed in Section 3.3 to smooth boundary curves. To guarantee that two patches can still stitch seamlessly, we employ a similar strategy used in subdivision schemes, namely, computing the new position of a boundary vertex with only its boundary proximate vertices. Finally, boundary polygons may be independently faired as noisy curves. To preserve features and prevent shrinkage while fairing a noisy polygon P , we also insert the midpoint of each edge of P into P to yield a new polygon P and perform a constrained smoothing operation on P , i.e., smoothing P with interpolation to the inserted vertices. Here, we also call the inserted vertices as new vertices of P and others old vertices. For a boundary vertex vi of P , by the same reason as selection of the computation stencil of Dv for an interior vertex v, we choose a five-point scheme (2-neighborhood) to define its variation as follows: a b Dvi ¼ ð1 a bÞvi þ ðvi;1 þ vi;1 Þ þ ðvi;2 þ vi;2 Þ vi ; 2 2 where vi;1 , vi;1 are two neighbors of vi on the boundary and vi;2 , vi;2 are, respectively, neighbors of vi;1 , vi;1 not equal to vi , as illustrated in Fig. 7. a and b are two parameters required to be determined. Notice that vi;1 , vi;1 are new vertices of P thus remain fixed while smoothing, namely, we always have vli;1 ¼ ðv0i þ v0i;2 Þ=2 and vli;1 ¼ ðv0i þ v0i;2 Þ=2. Therefore, we can derive b l a 0 lþ1 l l l l l 0 0 vi ¼ vi þ kDvi ¼ vi þ k ð a bÞvi þ ðvi;2 þ vi;2 Þ þ ðvi;2 þ 2vi þ vi;2 Þ : 2 4 Noting that for a boundary vertex vi , vi;2 and vi;2 are its two boundary adjacent vertices in the noisy mesh M, we can finally conclude a fairing iteration formula for vi " !# X b X l a lþ1 l l 0 0 2vi þ vi ¼ vi þ k ð a bÞvi þ v þ vk : ð10Þ 2 k2BN ðiÞ k 4 k2BN ðiÞ M
M
To determine a and b, we consider the subdivision approach for generating B-splines of degree 7 (B-splines of lower degree have no five-vertex stencil). Noting that weights used to update control vertices of B-splines of degree 7 are {1,28,70,28,1}/ 128, we put a ¼ 2 28=128 ¼ 7=16 and b ¼ 2 1=128 ¼ 1=64. Fig. 8 demonstrates an example for smoothing an open mesh with Eq. (10).
Fig. 7. A boundary vertex and its 2-neighborhood.
G. Li et al. / Graphical Models 66 (2004) 160–179
171
Fig. 8. Open mesh smoothing. The first two images are, respectively, original and noisy mashes. The rest are smoothing results using our method produced in 12 and 24 iteration steps, respectively.
4. Convergence discussion As our method has a very simple linear iteration form, it is not difficult to perform continuity analysis. Let us first consider Eq. (7). Without loss of generality, assume jMj ¼ m and denote vl ¼ ðvl1 ; vl2 ; . . . ; vlm ÞT , then from Eq. (7) we have vlþ1 ¼ Svl þ r;
ð11Þ
where S ¼ ðsij Þi;j¼1;...;m with sij defined by 8 i ¼ j; < 1 7k 16 k sij ¼ 16n j 2 NM ðiÞ; : 0 others T
and r ¼ ðr1 ; r2 ; . . . ; rm Þ with ri defined by ! 3k X 1 X 0 ri ¼ v 8n f 2F ðiÞ jf j k2f k M
is a constant vector. Therefore matrix S will be responsible for the convergence of Eq. (7). Note that all elements of S are positive for 0 < k < 16=7 and its row sum X 7k k 3k ¼1 : jsij j ¼ 1 þ 16 16 8 j¼1;...;m P So, for 0 < k < 1, we always have j¼1;...;m jsij j < 1. Hence, we can deduce the following conclusion. Property 1. For parameter 0 < k < 1, Eq. (7) is convergent. As special cases Eqs. (8) and (9) are also convergent.
5. Taubins inflation Though there is no shrinkage in iterative process with Eq. (7), regarding the face centers as points of the underlying surface generally introduces small shrinkage due to the fact of linear approximation. It is easy to understand that the coarser the noisy mesh, the more the shrinkage as the approximate error becomes large. This can be illustrated by the example in Fig. 9. In the figure, (A) is the original mesh while
172
G. Li et al. / Graphical Models 66 (2004) 160–179
Fig. 9. Example for inflation operation. From left to right: (A) The original mesh; (B) the noisy mesh; and (C and D) results created separately without and with inflation operation.
Fig. 10. Efficiency comparison between TaubinÕs algorithm and the proposed approach for smoothing the cow and the apple models. In the figures, Taubin Eq. (i) indicates that the weights defined by Eq. (1) are used for the kjl algorithm. T, the cost (seconds) for the iterations; N, the number of iterations.
(C) is the result using our method. To correct this defect, we perform TaubinÕs inflation operator one time to enlarge the mesh M after finishing the iterative process of Eq. (1) via looking at M as a triangle mesh. Namely, for all vertices of M, the following equation is executed k 1 X 1 X 0 vnew ¼ vold þ v; ð12Þ i i jPB k 1 n f 2F ðiÞ jf j k2f k M
where k=ðjPB k 1Þ is the relaxation factor of the inflation operation. Fig. 9D displays the smoothed mesh generated by performing our smoothing algorithm and the above inflation operation alternately with jPB ¼ 0:3. Figs. 11 and 12 show other two examples generated by combining Eqs. (7) and (12). The first example is an arbitrary polygonal mesh and the second a rectangular mesh. These examples also demonstrate the excellent capability of our approach to exactly preserve the model features. For instance, the ripples of the body and head of the original ant model are completely preserved after fairing as shown in Fig. 12.
6. Results and discussions In contrast to other smoothing algorithms, our approach has not only the advantage of being capable of fairing arbitrary polygonal meshes, but also strong ability in preserving features when used for smoothing triangular meshes. Here, we give an
G. Li et al. / Graphical Models 66 (2004) 160–179
173
Fig. 11. The Venus model. The first two images of the top row are the original model and the third one is the noisy mesh. The rest images are smoothing results with each derived by performing our method on the last smoothed mesh with 10 iteration steps. The relaxation factor used is 0.5.
example to compare our approach with Laplacian smoothing, Taubin smoothing, and volume-preserving smoothing of [19]. A cow model of 2094 vertices and 5804 faces is tested in the example. From Figs. 13E and F we can observe that the Laplacian smoothing suffers from serious shrinkage. With increasing of iterative times, cowÕs legs become more and more slender and almost disappear finally. On the other hand, LiuÕs method leads to over-smoothing with substantial distortion introduced by local volume transfer. Contrary to the Laplacian smoothing, clumsy cowÕs legs are generated using this approach as shown in Figs. 13G and H. For TaubinÕs kjl smoothing, three kinds of weights, respectively, defined by Eqs. (2), (4), and (5) are examined. It seems that the results of our approach have little difference from those of TaubinÕs algorithm using the first kind of weights (see Figs. 13I–N). However, when zooming to the model closely, we can clearly observe their difference as shown in Fig. 14. The results produced with our method more faithfully preserve the triangle shape of the original model than that with the kjl algorithm associated with weights of Eq. (2) that leads to over-smoothing (see Figs. 14A–C). The kjl algorithm of the last two types of weights could faithfully regenerate triangle shapes (see Figs. 14D and E), but results in instability in this example. This exhibits
174
G. Li et al. / Graphical Models 66 (2004) 160–179
Fig. 12. An Ant model. The first two images of the top row are, respectively, the original and the noisy meshes. The rest images are smoothing results with each derived by performing our method on the last smoothed mesh with 10 iteration steps.
Fig. 13. Comparison with other methods. For all approaches, 10 and 50 iteration steps are performed. (A and B) The original and the noisy meshes. (C and D) Smoothing results with the presented method; (E and F) Laplacian smoothing; (G and H) Liu et al.Õs method; (I and J), (K and L), and (M and N) are produced using TaubinÕs algorithm with weights of Eqs. (2), (4), and (5), respectively.
G. Li et al. / Graphical Models 66 (2004) 160–179
175
Fig. 14. Local zooming. (A) Original smooth model; (B–D) Locally zoomed views corresponding to Figs. 13D, J, L, and N. It demonstrates that the proposed method causes the least distortion compared to TaubinÕs algorithm.
Fig. 15. Smoothing of a scanned model. (A) The scanned data; (B) result generated using our method in 50 smoothing steps and 5 inflating steps (every 10 smoothing steps followed by 1-step inflation); (C) TaubinÕs algorithm in 100 iteration steps with weights of Eq. (2); and (D) the method of Liu et al. in 20 steps of iteration.
burrs in some regions for large numbers of iteration step. It is probably caused by the geometry-dependent weights that may behave ill-conditioned due to bad triangle shapes. Fig. 15 depicts an example for smoothing an apple model with 36,330 vertices and 72,664 faces. The mesh was extracted from data scanned using a CT-scanner. As a comparison, results produced using the kjl algorithm [30] and the volume-preserving approach [19] are also presented. Fig. 15B is generated after 50 steps of iteration and 5 steps of inflation using our method (10 steps of iteration followed by 1-step inflation). Figs. 15C and D illustrate the smoothing results using the approaches of [30] and [19], respectively. The volume-preserving approach of Liu et al. reaches the desired result with the least iteration steps (20 iterations). The proposed algorithm achieves the same visual effect in 50 steps. Though 100 steps have been performed, TaubinÕs algorithm still exhibits undulation on the apple surface. Fig. 16 illustrates another example for comparing the result of proposed approach with that of Goskov et al. [9]. Fig. 16C is generated using their software BEDA with maximal blur while Fig. 16B is the same as Fig. 15B. Employing geometry-independent weights, the presented approach is more efficient for each individual smoothing step than TaubinÕs kjl algorithm using the weights of Eqs. (4) and (5). In addition, the proposed approach uses one inflation for a group of iteration steps, but one inflation step should be added for each round of smoothing iteration for the latter approach in order to prevent shrinkage. It can therefore be conclude that TaubinÕs method is computationally more expensive when
176
G. Li et al. / Graphical Models 66 (2004) 160–179
Fig. 16. Comparison between the proposed method and Goskov et al.Õs algorithm: (A) the original mesh; (B) result generated with the proposed method; and (C) result generated with Goskov et al.Õs algorithm.
performing the same number of smoothing steps. Fig. 10 uses two examples to demonstrate this conclusion. On the other hand, among the three methods, the volume-preserving approach seems most computing-expensive. The computation cost is about 15 times more than that of the proposed approach for performing the same number of iterations. However, its rapidly convergent speed remedies the disadvantage to some extent. To evaluate the smoothing quality of the presented approach, we first define approximation errors. Suppose Mo be a smooth model and M be the mesh generated by adding white noise onto Mo while Ms be the smoothed mesh by denoising M. Considering that both the distance error and the shape error of triangles are desired between Mo and Ms , we linearly refine Mo and Ms one time, respectively, to obtain Mo1 ¼ ðVo1 ; Eo1 ; Fo1 ; Po1 Þ and Ms1 ¼ ðVs1 ; Es1 ; Fs1 ; Ps1 Þ. Noticing that both Mo1 and Ms1 have the same topology and their corresponding vertices have the same index, it is reasonable to define the following three types of errors between Mo and Ms : 1
Eave ðMo ; Ms Þ ¼
jVo j 1 1 X P ðiÞ P 1 ðiÞ; o s 1 jVo j i¼1
Emax ðMo ; Ms Þ ¼ maxPo1 ðiÞ Ps1 ðiÞ; i2Vo1
Emin ðMo ; Ms Þ ¼ minPo1 ðiÞ Ps1 ðiÞ: i2Vo1
Table 1 shows the results using these errors to test TaubinÕs kjl algorithm with three kinds of weights and the presented approach for the Cow model depicted in Fig. 13. In the table, TaubinÕs algorithm with different weights defined in Eqs. (2), (4), and (5), respectively, is examined while the proposed smoothing without and with two steps of TaubinÕs inflation is conducted. From the table we can observe that the presented approach with inflations achieves the best approximation errors for almost all cases in this example. For all the aforementioned examples, a PC with a Pentium 1.6 GHz CPU and 256 MB RAM running Windows XP operating system is used for testing and evaluation.
G. Li et al. / Graphical Models 66 (2004) 160–179
177
Table 1 Comparison of the approximation errors between the results of TaubinÕs algorithm and those of the present approach for the Cow model Iteration steps
25
Error types
ave.
max
min
ave.
max
min
ave.
max
min
Taubin Eq. (2) Taubin Eq. (4) Taubin Eq. (5) Our approach Our approach + 2 inflations
0.00633 0.00891 0.00893 0.00658 0.00640
0.04048 0.03597 0.02545 0.03214 0.02758
0.00027 0.00030 0.00049 0.00030 0.00022
0.00741 0.00926 0.00930 0.00658 0.00640
0.04409 0.04032 0.02847 0.03220 0.02765
0.00020 0.00062 0.00062 0.00030 0.00023
0.00877 0.00988 0.01030 0.00761 0.00714
0.04736 0.10305 0.04883 0.04111 0.03859
0.00030 0.00055 0.00046 0.00053 0.00044
50
100
7. Conclusions In literature, one can find many papers on mesh fairing. However, most of them are developed for smoothing triangular meshes. This paper presents a unified approach for fairing arbitrary two-manifold polygonal meshes. As the dual mesh of a mesh can roughly determined its geometric shape, our smoothing scheme is able to prevent shrinkage and preserve features by interpolating the vertices of the dual mesh. The scheme explores very simple weights, which are only dependent on mesh connectivity, for linear iteration and thus has a concise form. Moreover, we have carefully treated the boundary case in this paper. Mathematical discussion shows that our scheme is convergent for a large range of relaxation factor. As assuming that the face centers of the noisy mesh locate at the corresponding anticipated surface only leads to linear approximation to the original geometry and therefore may introduce shrinkage. To further eliminate this kind of shrinkage, a Taubin-like inflation operator is explored to enlarge the resulted mesh. Experiments have demonstrated the effectiveness and robustness of our approach. In TaubinÕs smoothing algorithm, the inflation operation is built as the dual operator of the Laplacian. We only used it to enlarge the intermediate meshes here and there is no meaning of signal processing. Therefore it is more difficult to determine the relaxation factor l so as to produce an appropriate scale rate. To find a more measurable inflation operator is our future work. In addition, it is also interesting to incorporate crease enhancement techniques to further preserve sharp features. At last, approaches based on signal processing and resampling are probably valuable for arbitrary polygonal meshes in further investigation.
Acknowledgments The work described in this paper is supported by Postdoctoral Foundation of China, Strategic Research Grants Project #7001074 from City University of Hong Kong, National Natural Science Foundations of China (Grant Nos. 60373034, 60021201, 60133020, and 69925204), and the 973 program of China (Grant No. 2002CB312102). We are grateful to Dr. Xinguo Liu for letting us use his program
178
G. Li et al. / Graphical Models 66 (2004) 160–179
of the volume-preserving algorithm. Models of the Venus, mannequin, and Beethoven bust came from http://www.mrl.nyu.edu/~biermann/subdivision/ and the ant mesh was downloaded from http://www.3dcafe.com/asp/ meshes.asp. References [1] D. Barash, A fundamental relationship between bilateral filtering, adaptive smoothing, and the nonlinear diffusion equation, IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (6) (2002) 844–847. [2] E. Catmull, J. Clark, Recursively generated B-spline surfaces on arbitrary topological meshes, Computer-Aided Design 10 (6) (1978) 350–355. [3] U. Clarenz, U. Diewald, M. Rumpf, Anisotropic geometric diffusion in surface processing, in: Proceedings of IEEE Visualization Õ00, 2000, pp. 397–405. [4] M. Desbrun, M. Meyer, P. Schr€ oder, A.H. Barr, Implicit fairing of irregular meshes using diffusion and curvature flow, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1999, pp. 317–324. [5] M. Desbrun, M. Meyer, P. Schr€ oder, A.H. Barr, Discrete differential-geometry operators in nD. Multi-res modeling group TR, Caltech, September 2000. Available from http://www.multires. caltech.edu/pubs/pubs.htm. [6] F. Durand, J. Dorsey, Fast bilateral filtering for the display of high-dynamic-range images, ACM Transactions on Graphics (SIGGRAPH 2001) 21 (3) (2002) 257–266. [7] S. Fleishman, I. Drori, D. Cohen-Or, Bilateral mesh denoising, ACM Transactions on Graphics (SIGGRAPH 2003) 22 (3) (2003). [8] K. Fujiwara, Eigenvalues of Laplacians on a closed Riemannian manifold and its nets, Proceedings of the AMS 123 (8) (1995) 2585–2594. [9] I. Guskov, W. Sweldens, P. Schr€ oder, Multiresolution signal processing, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1999, pp. 325–334. [10] M. Halstead, M. Kass, T. DeRose, Efficient, fair interpolation using Catmull–Clark surfaces, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1993, pp. 35–44. [11] H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, J. Schweitzer, W. Stuetzle, Mesh optimization, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1993, pp. 19–26. [12] T.R. Jones, F. Durand, M. Desbrun, Non-iterative, feature-preserving mesh smoothing, ACM Transactions on Graphics (SIGGRAPH 2003) 22 (3) (2003) 943–949. [13] L. Kobbelt, Discrete fairing, in: Proceedings of the Seventh IMA Conference on the Mathematics of Surface, 1997, pp. 101–131. [14] L. Kobbelt, J. Vorsatz, U. Labsik, H.-P. Seidel, A shrink wrapping approach to remeshing polygonal surfaces, Computer Graphics Forum (EUROGRAPHICS 1999) 18 (3) (1999) 119–130. pffiffiffi [15] L. Kobbelt, 3-Subdivision, in: Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 2000, pp. 103–112. [16] L. Kobbelt, G. Taubin, Digital mesh processing. (SIGGRAPH Course Notes#17, electronic version). Publications Department, ACM Inc., Los Angeles, CA, 2001. [17] A. Khodakovsky, P. Alliez, M. Desbrun, P. Schr€ oder, Near-optimal connectivity encoding of 2manifold polygon, Graphical Models 64 (3/4) (2002) 147–168. [18] X. Liu, H. Bao, P.-A. Heng, T.-T. Wong, Q. Peng, Constrained mesh fairing, Computer Graphics Forum 20 (2) (2001) 115–124. [19] X. Liu, H. Bao, H.-Y. Shum, Q. Peng, A novel volume constrained smoothing method for meshes, Graphical Models 64 (3/4) (2002) 169–182. [20] H.P. Moreton, C.H. Sequin, Functional optimization for fair surface design, in: Computer Graphics Proceedings, ACM SIGGRAPH, 1992, pp. 167–175. [21] Y. Ohtake, A.G. Belyaev, I.A. Bogaevski, Polyhedral surface smoothing with simultaneous mesh regularization, in: Proceedings of the Geometric Modeling and Processing, Hong Kong, 2000, pp. 229–237.
G. Li et al. / Graphical Models 66 (2004) 160–179
179
[22] Y. Ohtake, A.G. Belyaev, I.A. Bogaevski, Mesh regularization and adaptive smoothing, ComputerAided Design 33 (11) (2001) 789–800. [23] Y. Ohtake, A.G. Belyaev, H.-P. Seidel, Mesh smoothing by adaptive and anisotropic Gaussian filter applied to mesh normals, in: G. Greiner, H. Niemann, T. Ertl, B. Girod, H.-P. Seidel (Eds.), Proceedings of Vision, Modeling, and Visualization 2002, Erlangen, Germany, 2002, pp. 203–210. [24] J. Peng, V. Strela, D. Zorin, A simple algorithm for surface denoising, in: IEEE Visualization Õ01, 2001, pp. 107–112. [25] H. Prautzsch, G. Umlauf, A G2 -subdivision algorithm, in: G. Farin, H. Bieri, G. Brunnet, T. DeRose (Eds.), Geometric Modeling, Computing Suppl., vol. 13, Springer-Verlag, Berlin, 1998, pp. 217–224. [26] R. Schneider, L. Kobbelt, Generating fair meshes with G1 boundary conditions, in: Proceedings of Geometric Modeling and Processing, Hong Kong, 2000, pp. 251–261. [27] R. Schneider, L. Kobbelt, Mesh fairing based on an intrinsic PDE approach, CAD 33 (11) (2001) 767– 777. [28] R. Schneider, L. Kobbelt, Geometric fairing of irregular meshes for free-form surface design, Computer Aided Geometry Design 18 (4) (2001) 359–379. [29] T. Tasdizen, R. Whitaker, P. Burchard, S. Osher, Geometric surface smoothing via anisotropic diffusion of normals, in: Procceedings of IEEE Visualization Õ02, 2002, pp. 125–132. [30] G. Taubin, A signal processing approach to fair surface design, in: Computer Graphics Proceedings, ACM SIGGRAPH, 1995, pp. 351–358. [31] G. Taubin, Dual Mesh Resampling, Pacific Graphics, Tokyo, Japan, 2001, pp. 180–189. [32] G. Taubin, Linear anisotropic mesh filtering. IBM Research Report RC-22213, 2001. [33] C. Tomasi, R. Manduchi, Bilateral filtering for gray and color image, in: Proceedings of the 1998 IEEE International Conference on Computer Vision, 1998. [34] J. Vollmer, R. Mencl, H. Muller, Improved Laplacian smoothing of noisy surface meshes, Computer Graphics Forum (EUROGRAPHICS) 18 (3) (1999) 131–138. [35] W. Welch, A. Witkin, Variational surface modeling, in: Computer Graphics Proceedings, ACM SIGGRAPH, 1992, pp. 157–166. [36] W. Welch, A. Witkin, Free-form shape design using triangulated surfaces, in: Computer Graphics Proceedings, ACM SIGGRAPH, 1994, pp. 247–256. [37] G. Westgaard, H. Nowacki, A process for surface fairing in irregular meshes, Computer Aided Geometry Design 18 (7) (2001) 619–638. [38] Z. Wood, M. Desbrun, P. Schr€ oder, D. Breen, Semi-regular mesh extraction from volumes, in: Proceedings of IEEE Visualization, 2000, pp. 275–282. [39] H. Zhang, E. Fiume, Mesh smoothing with shape or feature preservation, in: J. Vince, R. Earnshaw (Eds.), Advances in Modeling, Animation, and Rendering, 2002, pp. 167–182.