Combining Interior Orientation Variables to Predict the Accuracy of Rpas–Sfm 3D Models

: Remotely piloted aerial systems (RPAS) have been recognized as an e ﬀ ective low-cost tool to acquire photogrammetric data of low accessible areas reducing collection and processing time. Data processing techniques like structure from motion (SfM) and multiview stereo (MVS) techniques, can nowadays provide detailed 3D models with an accuracy comparable to the one generated by other conventional approaches. Accuracy of RPAS-based measures is strongly dependent on the type of adopted sensors. Nevertheless, up to now, no investigation was done about relationships between camera calibration parameters and ﬁnal accuracy of measures. In this work, authors tried to ﬁll this gap by exploring those dependencies with the aim of proposing a prediction function able to quantify the potential ﬁnal error in respect of camera parameters. Predictive functions were estimated by combining multivariate and linear statistical techniques. Four photogrammetric RPAS acquisitions were considered, supported by ground surveys, to calibrate the predictive model while a further acquisition was used to test and validate it. Results are preliminary, but promising. The calibrated predictive functions relating camera internal orientation (I.O.) parameters with ﬁnal accuracy of measures (root mean squared error) showed high reliability and accuracy.


Introduction
Descriptions of the earth's morphology and environmental processes involve the adoption of a large amount of data from a wide range of branches of knowledge, such as geology [1], hydrology [2] and geomorphology [3]. All of them require a digital surface model (DSM) as the basis to describe surface morphometry. DSM resolution and accuracy strongly affect landforms mapping: in fact, geomorphic elements, located at the same position, can be differently interpreted depending on DSM resolution [4]. Consequently, DSM resolution must be consistent with the size of the investigated element to optimize results [5]. Along the years, several approaches were proposed to improve 3D models accuracy and detail level of survey. Among these, terrestrial laser scanner (TLS) was widely considered as the gold standard to meet such requirements [6], even if its application determines high costs of both technological equipment and acquisition/processing data time. Researchers looked for a valid low-cost alternative to overcome these operational limits. Rieke-Zapp et al. [7] demonstrated that photogrammetry is a convincing tool to generate 3D models having comparable accuracy and resolution in respect of those from TLS. Figure 1. Study area. Colored orthophoto was generated during this work from one of the available datasets. In red, ground control points (GCPs) locations. They were surveyed by network real time kinematic (NRTK) Global Navigation Satellite System (GNSS) and georeferenced in the RDN2008/UTM Zone 33N (NE) reference frame (EPSG: 6708). In the background (grayscale) an orthophoto dated 2016 is shown as provided by the WMS Service of SIT Puglia.
Subhorizontal surfaces, close to the reef, are mainly composed of calcarenitic lithotypes showing a peculiar resistance to mechanical erosion processes. At the global scale, stratified layers and fracturing systems, subparallel to the coastline (ONO ESE), determine favorable conditions for land instability due to hydric erosion, accelerating natural withdrawal of coast fractions ( Figure 2). Moreover, local lack of vegetation increases imperviousness and, consequently, facilitate erosion dynamics [32]. The area is extremely complex from the archaeological perspective, as well. It hosts archeological evidence of a Neolithic community that was pulled out by excavation campaigns carried out in the past; currently, the area is in an evident state of abandon ( [37,38]).  Subhorizontal surfaces, close to the reef, are mainly composed of calcarenitic lithotypes showing a peculiar resistance to mechanical erosion processes. At the global scale, stratified layers and fracturing systems, subparallel to the coastline (ONO ESE), determine favorable conditions for land instability due to hydric erosion, accelerating natural withdrawal of coast fractions ( Figure 2). Moreover, local lack of vegetation increases imperviousness and, consequently, facilitate erosion dynamics [32]. The area is extremely complex from the archaeological perspective, as well. It hosts archeological evidence of a Neolithic community that was pulled out by excavation campaigns carried out in the past; currently, the area is in an evident state of abandon ( [37,38]).
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 31 Figure 1. Study area. Colored orthophoto was generated during this work from one of the available datasets. In red, ground control points (GCPs) locations. They were surveyed by network real time kinematic (NRTK) Global Navigation Satellite System (GNSS) and georeferenced in the RDN2008/UTM Zone 33N (NE) reference frame (EPSG: 6708). In the background (grayscale) an orthophoto dated 2016 is shown as provided by the WMS Service of SIT Puglia.
Subhorizontal surfaces, close to the reef, are mainly composed of calcarenitic lithotypes showing a peculiar resistance to mechanical erosion processes. At the global scale, stratified layers and fracturing systems, subparallel to the coastline (ONO ESE), determine favorable conditions for land instability due to hydric erosion, accelerating natural withdrawal of coast fractions ( Figure 2). Moreover, local lack of vegetation increases imperviousness and, consequently, facilitate erosion dynamics [32]. The area is extremely complex from the archaeological perspective, as well. It hosts archeological evidence of a Neolithic community that was pulled out by excavation campaigns carried out in the past; currently, the area is in an evident state of abandon ( [37,38]).   Flight were operated by a commercial quadcopter DJI Inspire 1, mounting a consumer DJI Zenmuse X3 camera (focal length 3.61 mm, pixel size 1.56 µm, effective pixels 12.4 M) equipped with a 3−axis gimbal (to compensate accidental movements of the drone). The DJI Ground Station Pro app, proposed by the Chinese company DJI (Dà−Jiāng Innovations) [41] was also used to automatize the flight plan. It ensured that all the flights were achieved along the same path and under the same conditions, e.g., cruising speed (4.0 m/s) and altitude (100 m AGL, above ground level) ( Figure 3). Missions were planned to obtain an average Ground Sampling Distance (GSD) of 0.043 m/pix, a forward and side overlaps of 85% and 75%, respectively, as suggested by [42,43]. A total of 77 images per flight were acquired. Camera was set nadiral and the stop&go mode was applied to reduce collection of blurry images due to forward motion [42]. RPAS position was recorded by a low cost GNSS/INS positioning receiver and saved into a metadata file.

