Determining the Suitable Number of Ground Control Points for UAS Images Georeferencing by Varying Number and Spatial Distribution

: Currently, products that are obtained by Unmanned Aerial Systems (UAS) image processing based on structure-from-motion photogrammetry (SfM) are being investigated for use in high precision projects. Independent of the georeferencing process being done directly or indirectly, Ground Control Points (GCPs) are needed to increase the accuracy of the obtained products. A minimum of three GCPs is required to bring the results into a desired coordinate system through the indirect georeferencing process, but it is well known that increasing the number of GCPs will lead to a higher accuracy of the ﬁnal results. The aim of this study is to ﬁnd the suitable number of GCPs to derive high precision results and what is the e ﬀ ect of GCPs systematic or stratiﬁed random distribution on the accuracy of the georeferencing process and the ﬁnal products, respectively. The case study involves an urban area of about 1 ha that was photographed with a low-cost UAS, namely, the DJI Phantom 3 Standard, at 28 m above ground. The camera was oriented in a nadiral position and 300 points were measured using a total station in a local coordinate system. The UAS images were processed using the 3DF Zephyr software performing a full BBA with a variable number of GCPs i.e., from four up to 150, while the number and the spatial location of check points (ChPs) was kept constant i.e., 150 for each independent distribution. In addition, the systematic and stratiﬁed random distribution of GCPs and ChPs spatial positions was analysed. Furthermore, the point clouds and the mesh surfaces that were automatically derived were compared with a terrestrial laser scanner (TLS) point cloud while also considering three test areas: two inside the area deﬁned by GCPs and one outside the area. The results expressed a clear overview of the number of GCPs needed for the indirect georeferencing process with minimum inﬂuence on the ﬁnal results. The RMSE can be reduced down to 50% when switching from four to 20 GCPs, whereas a higher number of GCPs only slightly improves the results. This presented a series of tests to determine the most suitable number and the spatial distribution of GCPs needed to georeference a block of nadiral UAS images acquired with a low-cost platform, i.e., DJI Phantom 3 Standard over an urban area. The images are of high resolution having a GSD of 1.1 cm and they were taken in a cross-strip ﬂight line conﬁguration with high overlaps (80% / 60%), with their processing being performed by bundle block adjustment with self-calibration. The total RMSE calculated for the 150 ChPs improved with more than 50% when using 20 GCPs with respect to the four and eight GCPs scenarios and up to 80% in the case of 150 GCPs. GCPs 20 GCPs higher two control areas inside the GCPs and ChPs boundary, errors ﬁve times smaller for the 20 GCPs scenario with respect to four GCPs, in the case of both studied distributions. A decreasing tendency of the median absolute deviation (stdDevMAD) is visible for all cases and, the 150 GCPs scenario, the di ﬀ erences are in the range of millimetres compared with 20 GCPs scenario also expressed for the total RMSE calculated for ChPs.


Introduction
Unmanned Aerial Systems (UAS) extend the applications range of aerial imaging to various domains, due to efficient acquisition time, costs, and repeatability.

Materials and Methods
Analysing the results reported in [Error! Reference source not found.], we can see that the differences between the RMSE that were obtained for each number of GCPs by processing the projects with Pix4D software and Zephyr Pro software [Error! Reference source not found.], are similar with a range of a few millimetres up to 3 cm, except for the first three case: using three, four, and five GCPs. Moreover, when using the 3DF Zephyr Pro software, all images were orientated without making any additional steps. So, for the approach presented in this paper, all of the tests have been performed using only the 3DF Zephyr Pro software.
The workflow that was applied for this study area can be structured in the following parts: field work with measurements for the local GNSS network and for determining the coordinates of GCPs and ChPs, the placement of targets for materialisation of GCPs and ChPs, planning and execution of UAS flight and office work with the adjustment of GNSS network, determination of GCPs and ChPs coordinates by weighted indirect adjustment model [Error! Reference source not found.] implemented in MATLAB, choosing and creating different files with the desired GCPs and ChPs for each scenario, processing the UAS images corresponding to each scenario, and quality assessment ( Figure 2). The chosen study area reflects a common situation for UAS projects in urban area, such as including artificial surfaces (roof, façade, parking, roads) and natural surfaces (green area). Often, the natural environment can be considered to be an obstacle or impracticable and the placement of GCPs and ChPs can be challenging or even impossible.

Materials and Methods.
Analysing the results reported in [32], we can see that the differences between the RMSE that were obtained for each number of GCPs by processing the projects with Pix4D software and Zephyr Pro software [29], are similar with a range of a few millimetres up to 3 cm, except for the first three case: using three, four, and five GCPs. Moreover, when using the 3DF Zephyr Pro software, all images were orientated without making any additional steps. So, for the approach presented in this paper, all of the tests have been performed using only the 3DF Zephyr Pro software.
The workflow that was applied for this study area can be structured in the following parts: field work with measurements for the local GNSS network and for determining the coordinates of GCPs and ChPs, the placement of targets for materialisation of GCPs and ChPs, planning and execution of UAS flight and office work with the adjustment of GNSS network, determination of GCPs and ChPs coordinates by weighted indirect adjustment model [33] implemented in MATLAB, choosing and creating different files with the desired GCPs and ChPs for each scenario, processing the UAS images corresponding to each scenario, and quality assessment ( Figure 2).

