© Wuhan University 2026

0 Introduction

Point cloud data are widely used in numerous fields, such as intelligent grasping^[1], simultaneous localization and mapping (SLAM)^[2], 3D reconstruction^[3], and object recognition^[4]. Cloud registration^[5] is a very important step before reconstruction. However, due to continuous advancement of laser scanning technology in recent years, the amount of point cloud data has become increasingly large. Traditional registration algorithms often take a long time and are inefficient when processing large-scale point cloud data. To address this issue, point cloud data can be downsampled before registration.

Voxel grid filtering^[6], a standard downsampling protocol, partitions data into static cells, collapsing internal points to a single centroid. While computationally fast, this mechanism sacrifices critical geometric details in feature-rich regions. Curvature-based strategies^[7-8] target high-curvature features but fail to adapt to non-uniform densities. The flaw is systemic: they strip away essential key points in sparse regions while hoarding redundant data in dense areas, incurring unnecessary computational overhead. To mitigate such information loss, research shifted to virtual feature points^[9], a domain detailed in recent surveys^[10-11]. Nevertheless, standard methods often fail to account for non-uniform point densities, degrading robustness. Although Li et al^[9] introduced adaptive mechanisms, they disproportionately emphasize local geometric details while neglecting the impact of global density variations.

Registration pipelines generally follow a coarse-to-fine hierarchy: rough alignment establishes the baseline, followed by precise optimization. For the coarse stage, local feature matching is the dominant paradigm. Rusu et al^[12] introduced the Point Feature Histogram (PFH), which aligns scans by analyzing the histograms of value relationships between points. However, the $O(n\cdot k^{\mathrm{2}})$ complexity of PFH renders it impractical for real-time applications. The Fast Point Feature Histogram (FPFH)^[13] emerged to address this, slashing the complexity to $O(n\cdot k)$. Despite its speed advantage, standard FPFH relies on a fixed search radius. This approach oversamples dense regions, incurring computational redundancy, yet fails to extract key features in sparse areas. Alternatively, statistical methods utilizing global geometric distributions have been developed. Principal Component Analysis (PCA)-based algorithms^[14] achieve initial alignment by matching principal axes derived from covariance matrices, offering high computational efficiency. However, these global methods often perform poorly when facing partial overlaps or outlier interference. To address local feature description, 3D Shape Context (3DSC)^[15] captures the distribution of neighboring points using spherical bins, yet it remains computationally intensive and sensitive to point density variations. For fine registration, the Iterative Closest Point (ICP) algorithm^[16] remains the standard. However, it needs coarse registration to provide a good initial position. Otherwise, the registration accuracy is severely compromised.

To accelerate the registration of medium- and large-scale point clouds, an FPFH-based registration scheme integrating dynamic downsampling and an adaptive neighborhood approach is proposed. This involves an estimation of global point densities and a minimum feature threshold for dynamic point downsampling, achieving a dramatic reduction in the number of sampled points while maintaining local features in regions of low point density. Instead of employing a fixed search range, the size of the FPFH search boundary varies with point densities, resulting in enhanced efficiency of feature point calculations. Finally, the fine registration process is achieved via the ICP algorithm. Experimental results on two datasets demonstrate that the proposed algorithm improves efficiency with great registration accuracy.

1 Methodology

1.1 Adaptive Voxel Downsampling via Sampling Theorem

The traditional voxel grid filter^[6] represents the 3D space with cubic voxels and encodes the set of points in each voxel using the central point. Although useful for uniformly distributed point clouds, the approach is inadequate for datasets with non-uniformly distributed points. Fixed voxel sizes face a trade-off: large voxels eliminate details in sparse regions, while small voxels fail to sufficiently reduce redundancy in dense regions. Our scheme is based on a scale-adaptive approach described by the 3D sampling theorem, and the specific details are presented below.

1) Averaging the Euclidean distances among the k-nearest neighbors offers a precise estimate for the point density in the region. The global spacing measure, obtained by combining all local point density measures, provides a basic benchmark.

$D=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{d}_i=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}\left(\frac{\mathrm{1}}{k}{\displaystyle \sum _{j=\mathrm{1}}^k}\Vert p_i-p_{ij}\Vert \right)$(1)

where D denotes the global average point spacing, n represents the total number of points in the point cloud, and d_i is the local average distance for the i-th point p_i . Furthermore, k specifies the number of nearest neighbors, p_ij represents the j-th nearest neighbor of p_i , and $\Vert p_i-p_{ij}\Vert$ is the Euclidean distance between $p_i$ and $p_{ij}$ .The value sets a standard for density against which the downsampling process must compare.

