DIMENSIONALITY BASED SCALE SELECTION IN 3D LIDAR POINT CLOUDS

Abstract. This papers presents a multi-scale method that computes robust geometric features on lidar point clouds in order to retrieve the optimal neighborhood size for each point. Three dimensionality features are calculated on spherical neighborhoods at various radius sizes. Based on combinations of the eigenvalues of the local structure tensor, they describe the shape of the neighborhood, indicating whether the local geometry is more linear (1D), planar (2D) or volumetric (3D). Two radius-selection criteria have been tested and compared for ﬁnding automatically the optimal neighborhood radius for each point. Besides, such procedure allows a dimensionality labelling, giving signiﬁcant hints for classiﬁcation and segmentation purposes. The method is successfully applied to 3D point clouds from airborne, terrestrial, and mobile mapping systems since no a priori knowledge on the distribution of the 3D points is required. Extracted dimensionality features and labellings are then favorably compared to those computed from constant size neighborhoods.


INTRODUCTION
Point cloud data from airborne and terrestrial devices provide a direct geometrical description of the 3D space. Such information is reliable, of high accuracy but irregular and not dense. However, the underlying structures and objects may be detected among sets of close 3D points. The local geometry is estimated by the distribution of points in the neighborhood. Finding the best neighborhood for each point is a main issue for a large variety of common processes: data downsampling, template fitting, feature detection and computation, interpolation, registration, segmentation, or modelling purposes. The notion of neighborhood and its fundamental properties are fully described in (Filin and Pfeifer, 2005). The neighbors of a lidar point are traditionally retrieved by finding the k nearest neighbors or all the points included in a small restricted environment (sphere or cylinder) centered on the point of interest. The main problem stems from the fact that the k and environment radius values are (1) usually heuristically chosen, and (2) assumed to be constant for the whole point cloud, instead of being guided by the data. This does not ensure that all theses neighbors belong to the same object as the current point. Therefore, its local description may be biased when including several distinct structures, and provides erroneous feature descriptors. Moreover, the relative variation in the spatial extent of geometrical structures is ignored. For aerial datasets, problems will occur at the borders between objects and for objects which size is inferior or close to the neighborhood size. For terrestrial datasets, in addition, the point density may significantly fluctuate due to foreground object occlusion, dependence on distance and relative orientation of the objects (Soudarissanane et al., 2009), leading to data sparseness or irregular sampling. This paper aims at proposing a methodology to find the optimal neighborhood radius for each 3D point on a lidar point cloud. An "optimal" neighborhood is defined as the largest set of spatially close points that belong to the same object as the point of interest. The inclusion of points lying on different surfaces is prohibited. The context of the study is rather general: in order to be applicable both on terrestrial (TLS) and airborne (ALS) datasets, the method is simply based on the point location, without requiring knowledge on intensity, echo number or full waveforms. The topology resulting from the sequential acquisition of the data is also considered to be lost ("unorganized" point cloud), preventing the adoption of specific scan line grouping methods (Hadjiliadis and Stamos, 2010). Furthermore, the process is designed out of the scope of any application, even if the final goal is indeed to be beneficial to any application requiring a correct local description around each 3D point. The problem of scale selection has been mainly tackled for surface reconstruction and feature extraction of scanned opaque objects. Several approaches have therefore been developed for noisy point clouds and irregular sampling issues. Most of them are surface-based i.e., they try to fit a curve or a surface of some form to the 3D point cloud (Pauly et al., 2006). Finding the optimal group that well represents the local geometrical properties is performed using indicators such as the normal and/or the curvature (Hoppe et al., 1992;Zwickler et al., 2002;Dey and Sun, 2005;Belton and Lichti, 2006). For instance, it is retrieved by minimizing the upper bound on angular error between the true normal and the estimated one. Starting for the minimal possible subset around the point of interest, the neighborhood is iteratively increased until the angular variance reaches a predefined threshold (Mitra et al., 2004). Such works have been theoretically improved by Lalonde et al. (2005), no more requiring knowledge on the data distribution, and applied to mobile mapping datasets. An alternative work, based of the expression of the positional uncertainty, is introduced in (Bae et al., 2009) for processing TLS datasets. However, these methods are effective for smoothly varying surfaces and may not be adapted to anthropic surfaces acquired with various kinds of lidar systems. Furthermore, the model-based assumption does not hold when dealing with objects without predefined shapes (e.g., vegetated areas) or with noise stemming for relief high frequencies (e.g., chimneys and facades for ALS data or pedestrians and points inside buildings for TLS data). Consequently, in our context, a more suitable solution is to directly compute shape features (Gumhold et al., 2001;Belton and Lichti, 2006), i.e., low-level primitives that may capture the variability of natural environments. The proposed methodology is developed in Section 2. The shape features are described in Section 3. The computation of the optimal neighborhood radius embedded in a multi-scale framework is proposed in Section 4. Results on various laser scanning datasets are presented in Section 5, and conclusions are drawn in Section 6.

