# Refined Satellite Image Orientation in the Free Open-Source Photogrammetric Tools Apero/micmac

E. Rupnik, M. Pierrot Deseilligny, A. Delorme, Y. Klinger

► To cite this version:

E. Rupnik, M. Pierrot Deseilligny, A. Delorme, Y. Klinger. Refined Satellite Image Orientation in the Free Open-Source Photogrammetric Tools Apero/micmac. ISPRS Annals of the Photogrammetry, 2016, III1, pp.83-90. 10.5194/isprs-annals-III-1-83-2016 . insu-03581656

HAL Id: insu-03581656

https://insu.hal.science/insu-03581656

Submitted on 20 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
REFINED SATELLITE IMAGE ORIENTATION IN THE FREE OPEN-SOURCE PHOTOGRAMMETRIC TOOLS APERO/MICMAC

E. Rupnik<sup>a,b*</sup>, M. Pierrot Deseilligny<sup>a</sup>, Arthur Delorme<sup>b</sup>, Yann Klinger<sup>b</sup>

<sup>a</sup> ENSG, École Nationale des Sciences Géographiques, Marne la Vallée, France -
(ewelina.rupnik, Marc.Pierrot-Deseilligny)@ensg.eu

<sup>b</sup> IPGP - Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Univ Paris Diderot, CNRS, Paris, France -
(delorme, klinger)@ipgp.fr

Commission I, WG I/4

**KEY WORDS:** Pushbroom, Orientation, Bundle Adjustment, Epipolar Geometry, MicMac, Open-source

**ABSTRACT:**

This publication presents the RPC-based bundle adjustment implemented in the freeware open-source photogrammetric tool Ap-
ero/MicMac. The bundle adjustment model is based on some polynomial correction functions, enriched with a physical constraint
that introduces the notion of a global sensor rotation into the model. The devised algorithms are evaluated against two datasets consist-
ing of two stereo and a triplet pair of the Pleiades images. Two sets of correction functions and a number of GCPs configurations are
examined. The obtained geo-referencing accuracy falls below the size of 1GSD.

**1. INTRODUCTION**

The value of high-resolution optical satellite imaging in earth sci-
ences is indisputable. Among the most common applications are
the mapping of natural hazards, monitoring the effects of the cli-
mate change, as well as large-scale land classification and sur-
veying.

Quantification of the ground displacements caused by Earthquakes
with e.g. SPOT-2 (Puymbroeck et al., 2000), SPOT-5 (Binet and
Bollinger, 2005; Vallage et al., 2015), Quickbird-II (Rosu et al.,
2015) has proven to give supixel accuracy result. That is, dis-
placements of few cm are possible to resolve. Within the context
of natural disasters and crisis areas, satellite imagery is also of
interest in mapping changes occurring around residential areas,
be it e.g. the collapsed buildings (Klonus et al., 2012), flooded
settlements, or structural damages like harbors (Cho et al., 2014).
Besides the short term observations, the long term phenomena
are observable too, for instance, in the study of glacier melts
caused by the changing climate (Jóhannesson et al., 2013). Last
but least, the satellite optical sensors offer a variety of ground
sample distances (GSD) and frequent revisit times thus are ad-
vantageous for cost-effective land-cover classification (Schindler,
2012), country/world-wide digital surface model (DSM) creation,
orthophoto generation (d'Angelo, 2013) or even updating of the
city building models/cadastral maps (Guerin et al., 2014).

**1.1 Modelling of the pushbroom sensors**

Pushbroom sensors are time-dependent sensors represented by
the central perspective projection in the across-track direction,
and the orthogonal projection in the along-track direction. As
such, they can be described *a*) rigorously by a physical model
(sensor's position, attitude etc.), or *b*) empirically by some ap-
proximating functions. Since the physical modelling varies from
sensor to sensor, and is mathematically complex, the common
practice is to work with the approximating functions which are
purely empirical and generated from the rigorous model (or dense
network of ground control points), yet reasonably accurate within

a specified validity zone. Further advantage of the empirical ap-
proach is in maintenance of its format due the fact that changes in
the sensor geometry/configuration do not necessitate any satellite
processing software re-implementation.

Examples of the empirical (a.k.a. generalized) functions are the
rational polynomial functions (RPCs), affine functions, or sim-
ply object space GRIDs. All three provide an efficient means
to functionally relate the image and object space. The 3D affine
model approximates the image to ground relation with orthogonal
projection, which is full-fledged for the satellites' narrow field of
views and high flying altitudes. See Fraser and Yamakawa (2003)
and Tao and Hu (2001) as good references on the affine and ra-
tional function models.

**1.2 Refining the RPC-based sensor orientation**

The general RPCs – nowadays customarily delivered with Pleiades,
SPOT 6/7, QuickBird, WorldView, GeoEye, IKONOS, Cartosat
etc. – of a stereo or a triplet pair are said to guarantee accu-
racy of plus or minus 10pix (CE90)(depending on the satellite
employed), corresponding to up to 12m in the ground space (Oh
and Lee, 2015). The absolute geo-location misalignment is due
to limitations of the sensor’s direct spatial positioning, primarily
the attitude, the position and velocity (Fraser and Hanley, 2005).
Notwithstanding, for large image blocks, multi-temporal acquisi-
tions, and/or images captured from different orbits, the RPCs typ-
ically reveal inconsistencies between the scenes, hindering their
suitability for application demanding high-precision orientation
parameters, e.g. for DSM generation with the semi-global match-
ing algorithms (d’Angelo, 2013). The RPC positioning is hence
subject to some refinement.

Consequently, direct and indirect methods of refining the RPCs
were conceived. The direct methods modify the original RPCs
(for example with the help of GCPs), whereas the indirect meth-
ods introduce auxiliary correction functions defined in the image
space. The approved methodology is that of indirectly correct-
ing for the differences between true and measured (derived from
RPCs) image coordinates with polynomial functions. Grodecki

<footnote>*Corresponding author</footnote>
and Dial (2003) and Fraser and Hanley (2003) independently proposed the method and laid down the recipe on how to estimate the parameters of the functions in a bundle adjustment routine. For a good overview of all the available refinement approaches the reader is referred to Xiong and Zhang (2009).

## 1.3 Objectives

The principal objectives of the work presented here were to a) have a free open-source software that is capable of handling the push-broom sensor data, b) include the pushbroom sensor orientation refinement within one implementation framework of the existing tools of Apero thus allow for future fusion of frame and push-broom sensors, and c) conceive a generic refinement algorithm that can handle all of the nowadays available pushbroom sensors (Pleiades, SPOT, Worldview, Quickbird, Ikonos, Aster, etc.).

