HLO: Half-kernel Laplacian Operator for Surface Smoothing

TL;DR

This paper introduces the Half-kernel Laplacian Operator (HLO), a novel surface smoothing method that better preserves features and reduces artifacts compared to traditional Laplacian-based techniques, especially on noisy meshes.

Contribution

The paper proposes the HLO, a new Laplacian operator that improves feature preservation and artifact reduction in surface smoothing, with a simple iterative approach and only one parameter.

Findings

HLO outperforms traditional Laplacian smoothing in feature preservation.

HLO significantly reduces shrinkage artifacts.

HLO is effective on high-noise meshes.

Abstract

This paper presents a simple yet effective method for feature-preserving surface smoothing. Through analyzing the differential property of surfaces, we show that the conventional discrete Laplacian operator with uniform weights is not applicable to feature points at which the surface is non-differentiable and the second order derivatives do not exist. To overcome this difficulty, we propose a Half-kernel Laplacian Operator (HLO) as an alternative to the conventional Laplacian. Given a vertex v, HLO first finds all pairs of its neighboring vertices and divides each pair into two subsets (called half windows); then computes the uniform Laplacians of all such subsets and subsequently projects the computed Laplacians to the full-window uniform Laplacian to alleviate flipping and degeneration. The half window with least regularization energy is then chosen for v. We develop an iterative approach to apply HLO for surface denoising. Our method is conceptually simple and easy to use because it has a single parameter, i.e., the number of iterations for updating vertices. We show that our method can preserve features better than the popular uniform Laplacian-based denoising and it significantly alleviates the shrinkage artifact. Extensive experimental results demonstrate that HLO is better than or comparable to state-of-the-art techniques both qualitatively and quantitatively and that it is particularly good at handling meshes with high noise. We will make our source code publicly available.

Tables

Table 1. Table 1: Parameters of the state-of-the-art mesh smoothing methods.

Method	Number of Parameters	Parameters
BMF	1	$k_{iter}$: number of iterations.
UNF	3	$T$: threshold for controlling the averaging weights. $n_{iter}$: number of iterations for normal update. $v_{iter}$: number of iterations for vertex update.
BNF(Local)	3	$v_{iter}$: number of iterations for vertex update. ${\sigma }_s$: variance parameter for the spatial kernel. $n_{iter}$: number of iterations for normal update.
L0	6	$\lambda$: weight for the L0 term in the target function.. ${\alpha }_0,{\beta }_0$: initial values for $\alpha$ and $\beta$ ${\mu }_{\alpha },\mu$: update ratios for $\alpha$ and $\beta$. ${\beta }_{max}$: maximum value of $\beta$.
GNF	3	$v_{iter}$: number of iterations for vertex update. ${\sigma }_r$: variance parameter for the range kernel. $n_{iter}$: number of iterations for normal update.
NLLR	3	${\sigma }_M$: the noise variance. $N_k$: number of similar vertices. $v_{iter}$: number of iterations for vertex update.
Ours	1	$k_{iter}$: number of iterations.

Table 2. Table 2: Quantitative comparisons with the representative mesh smoothing methods.

