Network optimization and design in group-wise registration of terrain corrected satellite images

Abstract

Direct web-based access to ready-to-use huge archives of satellite images and cloud-based services for planetary-scale data processing (e.g., Google Earth Engine or Amazon S3) is making possible to analyze unprecedented amounts of remotely sensed images simultaneously. Multiple images can be exploited to improve traditional results achieved through on-premises (on-site) processing, coupling cloud offerings, and redundant image information.

This paper will introduce the concept of image network optimization for the case of registration problems based on groups of terrain-geocoded images. The particular case of multi-image registration will be discussed, notwithstanding the proposed approach can be extended to other practical issues, as illustrated in the paper. The concept of network design and optimization for satellite images is mathematically formulated and quantified with a multi-purpose objective function comprising precision, reliability, and cost. Results are illustrated with theory and numerical simulations carried out with a rigorous stochastic approach, in which the significance of the different input variables is estimated. The developed network-based approach allows one to reduce the number of external constraints, mainly focusing only on images and their increasing availability through web-services integrated by massive cloud computation capability.

Keywords

Cost design network optimization precision reliability

1. Introduction

1.1 Why use network-based analysis for remotely sensed data?

The availability of a massive amount of satellite images directly accessible through the Internet is changing the traditional processing approach of remotely sensed data (Zhu et al., 2009). Processing is moving from personal computers and workstations (local storage) to cloud computing services. For instance, Sentinel-2 archives organized into tiles using the Military grid system are available to anyone via Amazon S3. The Sentinel image browser is available at http://sentinel-pds.s3-website.eu-central-1.amazonaws.com/browser.html. Amazon Web services (AWS) made available every day Landsat 8 scenes acquired from 2015, often within hours of production. Another example is Google Earth (GE) Engine (https://earthengine.google.com/), which combines multi-petabyte catalogs of satellite imagery and geospatial datasets with planetary-scale analysis. Processing runs in parallel on the machine composing the Google cloud. GE Engine is a new tool to perform highly-interactive algorithm development at the global scale (Horowitz, 2015). Some applications carried out with GE Engine are described in Alonso et al. (2016), Goldblatt et al. (2016), Pulighe and Lupia (2017), and Warren et al. (2015).

The opportunity to work without downloading massive datasets allows users to deal with an unprecedented amount of images with variable spatial, temporal, radiometric, and spectral resolution. Improved data access is correlated to more exhaustive usages of data, in which the combined use of more images than those traditionally exploited is quite attractive. On the other hand, large datasets introduce relevant problems and require efficient and reliable processing algorithms. Big data has, therefore, become a favorite word in remote sensing and, at the same time, a practical challenge for the development of applications and services capable of delivering information for productive work (Ma et al., 2015). The well-known 6Vs (Volume, Variety, Velocity, Veracity, Validity, and Volatility) have to be translated from a technological problem to a set of reliable solutions for data processing (Nativi et al., 2015).

Recent work carried out in the field of close-range photogrammetry and computer vision (say in the latest 6–8 syears) has increased the number of users of image-based reconstruction. Software like PhotoScan, PhotoModeler, 3D Zephyr, ContextCapture, Pix4D Mapper (among the others) is now used by different operators of different fields. This has also proved that (i) facilitating data availability and access (in this case digital cameras) coupled with (ii) automated processing solutions (such as low cost 3D reconstruction software) open new opportunities not only for the specialists of the sector, but also for a wide community of users interested in digital reconstruction.

Another important consideration deserves to be mentioned. The availability of automated processing algorithms for 3D reconstruction has led to an increment of the number of image acquired and processed. This is due to two main reasons: automated orientation procedures (coined Structure from Motion in computer vision) require more images than those strictly needed for manual processing, where a human operator selects points with interactive measurements. Procedures for image matching and orientation based on scale invariant features require short baselines, increasing the final number of images. The second reason is less technical and depends on user’s experience. Operators not very expert in the field of image analysis tend to acquire several images for both simple and complex applications. The list of users includes geologists, engineers, archeologists, designers, architects, etc., who are probably interested in such techniques for productive purposes, but also casual snapshooters or professional photographers. Photogrammetrists and computer vision specialists are a part of the users of image reconstruction software, where the high level of automation achieved in image processing has allowed users to “ignore” basic rules for the efficient production of 3D models. The opportunity to upload images and run automatic data processing is surely quite tempting. Obviously, a lack of “photogrammetric experience” has also implications regarding metric accuracy, level of detail, completeness, processing time, etc. (Nocerino et al., 2014).

In the case of satellite images, the author’s intention is not to expect a massive usage among different operators with limited experience in remote sensing for the growing availability of new cloud-based services. While the production of a 3D model from a set of images is a tangible result usable in a variety of applications (from pure visualization purposes to advanced applications where the model is corrected and refined), remotely sensed data processing still require expert users who have to analyze and interpret the result to deliver reliable products. This is a remarkable difference compared to 3D modeling solution, in which a “tangible visual result” is directly available (like an orthophoto or a textured 3D model). On the other hand, it is clear that the availability of a huge number of images coupled with efficient techniques for data processing will foster a more productive use of space data, requiring new downstream services and applications where satellite images collected by different sensors will be processed to deliver products.

More attention is also required to the different phases of the image processing pipeline. Reliable products from images can be automatically generated when all the various processing steps have reached a significant technological maturity, facing the new challenges of a multi-sensor approach with images with variable geometric, radiometric, temporal, and spectral resolution.