Measurement of Ground Control Points (GCPs)
For these experiments, the chosen solution was to design a spatial geodetic network that was determined by GNSS technology. The spatial geodetic network was designed with four points (A, B, C, and D), two of which were grounded in the green area behind the Faculty of Hydrotechnical Engineering, Geodesy, and Environmental Engineering, and the other two were located on the roof in the northern part of the building. Moreover, for the planimetric rectangular coordinates (x,y) and

Measurement of Ground Control Points (GCPs)
For these experiments, the chosen solution was to design a spatial geodetic network that was determined by GNSS technology. The spatial geodetic network was designed with four points (A, B, C, and D), two of which were grounded in the green area behind the Faculty of Hydrotechnical Engineering, Geodesy, and Environmental Engineering, and the other two were located on the roof in the northern part of the building. Moreover, for the planimetric rectangular coordinates (x,y) and orthometric altitude (H) a local reference and coordinate system was implemented, called LOCAL-HGIM, as previously described in [13]. There were 300 uniformly distributed points grounded in the study area: 29 points in the green space (reinforced concrete poles to ensure different heights), two points of the GNSS network (concrete), 155 points in the parking lot (metallic bolts), 104 points on the roof (marked with paint), and 10 points on the building façade ( Figure 3).

Measurement of Ground Control Points (GCPs)
For these experiments, the chosen solution was to design a spatial geodetic network that was determined by GNSS technology. The spatial geodetic network was designed with four points (A, B, C, and D), two of which were grounded in the green area behind the Faculty of Hydrotechnical Engineering, Geodesy, and Environmental Engineering, and the other two were located on the roof in the northern part of the building. Moreover, for the planimetric rectangular coordinates (x,y) and orthometric altitude (H) a local reference and coordinate system was implemented, called LOCAL-HGIM, as previously described in [Error! Reference source not found.]. There were 300 uniformly distributed points grounded in the study area: 29 points in the green space (reinforced concrete poles to ensure different heights), two points of the GNSS network (concrete), 155 points in the parking lot (metallic bolts), 104 points on the roof (marked with paint), and 10 points on the building façade ( Figure 3). Each new point was double measured using a mini prism from the ends of a GNSS base. The weighted indirect adjustment model was applied for each new point for determining the final spatial coordinates of the detail points, obtaining a total accuracy of a few millimetres, as previously described in [13]. The points were marked using plexiglass plates that were drawn with two black and orange triangles, 3 mm thick, 40 cm × 40 cm, with a central hole of 5 mm diameter, to assure their visibility on the UAS images.

Aerial Image and Terrestrial Laser Scanner (TLS) Acquisition
The flight over the study area was executed at 28 m height above ground with a low-cost UAS DJI Phantom 3 Standard at the beginning of November. This platform has an integrated digital camera with 3.61 mm focal length and a pixel size of 1.624 µm, with the images having a resolution of 4000×3000 pixels.
The flight planning was made with Pix4D Capture software, choosing an end lap and side lap of 80% and 60%, respectively, with the camera being oriented in nadiral position. For the height of 28 m, the flight was made in cross-strips (Figure 4a), with 122 images being acquired with 1.1 cm GSD, providing a very good overlap over the study area ( Figure 4b). Figure 4c presents the camera positions for each strip and the number of visible cameras corresponding to each GCP and ChP.
The flight over the study area was executed at 28 m height above ground with a low-cost UAS DJI Phantom 3 Standard at the beginning of November. This platform has an integrated digital camera with 3.61 mm focal length and a pixel size of 1.624 μm, with the images having a resolution of 4000×3000 pixels.
The flight planning was made with Pix4D Capture software, choosing an end lap and side lap of 80% and 60%, respectively, with the camera being oriented in nadiral position. For the height of 28 m, the flight was made in cross-strips (Figure 4a), with 122 images being acquired with 1.1 cm GSD, providing a very good overlap over the study area (Figure 4b). Figure 4c presents the camera positions for each strip and the number of visible cameras corresponding to each GCP and ChP.  For all of the studied scenarios, the automatically derived products, point clouds, and mesh surfaces were evaluated based on reference data that were acquired by terrestrial laser scanning technology (TLS).
The measurements for the TLS point cloud acquisition were performed using the Maptek I-Site 8820 terrestrial laser scanner, from four different station points in the middle of August. Three of these points A, C, and D are part of the GNSS network, and GNSS technology determined the last one, point E ( Figure 5). For all of the studied scenarios, the automatically derived products, point clouds, and mesh surfaces were evaluated based on reference data that were acquired by terrestrial laser scanning technology (TLS).
The measurements for the TLS point cloud acquisition were performed using the Maptek I-Site 8820 terrestrial laser scanner, from four different station points in the middle of August. Three of these points A, C, and D are part of the GNSS network, and GNSS technology determined the last one, point E ( Figure 5). Artificial spheres have been used in order to georeference the TLS point cloud, which were levelled and centred with the help of tripods and spirit level bubbles ( Figure 6).

TLS Point Cloud Processing
The TLS point clouds processing was done using two different methods. For the point clouds that were taken from the A, C, and D station points the direct georeferencing method was used, and for the point cloud taken from E station point, situated in the green space behind the building, the indirect georeferencing method was used. In both cases, the coordinates of the sphere's centres were determined by fitting the TLS points measured on their surfaces with the mathematical shape while using the least square method implemented into CloudCompare software [Error! Reference source not found.]. The spatial geodetic coordinates (B, L, HE) determined by GNSS technology for the spheres centres and TLS station points have been transformed into the local coordinate system that was designed for this research using a MATLAB script. The direct georeferencing process was performed using the I-Site Studio software, which was specifically designed for Maptek scanners.
The Gauss-Helmert transformation was used, while considering three common points materialised with artificial spheres in order to bring in coincidence the E point cloud with the A+C+D point cloud defined in the local coordinate system (Figure 7). Artificial spheres have been used in order to georeference the TLS point cloud, which were levelled and centred with the help of tripods and spirit level bubbles ( Figure 6). Artificial spheres have been used in order to georeference the TLS point cloud, which were levelled and centred with the help of tripods and spirit level bubbles ( Figure 6).

