Consistent Point Orientation for Manifold Surfaces via Boundary Integration

Abstract.

This paper introduces a new approach for generating globally consistent normals for point clouds sampled from manifold surfaces. Given that the generalized winding number (GWN) field generated by a point cloud with globally consistent normals is a solution to a PDE with jump boundary conditions and possesses harmonic properties, and the Dirichlet energy of the GWN field can be defined as an integral over the boundary surface, we formulate a boundary energy derived from the Dirichlet energy of the GWN. Taking as input a point cloud with randomly oriented normals, we optimize this energy to restore the global harmonicity of the GWN field, thereby recovering the globally consistent normals. Experiments show that our method outperforms state-of-the-art approaches, exhibiting enhanced robustness to noise, outliers, complex topologies, and thin structures. Our code can be found at https://github.com/liuweizhou319/BIM.

1. Introduction

Point clouds with globally consistent normals have numerous downstream applications, such as point set upsampling (Huang et al., 2013; Feng et al., 2022), point cloud filtering (Zhang et al., 2020; Lu et al., 2022), surface reconstruction (Kazhdan et al., 2006; Huang et al., 2009; Calakli and Taubin, 2011; Kazhdan and Hoppe, 2013; Lu et al., 2018; Mescheder et al., 2019; Sellán and Jacobson, 2022, 2023), and segmentation (Khaloo and Lattanzi, 2017). As a fundamental problem in computer graphics, it remains challenging due to the complexity of object geometry, topology, and the influence of noise.

Most existing methods (Hoppe et al., 1992; Cazals and Pouget, 2003; Pauly et al., 2003; König and Gumhold, 2009; Schertler et al., 2017; Metzer et al., 2021) estimate normals perpendicular to the point cloud surface and then use propagation to flip inconsistent normals. As pointed out by Li et al. (2023a), propagation strategies are affected by the direction of distribution of the unoriented normal vectors and often lead to undesirable results. Moreover, neighborhood-based normal estimation methods tend to perform poorly when dealing with noise and geometries with small reach. In the past few years, there has been significant advancement in the field of globally consistent normal estimation. iPSR (Hou et al., 2022) reconstructs surfaces from unoriented point clouds by iteratively refining the normals obtained from the PSR solver (Kazhdan et al., 2006). While effective for dense point clouds, iPSR encounters difficulties with sparse inputs. In sparsely sampled regions, the potential generation of erroneous normals by iPSR can result in disconnected components in the reconstructed surface. There are also deep learning-based methods for point normal estimation or orientation (Guerrero et al., 2018; Ben-Shabat et al., 2019; Zhu et al., 2021; Li et al., 2023a, b, c, d). While these methods enhance the precision of predictions, they typically do not reverse the normals when the prediction of normal orientation is inconsistent.

The generalized winding number (GWN) (Jacobson et al., 2013) has proven to be a powerful tool for robust segmentation of 3D surfaces into distinct inside and outside regions. Its application extends to the orientation of triangle meshes and point clouds, showcasing its versatility in the field of geometric processing. Building on this foundational work, Takayama et al. (2014) proposed to orient triangle soups by minimizing the Dirichlet energy associated with the GWN field. Barill et al. (2018) developed a point cloud-based approach to estimating the GWN, linking point cloud normals to the generalized winding number. More recently, PGR (Lin et al., 2022) and GCNO (Xu et al., 2023) adopted point-based GWN for orienting point clouds. They regularize the values of the GWN both inside and outside the point cloud. Despite these advancements, these methods are usually limited to small-scale models due to their high computational cost. Additionally, they are not effective in handling models with complex topology and often fail to robustly address issues posed by thin structures and noise.

Our work aims at developing a robust approach for computing globally consistent normal orientations in point clouds sampled from manifold surfaces. We identify harmonicity and jump boundary conditions as two key factors in point cloud orientation. Leveraging the fact that the Dirichlet energy of the GWN field can be expressed as an integral over the boundary (Takayama et al., 2014), we introduce a new boundary energy formulation. In this framework, the normals associated with each point are variables. This boundary energy coincides with the Dirichlet energy only when the normals are globally consistent. By imposing constraints on the range of the GWN values, we observe an inverse correlation between boundary energy and Dirichlet energy: the boundary energy is small for randomly oriented normals, and becomes large when normals are globally consistent. This observation motivates us to maximize the boundary energy from initially random normals. Employing an efficient L-BFGS solver, we iteratively enhance the global harmonicity of the generalized winding number field, eventually yielding globally consistent normals. Extensive experiments demonstrate the effectiveness of our method, achieving high-quality orientation and reconstruction results.

2. Related Work

The problem of global point orientation dates back to the 1990s, in the seminal work of Hoppe et al. (1992). The authors calculated initial normals utilizing principal component analysis, and subsequently corrected inconsistent normals through propagation. This pioneering work has inspired many follow-up works, including (König and Gumhold, 2009; Jakob et al., 2019; Levin, 1998; Mitra and Nguyen, 2003; Alliez et al., 2007; Boulch and Marlet, 2012; Xu et al., 2022), among others. Although these methods perform well with simple models, their performance deteriorates in the presence of complex geometries and data imperfections, especially in the presence of noise and outliers.

Reconstructing 3D surfaces from unoriented point clouds is closely related to the problem of point orientation. A popular approach for 3D reconstruction involves using an implicit function to represent the target surface, where the gradients of this function at a particular level set serve as the point normals. Current approaches predict globally consistent normals by using signed/unsigned distance fields or by considering point clouds as dipoles to generate generalized winding number fields. Mullen et al. (2010) computed an unsigned distance field for the input point cloud, subsequently determining the sign by minimizing the Dirichlet energy. NGLO (Li et al., 2023a) enhanced normal accuracy by optimizing gradients based on the fitted distance field. NeuralGF (Li et al., 2023b) introduced a novel paradigm of learning gradients for fitting a signed distance field with neural implicit functions, enabling unsupervised point cloud normal estimation. The Dipole method (Metzer et al., 2021) maximizes the electric field gradient of the point cloud using a greedy strategy, aiming to approximate the generalized winding number. iPSR (Hou et al., 2022) iteratively refined the surface by feeding the normals obtained from the previous iteration into the PSR (Kazhdan and Hoppe, 2013) solver, gradually improving the surface quality. PGR (Lin et al., 2022) treats normals and surface element areas as unknown parameters, interpreting the indicator function as a member of a parametric function space using the Gauss formula, and recovering the normals by estimating the values of the indicator function. GCNO (Xu et al., 2023) incorporated the prior knowledge of the generalized winding number field as an approximation to the indicator function for point clouds with globally consistent normals. Despite significant advancements in the field of 3D reconstruction, existing methods for reconstructing 3D surfaces from unoriented point clouds still suffer from various issues. For example, NGLO’s dependency on a training set limits its ability to robustly handle diverse types of data and noise levels. Dipole may generate conflicting normals at thin structures. iPSR can lead to disconnected components in cases of sparse point clouds. PGR lacks robustness to noise and demands high computational resources. GCNO’s effectiveness is contingent upon an even distribution of Voronoi vertices both inside and outside the point cloud–a condition that, when unmet, significantly degrades the quality of the results.

