EVALUATING VOXEL ENABLED SCALABLE INTERSECTION OF LARGE POINT CLOUDS

Laser scanning has become a well established surveying solution for obtaining 3D geo-spatial information on objects and environment. Nowadays scanners acquire up to millions of points per second which makes point cloud huge. Laser scanning is widely applied from airborne, carborne and stable platforms, resulting in point clouds obtained at different attitudes and with different extents. Working with such different large point clouds makes the determination of their overlapping area necessary but often time consuming. In this paper, a scalable point cloud intersection determination method is presented based on voxels. The method takes two overlapping point clouds as input. It consecutively resamples the input point clouds according to a preset voxel cell size. For all non-empty cells the center of gravity of the points in contains is computed. Consecutively for those centers it is checked if they are in a voxel cell of the other point cloud. The same process is repeated after interchanging the role of the two point clouds. The quality of the results is evaluated by the distance to the pints from the other data set. Also computation time and quality of the results are compared for different voxel cell sizes. The results are demonstrated on determining he intersection between an airborne and carborne laser point clouds and show that the proposed method takes 0.10%, 0.15%, 1.26% and 14.35% of computation time compared the the classic method when using cell sizes of of 10, 8, 5 and 3 meters respectively.


INTRODUCTION
Light Detection And Ranging (LiDAR) mapping techniques enable to quickly acquire geometric information in 3D (Vosselman and Maas, 2010).Compared to airborne photogrammetry, laser scanning has the advantages of high speed and high quality acquisition of points (Baltsavias, 1999a, Rönnholm et al., 2007).In recent years there is a rapid development of different platforms for laser scanning, such as Airborne Laser Scanning (ALS), Terrestrial Laser Scanning (TLS) and Mobile Laser Scanning (MLS).Laser scanning has already been applied to situations, such as tunnel inspection and monitoring (Gosliga et al., 2006) and road engineering (Wang et al., 2013, Jaakkola et al., 2008, Kumar et al., 2013, Pu et al., 2011), vegetation and land cover classification (Bater et al., 2011, Antonarakis et al., 2008, Yunfei et al., 2008, Puttonen et al., 2010), and object extraction and reconstruction (Olsen et al., 2010, Brasington et al., 2012).However, due to the scanning mechanism of the different systems, the obtained point clouds have different properties (Vosselman and Maas, 2010), such as non-uniform point density and occlusion effects (Wehr andLohr, 1999, Forest andSalvi, 2002).To fully sample a complicated scan, it may be required to scan its objects in multiple runs, and from multiple platforms.
Since different platforms have their own advantageous perspective, multiple platforms are now combined to obtain complete coverage of objects (Zhou and Vosselman, 2012, Baltsavias, 1999b, Holopainen et al., 2013, Sithole and Vosselman, 2004).Once point clouds from multiple sources are used, the determination of area of intersection becomes a compulsory operation.2D and 3D intersection have been studied and discussed in the fields of com- * Corresponding author putational geometry, computer graphics and robotics.In (Shamos and Hoey, 1976), optimal algorithms were developed to determine the intersection of geometric objects in 2D.In (Dobkin and Kirkpatrick, 1983) low complexity methods for suitably preprocessed polyhedral intersections were distinguished based on convex.In (Mount, 1992), by constructing an enveloping triangulation called a scaffold two preprocessed simple polygons are detected whether they intersect one another.And in (Wang, 1996) the intersection curves of two parametric surfaces are presented based on the concept of normal projection.In (Manocha andCanny, 1991, Grandine andKlein, 1997) new approaches were designed by determining numerical solution of differential algebraic equations.In 3D, approaches were also developed for intersection detection.In (Chazelle, 1989, Dobkin andKirkpatrick, 1985), an intersection algorithm based on the checking of the relationship between vertices and planes was presented.In (Pan et al., 2012), intersection of objects was determined by detecting objects' bounding boxes, and then objects collision checking and proximity computations were performed.Those literature did give the concept of object intersection determination, but the objects all had a primitive geometry shape and the number of objects is limited.The latest scanners obtain up to millions of high precision points (Vosselman andMaas, 2010, Hyyppä, 2011), and the data sets generated by those laser scanning systems are huge.
To determine intersections either it would be needed to convert the discrete point clouds to standard geometry objects or using classic point based methods, which make the procedure computational expensive.Varela-González et al., 2013).Some open source softwares such as CloudCompare (CloudCompare, 2015) and LasTools (Isenburg, 2015), have segmentation tools accessible from a Graphical User Interface (GUI) but don't have data based intersection tools.Also the classic intersection determination based on point clouds neighborhood searching is computational expensive and unable to deal with huge point clouds.
Therefore in this paper a novel point cloud intersection method is proposed that systematically visits only part of the points in the input point clouds.
This study presents a method for automatic and scalable area of intersection determination from two huge point cloud data.The results are evaluated by computing the minimum distance between the two overlapping point clouds at different scales.Also, the quality and efficiency of this method are compared to the classic point based intersection method.

