Shape Operator Metric for Surface Normal Approximation

Summary. This work deals with the problem of practical mesh generation for surface normal approximation. Part of its contribution is in presenting previous work in a uniﬁed framework. A new algorithm for surface normal approximation is then introduced which improves upon existing ones in a number of aspects. In particular, it produces better approximations of surfaces both in practice and in the theoretical limit regime. Additionally, it resolves in a simple way some of the problems that previous methods for surface approximation suﬀered from.


Introduction
Computing high-quality approximating meshes from surfaces is an important problem in computational geometry, with many practical implications.Although the approximation criteria can vary greatly, often, approximating either surface position, or a surface's normal field can be a good criteria in practice.As has been argued elsewhere [7,14], approximating a surface while minimizing normal approximation error can be useful in many applications.
There is a considerable body of previous work that deals with the surface approximation problem.Some notable examples include -nets [6], for surface and normal approximation, the Quadric Error Metric algorithm (QEM) [9] for surface approximation, and Variational Shape Approximation (VSA) [14], for surface and normal approximation.
In this paper, it is first discussed how the above three methods can be interpreted from within a unified framework.In this interpretation, they are essentially all minimizing a k-means like energy, where only the distance metrics are different.Interestingly enough, in the limit, for smooth surfaces, these distance measures converge to each other.The other key difference that is explored is how these means are used by the different algorithms to produce the output triangulation.
Next, a novel distance measure is proposed (Shape Operator Metric, or SOM), with a corresponding algorithm, that fills a natural gap in this framework.In particular, like VSA, it is designed for normal approximation.But, like QEM, it does not require a region-triangulation step.Such a step can complicate the implementation and, as it is discussed, it introduces a constant factor of inefficiency close to 2, in the limit of approximation.

Framework
Considered here are meshing algorithms for surface approximation that try to either meet a uniform error bound, or minimize the average approximation error over a surface M (minimizing error in the L ∞ or L 2 sense respectively).These kinds of algorithms can be naturally described as an optimization problem: over both a set of means {p j } (points on the surface), and a corresponding partition {V j } of M composed of the Voronoi cells of {p j } with respect to a chosen distance function D X .Optimal Voronoi partitions have in all (except perhaps the rarest) cases neither the shape nor the topology of a triangle mesh.Some further step is generally necessary before producing a triangle mesh as output.In the sequel, a meshing algorithm is referred to as a primal algorithm if it discretizes the boundaries of Voronoi cells, triangulates their interior, and outputs this set of triangles, as in [14]; while an algorithm in which the means are instead directly connected using the dual topology of the partition {V j } to produce a triangle mesh, as in [6], will be denoted as a dual algorithm.
Apart from the added algorithmic complexity, primal algorithms have an inherent approximation inefficiency in the limit.Roughly speaking, in smooth surface regions, in the limit, the relative sizes and aspect ratios of the Voronoi regions are optimized by minimizing the above energies.These relative sizing and aspect ratios will be maintained under mesh duality.But these sizes and aspect ratios are altered (within a constant factor) when the Voronoi regions are triangulated.The limit regime is explored in more detail in Appendix B, while the non-limit case is discussed experimentally in Section 5.
The algorithms of [6,9,14], as well as the one introduced in Sec. 3, all fit into this framework.In particular, the method in [14] introduced the idea of directly optimizing energies with the above form using a k-means/Lloyd-Max type algorithm.It then applies a primal meshing step to the resulting partition.
In the work of [6] one finds a set of means {p j } with bounds on the energy of Eq. 1.The means are then connected in a dual triangulation.
The QEM method of [9] applies a sequence of edge collapses to the input mesh, which can essentially be interpreted as an attempt to minimize 2. In particular, upon completion of the QEM algorithm, each vertex of the output triangulation can be thought of as a mean with a region of the surface associated to it: the portion of the surface that it uses to evaluate its associated (quadric) error (with adjacent regions slightly overlapping).In this sense QEM can be considered a dual algorithm.Even though the connectivity of the triangle mesh is not directly related to the Voronoi regions of Eq. 1, its limit behavior is analogous.