This research contributes to the approved methodology of empirically correcting the sensor orientation parameters by adding a physical piece of information into the bundle adjustment (BA), be it a global sensor rotation. Consequently, it shall lead to the strengthening the BA when none of very few GCPs are available.

Section 2. presents the devised adjustment approach and familiarizes the reader with one evaluation method, namely the residual transverse parallax. The experiments and the commentary on the conducted tests for Limoge and Oregon datasets are presented in Subsections 3.1, 3.2. Subsection 3.3 examines in short the effect of weighting parameters on the outcome BA accuracy.

## 2. METHODS

When devising the orientation correction functions, the following was assumed,

(a) the satellite trajectory is precise,

(b) only the calibration and attitude are incorrect, thus

(c) the refinement can be modelled by a smooth 2D polynomial in x, and y of the image space.

The necessary algorithm input data are the images with their respective RPCs, Eqs. (1)-(3), and some homologous points between them. The GCPs are optional in the adjustment.

$$Sample = x(\phi, \lambda, h) = \frac{P_1(\phi, \lambda, h)}{P_2(\phi, \lambda, h)} \quad (1)$$

$$Line = y(\phi, \lambda, h) = \frac{P_3(\phi, \lambda, h)}{P_4(\phi, \lambda, h)} \quad (2)$$

$$P(\phi, \lambda, h) = \sum_{i=0}^{m} \sum_{j=0}^{n} \sum_{k=0}^{o} a_{ijk} \phi^i \lambda^j h^k \quad (3)$$

where $P_i$ are the nominator and denominator functions of the RPC equations, and $m, n, \in (0, 3)$.

### 2.1 Polynomial compensation with a physical constraint

Analytically, the problem of sensor orientation refinement is formulated in two condition (observations) and three constraint equations.

<figure><img src="image_2.png" /><figcaption>Figure 1: The relationship between the change in sensor attitude $\delta\alpha_i$, and the arising 2D pixel displacements $R_i$.</figcaption></figure>

**Conditions** The observations consist of the bespoken polynomial functions, so it is the well-established refinement method described in Grodecki and Dial (2003) and Fraser and Hanley (2003) (a.k.a. the *Bias-compensated RPC model*). For every homologous point, their image coordinates accounts to the sum of the measured coordinates and their adjustable correction functions, being the $(m+n)^{th}$ order polynomials; see Eqs. (4)-(7).

$$x = x(\phi, \lambda, h) + D_x(x, y) \quad (4)$$

$$y = y(\phi, \lambda, h) + D_y(x, y) \quad (5)$$

