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

HAL
open science

# Natural Color Satellite Image Mosaicking Using Quadratic Programming in Decorrelated Color Space

Rémi Cresson, Nathalie Saint-Geours

► To cite this version:

Rémi Cresson, Nathalie Saint-Geours. Natural Color Satellite Image Mosaicking Using Quadratic Programming in Decorrelated Color Space. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 10.1109/JSTARS.2015.2449233. [hal-01373314](https://hal.science/hal-01373314v1)

HAL Id: hal-01373314

https://hal.science/hal-01373314v1

Submitted on 28 Sep 2016

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_2.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Natural color satellite image mosaicking using quadratic programming in decorrelated color space

Rémi Cresson and Nathalie Saint-Geours*

Generating mosaics of orthorectified remote sensing images is a challenging task because of the colorimetric differences between adjacent images introduced by land use, surface illumination, atmospheric conditions, and sensor. Most of the existing color correction methods involve pairwise techniques, which are limited when the collection of images is large with numerous overlaps. Besides, available techniques do not operate in a color space suited for true-color processing. This paper presents a simple and robust method to perform the global colorimetric harmonization of multiple overlapping remote sensing images in natural colors (RGB). Our parameter-free method deals simultaneously with any number of images, with any spatial layout, and without any single reference image. It is based on the resolution of a quadratic programming optimization problem. It operates in the $l_αβ$ decorrelated color space, which is well suited for human vision of natural scenes. The results obtained from the mosaicking of 132 RapidEye color orthoimages over mainland France demonstrate good potential for performing colorimetric harmonization automatically and effectively.

## I. INTRODUCTION

There is a growing need in Earth Sciences for high resolution (HR) to very high resolution (VHR) satellite data (decametric to metric) covering extensive areas like regions, countries or even whole continents [1]. In such applications, remote sensor footprints are usually small compared to the extensive area to be covered: satellite image mosaic generation is thus required to combine a large number of small-footprints scenes into a single big mosaic [2]. In most cases, the numerous cloud-free scenes to be combined in a mosaic have been acquired over a broad time window, due to the low orbital revisit cycles and the probability of cloud cover. Therefore, these images most often exhibit a diversity of atmospheric and seasonal conditions, making image mosaic generation a challenging task.

Image mosaicking involves several steps: i) geometric correction of images (registration), ii) seamline detection in regions where adjacent images are overlapping [3–6], iii) color harmonization of images, and iv) composition of the final mosaic, which may involve techniques such as image blending or image stitching [7–9], which ensures that the transition between two neighbouring images along their seamline will be soft enough. The issue of image registration has long been debated by the remote sensing community—see [10] for a review—and many tools now allow for a fully automatic geometric correction of images. Color harmonization is a more arduous task, and in spite of its importance in the image mosaicking process, it is far from being as mastered as geometric correction. The goal of color harmonization is to reduce color differences between adjacent images in the mosaic. Indeed,

the radiometry of a scene depends on land use, surface illumination, atmospheric conditions, and sensor: all these factors may vary from one image to another, and may result in radiometric heterogeneity in the mosaic. Adjustment techniques are thus required to compensate for these radiometric differences.

Regarding this issue, two different purposes of remote sensing image mosaicking must be distinguished. The first use of satellite image mosaics is quantitative; in this case, the analyst will further process the image mosaic to compute thematic indices, such as vegetation or built-up indices, or other quantitative products. It is then of great importance that the mosaic preserves a number of radiometric signal characteristics. The main goal is thus to ensure the radiometric fidelity of the mosaic: this is what *radiometric adjustment* or *radiometric normalization* techniques are used for. The second purpose of image mosaicking is to produce good-looking mosaics for illustrative use only. Such image mosaics should be visually appealing, seamless, and should hide all the technical issues related to the acquisition of satellite imagery. Their quality can mainly be assessed through visual inspection—some quantitative measures (color and structure similarity measures) have also been proposed for that purpose [11]. In this case, we can simply speak of *color harmonization* techniques. While most radiometric normalization techniques can also be used as simple color harmonization techniques for illustrative use only, the reverse is generally not true.

In this paper, we focus on a new color harmonization technique to generate satellite image mosaics for illustrative use only. Color harmonization is usually achieved using a variety of empirical techniques. Mean-variance equalization consists in correcting one image with respect to a reference image so that both images have the same mean and the same variance in their overlapping part. Histogram matching [12] goes one step further: the corrected image will have the same complete pixel value distribution as the reference image. Relative radiometric normalization techniques [13–15] and radiometric normalization tech-

<footnote>* R. Cresson and N. Saint-Geours are with the Irstea, UMR TETIS, BP 5095, Montpellier, France (e-mail: remi.cresson@teledetection.fr; nathalie.saint-geours@teledetection.fr).; Manuscript submitted November, 21st, 2014; Revised manuscript submitted May, 11th, 2015.</footnote>
niques based on invariant features [16] usually involve
the selection of a sample of pixels in pairs of images,
and the computation of correction models based on
some statistical regression.

Nevertheless, up to our knowledge, most of the
available methods for color harmonization have a few
drawbacks. First, existing works most often present
color harmonization techniques for pairs of images
only: one image is corrected to match a reference im-
age. To generate a mosaic with *n* > 2 images, one usu-
ally proceed sequentially : a first image is corrected
with respect to some reference image, and the two
images are merged into a first mosaic. Then, a third
image is corrected to match the color characteristics
of the current mosaic, and a 3-images mosaic is assem-
bled, and so on. Colorimetric correction is thus prop-
agated from the first reference image to the whole mo-
saic. However, depending on the spatial layout of the
*n* images and their overlapping areas, several paths
may be taken to build the whole mosaic. [17] shows
how different results may be obtained when changing
this sequential path, and discuss the accumulation of
errors associated with intermediary scene corrections.
Pairwise techniques are thus limited when the collec-
tion of satellite images is large, with large amount of
image overlaps.

Besides, the available color harmonization tech-
niques operate in the RGB color space, or directly on
the original Digital Numbers from spectral bands of
satellite data. However, it is well known in the com-
puter vision community that other color spaces are
more suitable for color processing, but these develop-
ments have not reached the remote sensing community
yet. Regarding this subject, [18] recently provided a
extensive review of numerous *color transfer* or *colour
mapping* techniques, which aim to recolor a given im-
age by deriving a mapping between that image and
another image serving as a reference.

The goal of this paper is thus to present a simple
and robust method for color harmonization of remote
sensing images in natural colors, to produce large im-
age mosaics for illustrative purposes. The method
works (i) simultaneously (rather than sequentially) for
any number of images, (ii) for any spatial layout of
the images and their overlapping parts, (iii) without
any single reference image, (iv) in a decorrelated color
space ($l\alpha\beta$) which is well suited for human vision of
natural scenes. It is limited to natural color images.
Besides, it is parameter-free and fully automated. We
first give in Section II a detailed description of the
mosaicking process. Our approach is then successfully
tested on a collection of 132 RapidEye satellite images
over mainland France supplied by GEOSUD project
(Section III and IV). We finally discuss the main ad-
vantages and limits of our approach in Section V.

II. METHOD

A. Overview