PDEs with jump boundary conditions are widely used in surface reconstruction, rendering, and fluid simulation. In surface reconstruction, jump boundary conditions are used to infer the implicit function from the point cloud with consistent normals. The indicator function obtained from the PSR (Kazhdan et al., 2006; Kazhdan and Hoppe, 2013) solver effectively satisfies the jump boundary conditions. The generalized winding number generated by point clouds with globally consistent normals (Barill et al., 2018) can be regarded as an approximate solution to the Poisson equation with jump boundary conditions. The surface can be reconstructed by extracting the iso-surfaces from the generalized winding number. In rendering, Poisson vector graphics (Hou et al., 2020) and diffusion curves (Orzan et al., 2008) treat the input strokes as jump boundary conditions and accomplish color rendering by solving Laplace’s equation. In fluid simulation (Leal, 2007), jump boundary conditions are used to represent the physical properties at the interface between two fluids. Recently, Feng et al. (2023) proposed a method for computing the winding number for discrete surfaces. They approached the problem by considering the relationship between the winding number and harmonic functions and solved it by solving Laplace’s equation with jump boundary conditions.

3. Preliminaries

Given a solid object $\Omega\subset\mathbb{R}^{3}$, we denote its boundary surface by $\partial\Omega$. Let $\mathcal{P}=\{\mathbf{p}_{i}|\mathbf{p}_{i}\in\partial\Omega\}_{i=1}^{n}$ be a point cloud sampling the manifold surface. For each point $\mathbf{p}_{i}$, we denote by $\hat{\mathbf{n}}_{i}$ its outward unit normal. We denote the inner product in $\mathbb{R}^{3}$ by $\langle\cdot,\cdot\rangle$.

3.1. Poisson’s Equation with Jump Boundary Conditions

Poisson’s equation with jump boundary conditions finds extensive applications across physics, scientific computing, simulation, and geometric processing. The equation is formulated as:

(1) $$\Delta u(\mathbf{x})=f(\mathbf{x}),\;\;\mathbf{x}\in\;\mathbb{R}^{3}\setminus\partial\Omega,$$

subject to the Dirichlet boundary condition

(2) $$u^{+}(\mathbf{x})-u^{-}(\mathbf{x})=k,\;\;\mathbf{x}\in\;\partial\Omega,$$

and the Neumann boundary condition

(3) $$\partial u^{+}(\mathbf{x})/\partial\mathbf{n}_{\mathbf{x}}=\partial u^{-}(\mathbf{x})/\partial\mathbf{n}_{\mathbf{x}},\;\;\mathbf{x}\in\;\partial\Omega,$$

where $u:\mathbb{R}^{3}\to\mathbb{R}$ is the sought solution, $f:\mathbb{R}^{3}\to\mathbb{R}$ represents a source term, and $k$ is a constant. The term $\mathbf{n}_{\mathbf{x}}$ denotes the globally consistent outward normal at $\mathbf{x}$. We also define exterior and interior boundary values as the limiting values of $u$ at the boundary from the exterior and interior, respectively. That is, $u^{\pm}(\mathbf{x})=\lim_{\epsilon\to 0}u(\mathbf{x}\pm\epsilon\mathbf{n}_{\mathbf{x}})$ for $x\in\partial\Omega$. When $f=0$, the solution $u$ is harmonic, and in this context, $u$ also minimizes the Dirichlet energy functional $\int_{\mathbb{R}^{3}\setminus\partial\Omega}\|\nabla u\|^{2}$.

3.2. Boundary Integral Equation

The boundary element method is a powerful tool for solving partial differential equations. BEM transforms the differential equation within the computational domain into a boundary integral equation. The boundary is then discretized into elements, allowing the discretization of the boundary integral equation as a dense system of linear equations without requiring explicit discretization of the computational domain.

Equation (1) can be expressed in the form of a Boundary Integral Equation (Costabel, 1987) using the Poisson kernel and Green’s function as follows:

(4) $$u(\mathbf{x})=\int_{\mathbb{R}^{3}\setminus\partial\Omega}G(\mathbf{x},\mathbf{y})f(\mathbf{y})\mathrm{d}\mathbf{y}$$

$$+ \int_{\partial\Omega}P(\mathbf{x},\mathbf{z})\left[u^{+}(\mathbf{z})-u^{-}(\mathbf{z})\right]-G(\mathbf{x},\mathbf{z})\left[\frac{\partial u^{+}(\mathbf{z})}{\partial\mathbf{n}_{\mathbf{z}}}-\frac{\partial u^{-}(\mathbf{z})}{\partial\mathbf{n}_{\mathbf{z}}}\right]\mathrm{d}\mathbf{z}$$

$$= \int_{\Omega}G(\mathbf{x},\mathbf{y})f(\mathbf{y})\mathrm{d}\mathbf{y}+\int_{\partial\Omega}P(\mathbf{x},\mathbf{z})\left[u^{+}(\mathbf{z})-u^{-}(\mathbf{z})\right]\mathrm{d}\mathbf{z}$$

where $G(\mathbf{x},\mathbf{y})$ represents Green’s function satisfying $\Delta G(\mathbf{x},\mathbf{y})=\delta_{\mathbf{y}}(\mathbf{x})$. Here $\delta_{\mathbf{y}}(\mathbf{x})$ is Dirac delta function and $P(\mathbf{x},\mathbf{y})$ denotes the Poisson kernel. The final term in the second integral vanishes due to the Neumann boundary condition.

3.3. Generalized Winding Number

