# 3D Point Cloud Denoising using Graph Laplacian Regularization of a Low Dimensional Manifold Model

3D point cloud - a new signal representation of volumetric objects - is a discrete collection of triples marking exterior object surface locations in 3D space. Conventional imperfect acquisition processes of 3D point cloud - e.g., stereo-matching from multiple viewpoint images or depth data acquired directly from active light sensors - imply non-negligible noise in the data. In this paper, we adopt a previously proposed low-dimensional manifold model for the surface patches in the point cloud and seek self-similar patches to denoise them simultaneously using the patch manifold prior. Due to discrete observations of the patches on the manifold, we approximate the manifold dimension computation defined in the continuous domain with a patch-based graph Laplacian regularizer and propose a new discrete patch distance measure to quantify the similarity between two same-sized surface patches for graph construction that is robust to noise. We show that our graph Laplacian regularizer has a natural graph spectral interpretation, and has desirable numerical stability properties via eigenanalysis. Extensive simulation results show that our proposed denoising scheme can outperform state-of-the-art methods in objective metrics and can better preserve visually salient structural features like edges.

• 6 publications
• 26 publications
• 5 publications
• 14 publications
• 44 publications
02/20/2019

### Point cloud denoising based on tensor Tucker decomposition

In this paper, we propose an algorithm for point cloud denoising based o...
10/14/2020

### PointManifold: Using Manifold Learning for Point Cloud Classification

In this paper, we propose a point cloud classification method based on g...
02/29/2012

### Perturbation of the Eigenvectors of the Graph Laplacian: Application to Image Denoising

The original contributions of this paper are twofold: a new understandin...
04/27/2016

### Graph Laplacian Regularization for Image Denoising: Analysis in the Continuous Domain

Inverse imaging problems are inherently under-determined, and hence it i...
11/09/2021

### Graph-Based Depth Denoising Dequantization for Point Cloud Enhancement

A 3D point cloud is typically constructed from depth measurements acquir...
09/15/2021

### Graph skeletonization of high-dimensional point cloud data via topological method

Geometric graphs form an important family of hidden structures behind da...
07/22/2019

### Feature Graph Learning for 3D Point Cloud Denoising

Identifying an appropriate underlying graph kernel that reflects pairwis...

## I Introduction

The three-dimensional (3D) point cloud has become an important and popular signal representation of volumetric objects in 3D space [1]. 3D point cloud can be acquired directly using low-cost depth sensors like Microsoft Kinect or high-resolution 3D scanners like LiDAR. Moreover, multi-view stereo-matching techniques have been extensively studied in recent years to recover a 3D model from images or videos, where the typical output format is the point cloud [2]. However, in either case, the output point cloud is inherently noisy, which has led to numerous approaches for point cloud denoising [3, 4, 5, 6].

Moving least squares (MLS)-based [7, 8] and locally optimal projection (LOP)-based methods [9, 10] are two major categories of point cloud denoising approaches, but are often criticized for over-smoothing [5, 6] due to the use of local operators. Sparsity-based methods [5, 4] assume sparse representation of the point normals and achieve good performance [11]

. However at high noise levels, normal estimation can be so poor that it leads to over-smoothing or over-sharpening

[5].

Non-local methods generalize the ideas in the non-local means [12] and BM3D [13] image denoising algorithms to point cloud denoising, and are shown to achieve state-of-the-art performance and robust to high level of noise [3, 4]. These methods exploit the inherent self-similarity characteristic between surface patches in the point cloud to preserve structural details, but their computational complexity is too high to be practical [3].

Utilizing an assumed self-similarity characteristic in images has long been a popular strategy in image processing [12, 13]. Extending on these earlier works, a more recent work [14] proposed a new low-dimensional manifold model (LDMM)—similar image patches are samples of a low-dimensional manifold in high-dimensional space—and achieved state-of-the-art results in various inverse imaging applications, e.g., denoising, inpainting, superresolution, etc..

Inspired by the work in [14], we exploit self-similarity of the surface patches by assuming that the surface patches in the point cloud lie on a manifold of low dimension. In the optimization, the manifold dimension serves as a regularization prior, which is approximated by a discrete graph Laplacian regularizer due to discrete observation of the surface patches. The main contributions of our work are as follows:

1. By adopting the LDMM, we exploit the surface self-similarity characteristic and simultaneously denoise similar patches to better preserve sharp features;

2. The computation of the manifold dimension is approximated by a discrete graph Laplacian regularizer with low computational complexity; and