TLS Point Cloud Processing
The TLS point clouds processing was done using two different methods. For the point clouds that were taken from the A, C, and D station points the direct georeferencing method was used, and for the point cloud taken from E station point, situated in the green space behind the building, the indirect georeferencing method was used. In both cases, the coordinates of the sphere's centres were determined by fitting the TLS points measured on their surfaces with the mathematical shape while using the least square method implemented into CloudCompare software [Error! Reference source not found.]. The spatial geodetic coordinates (B, L, HE) determined by GNSS technology for the spheres centres and TLS station points have been transformed into the local coordinate system that was designed for this research using a MATLAB script. The direct georeferencing process was performed using the I-Site Studio software, which was specifically designed for Maptek scanners.
The Gauss-Helmert transformation was used, while considering three common points materialised with artificial spheres in order to bring in coincidence the E point cloud with the A+C+D point cloud defined in the local coordinate system (Figure 7).

TLS Point Cloud Processing
The TLS point clouds processing was done using two different methods. For the point clouds that were taken from the A, C, and D station points the direct georeferencing method was used, and for the point cloud taken from E station point, situated in the green space behind the building, the indirect georeferencing method was used. In both cases, the coordinates of the sphere's centres were determined by fitting the TLS points measured on their surfaces with the mathematical shape while using the least square method implemented into CloudCompare software [34]. The spatial geodetic coordinates (B, L, HE) determined by GNSS technology for the spheres centres and TLS station points have been transformed into the local coordinate system that was designed for this research using a MATLAB script. The direct georeferencing process was performed using the I-Site Studio software, which was specifically designed for Maptek scanners.
The Gauss-Helmert transformation was used, while considering three common points materialised with artificial spheres in order to bring in coincidence the E point cloud with the A+C+D point cloud defined in the local coordinate system (Figure 7). For the evaluation purpose of the registration process, an overlapping area between the above two point clouds was compared using the Hausdorff distance that is implemented into CloudCompare software. A histogram of the Hausdorff distances distribution was calculated ( Figure 8b) after computing the Hausdorff distances between each point from E point cloud and the mesh surface created based on A+C+D TLS points measured on the parking lot surface (Figure 8a). It can be noticed that most of the calculated distances are in range of -1 mm to -12 mm, the standard For the evaluation purpose of the registration process, an overlapping area between the above two point clouds was compared using the Hausdorff distance that is implemented into CloudCompare software. A histogram of the Hausdorff distances distribution was calculated (Figure 8b) after computing the Hausdorff distances between each point from E point cloud and the mesh surface created based on A+C+D TLS points measured on the parking lot surface (Figure 8a). It can be noticed that most of the calculated distances are in range of −1 mm to −12 mm, the standard deviation being ± 3 mm by fitting a Gauss distribution, according to the histogram of the Hausdorff distances. For the evaluation purpose of the registration process, an overlapping area between the above two point clouds was compared using the Hausdorff distance that is implemented into CloudCompare software. A histogram of the Hausdorff distances distribution was calculated ( Figure 8b) after computing the Hausdorff distances between each point from E point cloud and the mesh surface created based on A+C+D TLS points measured on the parking lot surface (Figure 8a). It can be noticed that most of the calculated distances are in range of -1 mm to -12 mm, the standard deviation being ± 3 mm by fitting a Gauss distribution, according to the histogram of the Hausdorff distances.

Establishing the GCPs and ChPs Distribution
For UAS image processing, a total of 300 points has been chosen in two different spatial distributions: systematic distribution and stratified random distribution, with half of them serving as ChPs.  For the evaluation purpose of the registration process, an overlapping area between the above two point clouds was compared using the Hausdorff distance that is implemented into CloudCompare software. A histogram of the Hausdorff distances distribution was calculated ( Figure 8b) after computing the Hausdorff distances between each point from E point cloud and the mesh surface created based on A+C+D TLS points measured on the parking lot surface (Figure 8a). It can be noticed that most of the calculated distances are in range of -1 mm to -12 mm, the standard deviation being ± 3 mm by fitting a Gauss distribution, according to the histogram of the Hausdorff distances.

Establishing the GCPs and ChPs Distribution
For UAS image processing, a total of 300 points has been chosen in two different spatial distributions: systematic distribution and stratified random distribution, with half of them serving as ChPs.