The work carried out in this paper tries to extend the geometric registration problem from a traditional “reference-to-sensed” approach to the simultaneous exploitation of a complete time series, in which images are registered in a single step. In fact, accurate image registration is still a significant preprocessing step in remote sensing applications (Dawn et al., 2010). The lack of accurate image-to-image registration can result in significant errors in practical applications, as demonstrated by Townshend et al. (1992). Previous work (Barazzetti et al., 2014) has shown that a multi-image registration approach is feasible, more precise and robust than basic pairwise registration approaches. The availability of more than two images collected over the same area is an opportunity to achieve better results.

The aim of this paper is not an evaluation of specific case studies in different areas with different images. The work carried out wants to highlight the potential (and limitations) for the case of a generic network of satellite images, generalizing the approach from a single image pair to a complete network with an arbitrary number of images. The work wants to present the concept of network optimization and design applied to satellite image registration problems, which can be defined as the estimation of network quality before processing real data.

Results through numerical simulations and theoretical explanations will be discussed to clarify the different aspects of the design problem, testing the improved precision when different parameters are varied, such as the number of images, image-to-image points, and image-to-control points, as well as their relative percentage and distribution. Outcomes will be available through rigorous stochastic approaches, in which the variance-covariance matrix of unknown parameters will provide a numerical solution to validate the different configurations.

Although the proposed work focuses on network registration problems, it can be extended to different applications with satellite images. Starting from set of images, it is possible to improve the mathematical relationships (i.e., the functional model) between observed values and unknowns to take into consideration a complete network. For instance, the work presented in this paper could be extended to the relative normalization problem with pseudo-invariant features, where usually the same point is matched on single image pairs.

1.2 The multi-purpose objective function in satellite network design for data registration

We can define the aim of the optimization process for network-based image registration problems as the maximization of a multi-purpose objective function $F$ defined as:

$\displaystyle F=\alpha\frac{1}{\textit{precision}}+\beta\cdot\textit{% reliability}+\gamma\frac{1}{\textit{cost}}=\textit{max}$ (1)

where coefficients $\alpha$ , $\beta$ , $\gamma$ are normalization factors that make precision, reliability, and cost directly comparable. The use of such multi-purpose objective function has been proposed by different authors in previous work mainly related to network design problems for geodetic networks (Grafarend & Sansò, 1985; Ogundare, 2015). The idea behind the optimization problem is to have a precise network with a high reliability, i.e. network geometry should guarantee the detection and removal of gross errors (Eshagh & Kiamehr, 2013). At the same time, cost cannot be ignored during network design. Usually, the better the precision and reliability, the higher the cost.

Precision can be evaluated from the variance-covariance matrix of unknown parameters ${\bm{C}}_{xx}$ , including both image registration parameters and point coordinates. Matrix ${\bm{C}}_{xx}$ can be investigated to separate the effect of image and point parameters and to define specific criteria for comparison. For instance, precision can be expressed with the eigenvalues ${\bm{\lambda}}_{i}$ of point coordinates and specific criteria could be based on A-optimality [sum( $\lambda_{i}$ ) $=$ min], D-optimality [det( ${\bm{C}}_{xx}$ ) $=$ min], and E-optimality [max( $\lambda_{i}$ ) $=$ min]. The mathematical formulation for network-based optimization used in this work will also provide a solution for the choice of the parameters to be minimized as well as the criteria for simulation comparison.

Reliability is intended as the ability to detect and remove in an automated way gross errors (outliers) that play a fundamental role in the case of fully automated data processing. Modern techniques for outlier detection are based on sequential approaches based on robust estimators (such as the work initially presented in Fischler & Bolles, 1981), where a geometric model is coupled with iterative adjustment techniques based on random sampling. The aim is to find a set of inliers for the final estimation via least squares.

An additional consideration deserves to be mentioned. Precision, reliability, and accuracy are numerical indicators of random, gross and systematic errors in the network. Although this work will mainly focus on precision as an overall quality parameter of the network (see next sections), precision is not more important than reliability and accuracy. The choice of precision will be motivated by the intrinsic opportunity to specify (a priori) precision criteria. This aim can be achieved by exploiting the variance-covariance matrix, that can be analytically expressed as ${\bm{C}}_{xx}=({\bm{A}}^{T}{\bm{WA}})^{-1}$ , where ${\bm{A}}$ is the design matrix, and ${\bm{W}}$ is a weight matrix. Precision is, therefore, a function of network geometry and can be estimated (simulated) without acquiring real measurements. Network geometry encapsulated into the design matrix ${\bm{A}}$ provide all the required information to get an estimate of ${\bm{C}}_{xx}$ .

The cost can be quantified by considering CPU time (or the number of floating-point operations per second) needed to complete data processing. CPU cost is a parameter that cannot be ignored in the case of a large dataset of images. The formulation of a specific cost function has to take into consideration the number of images and points used. An initial formulation of the cost function could be based on a quadratic form, that depends on the combination of images. For instance, a dataset of $m$ images can be split into ( $m^{2}-m)$ /2 independent pairs, where the parameter $m^{2}$ becomes predominant for large datasets. Optimization solutions were developed in the field of computer vision/image analysis to reduce CPU time through approaches that do not exploit all image combinations, leading to a better cost function $m\cdot$ log( $m$ ), which however does not rigorously utilize all the possible image combinations. For this reason, the initial cost function will be assumed as quadratic and will be further investigated to determine a better model, as illustrated in the next section.

Different statistical parameters can be used to quantify the reliability of a network. On the other hand, recent experiences in the field of close-range photogrammetry have proved that reliable methods for automated image processing are now available. This is the case of photogrammetric applications where most steps of the processing workflow are fully automatic. As mentioned, different scientific papers and commercial software are an indicator of the significant maturity achievable in automated close-range image processing, where gross errors are detected and discarded exploiting robust estimators. Such algorithms can be adapted to deal with satellite images, as demonstrated in different scientific papers (e.g. Qiaoliang et al., 2009; Teke & Temizel, 2010; Zhao et al., 2009; Bouchiha & Besbes, 2013).