Models	Methods	MSAE $(\times {10}^{-2})$	$E_v$ $(\times {10}^{-3})$	Parameters
Amadillo (${\sigma }_n=0.5l_e$) (Figure 8) $\left|V\right|$: 43243 $\left|F\right|$: 86482	BMF BNF CNR GNF UNF L0 NLLR OURS	47.54 53.58 34.28 44.90 49.20 9.59 32.18 9.33	3.34 3.11 3.39 3.10 3.07 3.52 3.03 3.57	($30$) $(10,0.35,10)$ $(synthetic,10)$ $(10,0.35,10)$ $(0.5,10,10)$ $\left(Default\right)$ $\left(Default\right)$ ($5$)
iHbunny (${\sigma }_n=0.5l_e$) (Figure 9 row 1) $\left|V\right|$: 34834 $\left|F\right|$: 69451	BMF BNF CNR GNF UNF L0 NLLR OURS	5.07 3.57 29.7 3.39 3.87 4.32 25.24 6.25	4.99 4.39 3.96 4.25 4.38 4.70 3.58 4.32	($30$) $(10,0.35,20)$ $(synthetic,10)$ $(10,0.35,20)$ $(0.4,20,10)$ $\left(Default\right)$ $\left(Default\right)$ ($5$)
Nicolo (${\sigma }_n=0.2l_e$) (Figure 9 row 2) $\left|V\right|$: 14846 $\left|F\right|$: 29437	BMF BNF CNR GNF UNF L0 NLLR OURS	16.20 8.97 8.56 9.42 9.77 4.44 7.13 3.35	0.96 0.92 0.96 0.94 0.93 0.93 0.93 0.89	($30$) $(10,0.35,20)$ $(synthetic,3)$ $(10,0.35,20)$ $(0.5,20,10)$ $\left(Default\right)$ $\left(Default\right)$ ($3$)
Nicolo (${\sigma }_n=0.5l_e$) (Figure 9 row 3) $\left|V\right|$: 14846 $\left|F\right|$: 29437	BMF BNF CNR GNF UNF L0 NLLR OURS	56.33 65.33 41.35 58.15 61.66 5.67 38.44 10.33	1.13 1.00 1.22 0.93 0.94 1.02 0.92 1.08	($50$) $(10,0.45,20)$ $(synthetic,10)$ $(10,0.45,20)$ $(0.45,20,10)$ $\left(Default\right)$ $\left(Default\right)$ ($5$)
Vaselion (${\sigma }_n=0.2l_e$) (Figure 9 row 4) $\left|V\right|$: 38728 $\left|F\right|$: 77452	BMF BNF CNR GNF UNF L0 NLLR OURS	38.33 33.76 45.25 44.75 39.07 35.84 35.00 14.49	34.47 14.20 14.49 14.40 13.84 15.30 12.51 14.80	($30$) $(10,0.35,10)$ $(synthetic,1)$ $(10,0.35,10)$ $(0.5,10,10)$ $\left(Default\right)$ $\left(Default\right)$ ($3$)
Vaselion (${\sigma }_n=0.8l_e$) (Figure 9 row 5) $\left|V\right|$: 38728 $\left|F\right|$: 77452	BMF BNF CNR GNF UNF L0 NLLR OURS	167.9 171.9 126.0 166.6 158.7 14.72 82.25 32.35	30.06 12.58 14.82 13.98 13.26 16.49 11.29 15.87	($50$) $(10,0.45,30)$ $(synthetic,10)$ $(10,0.45,30)$ $(0.5,20,10)$ $\left(Default\right)$ $\left(Default\right)$ ($10$)
Boy (kinect fusion) (Figure 11 row 1) $\left|V\right|$: 76866 $\left|F\right|$: 152198	BMF BNF CNR GNF UNF L0 NLLR OURS	10.03 9.67 9.70 8.50 9.42 11.21 9.61 8.06	5.00 4.89 4.95 5.06 4.90 4.90 4.91 4.86	($30$) $(10,0.35,20)$ $(kinectfusion,10)$ $(10,0.35,20)$ $(0.5,10,10)$ $\left(Default\right)$ $\left(Default\right)$ ($10$)
David (kinect fusion) (Figure 11 row 2) $\left|V\right|$: 51789 $\left|F\right|$: 101937	BMF BNF CNR GNF UNF L0 NLLR OURS	10.74 10.70 15.85 9.14 10.28 12.59 10.39 9.32	1.60 1.572 2.09 1.65 1.566 1.60 1.58 2.12	($30$) $(10,0.45,20)$ $(kinectfusion,5)$ $(10,0.45,20)$ $(0.45,10,10)$ $\left(Default\right)$ $\left(Default\right)$ ($10$)
Pyramid (kinect fusion) (Figure 11 row 3) $\left|V\right|$: 35261 $\left|F\right|$: 69611	BMF BNF CNR GNF UNF L0 NLLR OURS	6.70 6.98 5.41 5.40 7.43 5.59 6.79 6.62	3.68 3.45 3.45 3.28 3.43 3.29 3.29 3.58	($30$) $(10,0.45,20)$ $(kinectfusion,3)$ $(10,0.45,20)$ $(0.6,20,10)$ $(\sqrt{2},0.1,0.1)$ $\left(Default\right)$ ($5$)

