3D Point Cloud Super-Resolution via Graph Total Variation on Surface Normals

TL;DR

This paper introduces a graph-based super-resolution method for 3D point clouds that enhances resolution by optimizing surface normals and point positions using total variation and ADMM, outperforming existing methods.

Contribution

The paper proposes a novel local algorithm leveraging graph total variation and bipartite graph approximation for 3D point cloud super-resolution, with an effective optimization framework.

Findings

Outperforms competing super-resolution methods objectively.

Produces higher quality 3D surface reconstructions.

Effectively preserves original point coordinates during enhancement.

Abstract

Point cloud is a collection of 3D coordinates that are discrete geometric samples of an object's 2D surfaces. Using a low-cost 3D scanner to acquire data means that point clouds are often in lower resolution than desired for rendering on high-resolution displays. Building on recent advances in graph signal processing, we design a local algorithm for 3D point cloud super-resolution (SR). First, we initialize new points at centroids of local triangles formed using the low-resolution point cloud, and connect all points using a k-nearestneighbor graph. Then, to establish a linear relationship between surface normals and 3D point coordinates, we perform bipartite graph approximation to divide all nodes into two disjoint sets, which are optimized alternately until convergence. For each node set, to promote piecewise smooth (PWS) 2D surfaces, we design a graph total variation (GTV) objective for nearby surface normals, under the constraint that coordinates of the original points are preserved. We pursue an augmented Lagrangian approach to tackle the optimization, and solve the unconstrained equivalent using the alternating method of multipliers (ADMM). Extensive experiments show that our proposed point cloud SR algorithm outperforms competing schemes objectively and subjectively for a large variety of point clouds.

Tables

Table 1. Table 1 : C2C ( × 10 − 1 absent superscript 10 1 \times 10^{-1} ) of different models

Model	Low	APSS	RIMLS	Initial	Prop.
Bunny	1.47	1.31	1.22	1.27	1.14
Dragon	1.52	1.43	1.34	1.45	1.25
Armadillo	1.49	1.38	1.30	1.37	1.21
Buddha	1.55	1.42	1.32	1.46	1.27
Asian Drag.	1.56	1.41	1.39	1.44	1.31
Lucy	1.58	1.46	1.40	1.51	1.28

Table 2. Table 2 : C2P ( × 10 − 2 absent superscript 10 2 \times 10^{-2} ) of different models

Model	Low	APSS	RIMLS	Initial	Prop.
Bunny	3.92	2.40	2.28	2.49	1.87
Dragon	2.81	1.71	1.58	1.77	1.26
Armadillo	2.55	1.42	1.30	1.37	1.07
Buddha	3.11	1.93	1.81	1.88	1.38
Asian Drag.	2.76	1.62	1.47	1.57	1.31
Lucy	3.53	2.36	2.11	2.44	1.68

Equations

$$\mathbf{Cp}=\mathbf{q}.$$

$$w_{i,j}=\exp \left\{-\frac{||{\mathbf{p}}_i-{\mathbf{p}}_j|{|}_2^2}{{\sigma }_p^2}\right\}{\cos }^2{\theta }_{i,j},\text{vspace}-5pt$$

$$||\mathbf{n}|{|}_{\text{GTV}}=\sum _{i,j\in \mathcal{E}}{w}_{i,j}||{\mathbf{n}}_i-{\mathbf{n}}_j|{|}_1,$$

$${\mathbf{n}}_i={\mathbf{A}}_i{\mathbf{p}}_i+{\mathbf{b}}_i,$$

$$\mathbf{m}=\mathbf{B}\mathbf{p}+\mathbf{v},$$

$$\displaystyle\underbrace{\left[\begin{array}[]{cc}{\mathbf{I}}&-{\mathbf{B}}\\ 0&{\mathbf{C}}\end{array}\right]}_{{\mathbf{H}}}\underbrace{\left[\begin{array}[]{c}{\mathbf{m}}\\ {\mathbf{p}}\end{array}\right]}_{{\mathbf{s}}}=\underbrace{\left[\begin{array}[]{c}{\mathbf{v}}\\ {\mathbf{q}}\end{array}\right]}_{\mathbf{d}}$$

$$\begin{array}{rl}& \underset{\mathbf{s}}{\min }\sum _{i,j}{w}_{i,j}\Vert {\mathbf{m}}_{i,j}{\Vert }_1\end{array}$$