The generalized winding number is the number of times a point in space is enclosed by a surface. As shown in (Jacobson et al., 2013), the generalized winding number $w$, generated by globally consistent outward normals, satisfies Laplace’s equation. This condition represents a special case of Equation (1) when $k=-1$ and $f\equiv 0$. Consequently, we can express $w$ using Equation (4) as

(5) $$w(\mathbf{x})=\int_{\partial\Omega}P(\mathbf{x},\mathbf{z})d\mathbf{z}.$$

The Poisson kernel $P(\mathbf{x},\mathbf{y})$, which is the directional derivative of Green’s function in the direction $\mathbf{v}_{\mathbf{y}}$, is given by:

$$P(\mathbf{x},\mathbf{y})=\frac{\partial G(\mathbf{x},\mathbf{y})}{\partial\mathbf{v}_{\mathbf{y}}}=\frac{1}{4\pi}\frac{\langle\mathbf{n}_{\mathbf{y}},\mathbf{x}-\mathbf{y}\rangle}{\|\mathbf{x}-\mathbf{y}\|^{3}}.$$

Given an oriented point cloud $\mathcal{P}$, where each point $\mathbf{p}_{i}$ is associated with a globally consistent outward normal $\mathbf{n}_{i}$, Barill et al. (2018) approximates the GWN for an arbitrary query point $\mathbf{q}\in\mathbb{R}^{3}$ by

(6) $$w(\mathbf{q}):=\sum\limits_{i=1}^{n}{a_{i}}\frac{1}{4\pi}\frac{\langle\mathbf{n}_{i},\mathbf{p}_{i}-\mathbf{q}\rangle}{\|\mathbf{p}_{i}-\mathbf{q}\|^{3}}=\sum\limits_{i=1}^{n}a_{i}P(\mathbf{q},\mathbf{p}_{i}).$$

Here $a_{i}$ is the geodesic Voronoi area of point $\mathbf{p}_{i}$, representing the contribution of point $\mathbf{p}_{i}$ to the surface integral. Since computing geodesic Voronoi diagram on discrete surfaces is time consuming, Barill et al. (2018) suggested calculating $a_{i}$ as a 2D Voronoi area by projecting $\mathbf{p}_{i}$ and its $k$-nearest neighbors to the tangent plane of $\mathbf{p}_{i}$. Jacobson et al. (2013) showed that the GWN for points that are away from the surface is harmonic, satisfying $\Delta w=0$. Furthermore, as observed by Xu et al. (2023), the GWN produced by a point cloud with randomly oriented normals is nearly zero everywhere (see Figure 2).

4. Method

Our method takes a point cloud $\mathcal{P}=\{\mathbf{p}_{i}\}_{i=1}^{n}$, which is sampled from a closed, orientable manifold surface $\partial\Omega$ as input. It begins by assigning a random normal to each point in the cloud. The core of our method lies in the regularization of the generalized winding number $w$ by optimizing an objective function related to the Dirichlet energy of $w$. After orientation, we use sPSR (Kazhdan and Hoppe, 2013) to obtain the watertight surface.

Notations

Throughout the paper, we denote by $\hat{\mathbf{n}}=\{\hat{\mathbf{n}}_{i}\}_{i=1}^{n}$ the set of globally consistent outward unit normal of the surface $\partial\Omega$, and by $\mathbf{n}=\{\mathbf{n}_{i}\}_{i=1}^{n}$ an arbitrary unit normal field. Although the GWN for the input point cloud is computed in a discrete manner, we treat the GWN as a continuous field for the sake of clarity in our discussion. We define $\Omega^{+}$ (resp. $\Omega^{-}$) to be the exterior (resp. interior) of the surface. The GWN for points $\mathbf{x}\in\partial\Omega$ is defined as $w^{\pm}=\lim_{\epsilon\to 0}w(\mathbf{x}\pm\epsilon\hat{\mathbf{n}})$. We define the consistently oriented normals, $\hat{\mathbf{n}}^{+}$ and $\hat{\mathbf{n}}^{-}$, respectively, for the boundary of $\Omega^{+}$ and $\Omega^{-}$, where $\hat{\mathbf{n}}^{+}=-\hat{\mathbf{n}}$ and $\hat{\mathbf{n}}^{-}=\hat{\mathbf{n}}$.

4.1. Boundary Energy

Since the GWN field associated with a globally consistent orientation is a harmonic function away from the surface and thereby minimizes the Dirichlet energy, Takayama et al. (2014) proposed the minimization of the Dirichlet energy in approximating the GWN from polygonal soups. The Dirichlet energy is defined as an integral over the whole space, excluding the boundary surface $\partial\Omega$, and is given by:

$$\min\int_{\mathbb{R}^{3}\setminus\partial\Omega}\|\nabla w\|^{2}:=\int_{\Omega^{+}}\|\nabla w\|^{2}+\int_{\Omega^{-}}\|\nabla w\|^{2}.$$

A straightforward approach to compute the Dirichlet energy involves subdividing the finite volume around $\partial\Omega$ into sufficiently small elements to accurately capture the variations of the GWN. Subsequently, numerical integration is performed on this discretized space to calculate the Dirichlet energy. In order to calculate integrals accurately, a high-resolution and high-quality tetrahedral mesh is typically necessary, particularly when $\partial\Omega$ has complex geometry and topology. An alternative solution, as proposed by Takayama et al. (2014), is to convert the volume integral $\int_{\mathbb{R}^{3}\setminus\partial\Omega}\|\nabla w\|^{2}$ into boundary integrals using Green’s first identity

(7) $$\int_{\mathbb{R}^{3}\setminus\partial\Omega}\|\nabla w\|^{2}=\int_{\partial\Omega}(w^{-}-w^{+})\nabla_{\hat{\mathbf{n}}}w^{-},$$

where $\nabla_{\hat{\mathbf{n}}}$ denotes the directional derivative along normal direction $\hat{\mathbf{n}}$.

We extend the Dirichlet energy associated with globally consistent normals $\hat{\mathbf{n}}$ to an arbitrary normal vector field $\mathbf{n}$, and define the boundary energy as

(8) $$f(\mathbf{n}):=\int_{\partial\Omega}(w^{-}-w^{+})\nabla_{\mathbf{n}}w^{-}.$$

