Non-uniform Point Cloud Upsampling via Local Manifold Distribution

Optimization-based methods typically require geometric prior information of the point cloud, such as edges and normals, and their performance tends to degrade for complex shapes. Alexa et al. (Alexa et al., 2003) first introduced a point cloud upsampling algorithm, which adds new points at the vertices of the Voronoi diagram in the local tangent space to achieve point cloud upsampling. Subsequently, Lipman et al. (Lipman et al., 2007) proposed a non-parametric method based on the Local Optimal Projection (LOP) operator, which resamples points based on the L1 norm. While this method is robust to noise and outliers, its performance on complex geometries is relatively poor. Later, Huang et al. (Huang et al., 2009) developed an improved weighted LOP, which performs upsampling in an edge-aware manner, reducing the quality degradation at sharp edges and corners. Huang et al. (Huang et al., 2013) also proposed an Edge-Aware Resampling (EAR) method, which preserves edges during point cloud upsampling but heavily depends on the provided normal information and parameter tuning. Overall, these methods heavily rely on prior information and place significant demands on the original point cloud.

With the development of deep learning techniques and inspired by the success of PointNet (Qi et al., 2017a), many deep learning methods have been proposed for various point cloud processing tasks such as classification, completion, denoising, and other applications. For point cloud upsampling, PU-Net (Yu et al., 2018b) was the first to apply neural networks to point cloud upsampling, using PointNet++ (Qi et al., 2017b) to extract point features and extending those features through a multi-branch MLP to generate upsampled point clouds. However, this approach does not consider the spatial relationships between points, and thus cannot guarantee the generation of uniformly distributed points. EC-Net (Yu et al., 2018a) employs a joint loss of point-edge distances, allowing for the generation of sharp new points at edges. MPU (Yifan et al., 2019) is a progressive method for upsampling point patches, but it struggles to generate new points in missing regions for non-uniform point clouds. PU-GAN (Li et al., 2019) was the first to apply Generative Adversarial Networks (GANs) to point cloud upsampling, introducing uniformity loss and becoming the first method capable of generating uniformly distributed upsampled point clouds. PU-GCN (Qian et al., 2021) uses graph convolutional networks for feature extraction and expansion, proposing a novel point cloud upsampling module called NodeShuffle. PUGeo (Qian et al., 2020) is the first method to utilize geometric approaches for point cloud upsampling, improving upsampling results by employing local differential geometry constraints. Dis-PU (Li et al., 2021) defines two cascaded subnetworks that complete point cloud upsampling in stages. NP (Feng et al., 2022) introduces the use of neural networks to continuously represent geometric shapes, enabling sampling at arbitrary resolutions. Recently, GeoUDF (Ren et al., 2023) introduced a local geometric representation that approximates local shapes using quadratic surfaces. RepKPU (Rong et al., 2024) proposed a new paradigm, kernel-to-displacement generation, for point generation. However, existing learning-based methods often fail to achieve good performance when handling highly sparse and non-uniform datasets. Drawing inspiration from LGSur-Net (Xiao et al., 2024), our approach performs Gaussian fitting on the local patches of point clouds and applies further constraints on the statistical manifold constructed from the parameters. Our method generates more uniform results for sparse and non-uniform point clouds.

Point cloud data typically originate from the surfaces of objects in three-dimensional space, and these surfaces can locally be approximated as manifolds. A manifold is a topological space, where an n-dimensional topological manifold $\mathcal{M}$ is a locally Euclidean, second-countable Hausdorff space. In other words, every portion of the manifold can locally be homeomorphic to Euclidean space. Therefore, we represent each local neighborhood of a point cloud as a manifold. For the selection of basis functions for fitting the surface, the fitting capability of Gaussian functions has been widely proven. Due to the excellent linear combination properties of Gaussian basis functions, they can form complex functions through weighted summation, thus capable of expressing surfaces of various shapes. Moreover, the shape of a Gaussian function is well captured by its covariance matrix. Furthermore, inspired by (Yu and Turk, 2013), we incorporate anisotropic kernels to accommodate the irregular distribution of particles. Additionally, anisotropic kernel regression and covariance analysis enable the accurate capture of the local distribution characteristics of fluid surfaces. Considering that a point cloud can be viewed as a sample drawn from some probability distribution, i.e., the distribution of local point cloud data can be considered as a probability distribution, the Gaussian distribution is a natural choice. It can effectively and simply describe the central tendency and the direction of spread of the data through its mean and covariance matrix. Due to the uniformity of the basis functions, we construct a unified manifold representation for the entire point cloud. Given the distributional properties of Gaussian functions, we use a statistical model to describe the structure and geometric shape of the manifold. As a result, the statistical manifold can be defined by the properties of the Gaussian functions themselves.