Table 3. Table 3: Runtime comparisons of several state of the art methods with our method (in seconds). The vertex iterations of all the methods are set to 10.

Models	UNF	BNF	GNF	L0	NLLR	OURS
Amadillo (Figure 8) $\left|V\right|$: 43243 $\left|F\right|$: 86482	5.68	7.12	117.5	55.8	75.9	1.81
iHbunny (Figure 9 Row 1) $\left|V\right|$: 34834 $\left|F\right|$: 69451	3.80	3.15	35.6	102.4	21.4	1.39
Nicolo (Figure 9 Row 2) $\left|V\right|$: 14846 $\left|F\right|$: 29437	0.89	0.83	28.0	23.8	19.0	0.59
Vaselion (Figure 9 Row 5) $\left|V\right|$: 38728 $\left|F\right|$: 77452	6.52	7.23	84.6	54.2	27.8	1.54
Sphere (Figure 10) $\left|V\right|$: 10242 $\left|F\right|$: 20480	0.39	0.54	8.97	16.1	15.1	0.38

Equations

$${\delta }_i=v_i-\frac{1}{|NV(v_i)|}\sum _{v_k\in NV(v_i)}{v}_k,$$

$$\frac{\partial V(x,t)}{\partial t}=\lambda △V(x,t),$$

$$△V={\nabla }^{2}V,$$

$$V^{t+1}=V^t+\lambda dt{\nabla }^2{V}^t,$$

$$\delta =(d\cdot n)n,$$

$$\mathbb{E}({\delta }_i^t)=\text{Arrowvert}{\delta }_i^t\text{Arrowvert}+\text{Arrowvert}{v}_i^t-v_i^0\text{Arrowvert},$$

$$E_v=\sqrt{\frac{1}{3{\sum }_{k\in F}{A}_k}\sum _{i\in V}\sum _{j\in F_V(i)}{A}_{j}dist(x_i^{{}^{\prime }},T{)}^2},$$

$$\text{MSAE}=\frac{{\sum }_{k\in F}{\theta }_k^2}{N_F}$$

Full text

1 Introduction

Surface meshes acquired through scanning or sensing equipment are inevitably contaminated with noise. The user has to process them with smoothing techniques before applying to down stream applications, such as shape analysis, animation and rendering. The main technical challenge of surface smoothing is preserving features while removing noise and alleviating shrinkage. The design of robust surface smoothing methods is therefore of particular need in nowadays.

The discrete Laplacian operator on surface meshes has proven highly useful in various tasks of digital geometry processing, for example, mesh fairing, parameterization, reconstruction, editing and compressing [1, 2, 3, 4, 5, 6], just name a few. The differential surface representation encodes information about the local shape of a surface such as the curvature and the orientation. Despite that the uniform Laplacian operator can effectively smooth surfaces, it fails to preserve features and leads to shrinkage.

Among the feature-preserving methods, the majority of them typically use the normal information of surfaces, such as face or vertex normals. In contrast, very few works are based on vertex position [7, 8]. Some recent works extended these two types by adding more delicate steps [9, 10, 11, 12, 13, 8, 14, 15]. However, these methods generally consist of multiple steps and numerous parameters, which are difficult to use/tune especially for complex models. Moreover, the users especially those out of the field may find it difficult and tedious to tune the involved parameters.

To overcome the above issues, we propose a novel, robust approach for feature-preserving mesh smoothing in this paper. Our key idea is to construct a “half-kernel” uniform Laplacian operator (HLO) that can approximate the Laplacians at feature and non-feature points using half windows. Specifically, we first analyze the differential property of feature points and found it conflicts with the existence assumption of the uniform Laplacian. We then naturally propose a half-window algorithm to generate different pairs of subsets (half windows) for each vertex by pairing one immediate neighbor to the other unique neighbor. We compute the Laplacians of all subsets which are further projected to the full-window Laplacians to alleviate flipping and degeneration. The final half-kernel Laplacian is automatically determined as the one that incurs the smallest regularization energy. The surface is finally updated with the determined half-kernel Laplacians in an iterative way.

Taking a noisy surface mesh as input, our approach can automatically output a quality version with preserving features and resisting shrinkage. The main contributions of this work are:

mathematical analysis of the differential property at feature points;

• 2.

