3D Point Cloud Fusion Method Based on EMD Auto-Evolution and Local Parametric Network

Abstract

Although the development of high-resolution remote sensing satellite technology has made it possible to reconstruct the 3D structure of object-level features using satellite imagery, the results from a single reconstruction are often insufficient to comprehensively describe the 3D structure of the target. Therefore, developing an effective 3D point cloud fusion method can fully utilize information from multiple observations to improve the accuracy of 3D reconstruction. To this end, this paper addresses the problems of shape distortion and sparse point cloud density in existing 3D point cloud fusion methods by proposing a 3D point cloud fusion method based on Earth mover’s distance (EMD) auto-evolution and local parameterization network. Our method is divided into two stages. In the first stage, EMD is introduced as a key metric for evaluating the fusion results, and a point cloud fusion method based on EMD auto-evolution is constructed. The method uses an alternating iterative technique to sequentially update the variables and produce an initial fusion result. The second stage focuses on point cloud optimization by constructing a local parameterization network for the point cloud, mapping the upsampled point cloud in the 2D parameter domain back to the 3D space to complete the optimization. Through these two steps, the method achieves the fusion of two sets of non-uniform point cloud data obtained from satellite stereo images into a single, denser 3D point cloud that more closely resembles the true target shape. Experimental results demonstrate that our fusion method outperforms other classical comparison algorithms for targets such as buildings, planes, and ships, and achieves a fused RMSE of approximately 2 m and an EMD accuracy better than 0.5.

1. Introduction

Three-dimensional (3D) reconstruction is a critical step in the interpretation and processing of remote sensing images. Three-dimensional data can provide a more comprehensive, detailed, and accurate stereoscopic description, allowing for more intuitive and precise judgment and understanding. The 3D reconstruction of satellite remote sensing images has traditionally been within the research domain of photogrammetry, primarily relying on the results of satellite stereo image (SSI) matching combined with a rational function model to determine the spatial coordinates of object points. However, due to the resolution limitations of SSIs, the accuracy of reconstruction is often insufficient to meet the requirements for high-precision 3D descriptions of targets. With the development of deep learning technology and the use of large-scale training datasets in computer vision, it has become possible to use convolutional neural networks to reconstruct high-precision terrestrial targets in satellite remote sensing images. Sebastien et al. [1] proposed a method for 3D reconstruction of building targets based on SSIs. This method utilizes a U-net to extract the building’s polygonal outlines, combined with a digital terrain model and a digital elevation model, to accurately extract the building model. Yi et al. [2] introduced a machine learning-based method for automatic 3D building reconstruction and vectorization. This approach uses the digital surface model (DSM) and panchromatic images as inputs, combining conditional generative adversarial networks (cGANs) with semantic segmentation networks to extract building target models and construct roof polygons. Chen et al. [3] proposed a 3D extraction algorithm for object targets of SSIs, which leverages RFM constraint deformation reasoning networks and self-similar convolutions. This method integrates DSMs and SSIs to reconstruct the 3D structure of terrestrial targets.

Due to the limitations of the satellite observation perspective, the results obtained from a single pair of stereo images are often insufficient to effectively describe the complete 3D structure of the target. On the other hand, with the increasing volume of satellite image data, both the same and different satellites can perform multiple observations of the same object within a short period. This leads to multiple sets of 3D reconstruction results for the target structure. Therefore, it is necessary to develop an effective 3D point cloud fusion method, utilizing the multiple 3D reconstruction results obtained from different observations to achieve a more accurate fusion of the target’s 3D point cloud, ensuring the precision of the reconstructed 3D structure.

Three-dimensional data fusion methods are primarily divided into traditional methods and deep learning methods. Traditional methods mainly focus on 3D registration, where two sets of 3D point cloud are first unified in the same coordinate system and then fused using a certain measure or direct merging. The most classic algorithm in 3D point cloud registration is the iterative closest point (ICP) algorithm [4]. This algorithm works by finding corresponding points between the source and target point cloud, constructing a rotation and translation matrix based on these corresponding points, and transforming the source point cloud into the target point cloud’s coordinate system using the computed matrix. The error function between the transformed source point cloud and the target point cloud is then estimated. If the value of the error function exceeds a threshold, the above operations are iteratively performed until the specified error requirement is met. Based on this, Peng et al. [5] utilized ICP to register two cross-source point cloud from structure from motion (SfM) and light detection and ranging (LiDAR). However, their approach made several assumptions, including removing outliers and manually selecting targets. Zhou et al. [6] proposed a fast global registration (FGR) method to overcome the limitations of traditional ICP. This method provides precomputed fast point feature histograms (FPFHs) [7] between two surfaces, allowing dense alignment at each iteration without the need for initialization or correspondence searching while optimizing point registration. Tsin et al. [8] proposed a similar kernel-based matching approach where the registration between two sets of point cloud is formulated as maximizing kernel correlation (KC). Due to the weighted contributions of multiple linked dynamic points, the KC objective is smooth and exhibits unique convergence properties. However, this method is susceptible to noise and outliers within objects. Chen et al. [9] proposed a 3D fusion method for combining 3D data from laser scanning and SfM. By using a scale-based PCA-ICP algorithm, the scale differences between two viewpoints are eliminated, feature points are automatically extracted for precise registration, and multiple point cloud are registered using an optimized ICP method, achieving robust and accurate 3D scene reconstruction. Huang et al. [10] proposed a systematic method for registering cross-source point cloud using a scale normalization technique to eliminate scale issues and employing a new graph construction method to jointly reconstruct 3D structures from different sources. András et al. [11] introduced a point cloud blending and volumetric fusion method based on 3D tetrahedralization (3DT) ray projection, which unifies point cloud from SfM or multi-view stereo (MVS) with coarser but more complete point cloud from airborne 3D reconstruction. Hoegner et al. [12] fused point cloud reconstructed from ground-based laser scanners and RGB cameras with thermal infrared images mounted on robots for indoor 3D reconstruction.

Deep learning-based 3D data fusion methods are commonly used in 3D reconstruction across various scenarios. Typically, these methods generate the 3D structure of an entire scene by fusing multiple sets of reconstruction results (such as depth map, point cloud, or mesh) in the final step. Li et al. [13] proposed a dual-layer neural fusion network that leverages recent advances in neural implicit representation and neural rendering to fuse multi-view depth maps into a 3D structure. Jia et al. [14] developed a region attention mechanism for 3D reconstruction to learn high-quality correspondences between two coarse point cloud, which were then fused into a high-precision point cloud. Jia et al. [15] also introduced a gated network for point cloud data fusion that first predicts 3D point cloud from dual-view RGB images and then fuses two sets of point cloud obtained from different views. Liu et al. [16] proposed a method for generating 3D building models based on the fusion of 3D point cloud and mesh data. This method uses a designed multi-source 3D data quality evaluation network (MS3DQE-Net) to assess the quality of 3D meshes and point cloud and guides the fusion of 3D building models based on the evaluation results. Additionally, using deep learning to enhance the ICP algorithm for higher precision 3D registration and fusion has become a new hotspot in the field of 3D fusion in recent years. Wang et al. [17] proposed a deep closest point (DCP) network based on deep learning to address local optima and other challenges in ICP, which is used for point cloud registration. Hu et al. [18] presented a 3D point cloud enhancement optimization method that combines point cloud obtained from a binocular structured light scanner with those reconstructed by the Colmap algorithm. This method generates a more complete and higher-quality dense point cloud by integrating the coarse registration of the PointNetLK network with the fine registration of the ICP algorithm. Yookwan et al. [19] proposed a method for fusing point cloud extracted from depth values, which are directly measured by an infrared camera and estimated from RGB images of the same scene using an improved ResNet-50. They also introduced an information theory alignment strategy using cross entropy ICP for point cloud registration.

