Fairing Implicit of Irregular Meshes using Diffusion and Curvature Flow
Mathieu Desbrun Mark Meyer Caltech Peter Schr der o Alan H. Barr
Abstract
In this paper, we develop methods to rapidly remove rough features from irregularly triangulated data intended to portray a smooth surface. The main task is to remove undesirable noise and uneven edges while retaining desirable geometric features. The problem arises mainly when creating high-delity computer graphics objects using imperfectly-measured data from the real world. Our approach contains three novel features: an implicit integration method to achieve efciency, stability, and large time-steps; a scale-dependent Laplacian operator to improve the diffusion process; and nally, a robust curvature ow operator that achieves a smoothing of the shape itself, distinct from any parameterization. Additional features of the algorithm include automatic exact volume preservation, and hard and soft constraints on the positions of the points in the mesh. We compare our method to previous operators and related algorithms, and prove that our curvature and Laplacian operators have several mathematically-desirable qualities that improve the appearance of the resulting surface. In consequence, the user can easily select the appropriate operator according to the desired type of fairing. Finally, we provide a series of examples to graphically and numerically demonstrate the quality of our results.
1 Introduction
While the mainstream approach in mesh fairing has been to enhance the smoothness of triangulated surfaces by minimizing computationally expensive functionals, Taubin [Tau95] proposed in 1995 a signal processing approach to the problem of fairing arbitrary topology surface triangulations. This method is linear in the number of vertices in both time and memory space; large arbitrary connectivity meshes can be handled quite easily and transformed into visually appealing models. Such meshes appear more and more frequently due to the success of 3D range sensing approaches for creating complex geometry [CL96]. Taubin based his approach on dening a suitable generalization of frequency to the case of arbitrary connectivity meshes. Using a discrete approximation to the Laplacian, its eigenvectors become the frequencies of a given mesh. Repeated application of the resulting linear operator to the mesh was then employed to tailor the frequency content of a given mesh. Closely related is the approach of Kobbelt [Kob97], who considered similar discrete approximations of the Laplacian in the construction of fair interpolatory subdivision schemes. In later work this was extended to the arbitrary connectivity setting for purposes of multiresolution editing [KCVS98]. The success of these techniques is largely based on their simple implementation and the increasing need for algorithms which can process the ever larger meshes produced by range sensing tech {mathieu|mmeyer|ps|barr}@cs.caltech.edu.
(a) (b) Figure 1: (a): Original 3D photography mesh (41,000 vertices). (b): Smoothed version with the scale-dependent operator in two integration step with dt = 5 105 , the iterative linear solver (PBCG) converges in 10 iterations. All the images in this paper are at-shaded to enhance the faceting effect. niques. However, a number of issues in their application remain open problems in need of a more thorough examination. The simplicity of the underlying algorithms is based on very basic, uniform approximations of the Laplacian. For irregular connectivity meshes this leads to a variety of artifacts such as geometric distortion during smoothing, numerical instability, problems of slow convergence for large meshes, and insufcient control over global behavior. The latter includes shrinkage problems and more precise shaping of the frequency response of the algorithms. In this paper we consider more carefully the question of numerical stability by observing that Laplacian smoothing can be thought of as time integration of the heat equation on an irregular mesh. This suggests the use of implicit integration schemes which lead to unconditionally stable algorithms allowing for very large time steps. At the same time the necessary linear system solvers run faster than explicit approaches for large meshes. We also consider the question of mesh parameterization more carefully and propose the use of discretizations of the Laplacian which take the underlying parameterization into account. The resulting algorithms avoid many of the distortion artifacts resulting from the application of previous methods. We demonstrate that this can be done at only a modest increase in computing time and results in smoothing algorithms with considerably higher geometric delity. Finally a more careful analysis of the underlying discrete differential geometry is used to derive a curvature ow approach which satises crucial geometric properties. We detail how these different operators act on meshes, and how users can then decide which one is appropriate in their case. If the user wants to, at the same time, smooth the shape of an object and equalize its triangulation, a scale-dependent diffusion must be used. On the other hand, if only the shape must be ltered without affecting the sampling rate, then curvature ow has all the desired properties. This allows us to propose a novel class of efcient smoothing algorithms for arbitrary connectivity meshes.
2 Implicit fairing
In this section, we introduce implicit fairing, an implicit integration of the diffusion equation for the smoothing of meshes. We will demonstrate several advantages of this approach over the usual ex-
plicit methods. While this section is restricted to the use of a linear approximation of the diffusion term, implicit fairing will be used as a robust and efcient numerical method throughout the paper, even for non-linear operators. We start by setting up the framework and dening our notation.
can be constructed by integrating the diffusion equation with a simple explicit Euler scheme, yielding: X n+1 = (I + dt L )X n . (8)
2.1 Notation and denitions
In the remainder of this paper, X will denote a mesh, xi a vertex of this mesh, and ei j the edge (if existing) connecting xi to x j . We will call N1 (i) the neighbors (or 1-ring neighbors) of xi , i.e., all the vertices x j such that there exists an edge ei j between xi and x j (see Figure 9(a)). In the surface fairing literature, most techniques use constrained energy minimization. For this purpose, different fairness functionals have been used. The most frequent functional is the total curvature of a surface S :
With the umbrella operator, the stability criterion requires dt < 1. If the time step does not satisfy this criterion, ripples appear on the surface, and often end up creating oscillations of growing magnitude over the whole surface. On the other hand, if this criterion is met, we get smoother and smoother versions of the initial mesh as n grows.
2.3 Time-shifted evaluation
The implementation of this previous explicit method, called forward Euler method, is very straightforward [Tau95] and has nice properties such as linear time and linear memory size for each ltering pass. Unfortunately, when the mesh is large, the time step restriction results in the need to perform hundreds of integrations to produce a noticeable smoothing, as mentioned in [KCVS98]. Implicit integration offers a way to avoid this time step limitation. The idea is simple: if we approximate the derivative using the new mesh (instead of using the old mesh as done in explicit methods), we will get to the equilibrium state of the PDE faster. As a result of this time-shifted evaluation, stability is obtained unconditionally [PTVF92]. The integration is now: X n+1 = X n + dt L (X n+1 ). Performing an implicit integration, this time called backward Euler method, thus means solving the following linear system: (I dt L )X n+1 = X n . (9) This apparently minor change allows the user not to worry about practical limitations on the time step. Consequent smoothing will then be obtained safely by increasing the value dt. But solving a linear system is the price to pay.
E (S ) =
S
2 + 2 d S . 1 2
(1)
This energy can be estimated on discrete meshes [WW94, Kob97] by tting local polynomial interpolants at vertices. However, principal curvatures 1 and 2 depend non-linearly on the surface S . Therefore, many practical fairing methods prefer the membrane functional or the thin-plate functional of a mesh X:
Emembrane (X) = Ethin plate (X) =
1 2
1 2
2 2 Xu + Xv dudv
(2) (3)
2 2 2 Xuu + 2 Xuv + Xvv dudv.
Note that the thin-plate energy turns out to be equal to the total curvature only when the parameterization (u, v) is isometric. Their respective variational derivatives corresponds to the Laplacian and the second Laplacian:
2.4 Solving the sparse linear system
Fortunately, this linear system can be solved efciently as the matrix A = I dt L is sparse: each line contains approximately six non-zero elements if the Laplacian is expressed using Equ. (7) since the average number of neighbors on a typical triangulated mesh is six. We can use a preconditioned bi-conjugate gradient (PBCG) to iteratively solve this system with great efciency1 . The PBCG is based on matrix-vector multiplies [PTVF92], which only require linear time computation in our case thanks to the sparsity of the matrix A. We review in Appendix A the different options we chose for the PBCG in order to have an efcient implementation for our purposes.
L (X) = Xuu + Xvv L 2 (X) = L L (X) = Xuuuu + 2 Xuuvv + Xvvvv .
(4) (5)
For smooth surface reconstruction in vision, a weighted average of these derivatives has been used to fair surfaces [Ter88]. For meshes, Taubin [Tau95] used signal processing analysis to show that a combination of these two derivatives of the form: ( + )L L 2 can provide a Gaussian ltering that minimizes shrinkage. The constants and must be tuned by the user to obtain this non-shrinking property. We will refer to this technique as the | algorithm.
2.5 Interpretation of the implicit integration
Although this implicit integration for diffusion is sound as is, there are useful connections with other prior work. We review the analogies with signal processing approaches and physical simulation. 2.5.1 Signal processing In [Tau95], Taubin presents the explicit integration of diffusion with a signal processing point of view. Indeed, if X is a 1D signal of a given frequency : X = ei , then L (X) = 2 X. Thus, the transfer function for Equ. (8) is 1 dt2 , as displayed in Figure 2(a) as a solid line. We can see that the higher the frequency , the stronger the attenuation will be, as expected. The previous lter is called FIR (for Finite Impulse Response) in signal processing. When the diffusion process is integrated using implicit integration, the lter in Equ. (9) turns out to be an Innite Impulse Response lter. Its transfer function is now 1/(1 + dt2 ), depicted in Figure 2(a) as a dashed line. Because this lter is always in [0, 1], we have unconditional stability.
use a bi-conjugate gradient method to be able to handle non symmetric matrices, to allow the inclusion of constraints (see Section 2.7).
1 We
2.2 Diffusion equation for mesh fairing
As we just pointed out, one common way to attenuate noise in a mesh is through a diffusion process: X = L (X). t (6)
By integrating equation 6 over time, a small disturbance will disperse rapidly in its neighborhood, smoothing the high frequencies, while the main shape will be only slightly degraded. The Laplacian operator can be linearly approximated at each vertex by the umbrella operator (we will use this approximation in the current section for the sake of simplicity, but will discuss its validity in section 4), as used in [Tau95, KCVS98]:
L (xi ) =
1 m
jN1 (i)
x j xi
(7)
where x j are the neighbors of the vertex xi , and m = #N1 (i) is the number of these neighbors (valence). A sequence of meshes (X n )
Attenuation
1
Attenuation Explicit filter Implicit filter
1
Filter for ten explicit integrations Filter for ten implicit integrations
0.8
0.8
in the context of a xed mesh, though: amplifying frequencies requires renement of the mesh to offer a good discretization.
0.6
0.6
0.4 0.4 0.2 0.2 0 0 0.5 1 1.5 2 2.5 3 0 0 0.2 0.4 0.6 0.8 1
2.7 Constraints
We can put hard and soft constraints on the mesh vertex positions during the diffusion. For the user, it means that a vertex or a set of vertices can be xed so that the smoothing happens only on the rest of the mesh. This can be very useful to retain certain details in the mesh. A vertex xi will stay xed if we impose L (xi ) = 0. More complicated constraints are also possible [BW98]. For example, vertices can be constrained along an axis or on a plane by modifying the PBCG to keep these constraints enforced during the linear solver iterations. We can also easily implement soft constraints: each vertex can be weighted according to the desired smoothing that we want. For instance, the user may want to smooth a part of a mesh less than another one, in order to keep desirable features while getting a smoother version. We allow the assignment of a smoothing value between 0 and 1 to attenuate the smoothing spatially: this is equivalent to choosing a variable factor on the mesh, and happens to be very useful in practice. Entire regions can be spray painted interactively to easily assign this special factor.
Frequency
Frequency
(a) (b) Figure 2: Comparison between (a) the explicit and implicit transfer function for dt = 1, and (b) their resulting transfer function after 10 integrations. By rewriting Equ. (9) as: X n+1 = (I dt L )1 X n , we also note that our implicit ltering is equivalent to I + dt L + (dt)2 L 2 + ..., i.e., standard explicit ltering plus an innite sequence of higher order ltering. Contrary to the explicit approach, one single implicit ltering step performs global ltering. 2.5.2 Mass-spring network Smoothing a mesh by minimizing the membrane functional can be seen as a physical simulation of a mass-spring network with zerorest length springs that will shrink to a single point in the limit. Recently, Baraff and Witkin [BW98] presented an implicit method to allow large time steps in cloth simulation. They found that the use of an implicit solver instead of the traditional explicit Euler integration considerably improves computational time while still being stable for very stiff systems. Our method compares exactly to theirs, but used for meshes and for a different PDE. We therefore have the same advantages of using an implicit solver over the usual explicit type: stability and efciency when signicant ltering is called for.
2.8 Discussion
Even if adding a linear solver step to the integration of the diffusion equation seems to slow down the problem at rst glance, it turns out that we gain signicantly by doing so. For instance, the implicit integration can be performed with an arbitrary time step. Since the matrix of the system is very sparse, we actually obtain computational time similar or better than the explicit methods. In the following table, we indicate the number of iterations of the PBCG method for different meshes and it can be seen that the PBCG is more efcient when the smoothing is high. These timings were performed on an SGI High Impact Indigo2 175MHz R10000 processor with 128M RAM.
Mesh Horse Dragon Isis Bunny Buddha Nb of faces 42,000 42,000 50,000 66,000 290,000 dt = 10 8 iterations (2.86s) 8 iterations (2.98s) 9 iterations (3.84s) 7 iterations (4.53s) 5 iterations (13.78s) dt = 100 37 iterations (12.6s) 39 iterations (13.82s) 37 iterations (15.09s) 35 iterations (21.34s) 28 iterations (69.93s)
2.6 Filter improvement
Now that the method has been set up for the usual diffusion equation, we can consider other equations that may be more appropriate or may give better visual results for smoothing when we use implicit integration. We have seen in Section 2.1 that both L and L 2 have been used with success in prior work [Ter88, Tau95, KCVS98]. When we use implicit integration, as Figure 3(a) shows, the higher the power of the Laplacian, the closer to a low-pass lter we get. In terms of frequency analysis, it is a better lter. Unfortunately, the matrix becomes less and less sparse as more and more neighbors are involved in the computation. In practice, we nd that L 2 is a very good trade-off between efciency and quality. Using higher orders affects the computational time signicantly, while not always producing signicant improvements. We therefore recommend using (I + dt L 2 )X n+1 = X n for implicit smoothing (a precise denition of the umbrella-like operator for L 2 can be found in [KCVS98]).
1 0.8
(I-L) -1 (I-L2 ) -1 (I-L3 ) -1 -1 (I-L4 )
1.2
1
Implicit filter Constant filter Resulting convolution
0.6
0.8
0.6 0.4 0.4 0.2 0.2
0 0 0.5 1 1.5 2 2.5 3
0 0 0.5 1 1.5 2 2.5 3
(a) (b) Figure 3: (a): Comparison between lters using L , L 2 , L 3 , and L 4 . (b): The scaling to preserve volume creates an amplication of all frequencies; but the resulting lter (diffusion+scaling) only amplies low frequencies to compensate for the shrinking of the diffusion. We also tried to use a linear combination of both L and L 2 . We obtained interesting results like, for instance, amplication of low or middle frequencies to exaggerate large features (refer to [GSS99] for a complete study of feature enhancement). It is not appropriate
To be able to compare the results with the explicit method, one has to notice that one iteration of the PBCG is only slightly more time consuming than one integration step using an explicit method. Therefore, we can see in the following results that our implicit fairing takes about 60% less time than the explicit fairing for a ltering of dt = 100, as we get about 33 iterations compared to the 100 integration steps required in the explicit case. We have found this behavior to be true for all the other meshes as well. The advantage of the implicit method in terms of computational speed becomes more obvious for large meshes and/or high smoothing value. In terms of quality, Figure 4(b) and 4(c) demonstrate that both implicit and explicit methods produce about the same visual results, with a slightly better smoothness for the implicit fairing. Note that we use 10 explicit integrations of the umbrella operator with dt = 1, and 1 integration using the implicit integration with dt = 10 to approximate the same results. Therefore, there is a denite advantage in the use of implicit fairing over the previous explicit methods. Moreover, the remainder of this paper will make heavy use of this method and its stability properties.
3 Automatic anti-shrinking fairing
Pure diffusion will, by nature, induce shrinkage. This is inconvenient as this shrinking may be signicant for aggressive smoothing. Taubin proposed to use a linear combination of L and L L to amplify low frequencies in order to balance the natural shrinking. Unfortunately, the linear combination depends heavily on the mesh in practice, and this requires ne tuning to ensure both stable
(a) (b) (c) (d) Figure 4: Stanford bunnies: (a) The original mesh, (b) 10 explicit integrations with dt = 1, (c) 1 implicit integration with dt = 10 that takes only 7 PBCG iterations (30% faster), and (d) 20 passes of the | algorithm, with = 0.6307 and = 0.6732. The implicit integration results in better smoothing than the explicit one for the same, or often less, computing time. If volume preservation is called for, our technique then requires many fewer iterations to smooth the mesh than the | algorithm. and non-shrinking results. In this section, we propose an automatic solution to avoid this shrinking. We preserve the zeroth moment, i.e., the volume, of the object. Without any other information on the mesh, we feel it is the most reasonable invariant to preserve, although surface area or other invariants can be used. for fairing. Indeed, no parameters need be tuned to ensure stability or to have exact volume preservation. This is a major advantage over previous techniques. Yet, we retain all of the advantages of previous methods, such as constraints [Tau95] and the possibility of accelerating the fairing via multigrid [KCVS98], while additionally offering stability and efciency. This technique also dramatically reduces the computing time over Taubins anti-shrinking algorithm: as demonstrated in Figure 4(c) and 4(d), using the | algorithm may preserve the volume after ne tuning, but one iteration will only slightly smooth the mesh. The rest of this paper exploits both automatic anti-shrinking and implicit fairing techniques to offer more accurate tools for fairing.
3.1 Volume computation
As we have a mesh given in terms of triangles, it is easy to compute the interior volume. This can be done by summing the volumes of all the oriented pyramids centered at a point in space (the origin, for instance) and with a triangle of the mesh as a base. This computation has a linear complexity in the number of triangles [LK84]. For the readers convenience, we give the expression of the volume of a mesh in the following equation, where x1 , x2 and x3 are the three k k k vertices of the kth triangle: nbFaces 1 V= (10) gk Nk 6 k=1 where g = (x1 + x2 + x3 )/3 and Nk = x1 x2 x1 x3 k k k k k k k
4 An accurate diffusion process
Up to this section, we have relied on the umbrella operator (Equ. (7)) to approximate the Laplacian on a vertex of the mesh. This particular operator does not truly represent a Laplacian in the physical meaning of this term as we are about to see. Moreover, simple experiments on smooth meshes show that this operator, using explicit or implicit integration, can create bumps or pimples on the surface, instead of smoothing it. This section proposes a sounder simulation of the diffusion process, by dening a new approximation for the Laplacian and by taking advantage of the implicit integration.
3.2 Exact volume preservation
After an integration step, the mesh will have a new volume V n . We then want to scale it back to its original volume V 0 to cancel the shrinking effect. We apply a simple scale on the vertices to achieve this. By multiplying all the vertex positions by = (V 0 /V n )1/3 , the volume is guaranteed to go back to its original value. As this is a simple scaling, it is harmless in terms of frequencies. To put it differently, this scaling corresponds to a convolution with a scaled Dirac in the frequency domain, hence it amplies all the frequencies in the same way to change the volume back. The resulting lter, after the implicit smoothing and the constant amplication lter, amplies the low frequencies of the original mesh to exactly compensate for the attenuation of the high frequencies, as sketched on Figure 3(b). The overall complexity for volume preservation is then linear. With such a process, we do not need to tweak parameters: the anti-shrinking lter is automatically adapted to the mesh and to the smoothing, contrary to previous approaches. Note that hard constraints dened in the previous section are applied before the scaling and do not result in xed points anymore: scaling alters the absolute, but not the relative position. We can generalize this re-scaling phase to different invariants. For instance, if we have to smooth height elds, it is more appropriate to take the invariant as being the volume enclosed between the height eld and a reference plane, which changes the computations only slightly. Likewise, for surfaces of revolution, we may change the way the scaling is computed to exploit this special property. We can also preserve the surface area if the mesh is a non-closed surface. However, in the absence of specic characteristics, preserving the volume gives nice results. According to specic needs, the user can select the appropriate type of invariant to be used.
4.1 Inadequacy of the umbrella operator
The umbrella operator, used in the previous sections corresponds to an approximation of the Laplacian in the case of a specic parameterization [KCVS98]. This means that the mesh is supposed to have edges of length 1 and all the angles between two adjacent edges around a vertex should be equal. This is of course far from being true in actual meshes, which contain a variety of triangles of different sizes. Treating all edges as if they had equal length has signicant undesired consequences for the smoothing. example, For the Laplacian can be the same for two very different congurations, corresponding to different frequencies as depicted in Figure 5. This distorts the ltering signicantly, as high frequencies may be considered as low ones, and vice-versa. Nevertheless, the advantage of the umbrella operator is that it is normalized: the time step for integration is always 1, which is very convenient. But we want a more accurate diffusion process to smooth meshes consistently, in order to more carefully separate high from low frequencies.
3.3 Discussion
When we combine both methods of implicit integration and antishrinking convolution, we obtain an automatic and efcient method
(a) (b) Figure 5: Frequency confusion: the umbrella operator is evaluated as the vector joining the center vertex to the barycenter of its neighbors. Thus, cases (a) and (b) will have the same approximated Laplacian even if they represent different frequencies.
We need to dene a discrete Laplacian which is scale dependent, to better approximate diffusion. However, if we use explicit integration [Tau95], we will suffer from a very restricted stability criterion. It is well known [PTVF92] that the time step for a parabolic PDE like Equ. (6) depends on the square of the smallest length scale (here, the smallest edge length min(|e|)): dt min(|e|)2 2
y Samples on grid A
0.8
y Samples of grid B
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
2
x
0 0
2
x
(a)
Error
3.5e-06 3e-06 2.5e-06 2e-06 1.5e-06 1e-06 5e-07 0 0 0.5 1 1.5 2
(b)
Grid A, regular FD Grid B, regular FD Grid B, extended FD
X
This limitation is a real concern for large meshes with small details, since an enormous number of integration steps will have to be performed to obtain noticeable smoothing. This is intractable in practice. With implicit integration explained in Section 2, we can overcome this restriction and use a much larger time step while still achieving good smoothing, saving considerable computation. In the next two paragraphs we present one design of a good approximation for the Laplacian.
Xi Xi-1 Xi+1
2.5
3
x
ui-1
ui
u
ui+1
4.2 Simulation of the 1D heat equation
The 1D case of a diffusion equation corresponds to the heat equation xt = xuu . It is therefore worth considering this example as a test problem for higher dimensional ltering. To do so, we use Milnes test presented in [Mil95]. Milne compared two cases of the same initial problem: rst, the problem is solved on a regular mesh on [0, 1], and then on an irregular mesh, taken to consist of a 1 uniform coarse grid of cells on [0, 1] with each of the cells in [ 2 , 1] subdivided into two ne cells as depicted in Figure 6(a) and 6(b). With such a conguration, classical nite difference coefcients for second derivatives can be used on each cell, except for the middle one which does not have centered neighbors. Milne shows that if no particular care is taken for this peripheral cell, it introduces a noise term that creates large inaccuracies larger than if the mesh was represented uniformly at the coarser resolution! But if we t a quadratic spline at this cell to approximate the second derivative, then the noise source disappears and we get more accurate results than with a constant coarse resolution (see the errors created in each case in one iteration of the heat equation in Figure 6(c)). This actually corresponds to the extension of nite difference computations for irregular meshes proposed by Fornberg [For88]: to compute the FD coefcients, just t a quadratic function at the sample point and its two immediate neighbors, and then return the rst and second derivative of that function as the approximate derivatives. For three points spaced and apart (see Figure 6(d)), we get the 1D formula: 2 xi1 xi xi+1 xi + . + Note that when = , we nd the usual nite difference formula. (xuu )i =
(c) (d) Figure 6: Test on the heat equation: (a) regular sampling vs. (b) irregular sampling. Numerical errors in one step of integration (c): using the usual FD weight on an irregular grid to approximate second derivatives creates noise, and gives a worse solution than on the coarse grid, whereas extended FD weights offer the expected behavior. (d) Three unevenly spaced samples of a function and corresponding quadratic tting for extended FD weights. compute them initially using the current edges lengths and keep their values constant during the PBCG iterations. In practice, we have not noted any noticeable drawbacks from this linearization. We can even keep the same coefcients for a number of (or all) iterations: it will correspond to a ltering relative to the initial mesh instead if the current mesh. For the same reason as before, we also recommend the use of the second Laplacian for higher quality smoothing without signicant increase in computation time. As demonstrated in Figure 7, the scale-dependent umbrella operator deals better with irregular meshes than the umbrella operator: no spurious artifacts are created. We also applied this operator to noisy data sets from 3D photography to obtain smooth meshes (see Figure 1 and 12). The number of iterations needed for convergence depends heavily on the ratio between minimum and maximum edge lengths. For typical smoothing and for meshes over 50000 faces, the average number of iterations we get is 20. Nevertheless, we still observe undesired behavior on at surfaces: vertices in at areas still slide during smoothing. Even though this last formulation generally reduces this problem, we may want to keep a at area intact. The next section tackles this problem with a new approach.
5 Curvature ow for noise removal
In terms of differential equations, diffusion is a close relative of curvature ow. In this section, we rst explore the advantages of using curvature ow over diffusion, and then propose an efcient algorithm for noise removal using curvature ow.
4.3 Extension to 3D
The umbrella operator suffers from this problem of large inaccuracies for irregular meshes as the same supposedly constant parameterization is used (Figure 7 shows such a behavior). Surprisingly, a simple generalization of the previous formula valid in 1D corresponds to a known approximation of the Laplacian. Indeed, Fujiwara [Fuj95] presents the following formula:
5.1 Diffusion vs. curvature ow
The Laplacian of the surface at a vertex has both normal and tangential components. Even if the surface is locally at, the Laplacian approximation will rarely be the zero vector [KCVS98]. This introduces undesirable drifting over the surface, depending on the parameterization we assume. We in effect fair the parameterization of the surface as well as the shape itself (see Figure 10(b)). We would prefer to have a noise removal procedure that does not depend on the parameterization. It should use only intrinsic properties of the surface. This is precisely what curvature ow does. Curvature ow smoothes the surface by moving along the surface normal n with a speed equal to the mean curvature : xi = i ni . t (12)
L (xi ) =
2 E
jN1
x j xi , with E = |ei j | (i)
jN1 (i)
|ei j |.
(11)
where |ei j | is the length of the edge ei j . Note that, when all edges are of size 1, this reduces to the umbrella operator (7). We will then denote this new operator as the scale-dependent umbrella operator. Unfortunately, the operator is no longer linear. But during a typical smoothing, the length of the edges does not change dramatically. We thus make the approximation that the coefcients of the matrix A = (I dt L ) stay constant during an integration step. We can
area of all the triangles of the 1-ring neighbors as sketched in Figure 9(a). Note that this area A uses cross products of adjacent edges, and thus implicitly contains information on local normal vectors. The complete derivation from the continuous formulation to the discrete case is shown in Appendix B. We nd the following discrete expression through basic differentiation: n = (a) (b) (c) (d) Figure 7: Application of operators to a mesh: (a) mesh with different sampling rates, (b) the umbrella operator creates a signicant distortion of the shape, but (c) with the scale-dependent umbrella operator, the same amount of smoothing does not create distortion or artifacts, almost like (d) when curvature ow is used. The small features such as the nose are smoothed but stay in place. Other curvatures can of course be used, but we will stick to the mean curvature: = (1 + 2 )/2 in this paper. Using this procedure, a sphere with different sampling rates should stay spherical under curvature ow as the curvature is constant. And we should also not get any vertex sliding when an area is at as the mean curvature is then zero. There are already different approaches using curvature ow [Set96], and even mixing both curvature ow and volume preservation [DCG98] to smooth object appearance, but mainly in the context of level-set methods. They are not usable on a mesh as is. Next, we show how to approximate curvature consistently on a mesh and how to implement this curvature ow process with our implicit integration for efcient computations. 1 4A
jN1 (i)
(cot j + cot j )(x j xi )
(14)
where j and j are the two angles opposite to the edge in the two triangles having the edge ei j in common (as depicted in Figure 9(b)), and A is the sum of the areas of the triangles having xi as a common vertex.
Xi
Xi Xj-1
j
Xj
Aj
A j
e ij
Xj
j
Xj+1
(a) (b) Figure 9: A vertex xi and its adjacent faces (a), and one term of its curvature normal formula (b). Note the interesting similarity with [PP93]. We obtain almost the same equation, but with a completely different derivation than theirs, which was using energies of linear maps. The same remark stands for [DCDS97] since they also nd the same kind of expression as Equ. (14) for their functional, but using this time piecewise linear harmonic functions.
5.2 Curvature normal calculation
It seems that all the formulations so far have a non-zero tangential component on the surface. This means that even if the surface is at around a vertex, it may move anyway. For curvature ow, we dont want this behavior. A good idea is to check the divergence of the normal vector, as it is the denition of mean curvature ( = div n): if all the normals of the faces around a vertex are the same, this vertex should not move then (zero curvature). Having this in mind, we have selected the following differential geometry denition of the curvature normal n: A =n 2A (13)
5.3 Boundaries
For non-closed surfaces or surfaces with holes, we can dene a special treatment for vertices on boundaries. The notion of mean curvature does not make sense for such vertices. Instead, we would like to smooth the boundary, so that the shape of the hole itself gets rounder and rounder as iterations go. We can then use for instance Equ. (11) restricted to the two immediate neighbors which will smooth the boundary curve itself. Another possible way is to create a virtual vertex, stored but not displayed, initially placed at the barycenter of all the vertices placed on a closed boundary. A set of faces adjacent to this vertex and connecting the boundary vertices one after the other are also virtually created. We can then use the basic algorithm without any special treatment for the boundary as now, each vertex has a closed area around it.
where A is the area of a small region around the point P where the curvature is needed, and is the derivative with respect to the (x, y, z) coordinates of P. With this denition, we will have the zero vector for a at area. As proven in Figure 8, we see that moving the center vertex xi on a at surface does not change the surface area. On the other hand, moving it above or below the plane will always increase the local area. Hence, we have the desired property of a null area gradient for a locally at surface, whatever the valence, the aspect ratio of the adjacent faces, or the edge lengths around the vertex.
xi xi xi xi
5.4 Implementation
Similarly to Section 4, we have a non-linear expression dening the curvature normal. We can however proceed in exactly the same way, as the changes induced in a time step will be small. We simply compute the non-zero coefcients of the matrix I dtK, where K represents the matrix of the curvature normals. We then successively solve the following linear system: (I dtK) X n+1 = X n . We can use preconditioning or constraints, just as before as everything is basically the same except for the local approximation of the speed of smoothing. As shown on Figure 10, a sphere with different triangle sizes will remain the same sphere thanks to both the curvature ow and the volume preservation technique. In order for the algorithm to be robust, an important test must be performed while the matrix K is computed: if we encounter a face of zero area, we must skip it. As we divide by the area of the face, degenerate triangles are to be treated specially. Mesh decimation to eliminate all degenerate triangles can also be used as suggested in [PP93].
Figure 8: The area around a vertex xi lying in the same plane as its 1-ring neighbors does not change if the vertex moves within the plane, and can only increase otherwise. Being a local minimum, it thus proves that the derivative of the area with respect to the position of xi is zero for at regions. To derive the discrete version of this curvature normal, we select the smallest area around a vertex xi that we can get, namely the
(a) (b) (c) (d) Figure 10: Smoothing of spheres: (a) The original mesh containing two different discretization rates. (b) Smoothing with the umbrella operator introduces sliding of the mesh and unnatural deformation, which is largely attenuated when (c) the scale-dependent version is used, while (d) curvature ow maintains the sphere exactly.
5.5 Normalized version of the curvature operator
We can now write the equivalent of the umbrella operator, but for the curvature normal. Since the new formulation has nice properties, we can create a normalized version that could be used in an explicit integration for quick smoothing. The normalization will bring the eigenvalues back in [1, 0] so that a time step up to 1 can be used in explicit integration methods. Its expression is simply: ( n)normalized = 1 j (cot lj + cot rj )
we developed a curvature ow process. The same implicit integration is used for this new operator that now offers a smoothing only depending on intrinsic geometric properties, without sliding on at areas and with preserved curvature for constant curvature areas. The user can make use of all these different tools according to the mesh to be smoothed. We believ...