A METHOD FOR ESTIMATING THE ACCURACY OF THE TIE POINTS CLOUD: THE CASE STUDY OF PALAZZO DE FALCO

: In the context of cultural heritage documentation and conservation, an accurate 3D model is essential to properly describe an artifact of historical relevance and provide the expert user with verified metric data. In this paper we propose a workflow aimed at estimating errors in the image orientation steps, which is useful for improving the accuracy of the sparse cloud (composed of Tie Points). Some parameters provided as input to the processing software, often overlooked by the user in favour of the default ones, are analysed. Specifically, tests are carried out using the commercial software package Agisoft Metashape and involve the analysis of three parameters: (i) the accuracy of measuring the Ground Control Points (GCPs) provided by the user; (ii) the accuracy of the Tie Points (TPs) computed by the software; and (iii) the accuracy of collimating the GCPs on the image. The accuracy of the GCPs is estimated by implementing the formulas of error propagation in MATLAB, considering the three concurring sources of error: the coordinate measurement with the Total Station (TS), the coordinate transformation from the local system to the UTM/ETRF00 System, and the GNSS measurement for estimating the transformation parameters. The transformation parameters are calculated using the Bursa Wolf method. The collimation accuracy of GCPs on the image is estimated by the reprojection error of each single GCPs on the image plane. The accuracy of the GCPs will be estimated by the standard error per unit weight (SEUW), accurate GCPs being expected to have a SEUW value close to unity. SEUW is also the overall indicator for estimating the accuracy of the TP S cloud because it is influenced, in different ways, by all three parameters analysed. Test Area is a historic building located in the town of Fisciano, a few kilometres from the city of Salerno, known as Palazzo De Falco. The results highlight how these parameters, if not properly considered, can significantly affect the final 3D model.