$$\underset{\mathbf{s}}{\min }\sum _{i,j}{w}_{i,j}\Vert {\mathbf{m}}_{i,j}{\Vert }_1+{\mathbf{y}}^{\top }(\mathbf{Hs}-\mathbf{d})+\frac{\rho }{2}{\Vert \mathbf{Hs}-\mathbf{d}\Vert }_2^2,$$

$${\Vert \mathbf{Hs}-\mathbf{d}\Vert }_2^2={\Vert \mathbf{Bp}+\mathbf{v}-\mathbf{m}\Vert }_2^2+{\Vert \mathbf{Cp}-\mathbf{q}\Vert }_2^2.$$

$${\mathbf{y}}^{\top }(\mathbf{Hs}-\mathbf{d})={\mathbf{y}}_1^{\top }(\mathbf{Bp}+\mathbf{v}-\mathbf{m})+{\mathbf{y}}_2^{\top }(\mathbf{Cp}-\mathbf{q}),$$

$$\begin{array}{r}\underset{\mathbf{m},\mathbf{p}}{\min }\sum _{i,j}{w}_{i,j}\Vert {\mathbf{m}}_{i,j}{\Vert }_1+{\mathbf{y}}_1^{\top }(\mathbf{Bp}+\mathbf{v}-\mathbf{m})+\\ {\mathbf{y}}_2^{\top }(\mathbf{Cp}-\mathbf{q})+\frac{\rho }{2}{\Vert \mathbf{Bp}+\mathbf{v}-\mathbf{m}\Vert }_2^2+\frac{\rho }{2}{\Vert \mathbf{Cp}-\mathbf{q}\Vert }_2^2.\end{array}$$

$$\begin{array}{r}\rho ({\mathbf{B}}^{\top }\mathbf{B}+{\mathbf{C}}^{\top }\mathbf{C}){\mathbf{p}}^{k+1}={\mathbf{C}}^{\top }(\rho \mathbf{q}-{\mathbf{y}}_2)+\\ {\mathbf{B}}^{\top }(\rho {\mathbf{m}}^k-\rho \mathbf{v}-{\mathbf{y}}_1)\end{array}$$

$$\begin{array}{r}\underset{\mathbf{m}}{\min }{\mathbf{y}}_1^{k^{\top }}(\mathbf{B}{\mathbf{p}}^{\mathbf{k}+1}+\mathbf{v}-\mathbf{m})+\frac{\rho }{2}{\Vert \mathbf{B}{\mathbf{p}}^{\mathbf{k}+1}+\mathbf{v}-\mathbf{m}\Vert }_2^2+\\ \sum _{i,j}{w}_{i,j}{\Vert {\mathbf{m}}_{i,j}\Vert }_1\end{array}$$

$${\Delta }_{\mathbf{m}}({\mathbf{p}}^{k+1},\mathbf{m},{\mathbf{y}}_1^k)=-\rho ({\mathbf{Bp}}^{k+1}+\mathbf{v}-\mathbf{m})-{\mathbf{y}}_1^k.$$

$${\text{prox}}_{g,t}(\mathbf{m})=\text{arg}\underset{\theta }{\min }\left\{g(\theta )+\frac{1}{t}||\theta -\mathbf{m}|{|}_2^2\right\}$$

