Fast Resampling of 3D Point Clouds via Graphs

TL;DR

This paper introduces a graph-based randomized resampling method for 3D point clouds that efficiently preserves application-specific features while reducing data size, ensuring invariance to transformations.

Contribution

It proposes a novel optimal resampling distribution based on graph filters that is invariant to shift, rotation, and scale, applicable to various 3D processing tasks.

Findings

Resampling preserves features effectively across applications.

Method outperforms existing techniques in efficiency and accuracy.

Validated on visualization, registration, and shape modeling tasks.

Abstract

To reduce cost in storing, processing and visualizing a large-scale point cloud, we consider a randomized resampling strategy to select a representative subset of points while preserving application-dependent features. The proposed strategy is based on graphs, which can represent underlying surfaces and lend themselves well to efficient computation. We use a general feature-extraction operator to represent application-dependent features and propose a general reconstruction error to evaluate the quality of resampling. We obtain a general form of optimal resampling distribution by minimizing the reconstruction error. The proposed optimal resampling distribution is guaranteed to be shift, rotation and scale-invariant in the 3D space. We next specify the feature-extraction operator to be a graph filter and study specific resampling strategies based on all-pass, low-pass, high-pass graph filtering and graph filter banks. We finally apply the proposed methods to three applications: large-scale visualization, accurate registration and robust shape modeling. The empirical performance validates the effectiveness and efficiency of the proposed resampling methods.

Figures

Tables

Table 1. TABLE I: Proposed resampling strategy with low-pass graph filtering provides a robust shape modeling for a fitness ball. The first column is the ground truth. The relative error is shown in the parentheses. Best results are marked in bold.

	Original ball	Noisy ball	Uniform resampling	Denoised ball	Low-pass graph filtering based resampling
	(Figure 8 (a) )	(Figure 9 (a) )	(Figure 9 (b) )	(Figure 9 (c) )	(Figure 9 (d) )
Radius	$0.3182$	$0.3478\left(9.3023\%\right)$	$0.3520\left(10.6223\%\right)$	$0.3143\left(1.2256\%\right)$	$0.3199\left(0.5343\%\right)$
Center-$x$	$0.0833$	$0.0903\left(8.4034\%\right)$	$0.0975\left(17.0468\%\right)$	$0.0799\left(4.0816\%\right)$	$0.0849\left(1.9208\%\right)$
Center-$y$	$0.1903$	$0.2136\left(12.2438\%\right)$	$0.1794\left(5.7278\%\right)$	$0.1866\left(1.9443\%\right)$	$0.1783\left(6.3058\%\right)$
Center-$z$	$1.1725$	$1.3803\left(17.7228\%\right)$	$1.1530\left(1.6631\%\right)$	$1.1618\left(0.9126\%\right)$	$1.1613\left(0.9552\%\right)$

Table 2. TABLE II: Proposed high-pass graph filtering based resampling strategy provides an accurate registration for a sofa. Best results are marked in bold. High-pass graph filtering based resampling chooses key points and provides the best registration performance.

	RMSE	${\mathrm{Error}}_{\mathrm{shift}}$	${\mathrm{Error}}_{\mathrm{rotation}}$
All points	$4.22$	$8.76$	$2.30\times {10}^{-3}$
Uniform resampling	$4.27$	$9.38$	$3.76\times {10}^{-3}$
High-pass graph filtering	$1.49$	$0.01$	$4.29\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$
based resampling (16)

Equations

$$X=\left[\begin{array}{cccc}{\mathbf{s}}_1& {\mathbf{s}}_2& \dots & {\mathbf{s}}_K\end{array}\right]=\left[\begin{array}{c}{\mathbf{x}}_1^T\\ {\mathbf{x}}_2^T\\ ⋮\\ {\mathbf{x}}_N^T\end{array}\right]\in {\mathbb{R}}^{N\times K},$$

