3D Calibration Test-Field for Digital Cameras Mounted on Unmanned Aerial Systems (UAS)

: Due to the large number of technological developments in recent years, UAS systems are now used for monitoring purposes and in projects with high precision demand, such as 3D model-based creation of dams, reservoirs, historical monuments etc. These unmanned systems are usually equipped with an automatic pilot device and a digital camera (photo/video, multispectral, Near Infrared etc.), of which the lens has distortions; but this can be determined in a calibration process. Currently, a method of “self-calibration” is used for the calibration of the digital cameras mounted on UASs, but, by using the method of calibration based on a 3D calibration object, the accuracy is improved in comparison with other methods. Thus, this paper has the objective of establishing a 3D calibration ﬁeld for the digital cameras mounted on UASs in terms of accuracy and robustness, being the largest reported publication to date. In order to test the proposed calibration ﬁeld, a digital camera mounted on a low-cost UAS was calibrated at three different heights: 23 m, 28 m, and 35 m, using two conﬁgurations for image acquisition. Then, a comparison was made between the residuals obtained for a number of 100 Check Points (CPs) using self-calibration and test-ﬁeld calibration, while the number of Ground Control Points (GCPs) variedand the heights were interchanged. Additionally, the parameters where tested on an oblique ﬂight done 2 years before calibration, in manual mode at a medium altitude of 28 m height. For all tests done in the case of the double grid nadiral ﬂight, the parameters calculated with the proposed 3D ﬁeld improved the results by more than 50% when using the optimum and a large number of GCPs, and in all analyzed cases with 75% to 95% when using a minimum of 3 GCP. In this context, it is necessary to conduct accurate calibration in order to increase the accuracy of the UAS projects, and also to reduce ﬁeld measurements.


Introduction
The calibration process of a digital camera is necessary to obtain metric information of the three-dimensional world using two-dimensional images. The purpose of camera calibration is to describe the projection model that relates both coordinate systems, i.e., the world and image coordinate systems, and to identify intrinsic camera parameters so that it can be used as a measurement device [1].
The experiments conducted by Samper et al. [2] for synthetic and real point calibration models demonstrate that the model using a planar target produces a larger reconstruction error of the world system coordinates than the model using 3D calibration objects. As we previously demonstrated [3][4][5],

Related Work
There are many studies describing the pre-calibration of digital cameras used for UAS mapping purposes during the last decade, and only a few with the camera mounted on a UAS and at a long object distance during calibration process. Some of them have been calibrated using planar calibration objects [13][14][15][16], some using a 3D target, but at short object distance [17][18][19], and only a few tested the impact of camera pre-calibration on UAS mapping projects [10,20].
Test fields have been used to calibrate and test digital photogrammetric sensors for over two decades [21]. The test-field designed by Perez et al. [15] has a set of 67 target points distributed over a flat surface of approximately 25 m × 25 m, but has not been tested in an UAS photogrammetric project. Several control points attached on two facades in an inner courtyard and on the ground were used as test-fields for UAS camera calibration by Kraft et al. [22]; however, the resulting calibration parameters have not been replaced in real projects. A 3-dimensional test field was designed on a building wall by Cramer et al. [23] with an extension of approximately 14 m × 3 m × 3 m and approx. 500 spatially-distributed coded and non-coded targets, but only the values of focal length obtained by test-field and in-situ calibration were compared. They concluded that for unconventional block geometry, a best pre-calibrated camera must be requested, and that, possibly, it is helpful to use the pre-calibration parameters as an approximation for in-situ calibration. A 2500-square-meter calibration field was designed by Yusoff et al. [20] on a football field with 36 calibration targets. Different heights for calibration on different object distances where tested, but using only three check lines for image mapping measurement accuracy. They found that a different height for calibration process was the optimum for a given object distance. In Harwin et al. [10] the camera pre-calibration performed using CalibCam and PhotoScan software based on a 25 m × 25 m test-field and 13 GCPs and using Lens based on a Checker Board pattern didn't improved the results in comparison with self-calibration. They suggest that there are grounds for undertaking a similar study over a larger area in order to produce more 'scalable' rules for camera calibration and GCP distribution.

Proposed 3D Calibration Field
Our study builds on previous work by quantifying the impact of camera calibration on the accuracy of UAS reconstruction projects. Furthermore, our study is a first attempt to assess the impact of pre-calibration parameters on real UAS reconstruction projects, both with nadiral single and double grid images and with oblique convergent images, using also a large number of GCPs i.e., 100 instead of a few check points or distances. So, the actual question we want to answer is: using test-field calibration, how much we can improve the accuracy of UAS reconstruction projects with different flight configuration while reducing manual effort.
Plates made of plexiglass whose center have been manually measured have been used for the ground materialization of GCPs. If coded targets are used for control-point materialization, the size of the circular target must be 10 pixels on the ground to automatically identify the center, as recommended by [23][24][25]. As in this study, a digital camera with a 4 mm focal length was used, the size of the circular target at 35 m height must be 43 cm. Thus, taking into consideration the large number of control points of the proposed test-field and the large dimension of each control point, the non-coded targets option was chosen. This way, the user is able to watch and verify the correct location and numbering of each control point.
A local geodetic network was designed for the 3D calibration test-field to minimize the errors during the map conversion process, and a local reference coordinate system for the X, Y, and Z coordinates was used to determine the GCPs positions. In this way, a homogeneous and unitary precision of all three components of spatial positioning was obtained. So, this study is the first to undertake a calibration field of such dimensions and by using a local coordinate system to determine the GCPs coordinates with millimeter precision, independent of deformations caused by the global geodetic system.

Study Area
The study area of about 1 ha is located near the Faculty of Hydro-technical Engineering, Geodesy and Environmental Engineering from "Gheorghe Asachi" Technical University of Iasi, Romania, covering the building roof, the parking lot and the green area in the vicinity.