Given an arbitrary normal field $\bf n$, we decompose each normal vector $\mathbf{n}$ into two components relative to the globally consistent outward normal $\hat{\mathbf{n}}$ – a tangential component $\mathbf{n}_{t}$ and a normal component $\mathbf{n}_{n}$, that is, $\mathbf{n}=\mathbf{n}_{n}+\mathbf{n}_{t}$. Leveraging the linearity of directional derivatives, we can express the boundary energy for $\mathbf{n}$ as a sum of its components: $f(\mathbf{n})=f(\mathbf{n}_{n})+f(\mathbf{n}_{t})$.

Empirical analysis leads to two key observations: 1) When normals are oriented randomly, both the tangential and normal components of the boundary energy approach zero, yet the Dirichlet energy remains high (Figure 2, bottom). 2) Conversely, when normals are globally consistent, the tangential component remains negligible $f(\mathbf{n}_{t})\approx 0$, whereas the normal component is large, leading to $f(\mathbf{n})\approx f(\mathbf{n}_{n})$ (Figure 2, top). Under this condition, the Dirichlet energy is reduced, aligning with the boundary energy. These observations highlight the competing trends in Dirichlet and boundary energies during the process of transforming random normals to consistent normals. This serves as the inspiration for our approach to directly maximize the boundary energy $f(\mathbf{n})$. As a result, our optimization becomes

(9) $$\max f(\mathbf{n}):=\int_{\partial\Omega}\left(w^{-}-w^{+}\right)\nabla_{\mathbf{n}}w^{-}$$

Remark 1

It is important to note that the boundary energy $f(\mathbf{n})$, despite a superficial similarity to the Dirichlet energy, represents a fundamentally different concept. These two forms of energy coincide only when the normal associated to each point is perpendicular to the surface. For an arbitrary normal field $\bf n$, the two energies diverge. Furthermore, computing the Dirichlet energy for an arbitrary normal field $\bf n$ involves discretizing the space for volume integration - a computationally intensive task. Thus, minimizing the Dirichlet energy with randomly initialized normals is theoretically possible, but it is impractical for real-world applications. Our method overcomes this obstacle by converting volume integration to boundary integration when dealing with random normals, significantly reducing computational demands and eliminating the need for interior discretization.

4.2. Discretization

4.2.1. Sampling GWN Field

Approximating the boundary integral in Equation (9) using the input points $\{\mathbf{p}_{i}\}_{i=1}^{n}$ yields

$$f(\mathbf{n})\approx\sum\limits_{i=1}^{n}a_{i}\left(w^{-}(\mathbf{p}_{i})-w^{+}(\mathbf{p}_{i})\right)\nabla_{\mathbf{n}}w^{-}(\mathbf{p}_{i}).$$

Given that $w$ and its gradient can be analytically computed using Equation (6), the primary task is to discretize $w^{-}(\mathbf{p}_{i})$ and $w^{+}(\mathbf{p}_{i})$. In our approach, the discretization involves identifying two distinct points $\mathbf{p}^{+}_{i}$ and $\mathbf{p}^{-}_{i}$ for each input point $\mathbf{p}_{i}$. These points are strategically positioned such that $\mathbf{p}^{+}_{i}$ is located outside $\Omega$ and $\mathbf{p}^{-}_{i}$ inside $\Omega$. This arrangement ensures that the approximations $w^{+}(\mathbf{p}_{i})\approx w(\mathbf{p}_{i}^{+})$ and $w^{-}(\mathbf{p}_{i})\approx w(\mathbf{p}_{i}^{-})$ hold true.

(a) Normal displacement (b) Voronoi diagram

Figure 1. Illustration of GWN sampling strategies. Input points are represented as large grey circles, with Voronoi vertices shown as small yellow circles. The interior $\Omega^{-}$ is shaded blue and the exterior is in gray. A representative point $\mathbf{p}_{i}$ is highlighted in black, with its associated samples $\mathbf{p}_{i}^{+}$ and $\mathbf{p}_{i}^{-}$ shown as red and blue triangles, respectively. (a) The normal displacement strategy, though straightforward and easy to implement, is ineffective for models with thin structures. In this example, it positions both $\mathbf{p}_{i}^{+}~{}(=\mathbf{p}_{i}+\mathbf{n}_{i})$ and $\mathbf{p}_{i}^{-}~{}(=\mathbf{p}_{i}-\mathbf{n}_{i})$ on the same side of the target surface (specifically, both inside the surface). (b) Using Voronoi diagrams, we identify $\mathbf{p}_{i}^{\pm}$ by selecting two Voronoi vertices associated with $\mathbf{p}_{i}$ that best align with the current normal $\mathbf{n}_{i}$. This approach ensures that for the majority of the input points, their chosen $\mathbf{p}^{\pm}$ are consistently positioned on opposite sides of $\partial\Omega$, thereby providing a reliable sampling of the GWN field around the target surface.

A straightforward approach for identifying $\mathbf{p}_{i}^{\pm}$ is to set $\mathbf{p}_{i}^{\pm}=\mathbf{p}_{i}\pm\epsilon\mathbf{n}_{i}^{{}^{\prime}}$, where $\mathbf{n}_{i}^{{}^{\prime}}$ denotes the current estimate of the normal and $\epsilon~{}(>0)$ is a predetermined value. However, this strategy does not reliably ensure that $\mathbf{p}_{i}^{-}$ and $\mathbf{p}_{i}^{+}$ are located inside and outside the surface, respectively, especially in cases involving complex geometries and thin structures (see Figure 1 (a)).

To avoid these issues, we opt for Voronoi vertices to discretize $\mathbf{p}^{\pm}$. Given that many Voronoi cells are infinite, we truncate them using an enlarged bounding box of the input point cloud $\mathcal{P}$. Both the Voronoi vertices and the vertices along the clipped boundary are considered as potential candidates for $\mathbf{p}^{\pm}$. This strategy is supported by two key observations: Firstly, most Voronoi vertices, although not necessarily far away from the point cloud, tend to be on opposite sides of the target surface in practice, thereby providing suitable candidate positions for $\mathbf{p}^{\pm}$. Secondly, as demonstrated in (Xu et al., 2023), Voronoi vertices exhibit considerable robustness against noise so that we can use them to distinguish interior and exterior. Computational results further confirm that for models with low levels noise, employing Voronoi vertices yields a reliable sampling of the GWN field around the surface. For a more detailed discussion, refer to the supplementary material. In our implementation, for each $\mathbf{p}_{i}$, we identify $\mathbf{p}_{i}^{+}$ and $\mathbf{p}_{i}^{-}$ as the two Voronoi vertices of the Voronoi cell associated with $\mathbf{p}_{i}$ that best align with the current normal $\mathbf{n}_{i}$ (see Figure 1 (b)).