Description
The methodology aims at finding the optimal neighborhood radius for each lidar point, working directly and exclusively in the 3D domain, without relying on surface descriptors (such as normals) or structures (such as triangulations or polygonal meshes). It is composed on two main steps: 1. Computation of three dimensionality features for each point, between predefined minimal and maximal neighborhood scale. These features describe the distribution of the points in the 3D space, and more exactly, the matching between the local point cloud and each of the three dimensionalities (linear, planar or volumetric).
2. Scale selection: retrieval of the neighborhood radius for which one dimensionality is most dominant over the two others.
The three dimensionality features (a1D-a2D-a3D) are computed exhaustively, at each point and for each acceptable neighborhood scale, from the local covariance matrix. An isotropic spherical neighborhood, centered on the point of interest, is adopted for this purpose. These low-level features, as well as the automatic set up of the radius lower and upper bounds are described in Section 3. Then, the optimal radius is retrieved by comparing the behaviours of these three features between the minimal and maximal acceptable radius. Two radius-selection criteria are tested to evaluate each scale and find the most relevant value. This multi-scale analysis is performed in order to capture variation in shape when aggregating points for an object distinct from the object of interest (edge effect or outliers). It is also useful in case of significant density variation and lack of support data for gathering points over a large volume while ensuring the conservation of the favorite dimensionality. Finally, our method presents three interesting characteristics: • Definition of a confidence index of the saliency of one dimensionality over the two other ones.
• Multi-scale analysis and automatic set up of the bounding scales. • Labelling of each point according to its privileged dimensionality, providing an interesting basis for segmentation and classification algorithms.

Datasets
In order to assess the relevance of the proposed approach for various point densities, point distributions and points of view, three kinds of lidar datasets are tested: airborne, terrestrial static, and acquired with a mobile mapping system (named ALS, TLS, and MMS, respectively).
ALS : Three datasets are used. The first one (ALSG) has been acquired over Biberach (Germany), covering both residential and industrial areas as well as a city center with small buildings (point density of 5 pts/m 2 ). The second dataset (ALSR) concerns a residential area in Russia, with 5 pts/m 2 (Shapovalov et al., 2010). Finally, the third one covers the dense city center of Marseille (France), with high buildings, and thus sparse points on the building facades (ALSF ). Three parallel strips are present: the point density therefore varies between 2 and 4 pts/m 2 (for one strip and for the overlapping areas, respectively).
TLS : The terrestrial scans acquired over the Agia Sanmarina church (Greece) have been processed (Bae et al., 2009). The dataset first offers a large variety of structures of various sizes as well as sparse vegetation on the ground. Furthermore, the point density significantly varies with the orientation of the surfaces with respect to the scanner position.
MMS : Datasets over two urban areas (France and United States, respectively MMSF and MMSU ) from distinct mobile mapping systems have been processed (Munoz et al., 2009). Such datasets also include man-made objects of various sizes and shapes, with varying point densities. The two specificities of MMS datasets are (1) gaps in the point cloud due to the occlusion of foreground objects, and (2) vertical privileged directions in the point clouds due to the sequential acquisition by lines.