Marking the Ground Control Points (GCPs)
Firstly, 349 uniformly distributed points in the study area were marked: 104 points on the roof (marked with paint, using a plastic pattern), 55 points distributed on 11 rows on the building façade being marked using a construction black marker, 29 points in the green space (reinforced concrete poles to ensure different heights having a metal bold in the superior part), 2 points from the GNSS network (concrete), and 159 points in the parking lot (metallic bolts) (Figure 1).
During the flight, the points marked on the façade and the roof were marked with A3 laminated sheets (Figure 1a-c); after having drawn two black and white triangles, their intersection represented the mathematical point.

GNSS Measurements
The solution to design a spatial geodetic network determined by GNSS technology has been chosen, providing a homogeneous and unitary precision of all three components of spatial positioning, currently representing the most efficient approach in the development of geodetic networks in ground measurement works.
The spatial geodetic network was designed with 4 points (A, B, C, D), two of which (A, B) were grounded in the green area behind the Faculty of Hydrotechnical Engineering, Geodesy, and Environmental Engineering; the other two (C, D) were located on the rooftop, in the northern part of the building at a height difference of about 9 m (Figure 3). In the foregoing research, a local reference and coordinate system was adopted.
The framework of GNSS measurements, is represented by the World Geodetic System (WGS84). The WGS84 is an Earth-centered, Earth-fixed terrestrial reference system and geodetic datum, having three spatial rectangular axes defined in relation with a rotation ellipsoid determined exclusively using data, techniques and satellite technologies.
The GNSS measurements effectuated in the spatial geodetic network were made by static relative positioning method. Three GNSS South S82-V receivers with double frequency GPS + GLONASS + GALILEO + COMPASS + SBAS have been used, having an horizontal precision of positioning of 3 mm ± 0.3 ppm and a vertical precision of 4 mm ± 0.5 ppm, in static mode.
The duration of an observation session is determined by the length of the polygon sides, the number of visible satellites, the geometry of their space position (as assessed by PDOP-Position Dilution of Precision) and the accuracy requirements imposed when determining the network points. The points were situated in the parking lot and the green space, have been marked with plates made by orange plexiglass and two black triangles, 3 mm thickness, 40 cm × 40 cm dimension, with a hole in the center of 5 mm diameter.
A perspective view over the 3D calibration field of digital cameras mounted on UASs, can be seen in Figure 2.

