Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Voxel-based spatial elongation filtering method for airborne single-photon LiDAR data

Open Access Open Access

Abstract

A novel voxel-based spatial elongation filtering method is proposed, to reduce noise in airborne single-photon lidar (SPL) data. In this method, six additional points are generated adjacent to each point of interest in the SPL data. Then, we count the number of points within each voxel and discriminate signals from noise via a predefined threshold. A filter performance evaluation index (taking into account the false alarm and signal loss rates, and the average distance between the residual noise points and their nearest signal points) is introduced. We compare the proposed and voxel-based spatial filtering method. The average false alarm rate found with our method (3.5%) is 18.6% smaller than that of the voxel-based spatial filtering method (4.3%).

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Single-photon lidar (SPL) is expected to usher in a new era of mapping, owing to its high data collection efficiency [1,2]. The laser beam of the SPL is split into multiple beamlets, and the weak echo signals are detected using a photon-counting detector array. Unfortunately, the high-sensitivity detectors used in SPL are susceptible to interference from the ambient light, leading to noise points in SPL detection results. Therefore, an effective filtering method is needed to eliminate noise and improve the quality of the SPL point clouds.

Many filtering methods have already been proposed. In early researches, point clouds were projected onto a two-dimensional (2D) plane, which consisted of flight and elevation directions, and 2D algorithms were then applied to remove noise from this plane [38]. However, the loss of one dimension of information makes this method incapable of eliminating noise near the signal. As a consequence, histogram-based filtering methods were instead put forward [911], these divide point clouds into bins at the expense of resolution. In 2016, A. Swatantran and H. Tang proposed a voxel-based spatial filtering method (VBSFM) [12,13]. This method splits the 3D space into massive voxels, and the number of points in the 27 voxels adjacent to the voxel of interest is compared to a predefined threshold, to distinguish whether the voxel is a signal voxel. However, the residual noise was found to be unacceptable, especially in regions with vegetation [14]. An adaptive ellipsoid searching filter [15] and voxel-spherical adaptive ellipsoid searching [16] methods were also proposed to reduce noise more effectively. However, these two methods analyze every single point in the SPL lidar data to distinguish whether it is a signal point or not, which leads to low computational efficiency and makes them unsuitable for large-scale de-noising [17].

In this study, a novel voxel-based spatial elongation filtering method (VBSEFM) is presented, to simultaneously reduce noise and preserve the signal. Before splitting the 3D space into massive voxels, for each point in the SPL data, six additional points are generated at different positions relative to the point of interest: above, below, left, right, in front, and behind. A filter performance evaluation index considering the false alarm rate, the signal loss rate, and the average distance between the residual noise points and their nearest signal points, is introduced to comprehensively evaluate the filtering effect. Since we do not know in advance which points are signal points in the actual SPL data, it is impossible to evaluate the filtering effect accurately. Thus, simulated SPL data are used to evaluate the filtering effect, and the performance evaluation indexes of VBSFM and our method, under various conditions, are calculated and compared.

2. Dataset

2.1 Original data

In this research, lidar point clouds – from streak tube imaging lidar (STIL) – with a high signal-to-noise ratio are specified as the original signal. In a manner similar to certain previous works [6,8,9], we add random Poisson-distribution noise to the original signal, to simulate the SPL point clouds. More information about the STIL can be found in Refs. [18,19]. Although there are some differences between the point clouds of STIL and SPL, there are also notable similarities. The average point density of the STIL is 6.5 points·m−2, which is comparable to that of the SPL data of the High-Resolution Quantum Lidar System (HRQLS) [2]. The average range accuracy of STIL is 11 cm, which is also close to that of the SPL data of HRQLS.

The area of interest, as shown in Fig. 1(a), includes ground, roofs, trees, and other objects. Three roof regions and two tree regions are detected and marked in Fig. 1(a). Close-up views of the regions Roof1, Roof2, Tree1, and Tree2 are displayed in Fig. 1(b), 1(c), 1(d), and 1(e), respectively.

 figure: Fig. 1.