We consider a collection of *n* partly overlapping
satellite images, geometrically corrected, in natural
colors (RGB). Our goal is to build a single mosaic
from these images for illustration purposes— i.e., ra-
diometric fidelity is not sought after. We propose to
build on the *color transfer* technique developed by [19]
to define linear color correction models for each image,
in a decorrelated color space. Besides, while [19] per-
formed mean-variance equalization on pairs of images,
we formulate the color correction problem for mosaic
generation simultaneously from *n* images as a multi-
variate optimization problem, namely a quadratic pro-
gramming (QP) problem for which standard solvers
can be used.

Our method can be decomposed into four steps
(Fig. 1):

1. In the first step, we convert the RGB values of the original images to a decorrelated color space named $l\alpha\beta$;

2. Then, we apply a linear color correction model to each image; decorrelation lets us treat the three color channels *l*,*α* and *β* separately, simplifying the method. We search for the optimal parameters of linear color correction models to minimize a cost function. Two different cost functions are tested: one is based on the comparison of both the mean and standard deviation of neighboring pair of images in their overlapping areas, the other is based on the root mean squared error (RMSE) between pair of images in their overlapping areas;

3. We then proceed to mosaic composition. This step is not the focus of our work, thus we chose to use a simple and fast compositing technique which copies the last image over earlier ones in areas of overlap. Other more advanced methods could have been used to compose a seamless mosaic, such as the multi-band blending strategy [7] or a variational approach [9];

4. Finally, we transform the output mosaic image back to the RGB color space.

In the following subsections, we give more details on
each step of the process. We first display the conver-
sion equations from RGB to *l*αβ color space in II B.
We then describe the color correction models in II C,
and define the cost functions to be minimized in II D.
Finally, the formulation of the QP problem and its
solution are detailed in II E.

B. Conversion of RGB values to *l*αβ color space

The first step of color harmonization is to convert
RGB values into a new color space called *l*αβ. This
<figure><img src="image_3.png" /><figcaption>FIG. 1: Flowchart.</figcaption></figure>

color space was first introduced by [20], based on data-driven human perception research that assumes the human visual system is ideally suited for processing natural scenes. In the $l\alpha\beta$ color space, the $l$ axis represents an achromatic channel, while the $\alpha$ and $\beta$ channels are chromatic yellow-blue and red-green opponent channels. The main feature of the $l\alpha\beta$ color space is that, for natural scenes, the correlations between the three channels' values are low. [19] used this color space to transfer colors from a reference image to a target image, performing a simple mean-variance equalization in the $l\alpha\beta$ space — this approach was further developed by [21] and [22]. This straightforward technique offers great results when the composition of reference and target images is similar—e.g. same proportion of different land covers in both images.

As explained in the introduction, we focus here on satellite images with natural colors. As the composition of two overlapping images is usually very similar, the color transfer technique developed by [19] seems to be suited to our needs. We thus perform RGB-to-$l\alpha\beta$ conversion on each image, resulting in the computation of three new channels $l$, $\alpha$, and $\beta$, in which inter-channel correlations are expected to be low. The conversion of RGB values to $l\alpha\beta$ values follows three steps—see [19] for more details:

* first, RGB values are transformed to LMS cone space values (Long, Middle, Short wavelength cone), that describe human eye spectral sensitivities [Eqn. (1)] (Matrix values in Eqn. (1) slightly differ from those given by [19, Eqn (2)]), which were found to be erroneous by an anonymous reviewer. Correct values were directly computed from [19, Eqn (1)].);

$$ \begin{bmatrix} L \\ M \\ S \end{bmatrix} = \begin{bmatrix} 0.3811 & 0.5783 & 0.0406 \\ 0.1967 & 0.7244 & 0.0790 \\ 0.0241 & 0.1228 & 0.8531 \end{bmatrix} \cdot \begin{bmatrix} R \\ G \\ B \end{bmatrix} \quad (1) $$

* values in the LMS space are most often highly skewed; we partly eliminate this skew by converting values to logarithmic scale [Eqn. (2)];

$$ \begin{bmatrix} L_{log} \\ M_{log} \\ S_{log} \end{bmatrix} = \begin{bmatrix} \log L \\ \log M \\ \log S \end{bmatrix} \quad (2) $$

* finally, a linear transformation introduced by [20] is applied to convert the signal from LMS values to the decorrelated $l\alpha\beta$ color space [Eqn. (3)].

$$ \begin{bmatrix} l \\ \alpha \\ \beta \end{bmatrix} = \begin{bmatrix} \frac{1}{\sqrt{3}} & 0 & 0 \\ 0 & \frac{1}{\sqrt{6}} & 0 \\ 0 & 0 & \frac{1}{\sqrt{2}} \end{bmatrix} \cdot \begin{bmatrix} 1 & 1 & 1 \\ 1 & 1 & -2 \\ 1 & -1 & 0 \end{bmatrix} \cdot \begin{bmatrix} L_{log} \\ M_{log} \\ S_{log} \end{bmatrix} \quad (3) $$

Equations for the inverse transformation from $l\alpha\beta$ space to RGB space can be derived from Eqn. (1) to (3). They are needed in the last step of satellite image mosaicking, to move back from the $l\alpha\beta$ color space to the original RGB color space.

### C. Parametric model for color correction

The second step of color harmonization for image mosaicking is to define color correction models for each of the $n$ images to be included in the mosaic. Because the correlation of channels' values in the $l\alpha\beta$ color space is low, we can use independent correction models for each of the three axes $l$, $\alpha$ and $\beta$. In this study we considered two different linear color correction models: a general linear transformation [Eqn. (4a)] and a linear transformation with zero y-intercept [Eqn. (4b)]:

$$ v' = a \cdot v + b \quad (4a) $$

$$ v' = a \cdot v \quad (4b) $$

where $a$ and $b$ are real-valued parameters, $v$ (input value) is the $l$, $\alpha$ or $\beta$-value of a pixel in the original image, and $v'$ is the $l$, $\alpha$ or $\beta$-value of the same pixel after transformation. A different value of both parameters $a$ and $b$ must be set for each channel and for each satellite image to be included in the mosaic. If $n$ denotes the number of images, the total number of parameters is $6n$ for the general case of Eqn. (4a), and only $3n$ for Eqn. (4b).

### D. Cost functions

In order to select the optimal parameters of color correction models defined in Eqn. (4), we first need to define some criteria to assess the quality of a mosaic with respect to color harmonization. A usual quality
<figure><img src="image_4.png" /><figcaption>FIG. 2: Overlapping area between two images i and j.</figcaption></figure>

metric for image mosaics is based on the comparison of the histograms of neighboring images in their overlapping parts: the closer the histograms are, the better the mosaic is. In this study, we consider two alternative cost functions to assess the quality of a mosaic. A first cost function is based on two summary statistics derived from these histograms: mean and standard deviation of neighboring images in their overlapping parts. A second cost function is based on the root mean squared error (RMSE) between two images in their overlapping part.

Independent cost functions are computed for each channel (l, α, β). From now on, without loss of generality, we will focus on a single channel (l, α, or β). Let us consider a pair (i, j) of overlapping images (Fig. 2). The area of their overlapping part is $S_{ij}$. We denote with $\mu_{ij}$ and $\sigma_{ij}$ (resp. and $\mu_{ji}$ and $\sigma_{ji}$) the mean and standard deviation of image i (resp. j) in the overlapping part [Eqn. (5)]. We will also denote $r_{ij}$ the RMSE of images i and j in their overlapping part.