The typical 3D data fusion methods mentioned above still face the following challenges when applied to SSIs 3D data fusion:

• SSIs data often have low resolution and limited observation angles, making it difficult to achieve the level of point cloud accuracy and density required by the typical 3D fusion methods mentioned earlier. Directly applying these methods can lead to a disordered target structure after fusion due to the insufficient accuracy and density of the original data.
• Most of the typical fusion methods discussed above focus on constructing a measure to minimize the absolute spatial distance between the registration or fusion results and each source of data. However, minimizing absolute spatial distance does not necessarily ensure that the fused shape is closer to the true target. There may exist fusion results with a smaller absolute spatial distance but greater shape distortion. The objects considered in this paper, such as planes, ships, and buildings, have distinct inherent features in human cognition, and therefore require a higher degree of shape authenticity in the fusion results.
• The typical fusion methods described above do not significantly increase the density of the point cloud after fusion, resulting in a substantial gap between the fused 3D structure and the true target.

To address the aforementioned challenges, this paper proposes a 3D point cloud fusion method based on Earth mover’s distance (EMD) auto-evolution and local parameterization network by combining the advantages of traditional methods and deep learning networks. The proposed method merges two sets of non-equally spaced point cloud data into a single set of 3D point cloud that more closely resemble the actual target shape and have a higher point cloud density. Our approach consists of two main stages. In the first stage, to tackle the fusion measurement problem, we introduce the EMD measure as a key metric for evaluating the fusion results. By incorporating steps such as point cloud geometric structure inheritance, point cloud auto-evolution, adaptive point cloud weighting, and symmetry constraints, we formulate an optimization problem for the EMD auto-evolution-based point cloud fusion method. Then, using an alternating iterative technique, we sequentially optimize the variables to produce preliminary fusion results. The second stage addresses the point cloud density issue. Based on the preliminary fusion results, a local parameterized network is constructed to optimize the fused point cloud. A joint loss is designed to train the network in an end-to-end manner. Specifically, the preliminary fused point cloud is parameterized into a 2D domain, where multi-scale pointwise features of the 3D point cloud are extracted and connected. After upsampling in the 2D parameter domain, the features are mapped back to 3D space through a linear transformation, yielding the optimized point cloud fusion result.

The main contributions of our work can be summarized as follows:

(1)

We propose a 3D point cloud fusion method based on EMD auto-evolution, which achieves the fusion of non-equally spaced point cloud data while ensuring that the fusion results more closely resemble the actual target shape.

(2)

We propose a point cloud optimization method based on local parameterization network that enhances point cloud density while recovering more of the points lost during the initial fusion process.

The remainder of this paper is organized as follows. Section 2 provides a detailed de-scription of the proposed method. Section 3 presents comparative experiments and a discussion of public data sets. Conclusions are drawn in Section 4.

2. The Proposed Method

2.1. 3D Point Cloud Fusion Based on EMD Auto-Evolution

To ensure the effectiveness of 3D point cloud fusion, it is necessary to establish a target-oriented accuracy metric for 3D point cloud. Commonly used metrics for evaluating the accuracy of 3D reconstruction in photogrammetry include mean elevation deviation and root mean square error (RMSE). However, these metrics are often insufficient for effectively assessing the accuracy of the 3D structure of a target. To accurately evaluate the correctness of a target’s 3D structure, this paper introduces the EMD measure, a typical metric in the field of computer vision, as an indicator of the accuracy of 3D point cloud reconstruction. Unlike commonly used distance metrics such as Euclidean distance or Chebyshev distance, the EMD measure is suitable for evaluating dissimilarity in non-uniformly sized data. This means it can be applied to assess the accuracy of 3D reconstruction when there is an inconsistency in the number of points within the point cloud. Additionally, EMD provides a shape dissimilarity measurement that aligns with human perception, making it particularly well suited for evaluating the accuracy of the 3D structure of a target.

Considering both the EMD measure and the practical need for immediate data acquisition and precise reconstruction from satellite data, this paper first constructs a 3D point cloud fusion method based on EMD auto-evolution. The method uses the 3D reconstruction results obtained from SSIs as a basis to achieve the fusion and updating of the target 3D point cloud, thereby ensuring the accuracy of the target’s 3D reconstruction results.

This work uses $P^a=[P_1^a,\cdots ,P_{2n}^a]$ to represent one target 3D point cloud a reconstructed from a set of SSIs and $P^b=[P_1^b,\cdots ,P_{2m}^b]$ to represent the same target 3D point cloud b reconstructed from a different set of SSIs. The fused target 3D point cloud is denoted by $P=[P_1,\cdots ,P_k]$. For the fusion algorithm, it is essential not only to effectively reflect the geometric structure of the target at the current moment but also to preserve the original geometric structure of the target. Given that the EMD can effectively measure the dissimilarity in the geometric structure of point cloud, a point cloud dissimilarity measurement method based on EMD will be embedded into the fusion algorithm to ensure the structural similarity between the fused point cloud and the input point cloud. Additionally, due to the limited observation perspectives of SSIs, the target point cloud obtained from satellite observations may contain points with low confidence. To prevent low-confidence outlier points from adversely affecting the 3D fusion of the target, it is necessary to eliminate these outlier points as much as possible during the fusion process. Therefore, this paper proposes an EMD auto-evolution point cloud fusion method, which gradually approximates the geometric structure of the target point cloud during the evolution process and effectively reduces the influence of potential outlier points in the point cloud on the fusion result, based on the EMD measure. The final output is a preliminary fusion result with high confidence. The framework of proposed method is illustrated in Figure 1.

2.1.1. EMD Distance

The EMD was first introduced by Stanford University in 2000 [20], where it was demonstrated to effectively capture image dissimilarity in a manner consistent with human perception. EMD is suitable for comparing datasets of unequal length, which makes it particularly useful for assessing the similarity of point cloud when there is an inconsistency in the number of points. EMD is often used as a metric to measure the accuracy of 3D reconstruction. When applied to compare the dissimilarity between two target point cloud, it can be expressed as the following formula:

$$\begin{array}{l}F={\displaystyle \sum _{j=1}^{2m}{\displaystyle \sum _{i=1}^{2n}{s}_{i,j}\cdot w_{i,j}}}\\ s.t.{\displaystyle \sum _{j=1}^{2n}{s}_{i,j}\le 1}, {\displaystyle \sum _{i=1}^{2m}{s}_{i,j}\le 1}\\ {\displaystyle \sum _{j=1}^{2m}{\displaystyle \sum _{i=1}^{2n}{s}_{i,j}=\min (2m,2n)}}\\ s_{i,j}\ge 0\end{array}$$

(1)

where 2m and 2n denote the total number of points in each of the two point cloud, and the coefficient 2 is for ease of calculation. $w_{i,j}$ represents the cost of transporting the i-th point from the source point cloud to the j-th point in the target point cloud, and $s_{i,j}$ indicates whether the i-th point in the source point cloud corresponds to the j-th point in the target point cloud. Note that EMD distance does not assume $m=n$, meaning that the number of source point cloud and target point cloud may not be consistent, which makes it applicable for comparing non-equal length data dissimilarity. The goal of this formula is to find a flow that minimizes the objective function F.

2.1.2. Inheritance of Geometric Structure Information of Input Point Cloud

Considering that objects such as planes, ships, and buildings have distinct inherent characteristics in human cognition, simplified templates for these types of objects are constructed to ensure that the 3D data of objects obtained from different sources can be registered to the same coordinate system through 3D registration algorithms [4].