3. An efficient similarity measure for discrete -pixel patch pairs is designed for graph construction that is robust to noise.

In addition, we show that our graph Laplacian regularizer has a graph spectral interpretation, and has desirable numerical stability properties via eigen-analysis in the graph spectral domain [15]. Extensive simulation results show that our proposed method can outperform the state-of-the-art methods in objective metrics and can better preserve visually salient features like edges.

The rest of the paper is organized as follows. Section II overviews some existing works. Section III defines the patch manifold associated with the 3D point cloud. Section IV formulates the denoising problem by describing how the manifold dimension is computed and approximated with the graph Laplacian regularizer. The algorithm implementation is discussed in Section V with graph spectral analysis to interpret the algorithm and a numerical stability analysis. Finally, Section VI and Section VII presents experimental results and concludes the paper respectively.

## Ii Related Work

Previous point cloud denoising works can be classified into four categories: moving least squares (MLS)-based methods, locally optimal projection (LOP)-based methods, sparsity-based methods, and non-local similarity-based methods.

MLS-based methods. MLS-based methods approximate a smooth surface from the input samples and project the points to the resulting surface. To construct the surface, the method in [16] first finds the local reference domain for each point that best fits its neighboring points in terms of MLS, then defines a function above the reference domain by fitting a polynomial function to neighboring data.

Several extensions, which address the unstable reconstruction problem in the case of high curvature, e.g., algebraic point set surfaces (APSS) [7] and its variant in [17], or preserve the shape features, e.g., robust MLS (RMLS) [18] and robust implicit MLS (RIMLS) [8], have also been proposed. These methods can robustly generate a smooth surface from extremely noisy input, but are often criticized for over-smoothing [5, 6].

LOP-based methods. Unlike MLS-based methods, LOP-based methods do not compute explicit parameters for the surface. For example, LOP method in [9]

outputs a set of points that represent the underlying surface while enforcing a uniform distribution over the point cloud with a repulse term in the optimization. Its modifications include weighted LOP (WLOP)

[19], which provides a more uniformly distributed output by adapting the repulse term to the local density, and anisotropic WLOP (AWLOP) [10], which preserves sharp features by modifying WLOP to use an anisotropic weighting function. LOP-based methods also suffer from over-smoothing due to the use of local operators, or generate extra features caused by noise [5, 6].

Sparsity-based methods. Sparsity-based methods assume sparse representation of the point normals. A sparse reconstruction of the point normals is obtained by solving a global minimization problem with sparsity regularization. Then the point positions are updated by solving another global (or ) minimization problem based on a locally planar assumption. Some of the examples include [20] which utilizes regularization, and [5] and [4] which adopt regularization and achieve state-of-art performance [11]. However, when the noise level is high, normal estimation can be so poor that it leads to over-smoothing or over-sharpening [5].

Non-local methods. Non-local methods generalize the ideas in the non-local means [12] and BM3D [13] image denoising algorithms to point cloud denoising. These methods exploit the self-similarities between surface patches in the point cloud. Some early trials are [21] and [22], which utilize a non-local means algorithm, while more recent methods in [3, 4] inspired by the BM3D algorithm are shown to outperform the non-local means methods. However, the computational complexity is typically too high to be practical.

Our method belongs to the fourth category, the non-local methods. Similar to [3, 4], we also utilize the self-similarity among patches via the low-dimensional manifold prior [14]. However, unlike [14] that requires time-consuming fixed point integration, due to the use of graph Laplacian regularizer, our proposed approach can be efficiently implemented and capable of preserving sharp feature better than existing schemes.

## Iii Patch Manifold

We first define a notion of patch manifold given a point cloud , , which is a (roughly uniform) discrete sampling of a 2D surface of a 3D object. Let be the position matrix for the point cloud. Noise-corrupted can be simply modeled as:

 V=U+E, (1)

where contains the true 3D positions, is a zero-mean signal-independent noise (we assume Gaussian noise in our experiments), and . To recover the true position , we consider the low-dimensional manifold model prior (LDMM) [14] as a regularization term for the ill-posed problem.

### Iii-a Surface Patch

We first define a surface patch in a point cloud. We select a subset of points from as the patch centers, i.e., . Then, patch centered at a given center is defined as the set of nearest neighbors of in , in terms of Euclidean distance.

The union of the patches should cover the whole point cloud, i.e.,

. There can be different choices of patch centers, and the degree of freedom can be used to trade off computation cost and denoising performance. Let