Establishing the GCPs and ChPs Distribution
For UAS image processing, a total of 300 points has been chosen in two different spatial distributions: systematic distribution and stratified random distribution, with half of them serving as ChPs.
The systematic and stratified random distribution applies to both GCPs and ChPs; one of the conditions imposed was that one point could not be simultaneously a GCP and a ChP.
Regarding the structure arrangement for the stratified random distribution of GCPs and ChPs, their placement is random in each surface type, in order to have high statistical precision over the study area, being automatically obtained. An advantage of this sampling scheme is the avoidance of clustering while covering the entire area and its characteristics.
From the opposite perspective, the systematic distribution can be defined in this case as a methodical approach, where the points, both GCPs and ChPs, constitute a pattern. Therefore, two grids were designed to be a guideline to manually classify the GCPs and ChPs since a perfect regular grid is not possible to achieve in common situations for UAS projects, especially in urban areas, where artificial (roofs, façades, parking lots, roads) and natural surfaces can be found ( Figure 10).
For each approach, nine scenarios were tested where the number of GCPs varied i.e., 4,8,20,25,50,75,100,125, and 150 while the number of ChPs was constantly 150. Both approaches, systematic and stratified random, had one ChPs distribution that was applied for all nine specific scenarios. Thus, the ChPs for the systematic distribution were different from the ChPs for the stratified random distribution.
The first scenario, including four GCPs, considered four different scenarios, which are detailed in the section below. The second scenario with eight GCPs was selected to highlight the thumb rule that when high planimetric accuracy is to be obtained the points have to be around the edge of the area, whereas, for high altimetric accuracy, the points should be placed in the middle of the area. It was found that it is necessary to place GCPs around the edge of the study area to minimize planimetry errors and create a stratified distribution inside with a density of around 1.7 GCP/ha to minimize the altimetry errors [28]. The third scenario with 20 GCPs includes points on the edges and in the middle of the area and its placement takes the first grid with the rule of having at least one point inside the grid cell into account. With the increasing of GCP number, both of the grids offered guidelines for classifying points into GCPs and ChPs.
While considering the fact that ChPs distributions is set as constant along the entire processing, both resulted distributions, systematic, and stratified random, have a third of the ChPs in common while all 150 ChPs have good coverage and representativeness for all types of surfaces (roof, façade, parking lot, green area).
their placement is random in each surface type, in order to have high statistical precision over the study area, being automatically obtained. An advantage of this sampling scheme is the avoidance of clustering while covering the entire area and its characteristics.
From the opposite perspective, the systematic distribution can be defined in this case as a methodical approach, where the points, both GCPs and ChPs, constitute a pattern. Therefore, two grids were designed to be a guideline to manually classify the GCPs and ChPs since a perfect regular grid is not possible to achieve in common situations for UAS projects, especially in urban areas, where artificial (roofs, façades, parking lots, roads) and natural surfaces can be found ( Figure 10).
For each approach, nine scenarios were tested where the number of GCPs varied i.e. 4, 8, 20, 25, 50, 75, 100, 125, and 150 while the number of ChPs was constantly 150. Both approaches, systematic and stratified random, had one ChPs distribution that was applied for all nine specific scenarios. Thus, the ChPs for the systematic distribution were different from the ChPs for the stratified random distribution.
The first scenario, including four GCPs, considered four different scenarios, which are detailed in the section below. The second scenario with eight GCPs was selected to highlight the thumb rule that when high planimetric accuracy is to be obtained the points have to be around the edge of the area, whereas, for high altimetric accuracy, the points should be placed in the middle of the area. It was found that it is necessary to place GCPs around the edge of the study area to minimize planimetry errors and create a stratified distribution inside with a density of around 1.7 GCP/ha to minimize the altimetry errors [Error! Reference source not found.]. The third scenario with 20 GCPs includes points on the edges and in the middle of the area and its placement takes the first grid with the rule of having at least one point inside the grid cell into account. With the increasing of GCP number, both of the grids offered guidelines for classifying points into GCPs and ChPs.
While considering the fact that ChPs distributions is set as constant along the entire processing, both resulted distributions, systematic, and stratified random, have a third of the ChPs in common while all 150 ChPs have good coverage and representativeness for all types of surfaces (roof, façade, parking lot, green area).

UAS Image Processing
For the approach that was studied in this paper, the 3DF Zephyr Pro software was used for processing the UAS images acquired at 28 m height.
For UAS image processing, an initial project was created without specifying the type of camera and the calibration parameters. After importing the images, the software makes an approximation for the interior and exterior orientation parameters for each camera position based on the EXIF information. All 300 GCPs and ChPs were manually measured on each oriented image that they appeared on (minimum 3).
Based on this one, new projects were created while considering each time the same ChPs, while the number of GCPs varied: 4, 8, 20, 25, 50, 75, 100, 125, and 150, each time importing the same image measurements for the GCPs. In order to bring the results into the LOCAL-HGIM coordinate system, the bundle block adjustment (BBA) process was performed, including only the GCPs as constraints (the "Constraint" and "Control" boxes were checked), 150 points serving as check points (only the "Control" option was selected).
The confidence weight for the constraints was set by default at 50%, since increasing the GCPs degree of confidence leads to less confidence of the automatically generated point clouds [29]. After running the bundle adjustment process, the interior and exterior orientation parameters have been adjusted. It has to be specified that in order to adjust the radial and tangential distortions parameters, the boxes "adjust radial distortion" and "adjust tangential distortion" should be checked prior running the project.
For the accuracy assessment of the georeferencing process, the GCPs and ChPs coordinates were compared with the ones that were determined with high precision by the GNSS and total station technology the errors being expressed both, in pixels and in meters. Furthermore, the role of GCPs is to constrain the block of images, the smallest errors being concentrated in their proximity. Thus, to obtain a global perspective over the errors the root mean square error (RMSE) was calculated for the points that were not included in the BBA process, namely the ChPs. For all scenarios, a dense point cloud, together with a mesh surface, were automatically generated.

