---
abstract: 'Along with increasingly popular virtual reality applications, the three-dimensional (3D) point cloud has become a fundamental data structure to characterize 3D objects and surroundings. To process 3D point clouds efficiently, a suitable model for the underlying structure and outlier noises is always critical. In this work, we propose a hypergraph-based new point cloud model that is amenable to efficient analysis and processing. We introduce tensor-based methods to estimate hypergraph spectrum components and frequency coefficients of point clouds in both ideal and noisy settings. We establish an analytical connection between hypergraph frequencies and structural features. We further evaluate the efficacy of hypergraph spectrum estimation in two common point cloud applications of sampling and denoising for which also we elaborate specific hypergraph filter design and spectral properties. The empirical performance demonstrates the strength of hypergraph signal processing as a tool in 3D point clouds and the underlying properties.'
author:
- 'Songyang Zhang, Shuguang Cui, , and Zhi Ding, '
title: Hypergraph Spectral Analysis and Processing in 3D Point Cloud
---

3D point clouds, hypergraph signal processing, hypergraph construction, denoising, sampling.

Introduction {#intro}
============

Recent developments in depth sensors and softwares make it easier to capture the features and create a three-dimensional (3D) model for an object and its surroundings[@c1]. In particular, with the low-cost scanners such as light detection and ranging (LIDAR) and Kinect, a new data structure known as the point cloud has achieved significant success in many areas, including virtual reality, geographic information system, reconstruction of art document and high-precision 3D maps for self-driving cars [@c2]. A point cloud consists of 3D coordinates with attributes such as color, temperature, texture, and depth [@c3]. Owing to the easy access to scanning sensors and the huge need in describing the 3D features, the use of point clouds has attracted significant attentions in areas of computer vision, virtual reality, and medical science. How to process the point clouds efficiently becomes an important topic of research in many 3D imaging and vision systems.

To analyze the features of point cloud, the first step is to construct an analytical model to represent the 3D structures. The literature provides several different models. In [@c4], the 3D space is partitioned into several boxes or voxels, and the point clouds are then discretized therein. One disadvantage of voxels is that a dense grid is required to achieve fine resolution, leading to spatial inefficiency [@c3]. A spatially efficient approach [@c5; @c6] is the octree representation of point clouds. An octree is a tree data structure in which each node has exactly eight children. It can partition a 3D space recursively, and represent the point clouds with partitioned boxes. Although efficient, octree suffers from discretization errors [@c3]. The bd-tree is another spatial decomposition technique and is robust in highly cluttered point cloud dataset. However, compared to octree structures, bd-trees are more difficult to update.

Recently, graphs and graph signal processing (GSP) have found applications in modeling point clouds. For example, the authors of [@c3] construct a graph based on pairwise point distances. Some other works, such as [@c8], construct graphs based on the $k$-nearest neighbors, where each vertex (point) has an edge connection to its $k$ nearest neighbors. There are several clear connections between graph features and point cloud characteristics. For example, the smoothness over a graph can describe the flatness of surfaces in point clouds. GSP-based tools such as filters and graph learning methods can process the point clouds and have shown great success because of the graph model’s ability to capture the underlying geometric structures. However, graph-based methods still face some challenges such as limited orders and measurement inefficiency. In a traditional graph, each edge can only connect two nodes, constraining graph-based models to describe only pairwise relationships. However, a multilateral relationship among multiple nodes is far more informative as in a point cloud model. For example, the points (nodes) on the same surface of a point cloud exhibit a strong multilateral relationship, which cannot be easily captured by an edge of a traditional graph. In fact, construction of an efficient graph for a given dataset is always an open question. Thus, studies on point clouds can benefit from more general and efficient models.

To develop an efficient model for point clouds, we explore a high-dimensional graph model, known as hypergraph [@c9]. Hypergraph can be a useful model in processing 3D point clouds. A hypergraph $\mathcal{H}=\{\mathcal{V},\mathcal{E}\}$ consists of a set of nodes $\mathcal{V}=\{\mathbf{v}_1, \dots,\mathbf{v}_K\}$ and a set of hyperedges $\mathcal{E}=\{\mathbf{e}_1, \dots,\mathbf{e}_K\}$. Each hyperedge in a hypergraph can connect more than two nodes. For example, a 3D shape together with its hypergraph model are shown as Fig. \[hyper\]. Obviously, a normal graph is a special case of hypergraph, where each hyperedge degrades to connect two nodes exactly. The hyperedge in a hypergraph can characterize the multilateral relationship among several related nodes (e.g., on a surface), thereby making hypergraph a natural and intuitive model for point clouds. Moreover, advances in hypergraph signal processing (HGSP) [@c9] are providing more hypergraph tools, such as HGSP-based filters and spectrum analysis, for effective point cloud processing.

However, processing the point clouds based on hypergraph still poses several challenges. Similar to GSP, the first problem lies in the construction of hypergraph for point clouds. The traditional hypergraph construction method for a general dataset relies on data structure. For example, in [@c11], a hypergraph model is constructed according to the sentence structure in natural language processing. The $k$-nearest neighbor model is another method to construct the hypergraph. In [@c9], a hypergraph can be formed from the feature distances for an animal dataset to achieve clustering. However, such distance-based or structure-based model may be rather lossy in information preservation. For example, the structure-based method may not preserve the correlation of some irregular structures, whereas the $k$-nearest neighbor method may narrowly emphasize the distance information. In addition to hypergraph construction, another issue in analyzing point cloud with hypergraph tools is the computation complexity of the spectrum space. In the HGSP framework, spectrum-based analysis plays an important role but needs to compute the spectrum space. Usually, the computation of hypergraph spectrum is based on orthogonal-CP decomposition, which incurs high-complexity when there are many nodes. Another challenge in point cloud processing is the effect of noise and outliers. Since a hypergraph model is constructed from observed data, noise can distort the hypergraph and degrade the performances of HGSP. Thus, mitigating noise effect and robustly estimating the hypergraph model for point clouds pose a significant challenge.

This work addresses the aforementioned problems. We propose novel spectrum-based hypergraph construction methods for both clean and noisy point clouds. For clean point clouds, we first estimate their spectrum components based on the hypergraph stationary process and optimally determine their frequency coefficients based on smoothness to recover the original hypergraph structure. For noisy point clouds, we introduce a method for joint hypergraph structure estimation and data denoising. We shall illustrate the effectiveness of the proposed hypergraph construction and spectrum estimation in two point clould applications: sampling and denoising. Our experimental results clearly establish a connection between hypergraph frequencies and point cloud features. The performance improvement in both applications demonstrates the strength and power of hypergraph in point cloud processing and the practical value of our estimation methods.

We organize the rest of the paper as follows. In Section \[pre\], we lay the foundation with respect to the preliminaries and notations of point clouds, tensor basics and hypergraph signal processing. Next, we propose means in estimating hypergraph spectrum for basic point clouds in Section \[h1\] and further develop means for hypergraph structure estimation of noisy point clouds in Section \[h2\]. With the proposed estimation methods, we study two important application scenarios and establish the effectiveness of hypergraph signal processing in Section \[appli\]. Finally, we present the conclusion and future directions in Section \[con\].

Preliminaries and Notations {#pre}
===========================

In this section, we cover basic background with respect to point cloud, tensor basics and hypergraph signal processing.

Point Clouds
------------

A point cloud is a set of 3D points obtained from sensors, where each point is attributed with coordinates and other features, like colors [@c10]. Since the coordinates are basic features of a point cloud, in this work, we mainly focus on gray-scale point clouds, where each node is characterized by its coordinates. We consider a matrix representation of the gray-scale point clouds, where a point cloud with $N$ nodes is denoted by a location matrix $$\mathbf{s}=[\mathbf{X_1\quad X_2\quad X_3}]=
    \begin{bmatrix}
    \mathbf{s}_1^T\\
    \mathbf{s}_2^T\\
    \ddots\\
    \mathbf{s}_N^T
    \end{bmatrix}\in\mathbb{R}^{N\times 3},$$ where $\mathbf{X}_i$ denotes a vector of the $i$th coordinates of all the points, and $\mathbf{s}_i$ is the three coordinates of $i$th point. With the information of coordinates, different models, such as graphs [@c3] and octrees [@c5], can be constructed to analyze the point clouds, for which we will discuss more in Section \[appli\].

Tensor Basics
-------------

Tensor is a high-dimensional generalization of matrix. A tensor can be interpreted as multi-dimensional arrays. The order of tensor is the number of indices to label the components of arrays [@c17]. For example, a scalar is a zeroth-order tensor; a vector is a first-order tensor; a matrix is a second-order tensor; and an $M$-dimensional array is an $M$th-order tensor [@c18]. In this work, an $M$th-order tensor is denoted by $\mathbf{A}\in\mathbb{R}^{I_1
    \times I_2\times \cdots \times I_M}$, whose entry in position $(i_1,i_2,\cdots,i_M)$ is labeled as $a_{i_1\cdots i_M}$. Here, $I_k$ is the dimension of $k$th order.

Tensor outer product is a widely used operation to construct a higher-order tensor from lower-order tensors. The tensor outer product between an $M$th-order tensor $\mathbf{U}\in \mathbb{R}^{I_1\times I_2\times ...\times I_M }$ with entries $u_{i_1 ... i_M}$ and an $N$th-order tensor $\mathbf{V}\in \mathbb{R}^{J_1\times J_2\times ...\times J_N }$ with entries $v_{j_1 ... j_N}$ is denoted by $$\mathbf{W}=\mathbf{U} \circ \mathbf{V},$$ where the result $\mathbf{W}\in \mathbb{R}^{I_1\times I_2\times ...\times I_M \times J_1 \times J_2 \times ... \times J_N}$ is an $(M+N)$th-order tensor with entries $$w_{i_1 ... i_M j_1 ... j_N}= u_{i_1 ... i_M} \cdot v_{j_1 ... j_N}.$$

Hypergraph Signal Processing
----------------------------

Hypergraph signal processing (HGSP) is a tensor-based framework [@c9]. In the HGSP framework, a hypergraph with $N$ nodes and longest hyperedge connecting $M$ nodes, is represented by an $M$-th order $N$-dimension representing tensor $\mathbf{A}=(a_{i_1i_2\cdots i_M})\in\mathbb{R}^{N^M}$. The representing tensor can be adjacency tensor or Laplacian tensor in different purposes [@c16]. In this paper, we refer the adjacency tensor as the representing tensor, in which each entry $a_{i_1i_2\cdots i_M}$ indicates whether nodes $\{\mathbf{v}_1,\mathbf{v}_2,\cdots,\mathbf{v}_M\}$ are connected. The computation of the edge weight can be found in [@c9].

With the orthogonal-CP decomposition, the representing tensor can be decomposed via $$\label{decom}
\mathbf{A}=\sum_{r=1}^{N}\lambda_r\cdot\underbrace{\mathbf{f}_r\circ...\circ \mathbf{f}_r}_{\text{$M$ times}},$$ where $\mathbf{f}_r$’s are orthonormal basis called spectrum components and $\lambda_r$ are frequency coefficients related to the hypergraph frequency. All the spectrum components $\{\mathbf{f}_1,\cdots,\mathbf{f}_N\}$ construct the hypergraph spectral space. Each pair $(\mathbf{f}_r,\lambda_r)$ is called the spectral pair of the hypergraph.

Given an original signal $\mathbf{s}=[s_1\quad s_2\quad...\quad s_N]^{\mathrm{T}}$, the hypergraph signal is defined as the $(M-1)$ times tensor outer product of $\mathbf{s}$, i.e., $$\mathbf{s}^{[M-1]}=\underbrace{\mathbf{s\circ...\circ s}}_{\text{$M-1$ times}}.$$

The hypergraph frequency is ordered by the total variation of the spectrum component, which is defined as $$\mathbf{TV}(\mathbf{\mathbf{f}_r})=||\mathbf{f}_r-\frac{1}{\lambda_{max}}\mathbf{A}\mathbf{f}_r^{[M-1]}||_1,$$ where $\mathbf{A}\mathbf{f}_r^{[M-1]}$ is the contraction between representing tensor $\mathbf{A}$ and the hypergraph signal.

A spectrum component with larger total variation is a higher-frequency component, which indicates a faster propagation over the given hypergraph. Moreover, a supporting matrix $$\label{sup}
\mathbf{P_s}=\frac{1}{\lambda_{\max}}
\begin{bmatrix}
\mathbf{f}_1& \cdots& \mathbf{f}_N
\end{bmatrix}
\begin{bmatrix}
\lambda_1& & \\
&\ddots& \\
& &\lambda_N
\end{bmatrix}
\begin{bmatrix}
\mathbf{f}_1^{\mathrm{T}}\\
\vdots\\
\mathbf{f}_N^{\mathrm{T}}
\end{bmatrix},$$ can be defined to capture the overall spectral information of the hypergraph.

Instead of reviewing many properties of HGSP here, other aspects such as hypergraph Fourier transform, hypergraph filter design and sampling theory can be found in [@c9].

Hypergraph Spectrum Estimation for Point Clouds {#h1}
===============================================

To process the 3D point clouds, the first step is to construct an optimal hypergraph to model the point clouds. As we mentioned in the Section \[intro\], it is time-comsuming and inefficient to first construct a hypergraph structure before tensor decomposition to obtain the hypergraph spectrum. Instead, we propose to directly estimate the hypergraph spectral pairs based on the observed data, and then recover the original representing tensor with Eq. (\[decom\]). In this section, we first estimate the hypergraph spectrum components $\mathbf{f}_r$’s based on the hypergraph stationary process, and optimize the frequency coefficients $\lambda_r$’s based on the smoothness for original point clouds.

Estimation of Hypergraph Spectrum Components
--------------------------------------------

In this part, we propose a method to estimate the hypergraph spectral components based on the hypergraph stationary process.

### Hypergraph Stationary Process

Before providing details of the estimation, let us first introduce some new definitions and properties necessary for spectrum estimation.

Stationarity is a cornerstone property that facilities the analysis of random signals and observations in traditional signal processing [@c19]. It has equal importance in graph and hypergraph signal processing. Based on graph shifting introduced in [@c22], a definition of graph stationary process proposed in [@c19] can analyze the properties of the different observations of nodes, or the random signals over the graphs. Furthermore, [@c123] introduces a method to estimate the graph spectrum space and graph diffusion for multiple observations based on the graph stationary process. Similarly, the hypergraph stationary process can be defined to estimate hypergraph spectrum.

Now, let us introduce the definition of the hypergraph stationary process. In [@c9], a polynomial hypergraph filter based on supporting matrix is defined as $$\mathbf{s}'=\sum_{k=1}^{a}\alpha_k\mathbf{P}^k\mathbf{s},$$ where $\mathbf{P}=\lambda_{max}\mathbf{P_s}$.

Similarly, based on the supporting matrix, a $\tau$-step shifting operation is defined as $\mathbf{P}_{\tau}=\mathbf{P}^{\tau}$. Then, similar to the definition of the stationary process in traditional digital signal processing and graph signal processing, a strict-sense stationary process in HGSP can be defined as follows.

(Strict-Sense Stationary Process) A stochastic signal $\mathbf{x}\in\mathbb{R}^N$ is strict-sense stationary over the hypergraph with $\mathbf{P}_\tau$ if and only if $$\mathbf{x}\overset{d}{=}\mathbf{P}_{\tau}\mathbf{x}$$ holds for any $\tau$.

Since the strict-sense stationary is hard to achieve and analyze in the real datasets, we introduce the weak-sense stationary process similar to traditional digital signal processing.

(Weak-Sense Stationary Process) A stochastic signal $\mathbf{x}\in\mathbb{R}^N$ is weak-sense stationary over the hypergraph with $\mathbf{P}_\tau$ if and only if $$\label{mean}
    \mathbb{E}[\mathbf{x}]=\mathbb{E}[\mathbf{P}_{\tau}\mathbf{x}]$$ and $$\label{time}
    \mathbb{E}[(\mathbf{P}_{\tau_1}\mathbf{x})((\mathbf{P}^H)_{\tau_2}\mathbf{x})^H]=\mathbb{E}[(\mathbf{P}_{\tau_1+\tau}\mathbf{x})((\mathbf{P}^H)_{\tau_2-\tau}\mathbf{x})^H]$$ hold for any $\tau$, where $\mathbb{E}(\cdot)$ refers to the mean of observations and $(\cdot)^H$ is the Hermitian transpose.

From the definition of the weak-sense stationary process (WSS), Eq. (\[mean\]) implies that the mean function of the signal must be constant, which is the same condition as in traditional digital signal processing (DSP) [@c23]. From the definition of supporting matrix, the $(i,j)$-th entry of $\mathbf{P}$ is the same as the $(j,i)$-th entry of $\mathbf{P}^H$, which indicates that $\mathbf{P}^H$ is the shifting in the opposite direction of $\mathbf{P}$. Then, the condition in Eq. (\[time\]) indicates that the hypergraph covariance function $K_{\mathbf{xx}}({\tau_1},-\tau_2)=K_{\mathbf{xx}}({\tau_1}+\tau,\tau-\tau_2)=K_{\mathbf{xx}}({\tau_1}+\tau_2,0)$, which is also consistent with the definition in traditional DSP.

With the definition of the hypergraph stationary process, we have the following properties regarding the relationship between signals and hypergraph spectrum.

A stochastic signal $\mathbf{x}$ is WSS if and only if it has zero-mean and its covariance matrix has the same eigenvectors as the hypergraph spectrum basis, i.e., $$\label{s1}
    \mathbb{E}[\mathbf{x}]=\mathbf{0}$$ and $$\label{s2}
    \mathbb{E}[\mathbf{x}\mathbf{x}^H]=\mathbf{V}\Sigma_\mathbf{x}\mathbf{V}^{H},$$ where $\mathbf{V}=[\mathbf{f}_1,\mathbf{f}_2,\cdots,\mathbf{f}_N]\in\mathbb{R}^{N\times N}$ are the hypergraph spectrum.

Since the hypergraph spectrum basis are orthonormal, we have $\mathbf{VV}^T=\mathbf{I}$. Then, the $\tau$-step shifting based on supporting matrix can be calculated as $$\begin{aligned}
    \mathbf{P}_\tau&=\underbrace{\mathbf{V}\Lambda_\mathbf{P}\mathbf{V}^T\mathbf{V}\Lambda_\mathbf{P}\mathbf{V}^T\cdots\mathbf{V}\Lambda_\mathbf{P}\mathbf{V}^T}_{\tau\quad times}\\
    &=\mathbf{V}\Lambda_\mathbf{P}^\tau\mathbf{V}^T.\label{poly}
    \end{aligned}$$

Now, the Eq. (\[mean\]) can be written as $$\mathbb{E}[\mathbf{x}]=\mathbf{V}\Lambda_P^\tau\mathbf{V}^T\mathbb{E}[\mathbf{x}].$$ Since $\mathbf{V}\Lambda_P^\tau\mathbf{V}^T$ does not always equal to $\mathbf{I}$, Eq. (\[mean\]) holds for arbitrary supporting matrix and $\tau$ if and only if $\mathbb{E}[\mathbf{x}]=\mathbf{0}$.

Next we show the sufficiency and necessity of the condition in Eq. (\[s2\]). The condition in Eq. (\[time\]) can be written as $$\mathbf{P}_{\tau_1}\mathbb{E}[\mathbf{xx}^H]((\mathbf{P})^H_{\tau_2})^H=\mathbf{P}_{\tau_1+\tau}\mathbb{E}[\mathbf{xx}^H]((\mathbf{P})^H_{\tau_2-\tau})^H.$$ Considering Eq. (\[poly\]) and the fact that hypergraph spectrum is real [@c9], Eq. (\[time\]) is equivalent to $$\mathbf{V}\Lambda_\mathbf{P}^{\tau_1}\mathbf{V}^H\mathbb{E}[\mathbf{xx}^H]\mathbf{V}\Lambda_\mathbf{P}^{\tau_2}\mathbf{V}^H=\mathbf{V}\Lambda_\mathbf{P}^{\tau_1+\tau}\mathbf{V}^H\mathbb{E}[\mathbf{xx}^H]\mathbf{V}\Lambda_\mathbf{P}^{\tau_2-\tau} \mathbf{V}^H,$$ which can be written as $$\label{eig}
    (\mathbf{V}^H\mathbb{E}[\mathbf{xx}^H]\mathbf{V})\Lambda_\mathbf{P}^{\tau}=\Lambda_\mathbf{P}^{\tau}(\mathbf{V}^H\mathbb{E}[\mathbf{xx}^H]\mathbf{V}).$$ If Eq. (\[eig\]) holds for arbitrary $\mathbf{P}$, $(\mathbf{V}^H\mathbb{E}[\mathbf{xx}^H]\mathbf{V})$ should be diagonal, which indicates $\mathbb{E}[\mathbf{x}\mathbf{x}^H]=\mathbf{V}\Sigma_\mathbf{x}\mathbf{V}^{H}$. Thus, the sufficiency of the condition is proved.

Similarly, we can apply Eq. (\[s2\]) on both sides of Eq. ($\ref{time}$), we can establish the necessity of the condition in Eq. (\[s2\]).

This theorem can be used to estimate the hypergraph spectrum, given multiple observations of several signal points.

### Estimation of Spectrum Components for Point Clouds

Now, we can use the property of stationary process to estimate the hypergraph spectrum of point clouds. The three coordinates of a point can be interpreted as three observations of the point from different angles, which describe the underlying multilateral relationship. Thus, we can assume that the point cloud signals follow the stationary process over the estimated underlying hypergraph structure. If the point cloud signals $\mathbf{s}$ follow the hypergraph stationarity, it should satisfy Eq. (\[s1\]) and Eq. (\[s2\]). Thus, a spectrum estimation method can be based on hypergraph staionarity. The details of the algorithm is described as follows.

: Point cloud dataset $\mathbf{s}=[\mathbf{X_1\quad X_2\quad X_3}]\in\mathbb{R}^{N\times 3}$. Calculate the mean of each row in $\mathbf{s}$, i.e.,\
$\mathbf{\overline s}=(\mathbf{X_1+X_2+X_3})/3$; Normalize the original point cloud data as zero-mean in each row, i.e., $\mathbf{s}'=[\mathbf{X_1-{\overline s},X_2-{\overline s},X_3-{\overline s}}]$; Calculate the eigenvectors $\{\mathbf{f}_1,\cdots,\mathbf{f}_N\}$ for $R_{\mathbf{s}'}=\mathbf{s'}(\mathbf{s'}^T)$; : Hypergraph spectrum $\mathbf{V}=[\mathbf{f}_1,\cdots,\mathbf{f}_N]$.

With *Theorem 1*, we can directly obtain an estimation of the hypergraph spectrum based on the hypergraph stationarity. Note that, here, we assume all the observations are from a clean point cloud without noise. The case of noisy point clouds will be discussed later in Section \[h2\].

Estimation of Frequency Coefficients {#es}
------------------------------------

Next, we discuss how we estimate the hypergraph frequency coefficients with the spectrum components based on the hypergraph smoothness.

In real applications, the large-scale networks are usually sparse, which makes it meaningful to infer that most entries of the hypergraph representing tensor for real datasets are zero [@c24]. In addition, the smoothness of signals is a widely-used assumption when estimating the underlying structure of graphs and hypergraphs [@c25]. Thus, the estimation of the hypergraph representing tensor with known spectrum components for a given dataset $\mathbf{s}$ can be generally formulated as $$\begin{aligned}
&\min_{\mathbf{\boldsymbol{\lambda}}}\quad \alpha \mbox{Smooth}
(\mathbf{s,\boldsymbol{\lambda}},\mathbf{f}_r)+\beta||\mathbf{A}||_T^2\label{e1}\\
s.t.\quad &\mathbf{A}=\sum_{r=1}^{N}\lambda_r\cdot\underbrace{\mathbf{f}_r\circ...\circ \mathbf{f}_r}_{\text{M times}}. \label{dec}\\
&\quad\mathbf{A}\in \mathcal{A}.\label{cs1}\\
&||\mathbf{A}||_T=\sqrt{\sum_{i_1,i_2,\cdots, i_M=1}^N a_{i_1i_2\cdots i_M}^2}.\label{t_norm}\end{aligned}$$

The constraint set $\mathcal{A}$ in (\[cs1\]) includes the prior information of the representing tensor. For example, if the representing tensor is the adjacency tensor, its entries should be non-negative. In the constraint of $(\ref{t_norm})$, $||\mathbf{A}||_T$ is the tensor norm which controls the sparsity of the hypergraph structure. The smoothness function Smooth$(\mathbf{s,\boldsymbol{\lambda}},\mathbf{f}_r)$ can be designed for specific problems. Typical functions can be hypergraph Laplacian regularization, label ranking, and total variation [@c9]. For convenience, we use the quadratic-form total variation based on the supporting matrix to describe the hypergraph smoothness, i.e., $$\mathbf{TV}(\mathbf{s})=||\mathbf{s}-(1/\lambda_{max})\mathbf{P}\mathbf{s}||^2_2.$$ This form of smoothness function suggested in [@c9] can capture the differences between one node and its neighbors over hypergraph. Since the signals are smooth over the estimated hypergraph, observations are also smooth. Thus, the final smoothness function for point cloud $\mathbf{s}=[\mathbf{X_1\quad X_2\quad X_3}]$ is $$\begin{aligned}
\mbox{Smooth}(\mathbf{s,\boldsymbol{\lambda}},\mathbf{f}_r)&=\sum_{i=1}^3||\mathbf{X}_i-\mathbf{P_sX}_i||^2_2\nonumber\\
&=\sum_{i=1}^3||\mathbf{X}_i-\sum_r \sigma_r(\mathbf{f}_r^T\mathbf{X}_i)\mathbf{f}_r||^2_2\nonumber\\
&=\sum_{i=1}^3||\mathbf{X}_i-\mathbf{W}_i\boldsymbol{\sigma}||^2_2,\label{sm1}\end{aligned}$$ where $\mathbf{W}_i=[(\mathbf{f}_1^T\mathbf{X}_i)\mathbf{f}_1\quad(\mathbf{f}_2^T\mathbf{X}_i)\mathbf{f}_2\quad\cdots\quad(\mathbf{f}_N^T\mathbf{X}_i)\mathbf{f}_N]$, $\sigma_r=\lambda_r/\lambda_{max}$ and $\boldsymbol{\sigma}=[\sigma_1\cdots\sigma_N]$.

Moreover, the tensor norm of a given hypergraph has the following property with the frequency coefficients.

Given a representing tensor $\mathbf{A}=\sum_{r=1}^{N}\lambda_r\cdot\underbrace{\mathbf{f}_r\circ...\circ \mathbf{f}_r}_{\text{M times}}$, the tensor norm $||\mathbf{A}||^2_T=\sum_{i_1,i_2,\cdots,i_M=1}^{N}a_{i_1i_2\cdots i_M}^2$ can be written in the form of frequency coefficients as $$||\mathbf{A}||^2_T=\sum_{r=1}^N \lambda_r^2=\boldsymbol{\lambda}^T\boldsymbol{\lambda},$$ where $\boldsymbol\lambda=[\lambda_1 \quad\lambda_2 \quad\cdots\quad\lambda_N]$.

Since $\mathbf{A}=\sum_{r=1}^{N}\lambda_r\cdot\underbrace{\mathbf{f}_r\circ...\circ \mathbf{f}_r}_{\text{M times}}$, we have $$a_{i_1i_2\cdots i_M}=\sum_{r=1}^N \lambda_rf_{r,i_1}f_{r,i_2}\cdots f_{r,i_M},$$ where $f_{r,i}$ is the $i$th element of $\mathbf{f}_r$. Then, the tensor norm is $$\begin{aligned}
    ||\mathbf{A}||^2_T
    %&=\sum_{i_1,i_2,\cdots,i_M}a_{i_1i_2\cdots i_M}^2\nonumber\\
    &=\sum_{i_1,i_2,\cdots,i_M}(\sum_{r=1}^N \lambda_rf_{r,i_1}f_{r,i_2}\cdots f_{r,i_M})^2\nonumber\\
    &=\sum_{i_1,i_2,\cdots,i_M}(\sum_{r=1}^N \lambda_rf_{r,i_1}\cdots f_{r,i_M})(\sum_{t=1}^N \lambda_tf_{t,i_1}\cdots f_{t,i_M})\nonumber\\
    &=\sum_{i_1,i_2,\cdots,i_M}\sum_{r,t}\lambda_r\lambda_tf_{r,i_1}\cdots f_{r,i_M} f_{t,i_1}\cdots f_{t,i_M}\nonumber\\   
    &=\sum_{r,t}\lambda_r\lambda_t\sum_{i_1,i_2,\cdots,i_M=1}^N(f_{r,i_1}f_{t,i_1})\cdots(f_{r,i_M}f_{t,i_M})\nonumber\\
%   &=\sum_{r,t}\lambda_r\lambda_t(\sum_{i=1}^Nf_{r,i}f_{t,i})^M\nonumber\\
    &=\sum_{r,t}\lambda_r\lambda_t(\mathbf{f}_r^T\mathbf{f}_t)^M.
    \end{aligned}$$ Since $\mathbf{f}_r$ is orthogonal, $\mathbf{f}_r^T\mathbf{f}_t=1$ holds if $r=t$; otherwise, $\mathbf{f}_r^T\mathbf{f}_t=0$. Thus, we obtain $||\mathbf{A}||^2_T=\sum_{r=1}^N \lambda_r^2$.

This property can help us build a connection from the tensor norm to the frequency coefficients directly.

Now, if we consider the representing tensor as the adjacency tensor and each hyperedge consists of three nodes since at least three nodes are required to construct a surface, we optimize the normalized frequency coefficients $\boldsymbol{\sigma}=\frac{1}{\lambda_{max}}\boldsymbol{\lambda}=[\sigma_1\quad 
\sigma_2 \quad \cdots\quad \sigma_N]^T$ via $$\begin{aligned}
\label{target}
& \min_{\boldsymbol{\sigma}}\quad \alpha\sum_{i=1}^3||\mathbf{X}_i-\mathbf{W}_i\boldsymbol{\sigma}||^2_2+\beta{ \boldsymbol\sigma^T\boldsymbol\sigma}\\
\hspace{-3mm}s.\; t. \; &\;\;0\leq \sigma_r\leq \max_i \sigma_{i}=1,\label{non}\\
&\sum_{r=1}^N \sigma_r f_{r,i_1}f_{r,i_2}f_{r,i_3}\geq 0, \quad i_1,i_2,i_3=1,2,\cdots,N.\label{adj}\end{aligned}$$

The constraint (\[adj\]) limits the estimated representing tensor as the adjacency tensor. The constraint (\[non\]) is the nonnegative constraint on weight and the factor [@c26]. Clearly, the optimization is non-convex with the constraint $\max_i \sigma_{i}=1$. However, if the position of the maximal frequency is known, the optimization problem can be solved by tools such as cvx [@c27; @c28]. Thus, we can develop the following algorithm to estimate the frequency coefficients.

: Point cloud dataset $\mathbf{s}=[\mathbf{X_1,X_2,X_3}]\in\mathbb{R}^{N\times 3}$, hypergraph spectrum $\mathbf{V}=[\mathbf{f}_1,\cdots,\mathbf{f}_N]$. i=1,2,...,iter [**[do]{}**]{}: Set $\sigma_i=1$ as the maximal normalized eigenvalue. Solve the optimization problem in Eq. (\[target\]). Find the optimal $i$ to minimize the target function. The optimal coefficients is the solution of Eq. (\[target\]) correlated to the optimal $i$. : Frequency coefficients [$\boldsymbol\sigma$]{}.

Note that, since we consider clean point cloud without noise, we usually set parameter $\alpha\ll\beta$. Then, from the estimated spectrum pair $(\mathbf{f}_r,\sigma_r)$ under normalization, we can recover the original adjacency tensor as Eq. (\[dec\]). Hence, the hypergraph construction process for a clean point cloud can be summarized as Fig. \[fram1\]. The recovery of original adjacency tensor is not always necessary in practical applications since storing the representing tensor is less efficient than storing the spectrum pairs.

![Estimation of Hypergraph Spectral Pairs for Original Point Clouds[]{data-label="fram1"}](Framework1.jpg){width="3in"}

Joint Spectrum Estimation and Denoising {#h2}
=======================================

In practical 3D imaging, perturbations such as noises and outliers often exist when generate a point cloud of an unknown object. These noises may significantly affect the performance of point cloud processing since many existing algorithms require quality datasets [@c27]. Thus, denoising remains a vital issue in practical point cloud applications.

Usually, to denoise point sets with sharp features is difficult, especially when the noise is large, as such features are hard to distinguish from noise effect. Generally, smoothness-based methods are common. In [@c30], a method based on $L_0$ norm of differences between $k$-nearest neighbors is introduced. In [@c31], Laplacian regularization is used to describe smoothness and to denoise noisy point sets. Other works, such as [@c8; @c29], minimize the total variation over graphs to denoise the point sets. Although smoothness-based methods have achieved notable successes, how to interpret and define an effective smoothness function for a general point set remains open. Furthermore, for graph-based smoothness methods, the construction of graph model remains a critical problem, since traditional methods based on distance suffers from the imprecise location measurement. To this end, a more general definition of smoothness and a more efficient denoising method for arbitrary point clouds are highly desirable.

In this section, we introduce a joint method to simultaneously estimate the hypergraph structure and denoise noisy point clouds. In Section \[h1\], we already introduce an estimation method of spectral pair $(\mathbf{f}_r,\sigma_r)$ for clean point clouds. A similar construction process can be developed for the noisy point clouds. As the estimation of spectrum components only depends on the observed data, we need to denoise the noisy observations while optimizing the frequency coefficients. As already discussed, the problem of denoising a signal on a hypergraph can be written as a convex minimization problem with the constraints that denoised signals should be smooth over the hypergraph. Accordingly, the general process of hypergraph denoising and estimation can be summarized as the following steps:

-   Step 1: Estimate the approximated hypergraph spectrum components from the observed noisy point clouds;

-   Step 2: Jointly estimate frequency coefficients and denoise the noisy observations;

-   Step 3: Update the noisy observations as denoised data and repeat Step 1 until enough iterations.

To estimate hypergraph spectral components of noisy data, the process is the same as Algorithm 1 based on hypergraph stationary process. To jointly estimate the frequency coefficients to recover the original underlying structure and to denoise the noisy point clouds, we propose the following objective. Given $N$ noisy points $\mathbf{s}=[\mathbf{X_1\quad X_2\quad X_3}]$, the joint estimation task can be formulated as

$$\begin{aligned}
    \min_{\boldsymbol{\sigma},\mathbf{Y}}\quad& \sum_{i=1}^3[||\mathbf{X}_i-\mathbf{Y}_i||^2_2+\alpha ||\mathbf{X}_i-\mathbf{W}_i\boldsymbol{\sigma}||^2_2]+\beta ||\mathbf{A}||_2^2\label{esfram}\\
    s.t.\quad &\mathbf{A}=\sum_{r=1}^{N}\lambda_r\cdot\underbrace{\mathbf{f}_r\circ...\circ \mathbf{f}_r}_{\text{M times}} \in \mathcal{A},\nonumber\\
 &\;\;0\leq \sigma_r\leq \max_i \sigma_{i}=1, \nonumber\\
    &\mathbf{W}_i=[(\mathbf{f}_1^T\mathbf{X}_i)\mathbf{f}_1\quad(\mathbf{f}_2^T\mathbf{X}_i)\mathbf{f}_2\quad\cdots\quad(\mathbf{f}_N^T\mathbf{X}_i)\mathbf{f}_N].
    \nonumber\end{aligned}$$

The resulting $\mathbf{Y=[Y_1\quad Y_2 \quad Y_3]}$ is the denoised point clouds, and $(\alpha, \beta)$ are two positive regularization parameters. The first part in Eq. (\[esfram\]) lets the denoised point cloud maintain the observed structural features. The second part is the smoothness function derived from Eq. (\[sm1\]) which adjusts positions of noisy points. The third part is the tensor norm regularization to control hypergraph sparsity.

The optimization problem of Eq. (\[esfram\]) is not convex in $\mathbf{Y}$ and $\boldsymbol{\sigma}$. Therefore, similar to [@c25], we split the problem into two subproblems. For each subproblem, we fix one variable set to solve the other one. Upon convergence, the solution corresponds to a local minimum and not necessarily a global minimum.

We first initialize $\mathbf{Y}$ as the observed signals $\mathbf{X}$ and solve the following problem similar to that in Section \[h1\]. $$\label{target1}
\min_{\boldsymbol{\sigma}} \alpha\sum_{i=1}^3||\mathbf{X}_i-\mathbf{W}_i\boldsymbol{\sigma}||^2_2+\beta{ \boldsymbol\sigma^T\boldsymbol\sigma}$$ $$\begin{aligned}
s.t.\quad
&\;\;0\leq \sigma_r\leq \max_i \sigma_{i}=1, \nonumber\\
&\sum_{r=1}^N \sigma_r f_{r,i_1}f_{r,i_2}f_{r,i_3}\geq 0, \quad i_1,i_2,i_3=1,2,\cdots,N.\nonumber\end{aligned}$$ This problem can be solved similarly to the solution of clean point cloud with Algorithm 2.

Once the estimated frequency coefficients are found, we solve the subproblem of point cloud denoising $$\begin{aligned}
\label{ss2}
    \min_\mathbf{Y} \sum_{i=1}^3[||\mathbf{X}_i-\mathbf{Y}_i||^2_2+\alpha ||\mathbf{X}_i-\mathbf{W}_i\boldsymbol{\sigma}||^2_2],\end{aligned}$$ whose close-form solution for each coordinate is $$\label{ss3}
    \mathbf{Y}_i=[\mathbf{I}+\alpha(\mathbf{I-P_s})^T(\mathbf{I-P_s})]^{-1}\mathbf{X}_i.$$

Note that $\mathbf{P_s}$ is the supporting matrix. We then update the frequency components based on the denoised point clouds, and repeatedly carry out Step 1 to Step 3 until getting the final solution. In practice, we generally observe the convergence within only a few iterations. The complete algorithm is summarized in Algorithm 3 as shown in Fig. \[fram2\]. Unlike for clear point clouds, we emphasize more on the smoothness of signals over the hypergraph. The parameter $\alpha$ can be set larger than used when dealing with clean point clouds.

: Noisy observations of point clouds $\mathbf{s}=[\mathbf{X_1,X_2,X_3}]\in\mathbb{R}^{N\times 3}$. : Calculate the spectrum components $\mathbf{f}_r$’s from the observed point cloud $\mathbf{s}$ as Algorithm 1. i=1,2,...,iter [**[do]{}**]{}: Find the optimal $\boldsymbol{\sigma}$ for the first subproblem in Eq. (\[target1\]) with Algorithm 2. Solve the optimization problem in Eq. (\[ss2\]) with $\mathbf{Y}$ in Eq. (\[ss3\]). Update the observed signals as $\mathbf{Y}$ and recalculate the spectrum components $\mathbf{f}_r$’s. : Spectral pairs $(\mathbf{f}_r,\sigma_r)$’s, denoised point clouds $\mathbf{Y}$.

![Joint Hypergraph Estimation and Denoising for Noisy Point Cloud.[]{data-label="fram2"}](Framework2.jpg){width="3in"}

Application Examples {#appli}
====================

In this section, we examine two application examples to test the efficacy of the proposed method in estimating hypergraph structure for both clear and noisy point clouds.

Sampling
--------

Sampling is an important operation to facilitate analysis of very large point clouds. In this part, we consider different sampling strategies depending on different kinds of applications. Some interesting connections are found from the hypergraph frequency and point cloud features.

### Resampling using Harr-like Highpass Filtering

Filtering helps extract select features of a given dataset. In some applications such as boundary detection, accurate extraction of shape features of point clouds is important. Thus, an efficient sampling should retain the features of the original point cloud. In our estimation of hypergraph structure, smoothness is a significant feature to model point clouds. Ideally, smoothness over the original surface of a point cloud should correspond to smoothness over its hypergraph model. Therefore, we can also design a Harr-like high-pass filter to extract sharp features over the surfaces.

Let $\mathbf{I}$ be an identity matrix of appropriate size. Similar to that in GSP [@c3], a Haar-like high-pass filter is designed as $$\begin{aligned}
\mathbf{H}
&=\mathbf{I-P_s}\\
&=\mathbf{V}
\begin{bmatrix}
1-\sigma_1&0&\cdots&0\\
0&1-\sigma_2&\cdots&0\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\cdots&1-\sigma_N
\end{bmatrix}\mathbf{V}^T.\end{aligned}$$ The filtered signal is $$(\mathbf{Hs})_i=\mathbf{s}_i-\sum_j {P_s}_{(ij)}\mathbf{s}_j,$$ which reflects the differences between nodes and their neighbors over the hypergraph. Note that, the frequency coefficients together with their corresponding spectral components are ordered decreasingly here, i.e., $\sigma_{i}\geq\sigma_{i+1}$. From the definition of total variation, more smoothness corresponds to larger total variation. Thus, we can extract the sharp features over the point clouds by sampling the nodes with large value of $||\mathbf{s}_i-\sum_j {P_s}_{(ij)}\mathbf{s}_j||^2_2$.

To test this application, we estimate the spectral pairs for clean point clouds and filter the signals over several synthetic datasets. We randomly generate multiple points over the surfaces of basic graphics shown as Fig. \[ori\], and sample the point clouds using the high-pass filter (HPF) given in Fig. \[sam\]. From the test results, we can see that the sampled points of the surfaces in Fig. \[sur2\] mainly congregate near the corners and edges, which are the sharp parts of the point clouds. In addition, the sampled nodes for a cube shape are also crowded near edges and corners. On the other hand, the sampled nodes of a cylinder are mostly at the boundaries of the cylinder. Our test results show that the Harr-like HPF can extract sharp features from point cloud surfaces, which correspond to the least smooth parts of the estimated hypergraph. Moreover, since the total variation measures the order of frequency, sharp features over the point cloud correspond to high frequency components. Thus, the hypergraph model and the estimated spectral pairs are efficient when extracting features of 3D point clouds.

### Down-Sampling with Hypergraph Fourier Transform

Projecting signals into a suitable orthonormal basis is a widely-used sampling method [@c32]. The work of [@c9] develops a sampling theory based on hypergraph signal processing as follows:

-   Step 1: Order the spectrum components from low frequency to high frequency based on their total variations.

-   Step 2: Implement hypergraph Fourier tranform as $$\mathcal{F}(\mathbf{s})=[(\mathbf{f}_1^T\mathbf{s})^{M-1}\quad (\mathbf{f}_2^T\mathbf{s})^{M-1}\quad \cdots\quad (\mathbf{f}_N^T\mathbf{s})^{M-1}]^T.$$

-   Step 3: Use $C$ transformed signal components in the hypergraph frequency domain to represent $N$ signals in the original vertex domain.

More specifically, for a $K$-bandlimitted hypergraph signal, a perfect recovery is available with $K$ samples in hypergraph frequency domain. Similarly, we can sample the point clouds based on the hypergraph Fourier transform. To test the performance of the sampled signals, we implement hypergraph Fourier transform (HGFT) on each coordinates of the point clouds, i.e., $\mathcal{F}(\mathbf{X}_i)$ for all $i$. Then, we take the first $C$ transformed signals in all coordinates. Finally, we implement the inverse hypergraph Fourier transform (iHGFT) to obtain the sampled shapes of the original point clouds. Note that, perfect recovery happens with $C$ samples, if $(\mathcal{F}(\mathbf{X}_i))_{j+C}=0$ for $i,j\in {\cal Z}^+$.

We test the recovered point clouds for animal point datasets [@c33; @c34; @c35; @c36] with the GSP-based methods. For the GSP-based method, we construct the a graph adjacency matrix $\mathbf{W}$ with Guassian model, i.e., $$\label{adjj}
W_{ij}=\left\{
\begin{aligned}
\exp\left(-\frac{||\mathbf{s}_i-\mathbf{s}_j||^2_2}{\delta^2}\right),&\quad||\mathbf{s}_i-\mathbf{s}_j||^2_2\leq t;\\
0,&\quad \mbox{otherwise},
\end{aligned}
\right.$$ where $\mathbf{s}_i$ is the coordinates of the $i$th node. Then, we sample the point clouds using the signals after the graph Fourier transform (GFT).

The test point cloud is shown as Fig. \[sam2\]. We first compare the mean squared error (MSE) between the recovered point clouds and original point clouds shown as Fig. \[sam3\]. From the experimental results, we can see that the HGSP-based method has smaller error than the GSP downsampling method, clearly indicating hypergraph to be a better model. However, sometimes, MSE alone cannot tell the true story in terms of the performance for the recovered point clouds. To explore more, we compare the recovered point clouds directly in Fig. \[sam4\]. From the experimental results, we can see that HGSP-based method captures the overall structure of the point clouds with very few samples, whereas the GSP-based method requires more samples to get sufficient details. The MSE of GSP mainly stems from some outliers when taking more than 90 percent of the samples. The experiments show that HGSP-based method is a better tool for applications which need to recover an overall shape of point clouds from limited data storage. Our test shows hypergraph to be a suitable model for point clouds, and the estimated hypergraph spectral pairs capture the point cloud characteristics very well.

Denoising
---------

From estimated hypergraph spectral pairs from noisy point clouds, the performance of denoising is an intuitive metric of how good the estimates are. There are multiple methods developed to denoise noisy point clouds. The authors of [@c8] proposed a graph-based method to denoise based on total variation (GSP-TV). This method constructs a graph based on observed coordinates first before solving the denoising optimization $$\min_{\mathbf{Y}} ||\mathbf{X}-\mathbf{Y}||_2^2+\alpha \mathbf{TV}(\mathbf{Y,W}),$$ where $\mathbf{X}$ is the observed coordinates, and $\mathbf{W}$ is the adjacency matrix. Here, the graph total variation $\mathbf{TV}(\mathbf{Y,W})$ is applied in describing the smoothness over the graphs. In addition to total variation, Laplacian regularization (LR) has also been used in denoising with a basic formulation $$\min_{\mathbf{Y}} ||\mathbf{X}-\mathbf{Y}||_2^2+\alpha ||\mathbf{Y}^T\mathbf{LY}||^2_2,$$ where $\mathbf{L}$ is the Laplacian matrix. Developed from traditional Laplacian regularization methods, a mesh Laplacian smooth (MLS) method is given in [@c37].

                                  HGSP      GSP(TV)    MLS       LR     Noisy
  ---------------------------- ----------- --------- -------- -------- --------
   Uniform$\sim$U(-0.03,0.03)   **32.60**    45.94    56.63    48.86    63.84
   Uniform$\sim$U(0.08,0.16)    **98.36**   160.18    205.15   168.17   220.96
    Guassian$\sim$N(0,0.08)     **41.10**    42.36    49.41    64.00    76.54
   Guassian$\sim$N(0.02,0.08)   **73.43**    76.07    83.25    123.08   142.11
        Impulse (p=0.08)        **34.53**    45.45    50.89    40.53     60.5

  : Error in Dfferent Kinds of Noise[]{data-label="t1"}

\

To validate the performance of our denoising method, we compare with the aforementioned traditional methods using the Standford bunny dataset with 3595 points and sampled bunny with 397 points shown as Fig. \[den\]. We compare different methods in the sampled bunny dataset adding zero-mean Guassian noise with variance $\sigma^2$, and zero-mean Uniform noise with the interval $B-A$, respectively. We use the error denoted by $$Error=\sum_{i=1}^N\sum_{j=1}^3 |X_{ji}-Y_{ji}|,$$ where $X_{ij}$ and $Y_{ji}$ are the $j$th coordinates of observed and denoised point $i$, respectively, to measure the performance. We repeat the test on 1000 randomly generated noisy data. The error between the original dataset and the denoised dataset is shown in Fig. \[den1\]. The error of the noisy point clouds before denoising is also given as a reference in Fig. \[den1\]. From the test results, we can see that the HGSP-based method can achieve the lowest error, which demonstrates the effectiveness of the proposed denoising methods and estimated spectral pairs. The comparison in other types of noise is shown in Table. \[t1\]. More specifically, the methods based on total variation, i.e, HGSP and GSP-TV, have better performance than the methods based on Laplacian regularization, which indicates the total variation has a more efficient representation of the surface smoothness. The denoised bunny with 3595 samples is shown in Fig. \[den11\], using our proposed method to denoise the noisy bunny. The successful recovery of the bunny point cloud presents a strong evidence that our estimated spectral pairs and denoising method are powerful tools in processing noisy datasets.

Conclusions and Future Directions {#con}
=================================

In this work, we develop HGSP tools for effectively processing 3D point clouds. We first introduce a novel method to estimate hypergraph spectral components and presented an optimization formulation to optimally select frequency coefficients to recover the optimal hypergraph structure. We develop a HGSP algorithm to jointly estimate hypergraph spectrum pairs and denoise noisy point clouds. To test the practicality and efficacy of our proposed hypergraph tools, we study two point cloud application examples. Our results illustrate significant performance improvements for both sampling and denoising applications. Moreover, we establish a clear connection between hypergraph frequency components and features on point-cloud surface that can be exploited in future studies.

Our work establish hypergraph signal processing as an efficient tool in tackling high-dimensional interactions among multiple nodes. In addition to sampling and denoising, HGSP can find good applications in many other aspects of point clouds through estimation of spectral components and frequency coefficients. One direction is the design of filters to analyze the spectral properties and surface features of 3D point clouds. Another interesting problem is the recovery of point clouds from low dimensional samples. Beyond point clouds, HGSP can also effectively handle datasets with other complex underlying structure.

[1]{} A. Aldoma, Z. Marton, F. Tombari, W. Wohlkinger, C. Potthast, B. Zeisl, R. B. Rusu, S. Gedikli, and M. Vincze, “Tutorial: point cloud library: three-dimensional object recognition and 6 dof pose estimation," *IEEE Robotics and Automation Magazine*, vol. 19, no. 3, pp. 80-91, Sep. 2012.

R. B. Rusu and S. Cousins, “3d is here: Point cloud li-brary (pcl)," in *2011 IEEE International Conference onRobotics and Automation*, Shanghai, China, May 2011, pp. 1–4.

S. Chen, D. Tian, C. Feng, A. Vetro and J. Kovačević, “Fast Resampling of Three-Dimensional Point Clouds via Graphs," in *IEEE Transactions on Signal Processing*, vol. 66, no. 3, pp. 666-681, Feb, 2018.

J. Ryde, and H. Hu, “3D mapping with multi-resolution occupied voxel lists" *Autonomous Robots*, vol. 28, no. 2, pp. 169-185. Apr. 2010.

A. Hornung, K. M. Wurm, M. Bennewitz, C. Stachniss, and W. Burgard, “Octomap: An efficient probabilistic 3D mapping framework based on octrees," *Autonomous Robots*, pp. 189–206, Apr. 2013.

J. Peng and C. C. Jay Kuo, “Geometry-guided progressive lossless 3D mesh coding with octree (OT) decomposition," *ACM Trans. Graph. Proceedings of ACM SIGGRAPH*, vol. 24, no. 3, pp. 609–616, Jul. 2005.

S. Arya, D. M. Mount, N. S. Netanyahu, R. Silverman, and A. Y. Wu. “An optimal algorithm for approximate nearest neighbor searching fixed dimensions," *Journal of the ACM*, vol.45, np. 6, pp. 891–923, Nov. 1998.

Y. Schoenenberger, J. Paratte, and P. Vandergheynst, “Graph-based denoising for time-varying point clouds," in *2015 3DTV-Conference: The True Vision-Capture, Transmission and Display of 3D Video (3DTV-CON)*, Lisbon, Portugal, Aug. 2015, pp. 1-4.

S. Zhang, Z. Ding, and S. Cui, “Introducing hypergraph signal processing: theoretical foundation and practical applications," arXiv:1907.09203, Oct. 2019.

R. B. Rusu, Z. C. Marton, N. Blodow, M. Dolha, and M. Beetz, “Towards 3D point cloud based object maps for household environments," *Robotics and Autonomous Systems*, vol. 56, no. 11, pp. 927-941, Nov. 2008.

I. Konstas, and M. Lapata, “Unsupervised concept-to-text generation with hypergraphs," in *Proceedings of the 2012 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies*, Montreal, Canada, Jun. 2012, pp.752-761.

A. Banerjee, C. Arnab, and M. Bibhash, “Spectra of general hypergraphs,” *Linear Algebra and its Applications*, vol. 518, pp. 14-30, Dec. 2016.

S. Zhang, H. Zhang, H. Li and S. Cui, “Tensor-based spectral analysis of cascading failures over multilayer complex systems,” in *Proceedings of 56th Annual Allerton Conference on Communication, Control, and Computing (Allerton)*, Monticello, USA, Oct. 2018, pp. 997-1004.

T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” *SIAM Review*, vol. 51, no. 3, pp. 455-500, Aug. 2009.

A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Stationary graph processes and spectral estimation," *IEEE Transactions on Signal Processing*, vol. 65, no. 22, pp. 5911-5926, Aug. 2017.

B. Pasdeloup, V. Gripon, G. Mercier, D. Pastor, and M. G. Rabbat, “Characterization and inference of graph diffusion processes from observations of stationary signals," *IEEE Transactions on Signal and Information Processing over Networks*, vol. 4, no. 3, pp. 481-496, Aug. 2017.

G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro, *Connecting the dots: Identifying network structure via graph signal processing*, *IEEE Signal Processing Magazine*, vol. 36, no. 3, pp. 16-43, May 2019.

B. Pasdeloup, V. Gripon, G. Mercier, D. Pastor, and M. G. Rabbat, “Characterization and inference of graph diffusion processes from observations of stationary signals," *IEEE Transactions on Signal and Information Processing over Networks*, vol. 4, no. 3, pp. 481-496, Aug. 2017.

A. Ortega, P. Frossard, J. Kovacevic, J. M. F. Moura, and P. Vandergheynst, “Graph signal processing: overview, challenges, and applications," *Proceedings of the IEEE*, vol. 106, no. 5, pp. 808-828, Apr. 2018.

K. I. Park, *Fundamentals of probability and stochastic processes with applications to communications*. Springer International Publishing, 2018.

S. Segarra, A. G. Marques, G. Mateos, and A. Ribeiro, “Network topology inference from spectral templates," *IEEE Transactions on Signal and Information Processing over Networks*, vol. 3, no. 3, pp. 467-483, Sep. 2017.

X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning Laplacian matrix in smooth graph signal representations," *IEEE Transactions on Signal Processing*, vol. 64, no. 23, pp. 6160-6173, Dec. 2016.

M. Grant and S. Boyd. “Graph implementations for nonsmooth convex programs," *Recent Advances in Learning and Control (a tribute to M. Vidyasagar)*.

V. Blondel, S. Boyd, and H. Kimura, editors, pages 95-110, *Lecture Notes in Control and Information Sciences*, Springer, 2008. http://stanford.edu/ boyd/graph\_dcp.html.

A. Afshar, J. C. Ho, B. Dilkina, I. Perros, E. B. Khalil, L. Xiong, and V. Sunderam, “Cp-ortho: an orthogonal tensor factorization framework for spatio-temporal data," in *Proceedings of the 25th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems*, Redondo Beach, CA, USA, Jan. 2017, p. 67.

Y. Sun, S. Schaefer, and W. Wang, “Denoising point sets via l0 minimization," *Computer Aided Geometric Design*, vol. 35, pp. 2–15, May 2015.

J. Zeng, G. Cheung, M. Ng, J. Pang, and Yang, C. “3d point cloud denoising using graph laplacian regularization of a low dimensional manifold model," arXiv:1803.07252, 2018. C. Dinesh, G. Cheung, I. V. Bajic, and C. Yang, “Fast 3d point cloud denoising via bipartite graph approximation and total variation," arXiv:1804.10831, 2018.

A. Sandryhaila, and J. M. F. Moura, “Discrete signal processing on graphs: frequency analysis," *IEEE Transactions on Signal Processing,* vol. 62, no. 12, pp. 3042-3054, Apr. 2014. A. M. Bronstein, M. M. Bronstein, and R. Kimmel, “Numerical geometry of non-rigid shapes", Springer, 2008. ISBN: 978-0-387-73300-5.

A. M. Bronstein, M. M. Bronstein, and R. Kimmel, “Efficient computation of isometry-invariant distances between surfaces", *SIAM J. Scientific Computing*, vol. 28, no. 5, pp. 1812-1836, 2006.

A. M. Bronstein, M. M. Bronstein, U. Castellani, B. Falcidieno, A. Fusiello, A. Godil, L. J. Guibas, I. Kokkinos, Z. Lian, M. Ovsjanikov, G. Patané, M. Spagnuolo, and R. Toldo, “SHREC 2010: robust large-scale shape retrieval benchmark", *Proc. EUROGRAPHICS Workshop on 3D Object Retrieval (3DOR)*, 2010.

A. M. Bronstein, M. M. Bronstein, B. Bustos, U. Castellani, M. Crisani, B. Falcidieno, L. J. Guibas, I. Kokkinos, V. Murino, M. Ovsjanikov, G. Patané, I. Sipiran, M. Spagnuolo, and J. Sun, “SHREC 2010: robust feature detection and description benchmark", *Proc. EUROGRAPHICS Workshop on 3D Object Retrieval (3DOR)*, 2010.

S. W. Cheng, and M. K. Lau, “Denoising a point cloud for surface reconstruction," arXiv:1704.04038, 2017.