be the patch coordinates, composed of the points in .

### Iii-B Patch Manifold

Here we adopt the basic assumption in [14] that the patches sample a low-dimensional manifold embedded in , which is called the patch manifold associated with the point cloud . In order to evaluate similarity among patches, we first need to align the patches; i.e., the coordinates should be translated with respect to , so that lies on the origin point . Here we set to be the translated coordinates.

### Iii-C Low Dimensional Patch Manifold Prior

The LDMM prior assumes that the solution contains patches that minimize the patch manifold dimension. We can thus formulate a maximum a posteriori (MAP) problem with prior and fidelity terms as follows:

 minUdim(M(U))+λ∥V−U∥2F, (2)

where is a parameter that trades off the prior with the fidelity term, and is the Frobenius norm. Note that given a certain strategy of patch selection, the patches are determined by the point cloud , and the patches in turn define the underlying manifold . Hence we view as a function of .

The patches can be very different and sampled from different manifolds of different dimensions. For example, a flat planar patch belongs to a manifold of lower dimension than a patch with corners. The dimension of the patch manifold, becomes a function of the patch, and the integration of over is used as the regularization term,

 minU∫Mdim(M(U))(p)dp+λ∥V−U∥2F, (3)

where is the dimension of at . Here is a point on .

## Iv Problem Formulation

In this section, we discuss the calculation of the manifold dimension, approximate this computation with the graph Laplacian regularizer, and design a patch distance measure for the graph construction on the manifold.

### Iv-a Manifold Dimension Computation

To solve (3), we first need to compute the manifold dimension . Here we overview how the manifold dimension is computed in [14].

First, let , where , be the coordinate functions on the manifold embedded in , i.e.,

 αi(p)=pi,  ∀p=[p1,…,p3k]⊤∈M. (4)

According to [14], the dimension of at is given by:

 dim(M)(p)=3k∑i=1∥∇Mαi(p)∥2, (5)

where denotes the gradient of the function on at . Then the integration of over is given as,

 (6)

### Iv-B Discretization of Manifold Dimension

The formula in (6) is a sum of integrals on continuous manifold along different dimensions, but our observations of the manifold are discrete and finite. Hence, we approximate (6) using a discrete graph Laplacian regularizer.

Next, we first introduce the graph construction on a manifold, which induces the graph Laplacian regularizer. Then we discuss how the graph Laplacian regularizer approximates the manifold dimension.

#### Iv-B1 Constructing Graph on a Manifold

We construct a discrete graph whose vertex set is the observed surface patches lying on , i.e., . Let denote the edge set, where the edge between -th and -th patches is weighted as,

 wmn=(ρmρn)−1/γψ(dmn). (7)