a half-window algorithm to generate multiple half windows (subsets) for each vertex;

• 3.

a half-kernel uniform Laplacian operator (HLO) for feature-preserving and shrinkage-resisting surface smoothing.

Our method††footnotemark: is conceptually simple and easy to use since it involves only a single parameter (i.e., the number of vertex update iterations). We demonstrate that the proposed HLO substantially outperforms the uniform Laplacian in preserving features and resisting shrinkage. We evaluate our method on both synthetic and real-world models, and observe that our method can produce results with comparable or higher quality than the state-of-the-art methods. We also show that our method is faster than the state-of-the-art methods.

2 Related Work

There exists a large body of literature of mesh denoising. Due to space limit, we review only the mostly related works to ours and refer the readers to the comprehensive surveys [16, 17].

The classical Laplacian smoothing methods [18, 19] are simple and fast. However, its isotropic property leads to feature-wiping and shrinking artifacts. Taubin [20] proposed a non-shrinking, two-step smoothing method with positive and negative damping factors. Desbrun et al. [21] proposed a fairing method based on diffusion and curvature flow to process irregular meshes. Later, various isotropic smoothing methods have been introduced based on volume preservation, pass frequency controlling, differential properties, etc,. [22, 23, 24, 25, 26].

The above isotropic methods are effective to remove noise, however they also wipe out features. Various anisotropic methods have been proposed to preserve features. Representative methods are diffusion/differential-based methods [27, 28, 29, 30, 31, 32, 33, 34, 7], bilateral filters [35, 36, 37, 38, 39], methods combining normal filtering and vertex update [37, 38, 40, 41, 42, 43, 44, 45, 46, 35, 47, 48, 43, 49, 14, 15]. Bilateral filtering is initially used in image denoising [50], and it is successfully extended to mesh denoising. Regarding the feature-preserving techniques, the utilization of filtered normals has seen noticeable progress in recent years [37, 38, 40, 41, 42, 43, 44, 45, 46, 35, 47, 48, 43, 49, 14, 15, 51, 52]. Some recent methods have focused on vertex and face classification before mesh denoising [9, 10, 11, 12, 13, 53, 54]. Nevertheless, the classification results depends largely on the level of noise. Lu et al. [8, 14, 15] presented the ideas of pre-filtering before real smoothing and can robustly handle heavy noise with frequent flipped triangle faces. Arvanitis et al. [55] introduced a novel coarse-to-fine graph spectral processing approach for mesh denoising.

Another line of anisotropic methods focused on the sparse perspective. Compared with non-feature vertices, feature vertices in a mesh are usually sparse, which can be reconstructed or detected by solving a sparse problem [56]. Sparsity was introduced into mesh smoothing in some recent works [7, 57, 8, 58]. For instance, He et al. [7] developed a $L_{0}$ minimization framework with an area-based edge operator which is generally sparse in a surface mesh. It is, however, non-convex and difficult to solve. Zhao et al. [58] designed an improved alternating optimization strategy to solve the $L_{0}$ minimization which incorporates both vertex positions and face normals. Wang et al. [57] proposed a method to decouple noise and features by weighted $L_{1}$-analysis compressed sensing, and they prove that the pseudo-inverse matrix of the Laplacian of a mesh is a coherent dictionary for sparsely representing sharp feature signals on the shape. Recently, Lu et al. [8] detected features by introducing a novel $L_{1}$ minimization. Lu and his colleagues [59] proposed a low-rank matrix approximation approach for geometry filtering and demonstrated various geometry processing applications. More recently, the low-rank optimization was further extended to mesh denoising [60, 61].

3 Half-kernel Laplacian Operator

To better understand the proposed half-kernel Laplacian operator (HLO), we briefly introduce the uniform Laplacian on meshes and the Laplacian diffusion flow in the first place. We then analyze the differential property of feature points and finally explain how we construct the HLO.

3.1 Uniform Laplacian

Given a surface mesh $M=(V,E,F)$ with $N$ vertices, we have the set of vertices $V$, the set of edges $E$ and the set of faces $F$. The $i$-th vertex $v_{i}\in V$ is represented by the coordinates $v_{i}=(x_{i},y_{i},z_{i})$. The differential coordinates (i.e., $\delta$-coordinates) of vertex $v_{i}$ are defined as

