WEIGHTED ICP POINT CLOUDS REGISTRATION BY SEGMENTATION BASED ON EIGENFEATURES CLUSTERING

Dense point clouds can be nowadays considered the main product of UAV (Unmanned Aerial Vehicle) photogrammetric processing and clouds registration is still a key aspect in case of blocks acquired apart. In the paper some overlapping datasets, acquired with a multispectral Parrot Sequoia camera above some rice fields, are analysed in a single block approach. Since the sensors is equipped with a navigation-grade sensor, the georeferencing information is affected by large errors and the so obtained dense point clouds are significantly far apart: to register them the Iterative Closes Point (ICP) technique is applied. ICP convergence is fundamentally based on the correct selection of the points to be coupled, and the paper proposes an innovative procedure in which a double density points subset is selected in relation to terrain characteristics. This approach reduces the complexity of the calculation and avoids that flat terrain parts, where most of the original points, are de-facto overweighed. Starting from the original dense cloud, eigenfeatures are extracted for each point and clustering is then performed to group them in two classes connected to terrain geometry, flat terrain or not; two metrics are adopted and compared for k-means clustering, Euclidean and City Block. Segmentation results are evaluated visually and by comparison with manually performed classification; ICP are then performed and the quality of registration is assessed too. The presented results show how the proposed procedure seem capable to register clouds even far apart with a good overall accuracy.


INTRODUCTION
Dense point clouds can be nowadays considered the main product of UAV (Unmanned Aerial Vehicle) photogrammetric processing. They were largely studied in literature from several points of views such as algorithms and strategies for their generation, quality assessment, both in terms of accuracy and precision, segmentation techniques. Clouds registration is also a key aspect in case of blocks acquired apart; that usually happens when large datasets, covering wide areas, are acquired with several UAV mission or when time-series are considered. In the former case, data can be processed following a single block strategy and therefore a cloud is generated for each of them, in the latter one, clouds are obtained at different times; in either case, clouds must be registered to each other. Besides, similar issues can be observed in data acquired by Terrestrial and Aerial Laser Scanning systems where some misalignments can be present between scans. The authors have already worked on this topic and have recently published a study about . The work analysed the geometric consistency of two overlapping datasets, acquired with a multispectral Parrot Sequoia camera above some rice fields. The blocks were processed within Pix4D software package following different strategies. One of them concerned the Direct Georeferencing (DG) of each single photogrammetric block using the information registered by the camera GPS receiver. Since this is a navigation-grade sensor, the georeferencing information is affected by large errors and the so obtained dense point clouds are significantly far apart. The paper focused on their geometric consistency by exhaustively evaluating the distance between the generated point clouds using the Iterative Closes Point (ICP) technique.

