PREDICTING THE ACCURACY OF PHOTOGRAMMETRIC 3D RECONSTRUCTION FROM CAMERA CALIBRATION PARAMETERS THROUGH A MULTIVARIATE STATISTICAL APPROACH

Several tools have been introduced to generate accurate 3D models. Among these, Unmanned Aerial Vehicles (UAVs) are an effective low-cost tool to go beyond on-fields effort limits since they allow to fly over areas difficult to reach and to reduce the time needed to collect and process photogrammetric pictures as well. Combining their versatility with Structure from Motion (SfM) techniques efficiency has provided a widely accessible approach to generate accurate photogrammetric products. However, the outcome resolution and coherences also depend on sensor traits. Therefore, UAVs are usually equipped with low-cost non-metric cameras, with the consequent requirement for a calibration procedure to increase the final 3D models accuracy. Although several researchers have highlighted the strong impact of camera calibration parameters on the photogrammetric outcomes, their linkage has not been explored yet. This paper is aimed at investigating their relationship and to propose a novel predicting function of 3D photogrammetric reconstruction accuracy. Such function was estimated thanks to the application of the Principal Components Analysis (PCA) technique. Four photogrammetric UAV flight surveys provided the input data of PCA while an extra dataset was used to validate the results. Once PCA was completed, a synthetic index was proposed and the coefficient of determination was calculated between the index and error components. Synthetic indices values for the various datasets were applied as baseline to detect a predictive function able to assess the northern and eastern error components with a deviation of 0.005 m and 0.003 m, respectively. The proposed approach shows promising and satisfying results for predicting 3D models accuracy.


INTRODUCTION
Over the years, several tools have been developed to generate accurate three-dimensional textured models, useful in a wide range of professional and research applications. Among these, in the last decade, the Terrestrial Laser Scanner (TLS) has been considered as the reference standard to generate detailed terrain landform (Medjkane et al., 2018). Nevertheless, it is extremely expensive both in terms of technological equipment and acquisition/processing data time. Consequently, it cannot be efficiently applied in all conditions. Today, UAVs (Unmanned Aerial Vehicles) represent an effective low-cost alternative to overcome the limits imposed by TLS (Roşca et al., 2018;Manfreda et al., 2018). First of all, they are versatile and flexible, allowing flight quotas reduction with a consequent increment of input data resolution (Nex, Remondino, 2013;Capolupo et al., 2014); in addition, UAVs allow to achieve areas difficult to reach without endangering people (Capolupo et al., 2018); and, finally, their most eligible property: they are able to drastically reduce the time required to collect and process the input data (Nex, Remondino, 2013). The possibility to adapt the Ground Sample Distance (GSD) to the size of the object under investigation is a direct consequence of the above-mentioned properties (Capolupo et al., 2015). Moreover, their combined use with Structure from Motion (SfM) technique and Computer Vision (CV) has provided an essential approach to generate accurate Digital * Corresponding author: Alessandra Capolupo, Dept. of Civil, Environmental, Land, Construction and Chemistry (DICATECh), Politecnico di Bari, Via Orabona 4 -70125 Bari, Italyalessandra.capolupo@poliba.it Surface Models (DSM) and 3D Point Clouds. Neverteless, the resolution of the obtained 3D models is affected by other factors as well, such as the procedure applied to generate them and the characteristics of the integrated sensor of the camera mounted on UAVs (e.g. pixel size, focal length of the lens, etc.) (Cramer et al., 2017;Kraft et al., 2016). Several scientific contributions cooperated in designing a validated methodology suitable for facing the first issue (production of high-resolution outcomes) and generating products with a resolution comparable to the most established topographic procedures (Caroti et al., 2017;Saponaro et al., 2019b). For instance, Mesas-Carrascosa et al., (2015) and Manfreda et al., (2019) investigated the influence of the flight plan parameters defining the general requirements to be adopted to optimize the entire process chain. Conversely, Saponaro et al., (2019b), Padró et al., (2019), Rangel et al., (2017) explored the impact of georeferencing strategies on the obtained accuracy in order to minimize the data acquisition and processing efforts without losing quality in the final products. Lumban-Gaol et al., (2018) and Benassi et al., (2017), instead, analysed the importance of the rational parameterization during the orientation phase by applying a different dedicated environment and they evaluated the metric consistency of the achievable results. The second issue (cost of equipment and processing data time) was tackled by equipping the UAVs with a low-cost light digital camera, with the consequent requirement for a calibration procedure during the metric reconstruction phase in order to increase the final 3D models accuracy (Salvi et al., 2002). Although the camera calibration procedure is based on the detection of focal length (f), principal point offset (Cx, Cy) and lens distortion parameters (k1, k2, k3, k4, p1, p2, p3, p4), these variables are actually taken into account just in the field of precise photogrammetry. Radial (k1, k2, k3, k4) and decentring (p1, p2, p3, p4) distortions include the aberrations which influence images position (Fryer, 1996). In many other applications, such as CV works, just the focal length is considered (Remondino, Fraser, 2006). The calibration phase is really essential in the applications where aerial photogrammetric pictures are involved since they are strongly affected by the geometric instability and limited precision and accuracy of the camera (Warner, Carson, 1991). Such issues are even more relevant with digital commercial cameras, commonly mounted on the UAVs since, in their case, defining the interior orientation of the camera is extremely difficult and, often, impossible (Pérez et al., 2011). Therefore, calibration algorithms suitable for processing UAVs aerial photos have been introduced and implemented in the widespread photogrammetric software, i.e. Agisoft Photoscan Professional (Zhang, 2000). Therefore, as highlighted by Oniga et al., (2017), several scholars investigated the importance and the impact of camera calibration parameters on the photogrammetric outcomes, nevertheless the estimation of the model accuracy could be achieved has not been explored yet. Therefore, this paper proposes a novel algorithm predicting the accuracy of 3D models generated by UAVs -SfM technique, based on the analysis and on the combination of camera calibration parameters using multivariate and linear statistical techniques. Such methodology provides functions suitable for predicting the 3D model accuracy.