$$\Psi_{i,j}=\left\{\begin{array}[]{rl}1,&j=\mathcal{M}_{i};\\ 0,&\mbox{otherwise}.\end{array}\right.$$

$$\operatorname{W}_{i,j}=\left\{\begin{array}[]{rl}e^{-\frac{\left\|\mathbf{x}^{\rm(c)}_{i}-\mathbf{x}^{\rm(c)}_{j}\right\|_{2}^{2}}{\sigma^{2}}},&\left\|\mathbf{x}^{\rm(c)}_{i}-\mathbf{x}^{\rm(c)}_{j}\right\|_{2}\leq\tau;\\ 0,&\mbox{otherwise},\end{array}\right.$$

$$h(A)=\sum _{\ell =0}^{L-1}{h}_{\ell }{A}^{\ell }=h_{0}I+h_{1}A+\dots +h_{L-1}{A}^{L-1},$$

$$A=V\Lambda V^{-1},$$

$$\hat{\mathbf{s}}=V^{-1}\mathbf{s}.$$

$$D_{f(X)}(\Psi )={\Vert S{\Psi }^T\Psi f(X)-f(X)\Vert }_2^2,$$

$${\mathbb{E}}_{\Psi \sim \pi }\left({\Psi }^T\Psi f(X)\right)$$

$${\mathbb{E}}_{\Psi \sim \pi }\left(S{\Psi }^T\Psi f(X)\right)$$

$${\mathbb{E}}_{\Psi \sim \pi }{D}_{f(X)}(\Psi )=\mathrm{Tr}\left(f(X)Qf(X{)}^T\right),$$

$$D_f(\Psi )$$

$${\mathbb{E}}_{\Psi \sim \pi }\left(D_f(\Psi )\right)=c^2\mathrm{Tr}\left(FQ{F}^T\right)+\mathrm{Tr}\left(F{X}_{\mathrm{o}}Q(F{X}_{\mathrm{o}}{)}^T\right),$$

$${\pi }_i^{\ast }\propto {\Vert f_i(X)\Vert }_2,$$

$${\pi }_i^{\ast }\propto \sqrt{c^2{\Vert F_i\Vert }_2^2+{\Vert {\left(F{X}_{\mathrm{o}}\right)}_i\Vert }_2^2},$$

$$f(X)=h(A)X=\sum _{\ell =0}^{L-1}{h}_{\ell }{A}^{\ell }X,$$

$$h(A)=Vh(\Lambda )V^{-1}$$

$$h({\lambda }_i)=\sum _{\ell =0}^{L-1}{h}_{\ell }{\lambda }_i^{\ell }=c_i,$$

$${\pi }_i^{\ast }\propto \sqrt{c^2+{\Vert {\left(X_{\mathrm{o}}\right)}_i\Vert }_2^2}.$$

$$f_i(X)={\Vert {\left(h(A)X\right)}_i\Vert }_2^2,$$

$$h(A)1=\sum _{\ell =0}^{N-1}{h}_{\ell }{A}^{\ell }1=\sum _{\ell =0}^{N-1}{h}_{\ell }1=0,$$

$$h_{\mathrm{HH}}(A)$$

$${\left(h_{\mathrm{HH}}(A)X\right)}_i={\mathbf{x}}_i-\sum _{j\in {\mathcal{N}}_i}{A}_{i,j}{\mathbf{x}}_j.$$

$${\pi }_i^{\ast }$$

$$\mathrm{Tr}\left(X^{T}LX\right)=\sum _i\sum _{j\in {\mathcal{N}}_i}{W}_{i,j}{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }_2^2.$$

$$f_i(X)=\sum _{j\in {\mathcal{N}}_i}{W}_{i,j}{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert }_2^2.$$

$$\text{vspace}-5mm{h}_{\mathrm{IL}}(A)$$

$$h_{\rm IL}(\lambda_{i})\ =\ \left\{\begin{array}[]{rl}1,&i\leq b;\\ 0,&\mbox{otherwise}.\end{array}\right.$$

$${\pi }_i^{\ast }$$

$$h_{\mathrm{HL}}(A)=I+\frac{1}{|{\lambda }_{\max }|}A$$

$${\pi }_i^{\ast }$$

$${\mathbb{E}}_{\Psi \sim \pi }{\left({\Psi }^T\Psi f(X)\right)}_i={\mathbb{E}}_{\mathcal{M}}\left(\sum _{{\mathcal{M}}_j\in \mathcal{M}}{f}_{{\mathcal{M}}_j}(X){\delta }_{{\mathcal{M}}_j=i}\right)$$

Full text

I Introduction

With the recent development of 3D sensing technologies, 3D point clouds have become an important and practical representation of 3D objects and surrounding environments in many applications, such as virtual reality, mobile mapping, scanning of historical artifacts, 3D printing and digital elevation models [1]. A large number of 3D points on an object’s surface measured by a sensing device are called a 3D point cloud. Other than 3D coordinates, a 3D point cloud may also comprise some attributes, such as color, temperature and texture. Based on storage order and spatial connectivity between 3D points, there are two types of point clouds: organized point clouds and unorganized point clouds [2]. 3D points collected by a camera-like 3D sensor or a 3D laser scanner are typically arranged on a grid, like pixels in an image; we call those point clouds organized. For complex objects, we need to scan these objects from multiple view points and merge all collected points, which intermingles the indices of 3D points; we call those point clouds unorganized. It is easier to process an organized point cloud than an unorganized point cloud as the underlying grid produces a natural spatial connectivity and reflects the order of sensing. To make it general, we consider unorganized point clouds in this paper.

3D point cloud processing has become an important component in many 3D imaging and vision systems. It broadly includes compression [3, 4, 5, 6], visualization [7, 8], surface reconstruction [9, 10], rendering [11, 12], editing [13, 14] and feature extraction [15, 16, 17, 18, 19]. A challenge in 3D point cloud processing is how to handle a large number of incoming 3D points [20, 21]. In many applications, such as digital documentation of historical buildings and terrain visualization, we need to store billions of incoming 3D points; additionally, real-time sensing systems generate millions of data points per second. A large-scale point cloud makes storage and subsequent processing inefficient.

To solve this problem, an approach is to consider efficient data structures to represent 3D point clouds. For example, [22, 23] partitions the 3D space into voxels and discretizes point clouds over voxels; a drawback is that to achieve a fine resolution, a dense grid is required, which causes space inefficiency. [24, 25] presents an octree representation of point clouds, which is space efficient, but suffers from discretization errors. [26, 27] presents a probabilistic generative model to model the distribution of point clouds; drawbacks are that those parametric models may not capture the true surface, and it is inefficient to infer parameters in the probabilistic generative model.

Another approach is to consider reducing the number of points through mesh simplification. The main idea is to construct a triangular or polygonal mesh for 3D point clouds, where nodes are 3D points (need not be from the input points) and edges are connectivities between those points respecting certain restrictions (e.g., e.g. belonging to a manifold). The mesh is simplified by reducing the number of nodes or edges; that is, several nodes are merged into one node with local structure preserved. Surveys of many such methods can be found in [28, 29, 30]. Drawbacks of this approach are that mesh construction requires costly computation, and mesh simplification changes the positions of original points, which causes distortion.

In this paper, we consider resampling 3D point clouds; that is, we design application-dependent resampling strategies to preserve application-dependent information. For example, conventional contour detection in 3D point clouds requires careful and costly computation to obtain surface normals and classification models [31, 27]. We efficiently resample a small subset of points that is sensitive to the required contour information, making the subsequent processing cheaper without losing accuracy; see Figure 1 for an example. Since the original 3D point cloud is sampled from an object, we call this task resampling. This approach reduces the number of 3D points without changing the locations of original 3D points. After resampling, we unavoidably lose information in the original 3D point cloud.

The proposed method is rooted in rooted in graph signal processing, which is a framework to explore the interaction between signals and graph structure [32, 33]. We use a graph to capture local dependencies among points, representing a discrete version of the surface of an original object. The advantage of using a graph is to capture both local and global structure of point clouds. Each of the 3D coordinates and other attributes associated with 3D points is a graph signal indexed by the nodes of the underlying graph. We thus formulate a resampling problem as graph signal sampling. However, graph sampling methods usually select samples in a deterministic fashion, which solves nonconvex optimization problems to obtain samples sequentially and requires costly computation [34, 35, 36, 37]. To reduce the computational cost, we propose an efficient randomized resampling strategy to select a subset of points. The main idea is to generate subsamples according to a non-uniform resampling distribution, which is both fast and provably preserves application-dependent information in the original 3D point cloud.

We first propose a general feature-extraction based resampling framework. We use a general feature-extraction operator to represent application-dependent information. Based on this feature-extraction operator, we quantify the quality of resampling by using a simple, yet general reconstruction error, where we can derive the exact mean square error. We obtain the optimal resampling distribution by optimizing the mean square error. The proposed optimal resampling distribution is guaranteed to be shift/rotation/scale-invariant.

We next specify a feature extraction operator to be a graph filter and study the specific optimal resampling distributions based on all-pass, low-pass and high-pass graph filtering. In each case, we derive an optimal resampling distribution and validate the performance on both simulated and real data. We further combine all the proposed techniques into an efficient surface reconstruction system based on graph filter banks, which enables us to enhance features in a 3D point cloud.

We finally apply the proposed methods on three applications: large-scale visualization, accurate registration and robust shape modeling. In large-scale visualization, we use the proposed high-pass graph filtering based resampling strategy to highlight the contours of buildings and streets in a urban scene, which avoids saturation problems in visualization; in accurate registration, we use the proposed high-pass graph filtering based resampling strategy to extract the key points of a sofa, which makes the registration precise; in robust shape modeling, we use the proposed low-pass graph filtering based resampling strategy to reconstuct a surface, which makes the reconstruction robust to noise. The performances in those three applications validate the effectiveness and efficiency of the proposed resampling methods.

Contributions. This paper considers a widely-used task from a novel theoretical perspective. As a preprocessing step, resampling a large-scale 3D point cloud uniformly is widely used in many tasks of large-scale 3D point cloud processing and many commercial softwares; however, people treat this step heuristically. This paper considers resampling 3D points from a theoretical signal processing perspective. For example, our theory shows that uniform resampling is the optimal resampling distribution when all 3D points are associated with the same feature values. The main contributions of the paper are as follows: We propose

• •

a novel theoretical resampling framework for 3D point clouds with exact mean square error and optimal resampling distribution;

• •

a novel feature-extraction operator for 3D point clouds based on graph filtering;

• •

extensive empirical studies of the proposed resampling strategies on both simulated data and real point clouds.

This paper also points out many possible future directions of 3D point cloud processing, such as efficient 3D point cloud compression system based on graph filter banks, surface reconstruction based on arbitrary graphs and robust metric to evaluate the visualization quality of a 3D point cloud.

Outline of the paper. Section II formulates the resampling problem and briefly reviews graph signal processing. Section III proposes a resampling framework based on general feature-extraction operator and Section IV considers a graph filter as a specific feature-extraction operator. Three applications are presented in Section V. Section VI concludes the paper and provides pointers to future directions.

II Problem Formulation

In this section, we cover the background material necessary for the rest of the paper. We start with formulating a task of resampling a 3D point cloud. We then introduce graph signal processing, which lays a foundation for our proposed methods.

II-A Resampling a Point Cloud

We consider a matrix representation of a point cloud with $N$ points and $K$ attributes,

[TABLE]

where $\mathbf{s}_{i}\in\mathbb{R}^{N}$ represents the $i$th attribute and $\mathbf{x}_{i}\in\mathbb{R}^{K}$ represents the $i$th point. Depending on the sensing device, attributes can be 3D coordinates, RGB colors, textures, and many others. To distinguish 3D coordinates from other attributes, $\operatorname{\operatorname{X}_{\rm c}}\in\mathbb{R}^{N\times 3}$ represents 3D coordinates and $\operatorname{\operatorname{X}_{\rm o}}\in\mathbb{R}^{N\times(K-3)}$ represents other attributes.

The number of points $N$ is usually large. For example, a 3D scan of a building usually needs billions of 3D points. It is challenging to work with such a large-scale point cloud from both storage and data analysis perspectives. In many applications, however, we are interested in a subset of 3D points with particular properties, such as key points in point cloud registration and contour points in contour detection. To reduce the storage and computational cost, we consider resampling a subset of representative 3D points from the original 3D point cloud to reduce the scale. The procedure of resampling is to resample $M~{}(M<N)$ points from a point cloud, or select $M$ rows from the point cloud matrix $\operatorname{X}$. The resampled point cloud is $\operatorname{X}_{\mathcal{M}}\ =\ \Psi\operatorname{X}\ \in\mathbb{R}^{M\times K},$ where $\mathcal{M}=(\mathcal{M}_{1},\ldots,\mathcal{M}_{M})$ denotes the sequence of resampled indices, called resampled set, $\mathcal{M}_{i}\in\{1,\ldots,N\}$ with $|\mathcal{M}|=M$ and the resampling operator $\Psi$ is a linear mapping from $\mathbb{R}^{N}$ to $\mathbb{R}^{M}$, defined as

[TABLE]

The efficiency of the proposed resampling strategy is critical. Since we work with a large-scale point cloud, we want to avoid expensive computation. To implement resampling in an efficient way, we consider a randomized resampling strategy. It means that the resampled indices are chosen according to a resampling distribution. Let $\{\pi_{i}\}_{i=1}^{N}$ be a series of resampling probabilities, where $\pi_{i}$ denotes the probability to select the $i$th sample in each random trial. Once the resampling distribution is chosen, it is efficient to generate samples. The goal here is to find a resampling distribution that preserves information in the original point cloud.

The invariant property of the proposed resampling strategy is also critical. When we shift, rotate or scale a point cloud, the intrinsic distribution of 3D points does not changed and the proposed resampling strategy should not change.

Definition 1**.**

A resampling strategy is shift-invariant when a sampling distribution $\pi$ is designed for a point cloud, $\operatorname{X}=\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix}$, then the same sampling distribution $\pi$ is designed for its shifted point cloud, $\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}+{\bf 1}\mathbf{a}^{T}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix}$ with $\mathbf{a}\in\mathbb{R}^{3}$.