Remark 2

At the beginning of the optimization process with randomly initialized normals, some points $\mathbf{p}_{i}$ may experience inside-out flipping, where $\mathbf{p}_{i}^{+}$ is positioned inside and $\mathbf{p}_{i}^{-}$ outside. This inside-out configuration does not hinder our boundary energy maximization, as the samples $\mathbf{p}_{i}^{+}$ and $\mathbf{p}_{i}^{-}$ still remain on opposite sides of the target surface. Computational results show that the energy maximization process effectively reduces the frequency of such inside-out cases. For example, on the Bunny model, after 10 iterations, all $\mathbf{p}_{i}^{+}$ are consistently positioned outside and all $\mathbf{p}_{i}^{-}$ inside. See Figure 4 in the supplementary material.

4.2.2. Energy Maximization

Using the Voronoi vertices to sample the GWN field around $\partial\Omega$, we can express the maximization problem in Equation (9) as:

$$\max\sum\limits_{i=1}^{n}a_{i}\left(w(\mathbf{p}_{i}^{-})-w(\mathbf{p}_{i}^{+})\right)\nabla_{\mathbf{n}_{i}}w(\mathbf{p}_{i})$$

Figure 2. Visualization of the GWN field $w$ and its gradient $\|\nabla w\|^{2}$ for globally consistent normals (top) and random normals (bottom). The decomposed boundary energies $f(\mathbf{n}_{n})$ and $f(\mathbf{n}_{t})$ are (167.53, -3.21$\times 10^{-15}$) and ($4.77\times 10^{-6}$, $2.11\times 10^{-5}$), respectively.

Directly maximizing the above boundary energy often results in normals that, while perpendicular to the boundary surface, are not globally consistent. The inconsistency is mainly due to ambiguities associated with the jump boundary conditions. Figure 3 shows two common types of ambiguities encountered in 3D space: Case 1, where $w^{+}=-1$ and $w^{-}=0$, and Case 2, where $w^{+}=1$ and $w^{-}=2$. A direct approach to resolving the ambiguity in jump boundary conditions involves constraining the winding numbers to fall within the range of 0 and 1, that is:

$$0\leq w(\mathbf{p}_{i}^{\pm})\leq 1,\quad i=1,2,\ldots,n.$$

  Input: Point cloud $\mathcal{P}=\{\mathbf{p}_{i}\}_{i=1}^{n}$   Output: Globally consistent outward unit normal for each point   Compute the 3D Voronoi diagram with $\{\mathbf{p}_{i}\}$ as seeds   for each point $\mathbf{p}_{i}\in\mathcal{P}$ do      Assign a random normal $\mathbf{n}_{i}=(\cos u_{i}\sin v_{i},\sin u_{i}\sin v_{i},\cos v_{i})^{T}$      Calculate the geodesic Voronoi area $a_{i}$      Collect the Voronoi vertices $\{\mathbf{q}_{i}^{k}\}_{k=1}^{m_{i}}$   end for   while the L-BFGS solver does not converge do      for each point $\mathbf{p}_{i}\in\mathcal{P}$ do         Identify $\mathbf{p}_{i}^{+}=\arg\min\limits_{\{\mathbf{q}_{i}^{k}\}}\angle(\mathbf{q}_{i}^{k}-\mathbf{p}_{i},\mathbf{n}_{i}^{{}^{\prime}})$         Identify $\mathbf{p}_{i}^{-}=\arg\max\limits_{\{\mathbf{q}_{i}^{k}\}}\angle(\mathbf{q}_{i}^{k}-\mathbf{p}_{i},\mathbf{n}_{i}^{{}^{\prime}})$         Calculate $w(\mathbf{p}_{i}^{+})$ and $w(\mathbf{p}_{i}^{-})$ as in Equation (6)         Calculate $\nabla w(\mathbf{p}_{i})=\sum\limits_{j=1,j\neq i}^{n}a_{i}\nabla P(\mathbf{p}_{i},\mathbf{p}_{j})$      end for      Calculate the boundary energy as in Equation (10) and its gradient      Update $\{u_{i},v_{i}\}_{i=1}^{n}$ using the L-BFGS solver   end while

Algorithm 1 Globally consistent point orientation via optimizing the boundary energy of GWN

However, this method introduces $2n$ inequality constraints, significantly increasing the complexity of the optimization problem. Instead, we propose not to directly restrict the value of $w$ using the above inequality constraints. Instead, we 1) replace the GWN difference $w^{-}-w^{+}$ in Equation (9) with $|w^{-}|-|w^{+}|$, which encourages $w\geq 0$, thereby penalizing Case 1; and 2) add an additional term

$$g\left(w(\mathbf{p}_{i}^{\pm})\right)=\left(1-e^{\max\{w(\mathbf{p}_{i}^{+}),1\}}\right)+\left(1-e^{\max\{w(\mathbf{p}_{i}^{-}),1\}}\right),$$

which encourages $w\leq 1$, thereby penalizing Case 2. With these changes, the final energy becomes an unconstrained optimization problem:

(10) $$\max\sum\limits_{i=1}^{n}a_{i}\left(|w(\mathbf{p}_{i}^{-})|-|w(\mathbf{p}_{i}^{+})|\right)\nabla_{\mathbf{n}_{i}}w(\mathbf{p}_{i})+g\left(w(\mathbf{p}_{i}^{\pm})\right).$$

We solve this optimization using the L-BFGS solver. During each iteration, we examine the current normal $\mathbf{n}_{i}$. Through a linear scan, we determine $\mathbf{p}_{i}^{+}$ and $\mathbf{p}_{i}^{-}$, which correspond to the Voronoi vertices forming the smallest and largest angles, respectively, with $\mathbf{n}_{i}$. To force the normals to be unit length, we parameterize each normal as $\mathbf{n}_{i}=(\cos u_{i}\sin v_{i},\sin u_{i}\sin v_{i},\cos v_{i})^{T}$, where $u_{i}$ and $v_{i}$ serve as the variables in the objective function Equation (10).

