## 2 Methods
### Geometric Consistency Assessment
For each mesh, the associated 9 cells contain the planimetric coordinates of the central point; mesh size; mesh point density, where we counted the number of the cloud points lying in the mesh and divide this by the area; mesh planimetric bounding box; mesh edge, intended as the points defining the corresponding polygon; coordinates of the original cloud points lying in the mesh; interpolating plane's normal vector; parameters of the plane and its fitting goodness (the so-called Matlab gof (goodness of fit) data structure); and the descriptors of the skeleton points located inside. Once SPCA and SPCB were created, they were used to fit the ICP transformation. The procedure (Figure 5) is iterative; it starts with the original SPCA and SPCB and finally produces the parameters of a six-parameter 3D rigid transformation aligning SPCB to SPCA. Each iteration:  
* has in its input the running SPCB (produced by the previous iteration; in the first iteration, this is SPCB itself), the SPCA (remaining unaltered for all the process), and the running 3D rigid transformation determined so far;
* couples each point of the running SPCB with the closest point belonging to SPCA;
* performs outlier rejection based on the points' distances, the angle between the surface normals at the two points considered, and the norm of the difference between their RGB vectors;
* determines the parameters of a refinement of the 3D rigid transformation aligning running SPCB and SPCA, by solving a non-linear least squares (NLLS) problem defined by Equation (1); in plain words, the NLLS solver evaluates the distance between each B-point (belonging to running SPCB) and the plane passing through the paired A-point (belonging to SPCA) and the normal to SPCA; it determines the unknowns in order to minimize the sum of all the distances;  
Figure 4: The skeleton structure for sub-block 3. (**a**) Overview: the parts shown in light grey have a lower resolution of 1 pt/m\({}^{2}\); the darker parts have a higher resolution of 64 pt/m\({}^{2}\). (**b**) Detailed view: individual skeleton points are visible, and their different densities can be easily detected. * applies the determined 3D transformation to the running SPCB;
* composes the newly determined transformation with that received in the input;
* returns the updated running SPCB and coordinate transformation. The process is stopped when the latest estimated transformation is negligible. The above-described process was coded in a Matlab toolbox coded at the University of Pavia. Among other features, the NLLS problem is solved with a robust approach, based on the Huber method. Moreover, the procedure takes advantage of the k-d tree functionalities (k-dimensional tree) to speed up point coupling, which are available in the used Matlab environment [40]. The k-d tree engine is trained for SPCA, which was kept fixed throughout the procedure. The mathematical formulation [41; 42] of the estimation for the coordinate transformation is  
\[CPTopt\ =\ arg\min_{T}\sum\left(\left(CPT\cdot SPCB_{i}-SPCA_{i}\right)\cdot n _{i}\right)^{2} \tag{1}\]  
where \(SPCB_{i}\) is the generic point of the skeleton B; \(SPCA_{i}\) is the correspondent point of the skeleton A, derived by nearest neighbor searching; \(n_{i}\) is the normal vector at point \(SPCA_{i}\); \(CPT\) is the \(4\times 4\)  
Figure 5: Flowchart for ICP procedure. Note: K-D Tree, K-dimensional tree algorithm; MDL, KD tree model object; CTP, current transformation Parameters; TP, transformation parameters. 3D rigid-body transformation matrix estimated from previous iterations; and _CPTopt_ is the 4 \(\times\) 4 3D rigid-body transformation matrix estimated during the current iteration.