Study Area and Field Operations
The algorithm was benchmarked and tested on a coastline stretch of about 400 m located south of Bari, in Puglia (Southern Italy) ( Figure 1). The entire coastal strip is characterized by frequent cliff collapses, because of which the site appears highly vulnerable to any environmental phenomenon. In particular, hydrological risk is widespread and visible along the shoreline. Five flights were carried out over the experimental pilot on December 2018, January 2019, February 2019, March 2019 and October 2019, respectively. All of them were performed using a commercial quadcopter (DJI Inspire 1), equipped with DJI ZenMuse X3, a non-metric camera, which is characterized by a focal length of 3.61mm, a pixel size of 1.56μm with an effective number of pixels equal to 12.4M. Moreover, even a low-cost GNSS/INS positioning system was collocated on the UAV to accurately achieve the waypoints during the survey campaigns. The flight was planned using iOS app DJI Ground Station Pro, setting the flight quota equal to 100m Above Ground Level (AGL) and the cruising speed to 5.5m/s. Each plan was composed of ten strips, perpendicular to the shoreline, and 77 waypoints, ensuring a longitudinal and transversal overlap of 80% and 70%, respectively. Moreover, each survey was performed under clear sky conditions and following the same flight pattern. 77 pictures were totally acquired in each survey. Flight plan parameters were fixed in order to obtain an expected GSD of about 4.3cm/px, suitable for describing the landform of the selected experimental areas.
After the first flight campaign, thirty fixed points, uniformly distributed on the investigated area (Figure 1), were measured using Global Navigation Satellite System (GNSS) Leica Viva CS10/GS10 receivers, exploiting the Network Real-Time Kinematic (nRTK) mode. Thus, their accuracy was equal to 0.02 m along each axis. Subsequently these points were used as Ground Control Points (GCP) or Check Points (CPs) according to the selected georeferencing strategy. All the collected data were organized in 5 dataset and separately processed according to the operative workflow described in Figure 2. The image blocks acquired in December, January, February and October were used to calibrate the models, while the remaining dataset taken on March was utilized to validate it.

