Comparison of Different Algorithms to Orthorectify WorldView-2 Satellite Imagery

Due to their level of spatial detail (pixel dimensions equal to or less than 1 m), very high-resolution satellite images (VHRSIs) need particular georeferencing and geometric corrections which require careful orthorectification. Although there are several dedicated algorithms, mainly commercial and free software for geographic information system (GIS) and remote sensing applications, the quality of the results may be inadequate in terms of the representation scale for which these images are intended. This paper compares the most common orthorectification algorithms in order to define the best approach for VHRSIs. Both empirical models (such as 2D polynomial functions, PFs; or 3D rational polynomial functions, RPFs) and rigorous physical and deterministic models (such as Toutin) are considered. Ground control points (GCPs) and check points (CPs)—whose positions in the image as, well as in the real world, are known—support algorithm applications. Tests were executed on a WorldView-2 (WV-2) panchromatic image of an area near the Gulf of Naples in Campania (Italy) to establish the best-performing algorithm. Combining 3D RPFs with 2D PFs produced the best results.


Introduction
Very high-resolution satellite images (VHRSIs) are widely used in several application fields, at both scientific and commercial levels, because of their detailed information (pixel dimensions equal to or less than 1 m in panchromatic band).
The geometric use of these images requires careful orientation and orthorectification in order to georeference the images and correct the geometric deformations introduced during acquisition [1].
Remotely sensed images contain such significant geometric distortions that they cannot be used directly with map-based products in a geographic information system (GIS) [2].This characteristic is particularly pertinent to VHRSIs for several reasons (e.g., the reduced dimensions of the pixels and the off-nadir viewing).
Such 2D/3D empirical models can be used when the parameters of the acquisition systems or a rigorous 3D physical model are not available.In addition, they are based on different mathematical functions (2D or 3D polynomial functions, 3D rational functions) [15].
Physical models (also called rigorous models) are based on a standard photogrammetric approach, where the collinearity equations link the image and the ground coordinates, and the parameters involved have a physical meaning.In addition, knowledge regarding the sensor on the specific satellite and orbit characteristics is required [16].
These 2D/3D empirical models do not reflect the real physical condition of imaging, but they define it approximately.Three different sensor replacement models are defined in ISO 19130 "Geographic Information-Imagery Sensor Models for Geopositioning": (1) simple polynomial fitting; (2) rational function model-RFM; and (3) universal sensor model-USM [14].Simple polynomial fitting is useful only for low-and medium-resolution satellite images [17] and can be adopted for VHRSIs in small area mapping.In fact, we can also assume that the sensor moves linearly in space and that the attitude is almost unchanged [18].The relationship between the position of a point on the image and the corresponding point on the object is represented by polynomial functions.The coefficients of each polynomial are computed using an appropriate number of ground control points (GCPs) with known locations (X, Y) on the map and easy detections (X', Y') in the image [19].
An RFM matches the image and object spaces through a ratio of two polynomial functions to compute the image row, and a similar ratio to compute the image column [20].Also, in this case an appropriate number of GCPs is necessary to calculate the coefficients of each polynomial, named RPCs (rational polynomial coefficients).RPCs are widely adopted by image vendors and government agencies, who do not want to provide any satellite/sensor information with the image, and by commercial workstation suppliers; consequently, the various image vendors around the world provide third-order RF parameters [21].RFMs have been used to approximate physical sensor models for many years due to their capability to maintain the full accuracy of different physical sensor models, their unique characteristic of sensor independence, and real-time calculation [22].
With USMs, an image scene is divided into multiple subdivisions or sections and RFM is applied to each section; the modified polynomials are of a variable order (up to the fifth power of the two horizontal coordinates and up to the third power for elevation); the denominator polynomials are usually omitted; image positional correction tables can be incorporated; and the ground coordinates can be expressed in any ground coordinate system [23].
Several applications about the use of 2D/3D empirical models as well as 2D/3D physical and deterministic models to orthorectify VHRSIs are available in literature.Chmiel et al. [24] compared rigorous approach and RFM to achieve the best geometric accuracy in orthorectified imagery products obtained from IKONOS (Geo Ortho-kit), QuickBird (Standard OrthoReady), and EROS 1A level data.Wolniewicz [25] analyzed the influence of the number and distribution of GCPs on the planimetric accuracy reached in the final QuickBird-and IKONOS-corrected panchromatic image using different orthorectification methods.Afify and Zhang [26] assessed the attainable geometric accuracies of RFMs applied to IKONOS and Quickbird images when user adjustments are introduced.Aguilar et al. [27] compared results recorded for IKONOS (Geo Ortho-kit) and QuickBird basic imagery, applying two different sensor models: (1) a 3D rational function refined by the user with zero-or first-order polynomial adjustment and (2) the 3D Toutin physical model.The influence of the number of GCPs on the geometric accuracy reached with RFM applied to GeoEye-1 images was analyzed by Maglione et al. [28].This paper describes orthorectification that was conducted on WorldView-2 (WV-2) imagery using several algorithms which are available in common commercial and free software.Both rigorous and empirical approaches are considered and positional accuracies of the results are compared and discussed in order to define the best algorithm or the best combination.To our knowledge this, is the first work that compares several methods (including polynomial functions and rational functions, with variable number of GCPs) on WV-2 imagery.The paper is organized as follows.Section 2 describes the data and methods; various algorithms are applied to WV-2 imagery regarding the area around Vesuvius.Section 3 presents and discusses the results.Finally, Section 4 draws some conclusions.

