Neural Octahedral Field: Octahedral prior for simultaneous smoothing and sharp edge regularization

Reconstructing surface has been extensively researched in recent decades and continues to evolve rapidly with the advancement of deep learning. We focus on reviewing the approaches most related to ours and refer the interested readers to more detailed reviews (Huang et al., 2022; Sulzer et al., 2024).

Point Set Surface (Alexa et al., 2001) (PSS) fits polynomials to approximate the surface locally. The query points are evaluated using the Moving Least Square (MLS) (Levin, 2004), with the gradient of the polynomials indicating the surface normals (Alexa and Adamson, 2004). Algebraic Point Set Surface (APSS) (Guennebaud and Gross, 2007) fits the algebraic sphere to improve the stability of PSS over regions of high curvature. Robust Implicit MLS (Öztireli et al., 2009) integrates robust statistics with local kernel regression to handle outliers. Fleishman et al. (2005) apply forward search (Atkinson and Riani, 2012) to filter outliers that cannot be identified by robust statistics. Variational Implicit Point Set Surface (VPSS) (Huang et al., 2019) leverages the eikonal constraint to ease the need for input normals, but at the cost of cubic complexity. In contrast to local fitting, Poisson Surface Reconstruction (PSR) (Kazhdan et al., 2006) represents the surface globally as the level set of implicit indicator function, whose gradient matches the input normal. Screened Poisson Surface Reconstruction (SPSR) (Kazhdan and Hoppe, 2013) integrates positional constraints to balance precision and smoothness. Recently, Neural Kernel Surface Reconstruction (NKSR) (Huang et al., 2023) integrates aspects from both ends–it models surface as the zero level set of the global implicit function, but with its value determined by learned proximity (neural kernel) with respect to input samples.

Edge-Aware Point Set Resampling (EAR) (Huang et al., 2013) applies bilateral filtering to reliable fit normals away from edges, then progressively sample towards the edge to capture sharp features. Wei et al. (2023b) predict cross field and leverage its crease align nature Jakob et al. (2015); Huang and Ju (2016) to tessellate space for convolution. Sarkar et al. (2018) randomly samples square patches and applies RANSAC (Fischler and Bolles, 1981) stores local height in square matrices and performs denoising by to low rank matrix factorization. Zeng et al. (2019) sample overlapping patches and connect samples by projection distance to perform graph smoothing. (Wei et al., 2023a) use line processes to perform smoothing up to a optimizable threshold. Deep Geometric Prior (DGP) (Williams et al., 2019) parameterize local patch with MLP and rely on its smooth prior to denoise.

Neural implicit representation encodes the surface as the zero-level set of the distance function using coordinate MLP. We primarily review its application in fitting-based surface reconstruction. Readers interested in neural field in general should check out the comprehensive survey by Xie et al. (2021).

DeepSDF (Park et al., 2019) and Occupancy Networks (Mescheder et al., 2019) parameterize the signed distance function (SDF) or occupancy classification over spatial locations, which are capable of encoding geometry of arbitrary resolutions. SAL (Atzmon and Lipman, 2020) sample in the vicinity of the input point cloud and learns an unsigned distance function. SALD (Atzmon and Lipman, 2021) further improves the quality of reconstruction with normal supervision. IGR (Gropp et al., 2020) introduces the eikonal regularization term, which under stochastic optimization, the MLP converges to a faithful SDF. When normal is not available, the eikonal regularization is insufficient to guarantee a unique solution, resulting in ambiguities and artifacts. Park et al. (2023) model the SDF gradient as the unique solution to the $p$-Poisson equation and additionally minimize the surface area to improve hole filling. DiGS (Ben-Shabat et al., 2023) initializes the MLP with a geometric sphere and minimizes the divergence of its gradient to preserve the orientation consistency. Neural Singular Hessian (NSH) (Wang et al., 2023) encourages the Hessian of MLP to have zero determinant and produces a topologically faithful initialization. However, both DiGS and NSH require annealing the regularization weights to fit finer details, which are prone to noise corruption. Our idea is to ”snapshot” the initialization as a prior, in making up for the emerging of noise (Figure 2).

In this section, we give a brief summary of the Spherical Harmonic (SH) representation of octahedral frame by Huang et al. (2011), Ray et al. (2016) and Palmer et al. (2020). A full review is available in the section A of the supplementary.