Definition 2**.**

A resampling strategy is rotation-invariant when a sampling distribution $\pi$ is designed for a point cloud, $\operatorname{X}=\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix}$, then the same sampling distribution $\pi$ is designed for its rotated point cloud, $\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}\operatorname{R}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix}$, where $\operatorname{R}\in\mathbb{R}^{3\times 3}$ is a 3D rotation matrix.

Definition 3**.**

A resampling strategy is scale-invariant when a sampling distribution $\pi$ is designed for a point cloud, $\operatorname{X}=\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix}$, then the same sampling distribution $\pi$ is designed for its rotated point cloud, $\begin{bmatrix}c\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix}$, where constant $c>0$.

Our aim is to guarantee that the proposed resampling strategy is shift, rotation and scale invariant.

II-B Graph Signal Processing for Point Clouds

A graph is a natural and efficient way to represent a 3D point cloud because it represents a discretized version of an original surface. In computer graphics, polygon meshes, as a class of graphs with particular connectivity restrictions, are extensively used to represent the shape of an object [38]; however, mesh construction usually requires sophisticated geometry analysis, such as calculating surface normals, and the mesh representation may not be the most suitable representation for analyzing point clouds because of connectivity restrictions. Here we extend polygon meshes to general graphs by relaxing the connectivity restrictions. Such graphs are easier to construct and are flexible to capture geometry information.

Graph Construction. We construct a general graph of a point cloud by encoding the local geometry information through an adjacency matrix $\operatorname{W}\in\mathbb{R}^{N\times N}$. Let $\mathbf{x}^{\rm(c)}_{i}\in\mathbb{R}^{3}$ be the 3D coordinates of the $i$th point; that is, the $i$th row of $\operatorname{\operatorname{X}_{\rm c}}$. The edge weight between two points $\mathbf{x}^{\rm(c)}_{i}$ and $\mathbf{x}^{\rm(c)}_{j}$ is

[TABLE]

where variance $\sigma$ and threshold $\tau$ are parameters. Equation (3) shows that when the Euclidean distance of two points is smaller than a threshold $\tau$, we connect these two points by an edge and the edge weight depends on the similarity of two points in the 3D space. The weighted degree matrix $\operatorname{D}$ is a diagonal matrix with diagonal element $\operatorname{D}_{i,i}=\sum_{j}\operatorname{W}_{i,j}$ reflecting the density around the $i$th point. This graph is approximately a discrete representation of the original surface and can be efficiently constructed via a tree data structure, such as octree [24, 25]. Here we only use the 3D coordinates to construct a graph, but it is also feasible to take other attributes into account (3). Given this graph, the attributes of point clouds are called graph signals. For example, an attribute $\mathbf{s}$ in (1) is a signal index by the graph.

Graph Filtering. A graph filter is a system that takes a graph signal as an input and produces another graph signal as an output. Let $\operatorname{A}\in\mathbb{R}^{N\times N}$ be a graph shift operator, which is the most elementary nontrivial graph filter. Some common choice of a graph shift operator is the adjacency matrix $\operatorname{W}$ (3), the transition matrix $\operatorname{D}^{-1}\operatorname{W}$, the graph Laplacian matrix $\operatorname{D}-\operatorname{W}$, and many other structure-related matrices. The graph shift replaces the signal value at a node with a weighted linear combination of values at its neighbors; that is, $\mathbf{y}\ =\ \operatorname{A}\mathbf{s}\ \in\ \mathbb{R}^{N},$ where $\mathbf{s}\in\mathbb{R}^{N}$ is an input graph signal (an attribute of a point cloud). Every linear, shift-invariant graph filter is a polynomial in the graph shift [32]

[TABLE]

where $h_{\ell}(\ell=0,1,\ldots,L-1)$ are filter coefficients and $L$ is the length of this graph filter. Its output is given by the matrix-vector product $\mathbf{y}\ =\ h(\operatorname{A})\mathbf{s}\ \in\ \mathbb{R}^{N}.$

Graph Fourier Transform. The eigendecomposition of a graph shift operator $\operatorname{A}$ is [39]

[TABLE]

where the eigenvectors of $\operatorname{A}$ form the columns of matrix $\operatorname{V}$, and the eigenvalue matrix $\Lambda\in\mathbb{R}^{N\times N}$ is the diagonal matrix of corresponding eigenvalues $\lambda_{1},\,\ldots,\,\lambda_{N}$ of $\operatorname{A}$ ($\lambda_{1}\geq\lambda_{2}\geq\ldots,\,\geq\lambda_{N})$. These eigenvalues represent frequencies on the graph [39] where $\lambda_{1}$ is the lowest frequency and $\lambda_{N}$ is the highest frequency. Correspondingly, $\mathbf{v}_{1}$ captures the smallest variation on the graph and $\mathbf{v}_{N}$ captures the highest variation on the graph. $\operatorname{V}$ is also called graph Fourier basis. The graph Fourier transform of a graph signal $\mathbf{s}\in\mathbb{R}^{N}$ is

[TABLE]

The inverse graph Fourier transform is $\mathbf{s}\ =\ \operatorname{V}\widehat{\mathbf{s}}=\sum_{k=1}^{N}\widehat{s}_{k}\mathbf{v}_{k},$ where $\mathbf{v}_{k}$ is the $k$th column of $\operatorname{V}$ and $\widehat{s}_{k}$ is the $k$th component in $\widehat{\mathbf{s}}$. The vector $\widehat{\mathbf{s}}$ in (6) represents the signal’s expansion in the eigenvector basis and describes the frequency components of the graph signal $\mathbf{s}$. The inverse graph Fourier transform reconstructs the graph signal by combining graph frequency components.