The Area Considered
In this study, WV-2 imagery concerning the area around Vesuvius (Campania, Italy) was considered (Figure 1).The dataset was acquired on 22 April 2015 and supplied by the provider as Ortho Ready Standard product (ORS2A), which means it had already been rectified using the average terrain elevation of the scene (low accuracy) [33].The scene is georeferenced in the UTM/WGS84 33 N zone T coordinates system and extends from 40 discussed in order to define the best algorithm or the best combination.To our knowledge this, is the first work that compares several methods (including polynomial functions and rational functions, with variable number of GCPs) on WV-2 imagery.The paper is organized as follows.Section 2 describes the data and methods; various algorithms are applied to WV-2 imagery regarding the area around Vesuvius.Section 3 presents and discusses the results.Finally, Section 4 draws some conclusions.

The Area Considered
In this study, WV-2 imagery concerning the area around Vesuvius (Campania, Italy) was considered (Figure 1).The dataset was acquired on 22 April 2015 and supplied by the provider as Ortho Ready Standard product (ORS2A), which means it had already been rectified using the average terrain elevation of the scene (low accuracy) [33].The scene is georeferenced in the UTM/WGS84 33 N zone T coordinates system and extends from 40°46′22.06′′north to 40°51′46.53′′north in latitude and 14°20′2.74′′east to 14°24′40.70′′east in longitude.This area is particularly interesting due to both its anthropological and natural features.In fact, most of the area is within the Vesuvius National Park, founded in 5 June 1995, which houses 610 vegetable species and 229 vertebrates and invertebrates [34].The area also has one of the highest urban densities in Europe and includes the archeological remains of Pompeii, the ancient town that was destroyed by the eruption of Vesuvius in 79 AD.Thus, a high-performance VHRSI such as WV-2 would clearly support research activities in many fields, if corrected and georeferenced.