Field Data Campaigns and Operative Workflow
An extensive flight campaign was planned between 2018 and 2019. Flight area was selected far away from urban sites, to compliant Italian national regulations about RPAS operations, [39] and [40]. Five flights were programmed and performed in December 2018, January 2019, February 2019, March 2019 and October 2019, respectively (Table 1). Surveys were interrupted between April and September 2019 to respect operational restrictions related to touristic season [39].  [41] was also used to automatize the flight plan. It ensured that all the flights were achieved along the same path and under the same conditions, e.g., cruising speed (4.0 m/s) and altitude (100 m AGL , above ground level) ( Figure 3). Missions were planned to obtain an average Ground Sampling Distance (GSD) of 0.043 m/pix, a forward and side overlaps of 85% and 75%, respectively, as suggested by [42] and [43]. A total of 77 images per flight were acquired. Camera was set nadiral and the stop&go mode was applied to reduce collection of blurry images due to forward motion [42]. RPAS position was recorded by a low cost GNSS/INS positioning receiver and saved into a metadata file. To operate all along the tests with the same ground control points, a GNSS survey campaign was operated to position permanent natural elements having adequate size and color with respect to RPAS image features. Thirty homogeneously distributed points were consequently, surveyed by Leica Viva CS10/GS10 GNSS receiver and successively used as ground control points (GCP) or check points (CPs). Survey was operated in a network real time kinematic (NRTK) mode based on Leica To operate all along the tests with the same ground control points, a GNSS survey campaign was operated to position permanent natural elements having adequate size and color with respect to RPAS image features. Thirty homogeneously distributed points were consequently, surveyed by Leica Viva Remote Sens. 2020, 12, 2674 5 of 31 CS10/GS10 GNSS receiver and successively used as ground control points (GCP) or check points (CPs). Survey was operated in a network real time kinematic (NRTK) mode based on Leica SmartNet Italpos network and determined a final 3D accuracy of 0.02 m. Reference system was RDN2008/UTM zone 33N (NE) (EPSG: 6708).
Image datasets were separately processed according to the workflow of Figure 4. Image blocks of December, January, February and October, were used to calibrate the predictive model (see forward on), while the March dataset was used for validation.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 31 SmartNet Italpos network and determined a final 3D accuracy of 0.02 m. Reference system was RDN2008/UTM zone 33N (NE) (EPSG: 6708). Image datasets were separately processed according to the workflow of Figure 4. Image blocks of December, January, February and October, were used to calibrate the predictive model (see forward on), while the March dataset was used for validation.