III Resampling based on Feature Extraction

During resampling, we reduce the number of points and unavoidably lose information in a point cloud. Our goal is to design an application-dependent resampling strategy, preserving selected information depending on particular needs. Those information are described by features. When detecting contours, we usually need careful and intensive computation, such as calculating surface normals and classifying points [31, 27]. Instead of working with a large number of points, we consider efficiently sampling a small subset of points that captures the required contour information, making the subsequent computation much cheaper without losing contour information. We also need to guarantee that the proposed resampling strategy is shift/rotation/scale-invariant for robustness. We will show that some features naturally provide invariance and other may not. We will handle the invariance by considering a general objective function.

III-A Feature-Extraction based Formulation

Let $f(\cdot)$ be a feature-extraction operator that extracts targeted information from a point cloud according to particular needs; that is, the features $f(\operatorname{X})\in\mathbb{R}^{N\times K}$ are extracted from a point cloud $\operatorname{X}\in\mathbb{R}^{N\times K}$111 For simplicity, we consider the number of features to be the same as the number of attributes. The proposed method also works when the number of features and the number of attributes are different.. Depending on an application, those features can be edges, key points and flatness [16, 17, 18, 40, 19]. In this section, we consider feature-extraction operator at an abstract level and use graph filters to implement a feature-extraction operator in the next section.

To evaluate the performance of a resampling operator, we quantify how much features are lost during resampling; that is, we sample features, and then interpolate to get back original features. The features are considered to reflect the targeted information contained in each 3D point. The performance is better when the recovery error is smaller. Mathematically, we resample a point cloud $M$ times. At the $j$th step, we independently choose a point $\mathcal{M}_{j}=i$ with probability $\pi_{i}$. Let $\Psi\in\mathbb{R}^{M\times N}$ be the resampling operator (2) and $\operatorname{S}\in\mathbb{R}^{N\times N}$ be a diagonal rescaling matrix with $\operatorname{S}_{i,i}=1/\sqrt{M\pi_{i}}$. We quantify the performance of a resampling operator as follows:

[TABLE]

where $\left\|\cdot\right\|_{2}$ is the spectral norm. $\Psi^{T}\Psi\in\mathbb{R}^{N\times N}$ is a zero-padding operator, which a diagonal matrix with diagonal elements $(\Psi^{T}\Psi)_{i,i}>0$ when the $i$th point is sampled, and [math], otherwise. The zero-padding operator $\Psi^{T}\Psi$ ensures the resampled points and the original point cloud have the same size. $\operatorname{S}$ is used to compensate non-uniform weights during resampling. $\operatorname{S}\Psi^{T}$ is the most naive interpolation operator that reconstructs the original feature $f(\operatorname{X})$ from its resampled version $\Psi f(\operatorname{X})$ and $\operatorname{S}\Psi^{T}\Psi f(\operatorname{X})$ represents the preserved features after resampling in a zero-padding form. Lemma 1 shows that $\operatorname{S}$ aids to provide an unbiased estimator.

Lemma 1**.**

Let $f(\operatorname{X})\in\mathbb{R}^{N\times K}$ be features extracted from a point cloud $\operatorname{X}$. Then,

[TABLE]

where $\mathbb{E}_{\Psi\sim\pi}$ means the expectation over samples, which are generated from a distribution $\Pi$ independently and randomly, and $\odot$ is row-wise multiplication.

The proof is shown in Appendix A-A.

The evaluation metric $D_{f(\operatorname{X})}(\Psi)$ measures the reconstruction error; that is, how much feature information is lost after resampling without using sophisticated interpolation operator. When $D_{f(\operatorname{X})}(\Psi)$ is small, preserved features after resampling are close to the original features, meaning that little information is lost. The expectation $\mathbb{E}_{\Psi\sim\pi}\left(D_{f(\operatorname{X})}(\Psi)\right)$ is the expected error caused by resampling and quantifies the performance of a resampling distribution $\pi$. Our goal is to minimize $\mathbb{E}_{\Psi\sim\pi}\left(D_{f(\operatorname{X})}(\Psi)\right)$ over $\pi$ to obtain an optimal resampling distribution in terms of preserving features $f(\operatorname{X})$. We now derive the mean square error of the objective function (7).

Theorem 1**.**

The mean square error of the objective function (7) is

[TABLE]

where $\operatorname{Q}\in\mathbb{R}^{N\times N}$ is a diagonal matrix with $\operatorname{Q}_{i,i}=1/\pi_{i}-1$.

The proof is shown in Appendix A-B.

We now consider the invariance property of resampling. The sufficient condition for the shift/rotation/scale-invariance of a resampling strategy is that the evaluation metric (7) be shift/ rotation/scale-invariance. Recall that a 3D point cloud is $\operatorname{X}=\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix}$, where $\operatorname{\operatorname{X}_{\rm c}}\in\mathbb{R}^{N\times 3}$ represents 3D coordinates and $\operatorname{\operatorname{X}_{\rm o}}\in\mathbb{R}^{N\times(K-3)}$ represents other attributes.

Definition 4**.**

A feature-extraction operator $f(\cdot)$ is shift-invariant when the features extracted from a point cloud and its shifted version are same; that is, $f(\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix})=f(\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}+{\bf 1}\mathbf{a}^{T}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix})$ with shift $\mathbf{a}\in\mathbb{R}^{3}$.

Definition 5**.**

A feature-extraction operator $f(\cdot)$ is rotation-invariant when the features extracted from a point cloud and its rotated version are same; that is, $f(\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix})=f(\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}\operatorname{R}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix})$ with $\operatorname{R}\in\mathbb{R}^{3\times 3}$ is a 3D rotation matrix.

Definition 6**.**

A feature-extraction operator $f(\cdot)$ is scale-invariant when features extracted from a point cloud and its scaled version are same; that is, $f(\begin{bmatrix}\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix})=f(\begin{bmatrix}c\operatorname{\operatorname{X}_{\rm c}}&\operatorname{\operatorname{X}_{\rm o}}\end{bmatrix})$ with constant $c>0$.

When $f(\cdot)$ is shift/rotation/scale-invariant, (7) does not change through shifting, rotating or scaling, leading to a shift/rotation/scale-invariant resampling strategy and it is sufficient to minimize $\mathbb{E}_{\Psi\sim\pi}\left(D_{f(\operatorname{X})}(\Psi)\right)$ to obtain a resampling strategy; however, when $f(\cdot)$ is shift/rotation/scale-variance, (7) may change through shifting, rotating or scaling, leading to a shift/rotation/scale-variant resampling strategy.

To handle shift variance, we can recenter a point cloud to the origin before processing; that is, we normalize the mean of 3D coordinates to zeros. To handle scale variance, we can normalize the magnitude of the 3D coordinates before processing; that is, we normalize the spectral norm $\left\|\operatorname{X}_{c}\right\|_{2}=c$ with constant $c>0$. The choice of $c$ depends on users’ preference and we will show that $c$ is a trade-off between 3D coordinates and the values of other attributes. From now on, we first recenter a point cloud to the origin and then normalize its magnitude to guarantee the shift/scale invariance of any 3D point cloud.

To handle rotation variance of $f(\cdot)$, we consider the following evaluation metric:

[TABLE]

where constant $c=\left\|\operatorname{X}_{c}\right\|_{2}$ is the normalized spectral norm of 3D coordinates.

Unlike $D_{f(\operatorname{X})}(\Psi)$ (7), to remove the influence of rotation, the evaluation metric $D_{f}(\Psi)$ considers the worst possible reconstruction error caused by rotation. In (III-A), we consider 3D coordinates as variables due to rotation. We constrain the spectral norm of 3D coordinates because a rotation matrix is orthornormal and the spectral norm of 3D coordinates does not change during rotation. We then minimize $\mathbb{E}_{\Psi\sim\pi}\left(D_{f}(\Psi)\right)$ to obtain a rotation-invariant resampling strategy even when $f(\cdot)$ is rotation-variant.