For each sample on the surface, we want to minimize the deviation between the predicted octahedral frame and its closest normal-aligned projection. Although the surface normal is not provided, the $f_{\omega}$ is always differentiable, that we utilize the gradient of its zero level set as the boundary condition to define our octahedral field:

where $\mathbf{p}_{i}^{\prime}=\mathbf{p}_{i}-f_{\omega}(\mathbf{p}_{i})\cdot\frac{\nabla f_{\omega}(\mathbf{p}_{i})}{\|\nabla f_{\omega}(\mathbf{p}_{i})\|_{2}}$ is the projection of $\mathbf{p}_{i}$ on zero level set, $\Pi(\mathbf{q},\mathbf{n})$ is the rewritten of (4) by treating both SH coefficient vector and normal as arguments, $\phi$ is a similarity measurement function.

In practice, to avoid the costly projection step, we adapt the weighting scheme from Ma et al. (2023), which we use $\nabla f_{\omega}(\mathbf{p}_{i})$ as the surrogate and weight it by distance to the zero level set. The alignment loss can then be simplified as:

where $\beta(\mathbf{p}_{i})=\exp(-100\cdot|f_{\omega}(\mathbf{p}_{i})|)$, and we choose cosine similarity $\phi(\cdot,\cdot)=(1-\langle\cdot,\cdot\rangle)$. We observe a minor difference when compared to direct projection. It is worth pointing out that the weighting is important–frames distant away surface may deviate significantly due to the existence of singularity. Note that (5) is a function of both MLPs, so we stop the gradient propagation for $f_{\omega}$ and focus on aligning $u_{\theta}$.

The design of (5) yields two benefits: First, as discussed in Section 2.4.2, projection (4) by Palmer et al. (2020) alleviates the need for additional parameterization or quadratic equality, especially when such equality is defined element-wise. Additionally, the choice of cosine similarity relaxes the unit norm constraint, which is topologically unsatisfiable for the smooth octahedral field due to the existence of singularities (Solomon et al., 2017). This also makes the output of $u_{\theta}$ a Euclidean vector space in $\mathbb{R}^{9}$, which we will elaborate on its importance in the next section.

However, unlike the boundary condition in the discrete setting, (5) is inherently a soft constraint. Compared to the exact one (Sukumar and Srivastava, 2022), it is notoriously difficult to train and requires dense boundary samples to fulfill (Raissi et al., 2019; Barschkis, 2023). It is especially challenging for our case, as small hole-like structures are a common occurrence in real-world scans. Those holes resemble cylinder volumes that require multiple singularity curves to satisfy. Compounded by the difficulty in capturing the interior of those holes using range sensors, our octahedral field cannot fully align with those regions. We will discuss its implications in Section 3.3.

The next step is to encourage the octahedral field to be smooth. Since $u_{\theta}$ directly outputs the SH coefficient vector, its MLP output space is of Euclidean topology and is also where the smoothness of the octahedral frames can be measured. This allows the use of efficient Lipschitz regularization (Liu et al., 2022) as a smoothness proxy.

Under the 1-Lipschitz activation function (ReLU, ELU, $\sin$, tanh, etc.), the norm of output variation with respect to the input is bounded by that layer’s Lipschitz constants $c_{i}$

where $u_{\theta i}$ denotes the $i$th layer of $u_{\theta}$, $\mathbf{x}$ is its input, $\mathbf{W}_{i}$ is its weight matrix, $l$ is number of layers. The LipMLP (Liu et al., 2022) initializes the Lipschitz constant per layer as the $L_{\infty}$ norm of its randomly initialized weight matrix $c_{i}=\|\mathbf{W}_{i}\|_{\infty}$. During training, $c_{i}$ is used to rescale the maximum absolute row-wise sum of the weight matrix $\mathbf{W}_{i}$ so $c_{i}=\|\mathbf{W}_{i}\|_{\infty}$ still holds. By shrinking all Lipschitz constants, we minimize the MLP’s output variation with respect to input, hence encouraging it to be smooth:

where softplus is a reparameterization to force the bound to be positive ($\text{softplus}(c_{i})\approx c_{i},\ c_{i}\gg 1$), and the $p=\infty$ is chosen purely for computational efficiency. This loss is highly efficient because it only needs first-order derivatives, and all weight matrices can be updated in parallel.