Results
Four different scenarios were tested when only using four GCPs in the BBA ( Figure 11). The study area has been divided in four and the corners have been marked with a rectangle. In each of the four corners, one GCP has been selected with respect to the interior of the rectangular boundary and with respect to the exterior part of boundary. Subsequently, 10 different combinations have been tested for the systematic distribution and, respectively, for the stratified random distribution. The area covered by the 300 GCPs and ChPs was divided in four and the corners were defined by a rectangular boundary, as mentioned above. In each of the four areas, one GCP was selected while taking into consideration two spatial distribution: in the exterior part of the boundary and in the interior part, not too far out from the corners. Table 1 lists the RMSE calculated for the 150 ChPs when using only four GCPs in the BBA process.  Figure 12 graphically represents the errors along the three axes for each scenario. Analysing the The area covered by the 300 GCPs and ChPs was divided in four and the corners were defined by a rectangular boundary, as mentioned above. In each of the four areas, one GCP was selected while taking into consideration two spatial distribution: in the exterior part of the boundary and in the interior part, not too far out from the corners. Table 1 lists the RMSE calculated for the 150 ChPs when using only four GCPs in the BBA process.  Figure 12 graphically represents the errors along the three axes for each scenario. Analysing the errors, the ChPs distribution has an influence on the residuals obtained, decreasing in both cases: for systematic distribution from 40.2 cm to 31.7 cm and for stratified random distribution from 38.2 cm to 29.6 cm. Nevertheless, the constant difference for RMSE T of 8.5 cm between interior and exterior can be highlighted, for both cases, systematic and stratified random distribution, the differences in meters between the targets selected for GCPs exterior and interior distribution being approximately 6 m. The most suitable distribution of the four GCPs is obtained when the targets are located in the interior of the four corners. It is noticeable that the vertical errors are six times larger than planimetric errors, directly influencing the total RMSE, the results being in accordance with [Error! Reference source not found.]. When comparing with the research that was conducted by Tomaštík et al. [Error! Reference source not found.], a few parameters were in common with the present study: low cost UAS system, image resolution, side-lap and end-lap, and size of the study area.
For the projects in which 8, 20, 25, 50, 75, 100, 125, and 150 GCPs, respectively, where used in the process of BBA, the errors along each axis have been computed: RMSEX, RMSEY, RMSEZ, along with the planimetric error, RMSEXY, and the total error RMSET. Table 2 lists the RMSE calculated for the 150 ChPs when using a systematic distribution of the GCPs and the ChPs and Table 3 lists the RMSE for the stratified random distribution scenario.   It is noticeable that the vertical errors are six times larger than planimetric errors, directly influencing the total RMSE, the results being in accordance with [26]. When comparing with the research that was conducted by Tomaštík et al. [26], a few parameters were in common with the present study: low cost UAS system, image resolution, side-lap and end-lap, and size of the study area.
For the projects in which 8, 20, 25, 50, 75, 100, 125, and 150 GCPs, respectively, where used in the process of BBA, the errors along each axis have been computed: RMSE X , RMSE Y , RMSE Z , along with the planimetric error, RMSE XY , and the total error RMSE T . Table 2 lists the RMSE calculated for the 150 ChPs when using a systematic distribution of the GCPs and the ChPs and Table 3 lists the RMSE for the stratified random distribution scenario.  The planimetric values (RMSE X,Y ) ranged between 2.0-4.7 cm for systematic and between 2.0-4.4 cm for the stratified random distribution, with the highest error corresponding to the eight GCPs scenario. It is mentionable that the differences between the two distributions are approximately a few millimetres, with a maximum of 3 mm. The vertical errors (RMSE Z ) are varying for the systematic distribution between 3.7-29.1 cm, while for the stratified random one is between 3.6-19.1 cm, the highest error corresponding to the eight GCPs scenario.
Analysing the two distributions, the difference between the RMSE T is 1 dm. The values are lower for the stratified random distribution, because the points are located on both the periphery and in the centre of the area. Moreover, while for the systematic distribution, the difference between four GCPs and eight GCPs scenarios is insignificant, only 2.1 cm, for the stratified random distribution is 1 dm.
The error tendency along all the tested scenarios along with the obtained result for each axis are presented in Figure 13 for both distribution methods: systematic (a) and stratified random (b). scenario. It is mentionable that the differences between the two distributions are approximately a few millimetres, with a maximum of 3 mm. The vertical errors (RMSEZ) are varying for the systematic distribution between 3.7-29.1 cm, while for the stratified random one is between 3.6-19.1 cm, the highest error corresponding to the eight GCPs scenario. Analysing the two distributions, the difference between the RMSET is 1 dm. The values are lower for the stratified random distribution, because the points are located on both the periphery and in the centre of the area. Moreover, while for the systematic distribution, the difference between four GCPs and eight GCPs scenarios is insignificant, only 2.1 cm, for the stratified random distribution is 1 dm.
The error tendency along all the tested scenarios along with the obtained result for each axis are presented in Figure 13 for both distribution methods: systematic (a) and stratified random (b). The graphic representations (Figure 13a and 13b) highlight a deep drop of the RMSET between eight GCPs and 20 GCPs cases, being three times lower for systematic and two times lower for stratified random. In terms of GSD, for 20 GCPs, the accuracy is 3 GSD in planimetry, for both distribution types, and 7.9 GSD in elevation for systematic distribution and 7.4 GSD in elevation for the stratified random distribution. Increasing the number of GCPs with 5, the accuracy is the same for the systematic distribution and it improves with only 6 mm for stratified random. Adding 25 GCPs, the accuracy improves with approximately 2 cm for both distributions. The following scenarios, 75GCPs, 100 GCPs, 125 GCPs, 150 GCPs, follow a decreasing tendency of errors, but at a lower rate in terms of millimetres. By having a great amount of GCPs (i.e. 50 GCPs), the planimetric errors are almost the same while the vertical errors are lightly reduced. Figure 14 individually highlights the spatial representation of the total root-mean-square-error of the check points for each scenario. It is noticeable the placement of both GCPs and ChPs for the systematic distribution (Figure 14a) and the stratified random one (Figure 14b) offers a more detailed analysis. The graphic representations (Figure 13a,b) highlight a deep drop of the RMSE T between eight GCPs and 20 GCPs cases, being three times lower for systematic and two times lower for stratified random. In terms of GSD, for 20 GCPs, the accuracy is 3 GSD in planimetry, for both distribution types, and 7.9 GSD in elevation for systematic distribution and 7.4 GSD in elevation for the stratified random distribution. Increasing the number of GCPs with 5, the accuracy is the same for the systematic distribution and it improves with only 6 mm for stratified random. Adding 25 GCPs, the accuracy improves with approximately 2 cm for both distributions. The following scenarios, 75GCPs, 100 GCPs, 125 GCPs, 150 GCPs, follow a decreasing tendency of errors, but at a lower rate in terms of millimetres. By having a great amount of GCPs (i.e., 50 GCPs), the planimetric errors are almost the same while the vertical errors are lightly reduced. Figure 14 individually highlights the spatial representation of the total root-mean-square-error of the check points for each scenario. It is noticeable the placement of both GCPs and ChPs for the systematic distribution ( Figure 14a) and the stratified random one (Figure 14b) offers a more detailed analysis.
Analysing the Figure 14, it can be seen that, for the ChPs situated on the roof edge and on the building façade, the errors are still larger when compared with the remaining ones, when using 20 GCPs or more. This was also the case of the Rangel et al. research [22], where the errors tend to be clustered in areas with abrupt slope changes. They also mentioned that altimetric errors are located in areas with little texture and radiometric uniformity. This fact also represents a particularity of the current study, with the roof area having a uniform texture (light grey) extended on 0.13 ha.
In addition, the errors were classified according to their location in the four main segments of the study area, namely: roof, façade, parking lot, and green area ( Figure 15). Due to considerably large amount of generated data, based on the RMSE T , the descriptive statistics were also listed in Figure 15, taking all cases from four GCPs to 150 GCPs into the analysis. A decreasing central tendency of RMSE T is clearly presented, being directly proportional with the variability, which, for all cases, shows a similar pattern, as described by Figure 15. Analysing the Figure 14, it can be seen that, for the ChPs situated on the roof edge and on the building façade, the errors are still larger when compared with the remaining ones, when using 20 GCPs or more. This was also the case of the Rangel et al. research [Error! Reference source not found.], where the errors tend to be clustered in areas with abrupt slope changes. They also mentioned that altimetric errors are located in areas with little texture and radiometric uniformity. This fact also represents a particularity of the current study, with the roof area having a uniform texture (light grey) extended on 0.13 ha.
In addition, the errors were classified according to their location in the four main segments of the study area, namely: roof, façade, parking lot, and green area ( Figure 15). Due to considerably large amount of generated data, based on the RMSET, the descriptive statistics were also listed in Figure 15, taking all cases from four GCPs to 150 GCPs into the analysis. Furthermore, the ChPs with high errors are located in the entire area, in all types of surfaces, natural or artificial, low and high altitude. In addition, the ChPs determined in the parking lot, which is located in the middle of the area, are influenced by the only ground control point marked on the right edge of it. The errors are gradually increasing with the distance, reaching values around 70 cm.
Special attention can be drawn to the eight GCPs scenarios where the behaviour of the ChPs error is directly influenced by the GCPs location. Regarding the systematic distribution, by adding four peripheric GCPs to the four GCPs scenario, the errors are decreasing below 10 cm in the middle of the roof and around 20 cm on the edges. The error discrepancy is reduced by analysing the spatial variability for the remaining three areas (façade, parking lot, green area). The mean is 0.25 m and the median is 0.21 m. tendency of RMSET is clearly presented, being directly proportional with the variability, which, for all cases, shows a similar pattern, as described by Figure 15. Furthermore, the ChPs with high errors are located in the entire area, in all types of surfaces, natural or artificial, low and high altitude. In addition, the ChPs determined in the parking lot, which is located in the middle of the area, are influenced by the only ground control point marked The stratified random distribution includes GCPs in the interior of the area, resulting in the errors being largely reduced, with exception being made by one ChP located in the farthest North with a value of 55.7 cm. Overall, the highest errors are only on the façade and in the parking lot having less than 30 cm each. The mean and the median, are equally, both having a value of 0.17 m.
In the cases of 20 GCPs, the greatest errors encountered on the façade and on the area's edges do not exceed 30 cm and the general decreasing tendency of the errors is highlighted. The average values became smaller than 10 centimetres, expressing for the both distributions very close ranges, having a mean of 0.074 m and a median of 0.062 m in the systematic case, while for the stratified random one, the mean and the median are 0.071 m, respectively, 0.058 m.
The interquartile range is almost the same for all scenarios, starting with 25 GCPs. Evaluating the computed outliers, in the case of systematic distribution, the total amount of values that are included in this category is no more than 10%. In contrast, the stratified random distribution does not exceed 4%.
With regards to the quality assessment of the UAV obtained point clouds and meshes, the results of both types of GCPs-ChPs distributions were compared with the reference data that were obtained by TLS technology while using the "Distance-Cloud/Mesh Dist" function implemented in CloudCompare software being calculated as the Hausdorff distances between each point and the corresponding triangle surface. For the comparisons of the mesh surfaces, a point cloud was generated for the mesh to be compared by sampling each surface by one-million points. Table 4 summarizes the standard deviation of the distances calculated between each point of the point cloud and each sampled point on the meshes automatically generated after the BBA process using four GCPs, 20 GCPs, and 150 GCPs, and the reference TLS mesh. It is noticeable to mention that before computing the standard deviation the point cloud and the sampled mesh were manually filtered based on the histograms of the calculated distances removing outliers. The errors are divided by half when using 20 GCPs for both scenarios, systematic and stratified random, as highlighted by Table 4. For the 150 GCPs case, the standard deviation is decreasing in the range of millimetres for the stratified random distribution and for the systematic distribution in the range of centimetres.