The kernel is a thresholded Gaussian function

 ψ(dmn)={exp(−d2mn2ϵ2)dmn

and is the Euclidean distance between the two patches and ,

 dmn=||pm−pn||2. (8)

The term is the normalization term, where is the degree of before normalization. The graph constructed in these settings is an -neighborhood graph, i.e., no edge has a distance greater than . Here , and is a constant.

#### Iv-B2 Graph Laplacian Regularizer

With the edge weights defined above, we define the symmetric adjacency matrix , with the -th entry given by . denotes the diagonal degree matrix, where entry . The combinatorial graph Laplacian matrix is [15].

For the coordinate function on defined in (4), sampling at positions of leads to its discretized version, . The graph Laplacian induces the regularizer . It can be shown that

 SL(αi)=α⊤iLαi=∑(m,n)∈Ewmn(αi(pm)−αi(pn))2. (9)

#### Iv-B3 Approximation with Graph Laplacian Regularizer

We now show the convergence of the discrete graph Laplacian regularizer to the dimension of the underlying continuous manifold.

Assume that

is a Riemannian manifold with boundary, equipped with the probability density function (PDF)

describing the distribution of the vertices on , and that belongs to the class of -Hlder functions [23] on . Then according to the proof in [24], there exists a constant depending only on such that for and the weight parameter , where denotes the manifold dimension, such that111We refer readers to [23] for the uniform convergence result in a more general setting and its corresponding assumptions on , , and the graph weight kernel function .

 sup∣∣∣cM2γ−1ϵ4(1−γ)(M−1)SL(αi)−SΔ(αi)∣∣∣=O(M−κ2κ+2δ+δ2+δκ), (10)

where is induced by the -th weighted Laplace-Beltrami operator on , which is given as

 SΔ(αi)=∫M∥∇Mαi(p)∥22h(p)2(1−γ)dp. (11)

Assuming that the vertices are uniformly distributed on , then

 ∫Mh(p)dp=1,h(p)=1|M|, (12)

where is the volume of the manifold . For implementation, similar to the setting in [24], we set , then becomes

 SΔ(αi)=1|M|∫M∥∇Mαi(p)∥22dp. (13)

From (10) and (13), the convergence is readily obtained by weakening the uniform convergence of (10) to point-wise convergence:

 limM→∞,ϵ→0|M|3k∑i=1α⊤iLαi∼3k∑i=1∫M∥∇Mαi(p)∥22dp, (14)

where means there exists a constant depending on , and , such that the equality holds. The rate of convergence depends on and . As the number of samples increases and the neighborhood size shrinks, the graph Laplacian regularization approaches its continuous limit. Moreover, if the manifold dimension is low, we can ensure a good approximation of the continuous regularization functional even if the manifold is embedded in a high-dimensional space. Consequently, given a point cloud, one can approximate the dimension of with the ’s and the constructed graph Laplacian .

So far, we obtain the approximation of (6) with (14), but the above graph construction requires the patch coordinates to be ordered so that the ’s can be defined and the patch distance in (8) determines the patch similarity. For example, if the patches are image patches of the same size, then the patch coordinates are naturally ordered according to pixel location, i.e., the -th entry in is the pixel value at the -th location in the image patch. However, surface patches in the 3D point cloud are unstructured, and there is no natural way to implement global coordinate ordering for all patches.

In the following, we first argue that the computation of the regularization term can be accomplished based on local pairwise correspondence between connected patches, relieving the need for global ordering. Then, we propose a patch distance measure for unstructured point cloud that is efficiently computed and robust to noise.

### Iv-C Patch Distance Measure

#### Iv-C1 From Global Coordinate Ordering to Local Correspondence

We modify the manifold dimension formula in (14):

 3k∑i=1α⊤iLαi =3k∑i=1∑(m,n)∈Ewmn(αi(pm)−αi(pn))2 (15) =∑(m,n)∈Ewmn3k∑i=1(αi(pm)−αi(pn))2 (16) =∑(m,n)∈Ewmnd2mn, (17)

where (17) follows from (16) according to the definition of in (8). From (17) we see that is not necessary to compute the graph Laplacian regularizer, and hence global coordinate ordering is not required. Moreover, since is itself a function of via (7), we can instead find the local pairwise correspondence between neighboring patches to compute patch distance in (17).

#### Iv-C2 Distance Measure in Continuous Domain

To measure the distance between the -th patch and -th patch, ideally

the two patches can be interpolated to two continuous surfaces, and the distance is calculated as the integral of the surface distance over a local domain around the patch center.

To define the underlying surface, we first define a reference plane. In Fig. 1(a), we examine a 2D case for illustration. The reference plane is tangent to the center of patch (origin point) and perpendicular to the surface normal at . Then the surface distance for patch with respect to normal is defined as a function , where is a point on the reference plane, superscript indicates that the reference plane is perpendicular to , while the subscript indicates the function defines patch . is then the perpendicular distance from to surface with respect to normal . Surface is similarly defined as . Note that because the patches are centered, which is the origin, but their surface normals and are typically different.

The patch distance is then computed as

 d−−→mn=√1|Ωm|∫x∈Ωm(fmm(x)−fmn(x))2dx, (18)

where is the local neighborhood at . is the area of . denotes the distance measured with reference plane perpendicular to .

Note that different reference planes lead to different distance values, so we alternately use and to define the reference plane. Fig. 1(b) illustrates the computation of with reference plane perpendicular to .

 d−−→nm=√1|Ωn|∫x∈Ωn(fnm(x)−fnn(x))2dx, (19)

where functions and define surfaces and , respectively, and the local neighborhood at .

is then given as,

 dmn=√d2−−→mn+d2−−→nm2. (20)

#### Iv-C3 Distance Measure with Discrete Point Observation

Since we only have discrete observations of the points on the patches, we instead measure the sum of the distances between points with the same projection on the reference plane.

First, we compute , where reference plane is perpendicular to the surface normal at . Specifically, patch is composed of points , while patch is composed of . The surface normal is given by,

 minnmk∑i=1((vim)⊤nm)2. (21)

It can be shown via Principal Component Analysis

[25]

that the solution is the normalized eigenvector according to the smallest eigenvalue of the covariance matrix

given by,

 Q=1kk∑i=1vim(vim)⊤. (22)

The same normal estimation method is used in PCL Library [1]. We then project both and to the reference plane, and the projections are and respectively, where the second index in subscript indicates that the reference plane is perpendicular to . The distances between and give the surface distance , and the distances between and give .

Ideally, for any in patch , there exists a point in patch whose projection as shown in Fig. 2(a). However, in real dataset, usually does not have a match in patch with exactly the same projection, as illustrated in Fig. 2(b). In this case, we replace the displacement value of (the green point in Fig. 2(b)), which has the same projection as , with the value of its nearest neighbor in patch in terms of the distance between their projections and . Then is computed as,

 d−−→mn= ⎷1kk∑i=1(fmm(vim)−fmn(vin))2. (23)

Similarly, to compute , we define reference plane with , then compute the projections and and displacements , . For each in patch , we match it to the closest point in patch in terms of projection distance. Then is computed as,

 d−−→nm= ⎷1kk∑i=1(fnm(vim)−fnn(vin))2. (24)

The final distance is given by

 dmn=√d2−−→mn+d2−−→nm2. (25)

#### Iv-C4 Planar Interpolation

The pairwise correspondence is based on nearest neighbor replacement, though more accurate interpolation can be adopted. However, due to the large size of the point cloud, implementing interpolation for all the points can be expensive. Thus we use nearest-neighbor replacement when the distance between point pair is under a threshold . When the distance goes above , we apply the interpolation method described as follows.

As shown in Fig. 3(a), for a point , to find its corresponding interpolation on the other patch, we find the three nearest points (also in terms of projection distance) to form a plane, and the interpolation

is given by its projection along the normal vector

on the plane. It can be easily derived that the distance between and is where is the normal vector for the plane , and .

Note that we choose to use projection on the reference plane (e.g. in Fig. 2(a)) to find the correspondence instead of the point position (e.g. in Fig. 2(a)), because using projection is more robust to noise. For example in Fig. 3(b), the underlying surfaces for two patches are both planar, where the circle points belong to one patch and the star points belong to the other. The correct connections are between points along the vertical lines (in blue). This is accomplished by using projection on the reference plane. On the other hand, if the connection is decided by point position, then the resulting connections are erroneous (in orange) and thus lead to inefficient denoising. Therefore point connection based on projection is closer to the ground truth and more robust to noise.

### Iv-D Graph Construction

Based on the above patch distance measure strategy, the connection between -th and -th patch is implemented as follows. If no interpolation is involved, the points in -th patch are connected with the nearest points in -th patch in terms of their projections on the reference plane decided by surface normal of patch . Also, the points in -th patch are connected with the nearest points in -th patch in terms of their projections on the reference plane decided by surface normal of patch . The edges are undirected and assigned the same weight decided by in (7).

If interpolation is involved, for example in Fig. 3(a), the weight between and is given by,

where is the weight between patch and . Point lies on patch and points , , lie on patch . is distance between and , for , and for . To simplify the implementation, we set a search window to be patches centered at the -nearest patch centers, and evaluate patch distance between these -nearest patches instead of all the patches in the point cloud.

In this way, the point domain graph is constructed, giving the graph Laplacian where is total number of points in all patches.

## V Algorithm Development

### V-a Denoising Algorithm

Given , the optimization is reformulated as:

 minUP⊤xLpPx+P⊤yLpPy+P⊤zLpPz+μ∥V−U∥2F, (27)

where are the vectorized -, - and -coordinates of the points in patches. Let . can be combined as .

is related to denoised 3D samples as follows:

 P=SU−C, (28)

where is a sampling matrix to select points from point cloud to form patches of 3D points each, and is for patch centering. Hence, the objective function can be rewritten as:

 minUtr((SU−C)⊤Lp(SU−C))+μ∥V−U∥2F. (29)

The optimization is non-convex because of ’s dependency on patches in . To solve (29) approximately, we take an alternating approach, where in each iteration, we fix and solve for , then update given , and repeat until convergence.

In each iteration, graph Laplacian is easy to update using the previously discussed graph construction strategy. To optimize for fixed , each of the coordinate is given by,

 (S⊤LpS+μIq)Uq=μVq+S⊤LpCq, (30)

where is the index for coordinates, and

is the identity matrix of the same size as

. We iteratively solve the optimization until the result converges. The proposed algorithm is referred to as Graph Laplacian Regularized point cloud denoising (GLR).

### V-B Graph Spectral Analysis

To impart intuition and demonstrate stability of our computation, in each iteration we can compute the optimal -, - and -coordinates in (29) separately, resulting in the following system of linear equations:

 (S⊤LpS+μI)u∗x=μvx+S⊤Lpcx, (31)

where , and represent the -coordinates (first column) of , and , respectively. In Section III, we assume that union of all patches covers all points in the point cloud , hence we can safely assume that .

Because is a sampling matrix, we can define as a principal sub-matrix222A principal sub-matrix of an original larger matrix is one where the -th row and column of are removed iteratively for different . of . Denote by the eigenvalues of matrix . The solution to (31) can thus be written as:

 u∗x=ΦΣ−1ΦT(μvx+S⊤Lpcx), (32)

where is an eigen-decomposition333Eigen-decomposition is possible because the target matrix is real and symmetric. of matrix ; i.e., contains as columns eigenvectors , and is a diagonal matrix containing eigenvalues on its diagonal. In graph signal processing (GSP) [15], eigenvalues and eigenvectors of a variational operator— in our case—are commonly interpreted as graph frequencies and frequency components. is thus an operator (called

graph Fourier transform

(GFT)) that maps a graph-signal to its GFT coefficients .

Observing that in (32) is a diagonal matrix:

 Σ−1=diag(1/(λμ1+μ),…,1/(λμN+μ)), (33)

we can thus interpret the solution in (32) as follows. The noisy observation (offset by centering vector ) is transformed to the GFT domain via and low-pass filtered per coefficient according to (33)—low-pass because weights for low frequencies are larger than large frequencies , for . The fact that we are performing 3D point cloud denoising via graph spectral low-pass filtering should not be surprising.

### V-C Numerical Stability via Eigen-Analysis

We can also estimate the stability of the system of linear equations in (32) via the following eigen-analysis. During graph construction, an edge weight is computed using (7), which is upper-bounded by . Denote by the maximum degree of a node in the graph, which in general . According to the Gershgorin circle theorem [25], given a matrix , a Gershgorin disc has radius and center at . For a combinatorial graph Laplacian , the maximum Gershgorin disc radius is the maximum node degree multiplied by the maximum edge weight, which is . Further, the diagonal entry for positive edge weights, which equals . Thus all Gershgorin discs for a combinatorial graph Laplacian matrix have left-ends located at . By the Gershgorin circle theorem, all eigenvalues have to locate inside the union of all Gershgorin discs. This means that the maximum eigenvalue for is upper-bounded by twice the radius of the largest possible disc, which is .

Now consider principal sub-matrix of original matrix . By the eigenvalue interlacing theorem, large eigenvalue for is upper-bounded by of . For matrix , the smallest eigenvalue , because: i) shifts all eigenvalues of to the right by , and ii) is PSD due to eigenvalue interlacing theorem and the fact that is PSD. We can thus conclude that the condition number444Assuming -norm is used and the matrix is normal, then the condition number is defined as the ratio . of matrix on the left-hand side of (31) can be upper-bounded as follows:

 C≤2ρmax+μμ. (34)