GNSS Measurements
The solution to design a spatial geodetic network determined by GNSS technology has been chosen, providing a homogeneous and unitary precision of all three components of spatial positioning, currently representing the most efficient approach in the development of geodetic networks in ground measurement works.
The spatial geodetic network was designed with 4 points (A, B, C, D), two of which (A, B) were grounded in the green area behind the Faculty of Hydrotechnical Engineering, Geodesy, and Environmental Engineering; the other two (C, D) were located on the rooftop, in the northern part of the building at a height difference of about 9 m (Figure 3). In the foregoing research, a local reference and coordinate system was adopted.
The framework of GNSS measurements, is represented by the World Geodetic System (WGS84). The WGS84 is an Earth-centered, Earth-fixed terrestrial reference system and geodetic datum, having three spatial rectangular axes defined in relation with a rotation ellipsoid determined exclusively using data, techniques and satellite technologies.
The GNSS measurements effectuated in the spatial geodetic network were made by static relative positioning method. Three GNSS South S82-V receivers with double frequency GPS + GLONASS + GALILEO + COMPASS + SBAS have been used, having an horizontal precision of positioning of 3 mm ± 0.3 ppm and a vertical precision of 4 mm ± 0.5 ppm, in static mode.
The duration of an observation session is determined by the length of the polygon sides, the number of visible satellites, the geometry of their space position (as assessed by PDOP-Position

GNSS Measurements
The solution to design a spatial geodetic network determined by GNSS technology has been chosen, providing a homogeneous and unitary precision of all three components of spatial positioning, currently representing the most efficient approach in the development of geodetic networks in ground measurement works.
The spatial geodetic network was designed with 4 points (A, B, C, D), two of which (A, B) were grounded in the green area behind the Faculty of Hydrotechnical Engineering, Geodesy, and Environmental Engineering; the other two (C, D) were located on the rooftop, in the northern part of the building at a height difference of about 9 m (Figure 3). In the foregoing research, a local reference and coordinate system was adopted.
Due to distances of less than 50 m between the network points, the duration of a work session was at least 20 min, ranging from 26 min to 1 h and 50 min, i.e., the time period required to resolve ambiguities from phase measurements. After data collection, the receivers were switched off and a part of them were moved as in the measurements scheme on the following bases (A-B-D, A-D-C, D-C-B). The measurements continued afterward ( Figure 4).  The measured GNSS bases, fell within the range PDOP 2.2 ÷ 4.2 (<5), HDOP 1.1 ÷ 1.7 (<2) and VDOP 1.8 ÷ 3.9 (<4), which proves a good configuration of satellite geometry with effect on positioning accuracy of GNSS network points.

Processing the GNSS Measurements
Observations from all sessions were completed in one day (28 September 2017); this is why the precise ephemeris downloaded from the IGS web site were used for processing.
Four stations with unknown coordinates resulted from the network pre-analysis, (12 components of the three-dimensional rectangular coordinates), 9 GNSS bases (27 relative spatial coordinates), and 27 -12 = 15 network degrees of freedom.
In the case of static measurements, there are several GNSS loops interconnected. For each formed polygon, the algebraic sum of the GNSS base relative coordinates components (ΔX, ΔY, ΔZ) should be zero. A large non-closure in one of the polygons indicates that there is an error of one or more GNSS bases measured in that polygon. In the non-closure analysis report of polygons formed within the GNSS network (Table 1), it can be seen that they are less than 7 mm in planimetric position, and 3 cm in vertical position, with total relative errors in the range 30 ÷ 292 ppm.
The GNSS network adjustment was performed in the WGS84 three-dimensional coordinate system, under the condition of a free geodetic network, in which no element is considered to be fixed, so, a network in which only the measurements to determine the network geometry is established. The adjustment results are shown in Appendix A.
For the current project, a heterogeneous "LOCAL HGIM" reference and coordinate system (CRS) was adopted, consisting of: The framework of GNSS measurements, is represented by the World Geodetic System (WGS84). The WGS84 is an Earth-centered, Earth-fixed terrestrial reference system and geodetic datum, having three spatial rectangular axes defined in relation with a rotation ellipsoid determined exclusively using data, techniques and satellite technologies.
The GNSS measurements effectuated in the spatial geodetic network were made by static relative positioning method. Three GNSS South S82-V receivers with double frequency GPS + GLONASS + GALILEO + COMPASS + SBAS have been used, having an horizontal precision of positioning of 3 mm ± 0.3 ppm and a vertical precision of 4 mm ± 0.5 ppm, in static mode.
The duration of an observation session is determined by the length of the polygon sides, the number of visible satellites, the geometry of their space position (as assessed by PDOP-Position Dilution of Precision) and the accuracy requirements imposed when determining the network points. Due to distances of less than 50 m between the network points, the duration of a work session was at least 20 min, ranging from 26 min to 1 h and 50 min, i.e., the time period required to resolve ambiguities from phase measurements. After data collection, the receivers were switched off and a part of them were moved as in the measurements scheme on the following bases (A-B-D, A-D-C, D-C-B). The measurements continued afterward ( Figure 4).   The measured GNSS bases, fell within the range PDOP 2.2 ÷ 4.2 (<5), HDOP 1.1 ÷ 1.7 (<2) and VDOP 1.8 ÷ 3.9 (<4), which proves a good configuration of satellite geometry with effect on positioning accuracy of GNSS network points.

Processing the GNSS Measurements
Observations from all sessions were completed in one day (28 September 2017); this is why the precise ephemeris downloaded from the IGS web site were used for processing.
Four stations with unknown coordinates resulted from the network pre-analysis, (12 components of the three-dimensional rectangular coordinates), 9 GNSS bases (27 relative spatial coordinates), and 27 -12 = 15 network degrees of freedom.
In the case of static measurements, there are several GNSS loops interconnected. For each formed polygon, the algebraic sum of the GNSS base relative coordinates components (ΔX, ΔY, ΔZ) should be zero. A large non-closure in one of the polygons indicates that there is an error of one or more GNSS bases measured in that polygon. In the non-closure analysis report of polygons formed within the GNSS network (Table 1), it can be seen that they are less than 7 mm in planimetric position, and The measured GNSS bases, fell within the range PDOP 2.2 ÷ 4.2 (<5), HDOP 1.1 ÷ 1.7 (<2) and VDOP 1.8 ÷ 3.9 (<4), which proves a good configuration of satellite geometry with effect on positioning accuracy of GNSS network points.

Processing the GNSS Measurements
Observations from all sessions were completed in one day (28 September 2017); this is why the precise ephemeris downloaded from the IGS web site were used for processing.
Four stations with unknown coordinates resulted from the network pre-analysis, (12 components of the three-dimensional rectangular coordinates), 9 GNSS bases (27 relative spatial coordinates), and 27 -12 = 15 network degrees of freedom. In the case of static measurements, there are several GNSS loops interconnected. For each formed polygon, the algebraic sum of the GNSS base relative coordinates components (∆X, ∆Y, ∆Z) should be zero. A large non-closure in one of the polygons indicates that there is an error of one or more GNSS bases measured in that polygon. In the non-closure analysis report of polygons formed within the GNSS network (Table 1), it can be seen that they are less than 7 mm in planimetric position, and 3 cm in vertical position, with total relative errors in the range 30 ÷ 292 ppm.
The GNSS network adjustment was performed in the WGS84 three-dimensional coordinate system, under the condition of a free geodetic network, in which no element is considered to be fixed, so, a network in which only the measurements to determine the network geometry is established. The adjustment results are shown in Appendix A.
For the current project, a heterogeneous "LOCAL HGIM" reference and coordinate system (CRS) was adopted, consisting of: • CRS 1: planimetric rectangular coordinates (x,y) in double tangent stereographical projection "Local-HGIM" with WGS84 datum, • CRS 2: ortometric altitude (H) in local vertical datum with the horizontal reference plan which pases through origin, having the altitude equal with the ellipsoidal height (because of small distances).
The "Local-HGIM" stereographical projection parameters are presented in Table 2, and the adjusted coordinates of GNSS network points in the new heterogenious reference and coordinate system named "LOCAL HGIM", together with accuracy assessment, are presented in Appendix A.   Table 2. "Local-HGIM" stereographical projection parameters.

Measurement of Ground Control Points (GCPs)
The detail points where determined by topographic measurements using a Leica TCR 405 with an angular precision of 5" (0.5 mgon) and distance precision of 2 mm + 2 ppm in Fine-Mode (IR). Each GCP have been double measured with a miniprism verticalized with a buble (Figure 5a) from the ends of a GNSS base (A-B for the points situated in the parking lot and the green area (Exhibit 5b) and C-D for the points situated on the roof).
As a preliminary check, the differences were made between the two coordinates strings thus obtained, resulting in a mean absolute value of 9 mm in the horizontal position and 1 mm in the vertical position for the parking and green area points, and 3 mm in the horizontal position and 6 mm in vertical position for the roof points.
The detail points where determined by topographic measurements using a Leica TCR 405 with an angular precision of 5″ (0.5 mgon) and distance precision of 2 mm + 2 ppm in Fine-Mode (IR). Each GCP have been double measured with a miniprism verticalized with a buble (Figure 5a) from the ends of a GNSS base (A-B for the points situated in the parking lot and the green area (Exhibit 5b) and C-D for the points situated on the roof). As a preliminary check, the differences were made between the two coordinates strings thus obtained, resulting in a mean absolute value of 9 mm in the horizontal position and 1 mm in the vertical position for the parking and green area points, and 3 mm in the horizontal position and 6 mm in vertical position for the roof points.
For the GCPs coordinate adjustments, the preliminary reduction of slant distances and zenithal angles to the direction between points at ground level was necessary. For the GCPs coordinate adjustments, the preliminary reduction of slant distances and zenithal angles to the direction between points at ground level was necessary.
At the end of this stage, the initial coordinates of the GCPs were obtained as a simple arithmetic mean between the two above-mentioned determinations, to be included further in the adjustment process, applying the least squares method.

Geodetic Measurements Processing
The adjustment of the topographic data was made in an application that was specifically developed in the Matlab, having as inputs the GNSS base points coordinates, the initial coordinates of the GCPs, and the topographic elements measured on the ground. To obtain the adjusted spatial coordinates of the GCPs, for each individual point, the adjustment functional model of indirect weight method was applied [25]. The total spatial positioning error for the whole project was less than 5 mm.

UAS Image Acquisition over the Proposed 3D Calibration Field
For this study, a low-cost UAS "DJI Phantom 3 Standard" was used, with an integrated "GoPro" digital camera (12 Mega pixel) equipped with a 6.31748 mm by 4.73811 mm image sensor. The digital image has a resolution of 4000 × 3000 pixels with a 1.6245 µm pixel size. UAS images were taken over the calibration field from 3 different heights: 23 m, 28 m, and 35 m, using the "Pix4D Capture" software for flight planning and control of the UAS system during the flight mission. The "circular" mode for image acquisition was selected for the area of interest (Figure 6 a, At the end of this stage, the initial coordinates of the GCPs were obtained as a simple arithmetic mean between the two above-mentioned determinations, to be included further in the adjustment process, applying the least squares method.

Geodetic Measurements Processing
The adjustment of the topographic data was made in an application that was specifically developed in the Matlab, having as inputs the GNSS base points coordinates, the initial coordinates of the GCPs, and the topographic elements measured on the ground. To obtain the adjusted spatial coordinates of the GCPs, for each individual point, the adjustment functional model of indirect weight method was applied [25]. The total spatial positioning error for the whole project was less than 5 mm.

UAS Image Acquisition over the Proposed 3D Calibration Field
For this study, a low-cost UAS "DJI Phantom 3 Standard" was used, with an integrated "GoPro" digital camera (12 Mega pixel) equipped with a 6.31748 mm by 4.73811 mm image sensor. The digital image has a resolution of 4000 × 3000 pixels with a 1.6245 μm pixel size.
UAS images were taken over the calibration field from 3 different heights: 23 m, 28 m, and 35 m, using the "Pix4D Capture" software for flight planning and control of the UAS system during the flight mission. The "circular" mode for image acquisition was selected for the area of interest ( Figure  6 a,b).
The image configuration for the camera calibration is close to the one described in [9]. Here, it is mentioned that the test field should be imaged perpendicularly and obliquely, and each image should have a relative rotation of 90° around the optical axis, with eight images being sufficient. As the DJI system has an integrated camera which cannot be separated from the drone itself, taking rotated images in vertical plane around the optical axis is impossible. So, 17 images were taken in an oblique position, the camera being orientated with the optical axis towards the geometric center. Additionally, 4 images in nadiral position were acquired in manual mode with the camera rotated each time at 90° around the optical axis, approximatively above the geometric center, used primarily to determine the principal point coordinates and affinity parameters.
The calibration was carried out as part of a bundle block adjustment using the "Pix 4D Mapper" and "3DF Zephyr Pro" software.
The workflow of this study is presented in Figure 7. The image configuration for the camera calibration is close to the one described in [9]. Here, it is mentioned that the test field should be imaged perpendicularly and obliquely, and each image should have a relative rotation of 90 • around the optical axis, with eight images being sufficient. As the DJI system has an integrated camera which cannot be separated from the drone itself, taking rotated images in vertical plane around the optical axis is impossible. So, 17 images were taken in an oblique position, the camera being orientated with the optical axis towards the geometric center. Additionally, 4 images in nadiral position were acquired in manual mode with the camera rotated each time at 90 • around the optical axis, approximatively above the geometric center, used primarily to determine the principal point coordinates and affinity parameters.

UAS Images Processing for Camera Calibration
The camera calibration process was used to obtain intrinsic parameters such as focal distance (f ), optical point coordinates (u 0 ,v 0 ), correction of radial distortion (k 1 , k 2 , k 3 ), correction of decentering distortion (p 1 , p 2 ) [1] as well as the extrinsic parameters (r ij , X 0 , Y 0 , Z 0 ).
The calibration was carried out as part of a bundle block adjustment using the "Pix 4D Mapper" and "3DF Zephyr Pro" software.
The workflow of this study is presented in Figure 7. To process the UAS images in "Pix4D Mapper" software, the desired images were imported into a new project, and regarding the "Output/GCP Coordinate System", the "Arbitrary coordinate system" was selected, the "Image Geolocation" being taken from the EXIF information. With "GCP/MTP Manager", the file with the control points coordinates was imported, the horizontal accuracy was changed to 0.002, and the vertical accuracy (m) to 0.001, the "Type" of GCP being left as the default, namely "3D GCP". First, 9 GCP were manually measured on each image using "Basic Editor" for a preliminary image block orientation through BBA.
As an advanced processing option, the geometrically verified matching was utilized for the matching strategy, and the "standard" option was selected with all internal and external parameters optimization for the calibration method. After performing the BBA, the images have been oriented and the rest of the 340 GCPs were measured easily on each image as they appeared, the approximate position being indicated by the software after performing a ray intersection (Figure 8). Even so, the GCP marking remains time-consuming for the operator (approximately 8 h). To process the UAS images in "Pix4D Mapper" software, the desired images were imported into a new project, and regarding the "Output/GCP Coordinate System", the "Arbitrary coordinate system" was selected, the "Image Geolocation" being taken from the EXIF information. With "GCP/MTP Manager", the file with the control points coordinates was imported, the horizontal accuracy was changed to 0.002, and the vertical accuracy (m) to 0.001, the "Type" of GCP being left as the default, namely "3D GCP". First, 9 GCP were manually measured on each image using "Basic Editor" for a preliminary image block orientation through BBA.
As an advanced processing option, the geometrically verified matching was utilized for the matching strategy, and the "standard" option was selected with all internal and external parameters optimization for the calibration method. After performing the BBA, the images have been oriented and the rest of the 340 GCPs were measured easily on each image as they appeared, the approximate position being indicated by the software after performing a ray intersection (Figure 8). Even so, the GCP marking remains time-consuming for the operator (approximately 8 h).
A final BBA ran and a quality report was created, as well as a graphic representation (Figure 9). The nadiral images haven't been oriented automatically by the software, so a minimum of 7 GCP was measured on each unoriented image. All GCPs that appear on each nadiral image could be measured after processing again the project. matching strategy, and the "standard" option was selected with all internal and external parameters optimization for the calibration method. After performing the BBA, the images have been oriented and the rest of the 340 GCPs were measured easily on each image as they appeared, the approximate position being indicated by the software after performing a ray intersection (Figure 8). Even so, the GCP marking remains time-consuming for the operator (approximately 8 h). A final BBA ran and a quality report was created, as well as a graphic representation (Figure 9). The nadiral images haven't been oriented automatically by the software, so a minimum of 7 GCP was measured on each unoriented image. All GCPs that appear on each nadiral image could be measured after processing again the project. The principal elements of the quality report for each calibration process for the 3 different heights and the two configurations for image acquisition, 21 images and 17 images respectively, can be found in Table 3. The values of interior orientation parameters for each case are listed in Table 4, where f is the focal length, u0 and v0 are the coordinates of the principal point, k1, k2 and k3 are the coefficients of the radial distortion, and p1 and p2 are the coefficients of the tangential distortion The same measurements for GCPs image coordinates have been used into "3DF Zephyr Pro" software. The principal elements of the quality report for each calibration process for the 3 different heights and the two configurations for image acquisition, 21 images and 17 images respectively, can be found in Table 3.
The values of interior orientation parameters for each case are listed in Table 4, where f is the focal length, u 0 and v 0 are the coordinates of the principal point, k 1 , k 2 and k 3 are the coefficients of the radial distortion, and p 1 and p 2 are the coefficients of the tangential distortion.  The same measurements for GCPs image coordinates have been used into "3DF Zephyr Pro" software.
To model the radial distortion ∆r, the odd-order polynomial ∆r = K 1 r 3 + K 2 r 5 + K 3 r 7 [26][27][28], where r is the radial distance, was used and the radial distortion profiles for all three heights, and for the two image configuration, as shown in Figure 10. To model the radial distortion Δr, the odd-order polynomial Δr = K1r 3 + K2r 5 + K3r 7 [26,28], where r is the radial distance, was used and the radial distortion profiles for all three heights, and for the two image configuration, as shown in Figure 10. The decentering distortion is due to the lens elements decentering along the optical axis [28], and was graphically represented in an analogous manner to radial distortion, using the function: (1)

Camera Calibration Using "3DF Zephyr Pro" Software
For UAS image processing using the "3DF Zephyr Pro" software, the images were imported into the software and the project was processed without specifying the calibration parameters, so for the interior and exterior orientation parameters for each camera position, an approximation was made based on the EXIF information. After performing SfM process, the image coordinates of each GCP were imported, a preliminary check for each one being made based on the reprojection error to see if all GCPs have been imported correctly. Then, the 3D coordinates have been added. Finally, the "constraint" and "control" boxes have been checked for all 349 GCPs, and a bundle adjustment process was performed adjusting all the interior orientation parameters. After performing BBA, all images have been automatically oriented, and the parameters of interior and exterior orientation have been calculated (Figure 11). If an additional measurement for a GCP is to be made, the correct position on the image is found by a ray intersection ( Figure 12).
The principal elements of the quality report for each calibration process done for the 3 different heights in "3DF Zephyr Pro" software and the two configurations for image acquisition, 21 images and 17 images respectively, are listed in Table 5.
The values of interior orientation parameters for each case are listed in Table 6. The radial and decentering distortion profiles drawn using the coefficients calculated with "3DF Zephyr Pro" software are shown in Figure 13. The decentering distortion is due to the lens elements decentering along the optical axis [28], and was graphically represented in an analogous manner to radial distortion, using the function:

Camera Calibration Using "3DF Zephyr Pro" Software
For UAS image processing using the "3DF Zephyr Pro" software, the images were imported into the software and the project was processed without specifying the calibration parameters, so for the interior and exterior orientation parameters for each camera position, an approximation was made based on the EXIF information. After performing SfM process, the image coordinates of each GCP were imported, a preliminary check for each one being made based on the reprojection error to see if all GCPs have been imported correctly. Then, the 3D coordinates have been added. Finally, the "constraint" and "control" boxes have been checked for all 349 GCPs, and a bundle adjustment process was performed adjusting all the interior orientation parameters. After performing BBA, all images have been automatically oriented, and the parameters of interior and exterior orientation have been calculated (Figure 11). If an additional measurement for a GCP is to be made, the correct position on the image is found by a ray intersection ( Figure 12).
The principal elements of the quality report for each calibration process done for the 3 different heights in "3DF Zephyr Pro" software and the two configurations for image acquisition, 21 images and 17 images respectively, are listed in Table 5.
The values of interior orientation parameters for each case are listed in Table 6. The radial and decentering distortion profiles drawn using the coefficients calculated with "3DF Zephyr Pro" software are shown in Figure 13

UAS Image Acquisition for Testing the Proposed 3D Calibration Field
In order to test the obtained interior orientation parameters, UAS images were taken with the "DJI Phantom 3 Standard" platform over the calibration field at the same heights as in the camera calibration processes, immediately after the calibration flights. The flight planning was made with "Pix4D Capture" software ( Figure 14), choosing the longitudinal and transversal overlap of 80% and 60% respectively, the camera being oriented in nadiral position. For 28 m height, the flight was made in double grid and 122 images were acquired, while for 23 m and 35 m, the flight was made in single grid and 85 and 51 images were acquired respectively.

UAS Image Acquisition for Testing the Proposed 3D Calibration Field
In order to test the obtained interior orientation parameters, UAS images were taken with the "DJI Phantom 3 Standard" platform over the calibration field at the same heights as in the camera calibration processes, immediately after the calibration flights. The flight planning was made with "Pix4D Capture" software ( Figure 14), choosing the longitudinal and transversal overlap of 80% and 60% respectively, the camera being oriented in nadiral position. For 28 m height, the flight was made in double grid and 122 images were acquired, while for 23 m and 35 m, the flight was made in single grid and 85 and 51 images were acquired respectively.

UAS Image Processing for Testing the Proposed 3D Calibration Field
A total of 63 scenarios were tested. The assessed key variables were: (i) the influence of precalibration process on the accuracy of BBA in comparison with self-calibration; (ii) the inclusion or exclusion of nadiral images from the image blocks used for the calibration processes; (iii) the influence of height between the calibration process and the UAS reconstruction project; and (iv) the number of GCPs. In order to process the UAS images, the coordinates of 150 points were selected from the total of 349 targets representing the calibration field, 50 of them being considered GCPs and the rest of 100 as Check Points. The distribution of these points can be found in Figure 15. Before processing UAS images with 3DF Zephyr Pro different scenarios were tested using both "Pix 4D Mapper" and "3DF Zephyr Pro" software. In both cases, a minimum of 3 GCPs were

UAS Image Processing for Testing the Proposed 3D Calibration Field
A total of 63 scenarios were tested. The assessed key variables were: (i) the influence of precalibration process on the accuracy of BBA in comparison with self-calibration; (ii) the inclusion or exclusion of nadiral images from the image blocks used for the calibration processes; (iii) the influence of height between the calibration process and the UAS reconstruction project; and (iv) the number of GCPs. In order to process the UAS images, the coordinates of 150 points were selected from the total of 349 targets representing the calibration field, 50 of them being considered GCPs and the rest of 100 as Check Points. The distribution of these points can be found in Figure 15.

UAS Image Acquisition for Testing the Proposed 3D Calibration Field
In order to test the obtained interior orientation parameters, UAS images were taken with the "DJI Phantom 3 Standard" platform over the calibration field at the same heights as in the camera calibration processes, immediately after the calibration flights. The flight planning was made with "Pix4D Capture" software ( Figure 14), choosing the longitudinal and transversal overlap of 80% and 60% respectively, the camera being oriented in nadiral position. For 28 m height, the flight was made in double grid and 122 images were acquired, while for 23 m and 35 m, the flight was made in single grid and 85 and 51 images were acquired respectively.

UAS Image Processing for Testing the Proposed 3D Calibration Field
A total of 63 scenarios were tested. The assessed key variables were: (i) the influence of precalibration process on the accuracy of BBA in comparison with self-calibration; (ii) the inclusion or exclusion of nadiral images from the image blocks used for the calibration processes; (iii) the influence of height between the calibration process and the UAS reconstruction project; and (iv) the number of GCPs. In order to process the UAS images, the coordinates of 150 points were selected from the total of 349 targets representing the calibration field, 50 of them being considered GCPs and the rest of 100 as Check Points. The distribution of these points can be found in Figure 15. Before processing UAS images with 3DF Zephyr Pro different scenarios were tested using both "Pix 4D Mapper" and "3DF Zephyr Pro" software. In both cases, a minimum of 3 GCPs were Before processing UAS images with 3DF Zephyr Pro different scenarios were tested using both "Pix 4D Mapper" and "3DF Zephyr Pro" software. In both cases, a minimum of 3 GCPs were manually measured on each image in which they appeared for a preliminary image block orientation through BBA. Then, the rest of GCPs were measured easily, the approximate position being indicated by the software after performing a ray intersection. We have concluded that even if a point has the "type" changed to "Check Point", it is included in the BBA. As we wanted to have a set of control points independent of BBA, we tried to make two separate lists of points: one for GCPs, and one for check points (CPs). When using a minimum number of GCPs, the "Pix 4D Mapper" is absolutely wrong when computing the images orientation, so we further replaced the calibration parameters into "3DF Zephyr Pro".

Quality Assessment
The residuals of the BBA are computed by comparing the CPs coordinates with those determined with high precision followed by the determination of the root mean square error (RMSE). So, the software is using only the points marked as "control" in the BBA process. A different number was used for GCPs: the minimum of 3, the optimum of 14, as found in [29], and a large number of 50.
Although the recovered camera parameters with both software were quite similar, by analyzing the results listed in Table 7, we can see that the errors are smaller in the case of using the calibration parameters calculated with "3DF Zephyr Pro" software. Therefore, we have done further tests using only these parameters. Table 7. The residuals of 100 CPs for a 28 m nadiral flight when using self-calibration and pre-calibration at 28 m height done with "Pix 4D Mapper" respectively "3DF Zephyr". When using self-calibration with the minimum number of GCPs, the errors are very large because the image block is affected by the bowl effects (Figure 16a). By modifying the camera interior orientation parameters before running the SfM, the distortions are corrected (Figure 16b). manually measured on each image in which they appeared for a preliminary image block orientation through BBA. Then, the rest of GCPs were measured easily, the approximate position being indicated by the software after performing a ray intersection. We have concluded that even if a point has the "type" changed to "Check Point", it is included in the BBA. As we wanted to have a set of control points independent of BBA, we tried to make two separate lists of points: one for GCPs, and one for check points (CPs). When using a minimum number of GCPs, the "Pix 4D Mapper" is absolutely wrong when computing the images orientation, so we further replaced the calibration parameters into "3DF Zephyr Pro".

Quality Assessment
The residuals of the BBA are computed by comparing the CPs coordinates with those determined with high precision followed by the determination of the root mean square error (RMSE). So, the software is using only the points marked as "control" in the BBA process. A different number was used for GCPs: the minimum of 3, the optimum of 14, as found in [29], and a large number of 50.
Although the recovered camera parameters with both software were quite similar, by analyzing the results listed in Table 7, we can see that the errors are smaller in the case of using the calibration parameters calculated with "3DF Zephyr Pro" software. Therefore, we have done further tests using only these parameters. Table 7. The residuals of 100 CPs for a 28 m nadiral flight when using self-calibration and precalibration at 28 m height done with "Pix 4D Mapper" respectively "3DF Zephyr".

Height for Pre-Calibration (m)
No. of GCPs When using self-calibration with the minimum number of GCPs, the errors are very large because the image block is affected by the bowl effects (Figure 16a). By modifying the camera interior orientation parameters before running the SfM, the distortions are corrected (Figure 16b). In the case of 23 m nadiral block images, not all 100 CPs could be measured on a minimum of 2 images, so only 87 CPs were used to compute the residuals, as shown in Table 8. Also, only 44 GCPs have been used as constraints in the BBA process.
The residuals for the 28 m nadiral flight have been calculated using 100 CPs, and are listed in Table 9. In the case of the 35 m nadiral block images, only 99 CPs were used to compute the residuals, as shown in Table 10.
In Tables 8-10, the residuals have also been listed in relation with the value of GSD calculated as a relationship between the flight height and the focal length multiplied by the sensor size. So, for In the case of 23 m nadiral block images, not all 100 CPs could be measured on a minimum of 2 images, so only 87 CPs were used to compute the residuals, as shown in Table 8. Also, only 44 GCPs have been used as constraints in the BBA process.
The residuals for the 28 m nadiral flight have been calculated using 100 CPs, and are listed in Table 9. In the case of the 35 m nadiral block images, only 99 CPs were used to compute the residuals, as shown in Table 10.
In Tables 8-10, the residuals have also been listed in relation with the value of GSD calculated as a relationship between the flight height and the focal length multiplied by the sensor size. So, for the 23 m height, a value of 0.93 cm was obtained for the GSD; for the 28 m height, a value of 1.1 cm, and for the 35 m, a value of 1.4 cm was obtained.

Testing the Pre-Calibration Parameters with an Oblique Flight
In order to test the pre-calibration parameters on an oblique flight, we chose a UAS reconstruction project done for an old monastery two years before camera calibration. Five reference points were uniformly distributed around the monastery, allowing good visibility without any blockage from vegetation. Three of them have been used as GCPs in order to place the project into a desired coordinate system, i.e., "Stereographical on unique secant plane 1970" and two of them as check points (Figure 17).
In order to test the pre-calibration parameters on an oblique flight, we chose a UAS reconstruction project done for an old monastery two years before camera calibration. Five reference points were uniformly distributed around the monastery, allowing good visibility without any blockage from vegetation. Three of them have been used as GCPs in order to place the project into a desired coordinate system, i.e., "Stereographical on unique secant plane 1970" and two of them as check points (Figure 17). The reference points were made by plexiglass, having the center marked by the intersection of two black triangles with 30 cm × 30 cm dimensions. The coordinates were measured with high accuracy using the GNSS technology, the final coordinates resulting as an average of two determinations. For UAS image acquisition, the flight was made in manual mode a total of 30 images from 30 different camera positions being acquired in a circular sequence with the camera axis oriented approximately at 45 degrees to perform oblique acquisition for the facades. Taking into account that the flight was done in manual mode, the altitudes above ground for each camera positions are different, ranging from 25 m to 30 m, with an average of 28 m. Table 11 summarizes the RMSE for the BBA process, and Table 12 presents the individual residuals for each CP.  The reference points were made by plexiglass, having the center marked by the intersection of two black triangles with 30 cm × 30 cm dimensions. The coordinates were measured with high accuracy using the GNSS technology, the final coordinates resulting as an average of two determinations. For UAS image acquisition, the flight was made in manual mode a total of 30 images from 30 different camera positions being acquired in a circular sequence with the camera axis oriented approximately at 45 degrees to perform oblique acquisition for the facades. Taking into account that the flight was done in manual mode, the altitudes above ground for each camera positions are different, ranging from 25 m to 30 m, with an average of 28 m. Table 11 summarizes the RMSE for the BBA process, and Table 12 presents the individual residuals for each CP. Table 11. The residuals of 2 CPs resulted after BBA process for an oblique flight when using self-calibration and pre-calibration at 3 different heights: 23 m, 28 m and 35 m.

Discussion
For this study, a 3D calibration test-field was designed and tested, being the largest reported to date. A number of 349 non-coded targets was measured with millimeter precision in a designed local coordinate system. The camera was calibrated at three different flying heights with two configurations for image acquisition, and the tests were performed on nadiral and oblique flights done at three heights, i.e., 23 m, 28 m, and 35 m, the same as for camera calibration. The pre-calibration parameters were compared with self-calibration parameters, by interchanging the heights after performing a total of 63 scenarios using 3, 14, and 50 GCPs and 100 CPs.
Knowing the camera calibration parameters, after the BBA process, the positions of the measured geometric features are as well adjusted, so the less the RMSE of the CPs object coordinates, the better the distortion parameters were estimated.
The errors are very large in the case of self-calibration with a minimum number of GCPs. In the case of a double grid nadiral flight, the RMSE is approximately 60 cm, while in the case of a nadiral single grid, the RMSE is 93 cm for the 35 m flight and 1.84 m for 23 m flight. As mentioned in [29][30][31], the errors of a Digital Surface Model (DSM) outside the GCPs area can reach 3 m, so the results support these findings. For all camera blocks, it was found that when the number of control points increases, the accuracy improves, reaching cm-to dm-level positioning accuracy.
The accuracy is improved by 75% when using pre-calibration in the orientation of the double grid nadiral block image with a minimum number of GCPs. When the number of GCPs increase, the accuracy improves by more than 50% (Figure 18).

Discussion
For this study, a 3D calibration test-field was designed and tested, being the largest reported to date. A number of 349 non-coded targets was measured with millimeter precision in a designed local coordinate system. The camera was calibrated at three different flying heights with two configurations for image acquisition, and the tests were performed on nadiral and oblique flights done at three heights, i.e., 23 m, 28 m, and 35 m, the same as for camera calibration. The pre-calibration parameters were compared with self-calibration parameters, by interchanging the heights after performing a total of 63 scenarios using 3, 14, and 50 GCPs and 100 CPs.
Knowing the camera calibration parameters, after the BBA process, the positions of the measured geometric features are as well adjusted, so the less the RMSE of the CPs object coordinates, the better the distortion parameters were estimated.
The errors are very large in the case of self-calibration with a minimum number of GCPs. In the case of a double grid nadiral flight, the RMSE is approximately 60 cm, while in the case of a nadiral single grid, the RMSE is 93 cm for the 35 m flight and 1.84 m for 23 m flight. As mentioned in [29,30], the errors of a Digital Surface Model (DSM) outside the GCPs area can reach 3 m, so the results support these findings. For all camera blocks, it was found that when the number of control points increases, the accuracy improves, reaching cm-to dm-level positioning accuracy.
The accuracy is improved by 75% when using pre-calibration in the orientation of the double grid nadiral block image with a minimum number of GCPs. When the number of GCPs increase, the accuracy improves by more than 50% (Figure 18). When using pre-calibration parameters in the orientation of a single grid nadiral block image with a minimum number of GCPs, the accuracy is improved by 95% for the 23 m flight and 90% for the 35 m flight. So, using SfM technique to orient the images and only 3 GCPs in the BBA process, the When using pre-calibration parameters in the orientation of a single grid nadiral block image with a minimum number of GCPs, the accuracy is improved by 95% for the 23 m flight and 90% for the 35 m flight. So, using SfM technique to orient the images and only 3 GCPs in the BBA process, the distortion parameters are poorly-modelled, particularly if the camera network is comprised of traditional nadiral images in strips and blocks [10]. When the optimum and a large number of GCPs are used in the process of BBA of the 23 m flight, the self-calibration led to slightly better results than pre-calibration in the range of 1-5 mm. When the optimum of GCPs is used in the case of 35 m flight, the pre-calibration improved the results by 86%, and the errors have almost the same values when a large number of GCPs i.e., 50, was used. So, using 1 GCP/65 sq. m in the process of BBA, the results of self-calibration and pre-calibration converge.
Even if a large number of GCPs are used for orienting a double grid nadiral flight, pre-calibration improves the results with respect to self-calibration with more than 50% in comparison with single grid, because good coverage of the study area is assured. So, having a pre-calibrated camera along with good image configuration and an optimal number of GCPs, an accuracy of half decimeter can be achieved. The impact of pre-calibration when using a large number of GCPs isn't so notable in the case of single grid, taking into account the side overlap of 40% and the fact that the roof has a uniform texture and it covers approximately 15% of the study area.
Orthogonal roll angles must be used to image the test-field for calibration purpose in order to break the projective coupling between interior and exterior orientation parameters, but, using a strongly 3D object point array and higher convergence angles for the images, it is possible to decouple the interior and exterior orientation parameters as mentioned in [7].
Analyzing the residuals when using the two image configuration, 17 oblique images and 21 oblique and nadiral rolled images around the optical axis, we can see for the two single grid nadiral flights that they are slightly different, being better overall in the case of 17 oblique images. Looking at the residuals obtained for the double grid nadiral flight with 17 images camera pre-calibration, they are better in percentage by approximately 35% with respect to self-calibration when using only 3 GCP in the BBA process and 3 ÷ 10% when using the optimum and a large number of GCPs. So, the exclusion of nadiral images from the image blocks used for the calibration process is possible, in agreement with the findings of [7].
In the case of the oblique flight done for an old monastery 2 years before calibration, the accuracy when using the pre-calibration parameters is improved by approximately 60%. The time variation of the interior orientation parameters had no impact on the accuracy of the UAS project even after this long period of time. This can be explained, as Cramer et al. [23] found that the DJI Phantom low-cost UAS system is already close to the concept of the metric camera design after performing a short distance test-field calibration in different epochs.
The smallest residuals computed for the 23 m nadiral flight were obtained using the 23 m pre-calibration height, so the same height for pre-calibration as for the UAS reconstruction project. Analyzing the residuals listed in Table 11 for the 28 m nadiral flight, we can see that better results were obtained when using 35 m pre-calibration height. In addition, analyzing the residuals calculated for the 35 m nadiral flight, the smaller values can be found for the 35 m pre-calibration height. We can conclude that the same height for pre-calibration must be used as for the UAS reconstruction project.

Conclusions
Our study assessed the accuracy of image block orientation through BBA process with a focus on the influence of pre-calibrated camera. For this purpose, the calibration parameters for the digital camera integrated into the "DJI Phantom 3 Standard" UAS platform calculated using the proposed test-field calibration and the "3DF Zephyr Pro" software for three different heights, namely 23 m, 28 m, and 35 m, and for two image configurations only with oblique images and with a combination of nadiral and oblique images, have been tested.
For all tests done, in the case of the double grid nadiral flight, the parameters calculated with the proposed 3D field improved the results by more than 50% when using the optimum and a large number of GCPs, and by 75% to 95% when using a minimum of 3 GCP in all other analyzed cases.
If the UAS reconstruction project involves an area with a very rugged terrain, placing the GCPs and check points and measuring their position is a very time-consuming step of the survey, and may be impossible. Since UAS are built for this purpose, namely for imaging areas that are not easily accessible, such situations are not uncommon. Using the pre-calibration parameters, a minimum number of 3 GCP can be used for BBA for obtaining a decimeter precision.
The self-calibration can be used when a number of well distributed GCPs is available and the network geometry is favorable, but the advantages of using test-field calibration are that the accuracy is significantly improved and the number of GCPs is considerable reduced. This research has investigated the relationship between height for calibration process and the UAS reconstruction project.
It has been shown that the interior orientation parameters calculated for a digital camera integrated in an UAS platform using a test-field calibration can be transferred with success to other UAS reconstruction projects. Also, by testing the interior orientation parameters obtained by using a 3D calibration test-field on UAS projects done before camera calibration (monastery) respectively after camera calibration (nadiral flights), we showed that the UAS cameras can be pre-or post-calibrated.
Further research can be done by performing camera calibration using the OrienTal software developed by Department of Geodesy and Geoinformation from TU WIEN and other commercial software such as Agisoft and Australis.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A