SHAPE FEATURES
In this work, we focus on shape features computed on the neighborhood of a lidar point. The neighborhood V r P of a point P at scale r is defined as the set of points P k verifying : Thus, the neighborhood environments are 3D spherical volumes. They ensure isotropy and rotation invariance, such that the computed shape descriptors are not biased by the shape of the neighborhood. Furthermore, r is the single parameter to be optimized. A classical approach consists in performing a Principal Component Analysis (PCA) of the 3D coordinates of V r P (Belton and Lichti, 2006;Gross et al., 2007). This statistical analysis uses the first and second moments of V r P , and results in three orthogonal vectors centered on the centroid of the neighborhood. The PCA synthesizes the distribution of points along the three dimensions (Tang et al., 2004), and thus models the principal directions and magnitudes of variation of the point distribution around the center of gravity. These magnitudes are combined to provide shape descriptors for each of the three dimensions. More advanced features based on harmonics or spin images (Frome et al., 2004;Golovinskiy et al., 2009) are not necessary since the segmentation task, which indeed requires contextual knowledge, is not tackled in this paper.

Principal Component Analysis
,n xi the center of gravity of the n lidar points of V r P .
Since C is a symmetric positive definite matrix, an eigenvalue decomposition exists and can be expressed as C = RΛR T , where R is a rotation matrix, and Λ a diagonal, positive definite matrix, known as eigenvector and eigenvalue matrices, respectively. The eigenvalues are positive and ordered so that λ1 ≥ λ2 ≥ λ3 > 0. ∀j ∈ [1, 3], σj = λj, denotes the standard deviation along the corresponding eigenvector − → vj . Thus, the PCA allows to retrieve the three principal directions of V r P , and the eigenvalues provide their magnitude. The average distance, all around the center of gravity, can also be modeled by a surface. The shape of V r P is then represented by an oriented ellipsoid. The orientation and the size informations are divided between R and Λ : R turns the canonical basis into the orthonormal basis ( − → v1, − → v2, − → v3) and √ Λ transforms the unit sphere to an ellipsoid (σ1, σ2 and σ3 being the lengths of the semi-axes). As enhanced in Figure 1 and for instance in (Gross et al., 2007), such an ellipsoid reveals the linear, planar or volumetric behaviour of the neighborhood i.e., whether the point set is spread in one, two or three dimensions (blue, gray, and green ellipsoids in Figure 1, respectively). Figure 1: Three examples of ellipsoids computed over three areas of interest of distinct dimensionalities for the TLS dataset (a tripod over a low and sparse vegetation -height colored).
We chose µ = σ1 because then : a1D + a2D + a3D = 1, so that the three features can be considered as the probabilities of each point to be labelled as 1D, 2D, or 3D. It will help us to select the most appropriate neighborhood size by finding which radius favors the most one dimensionality (see Equation 3).

Optimal neighborhood radius
As presented in Section 2, the dimensionality features are computed for increasing radius values between a lower bound and an upper bound (rmin and rmax, respectively). Their set up is presented in Section 4.2. They represent the minimal and maximal acceptable neighborhood radius according to the area of interest and the sensor. The [rmin, rmax] space has been sampled in 16 values, which is a suitable trade-off between tuning accuracy and computing time. Since the radius of interest is usually closer to rmin than to rmax, the r values are not linearly increased but with a square factor. This allows to have more samples near the radius of interest and less when reaching the maximal values.
Two radius-selection criteria have been developed, namely the entropy feature (E f ) and the similarity index (Si). First, a measure of unpredictability is given by the Shannon entropy of the discrete probability distribution {a1D, a2D, a3D}: The lower E f (V r P ) is, the more one dimensionality prevails over the two other ones. The relevance of scale selection is demonstrated in Figure 3. This criterion allows to define an optimal radius r * E f that minimizes E f (V r P ) in the [rmin, rmax] space : A dimensionality labelling is then provided by : is defined as the ratio of neighbors P k which dimensionality labelling is the same as P at scale r: where 1 {.} is the Indicator function and n is the number of points within V r P . Si evaluates the homogeneity of the labelling within V r P by measuring the labelling similarity between neighboring points. Finally, we have: Si(V r P ).
Si aims at reducing the noise of the results since r * Si is the scale for which the labelling of the current point is the most similar to the labellings of its neighbors at the same scale. Both criteria have been tested for our datasets. Results and conclusions are presented in Section 5.1. Figure 3: Illustration of the relevance of the entropy feature for scale selection (building roof in ALSG). Left: 3D point cloud (height colored). Right: 3D point cloud colored with E f computed for one point of interest (pink dot). The color of each neighbor P k corresponds to the E f value, computed on the smallest neighborhood containing P k . The entropy first decreases until the neighborhood reaches the edge of the roof (blue circle -optimal size). Then, the entropy increases and becomes maximum when the neighborhood gathers distinct objects (red circle), here the ground and a chimney. Figure 4 provides an example of the behaviour of the shape fea-tures and the radius-selection criterion for two areas of the TLS dataset. Figure 4c corresponds to an area with higher noise than in Figure