That is, we map the manifold representation of each neighborhood in the point cloud to a unified manifold. Since each point on the manifold represents a probability distribution, it is referred to as a statistical manifold.

Considering the local manifold properties of point clouds, we perform a fitting representation of the local neighborhoods. Let the input point cloud be denoted as $P$. To ensure the overlap of neighborhood regions and prevent gaps during surface fitting, we first randomly select a certain number of points from the original point cloud as initial query points. This selection process ensures a uniform distribution, avoiding excessive local concentration or sparsity. Subsequently, each local patch is updated by incorporating points based on nearest-neighbor neighborhoods. Finally, the centroid of each patch is computed and designated as the final query point. For each query point $P_{i}$ , we construct a local coordinate system, and the projection plane is denoted as $D_{xy}$. The specific transformation formula for the coordinates on the projection plane is given below:

We treat the local neighborhood of point $P_{i}$ as a 2D manifold in 3D space. Since the manifold surface in 3D space is homeomorphic to a 2D region, for each point and its local neighborhood on the manifold, there exists a local coordinate system such that this neighborhood can be mapped to a 2D Euclidean space via a homeomorphic mapping. Therefore, a 3D local neighborhood patch can be isomorphic to the 2D parameter domain D. To eliminate the complexity caused by surface tilt in the original global coordinates, we construct a 2D parameter domain $\mathcal{D}$ for the projection plane. This means that we can establish a continuous mapping between the original local neighborhood and $\mathcal{D}$:

Considering the challenges of directly fitting the point cloud surface in three-dimensional space, and given that local surfaces are inherently two-dimensional, it is unnecessary to use three-dimensional Gaussian functions to describe the redundant dimension. Instead, we perform surface function fitting in a two-dimensional parameter domain and transform the parameters back to represent the original point cloud surface. In the local coordinate region, we define a series of Gaussian kernels centered at the base points.

To ensure the fitting accuracy for sparse and non-uniform point clouds, we design the number of Gaussian functions $t$, based on the number of points within the neighborhood. Furthermore, drawing inspiration from the data-adaptive kernel regression method proposed in (Takeda et al., 2007), we capture the local characteristics of point clouds by adaptively adjusting the shape and size of the kernel function. To reduce the optimization difficulty, we fix the Gaussian function centers and only optimize the covariance matrix for each function. Through our design, the manifold of each local neighborhood in the point cloud is described as a weighted sum of multiple Gaussian functions, with the local point cloud distribution represented by the covariance matrices of the Gaussian components. Considering the statistical properties of the covariance matrix, we adopt the concept of Gaussian splatting(Kerbl et al., 2023), to ensure the positive definiteness of the covariance matrix, we decompose it into the product of a rotation matrix $R$ and a scaling matrix $S$, as expressed by the following equation:

During the entire optimization process, we refine the rotation matrix $R$ and scaling matrix $S$ to maintain the positive definiteness of the covariance matrix. Finally, we fit the local neighborhood manifold of the point cloud in the 2D parametric domain, expressed as follows:

In the weighted combination of Gaussian basis functions, although the basis functions themselves define part of the shape of the local surface, they do not directly define the position of the surface. The weighted Gaussian basis functions in equation (3) capture the shape features of the surface, but they do not specify the exact position of these shapes in three-dimensional space. Therefore, we design the local surface fitting for a local patch near a query point $q\in\mathbb{R}^{3}$ as follows:

This is our local Gaussian representation. Although our Gaussian representation inherently possesses strong approximation capabilities, extracting a reasonable Gaussian representation from information-deficient and sparse point clouds presents a significant challenge. During the optimization process of the above formulation, we leverage a generative model to learn and shape a prior. Specifically, we incorporate modules such as EdgeConv and cross-attention to account for the geometric structure of local neighborhoods, facilitating feature aggregation. The covariance matrix is optimized to model the Gaussian distribution, while the optimized Gaussian weight matrix yields the best parameter approximation, ultimately generating a Gaussian representation of the local neighborhood. Our approach ensures a refined consideration of the interactions between shallow and deep features, leveraging the powerful capabilities of neural networks to enable a more nuanced and detailed representation.