Fig. 1. Point clouds of original data. (a) An overview of the original data. (b), (c), (d), and (e) The close-up views of regions Roof1, Roof2, Tree1, and Tree2, respectively. These sub-figures are colored by pseudo-elevation. (f) and (g) The simulated SPL points with 5 MHz noise generated from the original data of regions Roof1 and Tree1, with red points indicating signal and blue points indicating noise.

Download Full Size | PDF

2.2 Simulated SPL point clouds

Random Poisson-distribution noise with different point densities is added to the point cloud of STIL to simulate SPL point clouds under different atmospheric conditions, with values typical of nighttime conditions (0.5 MHz), clear daytime conditions (2 MHz), and hazy daytime conditions (5 MHz). As mentioned in the introduction, single-photon detector is susceptible to ambient light. The Noise rate means the respond number of the single-photon detector per unit time due to ambient light noise. Stronger ambient light results in a higher noise rate.

Figures 1(f) and 1(g) show the point clouds of Roof1 and Tree1 with 5 MHz noise. The red points are signal points, and the blue points are noise points. In Fig. 1(f), the signal and noise points fill the entire 135×134×30 m3 3D space. In Fig. 1(g), the signal points and noise points fill the whole 86×87×30 m3 3D space.

3. Methodology

3.1 The voxel-based spatial elongation filtering method

As described above, the VBSFM considers points in 27 voxels near the voxel of interest during threshold processing. Noise voxels (such as voxel (i) shown in Fig. 2(a)) near the signal voxel are likely to be considered as signal voxels, since the number of points within the 27 voxels may exceed the threshold. If we consider fewer voxels during threshold processing, the possibility of residual noise will be lower. However, it is also inadvisable to care about only one voxel during threshold processing, as in this case, the signal loss would be too large. Voxels containing few signal points and some noise points (such as voxel (iii) shown in Fig. 2(a)) can be distinguished as a noise voxel because the number of points in the voxel is less than the threshold.

 figure: Fig. 2.

Fig. 2. (a) Noise points and signal points in three voxels. Circular points are signal points, and cross points are noise points. (b) Original SPL points and additional points generated by our method. Solid points are original SPL points, and hollow points are additional ones.

Download Full Size | PDF

We propose a novel method that generates six additional points around the point of interest (above, below, left, right, in front and behind) before splitting up the 3D space, as shown in Fig. 2(b). Solid points represent the original SPL points, and hollow points are additional ones. Therefore, each voxel contains interior SPL points alongside the additional points generated by points inside and around the voxel. During the threshold processing, only the number of points inside the voxel of interest (including original SPL points and additional points) is compared to the threshold. In this situation, noise voxels (such as voxel (i) shown in Fig. 2(a)) will be removed, and signal voxels (such as voxel (iii) shown in Fig. 2(a)) will be retained. The additional points generated by the points in voxel (ii) may cause the number of points inside voxel (iii) to exceed the threshold. The whole flowchart of our method is shown in Fig. 3, and the detailed process is as follows.

  • (1) Generate six additional points next to each point in the SPL point cloud. The distances between the additional points and the one of interest are determined by the product of the elongation factor p and the size of the voxel. Assuming that the coordinates of the point of interest are (x0,y0,z0), the coordinates of the six additional points are (x0+p·a,y0,z0), (x0-p·a,y0,z0), (x0,y0+p·b,z0), (x0,y0-p·b,z0), (x0,y0,z0+p·c), and (x0,y0,z0-p·c). Here p is the elongation factor, and a, b, and c represent the length, width, and height of the voxel, respectively.
  • (2) Combine the SPL point clouds with all additional points, and generate composite point clouds.
  • (3) Split the 3D space into 3D voxels of a given size, and count the number of points in each voxel.
  • (4) Compare the count value of each voxel to a threshold. Voxels with a count value not less than the threshold are distinguished as signal voxels, and voxels with a count value less than the threshold are distinguished as noise voxels.
  • (5) Delete all additional points in the signal and noise voxels, the remaining points inside the signal voxels are signal points, and those inside the noise voxels are noise points.

 figure: Fig. 3.

Fig. 3. Flowchart of the voxel-based spatial elongation filtering method.

Download Full Size | PDF

3.2 Performance evaluation index of filtering effect —Fl