For simplicity, we perform derivation for only linear feature-extraction operators. A linear feature-extraction operator $f(\cdot)$ is of the form of $f(\operatorname{X})=\operatorname{F}\operatorname{X}$, where $\operatorname{X}$ is a 3D point cloud and $\operatorname{F}\in\mathbb{R}^{N\times N}$ is a feature-extraction matrix.

Theorem 2**.**

Let $f(\cdot)$ be a rotation-varying linear feature-extraction operator, where $f(\operatorname{X})=\operatorname{F}\operatorname{X}$ with $\operatorname{F}\in\mathbb{R}^{N\times N}$. The exact form of $\mathbb{E}_{\Psi\sim\pi}D_{f}(\Psi)$ is

[TABLE]

where $\operatorname{Q}\in\mathbb{R}^{N\times N}$ is a diagonal matrix with $\operatorname{Q}_{i,i}=1/\pi_{i}-1$.

The proof is shown in Appendix A-C.

III-B Optimal Resampling Distribution

We now derive the optimal resampling distributions by minimizing the reconstruction error. For a rotation-invariant feature-extraction operator, we minimize (8).

Theorem 3**.**

Let $f(\cdot)$ be a rotation-invariant feature-extraction operator. The corresponding optimal resampling strategy $\pi^{*}$ is,

[TABLE]

where $f_{i}(\operatorname{X})\in\mathbb{R}^{K}$ is the $i$th row of $f(\operatorname{X})$.

The proof is shown in Appendix A-D. We see that the optimal resampling distribution is proportional to the magnitude of features; that is, points associated with high magnitudes have high probability to be selected, while points associated with small magnitudes have small probability to be selected. The intuition is that the response after the feature-exaction operator reflects the information contained in each 3D point and determines the resampling probability of each 3D point.

For a rotation-variant linear feature-extraction operator, we minimize (10).

Theorem 4**.**

Let $f(\cdot)$ be a rotation-variant linear feature-extraction operator, where $f(\operatorname{X})=\operatorname{F}\operatorname{X}$ with $\operatorname{F}\in\mathbb{R}^{N\times N}$. The corresponding optimal resampling strategy $\pi^{*}$ is,

[TABLE]

where constant $c=\left\|\operatorname{\operatorname{X}_{\rm c}}\right\|_{2}$, $\operatorname{F}_{i}$ is the $i$th row of $\operatorname{F}$ and $\left(\operatorname{F}\operatorname{\operatorname{X}_{\rm o}}\right)_{i}$ is the $i$th row of $\operatorname{F}\operatorname{\operatorname{X}_{\rm o}}$.

The proof is shown in Appendix A-E. We see that the optimal resampling distribution is also proportional to the magnitude of features. The feature comes from two sources: 3D coordinates and the other attributes. The tuning parameter $c$ in (12) is the normalized spectral norm used to remove the scale variance. The choice of $c$ trade-offs the contribution from 3D coordinates and the other attributes.

IV Resampling based on Graph Filtering

The previous section studied resampling based on an arbitrary feature-extraction operator. In this section, we design graph filters to efficiently extract features from a point cloud. Let features extracted from a point cloud $\operatorname{X}$ be

[TABLE]

which follows from the definition of graph filters (4). Since a graph filter is a linear operator, the corresponding optimal resampling distribution follows from the results in Theorems 3 and 4 by replacing $\operatorname{F}=\sum_{\ell=0}^{L-1}h_{\ell}\operatorname{A}^{\ell}$. All graph filtering-based feature-extraction operators are scale-variant due to linearity. As discussed earlier, we can normalize the spectral norm of a 3D coordinates to handle this issue. We thus will not discuss scale invariance in this section. We will see that by carefully using the graph shift operator $\operatorname{A}$ and filter coefficients $h_{i}$s, a graph filtering-based feature-extraction operator may be shift or rotation varying.

Similarly to filter design in classical signal processing, we design a graph filter either in the graph vertex domain or in the graph spectral domain. In the graph vertex domain, for each point, a graph filter averages the attributes of its local points. For example, the output of the $i$th point, $f_{i}(\operatorname{X})=\sum_{\ell=0}^{L-1}h_{\ell}\left(\operatorname{A}^{\ell}\operatorname{X}\right)_{i}$ is a weighted average of the attributes of points that are within $L$ hops away from the $i$th point. The $\ell$th graph filter coefficient, $h_{\ell}$, quantifies the contribution from the $\ell$th-hop neighbors. We design the filter coefficients to change the weights in local averaging.

In the graph spectral domain, we first design a graph spectrum distribution and then use graph filter coefficients to fit this distribution. For example, a graph filter with length $L$ is

[TABLE]

where $\operatorname{V}$ is the graph Fourier basis and $\lambda_{i}$ are graph frequencies (5). When we want the response of the $i$th graph frequency to be $c_{i}$, we set

[TABLE]

and solve a set of linear equations to obtain the graph filter coefficients $h_{\ell}$. It is also possible to use the Chebyshev polynomial to design graph filter coefficients [41]. We now consider some special cases of graph filters.

IV-A All-pass Graph Filtering

Let $h(\lambda_{i})=1$; that is, $h(\operatorname{A})=\operatorname{I}$ is an identity matrix with $h_{0}=1$ and $h_{i}=0$ for $i=1,\ldots,L-1$. The intuition behind this setting is that the original point cloud is trustworthy and all points are uniformly sampled from an object without noise, reflecting the true geometric structure of the object. We want to preserve all the information and the features are thus the original attributes themselves. Since $f(\operatorname{X})=\operatorname{X}$, the feature-extraction operator $f(\cdot)$ is rotation-variant. Based on Theorem 4, the optimal resampling strategy is

[TABLE]

Here the feature-extraction matrix $\operatorname{F}$ in (11) is an identity matrix and the norm of each row of $\operatorname{F}$ is 1. When we only preserve 3D coordinates, we ignore the term of $\operatorname{\operatorname{X}_{\rm o}}$ and obtain a constant resampling probability for each point, meaning that uniform resampling is the optimal resampling strategy to preserve the overall geometry information.

IV-B High-pass Graph Filtering

In image processing, a high-pass filter is used to extract edges and contours. Similarly, we use a high-pass graph filter to extract contours in a point cloud. Here we only consider the 3D coordinates as attributes ($\operatorname{X}=\operatorname{\operatorname{X}_{\rm c}}=\mathbb{R}^{N\times 3}$), but the proposed method can be easily extended to other attributes.

A critical question is how to define contours in a 3D point cloud. We consider that contour points break the trend formed by its neighboring points and bring innovation. Many previous works need sophisticated geometry-related computation, such as surface normal, to detect contours [31]. Instead of measuring sophisticated geometry properties, we describe the possibility of being a contour point by the local variation on graphs, which is the response of high-pass graph filtering. The corresponding local variation of the $i$th point is

[TABLE]

where $h(\operatorname{A})$ is a high-pass graph filter. The local variation $f(\operatorname{X})\in\mathbb{R}^{N}$ quantifies the energy of response after high-pass graph filtering. The intuition behind this is that when the local variation of a point is high, its 3D coordinates cannot be well approximated from the 3D coordinates of its neighboring points; in other words, this point bring innovation by breaking the trend formed by its neighboring points and has a high possibility of being a contour point.

The following theorem shows that in general the local variation is rotation invariant, but shift variant.

Theorem 5**.**

Let $f(\operatorname{X})={\rm diag}\left(h(\operatorname{A})\operatorname{X}\operatorname{X}^{T}h(\operatorname{A})^{T}\right)\in\mathbb{R}^{N}$, where diag$(\cdot)$ extracts the diagonal elements. $f(\operatorname{X})$ is rotation invariant and shift invariant unless $h(\operatorname{A}){\bf 1}=\mathbb{0}\in\mathbb{R}^{N}$.