High reliability can be achieved by integrating several strategies for outlier rejection in the procedure. Previous work in close-range photogrammetry has demonstrated that outlier rejection strategies can rely on high breakdown point estimators (Christmann, 1996), which allow the detection of possible outliers in the observations. The idea of having robustness in the estimation is to have safeguards against deviations from the assumptions. Within robust estimators, gross errors are defined as observations that do not fit the stochastic model of the estimated parameters. Their efficiency depends on different factors but primarily on the percentage of outliers. Basically, there is a lack of repeatability because of their random way of selecting the points (random sampling). Compared to data snooping techniques (a statistical test of the normalized residuals), robust estimators do not provide a measure or a judgment about the quality of the found (or rejected) outliers.

For this reason, reliability can be assumed as a part of CPU cost because of the need to integrate specific, robust algorithms that increase CPU time. From this point of view, reliability can be integrated into (Eq. (1)) to become an effective part of cost:

$F=\alpha\frac{1}{\textit{precision}}+\gamma\frac{1}{\textit{cost}^{\ast}}=% \textit{max}$ (2)

where cost ${}^{\ast}$ includes the extra-cost of the robust analysis for outlier detection and removal. In the next sections, precision and cost* will be directly related to the number of images and points, considering the algorithms available for multi-image matching and the mathematical relationships for an integrated approach for registration based on image networks.

2. CPU cost in networks for satellite image registration

The design process of the network is related to the choice of images and points, as well as their location and distribution, the operator for image matching and the mathematical model for geometric registration. Different operators for image matching were developed for photogrammetric and remote sensing applications. An extensive review is presented in Gruen (2012). The mathematical model for image registration is the geometric relationship between corresponding points in different images (Brown, 1992; Le Moigne et al., 2011). In the case of image matching between two images only, commercial software packages are based on models with different geometric transformations (such as rotation $+$ translation $+$ scale or affine). Such models can provide good result for high-resolution satellite images (such as Landsat 8 and Sentinel-2), that are available as terrain corrected images. Similar considerations are valid for medium and low-resolution images. Different examples with different matching strategies are described in Castro and Morandi (1987), Goshtasby et al. (1986), Goshtasby (2005), Pluim et al. (2001), Price et al. (1985). The case of raw very high-resolution images (such as WorldView-3, Ikonos, Quickbird, etc.) is not taken into consideration in this work because of the different mathematical model based on rational polynomial functions (Poli & Toutin, 2012).

Let us consider a set of $m$ (high, medium, or low) resolution images. An affine transformation can describe the generic transformation between an image and a set of control points:

$\displaystyle\bar{x}_{i}=a_{j}x_{i}+b_{j}y_{i}+c_{j}$

(3) $\displaystyle\bar{y}_{i}=d_{j}x_{i}+e_{j}y_{i}+f_{j}$

where $j$ is the image that requires registration, and $i$ is the point index. This is the case of most commercial software, in which an image is assumed as reference (coordinates $\bar{x}_{i}$ , $\bar{y}_{i}$ are fixed quantities) and the sensed image is transformed to match the reference using the coefficients estimated from the set of corresponding points.

The idea is to extend the pairwise registration approach for a network of images. Networks require the extraction of image points from all image combinations, applying robust algorithms for outlier detection, and determining the registration parameters of the different images in a single step. This combined adjustment extends the traditional approach based on individual pairs and makes the analysis more exhaustive, incorporating all the data in a unique adjustment process. Indeed, in the case where points cannot be matched in some image pairs, the availability of more images allows one to exploit additional combinations. More images also provide higher data redundancy and therefore better results regarding precision and reliability.

Figure 1.
The scheme of the proposed network-based approach. Images points can have a direct connection to control points (triangles) or can be visible in multiple images (circles). Images without control points (image 2) can be registered by concatenating images into a network.

The extension from a pair of images to a network (Fig. 1) does not modify the equations used for registration. Equation (1) remains the mathematical model used. As all the images have to be registered with an affine transformation (with different numerical values of transformation parameters) between image points and control points (reference), we can assume that the geometric model between two generic images is still an affine transformation. This means that the same point can be matched in different images, requiring a unique “id” for the points and the computation of the corresponding multiplicity, i.e., the number of images where the same point is visible. Points visible in a single image (in this case corresponding control points are needed), but also in different configurations (such as pairs or triplets of images) play in a similar way of tie points in photogrammetric projects.

A generic block of $m$ images can be considered composed of ( $m^{2}-m)$ /2 combinations of pairs, which are analyzed independently for the identification of the correspondences, and they are then progressively combined. The *cost** for image registration can be quantified as O( $m^{2}-m)$ /2, in which a set of operations are independently repeated for the $n$ image points of a generic image pair (Snavely et al., 2007, 2008). Let us consider data processing carried out with scale invariant feature algorithms based on the detector/descriptor concept, such as the work proposed by Bay et al. (2008), Lowe (2004), Mikolajczyk and Schmid (2005), Mikolajczyk et al. (2005). In the case of photogrammetric applications carried out with a terrestrial camera, matching is carried out for all the images and comparing the descriptors associated with each point. Images have an important role during the detection phase, which can be quantified for the whole network has a linear cost that depends on the number of images O( $m$ ). The comparison through the detector can be then carried out by exploiting all image combinations and all point descriptor combinations, which could be properly combined in an overall cost written as O(( $m^{2}-m)(n^{2}-n)$ /4). The use of kd-tree-like algorithms (Beis & Lowe, 1997) for the comparison of points could improve data processing time because of the huge number of points expected from satellite images. However, the same idea is not here used because since the aim is to exploit all image combination, without losing precious information in a generic pair. The CPU cost can be therefore quantified as O(( $m^{2}-m)n$ log( $n)$ /2). On the other hand, satellite images are georeferenced, and the identification of corresponding points can be carried out by integrating point location as an additional constraint to limit the search in specific areas around the reference point. The comparison of descriptors becomes very rapid and efficient, resulting in short CPU time.