$${\text{prox}}_{g,t}(m_{i,j,r})=\{\begin{array}{ll}{m}_{i,j,r}-t{w}_{i,j}& \text{if}{m}_{i,j,r}>t{w}_{i,j}\\ 0& \text{if}|m_{i,j,r}|\le t{w}_{i,j}\\ m_{i,j,r}+t{w}_{i,j}& \text{if}{m}_{i,j,r}<-t{w}_{i,j},\end{array}$$

$${\mathbf{m}}^{k+1}={\text{prox}}_{g,t}({\mathbf{m}}^k-t{\Delta }_{\mathbf{m}}({\mathbf{p}}^{k+1},{\mathbf{m}}^k,{\mathbf{y}}_1^k)).$$

$${\mathbf{y}}^{k+1}={\mathbf{y}}^k+\rho ({\mathbf{As}}^{k+1}-b).$$

Full text

Index Terms—  graph signal processing, point cloud super-resolution, graph total variation, convex optimization

1 Introduction

Point cloud, acquired directly by off-the-shelf 3D scanners like Microsoft Kinect or estimated indirectly via stereo-matching algorithms [1], is a recently popular 3D visual signal representation for free viewpoint image rendering, and is investigated in industrial standards like MPEG111https://mpeg.chiariglione.org/standards/mpeg-i/point-cloud-compression. Unlike 3D meshes, a point cloud is an unstructured list of 3D coordinates, and low-level processing tasks like compression [2, 3, 4] and denoising [5, 6] are challenging. In this paper, we focus on the little-studied point cloud super-resolution (SR) problem: how to increase the density of the 3D points while preserving the geometry of the underlying 3D object structure?

Compared to the formidable body of literature for image SR [7, 8, 9, 10], point cloud SR has received relatively little attention. In [11], resolution of a point cloud is increased by combining similar low-res point cloud patches via simple local shifting and aggregation. This method relies on the existence of a large number of similar point cloud regions, which may not be true in practice. In [12, 13, 14], surface interpolation is accomplished using the Moving Least Squares (MLS) method. However, it is observed that these methods often over-smooth due to MLS interpolation. As one solution, [15] proposed an edge-aware point upsampling method by first sampling away from object boundaries, and then progress towards object edges and corners. However, [15] is ad-hoc in methodology and still suffers from over-smoothing, albeit to a lesser extent. Recently a deep learning based method [16] learns multilevel features per point and then upsamples the set of points via a multi-branch convolution unit implicit in the feature space. However, a large training dataset with similar geometric characteristics as the target point cloud is required, which may not be practical.

In this paper, leveraging on recent progress in graph signal processing (GSP) [17, 18] and our previous work on point cloud denoising [5, 6], we propose a fast local point cloud SR algorithm based on graph total variation (GTV). First, we initialize new points at centroids of local triangles formed via Delaunay triangulation222https //en.wikipedia.org/wiki/Delaunay_triangulation using the low-res point cloud. We connect all points using a $k$-nearest-neighbor graph. Then, to establish a tractable linear relationship between surface normals and 3D point coordinates (similarly done in [5]), we perform bipartite graph approximation [19] to divide all nodes into two disjoint sets, which are optimized alternately until convergence. For each node set, in order to promote piecewise smoothness (PWS) in reconstructed 2D surfaces, we design a GTV objective for nearby surface normals, under the constraint that 3D coordinates of the original points are preserved. We pursue an augmented Lagrangian approach [20] and solve the unconstrained equivalent using the alternating method of multipliers (ADMM) and proximal gradient decent [20, 21]. We demonstrate in extensive experiments that our proposed point cloud SR algorithm outperforms competing schemes for a large variety of 3D point clouds objectively and subjectively.

2 Preliminaries

3D Point Cloud: We define a point cloud as a set of (roughly uniform) discrete samples of 3D coordinates on an object’s 2D surface in 3D space. Denote by ${\mathbf{p}}=\left[{\mathbf{p}}_{1}^{\top}\;\ldots\;\mathbf{p}_{N}^{\top}\right]^{\top}\in\mathbb{R}^{3N}$ the position vector for a full-resolution point cloud, where ${\mathbf{p}}_{i}\in\mathbb{R}^{3}$ is the 3D coordinate of a point $i$, and $N$ is the number of points in the point cloud. Similarly, denote by ${\mathbf{q}}=\left[\mathbf{q}_{1}^{\top}\;\ldots\;{\mathbf{q}}_{M}^{\top}\right]^{\top}\in\mathbb{R}^{3M}$ the position vector for a low-res point cloud to the above one, where ${\mathbf{q}}_{i}\in\mathbb{R}^{3}$, and $M$, $M<N$, is the number of points in the low-res point cloud. We assume that the observed low-res ${\mathbf{q}}$ is noiseless. Full-resolution ${\mathbf{p}}$ and low-res ${\mathbf{q}}$ are related through the sampling matrix $\mathbf{C}\in\{0,1\}^{3M\times 3N}$; i.e.,

[TABLE]

Graph Construction from a 3D Point Cloud: We construct a $k$-nearest-neighbor ($k$-NN) graph ${\mathcal{G}}$ to connect 3D points in a point cloud [22], so that each point can be filtered with neighboring points under a graph-structured data kernel [17, 18]. Specifically, consider an undirected graph $\mathcal{G}=({\mathcal{V}},{\mathcal{E}})$ with node set $\mathcal{V}$ and edge set $\mathcal{E}$. ${\mathcal{G}}$ is specified by $(i,j,w_{i,j})$, where $i,j\in{\mathcal{V}}$, $(i,j)\in{\mathcal{E}}$, and $w_{i,j}\in\mathbb{R}^{+}$ is the weight of an edge that connects nodes $i$ and $j$. Each 3D point $i$ is represented by a node $i\in{\mathcal{V}}$ and is connected via edges to its $k$ nearest neighbors $j$ in Euclidean distance, with weights $w_{i,j}$ that reflect inter-node similarities [6], i.e.,

[TABLE]

where $\theta_{i,j}$ is the angle between surface normals $\mathbf{n}_{i}$ and $\mathbf{n}_{j}$ at respective nodes $i$ and $j$, and $\sigma_{p}$ is a parameter. In words, (2) states that weight $w_{i,j}$ is large (close to 1) if points $i$ and $j$ are physically close and associated normal vectors are similar, and small (close to 0) otherwise. The weight definition in (2) is similar to bilateral filter weights defined in [23] with domain and range filters.

Surface Normals: A surface normal ${\mathbf{n}}_{i}\in\mathbb{R}^{3}$ of node $i$ is a vector that is perpendicular to the tangent plane at point $i$. Typically coordinates of the $k$ nearest neighbors of $i$ are used to compute ${\mathbf{n}}_{i}$ [24, 25, 26, 27, 28]. The most popular method is to fit a local plane to point $i$ and its $k$ nearest neighbors, and take the perpendicular vector to that plane (see e.g. [24, 25, 26]). An alternative is to calculate the normal vector as the weighted average of the normal vectors of the triangles formed by point $i$ and different pairs of its neighbors (see e.g. [27, 28]).

3 Proposed Algorithm

3.1 Algorithm Overview

We begin by adding a set of new interior points to the low-res point cloud $\mathbf{q}$ to populate the target full-res point cloud $\mathbf{p}$. Specifically, we construct a triangular mesh using Delaunay triangulation using points in ${\mathbf{q}}$, and then insert new points at the centroids of those triangles. Given an interpolated point cloud with appropriate density, we assume a graph signal smoothness prior to fine-tune the newly added 3D coordinates described as follows.

Specifically, we assume that the underlying 2D surface of a 3D object is piecewise smooth (PWS): surface normals ${\mathbf{n}}_{i}$ and ${\mathbf{n}}_{j}$ of two neighboring 3D points $i$ and $j$ should differ minimally over the 2D surface, except at the object boundaries (e.g., sides of a rectangular box), where the normal difference is large. Hence, to promote small overall deviation with sparse but large changes, we design GTV to promote PWS333$l_{1}$ norm can better promote PWS than $l_{2}$ norm [29]. of surface normals as follows (similarly done in [5]):

[TABLE]

We then formulate our point cloud SR problem as a minimization of the defined GTV while enforcing the locations of observed points $\mathbf{q}$. Here, ${\mathbf{p}}$ is the optimization variable, and ${\mathbf{n}}_{i}$’s are functions of $\mathbf{p}$. Unfortunately, using state-of-art surface normal estimation methods, each ${\mathbf{n}}_{i}$ is a nonlinear function of ${\mathbf{p}}_{i}$ and its neighbors. Hence, it is difficult to formulate a clean convex optimization using the present form of GTV in (3).

To overcome this problem, following our previous work on point cloud denoising [5], we first partition 3D points of the point cloud into two classes (say red and blue). When computing the normal for a red point, we employ as reference only neighboring blue points, and vice versa. Towards this goal, a $k$-neighborhood graph ${\mathcal{G}}$ is initially constructed for all 3D points $\mathcal{V}$ in a point cloud. We next divide 3D points ${\mathcal{V}}$ in a point cloud into two disjoint sets—red nodes ${\mathcal{V}}_{1}$ and blue nodes ${\mathcal{V}}_{2}$ where ${\mathcal{V}}_{1}\cup{\mathcal{V}}_{2}={\mathcal{V}}$ and ${\mathcal{V}}_{1}\cap{\mathcal{V}}_{2}=\emptyset$—via bipartite graph approximation [19]; i.e., finding a bipartite graph that best approximates the original non-bipartite graph in terms of Kullback-Leibler (KL) divergence444https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence. An example for the bipartite graph approximation is shown in Fig. 1. Doing so means that each red node $i\in{\mathcal{V}}_{1}$ has in its neighborhood two or more blue nodes, using which node $i$’s surface normal ${\mathbf{n}}_{i}$ can be written as a linear function of its 3D coordinate ${\mathbf{p}}_{i}$, i.e.,

[TABLE]

where $\mathbf{A}_{i}\in\mathbb{R}^{3\times 3}$ and $\mathbf{b}_{i}\in\mathbb{R}^{3}$ are computed from coordinates of the neighbouring blue nodes. See [5] for further details.

3.2 Optimization Objective

We next construct a new $k$-NN graph $\mathcal{G}_{1}~{}=~{}(\mathcal{V}_{1},\mathcal{E}_{1})$ connecting only neighboring red nodes, where $\mathcal{E}_{1}$ is the set of edges of the new graph. Edge weights $w_{i,j}$ are computed using (2), where $\theta_{i,j}$ is computed using initial surface normal estimates. Denote by ${\mathbf{m}}_{i,j}$ the difference between normals $\mathbf{n}_{i}$ and $\mathbf{n}_{j}$ of two connected nodes $i$ and $j$ in $\mathcal{G}_{1}$, i.e., ${\mathbf{m}}_{i,j}={\mathbf{n}}_{i}-{\mathbf{n}}_{j}$, where $(i,j)\in\mathcal{E}_{1}$. Given that each normal ${\mathbf{n}}_{i}$ is a linear function of node $i$’s 3D coordinates ${\mathbf{p}}_{i}$ (4), we can write ${\mathbf{m}}=\left[\{{\mathbf{m}}_{i,j}\}\right]^{\top}$ also as a linear function of ${\mathbf{p}}=\left[\{{\mathbf{p}}_{i}\}\right]^{\top}$:

[TABLE]

Given the relationship between low-res ${\mathbf{q}}$ and full-res ${\mathbf{p}}$ through sampling matrix (1), we can summarize the two constraints compactly as a single equality constraint as follows:

[TABLE]

where ${\mathbf{s}}=[{\mathbf{m}}\;\;{\mathbf{p}}]^{\top}$ is the combined optimization variable.

Our objective is to minimize the GTV of surface normals:

[TABLE]

subject to linear constraint (12). In our optimization framework, we first solve (13) for red nodes while the positions of blue nodes are fixed. Then, using the newly solved red nodes’ positions, we initialize the normals for blue nodes and construct another $k$-NN graph for blue nodes. Then (13) is solved for blue nodes while the positions of red nodes fixed. The two node sets are alternately optimized until convergence. We next discuss how to solve (13) for one node set.

3.3 Optimization Algorithm

To solve the problem in (13), we pursue an augmented Lagrangian approach and rewrite the constrained problem (13) into the following unconstrained form:

[TABLE]

where ${\mathbf{y}}$ is a dual variable and $\rho$ is a Lagrange multiplier. The unconstrained version (14) can now be solved using ADMM and proximal gradient decent. According to (12), we can rewrite $\left\lVert\mathbf{Hs-d}\right\rVert_{2}^{2}$ as

[TABLE]

Similarly, according to (12), $\mathbf{y}^{\top}(\mathbf{Hs-d})$ can be rewritten as

[TABLE]

where $\mathbf{y}^{\top}=[\mathbf{y}_{1}^{\top}\hskip 8.0pt\mathbf{y}_{2}^{\top}]$. We can thus rewrite (14) as

[TABLE]

As typically done in ADMM approaches, we solve (17) by alternately minimizing $\mathbf{p}$ and $\mathbf{m}$ and updating $\mathbf{y}$ one at a time in turn until convergence.

$\mathbf{p}$** minimization**: To minimize $\mathbf{p}$ having $\mathbf{m}^{k}$ and $\mathbf{u}^{k}$ fixed, we take the derivative of (17) with respect to $\mathbf{p}$, set it to [math] and solve for the closed form solution $\mathbf{p}^{k+1}$:

[TABLE]

One can show that matrix $(\mathbf{B}^{\top}\mathbf{B}~{}+~{}\mathbf{C}^{\top}\mathbf{C})$ is positive definite (PD) but can be numerically unstable with a large condition number (the ratio of the largest eigenvalue to the smallest). To obtain a stable solution nonetheless for (18), we use the iterative refinement method proposed in [21].

$\mathbf{m}$** minimization**: Keeping $\mathbf{p}^{k+1}$ and $\mathbf{y}^{k}$ fixed, the minimization of $\mathbf{m}$ becomes:

[TABLE]

where the first two terms are convex and differentiable, and the third term is convex but non-differentiable. We can thus use proximal gradient[21] to solve (19). The first two terms have gradient $\Delta_{\mathbf{m}}$:

[TABLE]

We can now define a proximal mapping $\text{prox}_{g,t}(\mathbf{m})$ for a convex, non-differentiable function $g()$ with step size $t$ as:

[TABLE]

We know that for our weighted $l_{1}$-norm in (19), the proximal mapping is just a soft thresholding function:

[TABLE]

where $m_{i,j,r}$ is the $r$-th entry of $\mathbf{m}_{i,j}$. We can now update $\mathbf{m}^{k+1}$ as:

[TABLE]

We compute (23) iteratively until convergence.

$\mathbf{y}$-update: Finally, we can update $\mathbf{y}^{k+1}$ simply:

[TABLE]

$\mathbf{p}$, $\mathbf{m}$ and $\mathbf{y}$ are iteratively optimized in turn using (18), (23) and (24) until convergence.

4 Experimental Results

The computed SR point clouds are compared against point cloud models obtained using competing schemes APSS [13] and RIMLS [14], as well as initial point cloud (after adding new points to the centroids of triangles). Point cloud models we use are Bunny, Dragon, Armadillo, Happy Buddha, Asian Dragon, and Lucy, provided in [30]. All point cloud models are first rescaled, so that each tightly fits inside a bounding box with the same diagonal. Both quantitative and visual comparisons are presented.

For numerical comparisons, we measure the point-to-point (denoted as C2C) error and point to plane (denoted as C2P) error [31] between ground truth and SR results. In C2C error, we first measure the average of the squared Euclidean distances between ground truth points and their closest points in the SR cloud, and also that between the points of the SR cloud and their closest ground truth points. Then the larger value among these two measures is computed as C2C error. In C2P error, we first measure the average of the squared Euclidean distances between ground truth points and tangent planes at their closest points in the SR cloud, and also that between the points of the SR cloud and tangent planes at their closest ground truth points. Again, the larger value between these two measures is computed as the C2P error.

First, the ground truth point cloud are downsampled to 30% of the original points using Poisson Disk Sampling method that was implemented in MeshLab software tool [32]. Then, downsampled point clouds are upsampled to get the same number of points of their ground truths by using the proposed method and existing methods. Numerical results are shown in Table 1 (with C2C error) and Table 2 (with C2P error). We observe that our proposed method has the lowest C2C and C2P errors. In each experiment, the selected parameters are $\rho=5$, $t=0.1$, and $k=8$ when constructing $k-$NN graphs.

Apart from the numerical comparison, visual results for Bunny and Armadillo models are shown in Fig. 2 and 3, respectively. We see that APSS and RIMLS schemes generate overly smooth models compared to our proposed method. Further, existing methods result in distorted surfaces with some details lost. This can be clearly seen at the left ear of bunny model and and the fingers (legs and hands) of armadillo model. On the other hand, using our proposed method, the details are well preserved without over-smoothing.

5 Conclusion

We pursue a graph-based approach to tackle the point cloud super-resolution (SR) problem, where the density of the point cloud is increased while preserving piecewise smoothness (PWS) of the intended object’s 2D surface. Specifically, we first initialize new points at centroids of local triangles, and construct a $k$-nearest-neighbor graph to connect all 3D points. We then divide the nodes into two disjoint sets (red and blue) via bipartite graph approximation, which are subsequently optimized in an alternating manner. Surface normal for each red node can now be defined as a linear function of its 3D coordinate using neighboring blue nodes as reference. We can thus formulate a convex optimization problem, with a graph total variation (GTV) on surface normals as objective and a linear equality constraint to enforce original point coordinates. Experimental results show that our proposed method outperforms competing schemes objectively and subjectively for a variety of 3D point clouds.