The proof is shown in Appendix A-F.

To guarantee that local variation is naturally shift invariant without recentering a 3D point cloud, we simply use a transition matrix as a graph shift operator; that is, $\operatorname{A}=\operatorname{D}^{-1}\operatorname{W}$, where $\operatorname{D}$ is the diagonal degree matrix. The reason is that ${\bf 1}\in\mathbb{R}^{N}$ is the eigenvector of a transition matrix, $\operatorname{A}{\bf 1}=\operatorname{D}^{-1}\operatorname{W}{\bf 1}={\bf 1}$. Thus,

[TABLE]

when $\sum_{\ell=0}^{N-1}h_{\ell}=0$. A simple design is a Haar-like high-pass graph filter

[TABLE]

Note that $\lambda_{\rm max}=\max_{i}|\lambda_{i}|=1$, where $\lambda_{i}$ are eigenvalues of $\operatorname{A}$, because the graph shift operator is a transition matrix. In this case, $h_{0}=1,h_{1}=-1$ and $h_{i}=0$ for all $i>1$, $\sum_{\ell=0}^{N-1}h_{\ell}=0$. Thus, a Haar-like high-pass graph filter is both shift and rotation invariant. The graph frequency response of a Haar-like high-pass graph filter is $h_{\rm HH}(\lambda_{i})=1-\lambda_{i}$. Since the eigenvalues are ordered descendingly, we have $1-\lambda_{i}\leq 1-\lambda_{i+1}$, meaning low frequency response attenuates and high frequency response amplifies.

In the graph vertex domain, the response of the $i$th point is

[TABLE]

Because $\operatorname{A}$ is a transition matrix, $\sum_{j\in\mathcal{N}_{i}}\operatorname{A}_{i,j}=1$ and $h_{\rm HH}(\operatorname{A})$ compares the difference between a point and the convex combination of its neighbors. The geometry interpretation of the proposed local variation is the Euclidean distance between the original point and the convex combination of its neighbors, reflecting how much information we know about a point from its neighbors. When the local variation of a point is large, the Euclidean distance between this point and the convex combination of its neighbors is long and this point provides a large amount of variation.

We can verify the proposed local variation on some simple examples.

Example 1**.**

When a point cloud forms a 3D line, two endpoints belong to the contour.

Example 2**.**

When a point cloud forms a 3D polygon/polyhedron, the vertices (corner points) and the edges (line segment connecting two adjacent vertices) belong to the contour.

Example 3**.**

When a point cloud forms a 3D circle/sphere, there is no contour.

When the points are uniformly spread along the defined shape, the proposed local variation (14) satisfies Examples 1, 2 and 3 from the geometric perspective. In Figure 2 (a), Point 2 is the convex combination of Points 1 and 3, and the local variation of Point 2 is thus zero. However, Point 4 is not the convex combination of Points 3 and 5 and the length of the red line indicates the local variation of Point 4. Only Points 1, 4 and 7 have nonzero local variation, which is what we expect. In Figure 2 (b), all the nodes are evenly spread on a circle and have the same amount of variation, which is represented as a red line. Similar arguments show that the proposed local variation (14) satisfies Examples 1, 2 and 3.

The feature-extraction operator $f(\operatorname{X})=\left\|h_{\rm HH}(\operatorname{A})\operatorname{X}\right\|_{F}^{2}$ is shift and rotation-invariant. Based on Theorem 3, the optimal resampling distribution is

[TABLE]

where $\operatorname{A}=\operatorname{D}^{-1}\operatorname{W}$ is a transition matrix.

Note that the graph Laplacian matrix is commonly used to measure variations. Let $\operatorname{L}=\operatorname{D}-\operatorname{W}\in\mathbb{R}^{N\times N}$ be a graph Laplacian matrix. The graph Laplacian based total variation is

[TABLE]

where $\mathcal{N}_{i}$ is the neighbors of the $i$th node and the variation contributed by the $i$th point is

[TABLE]

The variation here is defined based on the accumulation of pairwise differences. We call (18) pairwise difference based local variation. The pairwise difference based local variation cannot capture geometry change and violates Example 2. We show a counter example in Figure 3. The points are uniformly spread along the faces of a cube and Figure 3 shows two faces. Each point connects to its adjacent four points with the same edge weight. The pairwise difference based local variations of all the points are the same, which means that there is no contour in this point cloud. However, the black arrow points to a point that should be a contour point.

Experimental Validations. Figure 4 compares the Haar-like high-pass graph filtering based local variation (14) (second column) with that computed from the difference of normals (DoN) method (first column)  [42] which is used to analyze point clouds for segmentation and contour detection. As a contour detection technique, DoN computes the difference between surface normals calculated at two scales. In each plot, we highlight the points that have top $10\%$ largest DoN scores or local variations. In Figure 4 (a), we see that DoN cannot find the boundary in the plane because the surface normal does not change. The performance of DoN is also sensitive to predeisgned radius. For example, the difference of normals cannot capture precise contours in the hinge. On the other hand, local variation captures all the contours precisely in Figure 4 (b). We see similar results in Figures 4 (c), (d), (e) and (f). Further, difference of normals needs to compute the first principle component of the neighboring points for each 3D point, which is computationally inefficient. The local variation only involves a sparse matrix and vector multiplication, which is computationally efficient.

Figure 5 shows the local variation based resampling distribution on some examples of the point cloud, including hinge, cone, table, chair and trash container. The first column shows the original point clouds; the second and third rows show the resampled versions with respect to two local variations: pairwise difference based local variation (18) and Haar-like high-pass graph filtering based local variation (14). Two resampled versions have the same number of points, which is $10\%$ of points in the original point cloud.

For two simulated objects, the hinge and the cone (first two rows), the pairwise difference based local variation (18) fails to detect contour and the Haar-like high-pass graph filtering based local variation (14) detects all the contours. For the real objects, the Haar-like high-pass graph filtering based resampling (14) also outperform the pairwise difference based local variation (18). In summary, the Haar-like high-pass graph filtering based local variation (14) shows the contours of objects by using only $10\%$ of points.

The high-pass graph filtering based resampling strategy can be easily extended to detect transient changes in other attributes. Figure 6 (a) simulates a hinge with two different textures. The points in black have the same texture with value 0 and the points indicated by a green circle have a different texture with value 1. We put the texture as a new attribute and the point cloud matrix $\operatorname{X}\in\mathbb{R}^{N\times 4}$, where the first three columns are 3D coordinates and the fourth column is the texture. We resample $10\%$ of points based on the high-pass graph filtering based local variation (14). Figure 6 (b) shows the resamped point cloud, which clearly detects both the geometric contour and the texture contour.

IV-C Low-pass Graph Filtering

In classical signal processing, a low-pass filter is used to capture rough shape of a smooth signal and reduce noise. Similarly, we use a low-pass graph filter to capture rough shape of a point cloud and reduce sampling noise during resampling. Since we use the 3D coordinates of points to construct a graph (3), the 3D coordinates are naturally smooth on this graph, meaning that two adjacent points in the graph have similar coordinates in the 3D space. When a 3D point cloud is corrupted by noises and outliers, a low-pass graph filter, as a denoising operator, uses local neighboring information to approximate a true position for each point. Since the output after low-pass graph filtering is a denoised version of the original point cloud, it is more appropriate to resample from denoised points than original points.

IV-C1 Ideal low-pass graph filter

A straightforward choice is an ideal low-pass graph filter, which completely eliminates all graph frequencies above a given graph frequency while passing those below unchanged. An ideal low-pass graph filter with bandwidth $b$ is

[TABLE]

where $\operatorname{V}_{(b)}$ is the first $b$ columns of $\operatorname{V}$, and the graph frequency response is

[TABLE]