In previous works, the root-mean-square error (RMSE), false alarm rate, and signal loss rate were used to evaluate the filtering effect. As shown in Fig. 4, all three indexes vary monotonically as the threshold increases. Therefore, it is impossible to find the extremums of these three indexes and to regard the extremums as optimal values. Besides this, all these indexes evaluate the filtering effect one-sidedly. RMSE denotes the effects of the residual noise but largely depends upon the quality of the original signal. The false alarm rate and the signal loss rate represent the percentages of false alarm and signal loss, while the effects of the residual noise cannot be evaluated. To evaluate the filtering effects more comprehensively, we introduce Fl, as shown in Eq. (1).

$${F_l} = \frac{{k{T_\textrm{f}} + {F_\textrm{t}}}}{N}\Delta \bar{l}.$$
Here N represents the total number of original signal points, and $\Delta \bar{l}$ denotes the average distance between the residual noise points and their nearest signal points. Tf and Ft represent the signal loss number and the false alarm number, respectively, and k is the weight of Tf. Due to the extremely high data collection efficiency of the SPL, the signal loss number is of less importance than the false alarm number. In this paper, k is set to 0.5, which means that the importance of the false alarm number is twice that of the signal loss number. Fl denotes the influence of residual noise points and lost signal points after the filtering process. A smaller Fl indicates more effective filtering.

 figure: Fig. 4.

Fig. 4. Plots of Fl, RMSE, false alarm rate, signal loss rate, and $\Delta \bar{l}$ versus threshold for VBSFM and voxel-based spatial elongation filtering method. (a) Plots for VBSFM when the voxel size is 1×1×0.25 m3 and the noise density is 0.5 MHz. (b) Plots for our method when the voxel size is 1×1×0.5 m3 and the noise density is 0.5 MHz.

Download Full Size | PDF

As shown in Fig. 4, there is a minimum value of the Fl, caused by the combined influences of the false alarm and signal loss rates, we choose the minimum value of Fl as being the optimal filtering effect under this condition. The Fl bears a resemblance to the evaluation index used in J. Zhang’s work [5]. Both indexes use the signal loss number and the false alarm number to evaluate the filtering effect. However, there are cases in which the signal loss and false alarm number are relatively small, and the residual noise points are far away from the signal points, which leads to a large adverse impact on the overall accuracy. To evaluate the effects of residual noise points, we calculate $\Delta \bar{l}$ and introduce it into the expression of Fl.

4. Results and discussion

4.1 Optimal elongation factor of VBSEFM

As shown in Fig. 3, there is a coefficient elongation factor in the VBSEFM, which greatly affects the filtering. To set the coefficient correctly, we calculate Fl under various voxel sizes, elongation factors, noise densities, and regions. Besides this, the minimum Fl under each condition is chosen to represent the optimal filtering effect. Regions Roof1 and Tree1 are selected to represent the roof and tree regions, respectively. We optimize the elongation factor based on the principle that Fl, under different noise densities, should be kept at a relatively low value. The determined optimal elongation factors for different voxel sizes are given in Table 1.

Tables Icon

Table 1. The optimal elongation factors of VBSEFM for different voxel sizes

4.2 Filtering effect comparison

The minimum Fl, obtained by the VBSFM and our method, for the Roof1 and Tree1 regions under different noise densities and voxel sizes are presented in Fig. 5. As can be seen from the figure, when the voxels sizes are large (e.g., 4×4×1m3, 3×3×1 m3, and 2×2×1 m3), the Fl of our method are less than those of the VBSFM, whereas when the voxel sizes are small (e.g., 1×1×0.25 m3 and 1×1×0.125 m3), the Fl of our method are greater than those of VBSFM.

 figure: Fig. 5.

Fig. 5. The Fl obtained by VBSFM and our method under different noise densities, voxel sizes, and regions: (a) Region Roof1; (b) Region Tree1.

Download Full Size | PDF

This is a result of the difference in the number of voxels considered during threshold processing. There is an optimal voxel size for each voxel-based method. When the voxel size is small enough, the number of points in each voxel can be 1 or 0. In this situation, it is impossible to distinguish between signal voxels and noise voxels. On the other hand, a larger voxel size results in poor spatial resolution, which is detrimental to the filtering effect.