Hence for sufficiently small , the linear system of equations in (31) has a stable solution, efficiently solved using indirect methods like conjugate gradient (CG).

### V-D Complexity Analysis

The complexity of the algorithm relies on two main procedures: one is the patch-based graph construction, and the other is in solving the system of linear equations.

For graph construction, for the patches, the -nearest patches to be connected can be found in time. Then for -point patch distance measure, each pair takes ; with pairs, the complexity is in total. For the system of linear equations, it can be solved efficiently with CG based methods, with complexity of . Details about the parameter setting are given in Section VI.

## Vi Experimental Results

The proposed scheme GLR is compared with existing works: APSS [7], RIMLS [8], AWLOP [10], and the state-of-the-art moving robust principal components analysis (MRPCA) algorithm [4]. APSS and RIMLS are implemented with MeshLab software [26], AWLOP is implemented with EAR software [10], and MRPCA source code is provided by the author. Models we use are Fandisk, Gargoyle, DC, Daratech, Anchor, and Lordquas provided in [3] and [4].

### Vi-a Implementation Details

To speed up the implementation, we take 50% of the points as the patch centers, with the farthest point sampling [27] to assure spatially uniform selection. The search window size is set to be 16, and the patch size is set to be 30. Parameters that require tuning include in the objective and in deciding the graph edge weight, which will be discussed in Section VI-B and VI-C in detail for particular models.