UAV-imagery Processing
The image block acquired in each survey was treated as an independent dataset and, thus, separately processed in the same workspace. The photogrammetric reconstruction was performed using Agisoft PhotoScan software platform (v.1.4.1) (Agisoft LLC, St. Petersburg, Russia). The quality of the images of each dataset was improved removing the blurry photos. Thus, the quality images index was estimated through the tool Estimate Image Quality implemented in Agisoft Photoscan environment, which returns an average value greater than 0.8 for each dataset. The quality images index can result in a value between 0 and 1 (Agisoft, 2014): the highest the value, the highest the quality. Therefore, the 0.8 is a satisfying quality value (Agisoft, 2014).
To ensure the comparability among the results generated by the different dataset, a common workspace was set and RDN2008/UTM reference system zone 33N (NE) (EPSG:6708) was assigned to the whole workspace. Similarly, the accuracy of the on-board UAV equipment and GCPs coordinates were fixed. Camera Accuracy (m) option was set equal to 0.05m, meanwhile 10 degrees and 0.01 m were assigned to the Camera Accuracy (deg) option and to the Marker Accuracy (m) parameter, respectively. Conversely, in Image Coordinates Accuracy folder, the accuracy in pixels of the markers was parameterized with a value of 0.5 pixels, while the Tie Point Accuracy option was set equal to 3 px, as suggested in Saponaro et al., (2019a). In GPS/INS Offset option, the Lever-Arm vector was set equal to [0.005, 0.01, 0.25]m, considering a relative accuracy of 0.01m for each axis. Moreover, camera parameters correction option was enabled in order to improve their estimation proportionally to the photogrammetric block adjustments. Agisoft PhotoScan adopts the Brown's model (Brown, 1971) for describing camera lens parametrization, which is based on 13 parameters: adding to all the variables commonly applied in photogrammetric analysis, such as focal length, lens distortion parameters and principal point offset, the Skew parameters (B1, B2). These last ones are related to the affinity value and to the non-orthogonality coefficient, respectively. In Batch Process, "High mode option" was selected for the initial step of image alignments to optimize it and improve the quality of the obtained sparse point clouds. As enhanced by Saponaro et al. (2019a), the obtained points clouds were manually filtered to remove estimated points with a Reprojection Error value higher than 0.4. This step improves the compliance of the block by optimizing the camera calibration parameters estimation and the impact on the final accuracy values can be neglected. The GCPs were implemented in Agisoft Photoscan Professional to georeference the sparse point clouds, optimize the results of such phase and reduce the block deformation and systematic error (Gruen, Beyer, 2001). Once the alignment phase was completed, the leave-one-out technique was applied to extract 31 chunks from each dataset , characterized by a variable number of GCPs, from a maximum of 30 to a minimum of 0. Once all Bundle Block Adjustment (BBA) processes were completed, calibration camera configuration, estimated during the metric reconstruction of each chunk, was extracted and aggregated in a data table. They were considered as dependent variables and they were the input data of the subsequent processing phase consisting in simplifying the original set of data reducing the redundant information and extracting just the most relevant ones. Simultaneously, features were added to the obtained dense points clouds and a 3D textured model was generated.

Principal Component Analysis
Principal component analysis (PCA) is the oldest and the most popular multivariate statistical technique and, consequently, it is widely applied in almost all scientific disciplines in order to drastically reduce the large amount of input data. In fact, based on the variance maximization principle, PCA investigates the relationship among the several input variables, which are, in general, intercorrelated (Wold et al., 1987). After identifying the most significant information, it expresses the original data as a new set of linearly independent and orthogonal vectors, called principal components. The first principal component has the largest possible variance since it must explain the largest part of the inertia of the input data, while the second one must be designed orthogonal to the first component and include the largest possible remaining variance. All the other components must be constructed likewise to the second one (Abdi, Williams, 2010). In such way, the PCA technique compresses the number of original input data, removing the redundant information and maximizing their effectiveness (Bro, Smilde 2014). Thus, the principal components are extracted through the application of the Singular Value Decomposition (SVD) method, a generalization of the eigen-decomposition approach (equation 1). Following the annotations assumed in the equations: where j 2 = inertia of column j xij = element i of the column j The inertia provides an essential information since it reflects the importance of a specific component. Once the principal components were extracted, Kaiser's criterion was applied for determining the number of meaningful components to retain (Kaiser, 1960) and, consequently, all of them with an eigenvalue less than 1.0 were dropped. Thus, the weights related to each selected principal components were extracted as well, and all the original data were synthetized by computing the weighted average of the picked elements. Therefore, the result of such step was considered as a "synthetic index" representing the all procedure. This step as well as the further ones were performed programming a specific code developed in the open-source statistical software R.

Predictive Function
The relationships between the synthetic index and the northern and eastern error components of the metric reconstruction were explored through the implementation of the coefficient of determination (R 2 ). Examining the proportion of the variance in the dependent variable, R 2 measures the ability of the synthetic index to replicate the error components. Therefore, this step was essential to investigate the performance of the applied procedure. Its outcome, ranging between 0 and 1, was interpreted according to the indications reported by Mutanga et al., (2005): • strong correlation: R 2 higher than 0.7; • moderate correlation: R 2 comprises between 0.5 and 0.7; • weak correlation: R 2 less than 0.5. After investigating the correlation value, three different interpolation functions (linear, exponential, logarithmic) were implemented in order to test their ability in predicting the eastern and northern error components using the synthetic index as depend variable.