METHODOLOGY
Intersection problems are fundamental to many aspects of geometric computing.In 2D, the intersection of objects is defined as the set of points or area shared by both objects (Shamos andHoey, 1976, Grandine andKlein, 1997).While in 3D, the intersection of two objects is the volume that shared by them (Dobkin andKirkpatrick, 1985, Dobkin andKirkpatrick, 1990).However, those definitions are from objects that have a definite boundary.
For raw point clouds acquired by laser scanning, there are no clear boundaries.So the proposed method is an approximation method.As shown in Figure 1, the overall methodology presented in this paper consists of four steps.

Point clouds resampling
The imported point clouds are re-sampled as voxel cells.As a P ixel is a 2D representation, a voxel is a cube with a predefined edge length in 3D.This edge length defines the resolution of the voxel resampling procedure.the points that are within the geometric space of this voxel cell; and (iii), coordinates of the 3D center of gravity of all the points stored in the point indices container.

Intersected area determination
We call one of the imported point clouds Source and the other T arget.After the data sets are resampled as voxel cells and the center of gravity of the points in all cells is computed, the area of intersection is determined at the voxel cell level.

Result evaluation
The evaluation of the results consists of three parts.First the method presented in this paper is implemented and applied at different scales, which means the point cloud data are resampled by voxel cells of different cell sizes.Secondly, point based intersection determination is carried out with different radii on the same test point clouds.Finally, the results of the two different methods are evaluated and compared with respect to their distances.

Test data sets
The test data consists of two point clouds from different platforms.One is from airborne laser scanning (ALS) and the other from mobile laser scanning (MLS).General information of the two point clouds is given in   (Fugro, 2015).the overlap of the two pint clouds is considerable and contains several roads and buildings of the TUD campus.

Scalable testing of the method
Corresponding to the intersection results for four different voxel cell sizes, cloud to cloud distances are computed between the Source and T arget point clouds.The distances between two point clouds represents how far of each point in the Source, the closest point in the T arget point cloud is.We employ this as one aspect of intersection quality evaluation.The bigger the distance, the coarser the intersection area determination, while the smaller, the better the results.As can be seen from Figure 7, the smaller the voxel cell size used for computing, the smaller the distances between the two intersected point clouds.This denotes that if one wants better results, a smaller voxel cell size shall be used.Beaulieu and Rajwani, 2004).The probability density function is given in Equation 1.Here µ is the mean and σ is the standard deviation of ln(x), also named location parameter and scale parameter respectively.
The statistical analysis of the distances computed from the four intersection determination results are shown in Figure 7.It can be noted from the figure that the quality of the point cloud intersection varies upon the voxel cell size used.The parameters of the four results are shown in Table 3.2.While the voxel cell size changes from big to small, the average distances of the results are getting small as well.As learnt from the precessing with those different voxel cell sizes, the mean of the distances are around 20% of the voxel cell size.Also the mean and standard deviation of the Log-Normal distribution follow the same trend.