Surface Approximation
The algorithms of [6,9,14] use the following distances when optimizing Eqs. 1 or 2: where q c II (γ (t); γ(t)) is the "convexified" (using the absolute value of the eigenvalues) second fundamental form at point γ(t) and applied in direction γ (t) ∈ T γ(t) M , and P (p j , p) is the set of all paths that connect p j to p on the surface.
As described in Appendix A, for smooth surfaces, it is possible to write: for all λ > 0, for all non-parabolic p j ∈ M there is an open neighborhood V of p j such that ∀p ∈ V : where the notation λ , borrowed from [6], implies tight approximation to any desired degree of accuracy.Note that the exponent 4 above arises from the fact that -nets minimize a form of Euclidean distance between the surface and the approximation, while QEM and sV SA minimize squared Euclidean distance.Equation 6 is valid only for elliptic points p j .For hyperbolic points, D QEM and D sVSA still converge to the same value, but (D II ) 4 is only an upper bound of D QEM and D sVSA , [Note that one could have defined D II using |q II | 1/2 in the integrand of 3, where q II is the second fundamental form.This would make (D II ) 4 be a tight approximation of D QEM and D sVSA everywhere non-parabolic on M , but would no longer be a Riemannian metric.]The distance D II , is too expensive to compute in practice (since each evaluation involves computing a shortest path under the q c II surface metric).In contrast, both D QEM and D sVSA are efficiently computed using only local information at the arguments p j and p.

Normal Approximation
The problem of computing a mesh that approximates the normal field of a surface is considered next.It is noted that for this problem, the normals of the approximating mesh are piecewise constant.However, instead of being inferred from the vertex positions, the normals of the output mesh are optimally assigned to triangles.This distinction is necessary to avoid difficulties like those described in [5,7] that can occur when triangles have large internal angles, even if they have the right limit shape and size.
A similar analysis to that of Sec.2.1 can be made in this case.Here, the two relevant algorithms that are considered are [6,14].They use the following distances to optimize Eqs. 1 and 2 respectively: where q III is the surface's third fundamental form, and n : M → S 2 is the Gauss map.Analogously as proven in Appendix A, it is, for p in an appropriate, small enough neighborhood of a non-parabolic p j :

Behavior
To aid in our discussion, three different kinds of regions on a surface will be considered, and the algorithms under consideration evaluated separately for each.The following distinct types of regions on surfaces are considered: In smooth and non-parabolic regions, it can be shown that, in the limit, the regions of the partitions generated by optimizing D sVSA D nVSA and D QEM have the proper aspect ratio [10,14], which is a necessary condition for optimality for their respective surface or normal approximation problem.It can also be shown that, for everywhere-elliptical surfaces, and in the limit, the method of -nets [6] using D II produces results that are within a constant factor of the globally optimal L ∞ minimizer for the surface surface approximation problem.As discussed in Appendix B, in the limit, primal algorithms such as VSA will need roughly twice as many triangles as compared to dual algorithms such as QEM and -nets.
Near sharp features, these algorithms behave quite differently.In particular, one can see that D QEM measures error "from the viewpoint" of the variable of integration p ∈ M , while D sVSA does so from the viewpoint of the mean p j .As a result, QEM places means at high-curvature points, and thus is suited as a dual algorithm, while sVSA (the surface approximation version of VSA [14]) places them at low -curvature points, making it better suited as a primal algorithm.nVSA also tends to place means at low-curvature points.It is less clear the authors how -nets behaves in this regard.
Parabolic/curved regions are places on a smooth surface near, or at parabolic points where there is significant higher-than-second-order bending (i.e. the surface is not locally well approximated by a quadratic patch), such as near the parabolic line on a torus.Algorithms often need special care in this case.Near parabolic/curved points, an optimization using D QEM undersamples regions near curved parabolic lines.The original QEM algorithm [9] deals with this case by introducing special rules to prevent flips before edge collapses (which strictly-speaking breaks the energy-minimization formulation of Eq. 2).For -nets, an additional isotropic term is added to the distance to cope with such regions.VSA deals with this case, in which the Voronoi cell boundaries are highly curved, by discretizing these boundaries and triangulating the cells finely enough as to avoid undersampling.