$$D_x(x, y) = \sum_{i=0}^{m} \sum_{j=0}^{n} a_{ij} \cdot x^i y^j \quad (6)$$

$$D_y(x, y) = \sum_{i=0}^{m} \sum_{j=0}^{n} b_{ij} \cdot x^i y^j \quad (7)$$

where $a_{ij}$ and $b_{ij}$ are the polynomial coefficients to be estimated.

**Constraints** Defining the orientation improvement in terms of polynomials, as described above, is risky for the arising over-parameterization. The constraint equations are introduced into the system so as to reduce that risk. The first and most significant is the rotation constraint in Eq. (8), stating that there exists a global 3D rotation of the sensor, which causes a displacement of the image pixels being equal the displacements caused by the 2D polynomial defined in Eqs. (6), (7). Fig. 1 renders a graphical representation of the constraint. The small angles $\delta\alpha_i$ along respective coordinate axes, induce a 2D pixel motion, displayed in the displacement fields maps $R_i$. The composition of the three motions caused by the three angles gives rise to the total pixel displacement and is modelled by a 2D deforming polynomial.

The second constraint, specified in Eq. (9) ascertains that the displacements remain close to zero. Ultimately, Eq. (10) controls the convergence rate, declaring the difference in displacements between subsequent iterations be small too.

The $pds_{rot}$, $pds_{nul}$, $pds_{last}$ functions in Eqs. (8)-(10) are non-parameterized weight functions, whereas the $pds_{Nb}$ in Eq. (11) limits the influence of this particular category of equations.
$$ \sum_{x=0}^{m} \sum_{y=0}^{n} \left( \sum_{k}^{\{x,y\}} D_k(x, y) - \sum_{l=0}^{2} \delta \alpha_l R_{x,y,l} \right)^2 \cdot pds_{rot} \cdot pds_{Nb} = 0 \quad (8) $$

$$ \sum_{x=0}^{m} \sum_{y=0}^{n} \sum_{k}^{\{x,y\}} D_k(x, y)^2 \cdot pds_{nul} \cdot pds_{Nb} = 0 \quad (9) $$

$$ \sum_{it=1}^{iter} \sum_{x=0}^{m} \sum_{y=0}^{n} \sum_{k}^{\{x,y\}} \left( D_k(x, y)^{it} - D_k(x, y)^{it-1} \right)^2 \cdot pds_{last} \cdot pds_{Nb} = 0 \quad (10) $$

$$ pds_{Nb} = \frac{\sum pds_{all}}{s \cdot t} \quad (11) $$

where s, t correspond to a grid defined over the image size, and $pds_{all}$ is the sum of weights of all observation equations (typically the tie points and GCPs). Note that MicMac defines separate $pds_{Nb}$ weights for each of its observations (conditions) and constraints; refer to the manual (MicMac, Apero) for the comprehensive equations. The bundle adjustment in MicMac internally recomputes the initial RPCs to convert and subsequently work on object coordinates defined in a chosen projection system.

## 2.2 Residual transverse parallax

The residual transverse parallax of a stereo pair is a means to evaluate the internal accuracy of the sensor's orientations (the combined interior and exterior orientation). For an ideal pair, irrespective of whether the orientations are known through RPCs or the physical parameters, the homologous points should lie along their respective epipolar curves. In other words, the distance of a point in the right image to the epipolar curve of the corresponding point in the left image, shall equal zero (equivalent to no transverse parallax).

**Semi-global matching** To quantify the remaining transverse parallax MicMac performs semi-global matching (Hirschmüller, 2008; Pierrot-Deseilligny and Paparoditis, 2006) in the direction perpendicular to the epipolar curve. The similarity measure used is the normalized cross-correlation computed over a window. The usual values for the full-resolution matching step is 0.025 or 0.05 pix, and at each step the matching occurs on an area of 0.05 x 0.2 pix in the along and transverse epipolar line directions.

**Epipolar curves generation** To ease the computational effort, prior to the matching, the images are resampled to their epipolar geometries. It is known that, unlike for frame cameras, the epipolar curves of pushbroom sensors a) do not follow straight lines, and b) globally there exist no epipolar curve pairs (Gupta et al., 1997). Nonetheless, locally – within some limited depths – the curves can be approximated as straight lines and as such their corresponding epipolar lines in the left and right images are conjugate (Oh, 2011).