The ideal low-pass graph filter $h_{\rm IL}$ projects an input graph signal onto a bandlimited subspace [34] and $h_{\rm IL}(\operatorname{A})\mathbf{s}$ is a bandlimited approximation of the original graph signal $\mathbf{s}$. We show an example in Figure 7. Figure 7 (b), (c) and (d) shows that the bandlimited approximation of the 3D coordinates of a teapot gets better when the bandwidth $b$ increases. We see that the bandwidth influences the shape of the teapot rapidly: with ten graph frequencies, we only obtain a rough structure of the teapot. Figure 7 (e) shows that the main energy is concentrated in the low-pass graph frequency band.

The feature-extraction operator $f(\operatorname{X})=\operatorname{V}_{(b)}\operatorname{V}_{(b)}^{T}\operatorname{X}$ is shift and rotation-varying. Based on Theorem 4, the corresponding optimal resampling strategy is

[TABLE]

where $\mathbf{v}_{i}\in\mathbb{R}^{b}$ is the $i$th row of $\operatorname{V}_{(b)}$.

A direct way to obtain $\left\|\mathbf{v}_{i}\right\|_{2}$ requires the truncated eigendecomposition (6), whose computational cost is $O(Nb^{2})$, where $b$ is the bandwidth. It is potentially possible to approximate the leverage scores through a fast algorithm [43, 44], where we use randomized techniques to avoid the eigendecomposition and the computational cost is $O(Nb\log(N))$. Another way to leverage computation is to partition a graph into several subgraphs and obtain leverage scores in each subgraph.

IV-C2 Haar-like low-pass graph filter

Another simple choice is Haar-like low-pass graph filter; that is,

[TABLE]

where $\lambda_{\rm max}=\max_{i}|\lambda_{i}|$ with $\lambda_{i}$ eigenvalues of $\operatorname{A}$. The normalization factor $\lambda_{\rm max}$ is presented to avoid the amplification of the magnitude. We denote $\operatorname{A}_{\rm norm}=\operatorname{A}/|\lambda_{\rm max}|$ for simplicity. The graph frequency response is $h_{\rm HL}(\lambda_{i})=1+\lambda_{i}/|\lambda_{\rm max}|$. Since the eigenvalues are ordered in a descending order, we have $1+\lambda_{i}\geq 1+\lambda_{i+1}$, meaning low frequency response amplifies and high frequency response attenuates.

In the graph vertex domain, the response of the $i$th point is $\left(h_{\rm HL}(\operatorname{A})\operatorname{X}\right)_{i}\ =\ \mathbf{x}_{i}+\sum_{j\in\mathcal{N}_{i}}(\operatorname{A}_{\rm norm})_{i,j}\mathbf{x}_{j},$ where $\mathcal{N}_{i}$ is the neighbors of the $i$th point. We see that $h_{\rm HL}(\operatorname{A})$ averages the attributes of each point and its neighbors to provide a smooth output.

The feature-extraction operator $f(\operatorname{X})=h_{\rm HL}(\operatorname{A})\operatorname{X}$ is shift and rotation-variant. Based on Theorem 4, the corresponding optimal resampling strategy is

[TABLE]

To obtain this optimal resampling distribution, we need to compute the largest magnitude eigenvalue $\lambda_{\rm max}$, which takes $O(N)$, and compute $\left\|\left(\operatorname{I}+\operatorname{A}_{\rm norm}\right)_{i}\right\|_{2}^{2}$ and $\left\|\left(\left(\operatorname{I}+\operatorname{A}_{\rm norm}\right)\operatorname{\operatorname{X}_{\rm o}}\right)_{i}\right\|_{2}^{2}$ for each row, which takes $O(\left\|{\rm vec}(\operatorname{A})\right\|_{0})$ with $\left\|{\rm vec}(\operatorname{A})\right\|_{0}$ the nonzero elements in the graph shift operator. We can avoid computing the largest magnitude by using a normalized adjacency matrix or a transition matrix as a graph shift operator. A normalized adjacency matrix is $\operatorname{D}^{-\frac{1}{2}}\operatorname{W}\operatorname{D}^{-\frac{1}{2}}$, where $\operatorname{D}$ is the diagonal degree matrix, and a transition matrix is obtained by normalizing the sum of each row of an adjacency matrix to be one; that is $\operatorname{D}^{-1}\operatorname{W}$. In both cases, the largest eigenvalue of a transition matrix is one, we thus have $\operatorname{A}=\operatorname{A}_{\rm norm}$.

Experimental Validations. We aim to use a low-pass graph filter to handle a noisy point cloud. Figure 8 (a) shows a point cloud of a fitness ball, which contains $62,235$ points collected from a Kinect device. In this noiseless case, the surface of the fitness can be modeled by a sphere. Figure 8 (b) fits a green sphere to the fitness ball 222Figure 8 (b) is generated from a public software CloudCompare.. The radius and the central point of this sphere is $0.318238$ and $\begin{bmatrix}0.0832627&0.190267&1.1725\end{bmatrix}$. To leverage the computation, we resample a subset of points and fit another sphere to the resample points. We want these two spheres generated by the original point cloud and the resampled point cloud to be similar.

In many real cases, the original points are collected with noise. To simulate the noisy case, we add the Gaussian noise with mean zeros and variance $0.02$ to each points. Figures 9 (a) and (b) show a noisy point cloud and its resampled version based on uniform resampling, respectively. Figures 9 (c) and (d) show a denoised point cloud and its resampled version, respectively. The denoised point cloud is obtained by the low-pass graph filtering (IV-C2) and the resampling strategy is based on (IV-C1). We fit a sphere to each of the four point clouds and the statistics are shown in Table I. The relative errors are shown in the parenthesis, which is defined as ${\rm Error}\ =\ \arrowvert(x-\hat{x})/x\arrowvert,$ where $x$ is the ground truth, and $\hat{x}$ is the estimation. The denoised ball and its resampled version outperform the noisy ball and the uniform-sampled version because the estimated radius and the central point is closer to the original radius and the central point. This validates that the proposed resampling strategy with low-pass graph filtering provides a robust shape modeling for noisy point clouds.

IV-D Graph Filter Banks

In classical signal processing, a filter bank is an array of band-pass filters that analyze an input signal in multiple subbands and synthesize the original signal from all the subbands [45, 46]. We use a similar idea to analyze a 3D point cloud: separate an input 3D point cloud into multiple components via different resampling operators, allowing us to enhance different components of a 3D point cloud. For example, we resample both contour points and noncontour points to reconstruct the original surfaces, but we need more contour points to emphasize contours.

Figure 10 shows a surface reconstruction system for a 3D point cloud based on graph filter banks. In the analysis part, we separate a 3D point cloud $\operatorname{X}$ into $k$ subbands. In each subband, the information preserved is determined by a specific graph filter and we resample a subset of 3D points according to (11) and (12). The number of samples in each subband is determined by a sampling ratio $\alpha$. We have flexibility to use either the original 3D points or the 3D points after graph filtering. In the synthesis part, we use the resampled points to reconstruct the surface. A literature review on surface reconstruction algorithms is shown in [47]. Since each surface reconstruction algorithm has its own specific set of assumptions, different surface reconstruction algorithms perform differently on the same set of 3D points.

We measure the overall performance of a surface reconstruction system by reconstruction error, which is the difference between the surface reconstructed from resampled points and the original surface. This leads to a rate-distortion like tradeoff: when we resample more points, we encode more bits and the reconstruction error is smaller; when we resample fewer points, we encode less bits and the reconstruction error is larger. The overall goal is: given an arbitrary tolerance of reconstruction error, we use as few samples as possible to reconstruct a surface by carefully choosing a graph filter and sampling ratio in each subband. Such a surface reconstruction system will benefit a 3D point cloud storage and compression because we only need to store a few resampled points. Since a surface reconstruction system is application-dependent, the design details are beyond the scope of this paper.

V Applications

In this section, we apply the proposed resampling strategies to accurate registration. In this task, we use the proposed resampling strategy (16) to make two point clouds registered efficiently and accurately.