[TABLE]

where $\delta_{i}=(\delta_{x},\delta_{y},\delta_{z})$, $NV(v_{i})$ is the set of the neighboring vertices of the vertex $v_{i}$, and $|NV(v_{i})|$ is the degree (number of neighbors) of $v_{i}$. It is also called the uniform Laplacian due to the equal weights. The direction of the differential coordinates approximates the local normal direction and the magnitude linearly approximates the local mean curvature $H(v_{i})$ of $v_{i}$ [4].

3.2 Laplacian Diffusion Flow

Fairing a mesh can be viewed as a filtering process on mesh signal. In the discretized setting, it is usually done by solving a heat-diffusion-like partial differential equation as follows.

[TABLE]

where $\lambda$ is the diffusion speed. The Laplacian operator is defined as

[TABLE]

where $\nabla^{2}V=\nabla\cdot\nabla$ is the divergence of the gradient on the vertices of a given mesh. The above equation (Eq. (2)) is usually solved in an iterative way

[TABLE]

where $dt$ is the time step length.

The surface smoothing based on the uniform Laplacian often causes the issues of shrinkage and smoothing out features (Figure 2a). Motivated by these issues, we attempt to first analyze the differential property of feature points, which will be described in Sec. 3.3.

3.3 Differential Analysis of Feature Points

Intuitively speaking, feature points can be classified as the points with abrupt normals (i.e., with large dihedral angles) [14, 7]. Denote by $f$ a surface manifold and $a$ a feature point at a sharp edge on the cube model (Figure 1(a)). We simply analyze the situation when $f$ crosses the sharp edge (Figure 1(a)). It is obvious that $f^{\prime}(a+\epsilon u)\neq f^{\prime}(a+\epsilon v)$, where $f^{\prime}$ means the first-order derivative of $f$ and $\epsilon>0$. Thus, $f$ is not differentiable at the feature point $a$, which means that the second-order derivative $f^{\prime\prime}(a)$ does not exist either. This reveals that it does not make sense to apply the uniform Laplacian (i.e., second-order derivative) operator to a feature vertex with all its neighbors (i.e., full window, see Figure 1(a)). Figure 2(a) shows the feature-wiping and shrinking results with the uniform Laplacian operator in a full-window sense.

Ideally, the approximations for locations ($a+\epsilon v$) and ($a+\epsilon u$) should come from the upper region and the front region, respectively. The approximation for $a$ should base on its two neighboring vertices on the edge since it is differentiable along the edge [57]. Corner, the intersection of several edges, is theoretically non-differentiable and usually fixed in mesh denoising [57]. Sharp edge points and corners are difficult to approximate because they are identified with some means which makes surface smoothing less robust due to sensitivity to the noise level. To overcome this difficulty, we attempt to approximate each (feature or non-feature) vertex with a half-kernel Laplacian operator which targets at computing the Laplacians using half-window neighbors. It will be discussed in Section 3.4.

3.4 The Half-kernel Laplacian Operator

The above observation and analysis indicate that the uniform Laplacian operator should be performed only on half windows (i.e., subsets of neighbors) of the feature point instead of its whole neighborhood. This key insight encourages us to introduce a half-kernel Laplacian operator (HLO) for surface meshes.

Our operator is to approximate the Laplacians at vertices with half windows. To produce half windows, each of the immediate neighboring vertex of the current vertex is paired with another neighboring vertex to partition the local neighborhood into two half windows (left and right). The other neighbor paired by the starting neighbor is determined as the neighboring vertex which has the shortest distance to the plane defined by the current vertex, the starting neighbor and the centroid of the neighborhood. Refer to Figure 4(a). Consider a vertex $v_{0}$ with a neighbor $v_{D}$. Denote by $v_{A}$ the starting neighbor to be paired. $v_{0}^{\prime}$ is the centroid of the neighborhood of $v_{0}$. Figure 1(b) also shows a half window of a feature vertex. As a result, each vertex $v_{i}$ and its neighborhood has $|NV(v_{i})|$ partition choices which further generate $2|NV(v_{i})|$ paired subsets (half windows). The half-window generation algorithm is shown in Algorithm 1.