To better generate upsampled point clouds with a uniform distribution, we first resample a sufficient number of coordinate points in the 2D parameter domain of each local neighborhood and obtain the upsampled point cloud through local Gaussian representation mapping. Additionally, current point cloud upsampling methods fail to leverage the distribution information of the point cloud to supervise the generation of upsampled points. From the statistical manifold constructed by the Gaussian representations, we propose a novel method for correcting the distribution of the point cloud.

Since each point on the manifold represents a probability distribution, to ensure that the upsampled point cloud better aligns with the true distribution, the problem can be reformulated such that the distance between corresponding points on the statistical manifold is minimized. Here, we choose the Fisher-Rao distance between two distributions. The Fisher-Rao distance between two probability distributions on the statistical manifold is the length of a geodesic, which is the local shortest distance between the two points. This geodesic is defined under the Fisher information metric between the distributions. The elements of the Fisher information matrix and the classic Fisher-Rao distance are defined as follows:

Since our local Gaussian representation consists of multiple Gaussian functions (normal distributions), the Fisher-Rao distance between two multivariate normal distributions can be expressed as:

Following the approach in (Costa et al., 2015), we compute the Fisher-Rao distance for a univariate Gaussian probability distribution. Since the hyperbolic distance between two points $\boldsymbol{x}=\left(x_{1},x_{2}\right)$ and $\boldsymbol{y}=\left(y_{1},y_{2}\right)$ in the Poincaré half-plane model is given by:

We derive and establish the relationship between the Fisher-Rao distance for a univariate Gaussian probability distribution and the Poincaré distance $d_{H}$ within the Poincaré half-plane model in hyperbolic geometry:

Thus, the final computation formula is given as follows:

Now, to constrain the distribution of the point cloud, the objective is to minimize:

In this section, we introduce the design of the loss function, which includes constraints on the quality of the point cloud, the similarity between distributions, and the smoothness of the surface shape. First, considering the importance of point cloud quality for point cloud upsampling, our method should encourage the generated points to be closer to the real point cloud. Therefore, a reconstruction loss function is needed to evaluate the similarity between the two point clouds. Here, the reconstruction loss uses the Earth Mover’s distance(EMD).

Then, during the minimization process of Eq. (14), to prevent interference from outliers in some distributions, we impose a first-order differential constraint on the surface function, as shown in Eq. (16). This constraint limits the gradient variation of the surface, promoting smooth transitions in adjacent regions and enhancing surface continuity. Additionally, smoothing the gradient also contributes to improving the uniformity of the point cloud distribution.

Finally, based on the above description, the final combined loss function is as follows:

Compared to the traditional PUGAN (Li et al., 2019)dataset, the PU1K dataset (Qian et al., 2021) is more challenging because it contains a larger volume of data and more diverse categories. This dataset includes large objects with complex shapes, as provided by PUGCN. The PU1K dataset is collected from both PUGAN and ShapeNetCore (Xu et al., 2019). Since our method is designed for upsampling non-uniform point clouds, we processed the PU1K dataset to account for non-uniformity. During the training phase, we first sampled a dense point cloud from the original mesh as ground truth. For each patch during training, we applied non-uniform sampling as the input data.

For the testing dataset, we selected different resolutions of the PU1K dataset for evaluation. Considering the testing requirements of our method for upsampling sparse, non-uniform point clouds, we applied non-uniform processing to the original dataset as well. Additionally, to validate the generalization of our method and assess performance at different upsampling scales, we also evaluated the Sketchfab dataset (Qian et al., 2020) and the KITTI dataset (Geiger et al., 2013).