The MicMac's strategy commences with the generation of the conjugate piecewise linear epipolar curves for a 2D image grid. Every grid point commits to an epipolar line originating from the intersection of this grid point with the minimum and maximum Z-coordinates in the global reference frame. The Z-coordinates are known from the provided RPC data, and the entire procedure is conducted symmetrically for the left and the right image. Typically, the global direction of all epipolar curves does not follow the direction of the image rows. To make it compliant with the epipolar condition, hence bring it parallel to the image rows, the curves are approximately rotated so that the points in the left and the right images are row-wise aligned. The last step does

the actual resampling of images with a polynomial function by default of 9th degree.

The recovery of epipolar curves in MicMac resembles that of (Oh, 2011), thus no further implementation details will be discussed herewith.

## 3. EXPERIMENTS

The orientation refinement is tested on two datasets acquired with the Pleiades satellites. The first dataset – Limoge – is composed of four images, i.e. along (B/H = 0.15 and B/H = 0.3) and across-track (B/H = 0.2) stereo pairs; the second dataset – Oregon – is a along-track image triplet (B/H = 0.1). Both datasets capture a scene occupied by urban and rural landscapes. Evaluation of the results is based on

(a) $\sigma_0$ – bundle adjustment (BA) standard deviation,

(b) $\mu_{\Delta X, \Delta Y, \Delta Z}$ – mean discrepancies between BA-estimated and true values of the check points (CPs),

(c) $P_{x2}$ – the residual transverse parallax computed on the total image (cf. Subsection 2.2).

<figure><img src="image_3.png" /><figcaption>Figure 2: Limoge; sample GCPs in (a) satellite, (b) aerial images.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Oregon; sample GCPs in (a) satellite, (b) LiDAR pc1.</figcaption></figure>
The GCPs on the Limoge dataset were obtained by triangulation on aerial images with a ground sampling distance (GSD) of 0.3m (cf. Fig. 2). The estimated planimetric point accuracy is about 0.5 GSD, and the altitude accuracy is about 1.0 GSD. The testing scenarios differ in terms of the number/distribution of the employed GCPs (1–6 denoted with letters A–E; cf. Fig. 5) and the degree of the deforming polynomial (0 and 1st degree denoted with the respective numbers and placed after the letters indicating the GCPs scenario, e.g. A1 for degree 1 in GCP scenario A). Scenarios X include no GCPs in the processing. The total number of the control points equals 36.

The GCPs on the Oregon dataset were obtained from a 2x2m LiDAR point cloud data available from the *OpenTopography*<fnref n="1" /> platform (cf. Fig. 3). Due to the problematic accurate point identification in both the satellite images and the point clouds, the evaluation proceeds solely on the Z-coordinates of control points. In analogy to the Limoge data, the testing scenarios differ in GCPs number/distribution (1–6 GCPs, cf. Fig. 6), and the degree of the deforming polynomial. The naming convention of particular scenarios follows that of the Limoge dataset, too. The total number of the control points equals 19.

Taking the physical viewpoint, the 0-degree polynomial absorbs the in-track and across-track errors hence it models an image by two parameters and the authors refer to it as the shift model. The extra four parameters furnished by the 1st-degree polynomial absorb the surplus gyro drift, radial ephemeris error and leftover interior orientation errors. It is referred to as the drift model (Grodecki and Dial, 2003). The employed weighting parameters were: $pds_{rot} = 0.02$, $pds_{last} = 1e-15$, $pds_{nul} = 0$, $s \cdot t = 20$ and remained unchanged throughout the processing. In the X scenarios, the $pds_{Nb}$ for the condition equations on all tie points was set to 1000, and for scenarios A–E it was equal to 10. Bundle adjustment on both datasets was carried out in the UTM projection system.

## 3.1 Results on Limoge

In Tab. 2 and Fig. 5 the BA results with different GCP configurations are depicted. Two correction parameter sets were tested: modelling the shift, and the drift. The geo-positioning accuracy prior to the adjustment is listed at the top of the table. Notably, the mean reprojection error in the range of 3 pixels, and the significant differences on CPs suggest that both the relative and absolute orientations are apt for refinement.

Scenarios $X_0$ and $X_1$, correspond to the adjustment result without any GCPs. In this case, only the relative orientation between images will be improved, whereas in absolute, the result will remain related by a spatial similarity transformation. This test is to see whether no control information in the solution leads to undesired deformations in object space, and to verify that the imposed constraint do work as anticipated. To remove the effect of the absolute shift, the residuals on control points reported in Tab. 2 and Fig. 5 were aligned with the global frame by a translation (exclusively in $X_{0,1}$ cases).

The 0-degree polynomial perfectly models the relative orientation and does not introduce any systematic artefacts with respect to the global frame (the residuals most optimistic). The 1-degree polynomial result furnishes slightly higher residuals, nonetheless, still more than satisfactory and comparable to the best case with GCPs.