In our case, minimizing the Lipschitz bound also has geometric meaning. Denote the output for any two input $\mathbf{x}_{0}$, $\mathbf{x}_{1}$ as $\mathbf{q}_{0}=u_{\theta}(\mathbf{p}_{0})$, $\mathbf{q}_{1}=u_{\theta}(\mathbf{p}_{1})$ respectively, we have:

The Lipschitz constants bound norm of variation of SH coefficient vector (1) with respect to input spatial variations, hence its minimization is analogously to minimizing the discrete Dirichlet energy (2). Note that it is only possible because we output SH coefficients directly–the output space for most rotation representations are not Euclidean(Zhou et al., 2019), and for those are, they are not where the smoothness for octahedral frames is measured. Thus, they require more expensive gradient norm minimization or finite-difference equivalents (Huang et al., 2021; Zhang et al., 2023).

The core of our method is to exploit the interpolation property of the octahedral frames (Section 2.4.3) to simultaneously smooth and accentuate the sharp edges of the companion implicit geometry (Figure 4). We do it by matching the gradient of SDF to one of the six representation vectors of the octahedral frame.

Recall that align loss (5) is a function of both MLPs–we can simply flip the gradient propagation by fixing $u_{\theta}$ while updating $f_{\omega}$, which resembles alternating optimization. In examining the loss manifold (Figure 6), we find $L_{1}$ or $L_{2}$ loss works better than cosine similarity, so we adjust it slightly:

the normalization is employed because the projection (4) always yields valid SH 4 coefficient vector hence of unit norm. Given that we work on the same samples in the alignment phrase, the normalized $u_{\theta}(\mathbf{p}_{i})$ is unlikely to deviate from the octahedral variety. We remove the distance weighting, as we empirically find it contributes little here, possibly due to the inherent ambiguity in regularization process–the cubic symmetry of octahedral frame means the matched SDF gradient is not unique. This is also why our method does not contribute to improving the normal consistency of unoriented point cloud reconstruction. Conversely, we require a good initialization, such as DiGS or NSH, to realize our method’s full potential.

Note that smoothing is at the cost of reducing the capacity of the model to represent high-frequency information (figure 5). Therefore, we use NSH as our backbone and follow their practice, we also sample close-surface samples and off-surface samples, denoted $\mathbf{p}_{\text{close}}$ and $\mathbf{p}_{\text{off}}$ respectively. The off-surface point cloud is sampled uniformly in the unit bounding box, while close samples are drawn from the normal distribution, with sigma being the maximum KNN distance of $k=51$, centered at the scanning input $\mathbf{p}$. Our final losses are:

where $\lambda$s are balancing weight. Note that evaluating our loss only requires input point cloud.

We fix $\lambda_{\text{NSH}}=3$, $\lambda_{\text{eikonal}}=50$, $\lambda_{\text{off}}=100$, $\lambda_{\text{align}}=100$, $\lambda_{\text{regularize}}=10$, $\lambda_{\text{lip}}=10^{-6}$, and adapt two scheduling schemes:

• •

  Low-noise: Set $\lambda_{\text{positional}}=7000$, annealing $\lambda_{\text{NSH}}$ to $3\times 10^{-4}$ after $10\%$ of training step, start alignment and smoothing at $40\%$, regularization at $60\%$.

• •

  High-noise: Set $\lambda_{\text{positional}}=3500$, annealing $\lambda_{\text{NSH}}$ to $3\times 10^{-3}$ after $10\%$ of training step, start alignment and smoothing at $20\%$, regularization at $40\%$.

This might seem counterintuitive because we change the schedule of losses rather than adjusting their weights. The rationale is that we want our octahedral field to capture the geometry before it is fully corrupted, while not overfitting the overly smooth initialization. Given that small holes are often not observed sufficiently during scanning and emerge later in the fitting process (Figure 5), we postpond scheduling to capture those small details. However, when the noise level is high, the noise is corrupted so quickly that we cannot afford to leave more to come, so we start our alignment shortly after its annealing stage.

Note that $\lambda_{\text{positional}}$ and $\lambda_{\text{NSH}}$ in the low noise scheme are the default values in their original implementation, the tuning in our high noise scheme aims to slow down the noise corruption. The resulting effect is illustrated as the fine-tuned one in Figure 2.