* Corresponding author
As is well-known in literature (Besl and McKay, 1992;Chen and Medioni, 1992;Toldo et al., 2010), ICP is a procedure aiming to align point clouds without requiring the identification of homologous points. It starts by associating each point of a cloud to its closest point belonging to another cloud. The obtained coupled points are then used to estimate a coordinate transformation, typically a roto-translation, having six parameters, also known as a rigid body transformation (Low, 2004). The procedure is iterated until the latest estimated transformation is negligible. The process is strongly affected by point cloud shape because certain types of geometry can lead to an instable solution which means that the minimizing transformation is not unique (Gelfand et al., 2003). A common example is constituted by two planes that are parallel to each other and to the xy-plane. Once the planes are aligned, there are still three degrees of freedom: a relative translation between them in the xy-plane and a rotation around the z axis. This example is not too far from a real case in which a flat area was surveyed with an UAV, like the rice fields of our test-site. The choice of a subset of points is commonly adopted for improving final solution stability since the selected features have characteristics suitable to solve this rank deficiency. However, the correct points selection heavily affects the final alignment accuracy, therefore this step represents a crucial task (Glira, 2015). In the previous paper, a double density points subset is proposed to solve this task where the choice of the density level is connected to the terrain characteristics: lower in flat areas and a higher in variable terrain where ditches or dirt roads embankments are present. This approach reduces the complexity of the calculation, avoids that flat terrain parts, where most of the original points, are de-facto overweighed and favour the ICP convergence. However, the construction of this double density structure was realized with a manual approach limiting the extensive use of the proposed procedure to large dataset. In the past years, several authors have instead proposed and compared automatic sampling strategies having different levels of complexity: Masuda and Yokoya (1995) have evaluated the benefits of a random sampling while Rusinkiewicz and Levoy (2001) have proposed a technique in which points are selected in such a way that the distribution of their normals, in angular space, is as uniform as possible. Gelfand et al. (2003) have extended the last sampling strategy (Rusinkiewicz and Levoy, 2001) considering the eigenvectors of the covariance matrix. Since each eigenvector corresponds to a main motion, that can be described as a rotation around an axis and a translation along that axis, the analysis of theirs values can allow to evaluate stability of the transformation. Based on the spatial information of all 3D points within the local neighbourhood, invariant moments representing geometric properties can be calculated for each of them (Maas and Vosselman, 1999). The eigenvalues can directly be used to describe the local 3D structure or, alternatively, further measures based on these eigenvalues can be derived which encapsulate special geometric properties such as linearity or planarity (West et al., 2004;Mallet et al., 2011). These geometric descriptors, called eigenfeatures, can then be used to identify points useful for ICP registration. This strategy has the advantage to be applicable to any type of clouds be it produced by photogrammetry, as in the present paper, of by laser systems. Within this work, the eigenfeatures are used to characterize all the points and, starting from these descriptors, relevant ones are extracted from the original clouds in order to be used for ICP registration. Once again, a double density points subset is created. In this new realise, the choice between low-and high-density areas is done not manually but fully automatically thanks the use of eigenfeatures. K-means technique is used to group them in two cluster according to their geometry, flat terrain or not; besides, two distance metrics are compared, Euclidean and City Block. Finally, ICP registration is performed using the so-obtained segmentation.

METHODOLOGY
The proposed methodology for the point cloud registration consists of four steps procedure: eigenfeature extraction and selection, k-means clustering, point cloud segmentation and weighted ICP registration. All the steps are implemented in Matlab, realise 2019b.

Eigenfeatures extraction and selection
The adequate choice of a neighbourhood for determining the eigenfeature values of each point, depends on the characteristics of the cloud data especially to its points density and 3D shape. The choice can be based on a-priori definition of the search area in terms of radius or number of points (Friedman et al., 1977;Arya et al., 1998), or adapting this parameter according with the local geometry of the point cloud (Weinmann et al., 2015c;Farella et al., 2019). While the former requires an empiric knowledge of the scene, the latter is more versatile because it is not restricted to a specific dataset. The procedure implemented in this paper follows the first strategy and fixes a constant search radius because the area is flat, almost comparable to a surface, and the density is substantially uniform. This means that for each point belonging to the cloud, a list of neighbours, falling in search radius, can be associated. Then, for each 3D point and its neighbours, the derived normalized eigenvalues with = 1,2,3 can be extracted using the Principal Component Analysis (PCA). These values, obtained from the covariance matrix, represent the variation of the points distribution along the three principal orthogonal directions. Eigenvalues can be combined to obtain some shape descriptors called eigenfeatures (Weinmann et al., 2015a;Farella et al., 2019)  Eigenfeatures extraction was easily implemented in our modules thanks to the use of geoFEX Matlab toolbox (GeoFEX toolbox), developed by the Institute of Photogrammetry and Remote Sensing in Karlsruhe, Germany (Weimann et al., 2015b). Even if all the eigenfeatures have been extracted, it must consider they may contain redundant information. Besides, some of them can be substantially irrelevant; linearity for instance, which expresses the local similarity of the cloud to linear elements, do not contribute significantly in case of flat terrain. In our specific case, eigenfeatures capable to identify the presence of elements of discontinuity, such the change of curvature, could contain the information useful for a reliable registration. As highlighted by some authors (Weimann et al., 2013;Roffo, 2016), it is often desirable to select a compact subset of the most relevant features which allows for classification/clustering without significant loss of information. As we have chosen an unsupervised approach, eigenfeatures selection is a difficult task due to the absence of class label that could guide it. Among unsupervised selection method, we have adopted the Laplacian Score (LS), proposed by He et al. (2005), which is based on the observation that data from the same class are often close to each other, therefore the importance of a feature is evaluated by its power of locality preserving.