Including one GCP in the shift model compensates for the misclosures in relative orientation (see the $\sigma_0$ of the BA) and improves

the geo-location. The $B_0$ configuration barely affects the result, yet a more indicative reduction would be expected, were the GCPs distributed symmetrically on the diagonal of the area. The remaining GCPs configurations render equivalent results, implying that if adopting the shift model, as many as 3 GCPs suffice to obtain a geo-location accuracy that falls below 1.0 GSD.

The drift model with one and two GCPs clearly performs inferior to the shift model. The systematic error introduced in the adjustment (cf. Fig. 5(g)) is owing to overparameterizing the solution (6 unknowns per an image and only 3 observations for a single GCP) and the rather rigid weighting attached to the GCPs. It is believed to ameliorate if the control information is loosely weighted in the BA. According to Fraser and Hanley (2005), such a solution can be thought of as equivalent to free-network adjustment with inner constraints. Interestingly, the GCPs configuration with 3 points is commensurate with the result $E_0$, i.e. the shift model with the total of 6 GCPs. As a general note, the drift model delivers better accuracies with the increased number of GCPs.

Tab. 1 and Fig. 4 report on the remaining transverse parallax between images of across-track stereo pair, when 4 GCPs are employed, in the shift ($D_0$) and drift ($D_1$) models. In either case the parallaxes place below the BA standard deviation ($\sigma_0$), nonetheless, the reduced figures in the drift model indicate it is a more suitable choice for the given dataset. The visible sinusoidal pattern is possibly due to the sensor's small frequency vibrations, or the limited approximating capabilities of RPCs.

<figure><img src="image_5.png" /><figcaption>Figure 4: Limoge; residual transverse parallaxes computed for an across-track stereo pair with the matching interval of 0.025 pix. Left: the $D_0$; right: the $D_1$, scenarios. The parallax's magnitudes are found in Tab. 1.</figcaption></figure>

<table><caption>Table 1: Limoge; $\mu_{Px2}$, $\sigma_{Px2}$ – the mean and standard deviation of the residual parallax computed for across-track stereo pair.</caption><thead><tr><th></th><th>μ<sub>Px2</sub></th><th>σ<sub>Px2</sub></th></tr></thead><tbody><tr><td>D<sub>0</sub></td><td>0.00</td><td>0.06</td></tr><tr><td>D<sub>1</sub></td><td>0.02</td><td>0.03</td></tr></tbody></table>

## 3.2 Results on Oregon