Discussion
For this study, different tests were performed for determining the suitable number and the spatial distribution of GCPs that are needed for the indirect georeferencing process of UAS images with the minimum impact on the accuracy of the final products, such as point clouds and mesh surfaces. One nadiral flight has been performed at 28 m height, over an area of about 1 ha while using a low cost UAS. The UAS image processing has been performed using the 3DF Zephyr Pro software, with the accuracy assessment of the indirect georeferencing process relying on computing the RMSE based on ChPs. Moreover, the accuracy of point clouds and mesh surfaces automatically generated based on UAS images was tested while using the Hausdorff distances [34,35] between each individual point and sampled points on each triangle surface, respectively, and the reference surface that is represented by a TLS mesh surface.
The research highlights that when using a systematic or a stratified random distribution of GCPs, the errors decrease significantly in the case of 20 GCPs with respect to the errors that were obtained for the eight GCPs scenario. The total RMSE decreases with only 2 cm in the case of 50 GCPs with respect to the RMSE obtained for 20 GCPs, reaching half of the error when using a very large number of GCPs (the maximum considered for this study), i.e., 4.2 cm for the systematic distribution and 4.1 cm for the stratified random distribution. The planimetric errors are slightly different in the range of 1 cm when using 20 GCPs and 150 GCPs, respectively. The stratified random distribution of GCPs and ChPs gave the best results, with the total errors being smaller at the decimetre level for the eight GCPs scenario, centimetre level for the 20, 25, and 50 GCPs, and millimetre level for the rest of scenarios.
The results are in accordance with the findings of Rangel et al. [22], where 18 to 20 GCPs with uniform distribution on the image block is the optimum number since introducing more GCPs there is an insensitive improvement inaccuracy, even if the accuracy assessment was made for DSM and orthophotos of an open pit mine of 270 ha. The magnitudes of the RMSEz that were found in this research when using 20 GCPs in both systematic and stratified random distribution are almost the same as those found by Tonkin et al. [20], which evaluated the vertical accuracy of DSMs with 530 ChPs with known heights.
For a better data comprehension, three control areas ( Figure 16) have been specifically selected, two inside the surface covered by the GCPs and ChPs locations (Control Area 1 and 2) representing a parking lot and the building roof, respectively, and one outside the surface that covers the parking lot situated on the right side of the building and a car access (Control Area 3). The particularity of these surfaces is that they are man-made structures with a low degree of roughness.
For the point clouds and sampled meshes that correspond to the first two control areas, the median absolute deviation (stdDevMAD), has been calculated to evaluate the accuracy for both distributions: systematic (Table 5) and stratified random (Table 6) while using OPALS software [36,37], because gross matching errors can lead to very big errors and they should not be considered in the computation of the standard deviation.
The results are in accordance with the findings of Rangel et al. [Error! Reference source not found.], where 18 to 20 GCPs with uniform distribution on the image block is the optimum number since introducing more GCPs there is an insensitive improvement inaccuracy, even if the accuracy assessment was made for DSM and orthophotos of an open pit mine of 270 ha. The magnitudes of the RMSEz that were found in this research when using 20 GCPs in both systematic and stratified random distribution are almost the same as those found by Tonkin et. al [Error! Reference source not found.], which evaluated the vertical accuracy of DSMs with 530 ChPs with known heights.
For a better data comprehension, three control areas ( Figure 16) have been specifically selected, two inside the surface covered by the GCPs and ChPs locations (Control Area 1 and 2) representing a parking lot and the building roof, respectively, and one outside the surface that covers the parking lot situated on the right side of the building and a car access (Control Area 3). The particularity of these surfaces is that they are man-made structures with a low degree of roughness. For the point clouds and sampled meshes that correspond to the first two control areas, the median absolute deviation (stdDevMAD), has been calculated to evaluate the accuracy for both distributions: systematic (Table 5) and stratified random ( Table 6) while using OPALS software [Error! Reference source not found.,Error! Reference source not found.], because gross matching Figure 16. The representative surfaces used for quality assessment located inside (Control area 1: parking lot, Control area 2: roof) and outside (Control area 3: parking lot) the GCPs-ChPs area. According to [36], the median absolute deviation (stdDevMAD) is a robust standard deviation estimator computed from the median of absolute deviations from the median of all values (Equation (1)).
Being a robust statistic, the data quantification is more resilient due to the fact that large deviations as outliers do not have greater influence on the final result.
As highlighted by Tables 5 and 6, the errors are approximately five times smaller when using 20 GCPs for both scenarios, systematic and stratified random. For the 150 GCPs case, the median absolute deviation is decreasing in the range of millimetres for the systematic distribution and for the stratified random distribution in range of centimetres only for Control Area 1.
An influential factor of the statistical variables results is also represented by the position of the control surfaces with respect to the area that is covered by GCPs-ChPs. Therefore, in the case of the systematic distribution, for the Control area 2 located at the border of the GCPs-ChPs area, the errors are almost five times smaller for the 20 GCPs scenario with respect to four GCPs, for point clouds and sampled meshes, while for the control area 1 located in the middle of the area the errors are even seven times smaller. For the stratified random distribution, the errors are almost five times smaller for 20 GCPs scenario with respect to four GCPs, for both control areas, but, even so, they are smaller or almost equal for the control area 1 with respect to control area 2 in most cases.
A decreasing tendency of the median absolute deviation (stdDevMAD) for all cases is visible and, regarding the 150 GCPs scenario, the differences are in the range of millimetres when compared with the 20 GCPs scenario also expressed for the residuals of 150 ChPs calculated in the Results section.
The errors corresponding to each point of the point cloud and each sampled point belonging to mesh surfaces obtained for the four GCPs, 20 GCPs, and 150 GCPs scenarios in the case of control areas 1 and 2 after the comparison with the TLS mesh, can be visualized in Appendix A for both distribution. Comparing the point clouds and sampled meshes on control area 3 ( Figure 16) with the reference surface will give extra information regarding the error's propagation. Therefore, Table 7 lists the median absolute deviation (stdDevMAD) for the point clouds and sampled meshes that were obtained by UAS image processing using four GCPs in the BBA process, 20 GCPs and 150 GCPs, respectively. In the case of the systematic distribution, it can be noticed that the errors found in this area are of the same magnitude as those encountered in control areas 1 and 2 when only using four GCPs, decreasing to the decimetre level when using the suitable number, but still being two times larger than those encountered in the GCPs area. It is concluded when globally investigating the errors calculated for this area by comparing the median absolute deviation with the total RMSE that they are two times larger than in the GCPs area when using the suitable and the maximum number of GCPs. From all of the computed scenarios, in the case of both distributions, we found that the suitable number of GCPs to be used in the BBA process without any effect on the accuracy of the final products is 20 GCPs. Moreover, the total RMSE was smaller than in the case of systematic distribution when using the stratified random distribution for the spatial location of GCPs and ChPs.