Algorithms for Orthorectification
In this section, several algorithms that are present in the software most commonly used for image orthorectification are described.
According to Toutin [2], whatever the mathematical functions used, the geometric correction method and processing steps are: acquisition of image; acquisition of the GCPs and check points (CPs) with image coordinates and map coordinates X, Y, Z; computation of the unknown parameters of the mathematical functions used for the geometric correction model; and image rectification with or without digital elevation model (DEM).
The two-dimensional (2D) polynomial model has been adopted for image rectification since the 1970s [35].It relates images coordinates (X', Y') to cartographic coordinates (X, Y) by the following equations: If n is the order of the equation, the following relations are valid: The 2D polynomial function model of the first order (6 parameters) corrects the image by a rotation, a translation in x and y, scaling along both axes and an obliquitous transformation.The 2D first-order polynomial function is also called affine transformation.The 2D polynomial function model of second order (12 parameters) adds to the previous transformation's correction of the torsion and convexity along both axes.A cubic polynomial function (20 parameters) further increases the correction flexibility, in a way that does not necessarily correspond to any physical reality of the image acquisition system [36].
By substituting the coordinate values of GCPs in Equation ( 1), polynomial coefficients (a lij and b lij ) can be calculated.Consequently, the new location of each pixel in the cartographic plane can be defined.The minimum number of GCPs needed is established by the order n of the polynomial functions according to the formula [37]: To increase the positional accuracy of the resulting image, a greater number of GCPs is chosen, which are distributed regularly.Thus differences between map coordinates and corrected image coordinates for the same features tend to be reduced.To evaluate the quality of the results, errors not only for GCPs but also for other points, named check points (CPs), are considered [2,38].
Using polynomial functions is a very uncomplicated method to orthorectify VHRSI: often considered outdated, they correct for basic planimetric distortion at the GCPs; because they do not take ground elevation into consideration, they are limited to small and flat areas [39].
3D rational polynomial functions (RPFs) may be more useful.They define a relationship between the image coordinates (x', y') and the 3D object coordinates (X, Y, Z) [23,40]. Specifically: where P n 1 , P n 2 , P n 3 , P n 4 are usually cubic polynomials.Each of these includes 20 coefficients and can be expressed as: where: Equations in (4) are known in the literature as Upward RFM, as they enable the image coordinates to be obtained starting from the 3D coordinates of a ground point [41].
Substituting P in (4) with the polynomials in ( 5) and eliminating the first coefficient in the denominator, there are 39 RPF coefficients in each equation: 20 coefficients in the numerator and 19 (and the constant 1) in the denominator [42].In order to solve the 78 coefficients, at least 39 GCPs are required [43].
The 78 coefficients are calculated by the data provider considering the position of the satellite at the time of image acquisition.They are then included in the RPC (rational polynomial coefficient) file.
However, they can also be calculated using GCPs [43], as in the case above for polynomial functions (PFs).At least 39 of them are needed.In fact, if the camera model is not available, the ground control points need to be selected in a conventional way; that is, through collimation of the homologous points on the cartography or DEM (digital elevation model) or through specific GPS (Global Positioning System) survey campaigns [44].Usually the Equations in (4) are resolved using a least squares iterative process, on the basis of the measurement of a large number of GCPs.
The accuracy of the results depends on the number and the distribution of GCPs.Several GCPs (more than 39) with a regular planimetric and altimetric distribution contribute to the high quality of the results [17,45,46].A DEM of the whole area is also required for RPFs.
By applying PFs as well as RPFs to the primary image, a matrix of "empty" cells is computed.To calculate the radiometric value to be assigned to each pixel, a resampling method, such as nearest neighbor, bilinear interpolation, or cubic convolution, is used [47,48].The nearest-neighbor algorithm is used for this application because it preserves the original radiometric values, which is fundamental for further image processes (calculation of vegetation indexes, classifications, etc.).
Several algorithms are described in the literature for VHRSI orthorectification using rigorous physical models.Based on the collinearity equations, they reconstruct the geometry of the scene during the image acquisition through the knowledge of parameters concerning both the platform and the sensor.
The traditional photogrammetric method of transforming from object space to image space using collinearity equations is highly suited to frame cameras, but not to pushbroom sensor products [23].In fact, unlike the traditional frame-based aerial photos, each line of the linear array image is collected in a pushbroom fashion at a different instant of time; therefore, the perspective geometry is only valid for each line whereas it is close to a parallel projection in along-track direction; in addition, for each line, there is a different set of (time-dependent) values for the exterior orientation elements [49].For 1-meter GSD pushbroom sensors such as IKONOS, fully parametrized camera models are extremely complex and difficult to implement; this is highlighted in the IKONOS System Geometric and Mathematical Model document that consists of 183 pages [50].
Other examples of rigorous models are mentioned in the literature, i.e., the model implemented in SISAR, which is a software for high-resolution satellite images at the Area di Geodesia e Geomatica-Sapienza Università di Roma.The approximate values of the physical parameters can be computed thanks to the information contained in the metadata file provided with each image and corrected by a least squares estimation process based on an appropriate number of GCPs [55].
In this work orthorectification processes were carried out on the WV-2 panchromatic image using PCI Geomatica OrthoEngine Version 2015.The following algorithms were applied: Both GCPs and CPs were homogeneously distributed over the study area and located on well-defined features.Clearly, the source, accuracy, distribution and number of GCPs and CPs are very important for a correct orthorectification [2].To obtain a reliable measurement of GCPs and CPs with adequate planimetric accuracy, their (x, y) coordinates in UTM-WGS84 were derived from orthophotos of Campania with 0.20 m resolution (nominal scale: 1:5000).It should be noted that errors may occur when manually selecting corresponding points.Thus, specific GCPs and CPs were detected in well-defined locations such as road junctions and corners of buildings or swimming pools.For any ambiguous points, Canny edge detection algorithm [56] was used to support the identification of the corresponding points.In accordance with the literature [25], this helps in ensuring that the accuracy of the GCP and CP identification on the imagery was definitely below one pixel.Figure 2 reports an example of CP identification.
Other examples of rigorous models are mentioned in the literature, i.e., the model implemented in SISAR, which is a software for high-resolution satellite images at the Area di Geodesia e Geomatica-Sapienza Università di Roma.The approximate values of the physical parameters can be computed thanks to the information contained in the metadata file provided with each image and corrected by a least squares estimation process based on an appropriate number of GCPs [55].
In this work orthorectification processes were carried out on the WV-2 panchromatic image using PCI Geomatica OrthoEngine Version 2015.The following algorithms were applied: Both GCPs and CPs were homogeneously distributed over the study area and located on well-defined features.Clearly, the source, accuracy, distribution and number of GCPs and CPs are very important for a correct orthorectification [2].To obtain a reliable measurement of GCPs and CPs with adequate planimetric accuracy, their (x, y) coordinates in UTM-WGS84 were derived from orthophotos of Campania with 0.20 m resolution (nominal scale: 1:5000).It should be noted that errors may occur when manually selecting corresponding points.Thus, specific GCPs and CPs were detected in well-defined locations such as road junctions and corners of buildings or swimming pools.For any ambiguous points, Canny edge detection algorithm [56] was used to support the identification of the corresponding points.In accordance with the literature [25], this helps in ensuring that the accuracy of the GCP and CP identification on the imagery was definitely below one pixel.Figure 2   GCP and CP elevations were obtained by DEMs of the study area with a different horizontal resolution (cell size: 1 m, 20 m, 75 m).The DEM with a cell size of 1 m derived from laser-scanning data was supplied by the local administration (Provincia di Napoli) which indicates residuals of 0.16 m between the heights from the grid and those of the 3D points acquired with the RTK (real-time kinematic) survey.The DEM with a cell size of 20 m was provided by IGM (Istituto Geografico Militare), which indicates residuals of 7-10 m between the heights from the grid and those of the 3D points in plane and hilly zones.The DEM with a cell size of 75 m was provided by ISPRA (Istituto Superiore per la Protezione e la Ricerca Ambientale) and presents a vertical accuracy of 16-22 m in the considered area.
The flowchart in Figure 3 shows the inputs (orthophotos, WV-2 panchromatic image, DEMs), the algorithms used (2D PFs-fourth-and fifth-order; 3D RPFs with and without RPCs; Toutin's rigorous model; and a combination of 3D RPFs and 2D PFs), and the most important steps (GCP and CP residuals calculation, positional accuracy evaluation based on minimum, maximum, mean, standard deviation, RMS calculation, orthorectification process, and output).The orthorectification in two steps was introduced to increase the geopositional accuracy of the final WV-2 panchromatic image.In the first step, the best resulting 3D RPFs model (75 GCPs-15 CPs, DEM = 1 m) was applied.For the second step, first-order PFs using 15 GCPs was performed on the orthoimage resulting from the first step.

Results and Discussion
The positional accuracy of the primary panchromatic image was tested using 15 CPs. Figure 4 reports the CP distribution in the WV-2 primary panchromatic image.Table 1 shows the statistical values (minimum, maximum, mean, standard deviation, RMS) of errors in these CPs.The low level of positional accuracy of the primary image (RMS = 84.589m) is evident compared with the pixel dimensions.
Algorithms 2016, 9, 67 8 of 16 The orthorectification in two steps was introduced to increase the geopositional accuracy of the final WV-2 panchromatic image.In the first step, the best resulting 3D RPFs model (75 GCPs-15 CPs, DEM = 1 m) was applied.For the second step, first-order PFs using 15 GCPs was performed on the orthoimage resulting from the first step.

Results and Discussion
The positional accuracy of the primary panchromatic image was tested using 15 CPs. Figure 4 reports the CP distribution in the WV-2 primary panchromatic image.Table 1 shows the statistical values (minimum, maximum, mean, standard deviation, RMS) of errors in these CPs.The low level of positional accuracy of the primary image (RMS = 84.589m) is evident compared with the pixel dimensions.The performance of each rectification was evaluated by the residuals obtained for GCPs and CPs in XY direction.Figure 5 illustrates the distribution of all GCPs datasets (5-10-15-45-60-75) in WV-2 primary panchromatic image.A homogenous distribution was attempted for each group of GCPs.
For all the tests performed, 15 CPs (the same as the previous test on the primary image) were introduced.Statistical values of these errors are reported in the following tables.
Table 2 shows how for both fourth-and fifth-order PFs, 60 GCPs produce better results than 45 GCPs while the introduction of another 15 GCPs leads to small increases in accuracy.However, the results are not suitable (RMS for CPs = 10-25 m) in terms of the image geometric resolution (0.5 m).The performance of each rectification was evaluated by the residuals obtained for GCPs and CPs in XY direction.Figure 5 illustrates the distribution of all GCPs datasets (5-10-15-45-60-75) in WV-2 primary panchromatic image.A homogenous distribution was attempted for each group of GCPs.Tables 3 and 4 show how the original RPC file supplied by image provider limits CP residuals using a low number of GCPs (i.e., when a DEM with a cell size of 1 m is used, only 5 GCPs mean that the RMS on CPs can be limited to 3.272 m).Although height accuracy is less important than planimetric accuracy, the influence of DEM resolution on elevation extraction is also evident (e.g., with 15 GCPs, the RMS residuals for CPs was 3.747 m using DEM with a 75 m cell size and 1.688 m using DEM with a 1 m cell size.For all the tests performed, 15 CPs (the same as the previous test on the primary image) were introduced.Statistical values of these errors are reported in the following tables.
Table 2 shows how for both fourth-and fifth-order PFs, 60 GCPs produce better results than 45 GCPs while the introduction of another 15 GCPs leads to small increases in accuracy.However, the results are not suitable (RMS for CPs = 10-25 m) in terms of the image geometric resolution (0.5 m).Tables 3 and 4 show how the original RPC file supplied by image provider limits CP residuals using a low number of GCPs (i.e., when a DEM with a cell size of 1 m is used, only 5 GCPs mean that the RMS on CPs can be limited to 3.272 m).Although height accuracy is less important than planimetric accuracy, the influence of DEM resolution on elevation extraction is also evident (e.g., with 15 GCPs, the RMS residuals for CPs was 3.747 m using DEM with a 75 m cell size and 1.688 m using DEM with a 1 m cell size.Table 5 shows that to achieve good results a great number of GCPs as well as a high-resolution DEM are required to calculate RPCs.Tables 6 and 7 show how Toutin's rigorous model obtains a good performance using a low number of GCPs and a high-resolution DEM (e.g., with 15 GCPs, RMS residuals for CPs was 4.129 m using DEM with a 75 m cell size and 1.737 m using DEM with a 1 m cell size).
Table 8 reports the results obtained using the combination of RPF and PF algorithms.The GCPs used in this case are illustrated in Figure 6.This approach achieves low residuals, comparable with the WV-2 cell size as well as avoiding other geometrical distortions associated with higher-order PFs (i.e., twist with second-order).A combination of RPF and PF algorithms gave the best results.Table 8 reports the results obtained using the combination of RPF and PF algorithms.The GCPs used in this case are illustrated in Figure 6.This approach achieves low residuals, comparable with the WV-2 cell size as well as avoiding other geometrical distortions associated with higher-order PFs (i.e., twist with second-order).A combination of RPF and PF algorithms gave the best results.

Conclusions
The tests presented in this paper for WV-2 data show that RPFs are suitable for VHRSI orthorectification.Using the RPC file supplied by the provider, they can also be applied without GCPs, however the results present considerable positional errors, and few GCPs are sufficient to achieve good results.If RPCs are calculated by the user, several GCPs are necessary to limit residual values for Ground and Check Points.
RPFs need 3D information: Planimetric coordinates are not sufficient for GCPs, and a DEM of the entire scene is indispensable to complete the orthorectification.However, Z accuracy can be lower than that of the X and Y coordinates.As a consequence, height values can also be extracted from the DEM if they are not available from and GPS surveys.Detailed maps (scale 1:5000 or greater) as well as remotely sensed images with the same resolution (or better) of the panchromatic image-but already orthorectified-are suitable for extracting planimetric coordinates of GCPs and CPs.
A low number of GCPs is sufficient to achieve a good performance with Toutin's rigorous model.The results obtained with WV-2 showed the analogous impact of DEM by increasing the geometric resolution for Z extraction for both RPFs as well as Toutin's rigorous model.
Although PFs are inadequate for high-resolution satellite images, also when a large number of GCPs are used, they improve the results of RPFs if applied as a second step in the orthorectification process based on the combination of these algorithms.Future research will be aimed at a comparison of similar algorithms extended to other VHRSIs such as WorldView-3, and GeoEye-2.

Figure 1 .
Figure 1.WorldView-2 (WV-2) image and its location in Campania.The coordinates reported on the map are referred to as UTM/WGS84 (33 N-zone T) and expressed in meters.

Figure 1 .
Figure 1.WorldView-2 (WV-2) image and its location in Campania.The coordinates reported on the map are referred to as UTM/WGS84 (33 N-zone T) and expressed in meters.

Figure 2 .
Figure 2. Check point (CP) identification: (a) CP in orthophoto; (b) corresponding CP in the WV-2 primary panchromatic image.The coordinates reported are referred to as UTM/WGS84 (33 N-zone T) and are expressed in meters.

Figure 2 .
Figure 2. Check point (CP) identification: (a) CP in orthophoto; (b) corresponding CP in the WV-2 primary panchromatic image.The coordinates reported are referred to as UTM/WGS84 (33 N-zone T) and are expressed in meters.
. The DEM with a cell size of 20 m was provided by IGM (Istituto Geografico Militare), which indicates residuals of 7-10 m between the heights from the grid and those of the 3D points in plane and hilly zones.The DEM with a cell size of 75 m was provided by ISPRA (Istituto Superiore per la Protezione e la Ricerca Ambientale) and presents a vertical accuracy of 16-22 m in the considered area.The flowchart in Figure3shows the inputs (orthophotos, WV-2 panchromatic image, DEMs), the algorithms used (2D PFs-fourth-and fifth-order; 3D RPFs with and without RPCs; Toutin's rigorous model; and a combination of 3D RPFs and 2D PFs), and the most important steps (GCP and CP residuals calculation, positional accuracy evaluation based on minimum, maximum, mean, standard deviation, RMS calculation, orthorectification process, and output).

Figure 4 .
Figure 4. 15 CP distribution in WV-2 primary panchromatic image.The coordinates reported are referred to as UTM/WGS84 (33 N-zone T) and are expressed in meters.

Figure 4 .
Figure 4. 15 CP distribution in WV-2 primary panchromatic image.The coordinates reported are referred to as UTM/WGS84 (33 N-zone T) and are expressed in meters.

Figure 6 .
Figure 6.The combination of 3D RPFs and 2D polynomial functions (PFs): (a) 75 GCPs used for 3D RPFs without RPCs; (b) 15 GCPs used for 2D PFs first order.The coordinates reported are referred to as UTM/WGS84 (33 N-zone T) and expressed are in meters.

Figure 6 .
Figure 6.The combination of 3D RPFs and 2D polynomial functions (PFs): (a) 75 GCPs used for 3D RPFs without RPCs; (b) 15 GCPs used for 2D PFs first order.The coordinates reported are referred to as UTM/WGS84 (33 N-zone T) and expressed are in meters.

Table 1 .
Residuals (in meters) obtained for CPs in the primary WV-2 panchromatic image.

Table 1 .
Residuals (in meters) obtained for CPs in the primary WV-2 panchromatic image.

Table 2 .
Residuals (in meters) obtained for ground control points (GCPs) and CPs using polynomial functions.