### Vi-B Demonstration of Iterations

Here we show the results using the model Fandisk in Fig. 4

(a). Gaussian noise with standard deviation of 0.005 is added to the 3D position of the point cloud shown in Fig.

4(b). The surface reconstruction of denoising results in the first two iterations are shown in Fig. 4(c) and (d). The result after iteration 2 converges, so we only show the first two iterations. The parameter which increases along the iterations, and for iteration 1 and for iteration 2.

A comparison with other schemes is shown in Fig. 5. The results of APSS in Fig. 5(b) and RIMLS in Fig. 5(c) are over-smoothed at the folding regions and corners, while AWLOP result in Fig. 5(d) has non-negligible noise. The state-of-the-art MRPCA in Fig. 5(e) has good denoising performance except the artifact at the thin area where the disk shrinks to a sheet. For the proposed scheme in Fig. 5(f), the result is visually better where the folding and corner structures are well preserved without shrinkage at the thin area.

### Vi-C Numerical Evaluation

Apart from visual comparison, we measure the mean square error (MSE) between ground truth and denoising results. Specifically, we first measure the average of the squared Euclidean distances between ground truth points and their closest denoised points, and also between the denoised points and their closest ground truth points, then take the average between the two measures to compute MSE. The MSE results are shown in Table I, where the proposed scheme has the lowest MSE.