Following the practice of LipMLP, we spatially scale the input by $100$–it is equivalent to premultiply first weight matrix by the same amount Sitzmann et al. (2020) to speed up convergence, so our octahedral field can ”snapshot” as quickly as possible–our method cannot fix a already corrupted geometry.

For quantitative evaluation, we follow DiGS (Ben-Shabat et al., 2023) and report Chamfer distance $(\times 10^{3})$ (Fan et al., 2016) and Hausdorff distance $(\times 10^{2})$ over 1M points randomly sampled on surface or point cloud, depending on the method’s output. Given that neural implicit representations are susceptible to floating artifacts (Figure 8), we additionally report the F-score (Knapitsch et al., 2017). It clamps the closest matching distance by a threshold and reports as a percentage and hence is less sensitive to the mismatching of surplus parts. Following Tatarchenko et al. (2019), we use the distance threshold of $0.5\%$. Note that we do not postprocess the output for any methods and leave the extracted meshes or point clouds as is.

We first compare our method with existing approaches based on neural implicit representation fitting, namely DGR (Williams et al., 2019), SIREN (Sitzmann et al., 2020), DiGS (Ben-Shabat et al., 2023), NSH (Wang et al., 2023). We use the Surface Reconstruction Benchmark (SRB) (Berger et al., 2013) dataset, which consists of 5 scans that exhibit triangulation-based scanning patterns. We fit the full point cloud for all methods and use our low-noise scheme for our fitting process. DGR (Williams et al., 2019) outputs a dense denoised point cloud, which we use SPSR for visualization purposes.

Our method achieves the best F-score, second best Chamfer distance with the smallest standard deviation (Table 1), which is aligned with our cleaner reconstruction quality (Figure 7).

We further evaluate our methods on two widely evaluated datasets, ABC (Koch et al., 2019) and Thingi10k (Zhou and Jacobson, 2016). The former consists of CAD models of sharp edges, whereas the latter is made up of more general shapes. We leverage the 100 test split from Points2Surf (Erler et al., 2020) for both datasets and further use Blensor (Gschwandtner et al., 2011) to simulate the time-of-flight (ToF) scanning process. Specifically, we set sensor resolution to $176\times 144$, focal length to 10mm and scanning each object spherically for 30 scans, resulting in point clouds of size ranging from 20k to 100k. We generate scanning data for two noise levels, $\mathcal{N}(0,0.01L)$ and $\mathcal{N}(0,0.005L)$ (so $400$ in total), where $L$ is the length of the maximum edge of the model’s bounding box. Our two schemes in Section 3.3 are tuned with respect to these two noise levels, respectively, although we have shown that the low-level scheme also works for SRB as well. For reference, data-driven methods (Erler et al., 2020; Huang et al., 2023) typically randomize noise with the sigma in between $0.01L$ and $0.05L$, so our noise level is their lower bound.

We compare with 9 methods, 2 axiomatic ones (APSS (Guennebaud and Gross, 2007), SPSR (Kazhdan and Hoppe, 2013)), 1 resampling (EAR (Huang et al., 2013)), 2 patch-based denoising (GLR (Zeng et al., 2019), LP (Wei et al., 2023a)), 3 implicit fitting (SIREN (Sitzmann et al., 2020), DiGS (Ben-Shabat et al., 2023), NSH (Wang et al., 2023)) and 1 with data prior (NKSR (Huang et al., 2023)). For point cloud visualization, we use SPSR if the method outputs normals, otherwise we use the Advancing Front (Zienkiewicz et al., 2013). The full results are shown in figure 10, 11, and the quantitative evaluation is provided in the supplementary. Although our method shrines the most with CAD models, it succeeds in denoisng more general shapes such as sculptures as well. However, with noise sigma $0.01L$, our method somethingswould output a low-poly-like reconstruction–the likely cause would be, when the noise level is large, the implicit geometry evolves quickly such that constant changing surface normal makes the octahedral field hard to converge.

In our experiment, we find NSH, so is our method, fails to reconstruct three cases (Figure 8) that gives erroneously large closest matching distances ($50\times$ larger than the second worst case). Therefore, we report both original metrics and those with failure cases removed.