Comparison with point based intersection results
The point-wise intersection method first needs a predefined cloud to cloud distance as a criterion that defines the area of intersection.The intersection is defined as the points that are closer than the 1.0 meter distance threshold from each other.To evaluate the influence of point numbers to the processing, the original point clouds are down-sampled to different scales.For each scale, the computation time of the two methods are compared.In this study, the processing time of different voxel cell sizes are also compared.
Both algorithms are implemented in C and run on a Dell desktop computer, with a Intel E5-1620 processor and a memory of 16GB.
During the tests, the original point clouds are down-sampled at different percentages.As described in Table 3.3, the original point clouds are uniformly down-sampled to three scales.
Scenario MLS pts number (%) ALS pts number(%)  The computational time of voxel based method is related to the complexity of the algorithm.Suppose the bounding box of the point cloud is D and the voxel cell size is d, then the number of voxel cell is N = ( D d ) 3 .For each cell, its maximum and minimum boundary points are checked whether this cell contains the center point or not.So their are at least 32 operations for each cell.Then the total cost is at least Nt = 2 × ( D d ) 3 .For an imported point cloud, its spatial extent is constant and then the computational time will only related to the voxel cell sizes.Letting di and dj are two different voxel cell sizes that applied to the same point cloud, the ratio of their total computational costs is Rij = N i N j = ( , where di and dj are the voxel cell sizes.The smaller the voxel cell size, the larger the number of occupied cells and the longer the computational time.However, as the smaller the voxel cell size, there will also more empty cells.This letting the increase ratio of the non-empty cells be smaller than ( So the upper bound of computational expenses related to different cell size is ( . The computational time of the testing data sets is illustrated in Figure 8.For the presented voxel based methods, as denoted by the bottom horizontal label and the blue line, the computation time of voxel based method on different sizes.For the point based method, the more of the points number, the longer of the computation time it takes, as illustrated by the red line in Figure 8.In Figure 9 the computation time of the two methods is illustrated.The red line is the computation time of the point based intersection determination.For each intersection the distances between overlapping points in the intersection were calculated.In Figure 10 the mean and standard deviation of the distances are shown.It can be noticed that there is only small variation in the mean and standard deviation of the four pairs of testing scenarios when using the same distance threshold.From the first pair to the fourth pair, the computation time varies from 8.2 seconds to 1255.8 seconds, as shown in Figure 8.
For voxel based method, focusing on the original point clouds, the computation time of the presented method is 0.10%, 0.15%, 1.26% and 14.35% of the time needed by point based methods for voxel cell sizes of 10, 8, 5 and 3 meters respectively.
For the voxel based methods, Figure 9 also shows that computation time is robust to the number of points of the two clouds.However, the voxel cell size, which is also interpreted as scale of detail, effects the computation time and quality of results.If the voxel cell size is smaller, which also means the results get more accurate, the time expense will be more.However, the computational expense will also get more an more.Afterwards the method is compared with a traditional point based method.The processing times are compared for different scenarios of down-sampled point cloud pairs.

ISPRS
Experimental results show that to meet the same quality of results, voxel based methods take much less time than point based method.And smaller voxel cell sizes will generate more accurate results, however, computational expenses will also get higher.
Experimental results show that the quality of tradition point based method is related to the distance threshold and is generally better than that of the proposed method.However, the point based method is highly dependent on the number of points.While the voxel based method has more flexible choices on quality and computation time.Results on quality and computation time are balanced by the voxel cell size.
During processing and studying the presented methods, we also learned some experiences.Firstly, to more efficiently determine high quality area of intersection using this method, a relatively larger voxel cell size can be applied in the first run to get coarse results.Then based on these results, finer results can be obtained with a smaller voxel cell size.Secondly, the voxel cell resampling is based on uniform space sub-division and locations having no points are also sampled as empty cells.This introduces redundant memory space.As a solution, an octree data structure can be implemented to improve memory efficiency.

Figure 1 :
Figure 1: Overall methodology of the algorithm (1) Two overlapping point clouds are imported as input.In the figure, the red and green points are the two point clouds respectively.(2) The point clouds are re-sampled as voxel cells and the center of gravity of the points in each cell is computed.(3) The geometric relation between two voxels is checked and the area of intersection is determined.(4) The quality of the results is evaluated by computing the point based distances between the two point cloud data sets restricted to their determined intersection.
Figure 2 is an illustration of a voxel cell in 3D, where Pmax is the upper right back point and Pmin is the lower left front point of the cell.
Figure 4 illustrates this strategy.The green rectangles represent the voxel cells, which are generated from Source where the blue dots are the center of gravity of the voxel cells from T arget.Once the center point of a cell from T arget is detected within a voxel cell from Source, this cell is then regarded as a cell that belongs to the intersection.The green shaded cells together denote the area of intersection.

Figure 4 :
Figure 4: Determination of the area of intersection of two overlapping point clouds

Figure 6 :
Figure 6: Two overlapping multiple platform point clouds.
The next step is to quantify the distance computation results.To do this, statistical analysis is conducted on the results of the four different voxel cell sizes.Suppose the distances are random.Since the values of distances are all non-negative, this variable is a typical Log-Normal distribution (Beaulieu and Xie, 2004, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-3/W5, 2015 ISPRS Geospatial Week 2015, 28 Sep -03 Oct 2015, La Grande Motte, France This contribution has been peer-reviewed.The double-blind peer-review was conducted on the basis of the full paper.Editors: S. Oude Elberink, A. Velizhev, R. Lindenbergh, S. Kaasalainen, and F. Pirotti doi:10.5194/isprsannals-II-3-W5-25-2015

Figure 7 :
Figure 7: Cloud to cloud distances and statistical analysis from four different voxel cell sizes.

Figure 8 :
Figure 8: Computation time for different voxel cell sizes and down-sampling scenarios.
Figure 9: Computation time of the two methods

Table 1 :
Table 3.1.Both of the two point clouds sample part of the campus of Delft University of Technology (TUD).Information on multiple platform point clouds

Table 2 :
Mean and Std.dev. of the distances

Table 3 :
Scenarios of down-sampled pointsFor the Source point cloud 7.31%, 22.9% and 58.67% of the original points are kept and for the T arget point cloud 6.86%, 25.21% and 69.77% of the points are kept for comparison.Based on those scenarios, the two algorithms are compared.