K-means clustering
K-means clustering is an unsupervised method that aims to subdivide observations into clusters; a cluster refers to a collection of observations aggregated together according to certain similarities. Each observation is allocated to each cluster through reducing the inner distances; in other words, k-means identifies centroids and allocates each observation to the nearest one, while keeping cluster as small as possible (Jain, 2010). The algorithm aims at minimizing an objective function, in this case a squared error function, that is (Bora et al., 2014): The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) (1) where and are the cluster and observation numbers, respectively, and is the centroid of the j-th cluster. The observations ( ) are vectors which can contain descriptors. For a point cloud, they can be composed, for instance, by six elements, 3D position and RGB values; in the proposed strategy, these vectors contain the set of selected eigenvalues (Section 2.1). For our specific dataset, information about colour is not used because data concerns some blocks acquired with a multispectral Parrot Sequoia camera (see Section 3 for more information); with this sensor, point clouds are produced using the green band only and this radiometric information does not contribute significantly to clustering. The objective function can be represented by different distance metrics among which Euclidean and City Block, both tested in the present paper. Euclidean, or squared Euclidean, distance between two observations, and , with dimension is calculated as (Deza and Deza, 2013;Bora et al., 2014) = √∑( − ) 2

=1
( 2) while City Block (Manhattan) distance is defined as (Deza and Deza, 2013;Bora et al., 2014) If we considered two points in the xy-plane, the shorter distance between them is along the hypotenuse, which represents the Euclidean distance. The City Block distance is instead calculated as the sum of the distances in x and y direction, which is similar to the way people move in a city, like Manhattan, where it is possible to walk around the buildings but not through them. Finally, the number of clusters is not always known in advance and some methods can be used to find it. Among them, there is the silhouette coefficient (Kaufmann and Rousseauw, 1990;Lletı et al. 2003) that is capable to measure how similar an observation is to its own cluster compared to the other ones; an high value indicates that the observation is well matched to its own group and the cluster numerosity is appropriate. In our case, k-means is used to subdivide the points in clusters that have the same shape characteristics. As explained in Section 1, data nature can significatively influences ICP registration: flat terrain could insert instability while discontinuities, such as ditches or escarpments, may instead facilitate algorithm convergence. Clustering procedure allows to separate these two typologies of terrain thanks to the use of eigenfeatures (that are shape descriptors); for our purposes, should then be set equal to 2. Nevertheless, as reported in Section 4.2, a preliminary silhouette analysis is performed in order to confirm the correctness of our hypothesis.