The algorithms for outlier rejection based on random sampling also require a small CPU cost, which is negligible if compared to the time for point detection in the images. Although robust estimators have to be applied to all image pairs, their contribution is relatively small when compared to the time needed to extract the point from all the images. Finally, the estimation of the solution via least squares is another operation which does not require a significant CPU cost, since optimized algorithms provide the solution vector without inverting the normal matrix.

The overall cost of data processing can be therefore computed as a linear function that depends only on the number of images, which is related to the detection phase of corresponding points thanks to the initial georeferencing parameters:

$\displaystyle F=\alpha\frac{1}{\textit{precision}}+\gamma\frac{1}{\textit{cost% }^{\ast}}=\alpha\frac{1}{\textit{precision}}+\gamma^{\prime}\frac{1}{m}=% \textit{max}$ (4)

This relationship becomes even more useful for large networks in which the number of images increases. In such case, the cost of point detection is more significant when compared to the other phases of data processing.

One may ask why the large cost of image resampling (after the computation of transformation parameters) is not taken into consideration. This is a linear cost that depends of the number of images. However, network resampling is not different than traditional resampling used in pairwise image registration with commercial software. The network-based approach can be used with more images than those necessary, and then resampling only the images that will be processed to obtain additional products. The comparison presented in this paper is carried out from image matching to parameter computation, where the network-based approach is different from traditional processing.
3. Datum definition for a network of satellite images

In the case of a network where all the images require registration, the lack of a reference system gives a rank deficiency during the computation of registration parameter. The observable (image coordinates) does not contain information about the reference system, and a set of fixed control points has to be included in the registration. Fixed points have to be included in the adjustment and can be (i) pixel values of an additional image assumed as a reference, (ii) an external set of pixel coordinates, or (ii) a set of map coordinates.

The datum problem (Zero Order Design, ZOD) involves the choice of an optimal reference system in which images will be registered (Fraser, 1982). The use of a common datum to establish the reference system is related to the need of georeferenced images in remote sensing applications. High, medium and low-resolution satellite images are (often) delivered with a reference system based on geographic or cartographic coordinates. Other solutions, such as free network adjustment, as described in Granshaw (1980), based on criteria for datum definition that improve parameter precision are not taken into account in this work to meet basic requirements of remote sensing applications.

Given a set of fixed coordinate (control points) corresponding to points visible in the images (at least a single image pair if the number and distribution of image-to-image matches guarantee a connection among all the images), Eq. (2) can be turned into a system of linear equations. Unknown parameters in the network are not only image transformation parameters. Additional unknowns are added to take into consideration points matched between two or more images, whose coordinates are not given in the control point list. Such additional unknowns correspond to the image-to-image coordinates projected in the control point reference system. Image-to-image points allow one to register images without direct visibility of control points and, at the same time, to run the adjustment process in a single step.

The linear functional and stochastic models based on Eq. (2) have the form:

$\displaystyle{\bm{A}}{\bm{x}}-{\bm{b}}={\bm{v}}$

(5) $\displaystyle{\bm{C}}_{bb}={\bm{I}}$

where ${\bm{A}}$ is the design matrix, ${\bm{x}}$ the unknown vector (including transformation parameters and ground coordinates of image-to-image points), ${\bm{b}}$ the observation vector, ${\bm{v}}$ the residual vector, and ${\bm{C}}_{bb}$ the variance-covariance matrix of observations. The identity matrix ${\bm{C}}_{bb}={\bm{I}}$ is here used because observable are image points usually matched with the same automatic operator.

As mentioned, the rank defect of ${\bm{A}}$ (i.e., no inverse for the singular normal equations) can be removed if at least three non-collinear ground points are included in the adjustment. The linear system of normal equations has the form:

$\displaystyle{\bm{N}}{\bm{x}}=\left[\begin{array}[]{cc}{\bm{A}}_{1}&\bar{\bm{N% }}\\ \bar{\bm{N}}^{T}&{\bm{A}}_{2}\\ \end{array}\right]\left(\begin{array}[]{c}x_{1}\\ x_{2}\end{array}\right)=\left(\begin{array}[]{c}n_{1}\\ 0\end{array}\right)={\bm{n}}$ (6)

where ${\bm{x}}_{1}$ contains the unknown transformation parameters ( $a_{j}$ , $b_{j}$ , $c_{j}$ , $d_{j}$ , $e_{j}$ , $f_{j})$ for the different images and ${\bm{x}}_{2}$ the image coordinates of the image-to-image point projected in the reference system provided by the fixed points. Redundancy $r$ depends on the number of images, i.e., six geometric parameters per every image, the number of control points, and the number of image points along with their multiplicity. Matrix ${\bm{N}}$ has a banded form, where ${\bm{A}}_{1}$ is a hyper-diagonal matrix with sub-matrices 6 $\times$ 6 corresponding to image transformation parameters, whereas ${\bm{A}}_{2}$ is made up of a sequence of 2 $\times$ 2 matrices for additional image-to-image points not matched in the list of control points.

