Consistently Orienting Facets in Polygon Meshes byMinimizing the Dirichlet Energy of Generalized Winding Numbers
Kenshi Takayama

Alec Jacobson

Ladislav Kavan

Olga Sorkine-Hornung

ETH Zurich

ETH Zurich

University of Pennsylvania

ETH Zurich

Abstract
Jacobson et al. [JKSH13] hypothesized that the local coherency of the generalized winding number
function could be used to correctly determine consistent facet orientations in polygon meshes. We reporton an approach to consistently orienting facets in polygon meshes by minimizing the Dirichlet energy ofgeneralized winding numbers. While the energy can be concisely formulated and efficiently computed, wefound that this approach is fundamentally flawed and is unfortunately not applicable for most handmade
meshes shared on popular mesh repositories such as Google 3D Warehouse.

1Introduction

The robust inside-outside segmentation method based on the generalized winding number proposed by
Jacobson et al. [JKSH13] assumes that the polygonal facets in the input mesh are consistently oriented.
If this assumption does not hold, which is often the case for hand-made meshes shared on the Internet(e.g., Google 3D Warehouse), it is necessary to first correct the facet orientations. In this technical report,
we present an approach to correcting inconsistent facet orientations based on the Dirichlet energy of thegeneralized winding number. Our motivation came from the fact that if the mesh is watertight and its facets
are consistently oriented, the generalized winding number is piecewise constant, and hence its Dirichlet
energy is zero. Our intuition was that for a mesh with various imperfections such as holes and intersections,more consistently oriented facets would result in a smaller Dirichlet energy of the corresponding generalized
winding number.

2Method

2.1

Patch extraction

As a preprocessing step, facets that are topologically connected and forming a two-manifold surface aregrouped together as a patch, similar to Borodin et al.’s method [BZK04]. If the mesh has non-manifoldedges, the mesh is split around these edges, and multiple patches are generated. All facets belonging to
the same patch share the same decision about whether they should be flipped or not in the output. Givenn patches, our goal is to seek an assignment for n binary variables representing whether the i-th patch Pi
should be flipped or not (1 ≤ i ≤ n).

2.2

Generalized winding numbers

Let us first review the concept of generalized winding numbers [JKSH13]. The generalized winding number(hereafter we drop the word “generalized”) at a point p ∈ R3 is defined asw(p) =1 X
θτ (p) ,4π τ ∈M(1)

where M is the input mesh (denoted as a set of facets), τ ∈ M is a facet in M, and θτ is the signed solid
angle subtended by τ viewed from p. Note that the sign of the angle is defined to be positive for points in
the half space corresponding to the back side of a facet. In particular, the winding number of points insideand outside a watertight mesh is 1 and 0, respectively.

2.3

Dirichlet energy of the winding number

Since the winding number of a watertight and consistently oriented mesh is piecewise constant and hence
its Dirichlet energy is zero, our intuition is that the Dirichlet energy of the winding number for an imperfectmesh should decrease as its facet orientations become more consistent. We give a mathematical formulationof the Dirichlet energy in the following.

1

For each i-th patch Pi ⊆ M (1 ≤ i ≤ n), we want to determine a binary variable si ∈ {+1, −1} called
sign, where si = −1 means that all facets in Pi should be flipped in the output. Now we can define thewinding number with respect to the given sign s = (s1 , · · · , sn ) asws (p) =n
1 X Xsiθτ (p) .
4π i=1 τ ∈P(2)i

Note that the orientation of each triangle τ is fixed to the original configuration when computing θτ . TheDirichlet energy of ws denoted as E(s) is a function of s and is defined as the integral of the gradientmagnitude squared of ws over the volumetric domain Ω = R3 \ M:
ZE(s) =k∇ws (p)k2 dV .
(3)p∈Ω

Reorganizing the integrand as2k∇ws (p)k2

===n
1 X X
si
∇θτ (p)4π i=1 τ ∈Pi

 
nn
XX
1 X X∇θµ (p)
∇θτ (p) · sjsi
16π 2 i=1 τ ∈P
j=1
µ∈P
ji

 nn
XX1 XX∇θτ (p) · si sj ∇θµ (p) ,16π 2 i=1 j=1
τ ∈Pµ∈Pi(4)

(5)

(6)

j

the energy can be written concisely in quadratic form as
E(s) = s> Q s ,(7)where Q is an n-by-n matrix with elements
Z1 X X∇θτ (p) · ∇θµ (p) dV .
Qij =16π 2 τ ∈P µ∈P p∈Ωi(8)jQ is fixed per model independently of s. Our goal is to find a set of facet orientations s∗ that gives thesmallest Dirichlet energy:s∗ = arg min E(s).
(9)
s∈{−1,+1}nSince s2i = 1, the energy can be further reorganized as
E(s) = 2n−1
XnX

si sj Qij +i=1 j=i+1

n
XQii .

(10)i=1| {z }
const.

Ignoring the diagonals of Q, we can also considerẼ(s) =n−1
XnX