Shape Operator Metric for Normal Approximation
An obvious missing piece in this description is an algorithm that converges to Eq. 9 in the limit, but places means at high curvature points away from it, making it most suitable as a dual algorithm.In some sense this algorithm would be to nVSA what QEM is to sVSA.Moreover it can be efficiently computed an optimized (as in Eq. 1 or 2), has high approximation efficiency in the sense of Appendix B, and it avoids heavy undersampling near curved parabolic lines.
To begin, consider the definition D nVSA (p j , p) = n(p j ) − n(p) 2 , which measures normal error from either p j or p, and, similarly as QEM, construct an approximation that only depends on p j but not on any higher-order local information at p j .To do this, a second-order Taylor expansion of D nVSA (p j , p) around p is constructed (note that the zero-th and first order terms vanish): where are the principal curvatures, and {e 1 , e 2 } the principal directions.Note that D SOM , like D QEM and D VSA , can be efficiently computed only from local information at the endpoints, and, as will be shown in Sec. 4, results in an energy of the type of Eq. 1 or 2 that can be efficiently minimized using standard algorithms [13,12].
The SOM algorithm then simply outputs the dual trianglulation of this computed surface partition.
It follows from the fact that this is a dual algorithm, whose distance converges in the limit to that of D nVSA , and from the discussion of Appendix B, that this algorithm has the desired favorable (limit) efficient approximation characteristics when compared to the primal algorithm of [14].It is also shown in Appendix C that this algorithm produces elements that conform to the theoretically optimal limit shape and orientation.
Unlike [6,9,14], curved parabolic regions are dealt with in a natural way, without special consideration, which adds to the simplicity of the algorithm.The side figure illustrates this point, where a mean p j is placed at a parabolic line (red).Because the parabolic line is curved, p j does not lie along the flat direction when viewed from the point of view of nearby points p.An SOM primal-region centered around p j thus cannot grow too much along the parabolic line if the parabolic line curves.
It is possible to see that minimizing Eq. 2 using Eq. 10 has the effect of placing means at points of high-curvature.Consider the closely-related problem of gradient approximation of a scalar function f defined on the plane, and an analogous distance , where H f is the Hessian of f .In an everywhere-isotropic region, D fSOM = k(p) p j − p 2 , which, used for L 2 minimization in a form analogous to Eq. 2 over the plane, is an instance of the weighted k-means problem, which is well-known to place means at points with high weight [1] (high-curvature in this case).The case where H is not isotropic behaves similarly, but the weight can be thought of as directionally-varying.

Implementation of SOM
The energy of Eq. 10 is minimized in a way very similar to the algorithm of [12], which uses a probabilistic seeding of means, followed by a Lloyd relaxation [13] and has theoretical guarantees of closeness to the global optimum.In this work, the probabilistic seeding is simply replaced by iteratively placing means at the surface point with maximum minimum-distance to the current set of means, similarly as the greedy algorithm for computing -nets of [11].This is also similar to the optimization method of [6], except that the seeding is followed by a Lloyd relaxation, and is also similar to [14].The shape operator matrix S of Eq. 10 is estimated using the algorithm of [3].
Once the seed means have been placed, the Lloyd relaxation has two stages.The first creates a distance-dependent Voronoi partition of the surface, and the second computes the new means' locations from the current partition.
To compute a Voronoi partition, all vertices (as opposed to input triangles) are tagged as belonging to some primal Voronoi region, and Voronoi region boundaries are computed by splitting input triangles that have vertices in different regions, as in the side figure below.A Voronoi region is thus not constrained to be a collection of faces, but can have a boundary that cuts across triangles, which may slightly improve accuracy in practice.Also, in this way, Voronoi regions can meet at most at 3-way junctions.These 3-way junctions naturally dualize into triangles.Note that this generalizes to higherdimensions, so that, by construction, it will only output simplicies.
Given a Voronoi partition of the surface, the new means' locations are computed.First, note that the energy of Eq. 2 for the distance D SOM can be written as and so it is quadratic in p j .It is possible to compute the minimizer p j of Eq. 12 by solving a small linear system, but this would return a mean p j which is not constrained to be on the surface.Instead, the constants in equation 12 are computed in a first pass: A j = Vj S(p) 2 dp and b j = Vj S(p) 2 • pdp, for each Voronoi region V j .Then, for each input triangle (or split triangle) inside region V j the barycentric coordinates (u, v) of the minimizer p j of Eq. 12 can be found by solving where R ∈ R 3×2 is some basis of the supporting plane of the triangle.The minimizer may fall outside the triangle, so it is necessary to look for it along triangle edges and vertices as well.The final mean is the minimum over all the minimizers on each triangle, guaranteeing that p j is a point on the surface.Finally, instead of outputting p j directly as a (dual) vertex, a quadric error metric [9] for its associated Voronoi region V j is first computed, and its minimizer along the line passing through p j in direction n(p j ) is output.This small perturbation slightly improves the approximation.