Figure 11 (a) shows a point cloud of a sofa, which contains $1,204,055$ points collected from a Kinect based SLAM system [40]. As shown in Figure 11 (a), we split the original point cloud into two overlapping point clouds marked in red and blue, respectively. We intentionally shift and rotate the red part. The task is to invert the process and retrieve the shift and rotation. We use the iterative closest point (ICP) algorithm to register two point clouds, which is a standard algorithm to rotate and shift different scans into a consistent coordinate frame [48]. The ICP algorithm iteratively revises the rigid body transformation (combination of shift and rotation) needed to minimize the distance from the source to the reference point cloud. Figures 11 (b) and (c) show the registered sofa and the details of the overlapping part after registration, respectively. We see that the registration process recovers the overall structure of the original point cloud, but still leaves some mismatch in a detailed level.

Since it is inefficient to register two large-scale point clouds, we want to resample a subset of 3D points from each point cloud and implement registration. We will compare the registration performance between uniformly resampled point cloud and high-pass graph filtering based resampled point cloud. Note that high-pass graph filtering based resampling can enhance the contours and key points. Figures 11 (d) and (g) show the resampled point clouds based on uniform resampling and high-pass graph filtering based resampling, respectively. Two resampled versions have the same number of points, which is $5\%$ of points in the original point cloud. We see that Figures 11 (g) shows more contours than Figures 11 (d). Based on the uniformly resampled version Figures 11 (d), Figures 11 (e) and (f) show the registered sofa and the details of the overlapping part after registration, respectively. Based on the contour-enhanced resampled version Figures 11 (g), Figures 11 (h) and (i) show the registered sofa and the details of the overlapping part after registration, respectively. We see that the registration based on high-pass graph filtering based resampling precisely recovers the original point cloud, even in a detailed level. The intuition is that the high-pass graph filters enhance the contours, which make sharper match between the sources and targets, and thus the registration becomes easier.

The quantitative results are shown in Table II, where the first column shows the root mean square error (RMSE); the second column shows the shift error; The third column shows the rotation error. Specially, ${\rm RMSE}=\sqrt{\sum_{i=1}^{N}\min_{j=1,\ldots,N}\left\|\widehat{\mathbf{x}}_{i}-\mathbf{x}_{j}\right\|_{2}^{2}}$, ${\rm Error}_{\rm shift}=\left\|\widehat{\mathbf{a}}-\mathbf{a}\right\|_{2}$ and ${\rm Error}_{\rm rotation}=\left\|\widehat{\operatorname{R}}-\operatorname{R}\right\|_{\rm Frobenius}$, where $\widehat{\mathbf{x}}_{i}$, $\widehat{\mathbf{a}}$ and $\widehat{\operatorname{R}}$ are the 3D coordinates of the $i$th point, recovered shift vector and recovered rotation matrix after registration, respectively; $\mathbf{x}_{i}$, $\mathbf{a}$ and $\operatorname{R}$ are the ground-truth 3D coordinates of the $i$th point, ground-truth shift vector and ground-truth recovered rotation matrix. We see that high-pass graph filtering based resampled point cloud uses $20$-times fewer points and achieves even better results than using all the points. The shift and rotation errors of using high-pass graph filtering based resampling are significantly smaller than those of using all the points or using uniform resampling.

VI Conclusions

In this paper, we proposed a resampling framework to select a subset of points to extract application-dependent features and reduce the subsequent computation in a large-scale point cloud. We formulated an optimization problem to obtain the optimal resampling distribution, which is also guaranteed to be shift/rotation/scale invariant. We then specified the feature extraction operator to be a graph filter and studied the resampling strategies based on all-pass, low-pass and high-pass graph filtering. A surface reconstruction system based on graph filter banks was introduced to compress 3D point clouds. Three applications, including large-scale visualization, accurate registration and robust shape modeling, were presented to validate the effectiveness and efficiency of the proposed resampling methods. This work also pointed out many possible future directions of 3D point cloud processing, such as efficient 3D point cloud compression system based on graph filter banks, surface reconstruction based on arbitrary graphs, robust metric to evaluate the quality of a 3D point cloud.

Appendix A Appendices

A-A Proof of Lemma 1

Proof.

For the nonweighted version, we have

[TABLE]

where $\left(\Psi^{T}\Psi f(\operatorname{X})\right)_{i}\in\mathbb{R}^{K}$ is the $i$th row of $\Psi^{T}\Psi f(\operatorname{X})$, $f_{i}(\operatorname{X})\in\mathbb{R}^{K}$ is the $i$th row of $f(\operatorname{X})$, $\delta_{i}$ denotes an indicator function and equality (a) follows from the independent and identically distributed random sampling.

For the reweighted version, we have

[TABLE]

∎

A-B Proof of Theorem 1

Proof.

We first split the error into the bias term and the variance term,

[TABLE]

where the first term is bias and the second term is variance. Lemma (1) shows that the bias term is zero. So, we only need to bound the variance term.

For each element in the variance term, we have

[TABLE]

We finally combine all the elements and obtain (8). ∎

A-C Proof of Theorem 2

Proof.

Based on Theorem 1, we have

[TABLE]

∎

A-D Proof of Theorem 3

Proof.

The optimal resampling strategy is the solution of the following optimization problem,

[TABLE]

The corresponding Lagrange function is

[TABLE]

where the equality follows from Theorem 1. The derivative to $\pi_{i}$ is

[TABLE]

By setting its derivative to zero, we have

[TABLE]

Due to the complementary slackness, we have

[TABLE]

Thus, either $\mu_{i}$ or $\left\|f_{i}(\operatorname{X})\right\|_{2}$ is zero. In both cases, $\pi_{i}\propto\left\|f_{i}(\operatorname{X})\right\|_{2}$. ∎

A-E Proof of Theorem 4

Proof.

The optimal resampling strategy is the solution of the following optimization problem,

[TABLE]

The corresponding Lagrange function is

[TABLE]

where $\operatorname{F}_{i}$ is the $i$th row of $\operatorname{F}$ and $\left(\operatorname{F}\operatorname{\operatorname{X}_{\rm o}}\right)_{i}$ is the $i$th row of $\operatorname{F}\operatorname{\operatorname{X}_{\rm o}}$. The derivative to $\pi_{i}$ is

[TABLE]

By setting its derivative to zero, we have

[TABLE]

Due to the complementary slackness, we have

[TABLE]

Thus, either $\mu_{i}$ or $\left\|f_{i}(\operatorname{X})\right\|_{2}$ is zero. In both cases, $\pi_{i}\propto\sqrt{c^{2}\left\|\operatorname{F}_{i}\right\|_{2}^{2}+\left\|\left(\operatorname{F}\operatorname{\operatorname{X}_{\rm o}}\right)_{i}\right\|_{2}^{2}}$. ∎

A-F Proof of Theorem 5

Proof.

We first show the rotational invariance. Let $\operatorname{X}$ be the 3D coordinates of an original point cloud and $\operatorname{R}\in\mathbb{R}^{3\times 3}$ be a rotation matrix. The point cloud after rotating is $\operatorname{X}\operatorname{R}$. The local variation of $\operatorname{X}\operatorname{R}$ is

[TABLE]

where $\left(h(\operatorname{A})\right)_{i}$ is the $i$th row of $h(\operatorname{A})$ and (a) follows from any rotation matrix $\operatorname{R}$ is orthonormal.

We next show the shift variance. Let $\mathbf{a}\in\mathbb{R}^{3}$ be the shift and the point cloud after shifting is $\operatorname{X}+{\bf 1}\mathbf{a}^{T}$. The local variation of $\operatorname{X}+{\bf 1}\mathbf{a}^{T}$ is

[TABLE]

Thus, $f_{i}(\operatorname{X}+{\bf 1}\mathbf{a}^{T})=f_{i}(\operatorname{X})$ only when $h(\operatorname{A}){\bf 1}=\mathbb{0}\in\mathbb{R}^{N}$. ∎