## 2 Methods
### Geometric Consistency Assessment  
The iterative closest point (ICP) methodology was adopted to register overlapping point clouds, evaluate their distance, and estimate their geometric consistency. As is well-known in literature [37; 38; 39], ICP is a procedure aiming to align point clouds without requiring the identification of homologous points. It starts by associating each point of cloud A to its closest point belonging to cloud B, then a coordinate transformation (typically a roto-translation, having six parameters, also known as a rigid body transformation in literature [40]) is estimated, based on the obtained coupled points, and applied to one point set. The procedure is iterated until the latest estimated transformation is negligible. A dedicated Matlab procedure was specifically developed at the University of Pavia, implementing ICP and including some unique features. Procedure flowcharts are reported in Figure 3 and Figure 5.  
Preliminarily (Figure 3), the common enveloped area for the two point clouds is determined. The clouds are then trimmed according to this area, adding a further precautionary buffer to avoid edge effects; the buffer value is set manually. Therefore, a subset of points is extracted (the so-called skeleton) from each of the point clouds that have a variable density. The skeleton is constituted by squared meshes with sides measuring 2 m and belonging to two classes. There are skeleton points inside meshes. Those lying on flat terrain contain 1 pt/m\({}^{2}\), and therefore the spacing is 1 m. The others, which lie where there are ditches and escarpments, contain 64 pt/m\({}^{2}\), with a spacing of 0.125 m. The skeleton was adopted to reduce the complexity of the calculation and to avoid flat terrain parts, where most of the original points are de-facto overweighed. The classification of each mesh of the skeleton was performed by selecting all the original cloud points lying in the mesh and estimating the interpolating plane. By imposing suitable criteria concerning the residuals (low residuals mean flat areas) and the deviation from the vertical of the plane normal, the two classes (flat and variable terrain) were decided quite effectively. In the current version, the thresholds used for skeleton classification must be tuned by the operator and inserted manually; an improved one is under development, based on machine-learning. An example of the skeleton structure is reported in Figure 4 for sub-block 3, showing the skeleton points for flat and steep parts of the terrain. Figure 3: Flowchart for data preparation. After defining the skeleton as described above, we used it to define a new point cloud. By this, we mean a set of points for which we know the 3D coordinates, the normal vector of the surface at their position, and the color, totaling 9 descriptors. Original point clouds A and B will be referred with the acronyms PCA and PCB, respectively, while point clouds obtained from skeletons will be named skeleton point clouds A and B, shortened to SPCA and SPCB, respectively. Each element of the skeleton point clouds has 9 descriptors, as stated before; some were known from the definition and some were calculated. For each skeleton mesh, the fitting plane was used to estimate the height of the interior skeleton points. The plane's cartesian equation was also used to obtain its normal vector. Finally, the color was determined for each skeleton point by picking that of the closest original point. To perform quality assessment, data filtering, and further analysis, data used for each plane estimation was stored in a complex data structure. It is named cell array in the Matlab environment and can be thought of as a matrix where each element can store any kind of data structure. We created a cell array containing as many rows as the meshes constituting the skeleton and with 9 columns (the same number as the descriptors associated with each point cloud, but purely by chance--there is no relationship).