For the fused 3D point cloud $P=[P_1,\cdots ,P_k]$, it needs to have a similar geometric structure to the two input point cloud $P^a=[P_1^a,\cdots ,P_{2n}^a]$ and $P^b=[P_1^b,\cdots ,P_{2m}^b]$. If EMD distance is used to describe the geometric structure similarity of the target, it is necessary to minimize $EMD(P,P^a)$ and $EMD(P,P^b)$, as shown in the following formula:

$$\begin{array}{l}F*={\displaystyle \sum _{j=1}^k{\displaystyle \sum _{i=1}^{2n}{s}_{i,j}^a\cdot w_{i,j}^a}}+\lambda \times {\displaystyle \sum _{j=1}^k{\displaystyle \sum _{i=1}^{2m}{s}_{i,j}^b\cdot w_{i,j}^b}}\\ s.t.{\displaystyle \sum _{i=1}^{2n}{s}_{i,j}^a=1},{\displaystyle \sum _{i=1}^{2m}{s}_{i,j}^b=1}\\ 1\le {\displaystyle \sum _{j=1}^k{s}_{i,j}^a+s_{i,j}^b}\le 2\\ s_{i,j}^a\ge 0,s_{i,j}^b\ge 0\end{array}$$

(2)

where $\lambda$ is the weight factor, $s_{i,j}^a$ represents the correspondence between the i-th point of the reconstruction result of group a and the j-th point after fusion, $s_{i,j}^b$ represents the correspondence between the i-th point of the reconstruction result of group b and the j-th point after fusion, $w_{i,j}^a$ represents the Euclidean distance between the i-th point of the reconstruction result of group a and the j-th point after fusion, $w_{i,j}^b$ represents the Euclidean distance between the i-th point of the reconstruction result of group b and the j-th point after fusion. The constraints ${\displaystyle \sum _{i=1}^{2n}{s}_{i,j}^a=1}$ and $s_{i,j}^a\ge 0$ are employed to ensure that the point cloud of reconstruction result group a establishes an injective relationship with the fused point cloud, thereby retaining as much of the original point cloud information as possible. Similarly, the constraints ${\displaystyle \sum _{i=1}^{2m}{s}_{i,j}^b=1}$ and $s_{i,j}^b\ge 0$ are employed to ensure that the point cloud of reconstruction result group b establishes an injective relationship with the fused point cloud, thereby maximizing the extraction of all the information from the newly obtained point cloud. The constraint condition $1\le {\displaystyle \sum _{j=1}^k{s}_{i,j}^a+s_{i,j}^b}\le 2$ is used to ensure that each point in the fused result has a mapping relationship with the input point cloud, and that any point cannot link different point from the original target 3D reconstruction results.

Based on the above formula, the similarity between the fused point cloud P and the two input reconstruction result point cloud $P^a$ and $P^b$ can be calculated. However, relying solely on the formula cannot determine the specific values of the fused point cloud P. To obtain a fused point cloud with a geometric structure similar to the input point cloud, this paper constructs an auto-evolution point cloud fusion strategy, which iteratively updates the fused point cloud to effectively inherit the geometric structure information of the input point cloud.

2.1.3. Auto-Evolution of Point Cloud

Define the coordinates of the i-th point in the first set of reconstructed point cloud $P^a=[P_1^a,\cdots ,P_{2n}^a]$, as $[P_{i,x}^a,P_{i,y}^a,P_{i,z}^a]$. The coordinates of the i-th point in the second set of reconstructed point cloud $P^b=[P_1^b,\cdots ,P_{2m}^b]$, as $[P_{i,x}^b,P_{i,y}^b,P_{i,z}^b]$. The coordinates of the j-th point in the fused point cloud $P=[P_1,\cdots ,P_k]$, as $[P_{j,x},P_{j,y},P_{j,z}]$. To obtain the precise coordinates of the fused point cloud P, an auto-evolution point cloud fusion problem is formulated based on the evaluation of geometric structure similarity between the fused point cloud and the input point cloud, as shown in the following expression:

$$\begin{array}{l}\underset{s^a,s^b,P}{\min }{\displaystyle \sum _{j=1}^k{\displaystyle \sum _{i=1}^{2n}{s}_{i,j}^a}}\cdot {‖[P_{i,x}^a,P_{i,y}^a,P_{i,z}^a]-[P_{j,x},P_{j,y},P_{j,z}]‖}_2\\ +\lambda \times {\displaystyle \sum _{j=1}^k{\displaystyle \sum _{i=1}^{2m}{s}_{i,j}^b}}\cdot {‖[P_{i,x}^b,P_{i,y}^b,P_{i,z}^b]-[P_{j,x},P_{j,y},P_{j,z}]‖}_2\end{array}$$

(3)

The self-evolutionary point cloud fusion problem includes the variables $s^a$, $s^b$, and P. The objective of this optimization is to search for the possible coordinates of point cloud P such that it maintains a geometric structure similar to that of the input point cloud. Given that this is a non-convex problem, an alternating iterative strategy is employed to solve it. During updating of the variable P, the variables $s^a$ and $s^b$ are kept fixed. The term auto-evolution refers to the process in which, during each iteration, the correspondence between the automatically searched point cloud P and the input point cloud is established for points with the same name, gradually converging to the optimal solution for the point cloud coordinates.

2.1.4. Weight Adaptation for Outlier Point Cloud

Due to errors such as incorrect stereo matching, local deformation of the target point cloud may occur. These local deformations, or outlier point cloud, cannot accurately reflect the geometric structure of the target. Therefore, in the process of point cloud fusion, smaller weights should be assigned to these outlier point cloud to mitigate their adverse effects. During the evolution of the point cloud, the correspondence between the fused point cloud and the input point cloud is established autonomously. Considering that the outlier point cloud deviate significantly from their correct coordinates, the Euclidean distance between the outlier point cloud and its corresponding merged point cloud will be higher than that between other corresponding point cloud. Based on this characteristic, an adaptive weighting technique is embedded in the evolution process of the point cloud to minimize the impact of outlier point cloud on the accuracy of the fused point cloud, as shown in the following equation:

$$\begin{array}{l}\underset{s^a,s^b,P,c^a,c^b}{\min }{\displaystyle \sum _{i=1}^{2n}{c}_i^a}{\displaystyle \sum _{j=1}^k{s}_{i,j}^a}\cdot {‖[P_{i,x}^a,P_{i,y}^a,P_{i,z}^a]-[P_{j,x},P_{j,y},P_{j,z}]‖}_2\\ +\lambda \times {\displaystyle \sum _{i=1}^{2m}{c}_i^b}{\displaystyle \sum _{j=1}^k{s}_{i,j}^b}\cdot {‖[P_{i,x}^b,P_{i,y}^b,P_{i,z}^b]-[P_{j,x},P_{j,y},P_{j,z}]‖}_2\\ s.t.c^{aT}{c}^a=1,c^{bT}{c}^b=1\\ c_i^a\ge 0,c_i^b\ge 0\end{array}$$

(4)

where $c^a$ and $c^b$ are adaptive weighting factors. As fusion variables, they are iteratively updated during the auto-evolution of the point cloud. These factors adaptively select outlier point cloud and assign them smaller weights to reduce their impact on the accuracy of the fused point cloud, ensuring effective fusion of the input point cloud.

2.1.5. Point Cloud Symmetry Constraint

Many artificial structures exhibit obvious symmetry, such as planes, certain ships, and building targets, which can be considered symmetric about their principal planes. To improve the accuracy of target fusion, symmetry constraints are introduced during the point cloud fusion process, ensuring that the fused 3D point cloud retains the inherent symmetric structure of the target. On the other hand, symmetry constraints can reduce the number of point cloud to be fused, thereby improving the execution speed of the point cloud fusion algorithm.