si sj Qij

(11)

i=1 j=i+1

as the energy instead.

2.4

Computing the Dirichlet energy

We need to compute Qij in order to evaluate E(s). A straightforward approach would be to tessellate
the finite volume around the input mesh with sufficiently large volume and fine resolution to capture thevariations of the winding numbers and then numerically compute the integral over this tessellation. However,this approach can be computationally expensive, especially when the mesh contains many small details anda fine tessellation is needed for achieving sufficient accuracy.
Fortunately, since the winding number function is harmonic (i.e., ∆w(p) = 0) [JKSH13], we can convert
the volumetric integral over the domain Ω into a surface integral over the domain boundary ∂Ω using Green’s
first identity. For a once differentiable scalar function ψ and a twice differentiable scalar function φ, Green’s
first identity statesZZZ∇ψ · ∇φ dV =ψ∆φ dV +ΩΩψ (∇φ · n) dS(12)∂Ω

where n is an outgoing normal at the domain boundary (the positional variable p is omitted here for brevity).
By substituting ψ = θτ and φ = θµ and recalling that winding numbers are harmonic, we getZ
Z: 0Z
θ
∆θ
dV+
∇θ·∇θdV=θτ (∇θµ · n)dS ,(13)
τ
µτ
µ

ΩΩ∂Ω

2