$$ \mu_{ij} = E(v_i) \quad \text{and} \quad \sigma_{ij} = \sqrt{E\left[(v_i - \mu_{ij})^2\right]} \quad \text{and} \quad r_{ij} = \sqrt{\sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} \cdot (r'_{ij})^2} \quad (5) $$

with $v_i$ (resp. $v_j$) denoting pixel values in image i (resp. j), and E denoting the mean value over all pixels of the overlapping part between images i and j. After linear color correction model is applied to each image, the mean and standard deviation are transformed to $\mu'_{ij}$ and $\sigma'_{ij}$ (resp. $\mu'_{ji}$ and $\sigma'_{ji}$), and the RMSE is transformed to $r'_{ij}$.

### 1. Mean-and-variance-based cost function

An elementary cost $c_{ij}$ is defined as the half squared-difference between means of the two neighboring images (i, j) after transformation, weighted by the overlapping area:

$$ c_{ij} = \frac{1}{2} \cdot S_{ij} \cdot (\mu'_{ij} - \mu'_{ji})^2 \quad (6) $$

For a set of n images, the total cost related to mean-equalization is equal to the sum of pairwise costs:

$$ C^{\mu} = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=i}^{n} S_{ij} \cdot (\mu'_{ij} - \mu'_{ji})^2 \quad (7) $$

In an analogous way, the total cost related to variance-equalization is equal to the sum of pairwise

half squared-differences between standard deviation of two neighboring images i and j after transformation, weighted by the overlapping area:

$$ C^{\mu\sigma} = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=i}^{n} S_{ij} \cdot (\sigma'_{ij} - \sigma'_{ji})^2 \quad (8) $$

The first cost function $C^{\mu\sigma}$ is obtained by summing mean-equalization cost and variance-equalization cost:

$$ C^{\mu\sigma} = C^{\mu} + C^{\sigma} \quad (9) $$

### 2. RMSE-based cost function

The second cost function we consider is related to the minimization of RMSE criterion. It is equal to the sum of squared pairwise RMSE of neighboring images i and j after transformation, weighted by their overlapping area:

$$ C^r = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=i}^{n} S_{ij} \cdot (r'_{ij})^2 \quad (10) $$

## E. Minimization of cost functions

The image mosaic generation can now be formulated as a multivariate optimization problem: our goal is to find an optimal set of parameters for color correction models that minimizes the selected cost function, i.e. either the total mean-and-variance-based cost function $C^{\mu\sigma}$, or the RMSE-based cost function $C^r$. We first rewrite cost functions $C^{\mu\sigma}$ and $C^r$ as quadratic functions of the parameters a and b of the color correction models, then we set a number of linear constraints on these parameters, and we finally detail how to find the optimal set of parameters that minimize the cost functions under the specified constraints.

### 1. Canonic formulation of cost functions

Let $\mathbf{x}$ denote the vector of parameters of the n color correction models: $\mathbf{x} = (a_1, \dots, a_n, b_1, \dots, b_n)$ with $a_i$ and $b_i$ the slope and intercept coefficients related to image i. Cost functions $C^{\mu\sigma}$ and $C^r$ can be written as quadratic functions of parameter vector $\mathbf{x}$ (Appendix VII for detailed calculations):

$$ C^{\mu\sigma}(\mathbf{x}) = \frac{1}{2} \mathbf{x}^{\top} Q^{\mu\sigma} \mathbf{x} \quad \text{and} \quad C^{r}(\mathbf{x}) = \frac{1}{2} \mathbf{x}^{\top} Q^{r} \mathbf{x} \quad (11) $$

with $Q^{\mu\sigma} = \begin{pmatrix} H^{\mu} + H^{\sigma} & K \\ K^{\top} & L \end{pmatrix}$ and $Q^r = \begin{pmatrix} H^r & K \\ K^{\top} & L \end{pmatrix}$ two $2n \times 2n$ symmetric matrices, with $H^{\mu}, H^{\sigma}, H^r, K$ and L five $n \times n$ matrices defined by:
with $E_0$ a $n \times 2$ matrix defined by:

$$
\left\{
\begin{array}{ll}
H_{ij}^{\mu} = -S_{ij}\mu_{ij}\mu_{ji} & \text{if } i \neq j, \text{ and } H_{ii}^{\mu} = \sum_{j=1}^{n} S_{ij}\mu_{ij}^{2} \\
H_{ij}^{\sigma} = -S_{ij}\sigma_{ij}\sigma_{ji} & \text{if } i \neq j, \text{ and } H_{ii}^{\sigma} = \sum_{j=1}^{n} S_{ij}\sigma_{ij}^{2} \\
H_{ij}^{r} = -S_{ij}E[v_i v_j] & \text{if } i \neq j, \text{ and } H_{ii}^{r} = \sum_{j=1}^{n} S_{ij}(\mu_{ij}^2 + \sigma_{ij}^2) \\
K_{ij} = -S_{ij}\mu_{ij} & \text{if } i \neq j, \text{ and } K_{ii} = \sum_{j=1}^{n} S_{ij}\mu_{ij} \\
L_{ij} = -S_{ij} & \text{if } i \neq j, \text{ and } L_{ii} = \sum_{j=1}^{n} S_{ij}
\end{array}
\right.
\qquad (12)
$$

In the restricted case where the color correc-
tion models are linear transformations with zero y-
intercept [Eqn. (4b)], the vector of parameters is lim-
ited to **a** = (a₁, ..., aₙ) with aᵢ the slope coefficient
related to the transformation of image i. Cost func-
tions Cμσ and C'r are then equal to (Appendix VII):

$$
E_0 = \begin{pmatrix} S_{11}\mu_{11} & \dots & S_{nn}\mu_{nn} \\ S_{11}\sigma_{11} & \dots & S_{nn}\sigma_{nn} \end{pmatrix}
$$

3. Problem solving

The problem of finding an optimal set of color cor-
rection models that minimizes the total cost function
$C^{\mu\sigma}(\mathbf{x})$ or $C^r(\mathbf{x})$ can be formulated as a quadratic
programming (QP) problem, that is, a problem of op-
timizing a quadratic function of several variables sub-
ject to linear constraints. The optimal set of parame-
ters is given by:

$$
C^{\mu\sigma}(\mathbf{a}) = \frac{1}{2} \mathbf{a}^{\top} (H^{\mu} + H^{\sigma}) \mathbf{a} \quad \text{and} \quad C^{r}(\mathbf{a}) = \frac{1}{2} \mathbf{a}^{\top} H^{r} \mathbf{a} \qquad \mathbf{x}^{*} = \underset{\mathbf{x} \in \mathbb{R}^{2n}}{\operatorname{argmin}} \left\{ \frac{1}{2} \mathbf{x}^{\top} Q \mathbf{x} \right\} \quad \text{subject to} \quad E\mathbf{x} = \mathbf{d} \quad (18)
$$

2. Equality constraints

To ensure that we find consistent and non-
degenerated solutions to the optimisation problem,
we must add constraints on the parameters of color
correction models. First, the color correction models
should preserve the intensity of the mosaic, i.e. the
mean value of the n images:

$$
\frac{1}{n} \sum_{i=1}^{n} S_{ii} \mu_{ii} = \frac{1}{n} \sum_{i=1}^{n} S_{ii} \mu'_{ii} \quad (14)
$$

Besides, it is also desirable that the mean spread of
pixel values is preserved:

$$
\frac{1}{n} \sum_{i=1}^{n} S_{ii} \sigma_{ii} = \frac{1}{n} \sum_{i=1}^{n} S_{ii} \sigma'_{ii} \quad (15)
$$

These two constraints can be written in matrix no-
tations (Appendix VIII for detailed calculations):

$$
E\mathbf{x} = \mathbf{d} \tag{16}
$$

with $E$ a $2n \times 2$ matrix and $\mathbf{d}$ a $2$-vector defined by:

$$
E = \begin{pmatrix} S_{11}\mu_{11} & \dots & S_{nn}\mu_{nn} & S_{11} & \dots & S_{nn} \\ S_{11}\sigma_{11} & \dots & S_{nn}\sigma_{nn} & 0 & \dots & 0 \end{pmatrix}
$$
and
$$
\mathbf{d} = \begin{pmatrix} \sum_{i=1}^{n} S_{ii}\mu_{ii} \\ \sum_{i=1}^{n} S_{ii}\sigma_{ii} \end{pmatrix}
$$

In the restricted case of linear color correction mod-
els with zero y-intercept [Eqn. (4b)], only the first n
components of **x** are involved, and the equality con-
straint is given by:

$$
E_0 \mathbf{a} = \mathbf{d} \tag{17}
$$

with $Q = Q^{\mu\sigma}$ or $Q = Q^r$ depending on which cost function is to be minimized. In the restricted case where the color correction models are linear transformations with zero y-intercept, the vector of parameters is limited to $\mathbf{a} = (a_1, \dots, a_n)$, and the QP problem is formulated as:

$$
\mathbf{a}^* = \underset{\mathbf{a} \in \mathbb{R}^n}{\operatorname{argmin}} \left\{ \frac{1}{2} \mathbf{a}^\top H \mathbf{a} \right\} \quad \text{subject to} \quad E_0 \mathbf{a} = \mathbf{d} \quad (19)
$$

with $H = H^\mu + H^\sigma$ or $H = H^r$ depending on which cost function is to be minimized.

As the total cost $C^{\mu\sigma}(\mathbf{x})$ or $C^r(\mathbf{x})$ is strictly positive for any set of non-harmonized set of images, we have $\mathbf{x}^\top Q \mathbf{x} > 0$. This makes $Q$ a symmetric positive definite matrix. Therefore the quadratic function is strictly convex, and has a unique finite global minimum: this particular case of quadratic programming is called strictly convex quadratic program. To solve it, we use the primal-dual interior point algorithm implemented in the OOQP library [24].

III. DATA

Our method for color harmonization in satellite im-
age mosaicking was applied to a collection of high re-
solution Rapid Eye images, which was acquired over
mainland France (6.7 × 10⁵ km²) between 2009 and
2010 within the EQUIPEX GEOSUD project[23]. All
available images from the annual collection were used
without any pre-selection. In a pre-processing phase,
cloud-free images were obtained using a dedicated cut-
line network. Table I displays the main characteristics
of the image collection, while Fig. 3 shows the spa-
tial layout of the images. The number of overlaps per
image ranges from 2 to 13, with a mean value of 5
overlaps per image. The minimum overlap ratio for
an input image is 8.6% of its total area. As a whole,
overlapping parts cover 22% of the total mosaic area.
<table><caption>TABLE I: Characteristics of the image collection</caption><thead><tr><th>Characteristic</th><th>Value</th></tr></thead><tbody><tr><td>Satellite</td><td>RapidEye</td></tr><tr><td>Acquisition dates</td><td>May to October 2009<br>+ May to October 2010</td></tr><tr><td>Azimuth angle</td><td>100°</td></tr><tr><td>Image swath</td><td>Approximately 77 km</td></tr><tr><td>Spatial reference system</td><td>EPSG 2154</td></tr><tr><td>Spectral bands used</td><td>Red, Green, Blue</td></tr><tr><td>Number of images</td><td>132</td></tr><tr><td>Total spatial extent</td><td>2.3 × 10<sup>6</sup> km<sup>2</sup></td></tr><tr><td>Pixel size (orthorectified)</td><td>5 m</td></tr><tr><td>Pixel values encoding</td><td>16 bits</td></tr><tr><td>Total input images pixel count</td><td>91 × 10<sup>9</sup> pix.</td></tr><tr><td>Mean pixel count per image</td><td>6.9 × 10<sup>8</sup> pix.</td></tr><tr><td>Size of images files</td><td>550 Gb</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>FIG. 3: Layout of the 132 cloud-free RapidEye ortho-images, with selected pair of images A and B (bold) and overlapped area (in red)</figcaption></figure>

## IV. RESULTS

The color harmonization method presented in section II was implemented in C++ using OTB [26], ITK [27] and OOQP [24] libraries. For sake of clarity, we first analyse detailed results obtained on a single pair of partly overlapping images A and B (section IVA), then we display results related to the complete mosaic over mainland France (section IVB), and finally discuss some computational issues (section IV C).

### A. Simple case analysis: two images

We insist on the fact that our approach is designed to deal with any number of overlapping images in the same area. However, for sake of clarity, a pair of neighbouring images denoted A and B was selected to give detailed results on a simplified test case

(Fig. 3). We first applied the color harmonization method with general linear transformation and mean-and-deviation-based cost function $C^{\mu\sigma}$ to the pair of images A and B, which are overlapping in a small area (Fig. 4). Without color harmonization, image A displays a near-orange tone, whereas image B has a more near-green tone (Fig 4b). After correction (Fig. 4c), images have the same tone, and the boundary zone is almost indistinguishable. Fig. 4 also clearly suggests that the final tone of the mosaic is more influenced by the original tone of image A than by the original tone of image B. This results from the equality constraint on global intensity preservation [Eqn. (14)]: because image A is much larger than image B ($S_{AA} \gg S_{BB}$), corrected pixels from image A are closer to their original values than corrected pixels from image B.

Fig. 5 shows how pixel values of images A and B were transformed in the RGB color space. It appears that the transformation is not linear: the explanation is that linear color correction models are applied in the $l\alpha\beta$ color space, rather than directly in the original RGB channels. Again, we notice that pixel values in image B undergo a tougher transformation than in image A, resulting in a correction function with stronger curvature. The same figure plotted in the $l\alpha\beta$ color space, rather than in the RGB color space, would show a perfectly linear relationship between original pixel values $v$ and transformed pixel values $v'$, resulting from linearity in Eqn. (4).

Finally, Fig. 6 displays pixels distribution in $l\alpha\beta$ color space for both images A and B in their overlapping area, before and after correction. It shows how mean-variance equalization leads to corrected histograms for A and B that are close to one another, yet not exactly matching. In particular, the distribution modes for images A and B—unimodal distribution for $l$ channel, bimodal distribution for $\alpha$ and $\beta$ channels—are closer after correction than before. Table II presents an analysis of the range of images A and B pixels values in $l\alpha\beta$ color space. It shows how correction brings closer means and standard deviations of images: without correction, the absolute difference between mean values of A and B in $l\alpha\beta$ color space is reduced by [-92%, -97%, -100%] ($l, \alpha, \beta$ channels respectively). Same convergence appears with standard deviations: the absolute difference between standard deviations is reduced by [-73%, -55%, -40%] after correction.

### B. The France mosaic

Our approach for color harmonization was then applied to the whole collection of 132 RapidEye images over mainland France described in Section III. The output mosaic (Fig. 7) is 221320 × 244320 pixels in 8 bits encoding, 5m resolution, in GeoTIFF format, for a total size of 162.2 Gb. In comparison, we also tried to generate the entire output mosaic using a popular commercial software but it failed to process the entire dataset, most probably because of its excessive size.
<figure><img src="image_6.png" /><figcaption>FIG. 4: Color harmonization for pair of overlapping images A and B.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 5: Corrected pixel values vs. original pixel values for both images A and B in RGB color space.</figcaption></figure>

<table><caption>TABLE II: Analysis of the range of values in $l\alpha\beta$ space after color transformation</caption><thead><tr><th></th><th>Image A<br>(original)</th><th>Image B<br>(original)</th><th>Image A<br>(corrected)</th><th>Image B<br>(corrected)</th><th>Abs. Difference<br>(original)</th><th>Abs. Difference<br>(corrected)</th><th>Relative reduction %</th></tr></thead><tbody><tr><td>mean (l)</td><td>2.19E+00</td><td>2.31E+00</td><td>2.21E+00</td><td>2.22E+00</td><td>1.18E-01</td><td>9.30E-03</td><td>-92%</td></tr><tr><td>mean (α)</td><td>-6.01E-02</td><td>-5.85E-02</td><td>-6.03E-02</td><td>-6.02E-02</td><td>1.53E-03</td><td>4.10E-05</td><td>-97%</td></tr><tr><td>mean (β)</td><td>-1.59E-02</td><td>-2.24E-02</td><td>-1.75E-02</td><td>-1.75E-02</td><td>6.48E-03</td><td>1.50E-05</td><td>-100%</td></tr><tr><td>std (l)</td><td>2.17E-01</td><td>1.72E-01</td><td>2.10E-01</td><td>1.98E-01</td><td>4.53E-02</td><td>1.24E-02</td><td>-73%</td></tr><tr><td>std (α)</td><td>4.70E-02</td><td>3.74E-02</td><td>4.43E-02</td><td>4.85E-02</td><td>9.58E-03</td><td>4.28E-03</td><td>-55%</td></tr><tr><td>std (β)</td><td>1.25E-02</td><td>1.07E-02</td><td>1.19E-02</td><td>1.30E-02</td><td>1.77E-03</td><td>1.07E-03</td><td>-40%</td></tr></tbody></table>

Five different output mosaics were produced (Fig. 7), using either a linear transformation with zero y-intercept [Eqn. (4b)] or a general linear transformation [Eqn. (4a)], and minimizing either the mean-based cost $C^{\mu}$, the mean-and-deviation-based cost $C^{\mu\sigma}$, or the RMSE-based cost $C^r$. They are compared to a mosaic generated without any color correction. Mosaic characteristics are summarized in Tab. III. All mosaics offer satisfying visual results at the scale of France. However, a zoom over a 200 × 165 km square proves that there are significant differences between them. As explained in the introduction, our main

goal is to produce visually appealing mosaics in a fully automated way: the assessment of the five output mosaics is thus mainly qualitative. However, we also tried to compare the quality of the mosaics using a quantitative metric: we computed the max Peak Signal-to-Noise Ratio (PSNR) over all overlapping parts of images. The **PSNR** is a popular performance evaluation metric used to compare two similar images. It is computed as $PSNR = 10 \cdot \log_{10}(peak^2/MSE)$ with peak = $2^8 \cdot \sqrt{3}$ for 8 bits RGB images, and **MSE** the mean squared error between pairs of overlapping images over the whole mosaic. Table III shows
<figure><img src="image_8.png" /><figcaption>FIG. 6: Pixel distribution for *l*, *α*, *β* channels for images *A* (black lines) and *B* (red lines) in their overlapping area, before (solid lines) and after (dashed lines) correction.</figcaption></figure>

<table><caption>TABLE III: Characteristics of France mosaics with quality metrics</caption><thead><tr><th>Mosaic id</th><th>Correction model</th><th>Cost function</th><th>PSNR</th></tr></thead><tbody><tr><td>#0</td><td>none</td><td>none</td><td>32.948 dB</td></tr><tr><td>#1</td><td>v' = a · v</td><td>C<sup>μ</sup></td><td>34.660 dB</td></tr><tr><td>#2</td><td>v' = a · v</td><td>C<sup>μσ</sup></td><td>34.722 dB</td></tr><tr><td>#3</td><td>v' = a · v + b</td><td>C<sup>μσ</sup></td><td>34.960 dB</td></tr><tr><td>#4</td><td>v' = a · v</td><td>C<sup>r</sup></td><td>34.969 dB</td></tr><tr><td>#5</td><td>v' = a · v + b</td><td>C<sup>r</sup></td><td>35.413 dB</td></tr></tbody></table>

PSNR for the five generated mosaics. Mosaic #5 gives the highest PSNR metric: it corresponds to a mosaic generated with a general linear transformation [Eqn. (4a)], minimizing the RMSE-based cost function $C^r$. However, PSNR metric is only an approximation to human perception of mosaic quality: although a higher PSNR generally indicates that the mosaic is of higher quality, one may prefer the visual appearance of mosaic #3, which was corrected using mean-and-deviation-based cost function $C^{μσ}$.

### C. Computation time

Computational issues are a key limiting factor when generating large size image mosaics. For this reason, we give here some detailed elements on computation

times. The complete process of generating a mosaic following our method for color harmonization can be decomposed into three stages: i) computing image statistics in overlapping areas ($μ_{ij}$, $μ_{ji}$, $σ_{ij}$, $σ_{ji}$ and $E[v_i v_j]$), ii) solving the QP optimization problem described in Section II, and iii) composing the final output mosaic.