As in the previous case, the Oregon dataset was tested with two sets of correction parameters – the shift and the drift, the results of which are shown in Tab. 3 and Fig. 6. Here, unlike in the Limoge dataset where 4 images were acquired from two different orbits, the triplet stereo locates on a single orbit and it manifests in a very good relative orientation, still before the BA (the points' reprojection error of 0.5pix). Inclusion of one GCP in the shift model enhances the Z-accuracy on CPs to ≈1GSD, reaching the best result when 3 GCPs are used. With 4 and 6 GCPs the Z-accuracy slightly deteriorates. The viable explanation for this happening is the hampered identification of corresponding points in the Pleiades images and the LiDAR point cloud data. For the

<footnote n="1">¹http://www.opentopography.org</footnote>
<table><caption>Table 2: Limoge; $\sigma_0$ – bundle adjustment standard deviation, $\mu_{\Delta X, \Delta Y, \Delta Z}^{gt}$ – mean discrepancies between BA-estimated and true values on GCPs, and $\mu_{\Delta X, \Delta Y, \Delta Z}^{cp}$ – on CPs. The Y-coordinate runs along-track, the X-coordinate runs across-track.</caption><thead><tr><th></th><th>#GCP/CP</th><th>σ<sub>0</sub><br>[pix]</th><th>μ<sup>gt</sup><sub>ΔX</sub><br>[m]</th><th>μ<sup>gt</sup><sub>ΔY</sub><br>[m]</th><th>μ<sup>gt</sup><sub>ΔZ</sub><br>[m]</th><th>μ<sup>cp</sup><sub>ΔX</sub><br>[m]</th><th>μ<sup>cp</sup><sub>ΔY</sub><br>[m]</th><th>σ<sup>cp</sup><sub>ΔZ</sub><br>[m]</th></tr></thead><tbody><tr><td>-</td><td>0/36</td><td>3.0</td><td>-</td><td>-</td><td>-</td><td>2.3</td><td>3.2</td><td>12.4</td></tr><tr><td>X<sub>0</sub></td><td>0/36</td><td>0.83</td><td>-</td><td>-</td><td>-</td><td>0.17</td><td>0.26</td><td>0.51</td></tr><tr><td>A<sub>0</sub></td><td>1/35</td><td></td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td><td>0.46</td><td>0.33</td><td>1.12</td></tr><tr><td>B<sub>0</sub></td><td>2/34</td><td></td><td>0.04</td><td>0.04</td><td>0.01</td><td>0.36</td><td>0.35</td><td>1.05</td></tr><tr><td>C<sub>0</sub></td><td>3/33</td><td>0.84</td><td>0.03</td><td>0.03</td><td></td><td>0.39</td><td>0.36</td><td>0.61</td></tr><tr><td>D<sub>0</sub></td><td>4/32</td><td></td><td>0.07</td><td>0.04</td><td>0.02</td><td>0.30</td><td>0.33</td><td>0.63</td></tr><tr><td>E<sub>0</sub></td><td>6/31</td><td></td><td>0.06</td><td>0.06</td><td></td><td>0.30</td><td>0.35</td><td>0.61</td></tr><tr><td>X<sub>1</sub></td><td>0/36</td><td>0.82</td><td>-</td><td>-</td><td>-</td><td>0.29</td><td>0.63</td><td>0.61</td></tr><tr><td>A<sub>1</sub></td><td>1/35</td><td></td><td></td><td></td><td></td><td>4.31</td><td>10.90</td><td>21.90</td></tr><tr><td>B<sub>1</sub></td><td>2/34</td><td></td><td>&lt;0.01</td><td>&lt;0.01</td><td></td><td>3.94</td><td>3.11</td><td>1.04</td></tr><tr><td>C<sub>1</sub></td><td>3/33</td><td>0.82</td><td></td><td></td><td>&lt;0.01</td><td>0.33</td><td>0.38</td><td>0.58</td></tr><tr><td>D<sub>1</sub></td><td>4/32</td><td></td><td>0.02</td><td>0.04</td><td></td><td>0.27</td><td>0.34</td><td>0.41</td></tr><tr><td>E<sub>1</sub></td><td>6/31</td><td></td><td>0.03</td><td>0.06</td><td></td><td>0.25</td><td>0.35</td><td>0.36</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Limoge; BA results in different test scenarios. Red triangular markers correspond with the GCPs, the black circular markers are the CPs, and the vectors indicate the discrepancies between the BA-estimated and true values of control points in object space. Note the different scales in (g).</figcaption></figure>

<table><caption>Table 3: Oregon; $\sigma_0$ – bundle adjustment standard deviation, $\mu_{\Delta Z}^{gt}$ – mean discrepancies between BA-estimated and true Z-values on GCPs, and $\mu_{\Delta Z}^{cp}$ – on CPs.; both derived from LiDAR point cloud.</caption><thead><tr><th></th><th>#GCP/CP</th><th>σ<sub>0</sub><br>[pix]</th><th>μ<sup>gt</sup><sub>ΔZ</sub><br>[m]</th><th>μ<sup>cp</sup><sub>ΔZ</sub><br>[m]</th></tr></thead><tbody><tr><td>-</td><td>0/19</td><td>0.50</td><td>0.30</td><td>11.58</td></tr><tr><td>A<sub>0</sub></td><td>1/18</td><td></td><td>&lt;0.01</td><td>0.65</td></tr><tr><td>B<sub>0</sub></td><td>2/17</td><td></td><td>&lt;0.01</td><td>0.47</td></tr><tr><td>C<sub>0</sub></td><td>3/16</td><td>0.57</td><td>0.01</td><td>0.40</td></tr><tr><td>D<sub>0</sub></td><td>4/15</td><td></td><td>&lt;0.01</td><td>0.41</td></tr><tr><td>E<sub>0</sub></td><td>6/14</td><td></td><td>&lt;0.01</td><td>0.46</td></tr><tr><td>A<sub>1</sub></td><td>1/18</td><td></td><td></td><td>56.05</td></tr><tr><td>B<sub>1</sub></td><td>2/17</td><td></td><td></td><td>14.63</td></tr><tr><td>C<sub>1</sub></td><td>3/16</td><td>0.57</td><td>&lt;0.01</td><td>0.47</td></tr><tr><td>D<sub>1</sub></td><td>4/15</td><td></td><td></td><td>0.48</td></tr><tr><td>E<sub>1</sub></td><td>6/14</td><td></td><td></td><td>0.50</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6: Oregon; BA results in different test scenarios. Red triangular markers correspond with the GCPs, the black circular markers are the CPs, and the blue circle radius' indicate the discrepancies between the BA-estimated and true values of LiDAR-derived control points. Note the different scales in (a) and (e).</figcaption></figure>

same reason the increased points' reprojection error ($\sigma_0$) is observed. The drift model – again, for the same reasons mentioned in the subsection on Limoge – diverges from the true solution if solely a GCP is incorporated in the processing. Two GCPs considerably improve the geo-location, albeit still exhibit systematic errors at the far ends of the scene (cf. Fig. 6(e)).

The result of dense matching in the along-the-transverse-parallax direction is presented in Tab. 4 and Fig. 7. The residual parallax – computed for the along-track stereo pair – shows a global zero, and fluctuates with $\approx 1/10$ pixel irrespective of the model employed. It implies that the shift model could be the preferred refinement method.

<figure><img src="image_8.png" /><figcaption>Figure 7: Oregon; residual transverse parallaxes computed for along-track stereo pair with the matching interval of 0.05 pix. Left: the $D_0$; right: the $D_1$, scenarios. The parallax's magnitudes are found in Tab. 4.</figcaption></figure>

<table><caption>Table 4: Oregon; $\muPx_2$, $\sigmaPx_2$ – the mean and standard deviation of the residual parallax computed for along-track stereo pair.</caption><thead><tr><th></th><th>μPx<sub>2</sub></th><th>σPx<sub>2</sub></th></tr></thead><tbody><tr><td></td><td>[pix]</td><td>[pix]</td></tr><tr><td>D<sub>0</sub></td><td>0.0</td><td>0.08</td></tr><tr><td>D<sub>1</sub></td><td>0.0</td><td>0.07</td></tr></tbody></table>

### 3.3 Tests on weighting parameters

Fig. 8 presents the variation of accuracy as a function of the weighting parameters $pds_{rot}$, and $pds_{last}$. In the performed tests, the $pds_{null}$ was set to 0 thus excluding the constraint on small deformations.

The datasets manifest most optimal results for low $pds_{last}$ weight, turning it into a soft constraint. The $pds_{rot}$ weight is more conservative, however, our experiences prove that if no control information is included in the bundle solution, its value should tend again toward a softer constraint.

<figure><img src="image_9.png" /><figcaption>Figure 8: Accuracy as a function of weighting parameters $pds_{rot}$, and $pds_{last}$, in scenario $D_0$. Top: Limoge; bottom: Oregon.</figcaption></figure>
## 4. CONCLUSIONS

This work presents the RPC-based bundle adjustment (BA) im-
plemented in the free open-source software Apero/MicMac. The
adjustment model estimates the parameters of some polynomial
correction functions, constraining their estimation by the global
sensor rotation. Two datasets consisting of 4 and 3 Pleiades im-
ages were used as the real cases to evaluate both, the metric ca-
pabilities of the imagery, and the devised algorithms themselves.

A number GCPs configurations scenarios were arranged to ex-
amine the behavior of the BA. For each scenario two sets of cor-
rection parameters were taken into account – considering only a
shift error in the RPCs, and a shift plus a drift. It has been shown
that as long as the acquisitions took place from the same orbit
(Oregon), the shift parameters were sufficient to refine the geo-
localization of the data. On the other hand, images captured from
neighboring orbits (Limoge) proclaimed higher accuracy posi-
tioning when more parameters entered the BA, i.e. with the drift
model.

The number of GCPs employed in the adjustment and their distri-
bution affects the quality of the outcome sensor’s geo-referencing.
Generally speaking, a minimum of 3 control points are required
to obtain accuracy in the range of 1GSD. The empirical findings
further signal that the drift model (1st-degree polynomial) should
not be adopted in scenarios with less than 3 GCPs as it leads
to overparameterization. It has also been shown that bundle ad-
justment without control information assures quality results, and
does not introduce unwanted deformations.

Future works include automated estimation of weighting param-
eters, as well as investigating the possibility of a combined treat-
ment of satellite and aerial optical imagery.

## ACKNOWLEDGEMENTS

The authors would like to thank CNES for funding this project
through the TOSCA program. We thank Gregoire Maillet and
Benjamin Ferrand from the Institut National de l'Information Géo-
graphique et Forestiére (IGN) for sharing the Limoge dataset.
The numerical computations were performed on the S-CAPAD
cluster processing facility at the Institute de Physique du Globe
de Paris.

## References

Binet, R. and Bollinger, L., 2005. Horizontal coseismic deforma-
tion of the 2003 Bam (Iran) earthquake measured from SPOT-
5 THR satellite imagery. Geophysical Research Letters 32(2),
pp. L02307.

Cho, K., Wakabayashi, H., Yang, C., Soergel, U., Lanaras, C.,
Baltsavias, E., Rupnik, E., Nex, F. and Remondino, F., 2014.
Rapidmap project for disaster monitoring. In: 35th Asian Con-
ference on Remote Sensing (ACRS), Nay Pyi Taw, Myanmar.

d'Angelo, P., 2013. Automatic orientation of large multitemporal
satellite image blocks. In: Proceedings of international sym-
posium on satellite mapping technology and application.

Fraser, C. and Yamakawa, T., 2003. Applicability of the affine
model for Ikonos image orientation over mountainous terrain.
In: Workshop on HRM from Space, pp. 6–8.

Fraser, C. S. and Hanley, H. B., 2003. Bias compensation in ratio-
nal functions for IKONOS satellite imagery. Photogrammetric
Engineering & Remote Sensing 69(1), pp. 53–57.

Fraser, C. S. and Hanley, H. B., 2005. Bias-compensated RPCs
for sensor orientation of high-resolution satellite imagery. Pho-
togrammetric Engineering & Remote Sensing 71(8), pp. 909–
915.

Grodecki, J. and Dial, G., 2003. Block adjustment of high-resolution satellite images described by rational polynomials. Photogrammetric Engineering & Remote Sensing 69(1), pp. 59–68.

Guerin, C., Binet, R. and Pierrot-Deseilligny, M., 2014. Automatic detection of elevation changes by differential DSM analysis: Application to urban areas. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of 7(10), pp. 4020–4037.

Gupta, R., Hartley, R. et al., 1997. Linear pushbroom cameras. Pattern Analysis and Machine Intelligence, IEEE Transactions on 19(9), pp. 963–975.

Hirschmüller, H., 2008. Stereo processing by semiglobal matching and mutual information. Pattern Analysis and Machine Intelligence, IEEE Transactions on 30(2), pp. 328–341.

Jóhannesson, T., Björnsson, H., Magnusson, E., Guðmundsson, S., Pálsson, F., Sigurðsson, O., Thorsteinsson, T. and Berthier, E., 2013. Ice-volume changes, bias estimation of mass-balance measurements and changes in subglacial lakes derived by lidar mapping of the surface of icelandic glaciers. Annals of Glaciology 54(63), pp. 63–74.

Klonus, S., Tomowski, D., Ehlers, M., Reinartz, P. and Michel, U., 2012. Combined edge segment texture analysis for the detection of damaged buildings in crisis areas. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of 5(4), pp. 1118–1128.

MicMac, Apero, Pastis and Other Beverages in a Nutshell, 2016.

Oh, J., 2011. Novel Approach to Epipolar Resampling of HRSI and Satellite Stereo Imagery-based Georeferencing of Aerial Images. PhD thesis, The Ohio State University.

Oh, J. and Lee, C., 2015. Automated bias-compensation of ratio-
nal polynomial coefficients of high resolution satellite imagery
based on topographic maps. ISPRS Journal of Photogramme-
try and Remote Sensing 100, pp. 14–22.

Pierrot-Deseilligny, M. and Paparoditis, N., 2006. A multireso-
lution and optimization-based image matching approach: An
application to surface reconstruction from SPOT5-HRS stereo
imagery. In: ISPRS Workshop On Topographic Mapping From
Space, Vol. 36(1), Ankara, Turkey.

Puymbroeck, N. V., Michel, R., Binet, R., Avouac, J.-P. and
Taboury, J., 2000. Measuring earthquakes from optical satellite
images. Appl. Opt. 39(20), pp. 3486–3494.

Rosu, A.-M., Pierrot-Deseilligny, M., Delorme, A., Binet, R. and
Klinger, Y., 2015. Measurement of ground displacement from
optical satellite image correlation using the free open-source
software MicMac. ISPRS Journal of Photogrammetry and Re-
mote Sensing 100, pp. 48–59.

Schindler, K., 2012. An overview and comparison of smooth
labeling methods for land-cover classification. Geoscience
and Remote Sensing, IEEE Transactions on 50(11), pp. 4534–
4545.

Tao, C. V. and Hu, Y., 2001. A comprehensive study of the rational function model for photogrammetric processing. Photogrammetric Engineering and Remote Sensing 67(12), pp. 1347–1358.
Vallage, A., Klinger, Y., Grandin, R., Bhat, H. and Pierrot-Deseilligny, M., 2015. Inelastic surface deformation during the 2013 Mw 7.7 Balochistan, Pakistan, earthquake. Geology 43(12), pp. 1079–1082.

Xiong, Z. and Zhang, Y., 2009. A generic method for rpc refinement using ground control information. Photogrammetric Engineering & Remote Sensing 75(9), pp. 1083–1092.