Point cloud segmentation
As reported in the Section 1, ICP is an algorithm capable to register overlapping points clouds through an iterative procedure that minimize an error metric at each step. The registration does not require the identification of any homologous points because it is based on the association of each point of a cloud to the closest point belonging to another cloud; a coordinate transformation is then estimated and applied to one set of points. The procedure is iterated until the mutual distance between the two clouds is minimized. However, certain types of geometry can lead to a rank deficiency which means that the minimizing transformation is not unique; this happens for instance when two parallel planes are taken into consideration such as UAV blocks acquired on flat areas. This corresponds to our intuitive notion of three degree of freedom: a planimetric translation and a rotation around the vertical axis. The choice of a suitable subset of points is commonly adopted for improving the stability of the solution. In the proposed strategy, the subset is chosen through a subsampling in compliance with the clusters obtained at the previous stage. Indeed, the clouds are segmented in two regions characterized by a different point density: lower in flat areas (cluster #1) and a higher in variable terrain (cluster #2). A structure, called skeleton in the following, is then obtained by overlapping a 2 m width grid to the clouds and down-sampling them according to k-means results: the cells lying on flat terrain is set to a density of 1 pt/m 2 , therefore the spacing is 1 m, while, the others, which lie where there are ditches and escarpments, is set to 64 pt/m 2 , with a spacing of 0.125 m. By imposing a suitable threshold on the ratio between the two cluster in each cell, the two classes, low and high density, are decided quite effectively. The points horizontal positions are established following a regular grid while their vertical components are estimated by interpolation. The use of such structures introduces a threefold advantage: improves the algorithm stability (that is the main task), decreases the complexity of the ICP calculation, avoids overweighting of flat areas, and reduces point clouds noise thanks to interpolation. The skeletons so produced are then used for next step, ICP registration.

Weighted ICP registration
Once the skeletons are created, they are used to estimate the ICP transformation. As we do not use the original point clouds, but a subset of them, we can say that we perform a weighted estimation because the points belonging to discontinuities have a larger weight, due to their numerosity, respect to those lying in flat area. This allows a larger stability to the final transformation. First, each point of one skeleton is coupled with the closest points belonging to the other one using a k-d tree approach; then, an outlier rejection is performed based on the points' mutual distances, on colour difference, and on angle between their normal difference. Therefore, the selected couples are used to estimate a 3D rigid-body transformation based on a point-toplane metric (Chen and Medioni, 1992;Low, 2004) whose formulation is where is the generic point of the skeleton ; is the correspondent point of the skeleton , derived by nearest neighbour searching; is the normal vector at point ; is the 4 × 4 transformation matrix estimated from previous iteration; is the 4 × 4 transformation matrix estimated during the current iteration. The process is stopped when the latest estimated transformation is negligible. The complete description and full flowchart for the implemented ICP is reported in .

The equipment
The dataset was acquired with the HEXA-PRO™ UAV, which is operated by the Laboratory of Geomatics of the University of The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) Pavia and is shown in Figure 1. The vehicle was made by a small Italian company named Restart® and has the following main characteristics: 6 engines (290 W each one), Arducoptercompliant flight controller, maximum payload of 1.5 kg (partly used by the gimbal, weighting 0.3 kg), flight autonomy of approximately 15 min. The UAV was equipped with a Parrot Sequoia camera (see Figure 1c). Sequoia has a high-resolution RGB camera with a 4608 × 3456 pixels sensor, a pixel size of 1.34 μm, and a focal length of 4.88 mm; the ground sampling distance (GSD) is 1.9 cm at 70 m height above ground level (AGL). Sequoia also has four monochrome cameras that are sensitive to the following spectral bands: green (G, 530-570 nm), red (R, 640-680 nm), red-edge (RE, 730-740 nm), and nearinfrared (NIR, 77-810 nm). Their resolution is 1280 × 960, with a pixel size of 3.75 μm and a focal length equal to 3.98 mm; the GSD is 6.8 cm at the 70 m flying height (AGL), which was adopted for the described survey.

The blocks structure
On September 13, 2017, a photogrammetric survey was performed on the Santa Sofia farmstead, near Pavia, Northern Italy. The test-site is a flat area of about 36 ha, used exclusively to cultivate rice. The whole acquisition was obtained by five flight missions, the outlines for which is shown in Figure 2, where the optical orthomosaic, which was used as background, was derived from a previous survey. In total, the project constituted about 1300 multispectral images, each composed of four bands. The AGL height was 70 m and image overlapping was 80% and 60% along-and across-track, respectively. The present paper will only focus on flights 3 and 4, as it has a methodological purpose. The former is composed by 293 images while its extension is 12 ha; the latter has 226 images and covers an area of about 9.5 ha. The outline of the overlapping area between the two flights is highlighted in red in Figure 2.