More models are included in the evaluation. Gaussian noise with standard deviations of 0.01, 0.02, 0.03, and 0.04 is added to the 3D positions of the point cloud. Numerical results are shown in Table II, III, IV, and V, where the proposed scheme is shown to have the lowest MSE. The models are around 50000 in size. The parameter is set to increase along the iterations, i.e., , where for Gaussian noise with = 0.01, 0.02, 0.03, 0.04, respectively. is set to be for iteration 1, for iteration 2, and for the iterations afterwards.

Illustrations of the Daratech model are shown in Fig. 6 and 7 from different perspectives. AWLOP denoising results are not showing competitive results, which is also seen in Table IV, thus not included in the figures. As shown in Fig. 6(c-e) and Fig. 7(c-e), the existing schemes are under-smoothed in plane area, while in thin area the surface shrinks and the details are lost. However, for the proposed method, the details are well preserved without under-smoothing.

In addition, illustrations of the Lordquas model are shown in Fig. 8. For the existing schemes, the cigarette is broken and the edges (for example, the ears) are over-smoothed, while the proposed method is visually better, and structural details are preserved without over-smoothing.

## Vii Conclusion

In this paper, we propose a graph Laplacian regularization based 3D point cloud denoising algorithm. To utilize the self-similarity among surface patches, we adopt the low-dimensional manifold prior, and collaboratively denoise the patches by minimizing the manifold dimension. For efficient implementation, we approximate the manifold dimension with a graph Laplacian regularizer, and construct the patch graph with a new measure for the discrete patch distance. The proposed scheme is shown to have graph spectral low-pass filtering interpretation and numerical stability in solving the linear equation system. Experimental results suggest that our proposal can outperform existing schemes with better structural detail preservation.

## References

• [1] R. B. Rusu and S. Cousins, “3D is here: Point cloud library (PCL),” in Robotics and Automation (ICRA), 2011 IEEE International Conference on.   IEEE, 2011, pp. 1–4.
• [2]

M. Ji, J. Gall, H. Zheng, Y. Liu, and L. Fang, “Surfacenet: An end-to-end 3D neural network for multiview stereopsis,”

2017 IEEE International Conference on Computer Vision (ICCV)