Statistics computing and mosaic composition are CPU time consuming, because both need to access pixel data of input images. In addition, mosaic composing involves saving data to disk. If *m* is the total number of pixels to process, i.e. the total number of pixels in all *n* input images, the computational complexity of these two steps is $O(m)$. The so-called streaming mechanism of ITK library is used to sequentially process sub-regions of the images, which allowed us to process large datasets without restriction on images sizes. Statistics computing and mosaic composing were fully multi-threaded in order to take advantage of the latest computing architectures.

To better analyse the computational performance of our approach, we decomposed run time into the three stages (statistics computing, optimization, mosaic composing) for several sets of adjacent input images randomly selected from the complete dataset of 132 RapidEye images, for increasing number *n* of images (Fig. 8). All tests ran on a personal computer Intel(R) Xeon(R) with CPU E5-1650 @ 3.20GHz. Results clearly show that total run time increases linearly with the total number of pixels to process. For the entire RapidEye dataset over mainland France, statistics
<figure><img src="image_9.png" /><figcaption>FIG. 7: Output mosaics over mainland France and a zoom on a 200 × 165 km square: mosaic #0 (a-b), mosaic #1 (c-d), mosaic #2 (e-f), mosaic #3 (g-h), mosaic #4 (i-j), mosaic #5 (k-l). Tab III summarizes mosaic characteristics.</figcaption></figure>

were computed in 9h and final mosaic was composed
in 35h.