The photogrammetric processing
The photogrammetric project is carried out with Pix4Dmapper Pro, version 4.4.9. Since the original paper  deals with precision farming applications, only the four multispectral channels were considered, having 6.8 cm GSD. The higher resolution RGB imagery is then disregarded, as it is recorded in the JPEG format with a high compression factor, and has low quality compared to photogrammetry requirements.
Since the paper faces a method for point clouds registration, only Direct Georeferencing (DG) is considered as it is the most disadvantageous configuration. No Ground Control Points (GCPs) are inserted and images geolocation is obtained from the GPS receiver integrated with the camera; information contained in the EXIF files is instead adopted as internal parameters. As we knew from the Pix4D technical support, the parameters delivered into the EXIF are individually determined for each item at the factory. Their reliability is good, as reported in (Fernández-Guisuraga et al., 2018), in which the changes between nominal and optimized camera parameters were as low as 0.01%. The processing followed the usual pipeline Casella et al., 2020) image alignment, tie point extraction, bundle block adjustment (BBA); the two blocks are processed independently. Dense point clouds are generated using half image size resolution and obtaining an average density between 11 to 14 points per m 3 . In a preliminary test, the original image size resolution was also evaluated, but higher point density did not significantly improve the generation of orthophotos and reflectance maps.

RESULTS
The section describes the results obtained using the proposed strategy for point clouds registration. Since the paper has a methodological approach, only a subset of the available photogrammetric flights is taken into consideration: Blocks 3 and 4 ( Figure 2). For the same reason, only point clouds generated using DG configuration are processed because it represents the most disadvantageous configuration. Since the Parrot Sequoia sensor has a navigation-grade receiver on-board, the georeferencing information is affected by large errors and the so obtained dense point clouds are significantly far apart. Figure 3 shows the two generated green band orthophotos in which a significant planimetric shift is clearly present between the two blocks. Another suitable method to evaluate clouds distance is the use of profiles: Figure 4 highlights the position of an East-West longitudinal profile in the overlapping area. The two extracted profiles, 1 m thick, are reported in Figure 5 where blue and red lines represent Block 3 and Block 4, respectively. The image depicts the existence of a z-shift and a significant rotation. Some preliminarily steps are done before ICP registration: the common enveloped area is determined, and the two clouds are trimmed according to it. Besides, a further precautionary buffer is added to avoid edge effects. Clouds are then processed The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) following the steps previously described (Section 2): eigenfeatures extraction and selection, K-means clustering, point cloud segmentation and weighted ICP registration.

Eigenfeatures extraction and selection
First step of our procedure is the extraction of the eight eigenfeatures for each point of the two clouds (Section 2.1). A constant search radius is fixed inside which neighbours are identified; radius is set equal to 1 m, according to area characteristics and point clouds density, obtaining an average of 51 neighbours per point. Covariance matrix and eigenvalues are then determined, and eigenfeatures are calculated using the formulas reported in Table 1. Eigenfeatures are then normalized to the interval [0,1] and stored in a matrix having as many rows as the number of points, and eight columns. Laplacian rank method (He et al., 2005) is then applied to select main relevant features. As the final choice must be suitable for both point clouds, the feature rank is estimated three times: for the first cloud (block 3), for the second one (block 4) and for the join of the two datasets (block 3-4). The so-obtained classifications are then averaged to draw up a global ranking, as shown in Figure 6; the bars represent the mean of the ranking positions, so the lower values indicate high positions in the meaningfulness rank. Two features, scattering and change in curvature, stand out to be more relevant than others, such as linearity or anisotropy. The outcome is substantially coherent with the characteristics of the test-site which is mainly a flat terrain with homogenous density and highlights those eigenfeatures that can help to find shapes useful to ICP convergence. Moreover, the distance between the top two eigenfeatures and others is large enough to assume that they contain all the information needed for the next stages; for this reason, only scattering and change of curvature are selected for k-means clustering. Figure 6. Mean rank of the eight eigenfeatures