Bounding scales
The optimal neighborhood radius is retrieved between predefined lower and upper bounds. They depend on various characteristics, and are therefore specific to each dataset. The choice of the lower bound is driven by (1) the point cloud noise; (2) the sensor specifications; and (3) computational constraints : 1. Local scattering may appear for linear or planar surfaces.
It may stem from sensor noise or unfavorable geometrical configurations (e.g., laser beam nearly parallel to a surface). Shape features are then biased and the point cloud may locally be erroneously labelled as "volumetric". The knowledge of the intensity of such noise is a prerequisite for solving this issue. Alternatively, the minimum size can be estimated applying the methods mentioned in Section 1 for this purpose.
2. The laser beam deflection system or the kind of mapping device may create point clouds with very irregular landscape sampling (low values in one direction and larger values in the orthogonal one). This is particularly true for MMS since datasets are acquired line by line with forward motion of the device. To cope with this issue, the rmin value is increased until V r P blends points belonging to at least two different scan lines.
3. A minimal number of points is required for the PCA. We consider relevant statistics start with 10 points. If the two first above-mentioned issues are not problematic, the point density directly provides the lower bound.
The selection of the upper bound is not critical and can be tuned with the knowledge of the lower frequencies of the relief i.e., the size of the largest objects in the scene. For TLS and MMS datasets, this corresponds to facades (typically 3 m) whereas for ALS data, ground regions and large buildings are involved. As E f and Si values may remain constant for these large planar areas, a cut-off value around 5 m works well, even if in practice 3-4 m are sufficient.

Comparison between radius-selection criteria
Both criteria E f and Si have been tested, and radius selection results can be significantly different. The entropy feature criterion E f has been preferred to the similarity index Si for two main reasons. First, as illustrated in Figure 4, E f and Si may provide distinct optimal radius that may correspond to distinct behaviours and privileged dimensionalities. For instance, in Figure 4c, Si erroneously selects a radius for which the volumetric dimensionality is predominant whereas E f allows to detect the correct planar behaviour of the underlying surface. Si favors the most homogeneous scale, but, in practice, several radii may provide the same result (Si = 1). The smallest one is kept which has been revealed to be problematic for the TLS dataset since the optimal radius corresponds to the noise level. For ALS and MMS, such problem of erroneous dimensionality labelling does not occur. However, small radius are favored for large objects since Si may saturate, preventing the selection of larger radius. Consequently, the choice of the lower bound rmin for the multi-scale analysis becomes critical. Second, for each 3D point at each radius of interest, Si depends on the dimensionality labelling of the neighbors at this neighborhood size (i.e., E f ). Two iterations are therefore necessary to retrieve the optimal radius which requires more computational time and memory.