V. DISCUSSION

A. An innovative color harmonization method
for large sets of remote sensing images

Our first goal was to provide a robust method for
color harmonization to produce large mosaics from re-
mote sensing images. The formulation we give of the
color harmonization problem as a Quadratic Program-
ming (QP) problem overcomes some of the limitations
of existing techniques. It allows to deal simultane-
ously, rather than sequentially, with any number *n*
of input images and any spatial layout of these im-
ages and their overlapping areas. The optimal pa-
rameters *a* and *b* of color correction models for all im-
ages are searched simultaneously in a 2*n*-dimensional
space. As discussed in the introduction, this over-
comes the limitation of color harmonization tech-

Finally, Fig. 9 displays run time for the QP solving
stage, for increasing number *n* of input images. This
optimization step is very fast thanks to an efficient im-
plementation in OOQP library that takes advantage
of sparse matrices. The run time is higher for a gen-
eral linear color correction model [Eqn. (4a)], due to
the fact that the corresponding QP has twice more un-
known variables (*a* and *b* are unknown) than the zero
*y*-intercept case where only *a* is unknown [Eqn. (4b)].
Total QP solving time for the complete dataset of 132
RapidEye images is only 3 ms.
<figure><img src="image_10.png" /><figcaption>FIG. 8: Total run time for increasing number of input pixels, decomposed into statistics computation (white bars) and output mosaic composition (black bars). Run time devoted to the solving of QP optimization does not appear here because it is much smaller.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9: Run time to solve QP optimization problem for increasing number of input images, for a general (boxes) or zero y-intercept (crosses) linear color correction model.</figcaption></figure>

niques that work on pair of images only, in which several different paths can be taken to build sequentially a large mosaic from a set of *n* partly overlapping images, leading to non-unique results. Note that the use of the QP framework is closely linked to the linearity of the color correction models we used in the *l*αβ color space [Eqn. (4)]. More sophisticated and non-linear color correction models would surely give better results on single pairs of overlapping images, but would not be suitable for computationally efficient QP solving.

The second innovative part of our method is the use of the *l*αβ decorrelated color space, which was defined by [19] based on the characteristics of the human eye. Our work was an attempt to adapt this concept from the computer vision community to the context of remote sensing images. Our results clearly show that the use of the *l*αβ decorrelated color space is well suited to RGB satellite imagery, even if the quality of the resulting mosaic can only be assessed qualitatively by visual inspection.

## B. A fully-automated, computationally efficient tool

Our research also sought to result in an operational, fully-automated tool for color harmonization of large image sets. We efficiently implemented our approach using open-source libraries. No tuning parameter has to be set by the user. Total running time increases linearly with the number of input pixels. The main time-consuming task is the computation of image statistics in their overlapping parts: this could be reduced by computing statistics on a small sample of pixels, either randomly selected, or selected from invariant features. The amount of overlapping areas between input images is a key factor that drives the output quality of the mosaic. In our case study over mainland France, 22% of total mosaic area was covered by overlapping images; further research is needed to find extremal or even optimal values for this overlap ratio.

The tool was successfully applied to generate a mosaic over mainland France from RapidEye imagery at 5m resolution. In the near future, it will be used to produce, on a yearly basis, a high resolution mosaic of mainland France from satellite imagery acquired through the GEOSUD platform and the THEIA Land Data Center[25].

## C. Limits and further research

It should be noted that the presented method also has limits, which may restrict its use in some cases. To start with, the definition of the cost function *C*, that the image mosaic must minimize, could be questioned. In this study, two cost functions were used: one based on the mean and variance of the *n* input images in their overlapping parts [Eqn. (7) and (8)], the other one based on the RMSE criterion. A possible way forward would be to move from this mean-variance or RMSE paradigm to a broader, non-parametric framework, that could for example build on the histogram matching techniques to define a cost associated to a pair (i, j) of images based on a complete comparison of their histogram in their overlapping part.

Another possible improvement would be to use more sophisticated color correction models, different from the one we used from [19]. Many other color transfer techniques could be investigated — see [18] for a review of such techniques. However, only linear correction models allow to formulate the multi-images color harmonization as a Quadratic Programming problem. Moving to non-linear color correction models would thus prevent the use of efficient QP solvers to find optimal model parameters.

Besides, a key drawback of our approach is that it is based on the assumption that the composition of two neighboring images *i* and *j* in their overlapping part is similar—e.g. same proportion of different land covers in both images. Color differences are assumed to be mainly due to atmospheric or seasonal conditions. Hence, the color harmonization process may be biased
by color differences from i to j which are due to: i) land use changes between dates of acquisition of images i and j (e.g., clearcut harvesting), or ii) presence of clouds in i or j. To avoid this problem, cloud masking algorithms should be applied to input images before color harmonization. Our approach could also be improved by first pre-processing all input images using segmentation and classification algorithms to group pixels of overlapping areas depending on their vegetation type of surface material. Different linear models could then be fitted to each group of pixels. However, such a pre-processing step would require large CPU time, and may not be fully automated.

Finally, our approach is only suitable for optical imagery of Earth in natural colors (RGB). Further work would be needed to extend it to hyperspectral imagery (e.g. false color mosaic).

## VI. CONCLUSION

This work was carried out with a view towards producing high-resolution mosaics from RGB satellite images for illustrative use only. We formulated the color harmonization process as a Quadratic Programming (QP) problem, which allows us to deal simultaneously with any number of input images, with any spatial layout of their overlapping parts, and to use efficient numerical solvers to find the optimal linear color correction models. Besides, we introduced the use of the $l\alpha\beta$ decorrelated color space, which is ideally suited to process Earth observation scenes. The method was successfully implemented using open source libraries and will be applied every year to produce a high-resolution mosaic of mainland France from images acquired through the GEOSUD platform and the THEIA Land Data Center. Further research could focus on improving the definition of the cost function in the QP problem, or on enhancing the statistics computation by searching for invariant features in overlapping images.

# APPENDIX

## VII. CANONIC FORMULATION OF COST FUNCTIONS $C^{\mu\sigma}$ AND $C^r$

### A. Canonic formulation of cost function $C^{\mu\sigma}$

Mean-and-deviation-based cost function $C^{\mu\sigma}$ is equal to the sum of mean-based cost function $C^\mu$ and deviation-based cost function $C^\sigma$.

1. The mean-based cost function $C^\mu$ writes

$$C^{\mu} = \sum_{i=1}^{n} \sum_{j=i}^{n} c_{ij} = \frac{1}{4} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij}(\mu'_{ij} - \mu'_{ji})^2$$

The mean value after image color correction is linked to the mean value before correction by $\mu'_{ij} = a_i\mu_{ij} + b_i$ with $a_i$ and $b_i$ the slope and intercept coefficients related to the image i. We thus obtain the expression of cost function $C^\mu$ as a function of vector $\mathbf{x} = (a_1, \dots, a_n, b_1, \dots, b_n) = (\mathbf{a}, \mathbf{b})$:

$$C^{\mu}(\mathbf{x}) = \frac{1}{4} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} \left[ (a_i \mu_{ij} + b_i) - (a_j \mu_{ji} + b_j) \right]^2$$

We develop the expression:

$$C^{\mu}(\mathbf{x}) = \frac{1}{4} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} \cdot \left[ (a_i \mu_{ij} + b_i)^2 - 2(a_i \mu_{ij} + b_i)(a_j \mu_{ji} + b_j) + (a_j \mu_{ji} + b_j)^2 \right]$$

By swapping i and j in the sums, we notice that $\sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} \cdot [(a_i \mu_{ij} + b_i)^2] = \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} \cdot [(a_j \mu_{ji} + b_j)^2]$, and we obtain:

$$C^{\mu}(\mathbf{x}) = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} \cdot \left[ (a_i \mu_{ij} + b_i)^2 - (a_i \mu_{ij} + b_i)(a_j \mu_{ji} + b_j) \right] = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} \cdot \left[ \mu_{ij}^2 a_i^2 + 2\mu_{ij} a_i b_i + b_i^2 - \mu_{ij} \mu_{ji} a_i a_j - \mu_{ij} a_i b_j - \mu_{ji} a_j b_i - b_i b_j \right]$$

Using the definition of the three $n \times n$ matrices $H^\mu$, $K$ and $L$ given in Eqn. (12), we obtain:

$$C^{\mu}(\mathbf{x}) = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} (H_{ij}^{\mu} a_i a_j + K_{ij} a_i b_j + K_{ji} a_j b_i + L_{ij} b_i b_j)$$

Using matrix notation, we get:

$$C^{\mu}(\mathbf{a}, \mathbf{b}) = \frac{1}{2} \begin{pmatrix} \mathbf{a} \\ \mathbf{b} \end{pmatrix}^{\top} \begin{pmatrix} H^{\mu} & K \\ K^{\top} & L \end{pmatrix} \begin{pmatrix} \mathbf{a} \\ \mathbf{b} \end{pmatrix} \quad (20)$$

When the color correction model is a linear transformation with zero y-intercept, the vector of parameters is limited to $\mathbf{a} = (a_1, \dots, a_n)$. In this restricted case, the cost function has the same expression as $C^\mu$ with all $b_i = 0$, and we directly obtain its final expression by putting $\mathbf{b} = \mathbf{0}$ in Eqn. (20):

$$C^{\mu}(\mathbf{a}) = \frac{1}{2} \mathbf{a}^{\top} H^{\mu} \mathbf{a}$$

2. Deviation-based cost function $C^\sigma$ writes

$$C^{\sigma} = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=i}^{n} S_{ij} (\sigma'_{ij} - \sigma'_{ji})^2$$

No matter which type of linear color correction model is used (Eqn. (4a) or Eqn. (4b)), the standard
deviation $\sigma'_{ij}$ of image $i$ in the overlap area with image $j$ after correction is linked to that before correction by $\sigma'_{ij} = |a_i|\sigma_{ij}$ — it does not depend on $b_i$ — with $a_i$ the slope coefficient related to image $i$ and with all $a_i > 0$. Similarly to $C^\mu$ we have:

$$C^\sigma(\mathbf{x}) = C^\sigma(\mathbf{a}) = \frac{1}{4} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} (a_i \sigma_{ij} - a_j \sigma_{ji})^2 \\ = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} (\sigma_{ij}^2 a_i^2 - \sigma_{ij} \sigma_{ji} a_i a_j)$$

Moving from the second to the third line above is obtained by noticing that $\sum_{i,j=1}^n S_{ij}\sigma_{ij}^2 a_i^2 = \sum_{i,j=1}^n S_{ij}\sigma_{ji}^2 a_j^2$ by swapping $i$ and $j$ indices in the double sum. We obtain:

$$C^\sigma(\mathbf{x}) = C^\sigma(\mathbf{a}) = \frac{1}{2} \mathbf{a}^\top H^\sigma \mathbf{a}$$

with $H^\sigma$ a $n \times n$ matrix defined in Eqn. (12).

### 3. Mean-and-deviation-based cost function $C^{\mu\sigma}$

by summing the expression of mean-based cost function $C^\mu$ and deviation-based cost function $C^\sigma$, we find the expression of mean-and-deviation-based cost function $C^{\mu\sigma}$:

$$C^{\mu\sigma}(\mathbf{x}) = C^{\mu}(\mathbf{x}) + C^{\sigma}(\mathbf{x}) = \frac{1}{2} \mathbf{x}^{\top} Q^{\mu\sigma} \mathbf{x}$$

with $Q^{\mu\sigma} = \begin{pmatrix} H^\mu + H^\sigma & K \\ K^\top & L \end{pmatrix}$ a $2n \times 2n$ symmetric matrix.

## B. Canonic formulation of cost function $C^r$

The RMSE-based cost function $C^r$ writes:

$$C^r = \frac{1}{4} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} (r'_{ij})^2$$

with $r'_{ij}$ the RMSE after transformation defined in Eqn. (5):

$$\begin{aligned} (r'_{ij})^2 &= E[(v'_{i} - v'_{j})^2] \\ &= E[( (a_i v_i + b_i) - (a_j v_j + b_j))^2 ] \\ &= a_i^2 (\sigma_{ij}^2 + \mu_{ij}^2) + 2a_i b_i \mu_{ij} + b_i^2 \\ &+ a_j^2 (\sigma_{ji}^2 + \mu_{ji}^2) + 2a_j b_j \mu_{ji} + b_j^2 \\ &- 2 \cdot (a_i a_j E[v_i v_j] + a_j b_i \mu_{ji} + a_i b_j \mu_{ji} + b_i b_j) \end{aligned}$$

We develop the expression of $C^r$:

$$C^r(\mathbf{x}) = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} S_{ij} \cdot \left[ \begin{aligned} &a_i^2 (\sigma_{ij}^2 + \mu_{ij}^2) + 2a_i b_i \mu_{ij} - 2a_i b_j \mu_{ij} \\ &- a_i a_j E[v_i v_j] + (b_i - b_j)^2 \end{aligned} \right]$$

Using the definition of the three $n \times n$ matrices $H^r$, $K$ and $L$ given in Eqn. (12), we obtain:

$$C^r(\mathbf{x}) = \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} (H_{ij}^r a_i a_j + K_{ij} a_i b_j + K_{ji} a_j b_i + L_{ij} b_i b_j)$$

Using matrix notation, we get:

$$C^r(\mathbf{a}, \mathbf{b}) = \frac{1}{2} \begin{pmatrix} \mathbf{a} \\ \mathbf{b} \end{pmatrix}^\top \begin{pmatrix} H^r & K \\ K^\top & L \end{pmatrix} \begin{pmatrix} \mathbf{a} \\ \mathbf{b} \end{pmatrix}$$

which writes:

$$C^r(\mathbf{x}) = \frac{1}{2} \mathbf{x}^\top Q^r \mathbf{x}$$

with $Q^r = \begin{pmatrix} H^r & K \\ K^\top & L \end{pmatrix}$ a $2n \times 2n$ symmetric matrix.

## VIII. EQUALITY CONSTRAINTS

The first equality constraint is that the global mean intensity of all the images should be preserved by the color correction models:

$$\frac{1}{n} \sum_{i=1}^{n} S_{ii} \mu_{ii} = \frac{1}{n} \sum_{i=1}^{n} S_{ii} \mu'_{ii}$$

Using $\mu'_{ij} = a_i \mu_{ij} + b_i$ with $a_i$ and $b_i$ the slope and intercept coefficients related to the image $i$, we obtain:

$$\sum_{i=1}^{n} S_{ii} \mu_{ii} = \sum_{i=1}^{n} S_{ii} (\mu_{ii} a_i + b_i)$$

which writes:

$$(S_{11}\mu_{11}, \dots, S_{nn}\mu_{nn}, S_{11}, \dots, S_{nn}) \cdot \begin{pmatrix} a_1 \\ \vdots \\ a_n \\ b_1 \\ \vdots \\ b_n \end{pmatrix} = \sum_{i=1}^{n} S_{ii} \mu_{ii}$$

The second equality constraint is that the mean spread of pixel values over all images should be preserved by the color correction models:

$$\frac{1}{n} \sum_{i=1}^{n} S_{ii} \sigma_{ii} = \frac{1}{n} \sum_{i=1}^{n} S_{ii} \sigma'_{ii}$$

The standard deviation value after image color correction is linked to the standard deviation value before correction by $\sigma'_{ij} = a_i \sigma_{ij}$ and we obtain:

$$\sum_{i=1}^{n} S_{ii} \sigma_{ii} = \sum_{i=1}^{n} S_{ii} \sigma_{ii} a_i$$
which writes:

$$ (S_{11}\sigma_{11}, \dots, S_{nn}\sigma_{nn}) \cdot \begin{pmatrix} a_1 \\ \vdots \\ a_n \end{pmatrix} = \sum_{i=1}^{n} S_{ii}\sigma_{ii} $$

Therefore the two equality constraints (intensity preservation and spread preservation) can be written at once in matrix notations as shown in Eqn. (16).

the *Investissements d'Avenir* program managed by the French National Research Agency. The authors would also like to thank the anonymous reviewers and Agnès Bégué for their valuable comments and suggestions which have greatly contributed in improving the manuscript. They also thank the French National Institute of Geographic and Forestry Information (IGN) for providing the seamline network.

## ACKNOWLEDGMENTS

This work was supported by public funds received through GEOSUD, a project (ANR-10-EQPX-20) of

[1] J. Cihlar. Land cover mapping of large areas from satellites: status and research priorities. *International Journal of Remote Sensing*, 21:1093–1114, 2000.

[2] R. Bindschadler, P. Vornberger, A. Fleming, A. Fow, J. Mullins, D. Binnie, S. J. Paulsen, B. Granneman, and D. Gorodetzky. The Landsat image mosaic of Antarctica. *Remote Sensing of Environment*, 112:4214–4226, 2008.

[3] Y. Afek and A. Brand. Mosaicking of orthorectified aerial images. *Photogrammetric Engineering & Remote Sensing*, 64:115–125, 1998.

[4] Y. Shiren, L. Li, and G. Peng. Two-dimensional seampoint searching in digital image mosaicking. *Photogrammetric Engineering & Remote Sensing*, 55:49–53, 1989.

[5] J. Pan, M. Wang, D. Li, and J. Li. Automatic generation of seamline network using area Voronoi diagrams with overlap. *IEEE Transactions on Geoscience and Remote Sensing*, 47:1737–1744, 2009.

[6] M. Kerschner. Seamline detection in colour orthomage mosaicking by use of twin snakes. *Journal of Photogrammetry and Remote Sensing*, 56:53–64, 2001.

[7] P. J. Burt and E. H. Adelson. A multiresolution spline with application to image mosaics. *ACM Transactions on Graphics*, 2:217–236, 1983.

[8] R. Szeliski. Image alignment and stitching: A tutorial. Technical report, Microsoft Research, 2004.

[9] M.K. Ng and W. Wang. A variational approach for image stitching II: using image gradients. *SIAM J. Imaging Sci.*, 6:1345–1366, 2013.

[10] B. Zitova and J. Flusser. Image registration methods: a survey. *Image and Vision Computing*, 21:977–1000, 2003.

[11] W. Xu and J. Mulligan. Performance evaluation of color correction approaches for automatic multi-view image and video stitching. In *IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 263–270, 2010.

[12] Victor Jung-Der Tsa and Ying-Ting Huang. Automated image mosaicking. *Journal of the Chinese Institute of Engineers*, 28:329–340, 2005.

[13] D. Yuan and C. D. Elvidge. Comparison of relative radiometric normalization techniques. *Journal of Photogrammetry & Remote Sensing*, 51:117–126, 1996.

[14] F.G. Hall, D.E. Strebel, J.E. Nickeson, and S.J. Goetz. Radiometric rectification: Toward a common radiometric response among multidate, multisensor images. *Remote Sensing of Environment*, 35:11–27, 1991.

[15] C. D. Elvidge, D. Yuan, R. D. Weerackon, and R. S. Lunetta. Relative radiometric normalization of Landsat multispectral scanner data using an automated scattergram-controled regression. *Photogrammetric Engineering & Remote Sensing*, 61:1255–1260, 1995.

[16] J. R. Schott, C. Salvaggio, and W. J. Volchok. Radiometric scene normalization using pseudoinvariant features. *Remote Sensing of Environment*, 26:1–16, 1988.

[17] B. Guindon. Assessing the radiometric fidelity of high resolution satellite image mosaics. *Journal of Photogrammetry and Remote Sensing*, 52:229–243, 1997.

[18] H. S. Faridul, T. Pouli, C. Chamaret, J. Stauder, A. Tremeau, and E. Reinhard. A survey of color mapping and its applications. *SIAM J. Imaging Sci.*, 6:1345–1366, 2014.

[19] E. Reinhard, M. Ashikhmin, B. Gooch, and Shirley P. Color transfer between images. *IEEE Computer Graphics and Applications*, 21:34–41, 2001.

[20] D. L. Ruderman, T. W. Cronin, and C. C. Chiao. Statistics of cone responses to natural images: Implications for visual coding. *Journal of the Optical Society of America*, 15:2036–2045, 1998.

[21] Y. Qian, D. Liao, and J. Zhou. Manifold alignment based color transfer for multiview image stitching. In *Int. Conf. Image Processing (ICIP)*, 2013.

[22] Y.-H. Shin, M.-G. Park, Y.-S. Jeon, Y.-S. Moon, S.-H. Lee, and K.-J Yoon. Tone correction with dynamic objects for seamless image mosaic. In *ECCV Workshops*, 2012.

[23] Matrix values in Eqn. (1) slightly differ from those given by [19, Eqn (2)], which were found to be erroneous by an anonymous reviewer. Correct values were directly computed from [19, Eqn (1)].

[24] E. M. Gertz and S. J. Wright. Object-oriented software for quadratic programming. *ACM Transactions on Mathematical Software*, 29:58–81, 2003.

[25] www.equipex-geosud.fr.

[26] J. Inglada and E. Christophe. The Orfeo Toolbox remote sensing image processing software. In *IEEE*
*International Geoscience and Remote Sensing Symposium (IGARSS)*, 2009.

[27] T. S. Yoo, M. J. Ackerman, W. E. Laorensen, W. Schroeder, V. Chalana, S. Aylward, D. Metaxas, and R. Whitaker. Engineering and algorithm design for an image processing API: A technical report on

ITK - The Insight Toolkit. In *IOS Press Amsterdam*
J. Westwood, ed., editor, *Proceedings of Medecine*
*Meets Virtual Reality*, 2002.

[28] http://www.theia-land.fr/en.