Accuracy assessment
The final step is aimed at: • assessing the accuracy of the generated photogrammetric models; • assessing the accuracy of the predictive functions; • validating the predictive functions. Photogrammetric models accuracy depends on various factors, such as the quality of the collected images, GSD and low-cost camera mounted on the UAV (Fabris, Pesci, 2005). Therefore, to meet such purposes, three main steps were implemented: assessment of acquired pictures quality, evaluation of obtained GSD and, lastly, examination of the metric reconstruction accuracy. The evaluation of acquired photos quality was carried out using the Image Quality tool implemented in Agisoft Photoscan Professional detecting an average value of 0.8, as described in the previous section. Such value is satysfying and, thus, they were adapted for photogrammetric purposes. Next, considering the non-blurry pictures only, the coherence between the programmed and the obtained GSD was performed. Conversely, the third issue involved the comparison between the GCPs, imported in the Agisoft Photoscan Professional during the alignment phase, and the coordinates assigned to those points during the reconstruction phase. This procedure was automatically and iteratively detected by Agisoft Photoscan Professional software, calculating the Root Mean Square Error (RMSE) between the observed and the estimated coordinates. RMSE provided the accuracy level expected in the photogrammetric outcome (Butler et al., 1998). Consequently, only RMSE value lower than 0.5 can be accepted. The user's ability and experience in locating GCPs on the pictures strongly affect the outcome of this step and, therefore, RMSE value can be iteratively improved repeating it until RMSE achieves satisfying value. This step was repeated for all the implemented chunks. Once the predictive functions were estimated, their performance were explored in order to assess their accuracy and to detect the optimal functions for satisfying the requirements. This phase was carried out through the calculation of RMSE between the values predicted using the proposed functions and the eastern and northern error components obtained from the photogrammetric outcomes accuracy analysis. Lastly, their predictive ability was tested using the additional dataset acquired in March since it was not involved in the modelling phase. Thus, it was applied to predict its final reconstruction accuracy just elaborating the camera calibration parameters and the RMSE was computed between them to evaluate the accuracy of the results.