Conclusions
For the UAS reconstruction project, placing and measuring of GCPs can be very difficult and even impossible in some situations, so knowing the suitable number and the spatial distribution of these points is of primary interest for researchers and practitioners. This article presented a series of tests to determine the most suitable number and the spatial distribution of GCPs needed to georeference a block of nadiral UAS images acquired with a low-cost platform, i.e., DJI Phantom 3 Standard over an urban area. The images are of high resolution having a GSD of 1.1 cm and they were taken in a cross-strip flight line configuration with high overlaps (80%/60%), with their processing being performed by bundle block adjustment with self-calibration.
The total RMSE calculated for the 150 ChPs improved with more than 50% when using 20 GCPs with respect to the four and eight GCPs scenarios and up to 80% in the case of 150 GCPs.
The metric evaluation of automatically generated point clouds and mesh surfaces was completed by comparing the points and the sampled triangles with a reference mesh, namely a TLS mesh computing the median absolute deviation (stdDevMAD) that was implemented into OPALS software, concluding that the errors outside the GCPs area are two times larger than those that are encountered in the GCPs area when using 20 GCPs or a higher number. Regarding the two control areas situated inside the GCPs and ChPs boundary, the errors are almost five times smaller for the 20 GCPs scenario with respect to four GCPs, in the case of both studied distributions. A decreasing tendency of the median absolute deviation (stdDevMAD) is visible for all cases and, regarding the 150 GCPs scenario, the differences are in the range of millimetres compared with 20 GCPs scenario also expressed for the total RMSE calculated for ChPs.
Even if there is a slightly difference between the errors that were obtained for stratified random distribution with respect to systematic distribution, in most cases they are smaller, concluding that the stratified random distribution can be used for placing the GCPs over the study area.
It can be concluded that, in order to obtain high accuracy of the final products, the following guidelines can be considered: -control points in the corners are essential, but should be placed not too far out in the corners of the area of interest; -eight GCPs are better than four GCPs, which is to be expected. However, placing them along the border of the block is not optimal, and interior control points are improving the accuracy significantly; -a stratified random placement of control points offers a similar accuracy and an even better one than a systematic placement; -an increase in the number of control points leads to improved accuracy, which is to be expected. The accuracy converges to two GSD in planimetry and three GSD in elevation; -up to 20 control points the accuracy improves strongly, but for higher number of control points the improvements are only marginal. Placing 20 GCPs over a surface of approximately 4000 m 2 an accuracy of 3 GSD in planimetry and 7 GSD in elevation was obtained. Keeping in mind the characteristics of the study area and the flight design, the suitable number of GCPs needed for obtaining accurate results may be specific for the presented study; and, -extrapolating beyond the area enclosed by control points it leads to lower accuracy, also if a very high number of control points is used.
The results demonstrated a clear overview of the number and distribution of GCPs needed for the indirect georeferencing process with minimum influence on the final results. Future research will concentrate on a larger area, with more buildings, by flying at different flight heights.   Appendix A Appendix A Appendix A Appendix A