Since the input point cloud has already been aligned to a simplified template point cloud, the input point cloud is symmetric with respect to the cross-sectional plane formed by the x-axis and z-axis. Consequently, the fused point cloud should also maintain symmetry with respect to the cross-sectional plane formed by the x-axis and z-axis. Thus, the fused point cloud should satisfy the constraints shown in the following equation.

$$[P_{i,x}^{left},P_{i,y}^{left},P_{i,z}^{left}]=[P_{i,x}^{right},-P_{i,y}^{right},P_{i,z}^{right}]$$

(5)

where $[P_{i,x}^{left},P_{i,y}^{left},P_{i,z}^{left}]$ represents the fused point cloud coordinates for y < 0, and $[P_{i,x}^{right},-P_{i,y}^{right},P_{i,z}^{right}]$ represents the fused point cloud coordinates for y > 0.

2.1.6. Point Cloud Fusion Problem Based on EMD Auto-Evolution

Combining the above aspects of geometric structure information inheritance, auto-evolution of point cloud, weight adaptation for outlier point cloud, and point cloud symmetry constraint, the point cloud fusion problem based on EMD auto-evolution is formulated as follows:

$$\begin{array}{l}F={\displaystyle \sum _{i=1}^n{c}_i^{a,left}}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,left}}\cdot {‖[P_{i,x}^{a,left},P_{i,y}^{a,left},P_{i,z}^{a,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2\\ +\lambda \times {\displaystyle \sum _{i=1}^m{c}_i^{b,left}}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,left}}\cdot {‖[P_{i,x}^{b,left},P_{i,y}^{b,left},P_{i,z}^{b,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2\\ +{\displaystyle \sum _{i=n+1}^{2n}{c}_i^{a,right}}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,right}}\cdot {‖[P_{i,x}^{a,right},P_{i,y}^{a,right},P_{i,z}^{a,right}]-[P_{j,x}^{right},P_{j,y}^{right},P_{j,z}^{right}]‖}_2\\ +\lambda \times {\displaystyle \sum _{i=m+1}^{2m}{c}_i^{b,right}}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,right}}\cdot {‖[P_{i,x}^{b,right},P_{i,y}^{b,right},P_{i,z}^{b,right}]-[P_{j,x}^{right},P_{j,y}^{right},P_{j,z}^{right}]‖}_2\end{array}$$

(6)

where $s^{a,left}$, $s^{a,right}$, $s^{b,left}$, and $s^{b,right}$, respectively, represent the correspondence between the left half (i.e., y < 0) of the fused point cloud and the input point cloud group a, the right half (i.e., y > 0) of the fused point cloud and the input point cloud group a, the left half (i.e., y < 0) of the fused point cloud and the input point cloud group b, and the right half (i.e., y > 0) of the fused point cloud and the input point cloud group b.

This formula comprehensively encapsulates the four mechanisms mentioned above. As a result, the fused point cloud not only fully inherits the geometric structure of the input point cloud but also effectively prevents the adverse effects of outlier point cloud on the fusion process while preserving the inherent symmetric structure of the target point cloud. By solving this equation, the coordinates of the fused point cloud can be obtained.

Since solving Equation (6) is a non-convex optimization problem, an alternating iterative technique is employed to solve the equation effectively by sequentially updating the variable $s^a,s^b,P,c^a,c^b$.

First, the variables $s^a$ and $s^b$ are computed while keeping the other variables fixed. At this stage, the fusion problem is transformed into a linear programming problem, which can be solved using the typical interior point method or simplex algorithm.

Next, the variable P is calculated while keeping the remaining variables fixed. The fusion problem at this stage is represented as shown in Equation (7).

$$\begin{array}{l}\underset{P}{\min }{\displaystyle \sum _{i=1}^n{c}_i^a}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,left}}\cdot {‖[P_{i,x}^{a,left},P_{i,y}^{a,left},P_{i,z}^{a,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2\\ +\lambda \times {\displaystyle \sum _{i=1}^m{c}_i^b}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,left}}\cdot {‖[P_{i,x}^{b,left},P_{i,y}^{b,left},P_{i,z}^{b,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2\\ +{\displaystyle \sum _{i=n+1}^{2n}{c}_i^a}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,right}}\cdot {‖[P_{i,x}^{a,right},P_{i,y}^{a,right},P_{i,z}^{a,right}]-[P_{j,x}^{left},-P_{j,y}^{left},P_{j,z}^{left}]‖}_2\\ +\lambda \times {\displaystyle \sum _{i=m+1}^{2m}{c}_i^b}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,right}}\cdot {‖[P_{i,x}^{b,right},P_{i,y}^{b,right},P_{i,z}^{b,right}]-[P_{j,x}^{left},-P_{j,y}^{left},P_{j,z}^{left}]‖}_2\end{array}$$

(7)

The above equation represents an unconstrained quadratic optimization problem, whose optimal solution can be directly obtained by setting the gradient of the objective function with respect to the variable to zero. The expression for the optimal solution of the variable P is given as follows:

$$\begin{array}{l}[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]=\\ ({\displaystyle \frac{1}{{\displaystyle \sum _{i=1}^n{c}_i^a}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,left}}+\lambda \times {\displaystyle \sum _{i=1}^m{c}_i^b}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,left}}+{\displaystyle \sum _{i=n+1}^{2n}{c}_i^a}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,right}}+\lambda \times {\displaystyle \sum _{i=m+1}^{2m}{c}_i^b}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,right}}}})\\ \times ({\displaystyle \sum _{i=1}^n{c}_i^a}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,left}}\cdot [P_{i,x}^{a,left},P_{i,y}^{a,left},P_{i,z}^{a,left}]+\lambda \times {\displaystyle \sum _{i=1}^m{c}_i^b}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,left}}\cdot [P_{i,x}^{b,left},P_{i,y}^{b,left},P_{i,z}^{b,left}]\\ +{\displaystyle \sum _{i=n+1}^{2n}{c}_i^a}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,right}}\cdot [P_{i,x}^{a,right},-P_{i,y}^{a,right},P_{i,z}^{a,right}]+\lambda \times {\displaystyle \sum _{i=m+1}^{2m}{c}_i^b}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,right}}\cdot [P_{i,x}^{b,right},-P_{i,y}^{b,right},P_{i,z}^{b,right}])\end{array}$$

(8)

After calculating the point cloud coordinates according to Equation (8), the variables $c^a$ and $c^b$ are updated while keeping the other variables fixed. At this point, the fusion problem has been transformed into Equation (9).