K-means clustering
K-means clustering uses an iterative algorithm that assigns observations to clusters so that the sum of distances from each of them to the centroid is a minimum; the procedure returns the cluster index for each observation. For our dataset, observations are constituted by the two eigenfeatures, scattering and change of curvature, selected at the previous stage; this information is then used to cluster the point cloud in groups suitable for clouds registration. Even if, theoretically, should be set equal to 2 in order to separate flat areas from terrain having shape characteristics useful for ICP registration, such as ditches or escarpments, a preliminary silhouette coefficient analysis is performed. This method determines how well each observation lies within its cluster providing a value for each of them. These values range from -1 to 1; a high silhouette value indicates that a point is well matched to its own cluster, and poorly matched to other clusters. The optimal number of clusters is the one that maximizes the average silhouette over a range of possible values for . Figure 7 shows the results obtained up to 10 clusters where best value is reached for = 2. Besides, the curve also shows as silhouette coefficient has a sudden decrease from two to three clusters and descends constantly from that point on. This behaviour supports our initially hypothesis of using two clusters. Once established the clusters number, K-means clustering is run for both point clouds using the two metrics, Euclidean and City Block. The next figures show the results of the clustering: red dots represent cluster #1 identify as the flat terrain; blue dots are instead points belonging to discontinuities such as ditches or dirt roads embankments (cluster #2). Figure 8 and Figure 9 report results obtained using Euclidean metric for Blocks 3 and 4, respectively; Figure 10 and Figure 11, those found for City Block distance.  Figure 12 shows instead the orthophoto produced using only the green band imagery for Block 3. In the figure it is possible to identify the location of the elements of interest: the dirt roads (the cars used for surveying are also distinguishable) and the irrigation canals. From a first visual comparison, it is evident that there is a good correspondence between these elements and cluster #2. However, the two tested metrics have a different capability to detect them: City Block shows a larger uniformity along these linear discontinuities while Euclidean distance seems to identify them more irregularly. Moreover, both metrics classify some points lying inside cultivated fields as belonging to cluster #2; this happens especially for City Block distance. This behaviour could be related to the presence of weeds inside the planted crops; indeed, this vegetation, taller, could be identified as element of interest and then classified in the second cluster; the result must not be necessarily considered like a mistake, if it contributes to ICP convergence. Nevertheless, in this case, it seems to be more connected to the presence of noise in the point clouds. If we compare, for instance, the upper left corner of Figure 10 and Figure 11, it will be evident as the two images, even if referred to the same area, presents a different clustering.
Remembering that the two blocks have been independently processed, using the positions coming from the camera on-board GPS receiver, this disparity could due to a different quality of the external solutions which would cause, for Block 3, a noisier point cloud. Even if less manifest, the same phenomenon is also present in Figure 8 and Figure 9, obtained with Euclidean distance. The best result obtained by City Block metric could be explained considering that the researched elements (dirt roads and canals) have a structured quite similar to urban streets; both features are characterized by linear elements, orthogonal to each other. In the next step, only the results obtain by k-means clustering, based on City Block distance metric, will be used as starting point to create the double density structure needed to ICP registration.