INTRODUCTION
The use of close-range Photogrammetry in the survey and documentation of Cultural Heritage is now widely documented in literature, as the technique enables the creation of an accurate 3D model, which is essential to better describe a heritage site and provide the end user, usually a technician, with metric data useful for many different purposes (Colomina and Molina, 2014;Jasińska et al., 2023;Linder, 2009;Shen et al., 2023, Themistocleous, 2020. The approaches based on photogrammetric techniques are in some respects considered to be more efficient than active sensor-based techniques, such as Laser Scanning, since they allow for better description and interpretation of linear surface features, especially for the representation of linear structural elements (Alshawabkeh, 2020), as well as providing an accurate cracking picture for concrete or masonry structures (Jahanshahi and Masri, 2013). In contrast to methods using active sensors, photogrammetric techniques require a well-lighted scene, very accurate optical sensors, and highly accurate ground survey campaign to scale and orient the photogrammetric model (Remondino et al., 2017). Close-range photogrammetry and its application in the field of digital surveying have grown considerably in recent years, mainly due to the implementation of highly automated algorithms for image processing and reality-based model restitution (El Ghazouali et al., 2022). These algorithms are able to estimate the interior and exterior orientation parameters of the frames, even if they come from uncalibrated cameras, by operating in this case a self-calibration (Xie et al., 2022). The accuracy of the photogrammetric model is directly dependent on both the accuracy of the Ground Control Points (GCPs) used for georeferencing the images, the quality of the calibration of the optics and the size of the object to be surveyed (Barba et al., 2019a).
In most applications of photogrammetry for cultural heritage documentation, the cameras involved are often commercial, uncalibrated cameras. In this case, unlike applications that require high accuracy, which need the use of a calibrated or even photogrammetric camera, the SfM (Structure from Motion) algorithms being used allow the automated calculation of calibration parameters (self-calibration) of cameras that have not been previously calibrated using only point correspondences between different images (Eltner and Sofia, 2020). SfM algorithms require the combined use of nadiral and oblique images, which is necessary to improve shape definition, surface continuity and better description of sub-vertical walls (Barba et al., 2019b), as well as to improve the self-calibration process (Fraser, 2013). The first step of the process, consisting of the acquisition of the images, which must be in a higher number than strictly necessary to avoid shadow areas in the model and to improve the matching process, and the subsequent step of aligning the images and measuring the GCPs are the most delicate phases of the entire process of producing an accurate model (Sanz-Ablanedo et al., 2018). Poor model accuracy can nullify the high resolution of the data and, consequently, the output of the derived 3D model generated. The sparse point cloud, consisting only of Tie Points (TPs), is the basis for producing the complete 3D model. The analysis and removal of low-quality TPs is advisable because their presence affects the results of the next steps, consisting of the recalculation of orientation parameters, and the creation of the final output, the dense cloud (Barba et al., 2022). This aspect is well known in literature; by applying rigorous workflows and specific modelling strategies, including TPs processing and filtering techniques, it is possible to achieve levels of accuracy in 3D modelling comparable to those that can be expected by applying conventional and standardized ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume X-M-1-2023 29th CIPA Symposium "Documenting, Understanding, Preserving Cultural Heritage: Humanities and Digital Technologies for Shaping the Future", 25-30 June 2023, Florence, Italy metrological survey techniques (Riveiro et al., 2013). If these new methods of calculation and data processing help to increase the applicability and versatility of the photogrammetric technique, at the same time they raise the demand to define an unambiguous, or at least widely accepted, way to assess the quality of the models produced. A prerequisite for achieving this goal, however, is the knowledge and control of all the data involved in Bundle Block Adjustment, in terms of accuracy. To overlook this aspect would foreclose the possibility of properly evaluating any output resulting from the process itself. The aim of our study is to address this specific requirement and to clear up the meaning and role of some parameters involved in commercial algorithms that are sometimes overlooked in some applications.

METHODS AND MATERIALS
The case study to test the proposed methodology is an historic building located in the municipality of Fisciano, a few kilometers from the city of Salerno (Campania Region, Italy), known as Palazzo De Falco. The main facade of the building was surveyed by close-range photogrammetry, while both GCPs and Check Points (CPs) were measured by Total Station (TS) (Figure 1), combined with GNSS receivers used to frame the model in the cartographic reference system. The images were acquired with a DSRL Nikon D800E camera mounting an uncalibrated 24mm focal length lens. A total of 115 images were acquired, of which 76 were acquired with very close nadiral shots and 39 with oblique shots. The average GSD (Ground Sample Distance) results in about 2mm and the average image scale is 1/345.

Interior and Relative Orientation of images
Interior orientation parameters were estimated using selfcalibration algorithms. No additional camera calibration parameters were used. The preliminary step for calculating the relative exterior orientation parameters (alignment phase, camera orientation) is to identify and extract the significant points present on each frame (Key Points, KPs). The next stage is the matching step, in which the software runs the process of correlating the homologous points to identify the TPs and calculate the "sparse" cloud. This is a similar approach to that scale-invariant feature transform (SIFT), but it uses optimized algorithms to get higher alignment quality (structure estimation) (Nex and Remondino, 2014). The parameters set for frame alignment are: Full resolution of images (7,360 x 4,912) to search KPs, imposing a limit of 100,000 per frame; TPs limit equal to 0. TPs will be filtered in the next step of optimization.

Ground Control Points survey
The next step is the georeferencing of the images, requiring the measurement of GCPs. In our application, 22 GCPs were measured with TS Leica TCRA 1102, of stated precision of 0.6 mgon (2") in angle measurement (azimuthal and zenithal) and of 2 mm + 2 ppm on distance without prism. GCPs were chosen corresponding to easily recognizable details on the frames, such as sharp edges, to minimize collimation errors. The location of GCPs on the facade is shown in Figure 1a. The station point of the TS was materialized with a topographic nail at a position facing the facade of the building (Figure 1b). No artificial targets were used because the methodology was developed for expeditious surveys. Targets are very often not used in this case for multiple reasons, including the difficult accessibility of the façades of historic buildings without the use of a forklift that would make the survey no longer expeditious or even because of the restrictions imposed by Superintendence for Architectural Heritage. Three points (100, 200 and 300 in Figure 1b) measured with dual-frequency GNSS receiver in fast-static mode were also materialized so that the survey could be framed into the UTM/ETRF00 coordinate system. To transform the coordinates from the Local System (TS) to the Cartographic system with orthometric heights (Geoid ITALGEO 2005), we applied 7parameter Bursa-Wolf transformation. The coordinates of the GCPs will be associated with errors from three sources: the coordinate system transformation (roto-translation and scaling), the GNSS measurement, and the TS measurement, estimated with the formulas of propagation of uncertainty as described in Section 2.3.1. The implemented computational approach is aimed at associating to each single GCP the accuracy derived not only from the measurement of its coordinates made with TS but also from the measurements of a GNSS network referred to a Global System. In contexts where the survey object is very large, the measurements of the GCPs with respect to the local System are inserted into an external (Global) System via a frame network. The vertices of the frame network can be measured with GNSS receivers and in turn will be characterized by an uncertainty that is a function of the measurement modality. The uncertainties computed and associated with each GCP will, therefore, be used to give different weights to the GCPs in the image georeferencing phase. Not all the 22 measured points on the facade were used as GCPs for georeferencing, but 6 of them were used as CPs to validate the results, taking care to ensure an even distribution of the two groups in the surveyed scene.

Tie Points cloud optimization
In this step, following the interior and relative exterior orientation (section 2.1), the software runs a filtering in order to increase the quality of the TPs cloud. The filtering allows for a more accurate estimation of the interior camera orientation and relative image orientation parameters. The parameters used for filtering the TPs are: (i) reconstruction uncertainty, (ii) projection accuracy, (iii) reprojection error. "Reconstruction Uncertainty" (RU) is a numerical representation of the uncertainty in the position of a TP based on the geometric relationship of the cameras from which that point was projected or triangulated, considering geometry and redundancy. RU can also be thought of as the ratio between the largest and smallest semi-axis of the error ellipse created when triangulating 3D point coordinates between two images. The value of this parameter is dependent on the angle between the homologous radii that generate the position of the TPs, and therefore on the base-to-height (B/H) ratio. All those TPs with a RU greater than 10 were removed. This value corresponds to an angle between the homologous rays of two frames of 23°, smaller would not be acceptable. The second filtering parameter of TPs removes points based on the "Projection Accuracy" which is a measure of the "Mean key point size". Key point size (in pixels) is the standard deviation value (in σ) of the Gaussian blur at the scale at which the key point was found; lower standard deviation values are more precisely located in space. In our test, a threshold of 5 was set, all TPs with associated projection accuracy values greater than 5 will be removed. Lastly, the "Reprojection Error" ε is a geometric error corresponding to the image distance between a projected point and a measured one. It is used to quantify how closely an estimate of a 3D point recreates the point's true projection. All TPs with an associated error projection greater than 0.3 pixels were removed. The step of filtering is an iterative procedure. Removal of lowquality TPs will improve the estimation of the interior and exterior (relative) orientation parameters. Each time the TPs are removed, the camera orientation parameters and relative orientation of the images are updated, the TPs are recalculated, and the associated accuracy parameters are also estimated again.

Analysis of the accuracy parameters of the Tie Points Cloud
In the photogrammetric process, special attention should be paid to some parameters that greatly affect the accuracy of the TPs cloud, these are (i) the accuracy of the GCPs survey; (ii) the accuracy of the GCPs collimation; (iii) the accuracy of the Tie Points. The interior and exterior orientation of the images is a function of the accuracy of the GCPs survey and, of course, the accuracy of the collimation phase of the GCPs themselves on the images. The goodness of the interior and exterior orientation, in turn, affects the accuracy of the TPs cloud. A correct setting of these parameters avoids misleading statistics, while their incorrect estimation generates unrepresentative error statistics, as the interior orientation parameters are very sensitive to these parameters, which are often replaced by the default ones. The estimation of an indicator that is representative of the overall quality of the TPs cloud and that considers these three main parameters is therefore critical and unavoidable. The overall parameter chosen is σ0, the Standard Error per Unit Weight (SEUW). The farther the SEUW is from 1, the less reliable the estimated accuracy of TPs will be.
The following sections will introduce the methodology used to estimate the three parameters that influence the global indicator of the TPs cloud accuracy.

Accuracy of the GCPs survey (i):
The uncertainty associated with the coordinates of the GCPs was estimated in the process of transformation from the local system to the exterior (cartographic) system. To estimate the 7 parameters (rototranslation and scaling) we implemented in MATLAB the Bursa Wolf model, which requires a minimum of three common points (9 equations in 7 unknowns). Given the coordinates of the three points (100, 200 and 300 in Figure 1b where εx, εy, εz are the small rotations around the axes of the system. The estimation of the 7 transformation parameters x = [Δx Δy Δz εx εy εz k] T was carried out using the least squares criterion. For solving the model (3), the Combined Least Squares method (Deakin, 2006) was used, which allows adding the weights to the coordinates of the points of the exterior system {W}. The weights to be assigned to the points are derived directly from the. variance-covariance matrix of the GNSS baselines. The least-squares solution is iterative; we first calculate the approximate values of the 7 roto-translation parameters, then the corrections δx to be given to the approximate values x 0 of the unknown parameters are calculated by solving the normal system. The variance-covariance matrix of the transformation parameters, also containing the errors due to the different precision of the points of the system {W}, is given by: where v is the vector of residuals, f 0 is the vector of known terms, the matrix P contains the weights of the observations of the system {W} and the matrix A contains the partial derivatives: (5) Point's Transformation Uncertainty is estimated according to the relation (5) proposed by (Ren et al., 2015): where: R0 and R' are the rotation matrices contained in equations (1) and (2), is the coordinate variance-covariance matrix in the local system {M} calculated with the Jacobian J, according to the relations: where x, y, z are the GCPs coordinates in the local system {M}, Q0 is the variance-covariance matrix of the measures acquired by the TS, comprising the uncertainty both in distance and in angular measurements and J is the Jacobian matrix. The variance-covariance matrix , in addition to considering the errors due to the 7-parameter transformation, is calculated by also considering the accuracy by considering the accuracy of the points measured in both local {M} and external {W} systems, with TS and GNSS instruments.

Accuracy of the CCPs collimation (ii):
The accuracy in collimating GCPs on the image (expressed in pixels) depends largely on their position on the images and is obviously correlated with the image resolution. In case of a highly accurate collimation, which is more easily achieved if artificial targets are used, this parameter may assume values smaller than 0.5 pixels (which is the default value). In other cases, such as our test site, if natural points are used, the reprojection error on the points, which can be visualized after collimation on the images, can be used to evaluate accuracy. After importing the coordinates of the GCPs with associated Standard Deviations, a first calculation cycle is run in which only a few camera calibration parameters are input: f, cx, cy, k1, k2, k3, p1, p2, while leaving out others (Fraser, 2013). In detail, the parameters considered are the Focal length (f) measured in pixels, the principal point coordinates (cx, cy), i.e., coordinates of lens optical axis intersection with sensor plane in pixels, the three radial distortion coefficients (k1, k2, k3) and the two tangential distortion coefficients (p1, p2). Other parameters were neglected in the first calculation cycle: Affinity and Skew (non-orthogonality) transformation coefficients (B1 and B2), parameters that are worth considering only in case of high values of the reprojection error and the radial distortion coefficient (K4). Input of all parameters may lead to a lower reprojection error, but not reliable because the model is too complex and does not reflect reality.
After the first cycle, the marker accuracy value (in pixels) will be set equal to the mean value of the reprojection error obtained in that calculation cycle.

TS and GNSS survey
The three points (100) (200) and (300) used for the transformation of coordinates from the local system to the UTM/ETRF00 cartographic system were measured in fast-static mode with a GNSS receiver, the occupation time on each point was programmed to be 30 minutes with a log. of 1 second. To calculate the baselines, the AVEL station of the Campania Region permanent station network, located about 10 km from the case study, was used. The area is near the historic downtown, characterized by narrow streets and tall buildings on both sides of the street. This configuration induces multipath problems, which greatly lower the accuracy of the received signal. The calculation of baselines gives RMS in the order of 10 mm in planimetry and 30 mm in height. The three points also measured in the local system with TS, were used to estimate the 7 transformation parameters of the Bursa Wolf model. Table 1 shows the 7 calculated transformation parameters.

Bursa Wolf Parameters Values
Δx (  . This result seems plausible given that the largest contribution to the error comes from the GNSS measurement itself. While the error on the measurement made with TS, given the close distances of the GCPs from the TS station (30 m max) and the technique used to measure the GCPs (polar method), has a slight influence on the overall error. Although it is a known good rule to distribute GCPs (and CPs) uniformly over the object to be surveyed, as their distribution greatly affects the accuracy of the final model, in some cases of architectural surveying it is forced by the position on the façade of details having characteristics of good collimability on both it and the images. Because of this the distribution of GCPs on the façade of Palazzo De Falco is not regular in the two directions, they are actually positioned along the vertical and this affects the results. The distribution of points, which is more compact along the East direction (Figure 1), results in larger SD values in that direction and smaller SD values in the North one.

Accuracy assessment
GCPs, with associated SDs, were used for georeferencing the images. The mean error on the GCPs is 3 mm, ranging between 1 mm and 6 mm. The mean error on CPs is 5 mm, the maximum and minimum error are 10 mm and 1mm, respectively. The TPs cloud resulting after the alignment process was more than 600000 points. After filtering, the number of points decreased to just over 250000, about 60% of TPs were removed.
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume X-M-1-2023 29th CIPA Symposium "Documenting, Understanding, Preserving Cultural Heritage: Humanities and Digital Technologies for Shaping the Future", 25-30 June 2023, Florence, Italy The collimation accuracy of GCPs on the image was estimated by means of the reprojection error of each individual GCPs on the image frame plane. After the filtering process was performed, an average reprojection error of 0.3 pixels was estimated, the maximum error is 0.41 pix on the GCPs and 0.53 pix on the CPs, the minimum error is 0.12 pix on the GCPs and 0.21 pix on the CPs. The accuracy of the TPs was estimated by the SEUW, the accurate TPs will have a SEUW value close to unity. The estimated accuracy value of the TPs is 0.2 pixels. This value produced a SEUW value of 0.98. The variance-covariance matrix associated with each individual TPs was calculated by running a python script executable directly in the Metashape software. The error ellipsoid with k = 3 was derived from it. The major semi-axis was used to study the frequency distribution and to construct a tolerance interval (Natrella, 2013) suitable for the resulting frequency distribution (Figure 3). Since the length of the semi-axis is a quantity defined as positive, we construct a one-sided tolerance interval. This is calculated by imposing a confidence coefficient γ of 0.95 and a population rate P of 0.95. In our case, a non-parametric approach was used, resulting in an upper tolerance limit of 28 mm. This value represents the accuracy of the resulting TPs cloud. It is useful to analyse the semi-axis lengths associated with each individual TPs graphically in colour scale (Figure 4). The figure shows that the central zone of the model is characterized by higher accuracy (8 mm), this is attributable to a greater number of frames covering this zone and the presence of an appropriate number of oblique frames. The most accurate area is at ground level, this is attributable to the images being taken at street level, at human height. This result is also in accord with the SDs calculated on the GCPs. On the edges of the TPs cloud, however, the values are around 20 mm, with maximum values of about 28 mm, probably due to the lower presence of nadiral and oblique images.  To better highlight how the calibration of accuracy parameters can affect the quality of the sparse point cloud, Figure 5 shows the TPs cloud with associated major semi-axis lengths of the ellipsoid calculated for each single TP without the procedure being implemented. Although the maximum error (15 mm) is smaller than in the previous configuration (28 mm), the graphical distribution of errors is not representative of the acquisition scheme used and of the total number of images; looking at the figure it is not possible to identify areas characterized by different levels of accuracy. It should be pointed out again that the accuracy of the model, as is known, increases as the number of images, both oblique and nadiral, increases, as in the central part of the façade (Figure 4, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume X-M-1-2023 29th CIPA Symposium "Documenting, Understanding, Preserving Cultural Heritage: Humanities and Digital Technologies for Shaping the Future", 25-30 June 2023, Florence, Italy blue area). Figure 4 shows a plausible error distribution in agreement with the distribution and asset of the images, which is not evident in Figure 5, which shows a homogeneous error distribution far from the real world. For such reasons, the proposed procedure of analyzing the parameters that are involved in the photogrammetric process is relevant to the correct estimation of the accuracy of TPs.

CONCLUSIONS
The aim of this contribution was to provide a workflow for obtaining an accurate TPs cloud. For this reason, an attempt was made to control as many parameters as possible that govern the entire photogrammetric process. For our application, we limited ourselves to the Structure-from-Motion (SfM) phase of the process, leaving out that this is followed by the reconstruction of the dense cloud (Multi-View Stereo). The results show how these parameters, if not properly considered, can significantly influence the dense cloud and consequently the final 3D model. The proposed methodology also allows the control of a global parameter of accuracy for the TPs, which is often not considered, mainly if the weights of the observations are not properly calibrated.