2) The downsampling process can be represented by a model based on the expansion of an interval. The magnification factor m represents the ratio of the new interval to the original interval.

$m=\frac{v}{D}$(2)

Equation (2) binds m directly to the voxel size v. Consequently, deriving the optimal v is the central theoretical objective.

3) Through the formulation of the problem of retaining features, the Shannon-Nyquist theorem is used in connection with signal reconstruction to attain accuracy. The guiding principle states that error-free signal reconstruction must be able to capture the highest geometric frequency $f_{\mathrm{m}\mathrm{a}\mathrm{x}}$, which is determined by the minimum feature size $F_{\mathrm{m}\mathrm{i}\mathrm{n}}$.

$v\le \frac{\mathrm{1}}{\mathrm{2}{f}_{\mathrm{m}\mathrm{a}\mathrm{x}}}\approx \frac{F_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{2}}$(3)

4) Considering the non-ideal nature of real-world point clouds, a robustness coefficient $\alpha (\mathrm{0.5}\le \alpha \le \mathrm{0.7})$ is introduced to derive the maximum safe sampling interval:

$v_{\mathrm{m}\mathrm{a}\mathrm{x}}=\alpha \cdot \frac{F_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{2}}$(4)

5) By substituting the voxel size formula (4) into the magnification definition (2), we obtain the final theoretical expression for the adaptive $m$:

$m=\frac{v}{D}=\alpha \frac{F_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{2}D}$(5)

This formula achieves adaptive parameter setting for downsampling. The intensity of downsampling is not arbitrary; instead, it is objectively determined by the intrinsic properties of the point cloud (global point spacing $D$) and the extrinsic task requirements (minimum feature size $F_{\mathrm{m}\mathrm{i}\mathrm{n}}$).

Comparison of uniform voxel and Nyquist-based adaptive downsampling is shown in Fig.1.

Fig. 1 Comparison of uniform voxel and Nyquist-based adaptive downsampling

1.2 Scale-Adaptive FPFH Feature Extraction

To address the issue of low computational efficiency of traditional methods on non-uniform point clouds, this section first elucidates the calculation principle of FPFH, analyzes its bottleneck under the traditional fixed-radius strategy, and derives the acceleration mechanism of the proposed scale-adaptive method.

FPFH is an efficient operator describing the local geometric topology of point clouds. Its calculation process mainly includes the following three steps:

1) Local Frame Construction: For a query point $p_s$ and any point $p_t$ in its neighborhood, a local coordinate system (u, v, w) is established based on the normal vector $n_s$:

$\{\begin{array}{l}u=n_s,\\ v=u\times \frac{p_t-p_s}{\Vert p_t-p_s{\Vert }_{\mathrm{2}}},\\ w=u\times v,\end{array}$(6)

where $\Vert .\Vert$₂ denotes the L2 norm, and $\times$ represents the cross product. The resulting $(u,v,w)$ forms the local Darboux frame for the point pair.

2) Feature Component Calculation: In the local frame, three angular features $(\alpha ,\varphi ,\theta )$ between $p_s$ and $p_t$ are calculated to form a Simplified Point Feature Histogram (SPFH):

$\{\begin{array}{l}\alpha =v\cdot n_t,\\ \varphi =u\cdot \frac{p_t-p_s}{d},\\ \theta =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}(w\cdot n_t,u\cdot n_t).\end{array}$(7)

3) Neighborhood Weighted Fusion: To enhance feature descriptiveness, FPFH fuses the SPFH of the query point with those of its k-neighbors:

$\mathrm{F}\mathrm{P}\mathrm{F}\mathrm{H}(p_s)=\mathrm{S}\mathrm{P}\mathrm{F}\mathrm{H}(p_s)+\frac{\mathrm{1}}{k}{\displaystyle \sum _{i=\mathrm{1}}^k}\frac{\mathrm{1}}{{\omega }_i}\cdot \mathrm{S}\mathrm{P}\mathrm{F}\mathrm{H}(p_i),$(8)

where ${\omega }_i$ denotes the Euclidean distance between query point $p_s$ and its neighbor point $p_i$.

Figure 2 shows the FPFH feature extraction algorithm schematic diagram. Despite its effectiveness, the traditional FPFH algorithm typically sets a globally fixed search radius $r_{\mathrm{f}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{d}}$. In dense regions, the number of neighbor points $N_{\mathrm{n}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{b}\mathrm{o}\mathrm{r}}$ within this fixed radius grows cubically:

$N_{\mathrm{n}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{b}\mathrm{o}\mathrm{r}}\approx \frac{\mathrm{4}}{\mathrm{3}}\mathrm{\pi }{r}_{\mathrm{f}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{d}}^{\mathrm{3}}\cdot {\rho }_{\mathrm{m}\mathrm{a}\mathrm{x}}.$(9)

Fig. 2 FPFH feature extraction algorithm schematic diagram

Since the single-point calculation time of FPFH is proportional to the number of neighbors, the total computational complexity $T_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{d}}$ becomes extremely high:

$T_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{d}}\propto N\cdot {\rho }_{\mathrm{m}\mathrm{a}\mathrm{x}}\cdot r_{\mathrm{f}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{d}}^{\mathrm{3}},$(10)

where N denotes the total number of points in the dataset, and ${\rho }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ represents the maximum local point density.

Our method utilizes the optimal geometric description scale-voxel size v obtained from the preceding adaptive downsampling step. We dynamically bind the FPFH feature search radius r to the local voxel scale v:

$r_{\mathrm{a}\mathrm{d}\mathrm{a}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}=k_f\cdot v,$(11)

where $k_f$ is a constant scaling factor used to adjust the search radius. The introduction of $r_{\mathrm{a}\mathrm{d}\mathrm{a}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}$ ensures that the search volume scales proportionally with the point resolution, preventing an explosion in the number of neighbor points in high-density regions and thus maintaining high computational efficiency.

This dynamic binding leads to the following properties:

1) Density Normalization: After adaptive downsampling, the local density ${\rho }^{\text{'}}$ in any region is inversely proportional to the current voxel size v:

${\rho }^{\text{'}}\approx \frac{\mathrm{1}}{v^{\mathrm{3}}}.$(12)

2) Constant Neighbor Count: When calculating features, the number of neighbors $N_{\mathrm{n}\mathrm{e}\mathrm{w}}$ contained within the adaptive radius $r_{\mathrm{a}\mathrm{d}\mathrm{a}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}}$ is derived as:

$N_{\mathrm{n}\mathrm{e}\mathrm{w}}\approx V_{\mathrm{s}\mathrm{p}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}}\cdot {\rho }^{\text{'}}=\left(\frac{\mathrm{4}}{\mathrm{3}}\mathrm{\pi }(k_f{\cdot v)}^{\mathrm{3}}\right)\cdot \left(\frac{\mathrm{1}}{v^{\mathrm{3}}}\right),$(13)

where $V_{\mathrm{s}\mathrm{p}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}}$ denotes the volume of the adaptive search sphere, and ${\rho }^{\text{'}}$ represents the point density of the downsampled point cloud. The variable $v^{\mathrm{3}}$ in the numerator and denominator cancels out, yielding:

$N_{\mathrm{n}\mathrm{e}\mathrm{w}}\approx \frac{\mathrm{4}}{\mathrm{3}}\mathrm{\pi }{k}_f^{\mathrm{3}}=C.$(14)

This implies that the neighbor count $N_{\mathrm{n}\mathrm{e}\mathrm{w}}$ degenerates into a constant C related only to the coefficient $k_f$.

3) Theoretical Speedup Analysis: Since the neighbor count no longer varies with the original point cloud density $\rho$, the total time complexity $T_{\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{s}}$ of our method becomes $O(N^{\text{'}}\cdot C)$. The theoretical speedup ratio can be expressed as:

$\mathrm{S}\mathrm{p}\mathrm{e}\mathrm{e}\mathrm{d}\mathrm{u}\mathrm{p}=\frac{T_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{d}}}{T_{\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{s}}}\approx \frac{{\rho }_{\mathrm{m}\mathrm{a}\mathrm{x}}\cdot r_{\mathrm{f}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{d}}^{\mathrm{3}}}{\frac{\mathrm{4}}{\mathrm{3}}\mathrm{\pi }{k}_f^{\mathrm{3}}}\cdot \frac{N}{N^{\text{'}}},$(15)

where $T_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{d}}$ and $T_{\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{s}}$ represent the computational time of the traditional FPFH algorithm and our method, respectively. Here, ${\rho }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ denotes the maximum local point density of the original point cloud, and N/N′ represents the downsampling ratio.The analysis indicates that, without sacrificing feature extraction quality, the improved algorithm achieves a clear gain in computational efficiency.

Figure 3 shows comparison of neighborhood search strategy between traditional fixed-scale FPFH and our adaptive-scale FPFH.

Fig. 3 Comparison of neighborhood search strategy between traditional fixed-scale FPFH and our adaptive-scale FPFH

1.3 ICP Registration

Fine registration is sensitive to noise and initial pose, typically relying on a high overlap ratio and local feature consistency, thus requiring a relatively accurate initial registration. After coarse registration, the ICP ^[16] algorithm is employed for fine registration. The procedure is as follows:

1) For each point in the source point cloud, the closest corresponding point in terms of geometric distance will be searched for within the target point cloud space.

2) Outlier correspondences with excessively large distances or those that violate geometric constraints are filtered out to enhance matching robustness.

3) Based on the remaining correspondences, the optimal rigid transformation is computed using Singular Value Decomposition (SVD).

4) The current transformation matrix is applied to align the source point cloud with the target point cloud.

5) The above process is performed iteratively until the change in pose falls below a predefined threshold or the maximum number of iterations is reached. The error function is defined as follows:

$E(R,t)={\displaystyle \sum _{i=\mathrm{1}}^N}\Vert q_i-(R{p}_i+t){\Vert }^{\mathrm{2}},$(16)

where $E(R,t)$ represents the objective error function, R and t denote the rotation matrix and translation vector respectively, N is the total number of point correspondences, while $q_i$ and $p_i$ are the i-th pair of corresponding points in the target and source point clouds.

2 Experiment

2.1 Experimental Environment

The dataset used in this study was acquired using a laboratory-built high-precision three-dimensional scanning system based on line-structured laser technology, as shown in Fig. 4. A Keyence LJ-30 scanning head was used for scanning the workpiece, and a resolution in the range of a micrometer was obtained. For motion control, command signals were issued by the host computer to drive the XYZ stepper motors.

Fig. 4 Point cloud scanning system

To test the registration results of the proposed algorithm, experiments were conducted on two datasets of industrial blades. The first one represents a turbopump blade with dimensions of 78.83 mm × 105.50 mm× 23.46 mm and 1 024 897 points. The other set represents a damaged compressor blade that has larger dimensions of 114.42 mm× 109.36 mm× 42.60 mm with 1 282 190 points.

The two datasets exhibit distinct feature distributions. For the turbopump blade, many points are densely distributed over non-contour regions, resulting in substantial redundancy. If conventional downsampling strategies are applied directly, this redundancy leads to unnecessary computational overhead. In contrast, the point cloud of the damaged compressor blade shows pronounced geometric features at the edges and within internal cavity regions. Under a fixed downsampling scheme, critical geometric details in these areas are prone to being lost. To address the issue, the main body of the point cloud is treated as a low-frequency signal in this study. The minimum sampling frequency is set to 0.02, which effectively removes redundant points in relatively flat regions while allocating computational resources to areas with more complex geometric features. Combined with the computed global point spacing, this strategy is applied to downsample both point cloud datasets.

Both registration experiments were conducted in an identical computational environment. The computer used was running Windows 11 with an Intel Core i5-12600KF CPU and 32 GB of RAM. All the above-mentioned algorithms were developed in Python programming using the PyCharm environment.

2.2 Qualitative Registration Results

To intuitively evaluate the registration performance, we compared the proposed method with FPFH+ICP, 3DSC+ICP, and PCA+ICP. The visualization results are presented below, where the source point cloud is rendered in red and the target point cloud is blue.

Figure 5 illustrates the registration results for the engine turbopump blades. As shown in the magnified view (top-right), the registration quality varies significantly among the algorithms. Both 3DSC+ICP and PCA+ICP exhibit visible misalignment at the curved blade tip, resulting in a distinct "ghosting" effect where the red source points and blue target points are clearly separated. This indicates that these methods failed to converge accurately on the smooth surface. In contrast, FPFH+ICP and Ours+ICP achieve a seamless alignment. The source and target point clouds are perfectly superimposed, creating a uniform overlap. This visual evidence demonstrates that our method effectively handles the high point redundancy in these smooth regions and attains high precision comparable to the classic FPFH algorithm.

Fig. 5 Registration results of engine turbopump blades

Figure 6 presents the qualitative registration results for the damaged engine compressor blades. Unlike the smooth turbopump blades, this dataset represents a more challenging scenario characterized by severe geometric discontinuities and jagged fracture edges. As observed in the magnified view (bottom-right) focusing on the fracture zone, the registration performance differs markedly among the algorithms: The 3DSC+ICP method fails to preserve the sharp geometric structure at the boundary. Due to its sensitivity to non-uniform point density, the algorithm produces a chaotic point distribution, where the source and target clouds appear scattered and unable to converge on the edge features. Similarly, the PCA+ICP algorithm exhibits clear systematic errors. A visible translational offset exists between the red source contours and the blue target contours, indicating that the global principal component analysis was insufficient to handle the local deformations caused by the damage. In sharp contrast, Ours+ICP method demonstrates superior robustness. As shown in the magnified view, our method tightly aligns with complex, irregular fracture boundaries without any "ghosting" artifacts. By dynamically adjusting the feature extraction radius, our approach effectively captures the intricate details of the damaged regions, achieving a high-precision fit that visually outperforms the comparative methods.

Fig. 6 Registration results of damaged engine compressor blades

2.3 Quantitative Comparison and Analysis

To quantitatively evaluate the performance of the proposed algorithm, we compared it with FPFH+ICP, 3DSC+ICP, and PCA+ICP. The detailed comparison results regarding registration accuracy (Root Mean Square Error, RMSE) and computational time for both datasets are summarized in Table 1.

Accuracy Analysis: For the engine turbopump blades, the FPFH+ICP algorithm achieved the highest accuracy with an RMSE of $\mathrm{6.878}\times {\mathrm{10}}^{-\mathrm{3}}$ mm. This is expected, as the standard FPFH algorithm, with its fixed search radius, tends to capture broader contextual information in smooth regions, leading to stable feature descriptors. However, our method followed closely with an RMSE of $\mathrm{8.226}\times {\mathrm{10}}^{-\mathrm{3}}$ mm. The marginal difference indicates that our adaptive downsampling strategy successfully retained the essential geometric information while filtering out redundant data in flat regions. Conversely, for the damaged compressor blades, the performance trend reversed. Our method outperformed FPFH+ICP, achieving the lowest RMSE of $\mathrm{2.355}\times {\mathrm{10}}^{-\mathrm{3}}$ mm compared with that of FPFH ($\mathrm{3.294}\times {\mathrm{10}}^{-\mathrm{3}}$ mm). This superior accuracy is attributed to the adaptive neighborhood mechanism. While FPFH employs a fixed search radius that may under-sample critical features at the sparse, jagged fracture edges, our method dynamically adjusts the search scale based on local density. This ensures that the complex geometry of the damage is accurately described, leading to a more robust registration in non-uniform regions.

Efficiency Analysis: In terms of computational efficiency, our method demonstrates a significant advantage. For the large-scale damaged blade dataset, our method required only 19.42 s, reducing the processing time by approximately 55% compared with FPFH+ICP (43.25 s). Compared with the 3DSC+ICP method, our approach requires less than 10% of the processing time. This speedup validates the theoretical analysis in Section 1: by converting the neighbor search complexity from a cubic growth function to a constant scale via density normalization, we achieve high-speed processing without sacrificing accuracy.

Comparative Performance: Regarding the other baselines, PCA+ICP was the fastest (16.62 s) but suffers from poor accuracy (magnitude of ${\mathrm{10}}^{-\mathrm{1}}$ mm), proving that global features alone are insufficient for high-precision industrial inspection. 3DSC+ICP performed the worst in both metrics, likely due to its sensitivity to point density variations, making it unsuitable for scanning data with irregular distributions.

3 Conclusion

Aiming at the problem of low efficiency of existing algorithms in processing medium and large-scale point cloud registration, we propose an improved FPFH registration method based on dynamic downsampling and adaptive neighborhood mechanism. This method first calculates the global average point spacing based on the local density characteristics of the point cloud and then performs dynamic voxel downsampling on the original data to obtain a streamlined point cloud. Then, the neighborhood radius of the FPFH feature is dynamically adjusted according to the global average point spacing to enhance the adaptability of the feature descriptor to density changes. Based on the optimized feature set, the ICP algorithm is used to complete fine registration. Experiments show that the proposed algorithm significantly shortens the processing time while maintaining good accuracy. Contrary to fixed-radius approaches that failed at the sparse edges of the damaged compressor blade, our density-aware scheme successfully retained the critical geometric features, proving that the adaptive radius prevents information loss in non-uniform regions. Future work will extend this method to handle severe occlusion or integrate deep learning-based semantics.

References

All Tables

All Figures

	Fig. 1 Comparison of uniform voxel and Nyquist-based adaptive downsampling
In the text

	Fig. 2 FPFH feature extraction algorithm schematic diagram
In the text

	Fig. 3 Comparison of neighborhood search strategy between traditional fixed-scale FPFH and our adaptive-scale FPFH
In the text

	Fig. 4 Point cloud scanning system
In the text

	Fig. 5 Registration results of engine turbopump blades
In the text

	Fig. 6 Registration results of damaged engine compressor blades
In the text