Remark 1. Theoretically, the full window neighbors can be split into more than two half windows. We simply select the choice of two half windows, because corners are typically far fewer than the edge features in a shape, and recognizing corners may introduce new parameters and result in less robustness. Also, this half-window scheme is sufficient for non-feature vertices.

Remark 2. In general, there might exist multiple candidates with the same shortest distances when determining the other neighbor. We randomly select one candidate in this case, which we found no noticeable differences in mesh denoising. Note that the plane defined by the current vertex, the starting neighbor and the centroid of the neighborhood may degenerate into a line. We simply determine the other neighbor using the distance to the line instead.

We apply the uniform Laplacian operator independently to each of these subsets and compute the half-window Laplacians for each vertex $v_{i}$. Because half-kernel Laplacians introduce a shift in tangential direction, directly using the half-kernel Lapalcians may result in inferior results like degenerating and flipping (Figure 3).

We thus project them onto their corresponding full-kernel Laplacians (i.e., using full neighbors) to remove the tangential component and obtain the final half-kernel Laplacians. This way ensures small Laplacian magnitudes for feature vertices and large Laplacian magnitudes for non-feature vertices, thus smoothing non-features and preserving features. Figure 4(b-c) shows two such examples. We define the projection in Eq. (5).

[TABLE]

where $d$ indicates the half-kernel Laplacian and $n$ is the orientation of the full-kernel Laplacian.

We can easily enumerate all the intermediate half-kernel Laplacians for $v_{i}$: $\dots,\delta_{L}(v_{k}),\delta_{R}(v_{k}),\dots$, where $\delta_{L}(v_{k})$ and $\delta_{R}(v_{k})$ respectively denote the half-kernel Laplacians for the left subset and right subset with the starting neighbor $v_{k}$ ($v_{k}\in NV(v_{i})$). We next need to determine the optimal Laplacian among the $2|NV(v_{i})|$ half-kernel Laplacians. The optimal one will be selected based on the regularization energy which may vary in different applications. In this work, we define the regularization energy of each vertex as the sum of the norm of the Laplacian and the the distance to the original vertex position, shown as below.

[TABLE]

where $\delta_{i}^{t}$ can be any among the calculated $2|NV(v_{i})|$ half-kernel Laplacians in the $t$-th iteration. $v_{i}^{t}$ is the position of vertex $v_{i}$ in the $t$-th iteration and $v_{i}^{0}$ is the initial position of $v_{i}$. The former term, which naturally uses the computed $\delta_{i}^{t}$ and is a feature-preserving term, represents the movement between two consecutive iterations, and the latter term, a typical data term, describes the distance between the position at the $t$-th iteration and the initial position of vertex $v_{i}$. This energy enables a closest walking from the initial positions and positions in the previous iteration, thereby preserving features and the original shape. We compute the energy $\mathbb{E}$ for each half-kernel Laplacian of vertex $v_{i}$ and select the optimal Laplacian that incurs the smallest energy. Algorithm 2 summarizes the proposed half-kernel Laplacian operator. To smooth surfaces, the final half-kernel Laplacians in the $t$-th iteration are used to update vertex positions in the $(t+1)$-th iteration via Eq. (4). $\lambda dt$ in Eq. (4) and the number of iterations work in a proportional way to each other. In other words, increasing $\lambda dt$ would generally induce a decrease of the number of update iterations, and vice versa. A larger $\lambda dt$ would also possibly cause instability. To make our method more robust, we empirically set $\lambda dt$ to 1, which works very well in all our experiments. Figure 2(b) shows an example for our HLO.

4 Experimental Results

We first describe the parameter setting, and then compare the proposed HLO with the uniform Laplacian operator (ULO). Finally, We compare our method with the state-of-the-art denoising methods on both synthetic and scanned data, visually and quantitatively.

4.1 Parameter Setting

The compared state-of-the-art techniques are: the bilateral mesh filter (BMF) [36], the unilateral normal filter (UNF) [45], the bilateral normal filter (BNF) [37], the $L_{0}$ minimization method ($L_{0}$) [7], the guided normal filter (GNF) [43], the cascaded normal regression (CNR) [49], and the non-local low-rank normal filtering method (NLLR) [60]. The source code of these methods are available or the authors provide test results. We summarize the parameters of these methods in Table 1. Table 2 shows the parameter values of all methods on most models. Interested readers are referred to the papers for more specific details.