Scale selection
The six datasets have been processed with the proposed approach. Scale selection results are displayed in Figures 4 and 5, using the E f criterion. The expected behaviours are reached showing the efficiency of the multi-scale analysis using E f : • Planar areas exhibit high values for ALS and TLS datasets (see Figures 4b and 5a). This is particularly true for the TLS dataset since the multi-scale analysis allows to overtake the noise level.
• Lowest scales correspond to border areas, allowing to retrieve building and vegetation edges for ALS datasets, windows in MMS datasets (Figures 5c and d), or small architectural elements in TLS.
• The algorithm allows to select high values when no predominant dimensionality is found for small and medium radius values. A sufficient number of 3D points is collected to retrieve a coherent dimensionality : 1D for wires and poles (Figures 5c and d), 2D for building facades in ALS datasets (e.g., highest values in Figure 5b which correspond to a few number of points lying on the facades of high buildings) or 3D for vegetated areas (all datasets).
• The method deals well with varying point densities. Such variation may come from the number of overlapping strips for ALS datasets (a lower radius size is necessary when two  The effectiveness of the scale selection process can be assessed by retrieving the predominant dimensionality that has been retrieved (1D, 2D or 3D ?) for various objects of interest. Several examples are presented in Figure 5, and comparisons with manually labelled ALS and MMS data have been performed (cf. Figure 6). For both datasets, planar objects are correctly retrieved, with few errors corresponding to building edges or small urban items, especially for ALS data. Besides, objects with one privileged dimension such as poles, wires or trunks are mainly labelled as 1D (M M SU dataset, in Figure 6b). However, since these objects are in fact cylindrical with a small width (e.g., trunks or traffic lights), they may look planar or volumetric for medium-sized neighborhoods. Conversely, volumetric objects such as trees may locally look planar. This happens for large objects low point densities or when no multiple scatterings are found with vegetated areas. Such phenomena is limited for genuine 3D acquisitions (TLS or MMS datasets), whereas it is clearly enhanced for 2.5D data such as for the ALSR data (Figure 6a). High vegetation areas are principally labelled as planar, which is due to the dense canopy cover. Nevertheless, as displayed in Figure 5a, such effect almost disappears when 3D points are acquired within tree canopies.

Comparisons with constant neighborhood size
The adaptive size strategy allows to retrieve a correct number of points to estimate the local dimensionality of the point set. Such strategy may be compared with the results achieved with neighborhoods of fixed size for a whole dataset. Firstly, the three shape features are computed with both strategies for four classes of interest of the ALSR dataset (Figure 7). One can see that the planar behaviour of ground and building points is improved with the adaptive size. Besides, for vegetated areas, the volumetric behaviour is slightly enhanced, however this happens in conjunction with an increase of the planar behaviour (dense canopies).
Figure 7: Improvements in shape feature computation using an optimal neighborhood radius per point. adD and a * dD are the shape feature for dimensionality d, computed with neighborhoods of constant size and adaptive size, respectively. The constant neighborhood size corresponds to the 30 nearest neighbors of each point.
Secondly, in addition to improved shape features, Figure 8 shows that the privileged dimensionality is also better retrieved with the adaptive strategy (MMSF dataset). This is particularly true for trees (1D→3D), small building facades elements (1D→2D), and tree trunks (2D→1D). Furthermore, since the geometrical analysis is less affected by the 3D point distribution, the labelling procedure is far less noisy.
Figure 8: Dimensionality labelling for the MMSF dataset using a constant neighborhood size for the whole point cloud (left), and the optimal value per point selected with E f (right). The constant neighborhood size corresponds to a 1 m radius.

CONCLUSIONS AND FUTURE WORK
The objective of this paper was to provide a simple and general method for finding the optimal neighborhood radius for each 3D point of a lidar point cloud, with the simple knowledge of the point locations. The most appropriate scale of analysis has been automatically selected among a range of potential values using three low-level geometrical shape features describing the "dimensionality" of the local point set (1D-2D-3D). Results on distinct lidar point clouds acquired with airborne, terrestrial and mobile mapping systems under various conditions showed the effectiveness of the proposed method. No a priori knowledge or assumption on the point cloud distribution, density, laser scan pattern, object size or underlying structures was required to achieve such results. Several experiments demonstrated how favorably our method performed compared to constant neighborhood sizes, and how relevant such approach is. Finally, the prevailing dimensionality and the scale of interest for the underlying object are found for each point, providing interesting cues for many potential subsequent segmentation or classification algorithms.
Improvements will focus on two main issues. Firstly, as illustrated in Figure 3, in addition to the global minimum of the entropy feature providing the optimal neighborhood radius, the other minima are likely to be relevant for multi-scale characterization. Secondly, each 3D point is analyzed separately resulting in noisy dimensionality labelled point clouds. Therefore, such results would be improved and fastened by taking into account for close points the scale already computed for optimal neighbors. Further work will consist in detecting and characterizing boundaries between objects by using and quantifying variations in the entropy feature. Finally, segmentation and classification issues will be tackled in order to fully benefit from the computation of non biased shape features and from the retrieval of a correct neighborhood for each point.