which shows that the volumetric integral in (8) is now converted to a surface integral.
To compute this surface integral, we make the following notes on the domain and its boundary. Thedomain Ω consists of the entire 3D space R3 with a hole of infinitesimally thin slab carved out by each facetτ ∈ M, and the domain boundary ∂Ω consists of such holes. Each hole corresponding to τ has its front andback sides denoted by τ + and τ − , respectively. Note that the domain’s outgoing normals on τ + and τ − arepointing toward the hole and are opposite to each other, see Figure 1. Now the domain boundary can bewritten as[ +∂Ω =τ ∪ τ −.
(14)
τ ∈M

For a facet τ ∈ M and a point on it p ∈ τ , let us denote the corresponding point on the front and backsides as p+ ∈ τ + and p− ∈ τ − , respectively. Then, it is easy to see that the following properties hold:θτ (p+ )=−θτ (p− ) = −2π(15)+=∇θτ (p− ) ,(16)∇θτ (p )where the latter is due to the derivative of the arctangent. Note that θτ and ∇θτ are undefined along theedges of τ ; we simply ignore the edge parts (i.e., parts connecting τ + and τ − ) in the surface integral, since
the measure of all the edges vanishes (as opposed to the measure of the facets’ interiors).
Normal on

Original facet
normalNormal on

Figure 1: Notations for the domain boundary.
For any other third facet γ 6= τ , θτ and ∇θµ take the same value on γ + and γ − . Therefore,ZZZθτ (∇θµ · n)dSθτ (∇θµ · n)dS +θτ (∇θµ · n)dS =γ−γ+γ + ∪γ −ZZθτ (∇θµ · (−n))dSθτ (∇θµ · n)dS +=(18)γ+γ+=(17)0(19)

where in (18) we replaced the integral over γ − with an integral over γ + with a negated surface normal.
Therefore, the integral over ∂Ω is reduced to an integral over τ + ∪ τ − , i.e.,
ZX Zθτ (∇θµ · n)dS =θτ (∇θµ · n)dS
(20)
∂Ωγ∈M

γ + ∪γ −Zθτ (∇θµ · n)dS .

=(21)

τ + ∪τ −Using the properties above, we arrive at:Zθτ (∇θµ · n)dS =τ + ∪τ −==ZZ(−2π)(∇θµ · n)dS +(2π)(∇θµ · n)dS−τ+ZZτ
−2π∇θµ · n dS + 2π
∇θµ · (−n) dS+τ+Zτ
−4π∇θµ · n dS .

(22)
(23)
(24)

τ+In summary, Qij is computed as
Qij =X XF (τ, µ)(25)

τ ∈Pi µ∈Pj

whereF (τ, µ) = −14πZ∇θµ (p) · n(p) dS .
p∈τ +3(26)2.5Convergence of the Dirichlet energy

By deriving an analytic form of F (τ, µ) using Maple, we found that it evaluates to infinity (either positiveor negative depending on relative orientation of τ and µ) when some edges of τ and µ coincide geometrically.
An intuitive explanation for this is that the magnitude of ∇θµ (p) becomes infinitely large when p ∈ τ + goesinfinitesimally close to the edges of τ that are geometrically coinciding with the edges of µ. In particular,
the diagonal entries Qii always evaluate to positive infinity. Although we can ignore the diagonals using(11), off-diagonal entries Qij can still evaluate to infinity if some boundary edges of Pi and Pj coincide
geometrically. In our experiments, we numerically evaluate the energy with a sufficiently fine discretizationof ∂Ω where the finite entries in Q are converged.

2.6

Finding the minimizer of the Dirichlet energy

The minimization problem in (9) is referred to as Quadratic Unconstrained Binary Optimization (QUBO).
QUBO is NP-hard, and various heuristic methods exist [BHT07]. In our experiments, we take a brute
force approach and exhaustively search for the optimal solution, so that we can examine the nature of ourapproach with small toy examples.

3Results

We implemented our approach in 3D as well as in 2D where the input consists of open or closed polylines(still referred to as “patches” in the following) and the same formulation applies. Although here we show
our experiments only in 2D for the ease of visualization, the same argument also holds in 3D. Figure 2 showscases where the input corresponds to the boundary of a solid shape with some missing pieces. In such cases,
our algorithm produces the expected results, i.e., it orients the facets consistently and correctly identifies
the shape’s interior. On the other hand, Figure 3 shows a case where our approach is fundamentally flawed:
the input has some small patches attached to larger patches, and those small patches are oriented inside-out.
The fundamental reason for this is that the Dirichlet energy only measures the variation of the function
while ignoring the sign, and hence there is no penalty for orienting facets inside-out.
Optimized

11.3

Optimized

6.83

16.4

6.94

Figure 2: Examples where our approach produces reasonable results. The input polylines are represented bywhite lines with their orientations indicated by yellow arrows. The number below each image is the Dirichlet
energy of the corresponding winding number function. The color maps represent the associated winding numbers, where higher winding numbers correspond to warmer colors (the actual mapping is slightly adjusted foreach image to better visualize the distribution).

Optimized

9.28

9.01

Figure 3: An example representing a torso with two small objects attached to the shoulder and the side. Our
approach results in those small attachments oriented inside-out.

4For an intuitive explanation of this behavior, it is useful to consider the closure of a patch P denotedas P̄ [JKSH13], where the union of P and P̄ bounds a finite volumetric domain ΩP (i.e., ∂ΩP = P ∪ P̄).
Denoting the winding number defined by a patch P as wP (p), the following holds:
1, p ∈ ΩPwP (p) + wP̄ (p) =.
(27)
0, otherwise
Since the right hand side of (27) is piecewise-constant, the following also holds:∇wP (p) = −∇wP̄ (p) .

(28)

Let us now consider a simple case with two patches P1 and P2 as shown in Figure 4. We can decompose
the gradient and apply the above observation as:∇wP1 ∪P2 (p)=∇wP1 (p) + ∇wP2 (p)

(29)

=

−∇wP¯1 (p) − ∇wP¯2 (p) .

(30)

It is clear in the figure that the current configuration of P¯1 and P¯2 yields smaller Dirichlet energy than theother configuration, and hence P2 is oriented this way in our approach.

=+=−−=−Figure 4: An example showing the reason of why small attachments are oriented inside-out in our approach.
Note again that the mapping between the winding number and the color is adjusted in each image to better
visualize the function.

4Conclusion and future work

We conclude that the presented approach has a fundamental flaw and is not applicable for most handmade
meshes. One direction for further exploration is to apply a similar approach to surface reconstruction froman unoriented point cloud [MDGD+ 10, LW10], where only the positions of the points are given and finding
good normals for the points is critical. In this case, the unknowns are the continuous normal directions,potentially making the problem easier. A key question, however, would be how to generalize the definitionof winding numbers to point clouds.

Acknowledgements
We are grateful to Peter Kaufmann, Wenzel Jakob, Nobuyuki Umetani, Melina Skouras, Ilya Baran, andRyan Schmidt for their advice and feedback. This work was supported in part by the ERC grant iModel

5

(StG-2012-306877), by an SNF award 200021 137879, by an NSF grant IIS-1350330 and by a gift from Adobe
Research. Kenshi Takayama’s work was funded by JSPS Postdoctoral Fellowships for Research Abroad. Alec
Jacobson’s work was supported by an Intel Doctoral Fellowship.

References[BHT07]

Endre Boros, PeterL. Hammer, and Gabriel Tavares. Local search heuristics for quadratic
unconstrained binary optimization (qubo). Journal of Heuristics, 13(2):99–132, 2007.

[BZK04]

Pavel Borodin, Gabriel Zachmann, and Reinhard Klein. Consistent normal orientation forpolygonal meshes. In Proceedings of the Computer Graphics International, CGI ’04, pages
18–25, 2004.

[JKSH13]

Alec Jacobson, Ladislav Kavan, and Olga Sorkine-Hornung. Robust inside-outside segmentation using generalized winding numbers. ACM Trans. Graph., 32(4):33:1–33:12, 2013.

[LW10]Shengjun Liu and Charlie C.L. Wang. Orienting unorganized points for surface reconstruction.
Computers & Graphics, 34(3):209–218, 2010.

[MDGD+ 10] Patrick Mullen, Fernando De Goes, Mathieu Desbrun, David Cohen-Steiner, and Pierre Alliez.
Signing the unsigned: Robust surface reconstruction from raw pointsets. Computer Graphics
Forum, 29(5):1733–1741, 2010.

6