4.2 Comparison with the Uniform Laplacian

Figure 2 shows results processed by the uniform Laplacian and our HLO, respectively. We listed the corresponding results after 1, 3, 5, and 15 iterations. It is obvious that our HLO significantly alleviates the shrinking effect. The feature places (e.g., upper and lower parts) of the vase are preserved by HLO.

Figure 5 visualizes the one-to-one vertex errors of both the uniform Laplacian and our HLO on the Julius model. Our method better preserves surface features and produces more accurate smoothing results than the uniform Laplacian operator. The average vertex errors and mean curvature energy [62, 4, 63] of Figure 5 are displayed in Figure 6. The HLO induces much lower average vertex errors than the uniform Laplacian operator, which means it preserves the shape better. The first iteration sees a remarkable decrease on the mean curvature energy, and the two methods have a similar decreasing trend in mean curvature led by diffusion flow.

We also compared with the cotangent Laplacian. Figure 7 shows the results of the uniform Laplacian, cotangent Laplacian and our HLO. The cotangent Laplacian is more geometry aware than the uniform Laplacian. However, the cotangent Laplacian tends to generate sharp “thorns” on the surface (please zoom in to observe).

4.3 Visual Results

Synthetic models. We compare our method with the selected state-of-the-art methods on various models corrupted with synthetic noise. Following state-of-the-art mesh smoothing techniques, we generate synthetic models by adding zero-mean Gaussian noise with standard deviation $\sigma_{n}$ to the corresponding ground truth. $\sigma_{n}$ is proportional to the mean edge length $l_{e}$ of the input mesh.

The uniform Laplacian operator can unfold the flipped triangles during surface smoothing. Thanks to this property, our approach is less sensitive to high-level noise than most of the compared methods, as shown in Figure 8 and 9 (Bunny: $\sigma_{n}=0.5l_{e}$, Nicolo: $\sigma_{n}=0.5l_{e}$, Vaselion: $\sigma_{n}=0.8l_{e}$ ). *The flipped triangles are rendered in red. See the close-up views. *

Figure 8 shows visual results over the Armadillo model contaminated by Gaussian noise with $\sigma_{n}=0.5e_{l}$. We can see from the blown-up windows that the state-of-the-art methods tend to oversmooth or oversharpen the fine details, or retain excessive noise in the model. We can observe from the mouth and eyes that our method well preserves the small-scale features which are usually lost to some extent in the results by the other methods. Figure 9 (Bunny:$\sigma_{n}=0.5l_{e}$ and Nicolo:$\sigma_{n}=0.5l_{e}$) demonstrates that our approach preserves better surface details. Some methods may oversharpen the ear in Figure 9 (Bunny:$\sigma_{n}=0.5l_{e}$) and nose in Figure 9 (Nicolo:$\sigma_{n}=0.5l_{e}$). In Figure 9 (Vaselion:$\sigma_{n}=0.8l_{e}$), there are many curved features which are particularly difficult to recover by most existing methods when removing high-level noise. The method [7] can robustly remove the noise while it results in oversharpening in certain areas and loses some fine details. By contrast, our method outputs better results, in terms of features and details preservation. In Figure 10, our method well recovers the sphere shape which has fewer folded triangles than other methods. Notice that CNR [49] may produce better results if the model is trained on surface meshes with large noise.

Undoubtedly, our method can also well handle small or medium noise. Figure 9 (Nicolo:$\sigma_{n}=0.2l_{e}$ and Vaselion:$\sigma_{n}=0.2l_{e}$) shows two such examples. We can observe from Figure 9 (Nicolo:$\sigma_{n}=0.2l_{e}$) that our method and NLLR [60] generates similar results which are the best among all results. NLLR [60] uses non-local information for mesh smoothing while our method uses local information only. Figure 9 (Vaselion:$\sigma_{n}=0.2l_{e}$) demonstrates that the result by our HLO is better than the results by other techniques, in terms of noise removal and features preservation.