The solution (and its precision ${\bm{C}}_{xx}$ is obtained with standard Least Squares techniques:

$\displaystyle\hat{\bm{x}}={\bm{N}}^{-1}{\bm{b}},\sigma_{0}^{2}=\frac{(N\hat{x}% -b)^{T}(N\hat{x}-b)}{r},{\bm{C}}_{xx}={\bm{\sigma}}^{2}_{0}{\bm{N}}^{-1}$ (7)

Equation (7) are used for a theoretical description of the mathematical problem. On the other hand, they are not the best choice for practical applications, especially when the number of images and points becomes relevant. Other methods (e.g., Gaussian decomposition) can be used to estimate $\hat{x}$ with a very limited CPU cost. The computation of $C_{xx}$ remains a time-consuming process.
4. The configuration problem for a network of satellite images

The problem (FOD, first-order design) involves the choice of the number of images ( $m$ ), image points ( $n$ ), control points ( $c$ ), and their distribution. The goal is to find a form for the design matrix ${\bm{A}}$ for which ${\bm{C}}_{xx}$ meets specific precision criteria. In this case, we assume that image points are matched with the same precision because matching is carried out with the same operator. Given a set of images and points, ${\bm{C}}_{xx}$ can be estimated with a simulation process based on the inversion of the normal matrix, assuming that $\sigma_{0}=$ 1 pixel and thus normalizing the effect of the operator for image matching. The aim is not to find better algorithms for image matching, but an optimal configuration for ${\bm{A}}$ based on images and points. In fact, ${\bm{A}}$ encapsulates the geometry of the whole network.

Different simulations are presented in the next sections, in which the parameters were varied to understand the role of the input variables. Simulations were carried out by using images of 10,000 pixels $\times$ 10,000 pixels, assuming that images were always acquired over the same area, such as in the case of change detection from a satellite time series.

The particular symmetric form of Eq. (2) provides the same result for the precision of pixel coordinates ( $\sigma_{x}=\sigma_{y}$ ). This consideration is valid also for the transformation parameters of images ( $a_{j}$ , $b_{j}$ , $c_{j}$ , $d_{j}$ , $e_{j}$ , $f_{j}$ ), for which $\sigma_{cj}=\sigma_{fj}$ and $\sigma_{aj}=\sigma_{bj}=\sigma_{dj}=\sigma_{ej}$ . As the aim of the registration is to find the parameters of the transformation between the different images, the criteria used to compare different simulations is not based on point coordinates (like in geodetic, leveling, GPS and photogrammetric networks) but on transformation parameter precision. This is a remarkable difference concerning other types of the design process: we are interested in image parameters rather than points. Results are illustrated for translation ( $c$ , $f$ ) and scale ( $a$ , $b$ , $d$ , $e$ ) with only two graphs thanks to the symmetric form of Eq. (2).

4.1 The influence of image points on network precision ( $n=$ variable, $m=$ 2, $c=n$ )

Let us consider the simplest case of a single image pair ( $m=$ 2). We assume that all image points are at the same time control points ( $m=c$ ). The number of equations is $m\cdot n=$ 2 $\cdot n$ , and the number of unknowns is 6 $\cdot m$ , which correspond to the 12 affine transformation parameters. The distribution of the used image points (and control points is such case) is a regular grid obtained by dividing the image into rectangular cells. This work aims to run simulations with a variable number of image points and check the precision of image transformation parameters. The chosen grid size was varied from 2,000 $\times$ 2,000 pixels (25 points per image) to 50 $\times$ 50 pixels (40,000 points per image), using some intermediate steps. The size of normal matrix ${\bm{N}}$ is 12 $\times$ 12 and does not depend on the number of points because image points are at the same time control points.

The form of Eq. (2) separates the coefficients for parameters (the two blocks ( $a_{j}$ , $b_{j}$ , $c_{j}$ ) and ( $d_{j}$ , $e_{j}$ , $f_{j}$ ), for which the 6 $\times$ 6 part of each image is split into two additional 3 $\times$ 3 matrices. This means that the two blocks ( $a_{j}$ , $b_{j}$ , $c_{j}$ ) and ( $d_{j}$ , $e_{j}$ , $f_{j}$ ) are uncorrelated. Moreover, the parameter of different images are uncorrelated.

The achieved precision of transformation parameters depends on $\sqrt{n}$ ., i.e., increasing the number of points $n$ improves precision of image parameters, whereas the number of unknowns remain unvaried. A huge number of points does not provide significant improvement of parameter precision. A precision on translation parameter of about 1/7 pixels is already reached with 400 points, whereas a value of about 1/20 pixels is obtained with 2500 points. This is an important result: although matching operators can extract a large number of points from large satellite images (in terms of number of pixels), points should be decimated to reach a uniform distribution in the images, without exceeding a regular grid of about 200 $\times$ 200 cells (i.e., 2500 points), from which the improvement is minimal.

4.2 The importance of control points: network analysis reduced to pairwise processing ( $n=$ constant, $m=$ variable, $c=n$ )

As mentioned in the previous section, a good choice regarding points is between 400 and 2500 points. In this example, $n=$ 400 was used. Data processing was carried out by varying the number of images, from 2 to 32 images. Control point configuration was set to match the number of image points. The size of the normal matrix still depends on the number of images only ( $6m\times 6m$ ) since all image points are at the same time control points. Considerations similar to those reported in the previous section can be extended in such cases: transformation parameters for $x$ and $y$ are uncorrelated, and parameters of different images are uncorrelated.

The precision of translation and scale parameters in the case of image processing with a variable number of images remains unvaried. The number of images does not provide any change in translation and scale parameters ( $\sigma=$ cost.). This means that the independent analysis of image pairs, i.e. a more traditional matching approach based on reference-to-sensed matching approach with a CPU cost $O(m)$ , or the proposed exhaustive network-based approach with a cost $O((m^{2}-m)/2)$ , provide the same result. This outcome could be apparently intended as a failure of the network-based approach. Why use a larger number of images if the precision is independent of the number of images? The first significant consideration is related to the achievable matching results with real images. As the number of image points is equal to the number of control points, we are in the case where automated matching provided the best result achievable. This means not only that we are working with an ideal algorithm for matching, but also that images do not have any disturbing elements (such as clouds or shadows) that could generate the loss of some control points. This is quite unrealistic when real data are used. Also, we are working with a complete set of control points, i.e., external constraints that are already available for all the images that require registration. This is also an entirely unrealistic situation.

4.3 Network analysis to reduce control points ( $m=$ variable, $n=$ const., $c=$ variable)

Let us consider a constant number of points $n=$ 2500 visible in every image. Starting from a set of control points, the aim is to check the precision of translation and scale parameters for different image configuration. The test was carried out by increasing the number of images $m$ (2, 4, 8, 16, 24), and reducing the number of control points with a decimation factor $d$ (1, 2, 3, 4, 5, 10, 20, 30, 40), which is a function of the number of image points $c=n$ / $d$ . The case $d=$ 1 provides a dataset in which all image points are also control points. The case $d>$ 1 introduces additional unknowns in the network-based processing, which corresponds of $x$ and $y$ coordinates of image-to-image points without connection to control points, which are transformed in the reference system of control points.

Figure 2.

Number of equations and unknowns for 2 and 24 images as a function of the decimation factor.

The choice of a constant value for $n$ is motived by the results obtained in Section 4.1, for which the improvement of precision is no longer significant with more than 2,500 image points. Here, it was decided to use $n=$ 2,500 (instead of 400 like in the previous section) to allow a significant decimation of control points, since $c$ depends on $n$ in the formulation used for simulation.

The number of equations and unknowns (i.e., rows and columns of matrix ${\bm{A}}$ ) can be estimated as:

$\displaystyle\textit{equations}=2mn$ (8) $\displaystyle\textit{uknowns}=6m+\left({n-\frac{n}{d}}\right)$

The previous equations highlight that the number of equations is directly proportional to the number of images. Adding other images will immediately result in a substantial improvement of data redundancy because the number of unknowns is only increased by 6 $m$ , whereas the number of equations by 2 $m n$ . Here, it is expected that $n>>m$ . The most significant parameter in the number of unknowns is related to the number of image points that do not have corresponding control points. Decimation of the control points of a factor two leads to a significant increase in the number of unknowns. Figure 2 shows the results for a network of 2 and 40 images. The increment regarding unknowns is quite limited, whereas the advantage in the number of equations is significant. Moreover, it does not depend on control points, resulting in a much better redundancy. It is also interesting that the number of unknowns rapidly changes with decimation of factor 2, and it becomes less significant for large decimation factors. It is also interesting to understand the form of normal matrix ${\bm{N}}$ for the case of 2 images with a decimation 2, which is made up of a 12 $\times$ 12 block corresponding to image parameters. The remaining diagonal part corresponds to the unknown coordinate of image-to-image points projected in the reference system used for image registration. Vertical and horizontal symmetrical bands complete the system.

Finally, the complete result with the simulation based both on a variable number of images and decimation factors are shown in Fig. 3. Results confirm that without control point decimation the number of images is not significant, whereas an analysis based on a network of images provide better results, with an improvement of a factor 2.7 between the case with 2 and 24 images, with control points decimated by 40 times. The second important outcome is the considerable variability between the solution with the full set of control points and the decimated dataset (40 times), for the case of 2 images. Results regarding precision are much worse with an estimated factor 4.1. In the case of an analysis based on 24 images, the degradation of precision corresponds to a factor 1.5. Increasing the number of images will provide a set of curves in Fig. 3, which tend to become horizontal lines, confirming a very limited degradation of precision notwithstanding the strong decimation factor. According to these considerations, we can state that increasing the number of images allows one to reduce the number of control points, giving the opportunity to focus on the images (now available in cloud-based services) rather than control points.

Figure 3.

The relationship between transformation parameter precision in the case of a variable number of images and control points. A similar graph is obtained for scale parameters.

4.4 Considerations and outcomes

The experiments in the previous section demonstrated that a huge number of image points is not significant, as precision improvement depends on $\sqrt{n}$ . This means that after extracting image points, they can be decimated to not exceed 2,500 points per image, reducing the size of design matrix ${\bm{A}}$ s.

We could say that the results regarding precision for a network of satellite images requiring registration, with a limited number of control points and several images, become equivalent to those of the same network with the full set of control points. At the same time, this allows one to quantify the precision with several images and a limited dataset of control points as those achievable with the full set of control points, leading to the following relationships for the precision:

$\displaystyle\sigma_{cj}=\sigma_{fj}=\sqrt{\frac{\left({\mathop{\sum}\nolimits% _{i}x_{i}^{2}}\right)^{2}-\left({\mathop{\sum}\nolimits_{i}x_{i}y_{i}}\right)^% {2}}{n\left[{\left({\mathop{\sum}\nolimits_{i}x_{i}^{2}}\right)^{2}-\left({% \mathop{\sum}\nolimits_{i}x_{i}y_{i}}\right)^{2}}\right]+2\left({\mathop{\sum}% \nolimits_{i}x_{i}}\right)^{2}\left({\mathop{\sum}\nolimits_{i}x_{i}y_{i}-% \mathop{\sum}\nolimits_{i}x_{i}^{2}}\right)}}$ (9) $\displaystyle\sigma_{aj}=\sigma_{bj}=\sigma_{dj}=\sigma_{ej}=\sqrt{\frac{n% \mathop{\sum}\nolimits_{i}x_{i}^{2}-\left({\mathop{\sum}\nolimits_{i}x_{i}y_{i% }}\right)^{2}}{n\left[{\left({\mathop{\sum}\nolimits_{i}x_{i}^{2}}\right)^{2}-% \left({\mathop{\sum}\nolimits_{i}x_{i}y_{i}}\right)^{2}}\right]+2\left({% \mathop{\sum}\nolimits_{i}x_{i}}\right)^{2}\left({\mathop{\sum}\nolimits_{i}x_% {i}y_{i}-\mathop{\sum}\nolimits_{i}x_{i}^{2}}\right)}}$

in which $\sum_{i}x_{i}^{2}=\sum_{i}y_{i}^{2}$ and $\sum_{i}x_{i}=\sum_{i}y_{i}$ for a simulated dataset in which a regular grid is used. Both grid and image have to be squared.

The precision of translation parameters provides an immediate evaluation of the overall accuracy achievable from a network with a large number of images. A graphic visualization of $\sigma_{cj}=\sigma_{fj}$ as a function of the number of image points $n$ is shown in Fig. 4, which reveals that precision is proportional to $\sqrt{n}$ , as already mentioned.

Figure 4.

The precision for a network of several images as a function of image points. When the number of images increases, the effect of control points on parameter precision becomes less significant.

The final cost function introduced in Section 2 for a large network of images becomes:

$\displaystyle F=\alpha\frac{1}{\textit{precision}}+\gamma\frac{1}{cost^{\ast}}% =\alpha^{\prime}\sqrt{n}+\gamma^{\prime}\frac{1}{m}=\textit{max}\ \text{(for % large $m$)}$ (10)

The relationships illustrate that the cost can be minimized by reducing the number of images. On the other hand, this is not feasible to obtain a sufficient precision with limited control points, as demonstrated in the previous sections. At the same time, precision can be improved by increasing the number of image points, and it does not significantly affect CPU time. A network of multiple images connected provides better registration results for images with a weak connection to the final reference system (e.g., few or even no image-to-control points) because the additional images can add new connections inside the network. Such result can also motivate the use of the large number of satellite images available in web-services, selecting more images than those strictly necessary for the following phases of the work.

5. Conclusion

This paper presented some considerations on the advantages and disadvantages of image network analysis applied to the registration of satellite images. Results achieved through simulation based on the variance-covariance matrix of least squares adjustment revealed that an increment of the number of images becomes significant when a complete dataset of control points is not available. The number of points matched is not significant after gathering a regular distribution with only 400 points.

Results regarding precision for a network of satellite images, in which images have to be registered with a limited number of control points, become more similar to those of the same network with the full set of control points. From the point of view of automation in satellite image registration, this is a significant point since most work can be directly carried out with the images, i.e., the data available in web-services. We could say that the network will be connected not only with control points (such as in independent processing image pairs), but additional image-to-image connections will balance the lack of a large dataset of control points.

The proposed approach was presented for the particular case of registration problems. On the other hand, it could be extended for other practical issues where image redundancy can improve automatic data processing. An example is the case of the mathematical formulation used for atmospheric normalization via corresponding image points, in which the extension from independent image pairs to a network of images can be carried out by introducing a mathematical model that encapsulates all the images. Future work will be carried out to extend the proposed approach to other applications.

References

Alonso,

Muñoz-Carpena,

Kennedy,

R.E.

& Murcia,

(2016). Wetland landscape spatio-temporal degradation dynamics using the new google earth engine cloud-based platform. Opportunities for non-specialists in remote sensing Transactions of the ASABE, 59(5), 1333-1344.

Barazzetti,

Gianinetto,

& Scaioni,

(2014). Automatic co-registration of satellite time series via Least Squares adjustment. European Journal of Remote Sensing, 47, 55-74.

Bay,

Ess,

Tuytelaars,

& Van Gool,

(2008). SURF: Speeded up robust features. Computer Vision Image Understanding, 110(3), 346-359.

Beis,

& Lowe,

, (1997). Shape indexing using approximate nearest neighbour search in high dimensional spaces. Proceedings of CVPR’07, June 1997, Puerto Rico (USA), 1000-1006.

Bouchiha,

& Besbes,

(2013). Automatic Remote-sensing Image Registration Using SURF. International Journal of Computer Theory and Engineering, 5(1), 88-92.

Brown,

G. L.

(1992). A survey of image Registration Techniques. ACM Computing Surveys, 24, 325-376.

Castro,

E.D.

& Morandi,

(1987). Registration of translated and rotated images using finite Fourier transform. IEEE Transactions on Pattern Analysis and Machine Intelligence, 9, 700-703.

Christmann,

, (1996). High Breakdown Point Estimators in Logistic Regression. IN Robust Statistics, Data Analysis, and Computer Intensive Methods, springer New York, 79-89.

Dawn,

Saxena,

& Sharma,

(2010). Remote Sensing Image Registration Techniques: A Survey. International Conference on Image and Signal Processing, ICISP 2010, 103-112.

10.

Eshagh,

& Kiamehr,

(2013). A strategy for optimum designing of the geodetic networks from the cost, reliability and precision views. Acta Geodaetica et Geophysica Hungarica, 3, 297-308.

11.

Fischler,

M.A.

& Bolles,

R.C.

(1981). Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM, 24(6), 381-395.

12.

Fraser,

C.S.

(1984). Network design considerations for non-topographic photogrammetry. Photogrammetric Engineering and Remote Sensing, 50(9), 1115-1126.

13.

Goldblatt,

You,

Hanson,

& Khandelwal,

A.K.

(2016). Detecting the boundaries of urban areas in India: A dataset for pixel-based image classification in google earth engine. Remote Sensing, 8(8), art. no. 634, 28 pages.

14.

Goshtasby,

Stockman,

G.C.

& Page,

C.V.

(1986). A region-based approach to digital image registration with subpixel accuracy. IEEE Transactions on Geoscience and Remote Sensing, 24, 390-399.

15.

Gosthtasby,

(2005). 2-D and 3-D image registration. John Wiley & Sons, Hoboken, NJ, U.S.A., 258 pages.

16.

Grafarend,

& Sansò,

(1985). Optimization and Design of Geodetic Networks. Springer-Verlag Berlin Heidelberg Publisher, 606 pages.

17.

Granshaw,

S.I.

(1980). Bundle adjustment methods in engineering photogrammetry. The Photogrammetric Record, 10(56), 181-207.

18.

Gruen,

(2012). Development and Status of Image Matching in Photogrammetry. The Photogrammetric Record, 27(137), 36-57.

19.

Horowitz,

F.G.

(2015). MODIS Daily Land Surface Temperature Estimates in Google Earth Engine as an Aid in Geothermal Energy Siting, in: Proceedings World Geothermal Congress 2015, Melbourne, Australia, 19-25 April 2015, 4 pages.

20.

Le Moigne,

Netanyahu,

N.S.

& Eastman,

R.D.

(2011). Image Registration for Remote Sensing. Cambridge University Press, U.K., 484 pages.

21.

Lowe,

(2004). Distinctive Image Features from Scale-Invariant Keypoints. International Journal of Computer Vision, 60(2): 91-110.

22.

Ma,

Wu,

Wang,

Huang,

Ranjan,

Zomaya,

& Jie,

(2015). Remote sensing big data computing: Challenges and opportunities, Future Generation Computer Systems, Volume 51, October 2015, pages 47-60.

23.

Mikolajczyk,

Tuytelaars,

Schmid,

Zisserman,

Matas,

Schaffalitzky,

Kadir,

& Van Gool,

(2005). A Comparison of Affine Region Detectors. International Journal of Computer Vision, 65(1/2), 43-72.

24.

Mikolajczyk,

& Schmid,

(2005). A performance evaluation of local descriptors. IEEE Transactions on Pattern Analysis and Machine Intelligence, 10(27), 1615-1630.

25.

Nativi,

Mazzetti,

Santoro,

Papeschi,

Craglia,

& Ochiai,

(2015). Big Data challenges in building the Global Earth Observation System of Systems. Environmental Modelling & Software, (68), 1-26.

26.

Nocerino,

Menna,

& Remondino,

(2014). Accuracy of typical photogrammetric networks in cultural heritage 3D modelling projects. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci, XL-5, 465-472.

27.

Ogundare,

J.O.

(2015). Survey Design and Analysis. In Precision Surveying, John Wiley & Sons, Inc, 209-236.

28.

Poli,

& Toutin,

(2012). Review of developments in geometric modelling for high resolution satellite pushbroom sensors. The Photogrammetric Record, 27(137), 58-73.

29.

Pluim,

J.P.W.

Maintz,

J.B.A.

& Viergever,

M.A.

(2001). Mutual information matching in multiresolution contexts. Image and Vision Computing, 19, 45-52.

30.

Price,

K.E.

(1985). Relaxation matching techniques – a comparison. IEEE Transactions on Pattern Analysis and Machine Intelligence, 7, 617-623.

31.

Pulighe,

& Lupia,

(2016). Mapping spatial patterns of urban agriculture in Rome (Italy) using Google Earth and web-mapping services. Land Use Policy, 59, 49-58.

32.

Qiaoliang,

Wang,

Liu,

& Chen,

(2009). Robust Scale-Invariant Feature Matching for Remote Sensing Image Registration. IEEE Geoscience and Remote Sensing Letters 6(2), 287-291.

33.

Snavely,

Seitz,

S.M.

& Szeliski,

(2007). Modeling the World from Internet Photo Collections. International Journal of Computer Vision, 80(2), 189-210.

34.

Teke,

& Temizel,

(2010). Multi-spectral Satellite Image Registration Using Scale-Restricted SURF. 20th International Conference on Pattern Recognition, Istanbul, 2310-2313.

35.

Townshend,

J.R.

Justice,

C.O.

Gurney,

& McManus,

(1992). The impact of misregistration on change detection. IEEE Transactions on Geosciences and Remote Sensing, 30(5), 1054-1060.

36.

Warren,

M.S.

Brumby,

S.P.

Skillman,

S.W.

Kelton,

Wohlberg,

. Mathis,

Chartrand,

Keisler,

& Johnson,

(2015). Seeing the Earth in the Cloud: Processing one petabyte of satellite imagery in one day.

37.

IEEE Applied Imagery Pattern Recognition Workshop, AIPR 2015, art. no. 7444536

38.

Zhao,

Liu,

Du,

& Li,

(2009). Feature-Based Geometric Registration of High Spatial Resolution Satellite Imagery. Joint Urban Remote Sensing Event, 1-5.

39.

Zhu,

Sun,

Liao,

Yang,

& Zhu,

(2009). Earth System Scientific Data Sharing Research and Practice. Journal of Geo-Information Science 28, 17-20.