Results
Some surfaces processed by the SOM algorithm are shown in figures 1 through 3.These meshes are compared with those produced by VSA, which are computed by exactly following [14].Note that, unlike SOM, VSA has a free parameter (the precision used to discretize the partition regions's boundaries), which has been tuned to improve VSA's output.These results are also compared with the output of QEM [9].Note that QEM optimizes (RMS) distance from the surface to the approximation, instead of normal error, and therefore the comparison is not strictly relevant; but it is included it as reference.Runtimes for SOM range from 5 sec.(bunny, input: 70k tris.) to 40 sec.(statue, input: 512k tris.), on a single core 2.0GHz CPU.Even though it is not necessarily what is being optimized for in this work, it is possible to consider (L ∞ ) Hausdorff, and RMS error in the sense of surface approximation.Note that, in most cases, QEM produces slightly better approximation of the surface than SOM, and significantly better than VSA.This is expected, as QEM optimizes surface approximation error (RMS error in the figures), while VSA and SOM both optimize normal error instead.Notice that, for smooth surfaces, and using (almost) the same number of triangles, SOM's approximation is appreciably finer than VSA's.On smooth surfaces, the approximation is significantly better for SOM at a given sampling rate.As can be seen in the primal partitions in figures 1 and 2, with an equal triangle budget, SOM is able to partition the surface into smaller regions that capture detail better.Note that the bunny is particularly troublesome for VSA, when compared with SOM, because its bumpy surface produces very curved regions that can output many triangles when their boundaries are discretized by the VSA algorithm.In general, in the above figures, triangles are elongated along the directions of minimum curvature, and tend to show very high anisotropy in places where this is possible: like the ears of the bunny or the statue's arms.Note that our algorithm offers no guarantees in terms of normal flips in the triangulation, which could show up occasionally in sparsely sampled regions.This behavior is similar to VSA and -nets, which also cannot guarantee to be free of flips.Figure 3(d-f) shows a surface composed of roughly flat parts separated by sharp features.On these kinds of surfaces VSA does particularly well, since it operates by locating roughly-flat patches and triangulating them.In particular, the region-triangulation phase of VSA is well-tuned to this problem, since the desired behavior in this case is to triangulate the flat polygons.SOM in this case naturally places means at sharp corners.But its connectivity is guided by the shape operator, which is almost everywhere degenerate here.This case is dealt with by computing the final mesh connectivity using a modified shape operator, which is set to a very high (isotropic) value in flat regions, effectively simulating a flat-polygon triangulation step (similar to the constrained Delaunay triangulation used in [14]).

Numerical validation
Unlike for surface approximation, there is, as far as the authors are aware, no standard way of measuring the surface normal approximation on a surface.If there was, away from the limit regime, a well-defined one-to-one correspondence between points on the surface and points on the approximation, then it would be possible to compute the (L 2 ) approximation error by integrating the distance between corresponding normals over the surface.However, this correspondence is not available.To analyze approximation error, the very closely-related problem of approximation of the gradient of a height field over  the plane is considered [it has optimal limit aspect ratio ξ 1 /ξ 2 , where ξ i are the eigenvalues of the heigh field Hessian [4,5]].Because both VSA and SOM only look at normals and shape operators, it is possible to naturally adapt both to the gradient approximation case by measuring distances between gradients, as opposed to normals, by computing a Hessian of the height field at each point, instead of a shape operator.Both algorithms must also be extended to force them to conform to the boundary of the domain.There is however, to our knowledge, no equivalent natural generalization of QEM [9] to the heigh field approximation case.Once again, the tunable parameter in VSA has been adjusted to the best results obtained.The input is a surface that is finely scan-converted on a squared piece of the plane (figure 5 topright corners.)Planar meshes obtained this way are shown in figure 4, while figure 5 shows the corresponding error plots for these two inputs, at several approximation levels.Notice how the mesh approximating the mask in 4.b more closely matches the features of the input than 4.a, even though both have the same number of triangles.The difference in RMS error is not always large for a fixed number of triangles, though it is significant.If, alternatively, an RMS error level is fixed, and the VSA and SOM approximations with that error are considered, it can be noted that the SOM mesh has significantly fewer triangles.

Summary and Conclusion
This work begins by placing some established algorithms for surface approximation into a common framework.From this analysis, it becomes apparent that a dual variational algorithm for surface normal approximation was previously missing.Such algorithm is introduced next, and its limit behavior shown to conform with the theoretical asymptotic aspect ratio (Appendix C).
It is further argued that this dual algorithm has several advantages over primal variational algorithms for surface normal approximation (such as VSA).
In particular, the limit approximation efficiency is discussed in Appendix B, which is shown to be approximately 1.75 times higher for a dual algorithm with the same (asymptotically optimal) limit aspect ratio.The approximation results of the proposed algorithm and established ones are also compared on practical data sets.While the primal VSA is still preferable for piecewise flat surfaces, where it successfully splits them into flat regions which are then triangulated, for general curved surfaces, the algorithm introduced in this paper is shown to perform better.This is further shown on quantitatively experiments, which are carried out on the very closely related problem of gradient approximation over the plane.

Fig. 4 .
Fig. 4. Gradient approximation using nVSA and our algorithm.Primal (top), algorithm output (middle), height field approximation (bottom).Red marks in the primal are locations of vertices in the dual.In both (a) and (b), the mask (left side) is approximated with 356 triangles, and the bunny (right side) with 468 triangles.