Raw scanned models. In addition to the synthetic models, we also tested the proposed HLO on scanned models corrupted with raw noise. Figure 11 show the smoothing results of all methods on three real scanned surface meshes. Our method produces very competitive results, in terms of preserving features. Taking the second row as example, our HLO and NLLR [60] produce the best results while other methods may oversharpen or oversmooth certain regions. Note that NLLR [60] uses non-local information while ours utilizes local information only.

There are several reasons for why our method outperforms or is comparable to state-of-the-art techniques. (1) Our method is position based which largely mitigates the issue of flipped triangles. By contrast, most normal-based approaches (e.g., [45, 37, 60]) cannot effectively overcome this issue, without using any pre-processing strategy like [14, 8]. [7] is also position based; however it may over-sharpen some regions unexpectedly, due to the edge-based operator on the whole shape. CNR [49] is a normal regression method that relies on immediate position update in each iteration; that is, it cannot directly predict the final normals. It thus easily smooths out features without a good mapping. (2) Multiple steps and parameters make those methods more complicated to search the solution space, thus leading to less automation and less robustness in handling different levels of noise, especially heavy noise. On the contrary, our method involves a single parameter, and is more simpler, automatic and robust in dealing with various noise.

4.4 Quantitative Evaluations

Besides the above visual comparisons, we also compare the state-of-the-art techniques with our approach from a quantitative perspective. Specifically, we employ $E_{v}$ and MSAE (mean square angular error) to respectively evaluate the positional error and normal error, as suggested by previous works [45, 37, 8]. These two metrics are calculated between the smoothing results and their corresponding ground truth.

According to  [45], $E_{v}$ is the $L^{2}$ vetex-based mesh-to-mesh error metric, and MSAE measures the mean square angular error between the face normals of the denoised mesh and those of the ground truth.

[TABLE]

where $A_{k}$ is the area of face $k$, and $dist(x^{{}^{\prime}}_{i},T)$ is the $L^{2}$ distance between the updated vertex $x^{{}^{\prime}}_{i}$ and a triangle of the reference mesh $T$ which is closest to $x^{{}^{\prime}}_{i}$.

[TABLE]

where $\theta_{k}$ is the angle between the $k$-th face normal of the denoised model and its corresponding normal in the ground-truth model, and $N_{F}$ is the number of faces in the 3D shape.

Table 2 summarizes the statistical numbers of $E_{v}$ and MSAE over most models for all the compared methods. As with previous works [45, 8], we also found that the visual comparison results might not necessarily agree with the $E_{v}$. This is because our method smooths surfaces via unfolding the flipped triangles in the noisy input, which is similar to the behavior of [7]. As a consequence, vertices may walk far from their ground-truth locations. Regarding MSAE, the $L_{0}$ minimization [7] and our method are usually the best two among all techniques. We also compare the runtime of our method with those methods and list them in Table 3. Our approach outperforms those methods in terms of speed. It has excellent speed performance in handling models with large sizes (e.g., Amadillo and Vaselion in Figure 10). It should be noted that our code has not been optimized or speeded up via parallelization.

5 Conclusion

In this work, we have introduced a novel method for feature-preserving surface smoothing. Motivated by the shrinking and feature-wiping issues of the uniform Laplacian operator, we analyzed the differential property at feature points. We developed a Half-kernel Laplacian Operator (HLO) which can preserve surface features and resist shrinkage. Various experiments show that our approach is better or comparable to state-of-the-art surface smoothing techniques, in terms of visual quality and quantitative evaluations. Our method is robust to high noise and it has a single parameter (the number of vertex update iterations) which is easy and intuitive to tune.

Our method still suffers from a few limitations. First, it has difficulty to recover smooth sharp edges for CAD-like models with a small number of vertices (Figure 13). This is because HLO is not an edge-based operator [7]. For real scanned CAD-like models with large number of vertices, our method can produce results with comparable quality to the state-of-the-art methods (Figure 11 row 3). As with previous works [37, 8], we also fixed the open boundaries of the input noisy mesh, which might sometimes lead to undesired results.

As the future work, we would like to incorporate new schemes to overcome the issue of recovering smooth sharp edges. We would also like to improve the performance via parallelization.

Acknowledgements

This work was supported in part by the National Natural Science Foundation of China under Grant 61907031. Xuequan Lu is supported by Deakin CY01-251301-F003-PJ03906-PG00447. Ying He is supported by MOE RG26/17.