Since our method considers only one voxel during threshold processing, the spatial resolution of our method is much better than that of VBSFM when the voxel size is relatively large. However, when the voxel size is relatively small, it is more likely for the VBSFM to distinguish between signal voxels and noise voxels.

To evaluate both methods, we compare the filtering effects of the two methods in their optimal voxel sizes. We choose the optimal voxel size based on the principle that the Fl in the optimal voxel size should be kept at a relatively low value. Finally, 1×1×0.25 m3 and 1×1×0.5 m3 are chosen as the optimal voxel sizes for the VBSFM and our method, respectively.

We calculate the Fl of both methods in their optimal voxel sizes, and the results are shown in Table 2. It can be seen that the Fl of both methods is close in roof regions, and our method has a smaller Fl in tree regions. When the noise density is greater, the Fl is also larger, and it has a greater chance of finding the difference in the filtering effect between the two filtering methods. Therefore, we calculate performance evaluation indexes in each region with 5 MHz noise to further compare the filtering effects of both methods. The results are displayed in Table 3.

Tables Icon

Table 2. The Fl of two filtering methods in different noise densities and regions

Tables Icon

Table 3. Performance evaluation indexes of two filtering methods in Roof1 and Tree1 regions, when noise density is 5 MHz. The false alarm rate is calculated by dividing the number of false alarm points with the number of original signal points. The signal loss rate is calculated by dividing the number of signal loss points with the number of original signal points.

Table 3 shows that the false alarm rates using our method in all regions (3.7%, 3.5%, 3.6%, 3.2%, 3.2%, 3.8%, 3.4%, and 3.3%) are less than those using the VBSFM (4.4%, 4.1%, 4.0%, 3.8%, 4.3%, 4.5%, 4.6%, and 4.3%). These results indicate that our method can remove noise more effectively. Besides this, the $\Delta \bar{l}$ of our method are smaller than those of VBSFM. This is because many of the residual noises of our method occur in voxel (iii), while many of the residual noises for VBSFM are in voxel (i), as shown in Fig. 2(a). The $\Delta \bar{l}$ of noise points in voxel (iii) is smaller than the $\Delta \bar{l}$ of noise points in voxel (i). Additionally, the signal loss rates of the tree regions are larger than those of roof regions. This is because of the complex structure of the point clouds in the tree region, and the low density of the point clouds.

Figure 6 shows the processed point clouds for the regions Roof1 and Tree2, after applying the VBSFM and our method. It can be seen that the processed point clouds of the VBSFM still feature many noise points near the signal. In contrast, there are fewer residual noise points in the processed point clouds of our method, and almost all features of the original signal are preserved, such as the door eaves and tree crowns in Fig. 6.

 figure: Fig. 6.

Fig. 6. Processed point clouds under different conditions. (a) Processed point clouds of region Roof1 with 5 MHz noise after using VBSFM, where the voxel size of VBSFM is 1×1×0.25 m3; (b) processed point clouds of Roof1 with 5 MHz noise after using our method, where the voxel size of our method is 1×1×0.5 m3; (c) processed point clouds of region Tree2 with 5 MHz noise after using VBSFM, where the voxel size of VBSFM is 1×1×0.25 m3; (d) processed point clouds of Tree2 with 5 MHz noise after using our method, where the voxel size of our method is 1×1×0.5 m3. Red points are signal points, and blue points are noise points.

Download Full Size | PDF

5. Conclusion

In this paper, a voxel-based spatial elongation filtering method is proposed, and the optimal elongation factors of this method under different voxel sizes are determined. Before splitting the 3D space into voxels, six points surrounding (above, below, left, right, in front, and behind) each point in the SPL data are generated. These additional points enable our method to make better use of the coordinates of each point in the SPL data.