Point cloud segmentation
A 2 meters width grid is superimposed on the dense clouds to perform segmentation. Each cell is down-sampled according to clustering results: the cells lying on flat terrain, cluster#1, is set to a density of 1 pt/m 2 , while, the others, which lie where there are ditches and escarpments (cluster #2), is set to 64 pt/m 2 . The two classes are established by imposing a threshold on the presence of the two clusters in each cell; if at least the 25% of points lying in a cell are classify as cluster #2, the cell owns to the high-density class otherwise to the lower one. The threshold is chosen is favour of high-density class in order to be sure not to disregard any possible area useful for ICP convergence. Figure 13 and Figure 14 the so obtained results for Block 3 and Block 4, respectively. Light grey represents 1 pt/m 2 density cells whereas the dark grey the 64 pt/m 2 ones. Dirt roads and canals are well represented in both cases even if the presence of noise in some flat areas of Block 3 is still present, as already discussed in Section 4.2. Figure 13. Skeleton structure for Block 3 obtained by point cloud segmentation The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) Figure 14. Skeleton structure for Block 4 obtained by point cloud segmentation This asymmetry in the skeleton structures for the two blocks could lead to an instability on the ICP algorithm and, for this reason, a combined analysis of the two segmentations is made. A cell comparison is then performed between the two skeletons in order to reach comparable structures. This operation could be considered like that one performed in photogrammetry for image alignment: after the extraction of points of interest, only the key ones (those connecting the different images) are considered for next processing stage. In the same way, after the creation of the skeletons of the two clouds separately, only the common cells are considered for ICP registration. This strategy allows to remove the noise present in Block 3, as show in Figure 15: the skeleton is now cleaner, and the researched shapes (roads and canals) are more manifest. To evaluate the quality of the automatic procedure proposed so far, a comparison with the segmentation achieved thanks to human interaction  is performed. Figure 16 reports the skeleton obtained with manual segmentation; the shapes are regular, and no black spots are present in the flat areas. The visual comparison between Figure 15 and Figure 16, for instance, shows a significant accordance between them; automatic procedure seems to bring to results substantially comparable to the manual ones.
To quantify the accordance, the confusion matrix is also calculated: each cell of the manual segmentation (true value) is compared to the corresponding cell of the automatic segmentation (predicted value). Both blocks are considered, and a unique confusion chart is produced summarizing the results ( Figure 17). In the main diagonal of the matrix it is possible to visualize how many cells are identified in the same way by both the manual and automatic approach. There is a good accordance in the low-density areas: more than 95% of them are correctly classify considering both user and product accuracy. The quality decreases around 75-77% considering the high-density cells. This difference in the results must be ascribed more to a misclassification of the human operator than to an error of the automatic algorithm. Comparing Figure 15 and Figure 16, especially where the dirt roads intersect to each other, it is evident as the two classifications, automatic and manual, have obtained different results: low density for former and high for the latter. Logically speaking, the upper part of a road is flat, so the algorithm has worked well. Human operator has simplified the reality introducing in this way an error source. Nevertheless, the result confirms the good agreement between manual-automatic segmentation and the overall accuracy is 92.5%. Figure 17. Confusion chart for manual/automatic segmentation; HD stands for High Density whereas LD for Low Density

Weighted ICP registration
The skeletons produced at the previous stage are then used to estimate 3D rigid-body transformation to register the two point clouds through ICP algorithm. It is known that ICP in not always reliable, especially when it is used to register almost flat clouds, as in our case. The registration can be performed in two directions, moving the first cloud toward the second one and vice versa, and the parameters of these transformations should coincide, aside from uncertainties. We use the comparison between the estimated and calculated inverse transformation to infer the precision of our estimation. Table 2 shows the results, where T is the transformation to align Block 3 to Block 4 (column 3), inv(T) is its inverse (column 4) and T' is the transformation to register Block 4 to Block 3 (column 5). Last column shows the difference between the two previous ones. There is a good accordance between the parameters: the worst residual for translation component is 10 cm while for orientation angle is 0.04 deg. A distance-equivalent error ( ) can be computed for the resulting angular residual ( ) by assuming a distance ( ) of 200 m, corresponding to the half-width of the considered test area. By applying the simple formula = • , where the angle is expressed in radiants, it can be found that = 14 . Now, we must consider the granularity of the datasets (i.e., the points' linear spacing). For the skeleton, the spacing is 12.5 cm for dense parts and 100 cm elsewhere. As residuals of the transformation substantially equal the discretization size, we consider the estimated transformations reliable and precise. Moreover, to evaluate the quality of the ICP transformation the distance between the two point clouds are evaluated before and after the transformation. All the results are shown in Table 3 in which the columns report: the identifier of the scenario, before and after ICP, the statistical figures considered, the differences for the horizontal and vertical components. The initial distance between the two clouds presents large values, as expected, while ICP produces a good alignment between them with a 3D rmse residual of about 16 cm.  Table 3. Summary statistics of the distance between overlapping point clouds Figure 18. The cloud 3D distance maps before the ICP registration, between Blocks 3 and 4 in the overlapping area, expressed in meters Figure 19. The cloud 3D distance maps after the ICP registration, between Blocks 3 and 4 in the overlapping area, expressed in meters Maps of the 3D distances between the overlapping point clouds are meaningful. Figure 18 and Figure 19 shows them before and after ICP, respectively; remarkably, the figures shown adopt the same colour map. Figure 18 highlights that the original clouds are quite far; the map of the distances is a sort of a ramp, meaning that the two clouds are not simply displaced but are also affected by a significant rotation (already reported in Figure 4). Figure 19 reports the distance between datasets after ICP registration; the two clouds are now close to each other and their distance is 15 cm, in the average. The result has the same order of magnitude of the clouds density (Section 3).
In the alignment step, Block 4 is moved toward Block 3 and a new profile can be then extracted for the former. Figure 20 shows once again the extracted longitudinal profiles: blue line represents the profile for Block 3, red line the profile extracted for Block4 before the ICP registration and green line that one obtained just after the alignment. The two clouds are now correctly positioned; nevertheless, it is present a slight deformation in the left side of the profile that cannot be reduced with a simply rigid-body transformation. The same phenomenon is also visible in the North-East area of Figure 19. In this case, only a more careful inner calibration could decrease this misalignment due probably to original cloud deformation.