Remark 3

Our Voronoi diagram strategy does not guarantee that the selected $\mathbf{p}^{\pm}$ will always be on opposite sides of the target surface, experimental results suggest that most input points do not face the issue of $\mathbf{p}^{\pm}$ being on the same side. In the rare occasions when $\mathbf{p}_{i}^{+}$ and $\mathbf{p}_{i}^{-}$ do end up on the same side, our optimization method effectively addresses these instances by penalizing them. Specifically, when $\mathbf{p}_{i}^{+}$ and $\mathbf{p}_{i}^{-}$ are on the same side, the value of $|w(\mathbf{p}_{i}^{-})|-|w(\mathbf{p}_{i}^{+})|$ reduces, leading to lower energy. Our objective is to maximize boundary energy, and with the normals updated, it is unlikely that these same-side samplings will recur in subsequent iterations.

Figure 3. The ambiguities of jump boundary conditions in Case 1 and Case 2. For each case, we plot the boundary energy (BE) and the Dirichlet energy (DE) with and without the constraints $0\leq w\leq 1$. By restricting the range of the winding numbers, DE and BE converge when the associated normals are globally consistent. We also visualize the GWN field using cut views and indicate the normals for the cross sections.

Remark 4

In our current implementation, we calculate $\nabla w$ and $w$ in a brute-force manner. Let $n$ be the number of input points in $\mathcal{P}$ and $m$ represent the number of Voronoi vertices. Each iteration begins with an $O(n+m)$ traversal of Voronoi vertices to identify all $\mathbf{p}^{\pm}$s. Subsequently, for each input point, computing $\nabla w$ and $w$ via the Poisson kernel requires $O(n)$ time. Therefore, the computational cost per iteration is $O(n^{2}+n+m)$. To enhance runtime efficiency, one possible strategy is to accelerate the computation using the fast multipole method (Greengard and Rokhlin, 1987). This enables querying neighboring points in $O(\log n)$ time, potentially reducing the overall time complexity per iteration to $O(n\log n)$.

Figure 4. Illustration of our algorithm on a 2D toy model. In different iterations, we visualize the normals, $\mathbf{p}^{+}$, $\mathbf{p}^{-}$ (row 1), winding number field (row 2), $w(\mathbf{p}^{+})$ and $w(\mathbf{p}^{-})$ (row 3). We highlight the regions with significant variations in normals between adjacent iterations using gray lines and red points.

4.3. Contrasting Our Method with Previous Approaches

Our method, while sharing similarities with the works of Takayama et al. (2014), PGR (Lin et al., 2022), and GCNO (Xu et al., 2023), particularly in their use of the generalized winding number field, diverges significantly in several key aspects.

Ours vs. Takayama et al. (2014).’s method

Although both methods aim at optimizing objective functions that are related to to the Dirichlet energy, there are two key differences that set our approach apart: 1) Takayama et al. primarily addressed scenarios where point normals are perpendicular to the boundary surface yet inconsistent. In such cases, we observe a significant increase in boundary energy (Figure 2). However, when normals are initialized randomly, the boundary energy tends to diminish to nearly zero, indicating a fundamental difference in the objective functions. 2) Their approach relies on a brute-force solver, exhaustively searching for the optimal normal configuration. As a result, its applicability is confined to small and simple models. In contrast, our method utilizes an L-BFGS solver.

Ours vs. GCNO (Xu et al., 2023)

Our approach fundamentally differs from GCNO in terms of methodology, discretization techniques and time complexity. 1) Methodology: Unlike GCNO, which regularizes the GWN of Voronoi vertices to exact values of 0 and 1, our method focuses on the boundary energy by evaluating the GWN difference $|w(\mathbf{p}_{i}^{-})|-|w(\mathbf{p}_{i}^{+})|$ and the GWN gradient. This strategy does not depend on precise GWN values, offering a more flexible approach. 2) Discretization: GCNO assumes a uniform distribution of Voronoi vertices both inside and outside the object—an assumption that is less reliable for models with complex geometries. Figure 9 shows an example where GCNO is unable to optimize its objective function when the assumption does not hold. Our method, by contrast, does not rely on such an assumption, making it adaptable to a broader range of 3D models. 3) Time complexity: GCNO requires computing $w$ for all Voronoi vertices associated with each point in the input point cloud. Our method, however, requires computing $\nabla w$ for $\mathbf{p}$ and $w$ for two specific points, $\mathbf{p}^{\pm}$, within the Voronoi cell. This leads to significantly lower computational demands in comparison to GCNO. Specifically, GCNO has a time complexity $O(nm)$ for each iteration, while ours is $O(n^{2}+n+m)$ per iteration. Given that the number of input points $n$ is typically much lower than the number of Voronoi vertices $m$ ($n\ll m$), as exemplified by the Bunny model ($n=8,830$, $m=224,085$), our method outperforms GCNO in terms of scalability and runtime performance. Experimental results indicate that GCNO is suited for models with up to 10K points, typically generating within an hour. However, for larger models containing 40K points, GCNO fails to complete within 24 hours. In contrast, our method can deliver results in just three hours. See the supplementary material for the details of the scalability test.

Ours vs. PGR (Lin et al., 2022)

PGR focuses on regularizing the GWN directly at the input points. In the context of noisy point clouds, achieving a GWN value of exactly 0.5 for these points is challenging. Unlike PGR, our method does not constrain GWN values. Instead we regularize the differences and gradients of the GWN both inside and outside of the input model, significantly enhancing our method’s robustness to noise. Furthermore, PGR requires solving a dense linear system, restricting its applicability to small models. In contrast, our method employs an L-BFGS solver to optimize the boundary energy, offering a solution that is both fast and memory-efficient.