$$\begin{array}{l}\underset{c^{a,left},c^{b,left},c^{a,right},c^{b,right}}{\min }{\displaystyle \sum _{i=1}^n{c}_i^{a,left}}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,left}}\cdot {‖[P_{i,x}^{a,left},P_{i,y}^{a,left},P_{i,z}^{a,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2\\ +\lambda \times {\displaystyle \sum _{i=1}^m{c}_i^{b,left}}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,left}}\cdot {‖[P_{i,x}^{b,left},P_{i,y}^{b,left},P_{i,z}^{b,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2\\ +{\displaystyle \sum _{i=n+1}^{2n}{c}_i^{a,right}}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,right}}\cdot {‖[P_{i,x}^{a,right},P_{i,y}^{a,right},P_{i,z}^{a,right}]-[P_{j,x}^{left},-P_{j,y}^{left},P_{j,z}^{left}]‖}_2\\ +\lambda \times {\displaystyle \sum _{i=m+1}^{2m}{c}_i^{b,right}}{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,right}}\cdot {‖[P_{i,x}^{b,right},P_{i,y}^{b,right},P_{i,z}^{b,right}]-[P_{j,x}^{left},-P_{j,y}^{left},P_{j,z}^{left}]‖}_2\end{array}$$

(9)

where $c^a=[c^{a,left},c^{a,right}]$ and $c^b=[c^{b,left},c^{b,right}]$. The fusion problem in Equation (9) can be solved using the Lagrangian dual function and the Karush–Kuhn–Tucker (KKT) conditions. The expression for the optimal solution is given in Equation (10).

$$\begin{array}{l}{c}_i^{a,left}={\displaystyle \frac{{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,left}}\cdot {‖[P_{i,x}^{a,left},P_{i,y}^{a,left},P_{i,z}^{a,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2}{\sqrt{{\displaystyle \sum _{i=1}^n{\left({\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,left}}\cdot {‖[P_{i,x}^{a,left},P_{i,y}^{a,left},P_{i,z}^{a,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2\right)}^2}}}}\\ c_i^{b,left}={\displaystyle \frac{{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,left}}\cdot {‖[P_{i,x}^{b,left},P_{i,y}^{b,left},P_{i,z}^{b,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2}{\sqrt{{\displaystyle \sum _{i=1}^m{\left({\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,left}}\cdot {‖[P_{i,x}^{b,left},P_{i,y}^{b,left},P_{i,z}^{b,left}]-[P_{j,x}^{left},P_{j,y}^{left},P_{j,z}^{left}]‖}_2\right)}^2}}}}\\ c_i^{a,right}={\displaystyle \frac{{\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,right}}\cdot {‖[P_{i,x}^{a,right},P_{i,y}^{a,right},P_{i,z}^{a,right}]-[P_{j,x}^{left},-P_{j,y}^{left},P_{j,z}^{left}]‖}_2}{\sqrt{{\displaystyle \sum _{i=1}^n{\left({\displaystyle \sum _{j=1}^k{s}_{i,j}^{a,right}}\cdot {‖[P_{i,x}^{a,left},P_{i,y}^{a,right},P_{i,z}^{a,right}]-[P_{j,x}^{left},-P_{j,y}^{left},P_{j,z}^{left}]‖}_2\right)}^2}}}}\\ c_i^{b,right}={\displaystyle \frac{{\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,right}}\cdot {‖[P_{i,x}^{b,right},P_{i,y}^{b,right},P_{i,z}^{b,right}]-[P_{j,x}^{left},-P_{j,y}^{left},P_{j,z}^{left}]‖}_2}{\sqrt{{\displaystyle \sum _{i=1}^m{\left({\displaystyle \sum _{j=1}^k{s}_{i,j}^{b,right}}\cdot {‖[P_{i,x}^{b,right},P_{i,y}^{b,right},P_{i,z}^{b,right}]-[P_{j,x}^{left},-P_{j,y}^{left},P_{j,z}^{left}]‖}_2\right)}^2}}}}\end{array}$$

(10)

By alternately iterating and updating the variables, the point cloud fusion algorithm terminates when the number of iterations exceeds a predefined threshold or when the value of the objective function converges earlier. The optimal fused point cloud coordinates P are then output.

$${‖F_t(s^a,s^b,c^a,c^b,P)-F_{t-1}(s^a,s^b,c^a,c^b,P)‖}_2\le \epsilon$$

(11)

For asymmetric targets, the point cloud symmetry constraint can be ignored. The variable $s^a,s^b,P,c^a,c^b$ is still updated sequentially, and the algorithm terminates when the number of iterations exceeds the threshold or when the objective function value converges early, outputting the optimal fused point cloud coordinates P.

2.2. Point Cloud Optimization Based on Local Parameterized Network

The aforementioned method enables the fusion of two sets of point cloud. However, the point cloud obtained before fusion are derived from SSIs and inherently exhibit issues such as sparsity. As a result, the preliminary point cloud obtained through EMD auto-evolution still suffers from low density, losing many points that effectively describe the target’s features. Additionally, as a simplified representation of 3D data, point cloud have limitations in visual quality compared to refined 3D products like meshes. Therefore, it is necessary to optimize the initial fused point cloud obtained in the previous step, enhancing its density and constructing mesh connections.

For targets with simple structures, local geometric optimization methods [21,22,23] can achieve good results on smooth surfaces with few features. However, for more complex targets such as planes, ships, and buildings, the unstructured and irregular nature of point cloud makes direct point cloud upsampling unfeasible. Currently, common approaches involve extracting point cloud features and performing upsampling at the feature level [24,25,26]. Although these methods can produce better results, they are largely influenced by image processing and rarely consider the target’s 3D geometric shape. As a result, various artifacts can be observed in their outputs. Moreover, these methods cannot directly obtain the mesh of optimized point cloud; a mesh model must be constructed using other mesh generation methods, which can introduce additional errors.

To address these issues, this paper proposes a point cloud optimization method based on local parameterization network. This approach enhances point cloud density while ensuring that the points lost during the initial fusion are largely recovered, enabling normal information extraction for mesh construction.

Given a point cloud $P=[P_1,\cdots ,P_i,\cdots ,P_k]$ with k points obtained from EMD auto-evolution evolution fusion and a specified upsampling factor R. Our goal is to generate a dense, uniformly distributed point cloud $P^R=[P_1,\cdots ,P_{1+R},\cdots ,P_i,\cdots ,P_{i+R},\cdots ,P_k,\cdots ,P_{k+R}]$ with corresponding normals for each point. This point cloud should contain more geometric details and as many points from $P^a$ and $P^b$ as possible. Our method is based on local parameterization network, comprising four steps: parameterizing the 3D point cloud into a 2D domain, extracting and connecting multi-scale hierarchical point wise features of the 3D point cloud, upsampling the point cloud in the 2D parameter domain, and mapping the 2D samples back to 3D space through a linear transformation. Local parameterization involves dividing the input sparse point cloud into multiple paths using a farthest point sampling algorithm. The fundamental theorem of surface local theory states that the local neighborhood of points on a regular surface can be fully determined by the first fundamental form and the second fundamental form [27]. Therefore, our key idea is to learn local parameterization for each point rather than computing and learning expensive global parameterizations [28].

Firstly, the local neighborhood of a point $P_i$ in 3D space is parameterized into a 2D domain using the differential mapping $\varphi :{ℝ}^2\to {ℝ}^3$, where $\varphi (0,0)=P_i$ and $(u,v)$ represent a point in the 2D domain (as shown in Figure 2). The Jacobian matrix $J_{\varphi }=[{\varphi }_u,{\varphi }_v]$ provides the best first-order approximation $\varphi (u,v)=\varphi (0,0)+[{\varphi }_u,{\varphi }_v]\cdot {(u,v)}^{\mathrm{T}}+O(u^2+v^2)$ of the mapping $\varphi$, where ${\varphi }_u$ and ${\varphi }_v$ are the tangent vectors defining the first fundamental form [27]. The normal at point $P_i$ can be calculated using the cross product $N_i={\varphi }_u(u,v)\times {\varphi }_v(u,v)={\varphi }_u(0,0)\times {\varphi }_v(0,0)$. For any point $\widehat{P}=P_i+J_{\varphi }\cdot {(u_1,v_1)}^{\mathrm{T}}$, since $(\widehat{P}-P_i)\cdot N_i=0$, it can be proven that point $\widehat{P}$ lies on the tangent plane of point $P_i$. This allows a point $(u_1,v_1)$ in the 2D domain to be mapped to a point $\widehat{P}$ on the tangent plane of a point $P_i$ in 3D space. Furthermore, we use the augmented Jacobian matrix $\widehat{J}=[{\varphi }_u,{\varphi }_v,{\varphi }_u\times {\varphi }_v]$ to calculate the normal $N_i=\widehat{J}\cdot {(0,0,1)}^{\mathrm{T}}$ and point $\widehat{P}=P_i+\widehat{J}\cdot {(u_1,v_1,0)}^{\mathrm{T}}$. Finally, the distance between a point $\widehat{P}$ on the tangent plane and its corresponding point $P$ in 3D is given by $\delta =‖P-\widehat{P}‖={\displaystyle \frac{{\kappa }_1{u}_1^2+{\kappa }_2{v}_1^2}{2}}+O(u_1^3,v_1^3)$, where ${\kappa }_1$ and ${\kappa }_2$ are the principal curvatures at point $\varphi (0,0)$, which are the eigenvalues of the second fundamental form. With this, we can reconstruct the local geometric relationship of point $P_i$ in the 2D parameter domain. In other words, based on the local neighborhood parameterization of $P_i$, any point $(u,v)$ in the 2D space can be correspondingly mapped to a point in the 3D space.

Specifically, as shown in Figure 3, given a local path $P=[P_1,\cdots ,P_i,\cdots ,P_k]$ of a 3D point cloud, the method first projects the local 3D point cloud onto a 2D parameter domain. In this 2D domain, new points $u=[u_i^1,\cdots ,u_i^r,\cdots ,u_i^R]$ and $v=[v_i^1,\cdots ,v_i^r,\cdots ,v_i^R]$ are generated for each point in the path. The normal $N_i={\widehat{J}}_i\cdot {(0,0,1)}^{\mathrm{T}}$ of the path is then calculated. The newly generated points in the 2D parameter domain are subsequently mapped onto the tangent plane of $P_i$ using ${\widehat{P}}_i^r=P_i+{\widehat{J}}_i\cdot {(u_i^r,v_i^r,0)}^{\mathrm{T}}$. Finally, by calculating the displacement distance ${\delta }_i^r$ along the normal direction, the points ${\widehat{P}}_i^r$ are projected back onto the object in 3D space.

2.2.1. Multi Scale Hierarchical Feature Extraction

We apply DGCNN [29] to extract hierarchical point-wise features, which can encode both local and global intrinsic geometric information of the input blocks. The hierarchical feature learning module extracts features from low to high levels. Intuitively, as the receptive field increases, skip connections [30], which are widely used in 2D vision tasks to improve feature quality and convergence speed, help retain details at different levels. Additionally, instead of direct feature concatenation, feature reweighting through a self-gating unit is applied to enhance the features.

Let $f_i^l$ denote the features extracted for point $P_i$ at level l ($l=1,\cdots ,L$). We first concatenate all features across L layers, i.e., $\widehat{f_i}=Concat(f_i^1,\cdots ,f_i^L)$, where $Concat(\cdot )$ represents the concatenation operator. The directly concatenated features are then fed into a small multi-layer perceptron (MLP) $h_r(\cdot )$ to obtain the logits $logi{c}_i=(logi{c}_i^1,\cdots ,logi{c}_i^L)$, as follows:

$$logi{c}_i=h_r(\widehat{f_i})$$

(12)

These weights are further fed into a softmax layer to produce new weighted coefficients $w_i=(w_i^1,\cdots ,w_i^L)$, as given by:

$$w_i^l={\displaystyle \frac{e^{f_i^l}}{{\displaystyle \sum _{j=1}^L{e}^{f_i^j}}}}$$

(13)

Finally, the multi-scale features are represented as a weighted concatenation:

$$f_i=Concat(w_i^1\cdot c_i^1,\cdots ,w_i^L\cdot c_i^L)$$

(14)

2.2.2. Point Cloud Upsampling Based on 2D Parameter Plane

In this module, we expand the input point cloud by R times. By regressing the obtained multi-scale features, a coarse dense point cloud and the corresponding coarse normals are generated. Specifically, the expansion process consists of two steps: first, adaptive upsampling is learned in the 2D parameter domain; then, the results are projected into the 3D tangent space through a learned linear transformation. The process is illustrated in Figure 4.

For each point $P_i$, we apply an MLP $g_1(\cdot )$ to its local surface features $f_i$ to generate R new 2D coordinate points $(u_i^r,v_i^r)$, as follows:

$$(u_i^r,v_i^r)=g_1(f_i)$$

(15)

where $r=1,\cdots ,R$. Furthermore, to fully utilize the information from the two pre-fusion point cloud $P^a$ and $P^b$ discussed in the previous section, the generation of the 2D coordinate points incorporates the corresponding 2D coordinates from $P^a$ and $P^b$ based on the adaptive weighting factors $c^a$ and $c^b$.

This approach allows the incorporation of points from the original point cloud before fusion that contribute to refining the target shape into the fusion result while ensuring the new 2D points are generated uniformly in the 2D domain. However, due to the transformation relationship between the 2D parameter domain and the 3D point cloud space, these newly upsampled points are non-uniformly distributed on the 3D surface, which more closely aligns with real-world scenarios.

For each point $P_i$ in the 2D parameter domain, the local surface features $f_i$ are used to predict a linear transformation matrix $\widehat{J}$, given by:

$$\widehat{J}=g_2(f_i)$$

(16)

where $g_2(\cdot )$ denotes the MLP. The $\widehat{J}$ is multiplied by the 2D coordinates $(u_i^r,v_i^r)$ to map the points onto the tangent plane of $P_i$, resulting in a coarse point cloud. These points lie on the tangent plane rather than the target surface.

$${\widehat{P}}_i^r=P_i+{\widehat{J}}_i\cdot {(u_i^r,v_i^r,0)}^{\mathrm{T}}$$

(17)

At the same time, we also estimate a coarse normal, specifically the normal $N_i$ of the tangent plane for each input point. To achieve this, we multiply the linear transformation matrix $\widehat{J}$ by a predefined normal (0, 0, 1) that is perpendicular to the 2D parameter domain:

$$N_i=\widehat{J}\cdot {(0,0,1)}^T$$

(18)

2.2.3. Mapping 2D Samples Back to 3D Space

Since the locally expanded point cloud on the 2D parameter plane is located on the tangent plane, it needs to be moved onto the target surface (typically a curved surface) and its normals updated. Specifically, each point ${\widehat{P}}_i^r$ in the point cloud is displaced along its normal $N_i$ by a distance ${\delta }_i^r$. This distance ${\delta }_i^r$ is calculated by regressing the features of each point connected with the coarse coordinates, as follows:

$${\delta }_i^r=g_3(Concat({\widehat{P}}_i^r,f_i))$$

(19)

where $g_3(\cdot )$ denotes the process of the MLP (multi-layer perceptron). The updated coordinates are then computed as:

$$P_i^r={\widehat{P}}_i^r+{\widehat{J}}_i\cdot {(0,0,{\delta }_i^r)}^{\mathrm{T}}$$

(20)

We update the normals in a similar manner, and the normal displacement $\Delta N_i^r$ of point $P_i^r$ is regressed as:

$$\Delta N_i^r=g_4(Concat({\widehat{P}}_i^r,f_i))$$

(21)

This value is further added to the corresponding coarse normal to obtain the refined normal:

$$N_i^r=\Delta N_i^r+N_i$$

(22)

where $g_4(\cdot )$ represents the MLP process.

2.2.4. Losses

In terms of loss, we designed a combined loss function to train the model end-to-end. Specifically, $Q^D$ is formed by selecting D points from the corresponding target high-density 3D point cloud ground truth, serving as the ground truth for $P^R$. During training, we use the Chamfer distance (CD) to measure the coordinate error between $P^R$ and $Q^D$, which is defined as:

$$L_{CD}={\displaystyle \frac{1}{D}}\left({\displaystyle \sum _{P_i^r\in P^R}{‖P_i^r-\eta (P_i^r)‖}_2}+{\displaystyle \sum _{Q_d\in Q^D}{‖Q_d-\psi (Q_d)‖}_2}\right)$$

(23)

where $\psi (Q_d)=argmi{n}_{P_i^r\in P^R}{‖P_i^r-Q_d‖}_2$ and $\psi (Q_d)=argmi{n}_{P_i^r\in P^R}{‖P_i^r-Q_d‖}_2$.

For the normals, the ground truth normals of points $P_i$ in the 3D point cloud are denoted as $\widetilde{N_i}$, and the normals of points $P^R$ are denoted as $\widehat{N_i}$. These are considered as the ground truths for the coarse normals $N_i$ and refined normals $N_R$, respectively. During training, we consider the error between both the coarse and refined normals:

$$L_{COARSE}={\displaystyle \sum _{i=1}^{D}L(N_i,\widetilde{N_i})}$$

(24)

$$L_{REFINED}={\displaystyle \sum _{i=1}^D{\displaystyle \sum _{r=1}^{R}L(N_i^r,\widehat{N_i})}}$$

(25)

where $L(N_i,\widetilde{N_i})=max\{{‖N_i-\widetilde{N_i}‖}_2^2,{‖N_i+\widetilde{N_i}‖}_2^2\}$ measures the unoriented difference between two normals. Finally, the combined loss function can be expressed as:

$$L_{ALL}={\lambda }_1{L}_{CD}+{\lambda }_2{L}_{COARSE}+{\lambda }_3{L}_{REFINED}$$

(26)

where ${\lambda }_1,{\lambda }_2,{\lambda }_3$ represents three positive parameters.

3. Experimental Results and Discussion

3.1. Data Sets, Metrics, and Implementation Details

• Dataset

The testing dataset A and all training data used in this paper were self-constructed through purchases and internet collection. The 3D data used for fusion consist of 3D point cloud data reconstructed from SSIs using a target 3D reconstruction algorithm [3]. The SSIs were obtained from the SuperView-1 and Jilin-1 satellite, covering four regions: China (Harbin), Japan, the United States, and Qatar. The 3D ground truth models for evaluating the fusion results were acquired through UAV oblique photogrammetry, Google Earth 3D, and internet collection.

Since there is no publicly available dataset for satellite-based 3D data fusion, we indirectly validated the fusion performance of our proposed method using a set of publicly available datasets for stereo matching and DSM fusion. Dataset B is a public dataset from the 2019 IGRASS Data Fusion Contest and contains World-View satellite data with a resolution of 0.35 m and DSM ground truth values obtained by LiDAR.

• Metrics

The experiment uses the EMD to evaluate the similarity of the fusion result to the actual target shape. A smaller EMD value indicates that the fusion result is closer to the true shape. For both the real 3D point cloud model and the fused 3D point cloud, the length between any two points in the point cloud is calculated, and half of the maximum length is used as the normalization radius. The center point coordinates of these two points are taken as the sphere center coordinates. All coordinates of the point cloud are then subtracted by the sphere center coordinates and divided by the normalization radius to obtain the normalized 3D point cloud. The EMD is then calculated as follows:

$$EMD=\underset{\phi :{\widehat{P}}^R\to {\widehat{P}}^{truth}}{\min }{\displaystyle \sum _{{\widehat{P}}^R}\frac{1}{k\times R}{‖{\widehat{P}}^R-{\widehat{P}}^{truth}‖}_2}$$

(27)

where $P^R$ represents the normalized fused point cloud, $P^{truth}$ represents the normalized point cloud ground truth, $\phi$ represents the unidirectional mapping formed by each point pair between $P^R$ and $P^{truth}$, and $k\times R$ denotes the number of points in the fused point cloud.

While EMD can only evaluate the accuracy of the target shape, the 3D data of ground targets obtained from satellite imagery in this study have clear physical attributes. Therefore, in addition to EMD, the root mean square error metric is added to describe the spatial positional accuracy of the fusion result.

The calculation of the RMSE metric requires selecting all N points from the fused 3D point cloud as test points and comparing them with the corresponding points in the ground truth. The calculation method is shown below:

$$RMSE={\displaystyle \frac{1}{\sqrt{N}}}{‖{\tilde{P}}^R-{\tilde{P}}^{truth}‖}_2$$

(28)

where ${\tilde{P}}^R$ and ${\tilde{P}}^{truth}$ represent the coordinates of a total of N test points in the fused point cloud and their corresponding nearest neighbor points in the ground truth point cloud, respectively. When calculating point cloud accuracy, P is the 3D coordinate. When calculating DSM accuracy, P is the elevation value.

• Implementation Details

The experiments were conducted on a PC equipped with an Intel Core i7-10870H CPU, 16GB RAM, and an Nvidia RTX 2080 GPU, without using any parallel computing programs or other specialized hardware. The values of the three parameters in the combined loss function were empirically set to 100, 1, and 1, respectively. In our experiments, we focused exclusively on three types of targets: buildings, planes, and ships.

3.2. Results and Comparisons

3.2.1. Dataset A

To comprehensively demonstrate the effectiveness of the proposed algorithm, several recent algorithms similar to the proposed method were selected for comparison experiments, including two traditional methods (ICP [4] and 3DT [11]) and two deep learning methods (DCP [17] and MS3DQE-Net [16]). For ICP and DCP, the fusion results were obtained by directly merging the registration 3D point cloud. The parameters for all comparison methods were set according to the recommendations in their respective papers, and the network models were implemented using their original architectures. Some experimental results are shown in Figure 5. Since some of the comparison methods are specifically designed for point cloud and this paper only evaluates point cloud accuracy, and considering that 3D models in mesh form are not convenient for observing point cloud positions, only the point cloud are displayed in the result images for visual clarity.

From the experimental results shown in the figures, it is evident that the proposed method achieves better fusion performance on targets such as planes, ships, and buildings, outperforming other comparison methods in terms of point cloud density and target completeness. The quantitative evaluation in

Table 1. Comparison to different 3D point fusion methods in dataset A.

Method	EMD			RMSE (m)
	Plane	Ship	Building	Plane	Ship	Building
3D point cloud reconstructed from SSIs [3]	0.14	0.48	0.27	1.69	2.94	1.71
ICP	0.12	0.43	0.23	1.58	2.70	1.58
3DT	0.13	0.46	0.23	1.62	2.81	1.64
DCP	0.11	0.41	0.23	1.57	2.69	1.59
MS3DQE-Net	0.11	0.45	0.24	1.49	2.86	1.57
Proposed	0.09	0.39	0.21	1.42	2.24	1.50

shows that the proposed method achieves higher accuracy in both EMD and RMSE compared to other methods, indicating that the fusion results are closer to the real target in both spatial position and shape.

Figure 6 shows an enlarged version of some results from Figure 5. The first two sets of data in Figure 5 demonstrate the fusion results for two plane targets. Due to the relatively simple structure of plane targets, all methods achieved satisfactory fusion results. However, the proposed method shows superior fusion performance in the wing and tail areas, with fewer noise points, as shown by the black circle in the first set of results in Figure 6. The third and fourth sets of data illustrate the fusion results for two ship targets. Ship targets are large and structurally complex, particularly the types of ships discussed in this paper, which have a large and prominent bridge structure. Limited by the resolution of SSIs, it is difficult to completely reconstruct the overall structure of the ships using only satellite data, resulting in relatively low-quality source data for fusion. As seen by the black circle in Figure 6, the proposed method ensures that the bridge structure remains prominently raised after fusion, while some methods lose this critical detail. The last two sets of data show the fusion results for two building targets. The proposed method achieves good fusion results in flat areas such as rooftops and sides of buildings, with a significant reduction in noise points. Additionally, it effectively reconstructs unique roof structures and restores the shapes of ornamental structures more completely than other methods. This can be seen from red and black circles in the third set of results in Figure 6.

Regarding the comparison methods involved in the experiments, both ICP and DCP, despite using different approaches, fundamentally aim to register two sets of point cloud data to achieve overlapping fusion, resulting in similar accuracy in the experiments. The 3DT method was primarily proposed to handle 3D data obtained in urban areas with complementary information through different means. However, the 3D data in this paper come from the same data source, even though they are obtained from remote sensing images taken at different times. There is relatively little complementary information in the source data, which does not meet the ideal application scenario for 3DT. MS3DQE-Net is mainly designed for the fusion of LiDAR point cloud data in large-scale areas, and it is prone to errors when dealing with single targets where target information is incomplete.

3.2.2. Dataset B

Since there is no publicly available dataset for satellite-based 3D data fusion, we indirectly validated the fusion performance of our proposed method using a set of publicly available datasets for stereo matching and DSM fusion. Specifically, as shown in Figure 7, we perform stereo matching and coordinate calculation on two sets of stereo image pairs to obtain DSMs. Then, we extract the 3D models of building targets from each DSM and fuse the point clouds. The fused result is reprojected into a DSM and compared with the DSM ground truth provided in the dataset. This method allows us to indirectly validate the effectiveness of the fusion method by evaluating the accuracy of the fused point cloud in the height dimension. The stereo matching and coordinate calculation are achieved through the method proposed in reference [31].

The comparison methods used in the experiments with dataset B are the same as those for dataset A. The experimental results are shown in Figure 8, and the DSM accuracy is presented in

Table 2. Comparison to different 3D point fusion methods in dataset B.

Method	RMSE
DSM from SSIs [31]	2.66
ICP	2.39
3DT	2.42
DCP	2.25
MS3DQE-Net	2.18
Proposed	2.07

. In Figure 8, it can be observed that the proposed method effectively fuses the auxiliary structures on top of the buildings, and as indicated by the red circles, it is closer to the ground truth. In Table 2, the proposed method achieves better accuracy compared to other methods. Although the RMSE for dataset B is calculated only in the elevation dimension, its performance is worse than in dataset A. This is because the DSM in dataset B was derived from [31], and its data quality cannot be compared to that of dataset A, which was obtained from [3]. This also demonstrates that the proposed method can achieve good results even when fusing low-quality data.

3.3. Ablation Study

To verify the effectiveness of the key modules in the proposed fusion method, an ablation study was conducted by altering or removing certain components to compare the fusion performance. The experimental results are shown in

Table 3. Comparison results with different key modules.

EMD Auto-Evolution	Multi Scale Hierarchical Feature Extraction	Point Cloud Upsampling Based on 2D Parameter Plane	Mapping 2D Samples Back to 3D Space	EMD	RMSE (m)
	✓	✓	✓	0.32	1.87
✓		✓	✓	0.29	1.72
✓	✓		✓	0.35	1.95
✓	✓	✓		0.31	1.79
✓	✓	✓	✓	0.25	1.68

. When there is no “EMD auto-evolution” module, two sets of 3D point clouds are directly merged through registration.

From Table 3, we can draw the following conclusions: (1) EMD auto-evolution significantly improves fusion performance compared to directly merging point cloud. (2) The absence of point cloud upsampling on the 2D parameter plane not only affects the upsampling accuracy and introduces errors but also fails to recover the point cloud from the source data, resulting in greater errors. (3) While other modules do not have a decisive impact on accuracy, they all contribute to the overall performance.

3.4. Robustness Analysis

3.4.1. Noisy Data

To verify the robustness of the algorithm, Gaussian noise at different noise levels was added to the point cloud obtained from SSIs, referred to as point cloud 1 and point cloud 2. The proposed fusion algorithm was then applied to fuse the two noisy point cloud. The fusion results are shown in the figure below.

As shown in Figure 9, despite the addition of noise to the point cloud, the fusion results still maintain a good shape. The proposed algorithm remains highly effective even with such challenging data, demonstrating its robustness to noise.

3.4.2. Fusion Data Obtained from Different Data Sources

Although the number of satellites is continuously increasing, it can sometimes be challenging to obtain multiple sets of satellite data for the same region within a short period. Therefore, supplementing with other data sources can effectively enhance fusion efficiency. To this end, the point cloud obtained from optical satellite remote sensing images were fused with a coarse 3D model of the target obtained from SAR data to evaluate the effectiveness of the proposed fusion method. The optical and SAR remote sensing images used are shown in Figure 10.

Achieving high precision 3D reconstruction of single targets using SAR data alone remains challenging. Therefore, this study extracts 3D information using a building height estimation method based on SAR layover to directly obtain the heights of flat and simple sloped roofs [32]. Subsequently, 3D models are manually generated using the roof heights and building outlines. This method is suitable only for a small number of regularly shaped buildings. Considering that these simple buildings are not easy to observe when displayed as point cloud, the buildings in Figure 11 are presented in mesh form.

The experimental results in Figure 11 and

Table 4. Accuracy of fusion results from optical and SAR data sources.

Data Source	EMD				RMSE (m)
	Building 1	Building 2	Building 3	Building 4	Building 1	Building 2	Building 3	Building 4
Optical SSIs data	0.13	0.14	0.14	0.17	1.75	1.82	1.84	1.65
SAR data	0.19	0.24	0.23	0.25	1.80	1.86	1.91	1.93
Fusion results	0.09	0.08	0.12	0.14	1.75	1.78	1.79	1.57

show that the fused building targets retain their inherent structural characteristics while providing more regular and higher confidence building structures, resulting in a significant improvement in accuracy. The proposed method effectively addresses the fusion of 3D data obtained from different data sources.

4. Conclusions

In this work, we propose a 3D point cloud fusion method based on EMD auto-evolution and local parameterization network. The method consists of two stages. The first stage addresses the fusion measurement problem by introducing EMD as a key metric for evaluating the fusion results. We construct an optimization problem for the point cloud fusion method based on EMD auto-evolution by incorporating input point cloud geometric structure inheritance, auto-evolution of the point cloud, adaptive point cloud weighting, and symmetry constraints. An alternating iterative technique is used to sequentially optimize the variables, producing an initial fusion result. This ensures that the fusion result more closely approximates the true target shape while allowing for the fusion of point cloud with non-uniform lengths. The second stage addresses the point cloud density issue. Based on the initial fusion results, we construct a local parameterization network and design a combined loss function to train it end to end. By parameterizing the initial fused point cloud into a 2D domain, extracting and connecting multi-scale hierarchical point-wise features of the 3D point cloud, and resampling in the 2D parameter domain before mapping back to 3D space via a linear transformation, we produce an optimized point cloud fusion result. This approach increases point cloud density while ensuring the recovery of more points lost during the initial fusion. Experimental results demonstrate that our fusion method outperforms other classical comparison algorithms for targets such as buildings, airplanes, and ships. The RMSE of the fusion results can reach approximately 2 m, and the EMD accuracy is better than 0.5. Nevertheless, the proposed method still has the following issues: a) The fusion accuracy depends on the quality of the original data before fusion. b) The fusion method does not take into account the physical characteristics of the target, lacking the use of structural information about the target. These issues will be explored in future research.