All of our experiments are implemented using PyTorch. We set the cross-attention and graph convolution modules to three layers and the number of Gaussian factors per patch to $t=\left\lfloor\mathcal{N}(p_{i})/8\right\rfloor$. We use the Adam optimizer to train our model on an RTX 3090 GPU, with a batch size of 64, an upsampling factor of 4, and for 400 epochs. The initial learning rate is set to 0.001, decaying by a factor of 0.7 every 80 iterations. Additionally, to eliminate unnecessary degrees of freedom in the input data space and reduce the learning difficulty, we normalize each point’s coordinates by its patch radius and rotate the points into a local coordinate system defined by PCA. To avoid overfitting during training, we augment the network inputs with random rotations, scaling, and Gaussian noise perturbations.

Similar to recent point cloud upsampling works, for quantitative evaluation, we use three common metrics: Chamfer Distance ($\text{CD}\times 10^{-3}$), Hausdorff Distance ($\text{HD}\times 10^{-4}$), and Point-to-Surface Distance ($\text{P2F}\times 10^{-3}$). These metrics reflect the reconstruction quality of the upsampled point cloud. Additionally, since our method considers the distribution similarity of each patch, we include the Jensen-Shannon Divergence ($\text{JSD}\times 10^{-3}$) metric, which measures the similarity between point cloud distributions. Finally, we introduce the Uniformity ($\text{UNI}\times 10^{-2}$) metric proposed in PU-GAN (Li et al., 2019), which reflects the uniformity of the point cloud. For all metrics, smaller values indicate better performance.

We compare the method proposed in this paper with five existing point cloud upsampling methods, including PUGeo (Qian et al., 2020), PUCRN (Du et al., 2022), APUNet(Zhao et al., 2023), Grad-PU(He et al., 2023), and RepKPU(Rong et al., 2024). For a fair comparison, we retrained these methods on our dataset in the same experimental environment, using the official code and recommended settings.

The quantitative comparison results are shown in Table 1. We performed a quantitative analysis on sparse, non-uniform point clouds with 256, 1024, and 4096 input points. From the table, we can see that our method achieves superior results in upsampling sparse, non-uniform point clouds compared to most other methods. Although, for 256 input points, some metrics like HD and P2M are not as good as RepKPU, the visual results in Figure 3 show that RepKPU produces highly uneven results. This is because when dealing with extremely sparse point clouds, anomalies can occur during point generation, affecting some metrics significantly. In contrast, other methods tend to generate highly clustered points, which may yield better results for some metrics. However, the visual results clearly demonstrate that our method generates points that most closely align with the original features.

In addition to the quantitative results, we also present point cloud upsampling and surface reconstruction results for some models in Figure 2. Here, we use Poisson reconstruction (Kazhdan and Hoppe, 2013). As seen in the figure, when dealing with complex data, such as the elephant, our method produces smoother and more uniformly distributed points compared to other methods, providing excellent guidance for subsequent reconstruction. Our upsampling results effectively fill in gaps, while other methods tend to generate more noise and uneven point sets.

Additionally, to evaluate runtime performance, we compare these methods in terms of floating point operations (FLOPs), as well as the training and the inference time. As shown in Table 2, our method exhibits the shortest training time while also achieving competitive performance in FLOPs and inference efficiency.

To validate the robustness of our method against noisy results, we performed upsampling on point clouds with noise levels of 0.5%, 1%, and 2%. As shown in Table 2, our method provides the best results across different noise levels. From Figure 4, we can observe that our method effectively constrains the noisy outliers to be closer to the surface, and our method also achieves the best performance in terms of uniformity.

To validate the effectiveness of our method in real-world scenarios, we conducted experiments on the KITTI dataset. Since there is no ground truth point cloud, we provide a visual comparison in Figure 5. As seen, even when facing real scanned point clouds, our method is still able to fill in gaps and output a more uniform point distribution.

To verify the effectiveness of each module, we conduct an ablation study to demonstrate how each component impacts the final results. Specifically, we use the PU1K dataset as a benchmark. We focus on the selection of Gaussian functions(I), the distribution constraints on the manifold(II), and the design of the first-order derivative loss(III). For the function selection, we adopt the quadratic term representation of surfaces provided in GeoUDF (Ren et al., 2023). The quantitative results of the ablation study are shown in Table 3. It can be observed that each component we designed significantly improves the experimental performance. From Figure 6, it is evident that without the manifold distribution constraint, the points become more clustered. The absence of the first-order derivative loss results in interference from anomalous points. The experimental results strongly support the effectiveness and necessity of each strategy we employed.