Besides this, a comprehensive performance index Fl (including the false alarm rate, the signal loss rate and the average distance between the residual noise point and their nearest signal point), is introduced to evaluate the filtering effect. The Fl of VBSFM and our method under various voxel sizes, noise densities, and regions are calculated. Comparisons of the filtering effects of both methods in their optimal voxel sizes (1×1×0.25 m3 for the VBSFM and 1×1×0.5 m3 for our method) are carried out. The results show that our method has smaller Fl than those of VBSFM in tree regions, and similar Fl to VBSFM in roof regions. Moreover, the false alarm rates of our method in all these regions (3.7%, 3.5%, 3.6%, 3.2%, 3.2%, 3.8%, 3.4%, and 3.3%) are less than those of the VBSFM (4.4%, 4.1%, 4.0%, 3.8%, 4.3%, 4.5%, 4.6%, and 4.3%). In other words, the average false alarm rate using our method (3.5%) is 18.6% smaller than that using VBSFM (4.3%). The results show that our method can remove noise more effectively than the VBSFM whilst still retaining the signal. This method can be used to remove noise in SPL data and therefore improve the data quality, especially in areas containing vegetation.

Funding

National Key Scientific Instrument and Equipment Development Projects of China (2012YQ040164).

Disclosures

The authors declare no conflicts of interest.

References

1. J. J. Degnan, “The impact of single photon SLR technology on large scale topo-bathymetric mapping,” International Workshop on Laser Ranging, (2016).

2. J. J. Degnan, “Scanning, multibeam, single photon lidars for rapid, large scale, high resolution, topographic and bathymetric mapping,” Remote Sens. 8(11), 958 (2016). [CrossRef]  

3. X. Wang, Z. Pan, and C. Glennie, “A novel noise filtering model for photon-counting laser altimeter data,” IEEE Geosci. Remote. Sens. Lett. 13(7), 947–951 (2016). [CrossRef]  

4. M. Awadallah, S. Ghannam, L. Abbott, and A. Ghanem, “Active contour models for extracting ground and forest canopy curves from discrete laser altimeter data,” 13th International Conference on Lidar Applications for Assessing Forest Ecosystems, (2013).

5. J. Zhang and J. Kerekes, “An adaptive density-based model for extracting surface returns from photon-counting laser altimeter data,” IEEE Geosci. Remote. Sens. Lett. 12(4), 726–730 (2015). [CrossRef]  

6. M. S. Moussavi, W. Abdalati, T. Scambos, and A. Neuenschwander, “Applicability of an automatic surface detection approach to micro-pulse photon-counting lidar altimetry data: implications for canopy height retrieval from future ICESat-2 data,” Int. J. Remote Sens. 35(13), 5263–5279 (2014). [CrossRef]  

7. M. S. Awadallah, A. L. Abbott, V. A. Thomas, R. H. Wynne, and R. F. Nelson, “Estimating forest canopy height using photon-counting laser altimetry,” 13th International Conference on Lidar Applications for Assessing Forest Ecosystems, (2013).

8. L. A. Magruder, M. E. Wharton, K. D. Stout, and A. L. Neuenschwander, “Noise filtering techniques for photon-counting ladar data,” Proc. SPIE, (2012).

9. U. C. Herzfeld, B. W. Mcdonald, B. F. Wallin, T. A. Neumann, T. Markus, A. Brenner, and C. Field, “Algorithm for detection of ground and canopy cover in micropulse photon-counting lidar altimeter data in preparation for the ICESat-2 mission,” IEEE Trans. Geosci. Remote Sens. 52(4), 2109–2125 (2014). [CrossRef]  

10. K. H. Horan and J. P. Kerekes, “An automated statistical analysis approach to noise reduction for photon-counting lidar systems,” 2013 IEEE Geosci. Remote Sens. Symp. (2013).

11. D. Gwenzi and M. A. Lefsky, “Prospects of photon counting lidar for savanna ecosystem structural studies,” Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. XL-1(1), 141–147 (2014). [CrossRef]  

12. A. Swatantran, H. Tang, T. Barrett, P. Decola, and R. Dubayah, “Rapid, high-resolution forest structure and terrain mapping over large areas using single photon lidar,” Sci. Rep. 6(1), 28277 (2016). [CrossRef]  

13. H. Tang, A. Swatantran, T. Barrett, P. Decola, and R. Dubayah, “Voxel-based spatial filtering method for canopy height retrieval from airborne single-photon lidar,” Remote Sens. 8(9), 771 (2016). [CrossRef]  