Models (# points, Fig.)	$\text{Dipole}^{*}$	$\text{PGR}^{*}$	iPSR	$\text{NeuralGF}^{*}$	GCNO	Ours
BS (4000, supp.)	4.5/5.7/0.17/28/2.0	17/21/0.12/2.0/0.50	3.7/4.5/0.11/9.7/0.094	5.1/6.0/0.15/640/1.9	8.1/9.9/0.21/810/0.067	3.6/4.3/0.11/95/0.075
Bird (8187, Fig. 10)	12/37/1.7/50/1.9	7.4/11/0.063/8.1/0.80	2.2/5.1/0.065/14/0.10	5.2/8.8/0.069/640/1.9	89/98/28/2800/0.14	2.0/3.4/0.060/360/0.084
linkCupTop (8602, Fig. 10)	89/120/3.6/60/2.2	23/28/0.27/6.7/0.85	24/28/0.20/25/0.11	87/110/3.4/640/1.9	89/110/3.3/2200/0.14	12/15/0.19/400/0.091
BunnyPeel (8648, Fig. 6)	51/82/25/53/2.0	22/26/0.099/10/0.85	12/17/0.099/18/0.10	48/70/11/640/1.9	45/63/36/2400/0.14	7.2/12/0.094/400/0.092
Horse (8657, supp.)	8.2/23/0.11/51/2.4	13/17/0.065/7.1/0.85	3.1/5.0/0.070/14/0.10	6.6/9.5/0.065/640/1.9	8.2/9.9/0.072/2500/0.14	2.9/4.2/0.065/400/0.085
397horse (8784, supp.)	8.9/28/0.19/50/2.4	15/21/0.066/7.0/0.86	3.1/5.0/0.063/15/0.10	8.1/17/0.45/640/1.9	8.2/9.9/0.064/2700/0.14	2.9/4.2/0.062/400/0.086
Bunny (8830, supp.)	6.7/9.2/0.13/62/2.1	15/20/0.14/7.3/0.86	4.9/6.9/0.11/13/0.11	5.4/6.6/0.17/640/1.9	8.9/11/0.16/2700/0.14	3.1/4.1/0.12/410/0.094
Cup (10130, supp.)	4.4/19/0.14/60/2.0	11/14/0.13/9.7/1.0	1.5/2.5/0.15/11/0.11	3.9/4.5/0.13/640/1.9	7.4/9.5/0.13/3700/0.16	1.7/2.2/0.13/550/0.088
Botijo (11305, supp.)	14/37/0.23/63/2.4	16/20/0.080/13/1.1	3.0/4.9/0.077/15/0.12	15/38/2.3/640/1.9	35/64/54/4100/0.19	2.6/3.7/0.077/660/0.086
30cup (11311, supp.)	4.4/15/0.090/65/2.4	14/16/0.091/14/1.1	1.7/3.1/0.093/13/0.12	32/64/97/640/1.9	9.2/11/0.097/4400/0.17	2.1/3.0/0.086/670/0.095
Candle (12221, Fig. 6)	65/93/1.12/67/2.4	42/53/0.97/18/1.2	18/25/0.099/34/0.14	30/48/11/640/1.9	50/72/36/7900/0.20	14/26/0.094/790/0.089
Elk (12541, Fig. 6)	17/48/7.8/77/2.1	19/30/0.17/22/1.3	3.0/6.8/0.12/18/0.12	17/43/5.9/640/1.9	48/83/20/7400/0.20	2.9/4.6/0.11/1000/0.10
Lion (13138, Fig. 6)	14/28/0.57/71/3.2	21/28/0.061/16/1.4	9.1/14/0.068/14/0.11	12/17/0.068/640/1.9	19/31/0.20/5100/0.21	6.5/11/0.057/900/0.086
Vase (13375, Fig. 10)	50/84/8.5/81/2.5	26/33/0.26/18/1.4	12/16/0.23/39/0.13	77/104/5.3/640/1.9	88/114/11/6200/0.21	5.8/8.1/0.23/920/0.099
Fertility (13802, supp.)	52/94/120/76/3.2	15/20/0.073/20/1.5	2.3/4.0/0.10/19/0.14	5.0/6.3/0.081/640/1.9	8.0/11/0.091/10000/0.22	2.2/3.2/0.075/970/0.088
Fandisk (13854, supp.)	4.7/9.9/0.13/79/2.2	18/25/0.12/22/1.5	2.6/5.3/0.12/17/0.16	4.6/6.0/0.18/640/1.9	9.5/12/0.11/5500/0.20	2.8/4.3/0.11/990/0.092
Felinei (15057, supp.)	36/72/25/81/3.0	19/24/0.092/23/1.6	5.7/10/0.088/18/0.11	30/57/16/640/1.9	89/97/17/8500/0.24	4.1/7.2/0.084/1200/0.092
Pulley (16116, supp.)	13/31/0.12/86/3.1	26/33/0.092/23/1.8	5.0/7.5/0.14/16/0.16	12/17/0.090/640/1.9	7.4/8.8/0.12/12000/0.22	3.8/5.3/0.083/1300/0.092
Average	25/46/11/64/2.4	19/25/0.16/14/1.1	6.5/9.5/0.11/18/0.12	22/35/8.0/640/1.9	35/46/9.5/5100/0.17	4.6/7.0/0.11/690/0.090

Table 1. Quantitative results on clean point clouds. We present five metrics: mean and standard deviation of angular discrepancies (degree), CD ($10^{-3}$), runtime (seconds), and GPU or CPU memory utilization (MB). * indicates GPU-based. The best results are emphasized in bold. See results on noisy point clouds in the supplementary material.

Figure 5. Reconstruction results for the 18 test models.

5. Experimental Results

Experimental setup

We perform our experiments on a PC equipped with an i7-12700H CPU and 16GB of RAM. The GPU-based approaches (Lin et al., 2022; Metzer et al., 2021; Li et al., 2023b) were executed on a single RTX 3090. We evaluate our method using 18 models, with diverse geometry and topology (see Figure 5). All point clouds are uniformly scaled to fit within the cube $[-1,1]^{3}$. Noise was synthetically generated using a standard Gaussian distribution with $\mu=0.0$ and $\sigma^{2}=1.0$, introducing random displacements to each point. We modulate the noise intensity to 0.5% and 0.75% of the diagonal length of the bounding box. Our method involves a single parameter: the neighborhood size, which is used for approximating the Voronoi area, $a_{i}$, in the discretization of the GWN. Despite the diversity in model geometries, we use a neighborhood size of 15 across all tests. We use the L-BFGS solver for optimizing the unconstrained boundary energy. After establishing the point orientations, we use screened Poisson surface reconstruction (Kazhdan and Hoppe, 2013) on the oriented point clouds, generating watertight surfaces. We compare our method with five state-of-the-art approaches GCNO (Xu et al., 2023), Dipole (Metzer et al., 2021), PGR (Lin et al., 2022), NeuralGF (Li et al., 2023b), and iPSR (Hou et al., 2022), utilizing the official implementations provided by the authors along with their default parameter settings. We assess the quality of orientation and reconstruction by analyzing the angular discrepancy between the predicted normals and the ground-truth normals, along with the Chamfer Distance (CD), which measures the closeness between the ground-truth surface and the reconstructed surface. To quantitatively describe this distribution, we report the mean $\mu$ and the standard deviation $\sigma$ of the angle differences in degree. The smaller the metrics, the higher the consistency of the predicted normals. We also document the running time and the peak memory consumption for computational efficiency.

Orientation & reconstruction results

We present the performance of various methods across 18 noise-free models in Table 1, and provide results under noisy conditions in the supplementary material. In the noise-free scenarios, our method achieves the best scores in $\mu$, $\sigma$ and CD for 83%, 88%, and 88% of the models, respectively. Our method demonstrates resilience to low level noise, as it is a global method that can capture the normal consistency from a global perspective. See Figures 6 and 7 for the comparison of the reconstructions. Our method is particularly effective in processing models characterized by complex topology and thin structures, areas where competing methods frequently encounter difficulties.

Runtime

As analyzed in Section 4.2, our method’s time complexity is $O(n^{2}+n+m)$ per iteration. Typically, $m$ increases linearly with $n$ for common models. However, it is important to note that in the worst-case scenario, $m$ can scale quadratically with $n$. Consequently, our method is faster than GCNO, which operates with a time complexity of $O(nm)$ per iteration. Computational results demonstrate that our method is 5-10 times faster than GCNO across 18 test models, with an average speed-up factor of 7.3.

Space complexity

Our algorithm has space complexity of $O(n+m)$, since it only needs to store the input point clouds and its associated Voronoi diagram. Both PGR (Lin et al., 2022) and VIPSS (Huang et al., 2019) necessitate solving dense linear systems, resulting in a space complexity of $O(n^{2})$. This quadratic space complexity constrains them to smaller-scale models, generally limited to those with up to 10K points. In contrast, our method scales up to 100K points.

Comparison with supervised methods

We also compare our method with recent supervised approaches, NGLO (Li et al., 2023b) and SHS-Net (Li et al., 2023c). NGLO operates in two stages: The initial phase, which is unsupervised, produces a set of normals. While globally consistent, these may not always be perpendicular to the surface. This is followed by a supervised module aimed at refining the normals’ accuracy. However, since the second stage is mainly aimed at fine-tuning, NGLO may face issues with normal flipping if the initial stage fails to achieve globally consistent normals. SHS-Net is an end-to-end network that extracts both local, dense features and global, sparse features from the input point cloud. It orients points by penalizing discrepancies between these features in the latent space. Utilizing both local and global features is a double-edged sword. On one hand, it renders SHS-Net highly effective for simple models, where these features can be well aligned. On the other hand, it struggles with input models characterized by complex geometry and topology, primarily due to the significant discrepancies between local and global features. In Figure 10, we compare our method to NGLO and SHS-Net for four models featuring thin structures and complex topology. To assess the robustness of these methods, we intentionally introduce noise or outliers to each test model. We observed that both NGLO and SHS-Net reconstruct the chick body well, which has a relatively simple geometry. However, they encounter difficulties with the thin structures of the chick’s feet and fail to capture the complex topology of the other three models.

Limitations

Despite the demonstrated robustness and accuracy of our method, it is important to acknowledge its current limitations. As highlighted in (Takayama et al., 2014), the Dirichlet energy of GWN only measures the magnitude of the oscillation of a function and does not explicitly indicate the inside-outside orientation on non-manifold surfaces. As the proposed boundary energy is derived from the Dirichlet energy, our method is suitable only for point clouds sampled from manifold surfaces. Another limitation is the scalability issue due to the quadratic time complexity of our method. While iPSR (Hou et al., 2022) can handle large-scale models with millions of points, our method is constrained to smaller-scale models, with around 100K points.

6. Conclusions

We introduce a new method for globally consistent orientation for point clouds. Our method initializes each point with a randomly assigned normal, and then iteratively updates its normal by maximizing the boundary energy of the generalized winding number field associated with the normals. When the algorithm terminates, we restore the harmonicity of the GWN field, whose gradients are the globally consistent orientations. Computational results show that our method is particularly effective in handling models with complex topology and thin structures, areas that existing methods have trouble with.

Figure 6. Comparison with existing methods. Rows 1 & 2: Sparse point clouds (Kitten: 500 points; Hilbert Cube: 2,500 points). Row 3 & 4: The Elk model featuring thin structures with 0.5% and 0.75% noise level. Rows 5-8: Models with 0.75% noisy level. Rows 9 & 10: Two challenging models featuring thin plates, slender tubes, and closely adjacent faces.

Figure 7. The reconstruction results of several state-of-the-art point orientation methods (Dipole (Metzer et al., 2021), NeuralGF (Li et al., 2023b), GCNO (Xu et al., 2023), iPSR (Hou et al., 2022)) and our method. Our method is more effective than the other methods at handling point clouds with complex geometries and topologies. The quantitative numbers below each model indicate the mean and the standard deviation of angle differences (degree), and Chamfer distances ($10^{-3}$), respectively.

Figure 8. Comparison with VIPSS on the horse model with 4k points under 3 different noise levels: 0.0%, 0.5% and 0.75%. The last scenario is the results with 1% outliers but no noise. For each scenario, we show the input point cloud, the results from VIPSS, and our results. While VIPSS handles clean data and low-level noise effectively, it struggles with high-level noise and outliers.

Figure 9. The efficacy of GCNO heavily relies on a uniform distribution of Voronoi vertices both inside and outside the point cloud. This condition is often unattainable for models characterized by complex geometry and topology, leading to inaccuracies in point orientation and, consequently, subpar reconstruction outcomes. In the left, we plot the the objective functions and gradient norms for GCNO (red) and our method (blue) for an example where GCNO fails due to stopping optimizing its objective at an early stage. In the right, we show our results (top) and GCNO’s (bottom).

Figure 10. Comparison with the SOTA deep-learning approaches on point clouds subjected to 0.75% noise levels (top) and 1% of outliers (bottom). The 3-tuple in each figure is the error metrics, $\mu$ (degree), $\sigma$ (degree) and CD ($10^{-3}$).