Photogrammetric Products Generation
This study relies on the workflow proposed by [11,[44][45][46]. The work flowchart is shown in Figure 5; steps are detailed in the next sections. Agisoft PhotoScan (v.1.4.1, Agisoft LLC −St. Petersburg, Russia) software, currently known as Metashape, was used during the work to photogrammetrically process the data. This first step was aimed at properly setting workspace and removing blurry images, possibly compromising final outcomes. Five chunks (image blocks) of the same scenario were created and, for each of them, the same processing parameters were set to guarantee results comparability (Table 2).

Photogrammetric Products Generation
This study relies on the workflow proposed by [11,[44][45][46]. The work flowchart is shown in Figure 5; steps are detailed in the next sections. Agisoft PhotoScan (v. 1 SmartNet Italpos network and determined a final 3D accuracy of 0.02 m. Reference system was RDN2008/UTM zone 33N (NE) (EPSG: 6708). Image datasets were separately processed according to the workflow of Figure 4. Image blocks of December, January, February and October, were used to calibrate the predictive model (see forward on), while the March dataset was used for validation.

Photogrammetric Products Generation
This study relies on the workflow proposed by [11,[44][45][46]. The work flowchart is shown in Figure 5; steps are detailed in the next sections. Agisoft PhotoScan (v.1.4.1, Agisoft LLC −St. Petersburg, Russia) software, currently known as Metashape, was used during the work to photogrammetrically process the data. This first step was aimed at properly setting workspace and removing blurry images, possibly compromising final outcomes. Five chunks (image blocks) of the same scenario were created and, for each of them, the same processing parameters were set to guarantee results comparability ( Table 2).  This first step was aimed at properly setting workspace and removing blurry images, possibly compromising final outcomes. Five chunks (image blocks) of the same scenario were created and, for each of them, the same processing parameters were set to guarantee results comparability ( Table 2). Camera positioning accuracy was set equal to 3 m to make it consistent with the average 3D positioning accuracy value of the RPAS GNSS receiver (approximately 2.54 m). This value is known to depend on GNSS epoch recording frequency (Hz) and, consequently, on RPAS speed [46]. Image attitude accuracy was set equal to the software default value (10 degrees) since no information was available about attitude measurements from RPAS-integrated IMU (inertial measurement unit). This value was considered precautionary. Agisoft PhotoScan allows weights parameters during bundle block adjustment (BBA): tie points estimation are generally three times less relevant than GCPs accuracy during the reconstruction phase [17]. Thus, a value equal to 3 was set in its parameterization. All of these parameters are crucial being directly involved in collinearity equations [47][48][49] and, consequently, they heavily affect final accuracy of solution. Software default values were instead accepted as weights for attitude parameters (yaw, tilt, roll).
A quality assessment was performed to detect and remove "bad" images. This step is mandatory to improve final accuracy. A large number of issues conditioning image quality depends on the acquisition mode. For instance, the adoption of RPAS to image a complex scenario, like the study area, involves several problems due to the interaction between environmental conditions and equipment; in particular, heat and magnetic sources can impact on inertial measurement unit (IMU) and GNSS [8]. The "estimate image quality" procedure available within Agisoft PhotoScan was used to assess image characteristics, providing information about sharpness and detecting blurring and distortion. This procedure returns a score ranging between 0 and 1: the higher the value, the better the quality [50,51]. Moreover, to balance colors of final products, all images were homogenized in terms of brightness and contrast. Radiometric adjustments are known to do not condition the efficiency of SIFT algorithm during tie point detection.
A scale-invariant feature transform (SIFT) approach was used to minimize projection errors [52,53] and facilitate extraction of homologous points independently from brightness conditions, that, in the datasets, were strongly influenced by the presence of sea.

Second
Step: Image Block Orientation The II step was aimed at automatically collecting tie points (sparse points cloud) and solving image orientation [54,55]. It was separately run for each chunk using the 'High' accuracy mode and setting a threshold = 0 for the "limits of key points and tie points" parameter. This choice was preferred to avoid uncontrolled filtering of measured points.
Block Bundle Adjustment (BBA) was performed including I.O. parameters within the model unknowns to be estimated using the camera optimization panel available in Agisoft PhotoScan. BBA outputs correspond to the estimates of tie-points coordinates in the object space, I.O. and external orientation (E.O.) parameters, GNSS lever arm offset. Moreover, BBA scales the entire photogrammetric block proportionally with respect to GCPs. Camera I.O. parameters estimates from all the chunks were recorded and organized into a table to be successively analyzed. Sparse point clouds (tie-point positions in the object space) were generated and systematic errors, mainly caused by nonlinear distortions of lens [54], estimated; measured points were manually filtered to optimize results and minimize image block distortions. Three criteria implemented in the "gradual selection" tool were considered: (a) photogrammetric restitution uncertainty; (b) projection accuracy; (c) reprojection error.
(a) is intended to remove points with low base-height ratios [49], i.e., all those points located at the edges of images, generally characterized by a higher degree of restitution uncertainty, that majorly depends on a too small overlapping among pictures. Although this option does not affect final accuracy, it is useful for thinning clouds [56]. Conversely, (b) is aimed at detecting and cleaning out less reliable tie points [49]. A threshold value equal to 3 was used to exclude tie points with an uncertainty 3 times higher than the minimum one; (c) was intended to remove all points with a large residual value in order to decrease drastically restitution errors improving orientation parameters estimates [49]. It is worth to remind that, residuals directly impact on representativeness of the root mean square error (RMSE, [40]) of GCPs and CPs, possibly making it not suitable to define the actual final accuracy of measurements [51].
The three above-mentioned criteria permitted to remove the most of inaccurate points from clouds, thus improving consistency between model and reality. Threshold values adopted for each criterion were defined accordingly to previous works [42,43]. After filtering, about 20% of originally measured points were removed. To further improve geometric modeling of the area, 30 GCPs were used [5], ensuring a marker reprojection error less than 0.5 pixels. Once the alignment step was accomplished, a "progressive" cross-validation analysis was achieved to test actual accuracy of final measurements, as explained in Section 2.3.4.

Fourth Step: Progressive Cross-Validation
With respect to the 30 surveyed points, a cross-validation was run in a progressive mode, i.e., progressively migrating 1 point a time from the training (GCPs) to the validation (CPs) set. The choice of those points to be progressively moved from GCPs to CPs was accomplished in a balanced way; peripheral and central points were alternatively migrated taking care of maintaining a proper spatial distribution of remaining GCPs ( [22,57]). Thirty one chunks of CPs were finally obtained varying from 0 (all surveyed points were used as GCPs) to 30 points (all surveyed points were used as CPs and the image orientation completely relied on a direct georeferencing approach) [50,58]. We refer to the 0 CPs situation as complete indirect georeferencing (CIG) and to the 30 CPs situation as direct georeferencing (DG) [21,50]. To make the procedure repeatable, points were imported in the same order for all analyzed subdatasets. Table 3 summarizes the order that was followed while migrating a point from GCP to CP. Measures accuracy was estimated with reference to RMSE computed for both GCPs and CPs. RMSE was used to quantify error components (definitions are given forward on), with the hypothesis that both random and systematic errors were Gaussian distributed [59]. It is worth to remind that RMSE of GCPs only provide information about the goodness of fitting of calibrated equations, with no concern about the model capability of generalization.

Analyses of I.O. Parameters Estimates
Camera calibration by SfM implies that camera I.O. parameters (including lens distortion) are precisely known to minimize errors possibly affecting restitution step [60]. Such parameters are image independent and, consequently, they do not depend on position and attitude of the camera [54]. Several approaches have been proposed over the years to get appropriate estimates/measures of these parameters; camera self-calibration guarantees several benefits when working with low-cost cameras. It relies on the numeric estimate of I.O. parameters that are included among the unknowns in the equation system that is solved during BBA. The solution is consistent with the data [61], with the assumption that I.O. parameters continuously vary in time and, consequently, need to be estimated for each specific image block.
where f is the focal length, ∆x and ∆y represents the image corrections, ∆ f is the correction to the initial principal distance value, x and y are the coordinates of a general poin, K i are the radial distortion coefficients, P i are the tangential distortion coefficients, B i are the in-plane correction parameters for differential scaling between the horizontal and vertical pixel spacing and non-orthogonality (axial skew) between x and y axes, r is the image radial distance estimated using Equation (3): where x p and y p represents the principal point coordinates.
The full 10-parameter model, Equations (1)- (3), is adopted as the default one when using the fully automatic camera calibration procedure. Although, practically, it could appear as the best choice, Remote Sens. 2020, 12, 2674 9 of 31 from an accuracy point of view, in many cases, it does not represent the optimal solution, because of the different role and significance of those parameters within the BBA numeric solution [55]. For example, a high correlation was found between P i coefficients and principal point coordinates. Consequently, when removing P i from the unknowns, x p and y p somehow can absorb associated variation. In other words, users could not consider P i in the calibration phase and obtain similar results. With these premises, a preventive analysis was achieved to figure out eventual correlation existing among I.O. parameters and trying to minimize it.

Accuracy Assessment
Accuracy assessment was aimed at: (a) testing the quality of the original dataset; (b) testing the accuracy of the photogrammetric measurements.
Triggs et al. [54] highlighted that accuracy of photogrammetric products depends on many factors: quality of the processed images, GSD and camera type. Consequently, all these factors were considered. The first factor was taken into account (trying to minimize its effects) using the image quality tool and selecting a scale-invariant feature transform (SIFT) approach while running BBA by SfM (see Section 2.3.1).
Errors affecting final measures from oriented image blocks were evaluated with reference to the RMSE. It was computed separately for all the error components and for both GCPs and CPs. RMSE E for the East coordinate, RMSE N for the North coordinate, RMSE H for the height coordinate, RMSE T as the total 3D error and, finally, RMSE I for the positioning error (total) in the image space.
Agisoft PhotoScan Professional software can automatically and iteratively compute this parameter [48]. This software feature is important, permitting reiteration of trials that can be run selectively tuning all the involved parameters, including GCPs collection refinement.

PCA and Synthetic Index Generation
I.O. estimates and errors computed during BBA were compared to test their reciprocal relationship. I.O. estimates were preventively preprocessed by a self-developed R routine [66] aimed at extracting the most significant information by principal component analysis (PCA), based on the variance maximization principle [67]. PCA, probably the most popular multivariate statistical method, is intended for dimensionality reduction, obtained by removing redundant information of a multivariate dataset where variables can be intercorrelated [67]. After detecting the most relevant components, it converts the original dataset into a new one consisting of independent and orthogonal vectors, called principal components (PCs, [67]). The first component provides the most of information, describing the largest part of the input data inertia, absorbing (explaining) the most of data variance; the second component is orthogonal to the first one and absorbs (explain) the most of the remaining variance [68]. The same principle is used to obtain all the other components that, necessarily, will represent a decreasing level of information while incrementing their position within the transformation [69]. Consequently, first components compress most of the native information making possible to drastically reduce dimensionality of data by removing redundant content [62].
The singular value decomposition (SVD) approach [69] was applied to compute principal components. The number of explanatory PCs was selected by the Kaiser's criterion [70] that suggest setting an eigenvalue threshold = 1.0. Selected PCs were weighted and linearly combined in a "synthetic index" (hereinafter called SI). SI was obtained as the weighted average of all significant components (Equation (4)). Weights were directly extracted from the PCA procedure: where w i and Dim i are the weights and the principal components, respectively.

Predicting Accuracy of Measurements: Model Definition
Pearson's coefficient (R) was computed between SI and error components affecting final measurements from oriented blocks. The expectation was that SI could be a predictor of error components (generically called RMSE j ). According to the obtained R values (only high or moderate correlation was considered) an interpolation function was calibrated to predict the following errors: RMSE E , RMSE N , RMSE H , RMSE T , RMSE I . SI was assumed as independent variable (x, predictor) of calibrated functions. A 2nd order polynomial (Equation (5)) was found to well fit all significant relationships. Model parameters (a, b and c) were estimated by ordinary least squares for each investigated error (y). Goodness of fitting was tested with reference to the coefficient of determination (R 2 ):

Accuracy of Photogrammetric Measurements
Suitability of available images to generate a reliable photogrammetric product was assessed through the application of the image quality tool in Agisoft PhotoScan. All images showed a satisfying value presenting an average image quality index equal to 0.8 (higher than the threshold value set to 0.5); consequently, no image was removed.
Accuracy of photogrammetric measurements were assessed with reference to RMSE T for both GCPs and CPs. They were tested against the number of GCPs used during BBA (Figures 6 and 7). Highest RMSE values (for both GCPs and CPs) correspond to the solution obtained using 1 or 2 GCPs. This confirms that, minimally 3 points are needed to properly orient a 3D model in a 3D space to recover translation, rotation and scaling values [8]. Nevertheless, the SfM/SIFT based approach can provide reasonable solutions of BBA even when GCPs number is lower than the theoretically minimum one, i.e., 3, just exploiting position and attitude data from RPAS GNSS/IMU low quality system. In spite of this interesting consideration characterizing the new digital photogrammetry, what is evident is that accuracy and processing time increases with GCPs number.    In this case, the CIG method is not reported. As previously, highest RMSE values were obtained when GCPs number was lower than 3, demonstrating that GNSS and IMU measures from RPAS system, in spite of their low quality, can drive to a reasonable solution (0.5-1 m) that could be accepted for many applications. This appeared to be similar for both CGPs and CPs, making those RMSE values sufficiently representative of the actual accuracy condition when working in DG mode. For the goals of this work, this fact is secondary and certainly must be overcame focusing the attention on the best expected solutions. In general, one can say that a good solution for BBA is the one when RMSE GCP and RMSE CP are consistent each other. With respect to Figures 6 and 7 this situation occurs reaching a number of well distributed GCPs around 14. More GCPs would appear as useless.  In this case, the CIG method is not reported. As previously, highest RMSE values were obtained when GCPs number was lower than 3, demonstrating that GNSS and IMU measures from RPAS system, in spite of their low quality, can drive to a reasonable solution (0.5−1 m) that could be accepted for many applications. This appeared to be similar for both CGPs and CPs, making those RMSE values sufficiently representative of the actual accuracy condition when working in DG mode. For the goals of this work, this fact is secondary and certainly must be overcame focusing the attention on the best expected solutions. In general, one can say that a good solution for BBA is the one when RMSEGCP and RMSECP are consistent each other. With respect to Figure 6 and 7 this situation occurs reaching a number of well distributed GCPs around 14. More GCPs would appear as useless.  Table 4 reports the main statistics computed for all the estimated I.O. parameters with respect to all processed dataset. It can be noted that I.O. parameters can be grouped in two clusters: one including x p , y p , B 1 , B 2 , P 3 and P 4 is characterized by a great variability; another one, including remaining parameters, shows pretty similar values for all the computed metrics. Since all trials generated similar estimates of radial and tangential distortion parameters, these can be retained strictly dependent on the camera with no conditioning by other factors. Moreover, P i coefficients are known to be less significant than radial ones (one or two orders of magnitude smaller [59]). Boxplots of Figures 9-13 confirm this fact. B 2 -skew coefficients; K 1 , K 2 , K 3 , K 4 -radial distortion coefficients; P 1 , P 2 , P 3 , P 4 -decentering distortion coefficients.

Survey
Stat A further synthetic analysis can come from the coefficient of variation (CV = standard deviation/mean × 100.0, %) of I.O. parameters estimates. Results are reported in Figure 8. It can be noted that focal length f and K i coefficients are the most stables parameters. Conversely, x p , y p , B i and P i coefficients present a wide variability depending on operational conditions. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing A further synthetic analysis can come from the coefficient of variation (CV = standard deviation/mean x 100.0, %) of I.O. parameters estimates. Results are reported in Figure 8. It can be noted that focal length f and Ki coefficients are the most stables parameters. Conversely, xp, yp, Bi and Pi coefficients present a wide variability depending on operational conditions.  According to the boxplots of Figures 9-13 a similar trend can be recognized affecting all I.O. parameters. Specifically, the October dataset ( Figure 13) appears to not contain any outlier for the most of parameters (f, K 1 -K 2 -K 3 -K 4 -P 3 -P 4 ). Oppositely, in the other datasets no outlier was found for B 1 , P 1 , P 2 in the December dataset ( Figure 9); y p and x p did not present outliers in the January and February datasets (Figures 10 and 11); f, P 4 and K 2 had no outliers in the March dataset ( Figure 12). Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 31      . Positive values mean that a direct relationship between variables is present; negative values an inverse linkage. The higher the values, the stronger the correlation. All the datasets showed pretty similar values, that can be synthesized as it follows: focal length was highly correlated with radial aberrations coefficients; x p was moderately correlated with y p and, in the most of cases, with P i coefficients; y p was poorly/moderately correlated with all the other parameters (March dataset excluded); skew parameters (B 1 ,B 2 ) were poorly/moderately correlated with the other parameters (March dataset excluded); K i coefficients showed to be internally correlated each other and, externally, with the focal length; a moderate correlation was found between y p and P 2 . January and February dataset showed low correlation between K i with y p and P 2 . In spite of all specific situations, these results suggested that I.O. parameters contain redundant information. Consequently, a PCA was operated showing that the most of the decorrelated information could be explained by few PCs, whose identification was operated by Kaiser's criterion.  Table 9. Correlation matrix of camera I.O. parameters as estimated from the October 2019 dataset (f-focal length; x p and y p coordinates of the principal point offset; B 1 , B 2 -skew parameters; K 1 , K 2 , K 3 , K 4 -radial distortions; P 1 , P 2 , P 3 , P 4 -components of the decentering distortions). The errors affecting final measures from oriented image blocks were evaluated with reference to the above mentioned RMSE j . Analysis was separately performed for both GCPs and CPs and for all the datasets (Table 10). It showed similar values for all the datasets. RMSE I defines the reprojection error computed as the mean distance (in the image space) between the position expected for a tie point that participated to solve block orientation and the one resulting by reprojection of the correspondent 3D object point after image resection. To minimize alignment issues, the maximum values of error should be <1, according to [44]. This condition was respected for all the blocks: the maximum detected value was 0.31 and 0.36 for GCPs an CPs (Table 11), respectively. Conversely, RMSE N , RMSE H , RMSE T define the distance (in the object space) between the expected position of a GCP (or CP) and the one determined by photogrammetric measurement in the object space. The smaller this value, the higher the final accuracy. These parameters are known to be strongly influenced by user's experience in recognizing the proper markers (shape, color, stability, etc.) within image and operate the correspondent collimation. The highest RMSE T values were detected in the December (0.98 m) and March (1.01 m) datasets.
The correlation matrices were computed to analyze possible relationships among error components for both GCPs and CPs (Figures 10 and 11). Results showed a high direct correlation among all the variables, reprojection error excluded (RMSE I ).

Calibrating Predictive Models
Camera I.O. parameters estimates were analyzed by PCA to detect and remove redundant information. As usual, December, January, February and October blocks were separately processed to identify those PCs that, for each dataset, could synthesize the most of the original information. This was done with reference to the correlation plots relating I.O. parameters with PC ( Figure 14). Application of Kaiser's criterion showed that the strongest three PCs were enough to describe the most of I.O. parameters variance for the December (Figure 14a), January ( Figure 14b) and October (Figure 14d) datasets. Conversely, five PCs were needed to explain the most of information resident in I.O parameters estimates from the February dataset (Figure 14c)

Calibrating Predictive Models
Camera I.O. parameters estimates were analyzed by PCA to detect and remove redundant information. As usual, December, January, February and October blocks were separately processed to identify those PCs that, for each dataset, could synthesize the most of the original information. This was done with reference to the correlation plots relating I.O. parameters with PC ( Figure 14). Application of Kaiser's criterion showed that the strongest three PCs were enough to describe the most of I.O. parameters variance for the December (Figure 14a  The most significant PCs, as resulting from the previous analysis, were aggregated by eq. 4 to compute SI that was used as predictor within the predictive functions of RMSE. SI values are reported in Table 12. Before calibrating predictive functions possibly relating SI to RMSE the Pearson's coefficient was computed between SI and all the available RMSE for both GCPs and CPs (Table 13). The most significant PCs, as resulting from the previous analysis, were aggregated by Equation (4) to compute SI that was used as predictor within the predictive functions of RMSE. SI values are reported in Table 12. Before calibrating predictive functions possibly relating SI to RMSE the Pearson's coefficient was computed between SI and all the available RMSE for both GCPs and CPs (Table 13). Only situations showing a moderate (0.5 < R < 0.7) or a high Pearson's coefficient (>0.7) were modeled in the following step [71]. Consequently, since SI showed a moderate and high correlation with GCPs RMSE E and RMSE N , respectively these relationships were modeled. Conversely, since CPs RMSE E showed a low correlation with SI it was excluded from modeling.
A 2nd order polynomial showed to better fit the data. Results concerning GCPs and CPs are shown in Figures 15 and 16, respectively. Graphs also report the coefficient of determination (R 2 ). For GCPs, the lowest R 2 value was 0.68 (RMSE E ), while CPs RMSE N showed the highest R 2 (0.99).
Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 31 Only situations showing a moderate (0.5 < R < 0.7) or a high Pearson's coefficient (>0.7) were modeled in the following step [71]. Consequently, since SI showed a moderate and high correlation with GCPs RMSEE and RMSEN, respectively these relationships were modeled. Conversely, since CPs RMSEE showed a low correlation with SI it was excluded from modeling.  Figures 15 and 16, respectively. Graphs also report the coefficient of determination (R 2 ). For GCPs, the lowest R 2 value was 0.68 (RMSEE), while CPs RMSEN showed the highest R 2 (0.99).  Table 14 reports the coefficients of the significant 2nd order polynomial predictive functions for both CPs and GCPs.   Table 14 reports the coefficients of the significant 2nd order polynomial predictive functions for both CPs and GCPs.

Validating Predictive Models
Reliability and accuracy of the proposed predictive method were tested by applying the calibrated models to all the available datasets, March included. With respect to March dataset specifically, PCA analysis was applied to recognize significant PCs able to explain the most of I.O. parameters variance. Two PCs were found to satisfy the Kaiser's criterion and, consequently, they were used to compute March SI (3.060). March SI value was then used to predict RMSEj according to eq. 5 applied using the coefficients from Table 14. To summarize performances of models, all the RMSEj estimates were compared, by differencing, with the correspondent values from the BBA solutions (reference ones). RMSE was then computed for all the tested differences and estimated RMSEj. Results are reported in Table 15.

Discussion
This research was intended to explore dependency of the accuracy of measures from RPASbased SfM 3D models from the camera I.O. parameters estimates. Several researchers investigated these issue, detecting the strong influence of camera parameters on final accuracy of obtainable measures. Nevertheless, no proposal came for an operative procedure able to predict the potential achievable accuracy once camera parameters were known. In this work, authors tried to fill this gap by developing and proposing a simple method that, in their preliminary tests, provided promising

Validating Predictive Models
Reliability and accuracy of the proposed predictive method were tested by applying the calibrated models to all the available datasets, March included. With respect to March dataset specifically, PCA analysis was applied to recognize significant PCs able to explain the most of I.O. parameters variance. Two PCs were found to satisfy the Kaiser's criterion and, consequently, they were used to compute March SI (3.060). March SI value was then used to predict RMSEj according to Equation (5) applied using the coefficients from Table 14. To summarize performances of models, all the RMSEj estimates were compared, by differencing, with the correspondent values from the BBA solutions (reference ones). RMSE was then computed for all the tested differences and estimated RMSEj. Results are reported in Table 15.

Discussion
This research was intended to explore dependency of the accuracy of measures from RPAS-based SfM 3D models from the camera I.O. parameters estimates. Several researchers investigated these issue, detecting the strong influence of camera parameters on final accuracy of obtainable measures. Nevertheless, no proposal came for an operative procedure able to predict the potential achievable accuracy once camera parameters were known. In this work, authors tried to fill this gap by developing and proposing a simple method that, in their preliminary tests, provided promising results. The proposed approach integrates uni-and multivariate statistics to investigate and removing correlated information resident in I.O. parameters of the camera as estimated during BBA.
The study was based on five photogrammetric surveys by RPAS operated in December 2018 and January, February, March and October of 2019 (Table 1). A ground survey campaign was also done to position 30 well distributed GCPs. All flight missions were performed by DJI Inspire 1 drone, equipped with DJI ZenMuse X3 camera. To ensure comparability of datasets all missions were operated according to the same flight plan (path, speed and height) and BBA was performed setting the same parameters within Agisoft PhotoScan. The same skilled user was in charge of processing data and optimizing the photogrammetric solution.
A first step was aimed at testing image blocks quality and removing "bad" images. A second step was aimed at solving BBA iteratively changing the number of GCPs. A total of 155 chunks (31 chunks for each mission) were analyzed. It is worth to remind that this step was not aimed at selecting the best spatial strategy for positioning GCPs; differently, it was devoted to evaluating the achievable accuracy of final measures and the accuracy dependence on GCPs number. As shown in Figures 6 and 7, highest RMSE T correspond to those solutions were direct georeferencing plays the leading role, being the number of GCPs lower than the theoretical minimum value (3, [8]). By comparing GCPs and CPs RMSE T trend with the number of GCPs it was found that the optimal minimum number of GCPs was around 14 (Figures 6 and 7). All the datasets showed the same GCPs and CPs RMSE T trend. Just few dissimilarities were found, probably due to different environmental conditions (e.g., lighting, weather conditions, vegetation phenology). After this preliminary investigation that ensured about comparability of processed datasets, a deeper investigation concerned errors components separately (RMSE j ). Some basic statistics (e.g., maximum, minimum, mean and standard deviation) of I.O. parameters estimates were also computed for each processed dataset ( Table 10). All of them showed similar statistics. Results confirmed what reported in previous works: Agisoft PhotoScan, as well as other photogrammetric software, in general, cannot estimate stable I.O. parameters while initial conditions (e.g., GCPs number) change [29,59]. In fact, stats showed a high variability of solutions. Nevertheless, the order of magnitude remained the same in all dataset (Table 4) and a similar trend, as shown in the boxplots of Figures 9-13. Moreover, the obtained order of magnitude is coherent with the one obtained by other researchers [59]; they showed that P i parameters are smaller by one or two orders of magnitude than radial ones and that K i present the most significant deviations.
Correlations among I.O. parameters was then investigated finding a high degree of intracorrelation. Correlation values proved to be consistent with those reported in literature (Tables 5-9, [59]). Correlated information was aggregated by PCA [60]; it was found that from two to five PCs are generally enough to explain the most of variance of I.O. parameters estimates.
With these premises, an index (SI) was defined to summarize the decorrelated information that the first PCs were able to aggregate. SI was assumed as predictor of RMSE j and correspondent predictive function calibrated (Figures 15 and 16). Models calibration was obtained with reference to the December, January, February and October datasets. March dataset was differently used to validate predictive models.
In spite of the few observations used for calibration, the proposed predictive functions showed satisfying results when applied to the validation set generating RMSE j estimates very close to the actual values as computed during BBA by Agisoft PhotoScan. Results are must be assumed as preliminary, but certainly encouraging.

Conclusions
The present study was aimed at exploring the impact of the camera I.O. parameters on the accuracy of final photogrammetric 3D models. The possibility of calibrating predictive models for error estimates once given I.O. parameters estimates was the focus point of this research. The proposed procedure was tested on the study area of Torre a Mare (Apulian Region) through the acquisition of five photogrammetric datasets operated at different dates along the year.
After investigating image quality and achieved BBA, resulting camera I.O. parameters were explored to test their eventual intracorrelation and potential relationships with the accuracy of final photogrammetric measures. A high correlation among the most of parameters was found and a high level of information reduction can be achieved by PCA. Tests proved that from two up to five PCs are enough to explain the most of I.O. parameters variance. PCs selection was operated according to the Kaiser's criterion. The strongest PCs-if properly aggregated along SI-showed that they are able to predict final errors in photogrammetric measures. The SI was, in fact, found to be a good predictor of errors when included, as independent variable, within a 2nd order polynomial function. Predictive functions dependent on SI were applied to all datasets included the validation one of March, obtaining satisfying predicted results as shown in Table 15.
Although the proposed method seems promising and predictive functions estimates satisfying, further investigations must be done, with the man aim of improving generalization capability of models and testing their dependencies on other operational conditions like different areas, oblique acquisitions and GCP spatial distribution.