RESULTS AND FURTHER ACTIVITIES
The paper concerns the evaluation of an automatic procedure for point clouds registration. The procedure can be suitable to align datasets acquired by several surveying techniques such as UAV or laser scanners, terrestrial or aerial. In the paper two overlapping datasets, acquired with a multispectral Parrot Sequoia camera above some rice fields, are analysed in a single block approach. Since the sensors is equipped with a navigationgrade receiver, the georeferencing information is affected by large errors and the so obtained dense point clouds are significantly far apart. An innovative procedure, based on ICP algorithm, is then proposed for their registration. Eigenfeatures are extracted for all points and the most significant are then selected; they contain the geometric information needed to identify the terrain discontinuities useful for ICP convergence. Eigenfeature are then used to clusterize the clouds in two groups: flat area and variable terrain (dirt roads embankments and ditches). As the point clouds used in this paper derive from a multispectral camera, only the green band has contributed to their generation.
For this reason, no information about colour has been inserted into k-means clustering. In the further activities new datasets will be tested considering also images acquired by optical cameras in order to understand if these additional descriptors (red, green and blue channels) can influence and improve the results.
Two metrics are tested for clustering, Euclidean and City Block, but only the second one has shown results suitable for the next steps. This could be explained considering the characteristics of the researched elements that are mainly linear and orthogonal to each other like the streets of a big city. This hypothesis will be further evaluated in next activities taking into consideration other dataset having different terrain morphology. The clouds are then segmented into 2 meters width cells having two different point densities: 1 pt/m 2 and 64 pt/m 2 . The choice between the two classes is based on the clustering results: lower density for areas characterized by points belonging to flat terrain cluster, higher density elsewhere (cluster containing discontinuities). Finally, the so obtained structures are inserted as input in ICP algorithm and used to determine the transformation parameters capable to align the two clouds. The final alignment between the two clouds is good proving that the proposed procedure is reliable also to register datasets having challenging morphologies such as flat terrain. In the present work, the average 3D distance between the two clouds passes from more than 3 m to approximately 15 cm (equal to the discretization size). At the moment, the process is fully automatic except for radius search of eigenfeatures, sets manually equal to 1 m in the present paper. Further activities will also focus on this aspect evaluating different automatic strategies such as variable search radiuses or eigen-entropy. Finally, the described process is composed basically by four steps (the same used to identify the various subparts of Section 2 and 3): eigenfeature extraction and selection, k-means clustering, point cloud segmentation and ICP registration. The procedure, although has worked well, could be in our opinion furthered improved by simplifying the workflow.
A new procedure will be tested in which the clouds are firstly subdivided into cells and eigenfeatures are then extracted according to them. This approach would allow to clusterize directly the cells avoiding the segmentation stage.