14. J. M. Stoker, Q. A. Abdullah, A. Nayegandhi, and J. Winehouse, “Evaluation of single photon and Geiger mode lidar for the 3D elevation program,” Remote Sens. 8(9), 767 (2016). [CrossRef]  

15. X. Wang, C. Glennie, and Z. Pan, “An adaptive ellipsoid searching filter for airborne single-photon lidar,” IEEE Geosc. Remote Sens. Lett. 14(8), 1258–1262 (2017). [CrossRef]  

16. X. Wang, C. Glennie, and Z. Pan, “Adaptive noise filtering for single photon lidar observations,” IEEE Geosc. Remote Sens. Symp.Fort Worth, USA, (2017).

17. X. Wang, Z. Pan, and C. Glennie, “Weak echo detection from single photon lidar data using a rigorous adaptive ellipsoid searching algorithm,” Remote Sens. 10(7), 1035 (2018). [CrossRef]  

18. G. Ye, R. Fan, Z. Chen, W. Yuan, D. Chen, and P. He, “Range accuracy analysis of streak tube imaging lidar systems,” Opt. Commun. 360, 7–14 (2016). [CrossRef]  

19. Z. Chen, R. Fan, G. Ye, T. Luo, J. Guan, Z. Zhou, and D. Chen, “Depth resolution improvement of streak tube imaging lidar system using three laser beams,” Chin. Opt. Lett. 16(4), 041101 (2018). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (6)

Fig. 1.
Fig. 1. Point clouds of original data. (a) An overview of the original data. (b), (c), (d), and (e) The close-up views of regions Roof1, Roof2, Tree1, and Tree2, respectively. These sub-figures are colored by pseudo-elevation. (f) and (g) The simulated SPL points with 5 MHz noise generated from the original data of regions Roof1 and Tree1, with red points indicating signal and blue points indicating noise.
Fig. 2.
Fig. 2. (a) Noise points and signal points in three voxels. Circular points are signal points, and cross points are noise points. (b) Original SPL points and additional points generated by our method. Solid points are original SPL points, and hollow points are additional ones.
Fig. 3.
Fig. 3. Flowchart of the voxel-based spatial elongation filtering method.
Fig. 4.
Fig. 4. Plots of Fl, RMSE, false alarm rate, signal loss rate, and $\Delta \bar{l}$ versus threshold for VBSFM and voxel-based spatial elongation filtering method. (a) Plots for VBSFM when the voxel size is 1×1×0.25 m3 and the noise density is 0.5 MHz. (b) Plots for our method when the voxel size is 1×1×0.5 m3 and the noise density is 0.5 MHz.
Fig. 5.
Fig. 5. The Fl obtained by VBSFM and our method under different noise densities, voxel sizes, and regions: (a) Region Roof1; (b) Region Tree1.
Fig. 6.
Fig. 6. Processed point clouds under different conditions. (a) Processed point clouds of region Roof1 with 5 MHz noise after using VBSFM, where the voxel size of VBSFM is 1×1×0.25 m3; (b) processed point clouds of Roof1 with 5 MHz noise after using our method, where the voxel size of our method is 1×1×0.5 m3; (c) processed point clouds of region Tree2 with 5 MHz noise after using VBSFM, where the voxel size of VBSFM is 1×1×0.25 m3; (d) processed point clouds of Tree2 with 5 MHz noise after using our method, where the voxel size of our method is 1×1×0.5 m3. Red points are signal points, and blue points are noise points.

Tables (3)

Tables Icon

Table 1. The optimal elongation factors of VBSEFM for different voxel sizes

Tables Icon

Table 2. The Fl of two filtering methods in different noise densities and regions

Tables Icon

Table 3. Performance evaluation indexes of two filtering methods in Roof1 and Tree1 regions, when noise density is 5 MHz. The false alarm rate is calculated by dividing the number of false alarm points with the number of original signal points. The signal loss rate is calculated by dividing the number of signal loss points with the number of original signal points.

Equations (1)

Equations on this page are rendered with MathJax. Learn more.

F l = k T f + F t N Δ l ¯ .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.