RESULTS AND DISCUSSION
This paper is intended to investigate the relationship between the camera calibration parameters and the maximum metric reconstruction accuracy that could be achieved by an expert user. Several research works have demonstrated that the model applied to calibrate a camera and the obtained results of this procedure strongly affects the photogrammetric outcomes accuracy (Warner & Carson, 1991;Oniga et al., 2017). Conversely, none explored the potential accuracy that could be reached. Thus, this paper is aimed at filling this scientific gap. Therefore, five flight surveys were carried out on the experimental site of Torre a Mare, a coastline stretch in Puglia, using the UAV DJI Inspire 1, equipped with a metric-camera, DJI ZenMuse X3. Such camera was subjected to a calibration procedure through the Brown's model implemented in Agisoft Photoscan Professional. This software was also used to investigate the images quality as well as to metrically reconstruct the scene. Before aligning the photos, their quality was analysed through the index quality tool and an average value equal to 0.8 was pulled out, demonstrating their suitability for metric reconstruction purposes. Thus, they were subjected to the further photogrammetric steps in Agisoft environment and highresolution 3D models were extracted. A GSD of 4.11cm/px, 4.73cm/px, 4.82cm/px, 4.15cm/px, 4.29cm/px were obtained from the surveys carried out in December, January, February, March and October, respectively. A GSD of about 4.3cm was expected. Thus, matching the value of expected and obtained GSD, the reliability of the generated outcomes was shown. Moreover, their accuracy was evaluated as well by computing RMSE value according to the amount of GCPs implemented in the reconstruction (Figure 3). Although Figure 3 shows the optimal number of GCPs to be set to generate highly accurate photogrammetric products, it is reported to show the absolute accuracy obtained in the 31 chunks processed for each survey. The whole reconstruction shows the same trend: the optimal volume of GCPs is equal to 3 and the accuarcy values are relatively close. This means that the traits of the 3D photogrammtric model are similar. Afer investigating the goodness of the outcomes of photogrammetric steps, the calibration parameters for each chunk were extracted, organized in a data table and used as input information of the PCA technique. Each dataset was separately investigated. Such method has the invaluable property to compress the size of the original dataset extracting the most important information only (Wold et al., 1987;Bro, Smilde 2014). Figures 4-8 show the correlation plot between camera calibration parameters and principals componentes. They enhance that just three variables are needed to describe the whole original dataset taking into account the survey performed in December, January and October. On the contrary, two and five components should be considered for the survey carried out in March and February, respectively. As already discussed in the previous section, Kaiser's criterion (Kaiser,1960) parameters; k1, k2, k3, k4: radial distortions; p1, p2, p3, p4: components of the decentring distortions; Dim.1-Dim.13: principal components.  The computed synthetic index were listed in Table 1. The correlation between them and the eastern and northern error components were computed as well and also reported in Table 1.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII- B2-2020, 2020XXIV ISPRS Congress (2020 According to the indications reported by Mutanga et al., (2005), the Syntethic index shows a strong (-0.8) and a moderate (0.5) correlation with the northern and eastern error components, respectively. Nevertheless, while a positive correlation was detected with the eastern component, a negative one was identified in the other case. Thus, the former are directly related; on the contrary, the latter show an inverse dependency. Although these results demonstrate that the synthetic index can be used as a depend variable of a predictive function, we are not able to define a priori interpolating function to be used. Therefore, five functions were tested to identify the best one: linear, logarithimic, exponential, power and polinomial. Figures 9-13 show the outcomes of the five interpolation functions to predict the northern error component. The coeffiecient of determination (R 2 ) was calculated for each functions in order to detect those ones showing the best performance. Although all the implemented functions show a good value of R 2 , higher than 0.51, the polinomial equation seams the best one since the coefficient of determination was equal to 0.78. Figure 9. Predictive function extracted linearly interpolating the depend (synthetic index) and independent variables (northern error component). R 2 : coefficient of determination. Figure 10. Predictive function extracted through a logarithmic interpolation between depend (synthetic index) and independent variables (northern error component). R 2 : coefficient of determination. Figure 11. Predictive function extracted through an exponential interpolation between depend (synthetic index) and independent variables (northern error component). R 2 : coefficient of determination. Figure 12. Predictive function extracted through a polynomial interpolation between depend (synthetic index) and independent variables (northern error component). R 2 : coefficient of determination. Figure 13. Predictive function extracted through a power interpolation between depend (synthetic index) and independent variables (northern error component). R 2 : coefficient of determination.
Conversely, the outcomes produced considering the eastern error components (Figures 14-18) has a different trend, indeed, just the polynomial functions appear reliable, showing a high coefficient of determination (0.68). Figure 14. Predictive function extracted through a polynomial interpolation between depend (synthetic index) and independent variables (eastern error component). R 2 : coefficient of determination. Figure 15. Predictive function extracted through a logarithmic interpolation between depend (synthetic index) and independent variables (eastern error component). R 2 : coefficient of determination.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) Figure 16. Predictive function extracted through an exponential interpolation between depend (synthetic index) and independent variables (eastern error component). R 2 : coefficient of determination. Figure 17. Predictive function extracted through a linear interpolation between depend (synthetic index) and independent variables (eastern error component). R 2 : coefficient of determination. Figure 18. Predictive function extracted through a power interpolation between depend (synthetic index) and independent variables (eastern error component). R 2 : coefficient of determination.
Moreover, to test the accuracy of the implemented functions, the RMSE between the values predicted by the introduced equations and the eastern and northern error components obtained from the dataset acquired in March were computed, since it was not included in their definition. For brevity, the average of RMSE related to polynomial equations are reported since they were recognized as the best one: 0.005 and 0.003 for northern and eastern error components, respectively.

CONCLUSION
The relationship between camera calibration parameters and the accuracy of the photogrammetric outcomes has not been deeply examined yet in literature. In fact, although several researchers have demonstrated the importance of calibration steps and the strong impact of applying different methods on the results, there are no studies quantifying the potential accuracy that should be achieved by an expert user. Therefore, this paper is intended to make a first step in that direction, proposing a methodology to investigate such potential and to define predictive functions of error components. The resulting equation was computed through the implementation of multivariate and linear statistics techniques. Their combination shows promising results for predicting 3D models accuracy; indeed, the northern and the eastern error components were detected with a deviation value equal to 0.005 m and 0.003 m, respectively. Further works are required to test the methodology on the other error component too and on different landscape scenarios as well.