, pp. 2326–2334, 2017.
• [3] G. Rosman, A. Dubrovina, and R. Kimmel, “Patch-collaborative spectral point-cloud denoising,” Computer Graphics Forum, vol. 32, no. 8, pp. 1–12, 2013.
• [4] E. Mattei and A. Castrodad, “Point cloud denoising via moving rpca,” Computer Graphics Forum, pp. 1–15, 2016.
• [5] Y. Sun, S. Schaefer, and W. Wang, “Denoising point sets via l minimization,” Computer Aided Geometric Design, vol. 35, pp. 2–15, 2015.
• [6] Y. Zheng, G. Li, S. Wu, Y. Liu, and Y. Gao, “Guided point cloud denoising via sharp feature skeletons,” The Visual Computer, pp. 1–11, 2017.
• [7] G. Guennebaud and M. Gross, “Algebraic point set surfaces,” ACM Transactions on Graphics (TOG), vol. 26, no. 3, p. 23, 2007.
• [8] A. C. Öztireli, G. Guennebaud, and M. Gross, “Feature preserving point set surfaces based on non-linear kernel regression,” Computer Graphics Forum, vol. 28, no. 2, pp. 493–501, 2009.
• [9] Y. Lipman, D. Cohen-Or, D. Levin, and H. Tal-Ezer, “Parameterization-free projection for geometry reconstruction,” ACM Transactions on Graphics (TOG), vol. 26, no. 3, p. 22, 2007.
• [10] H. Huang, S. Wu, M. Gong, D. Cohen-Or, U. Ascher, and H. R. Zhang, “Edge-aware point set resampling,” ACM Transactions on Graphics (TOG), vol. 32, no. 1, p. 9, 2013.
• [11] X.-F. Han, J. S. Jin, M.-J. Wang, W. Jiang, L. Gao, and L. Xiao, “A review of algorithms for filtering the 3D point cloud,” Signal Processing: Image Communication, 2017.
• [12] A. Buades, B. Coll, and J.-M. Morel, “A non-local algorithm for image denoising,” in

Computer Vision and Pattern Recognition (CVPR), 2005 IEEE Computer Society Conference on

, vol. 2.   IEEE, 2005, pp. 60–65.
• [13] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by sparse 3-D transform-domain collaborative filtering,” IEEE Transactions on Image Processing, vol. 16, no. 8, pp. 2080–2095, 2007.
• [14] S. Osher, Z. Shi, and W. Zhu, “Low dimensional manifold model for image processing,” SIAM Journal on Imaging Sciences, vol. 10, no. 4, pp. 1669–1690, 2017.
• [15]

D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,”

IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 83–98, May 2013.
• [16] M. Alexa, J. Behr, D. Cohen-Or, S. Fleishman, D. Levin, and C. T. Silva, “Computing and rendering point set surfaces,” IEEE Transactions on Visualization and Computer Graphics, vol. 9, no. 1, pp. 3–15, 2003.
• [17] G. Guennebaud, M. Germann, and M. Gross, “Dynamic sampling and rendering of algebraic point set surfaces,” Computer Graphics Forum, vol. 27, no. 2, pp. 653–662, 2008.
• [18] R. B. Rusu, N. Blodow, Z. Marton, A. Soos, and M. Beetz, “Towards 3D object maps for autonomous household robots,” in Intelligent Robots and Systems (IROS), 2007 IEEE/RSJ International Conference on.   IEEE, 2007, pp. 3191–3198.
• [19] H. Huang, D. Li, H. Zhang, U. Ascher, and D. Cohen-Or, “Consolidation of unorganized point clouds for surface reconstruction,” ACM Transactions on Graphics (TOG), vol. 28, no. 5, p. 176, 2009.
• [20] H. Avron, A. Sharf, C. Greif, and D. Cohen-Or, “l-sparse reconstruction of sharp point set surfaces,” ACM Transactions on Graphics (TOG), vol. 29, no. 5, p. 135, 2010.
• [21] T. Guillemot, A. Almansa, and T. Boubekeur, “Non local point set surfaces,” in 3D Imaging, Modeling, Processing, Visualization and Transmission (3DIMPVT), 2012 Second International Conference on.   IEEE, 2012, pp. 324–331.
• [22] J. Digne, “Similarity based filtering of point clouds,” in Computer Vision and Pattern Recognition Workshops (CVPRW), 2012 IEEE Computer Society Conference on.   IEEE, 2012, pp. 73–79.
• [23] M. Hein, “Uniform convergence of adaptive graph-based regularization,” Lecture Notes in Computer Science, vol. 4005, p. 50, 2006.
• [24] J. Pang and G. Cheung, “Graph Laplacian regularization for image denoising: Analysis in the continuous domain,” IEEE Transactions on Image Processing, vol. 26, no. 4, pp. 1770–1785, 2017.
• [25] R. Horn and C. Johnson, Matrix Analysis.   Cambridge University Press, 2012.
• [26] P. Cignoni, M. Callieri, M. Corsini, M. Dellepiane, F. Ganovelli, and G. Ranzuglia, “Meshlab: An open-source mesh processing tool.” in Eurographics Italian Chapter Conference, vol. 2008, 2008, pp. 129–136.
• [27] Y. Eldar, M. Lindenbaum, M. Porat, and Y. Y. Zeevi, “The farthest point strategy for progressive image sampling,” IEEE Transactions on Image Processing, vol. 6, no. 9, pp. 1305–1315, 1997.