R Cresson
email: remi.cresson@teledetection.fr

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 overlaping 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 [START_REF] Cihlar | Land cover mapping of large areas from satellites: status and research priorities[END_REF]. 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 [START_REF] Bindschadler | The Landsat image mosaic of Antartica[END_REF]. 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 [START_REF] Afek | Mosaicking of orthorectified aerial images[END_REF][START_REF] Shiren | Two-dimensional seampoint searching in digital image mosaicking[END_REF][START_REF] Pan | Automatic generation of seamline network using area Voronoi diagrams with overlap[END_REF][START_REF] Kerschner | Seamline detection in colour orthoimage mosaicking by use of twin snakes[END_REF], iii) color harmonization of images, and iv) composition of the final mosaic, which may involve techniques such as image blending or image stitching [START_REF] Burt | A multiresolution spline with application to image mosaics[END_REF][START_REF] Szeliski | Image alignment and stitching: A tutorial[END_REF][START_REF] Ng | A variational approach for image stitching II: using image gradients[END_REF], 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 [START_REF] Zitova | Image registration methods: a survey[END_REF] 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 heterogeneousness in the mosaic. Adjustement 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 [START_REF] Xu | Performance evaluation of color correction approaches for automatic multi-view image and video stitching[END_REF]. 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 [START_REF] Victor | Automated image mosaicking[END_REF] goes one step further: the corrected image will have the same complete pixel value distribution as the reference image. Relative radiometric normalization techniques [START_REF] Yuan | Comparison of relative radiometric normalization techniques[END_REF][START_REF] Hall | Radiometric rectification: Toward a common radiometric response among multidate, multisensor images[END_REF][START_REF] Elvidge | Relative radiometric normalization of Landsat multispectral scanner data using an automated scattergram-controled regression[END_REF] and radiometric normalization tech-niques based on invariant features [START_REF] Schott | Radiometric scene normalization using pseudoinvariant features[END_REF] 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 image. To generate a mosaic with n > 2 images, one usually 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 assembled, and so on. Colorimetric correction is thus propagated from the first reference image to the whole mosaic. However, depending on the spatial layout of the n images and their overlapping areas, several paths may be taken to build the whole mosaic. [START_REF] Guindon | Assessing the radiometric fidelity of high resolution satellite image mosaics[END_REF] 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 collection of satellite images is large, with large amount of image overlaps.

Besides, the available color harmonization techniques 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 computer vision community that other color spaces are more suitable for color processing, but these developments have not reached the remote sensing community yet. Regarding this subject, [START_REF] Faridul | A survey of color mapping and its applications[END_REF] recently provided a extensive review of numerous color transfer or colour mapping techniques, which aim to recolor a given image 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 image 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αβ) 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 successfuly 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 advantages 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., radiometric fidelity is not sought after. We propose to build on the color transfer technique developed by [START_REF] Reinhard | Color transfer between images[END_REF] to define linear color correction models for each image, in a decorrelated color space. Besides, while [START_REF] Reinhard | Color transfer between images[END_REF] performed mean-variance equalization on pairs of images, we formulate the color correction problem for mosaic generation simultaneously from n images as a multivariate optimization problem, namely a quadratic programming (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αβ;

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 [START_REF] Burt | A multiresolution spline with application to image mosaics[END_REF] or a variational approach [START_REF] Ng | A variational approach for image stitching II: using image gradients[END_REF];

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 conversion 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. color space was first introduced by [START_REF] Ruderman | Statistics of cone responses to natural images: Implications for visual coding[END_REF], based on datadriven human perception research that assumes the human visual system is ideally suited for processing natural scenes. In the lαβ color space, the l axis represents an achromatic channel, while the α and β channels are chromatic yellow-blue and red-green opponent channels. The main feature of the lαβ color space is that, for natural scenes, the correlations between the three channels' values are low. [START_REF] Reinhard | Color transfer between images[END_REF] used this color space to transfer colors from a reference image to a target image, performing a simple mean-variance equalization in the lαβ space -this approach was further developed by [START_REF] Qian | Manifold alignment based color transfer for multiview image stitching[END_REF] and [START_REF] Shin | Tone correction with dynamic objects for seamless image mosaic[END_REF]. 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 developped by [START_REF] Reinhard | Color transfer between images[END_REF] seems to be suited to our needs. We thus perform RGB-tolαβ conversion on each image, resulting in the computation of three new channels l, α, and β, in which inter-channel correlations are expected to be low. The conversion of RGB values to lαβ values follows three steps-see [START_REF] Reinhard | Color transfer between images[END_REF] for more details:

• first,
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ L M S ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣
0.3811 0.5783 0.0406 0.1967 0.7244 0.0790 0.0241 0.1228 0.8531

⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⋅ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ R G B ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (1)
• values in the LMS space are most often highly skewed; we partly eliminate this skew by converting values to logarithmic scale [Eqn. [START_REF] Bindschadler | The Landsat image mosaic of Antartica[END_REF]];

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ L log M log S log ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ log L log M log S ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (2)
• finally, a linear transformation introduced by [START_REF] Ruderman | Statistics of cone responses to natural images: Implications for visual coding[END_REF] is applied to convert the signal from LMS values to the decorrelated lαβ color space [Eqn. ( 3)].

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ l α β ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 √ 3 0 0 0 1 √ 6 0 0 0 1 √ 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⋅ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 1 1 1 1 -2 1 -1 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⋅ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ L log M log S log ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (3)
Equations for the inverse transformation from lαβ space to RGB space can be derived from Eqn. ( 1) to [START_REF] Afek | Mosaicking of orthorectified aerial images[END_REF]. They are needed in the last step of satellite image mosaicking, to move back from the lαβ 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αβ color space is low, we can use independent correction models for each of the three axes l, α and β. 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 ⋅ v + b (4a) v ′ = a ⋅ v (4b)
where a and b are real-valued parameters, v (input value) is the l, α or β-value of a pixel in the original image, and v ′ is the l, α or β-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 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 µ ij and σ ij (resp. and µ ji and σ ji ) the mean and standard deviation of image i (resp. j) in the overlapping part [Eqn. [START_REF] Pan | Automatic generation of seamline network using area Voronoi diagrams with overlap[END_REF]]. We will also denote r ij the RMSE of images i and j in their overlapping part.

µ ij = E (v i ) and σ ij = E (v i -µ ij ) 2 and r ij = E (v i -v j ) 2 (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 µ ′ ij and σ ′ ij (resp. µ ′ ji and σ ′ ji ), and the RMSE is transformed to r ′ ij .

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 = 1 2 ⋅ S ij ⋅ µ ′ ij -µ ′ ji 2 (6)
For a set of n images, the total cost related to meanequalization is equal to the sum of pairwise costs:

C µ = 1 2 n i=1 n j=i S ij ⋅ µ ′ ij -µ ′ ji 2 (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 σ = 1 2 n i=1 n j=i S ij ⋅ σ ′ ij -σ ′ ji 2 (8)
The first cost function C µσ is obtained by summing mean-equalization cost and variance-equalization cost:

C µσ = C µ + C σ (9)

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 = 1 2 n i=1 n j=i S ij ⋅ r ′ ij 2 (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 µσ , or the RMSE-based cost function C r . We first rewrite cost functions C µσ 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.

Canonic formulation of cost functions

Let x denote the vector of parameters of the n color correction models: x = (a 1 , . . . , a n , b 1 , . . . , b n ) with a i and b i the slope and intercept coefficients related to image i. Cost functions C µσ and C r can be written as quadratic functions of parameter vector x (Appendix VII for detailed calculations):

C µσ (x) = 1 2 x ⊺ Q µσ x and C r (x) = 1 2 x ⊺ Q r x (11)
with

Q µσ = H µ + H σ K K ⊺ L and Q r = H r K K ⊺ L two 2n × 2n symmetric matrices, with H µ , H σ , H r , K
and L five n × n matrices defined by:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ H µ ij = -S ij µ ij µ ji if i ≠ j, and H µ ii = ∑ n j=1 j≠i S ij µ 2 ij H σ ij = -S ij σ ij σ ji if i ≠ j, and H σ ii = ∑ n j=1 j≠i S ij σ 2 ij H r ij = -S ij E[v i v j ] if i ≠ j, and H r ii = ∑ n j=1 j≠i S ij (µ 2 ij + σ 2 ij ) K ij = -S ij µ ij if i ≠ j, and K ii = ∑ n j=1 j≠i S ij µ ij L ij = -S ij if i ≠ j, and L ii = ∑ n j=1 j≠i S ij (12)
In the restricted case where the color correction models are linear transformations with zero yintercept [Eqn. (4b)], the vector of parameters is limited to a = (a 1 , . . . , a n ) with a i the slope coefficient related to the transformation of image i. Cost functions C µσ and C r are then equal to (Appendix VII):

C µσ (a) = 1 2 a ⊺ (H µ + H σ ) a and C r (a) = 1 2 a ⊺ H r a (13)

Equality constraints

To ensure that we find consistent and nondegenerated 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:

1 n n i=1 S ii µ ii = 1 n n i=1 S ii µ ′ ii (14)
Besides, it is also desirable that the mean spread of pixel values is preserved:

1 n n i=1 S ii σ ii = 1 n n i=1 S ii σ ′ ii (15)
These two constraints can be written in matrix notations (Appendix VIII for detailed calculations):

Ex = d (16)
with E a 2n × 2 matrix and d a 2-vector defined by:

E = S
= ⎛ ⎜ ⎜ ⎝ n ∑ i=1 S ii µ ii n ∑ i=1 S ii σ ii ⎞ ⎟ ⎟ ⎠
In the restricted case of linear color correction models with zero y-intercept [Eqn. (4b)], only the first n components of x are involved, and the equality constraint is given by:

E 0 a = d (17)
with E 0 a n × 2 matrix defined by: E 0 = S 11 µ 11 . . . S nn µ nn S 11 σ 11 . . . S nn σ nn

Problem solving

The problem of finding an optimal set of color correction models that minimizes the total cost function C µσ (x) or C r (x) can be formulated as a quadratic programming (QP) problem, that is, a problem of optimizing a quadratic function of several variables subject to linear constraints. The optimal set of parameters is given by:

x ⋆ = argmin x∈R 2n 1 2 x ⊺ Qx subject to Ex = d ( 18
)
with Q = Q µσ 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 a = (a 1 , . . . , a n ), and the QP problem is formulated as:

a ⋆ = argmin a∈R n 1 2 a ⊺ Ha subject to E 0 a = d (19)
with H = H µ + H σ or H = H r depending on which cost function is to be minimized.

As the total cost C µσ (x) or C r (x) is strictly positive for any set of non-harmonized set of images, we have x ⊺ Qx > 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 [START_REF] Gertz | Object-oriented software for quadratic programming[END_REF].

III. DATA

Our method for color harmonization in satellite image mosaicking was applied to a collection of high resolution Rapid Eye images, which was acquired over mainland France (6.7 × 10 5 km 2 ) between 2009 and 2010 within the EQUIPEX GEOSUD project [START_REF]Matrix values in Eqn. (1) slightly differ from those given by [19, Eqn (2)], which were found to be erroneous by an anonymous reviewer[END_REF]. 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 cutline network. Table I displays the main characteristics of the image collection, while Fig. 3 shows the spatial 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. Mean pixel count per image 6.9 × 10 8 pix.

Size of images files 550 Gb 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)

IV. RESULTS

The color harmonization method presented in section II was implemented in C++ using OTB [START_REF] Inglada | The Orfeo Toolbox remote sensing image processing software[END_REF], ITK [START_REF] Yoo | Engineering and algorithm design for an image processing API: A technical report on ITK -The Insight Toolkit[END_REF] and OOQP [START_REF] Gertz | Object-oriented software for quadratic programming[END_REF] libraries. For sake of clarity, we first analyse detailed results obtained on a single pair of partly overlapping images A and B (section IV A), then we display results related to the complete mosaic over mainland France (section IV B), 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 meanand-deviation-based cost function C µσ 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. [START_REF] Hall | Radiometric rectification: Toward a common radiometric response among multidate, multisensor images[END_REF]]: because image A is much larger than image B (S AA ≫ 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αβ 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αβ 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. [START_REF] Shiren | Two-dimensional seampoint searching in digital image mosaicking[END_REF].

Finally, Fig. 6 displays pixels distribution in lαβ 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 α and β channels-are closer after correction than before. Table II presents an analysis of the range of images A and B pixels values in lαβ 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αβ color space is reduced by [-92%, -97%, -100%] (l, α, β 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. Five different output mosaics were produced (Fig. 7), using either a linear transformation with zero yintercept [Eqn. (4b)] or a general linear transformation [Eqn. (4a)], and minimizing either the meanbased cost C µ , the mean-and-deviation-based cost C µσ , 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 appealling 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 ⋅ log 10 (peak 2 MSE) with peak = 2 8 ⋅ √ 3 for 8 bits RGB images, and MSE the mean squared error between pairs of overlapping images over the whole mosaic. Table III shows

#1 v ′ = a ⋅ v C µ 34.660 dB #2 v ′ = a ⋅ v C µσ 34.722 dB #3 v ′ = a ⋅ v + b C µσ 34.960 dB #4 v ′ = a ⋅ v C r 34.969 dB #5 v ′ = a ⋅ v + b C r 35.413 dB
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-anddeviation-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 compu-tation 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

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)
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.

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

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 implementation in OOQP library that takes advantage of sparse matrices. The run time is higher for a general linear color correction model [Eqn. (4a)], due to the fact that the corresponding QP has twice more unknown 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.

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 remote sensing images. The formulation we give of the color harmonization problem as a Quadratic Programming (QP) problem overcomes some of the limitations of existing techniques. It allows to deal simultaneously, rather than sequentially, with any number n of input images and any spatial layout of these images and their overlapping areas. The optimal parameters a and b of color correction models for all images are searched simultaneously in a 2n-dimensional space. As discussed in the introduction, this overcomes the limitation of color harmonization tech- 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. [START_REF] Shiren | Two-dimensional seampoint searching in digital image mosaicking[END_REF]]. More sophisticated and nonlinear 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 [START_REF] Reinhard | Color transfer between images[END_REF] 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 timeconsuming 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 neeeded to find extremal or even optimal values for this overlap ratio.

The tool was succesfully 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. [START_REF] Burt | A multiresolution spline with application to image mosaics[END_REF] and [START_REF] Szeliski | Image alignment and stitching: A tutorial[END_REF]], 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 [START_REF] Reinhard | Color transfer between images[END_REF]. Many other color transfer techniques could be investigated -see [START_REF] Faridul | A survey of color mapping and its applications[END_REF] for a review of such techniques. However, only linear correction models allow to formulate the multiimages 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αβ decorrelated color space, which is ideally suited to process Earth observation scenes. The method was succesfully 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. deviation σ ′ ij of image i in the overlap area with image j after correction is linked to that before correction by σ ′ ij = a i σ 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 µ we have:

C σ (x) = C σ (a) = 1 4 n i=1 n j=1 S ij (a i σ ij -a j σ ji ) 2 = 1 2 n i=1 n j=1 S ij (σ 2 ij a 2 i -σ ij σ ji a i a j )
Moving from the second to the third line above is obtained by noticing that

∑ n i,j=1 S ij σ 2 ij a 2 i = ∑ n i,j=1 S ij σ 2
ji a 2 j by swaping i and j indices in the double sum. We obtain:

C σ (x) = C σ (a) = 1 2 a ⊺ H σ a
with H σ a n × n matrix defined in Eqn. [START_REF] Victor | Automated image mosaicking[END_REF].

r ′ ij 2 = E[(v ′ i -v ′ j ) 2 ] = E[((a i v i + b i ) -(a j v j + b j )) 2 ] = a 2 i (σ 2 ij + µ 2 ij ) + 2a i b i µ ij + b 2 i + a 2 j (σ 2 ji + µ 2 ji ) + 2a j b j µ ji + b 2 j -2 ⋅ (a i a j E[v i v j ] + a i b j µ ij + a j b i µ ji + b i b j )
We develop the expression of C r :

C r (x) = 1 2 n i=1 n j=1 S ij ⋅ a 2 i (σ 2 ij + µ 2 ij ) + 2a i b i µ ij -2a i b j µ ij -a i a j E[v i v j ] + (b i -b j ) 2
Using the definition of the three n × n matrices H r , K and L given in Eqn. [START_REF] Victor | Automated image mosaicking[END_REF], we obtain:

C r (x) = 1 2 n i=1 n j=1
H r ij 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 (a, b) = 1 2 a b ⊺ H r K K ⊺ L a b
which writes:

C r (x) = 1 2 x ⊺ Q r x
with Q r = H r K K ⊺ L a 2n × 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: The second equality constraint is that the mean spread of pixel values over all images should be preserved by the color correction models:

1 n n i=1 S ii σ ii = 1 n n i=1 S ii σ ′ ii
The standard deviation value after image color correction is linked to the standard deviation value before correction by σ ′ ij = a i σ ij and we obtain: Therefore the two equality constraints (intensity preservation and spread preservation) can be written at once in matrix notations as shown in Eqn. [START_REF] Schott | Radiometric scene normalization using pseudoinvariant features[END_REF].

FIG. 2 :

 2 FIG.2: Overlapping area between two images i and j.

FIG. 4 :

 4 FIG. 4: Color harmonization for pair of overlapping images A and B.

FIG. 5 :

 5 FIG. 5: Corrected pixel values vs. original pixel values for both images A and B in RGB color space.

FIG. 6 :

 6 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.

FIG. 8 :FIG. 9 :

 89 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.

3 .S ij r ′ ij 2

 32 Mean-and-deviation-based cost function C µσ by summing the expression of mean-based cost function C µ and deviation-based cost function C σ , we find the expression of mean-and-deviation-based cost function C µσ :C µσ (x) = C µ (x) + C σ (x) = 1 2 x ⊺ Q µσ x with Q µσ = H µ + H σ K K ⊺ L a 2n × 2n symmetric matrix.B. Canonic formulation of cost function C rThe RMSE-based cost function C r writes: with r ′ ij the RMSE after transformation defined in Eqn. (5):

S

  ii µ ′ ii Using µ ′ ij = a i µ ij + b i with a i and b i the slope and intercept coefficients related to the image i, we obtain: n i=1 S ii µ ii = n i=1 S ii (µ ii a i + b i ) which writes: S 11 µ 11 , . . . , S nn µ nn , S 11 , . . . , S nn ⋅

S

  ii σ ii = n i=1 S ii σ ii a iwhich writes: S 11 σ 11 , . . . , S nn σ nn ⋅

  B. Conversion of RGB values to lαβ color spaceThe first step of color harmonization is to convert RGB values into a new color space called lαβ. This

	Input images	1	RGB to lαβ color space	Input images	in overlaps Image statistics
				QP problem
					solving
			2	Color	Optimal linear
				correction	correction models
				Corrected
				images
			3	Mosaick
				composition
	Output mosaic	4	color space lαβ to RGB	Output mosaic
	RGB color space			lαβ color space
			FIG. 1: Flowchart.

  11 µ 11 . . . S nn µ nn S 11 . . . S nn S 11 σ 11 . . . S nn σ nn 0 . . . 0 and d

TABLE I :

 I Characteristics of the image collection

	Characteristic	Value
	Satellite	RapidEye
	Acquisition dates	May to October 2009
		+ May to October 2010
	Azimuth angle	100 ○
	Image swath	Approximately 77 km
	Spatial reference system	EPSG 2154
	Spectral bands used	Red, Green, Blue
	Number of images	132
	Total spatial extent	2.3 × 10 6 km 2
	Pixel size (orthorectified)	5 m
	Pixel values encoding	16 bits
	Total input images pixel count 91 × 10 9 pix.

TABLE II :

 II Analysis of the range of values in lαβ space after color transformation

		Image A	Image B	Image A	Image B	Abs. Difference	Abs. Difference	Relative
		(original)	(original)	(corrected)	(corrected)	(original)	(corrected)	reduction %
	mean (l) 2.19E+00 2.31E+00 2.21E+00 2.22E+00	1.18E-01	9.30E-03	-92%
	mean (α) -6.01E-02 -5.85E-02 -6.03E-02 -6.02E-02	1.53E-03	4.10E-05	-97%
	mean (β) -1.59E-02 -2.24E-02 -1.75E-02 -1.75E-02	6.48E-03	1.50E-05	-100%
	std (l)	2.17E-01 1.72E-01 2.10E-01	1.98E-01	4.53E-02	1.24E-02	-73%
	std (α)	4.70E-02 3.74E-02 4.43E-02	4.85E-02	9.58E-03	4.28E-03	-55%
	std (β)	1.25E-02 1.07E-02 1.19E-02	1.30E-02	1.77E-03	1.07E-03	-40%

TABLE III :

 III Characteristics of France mosaics with quality metrics

	Mosaic id Correction model Cost function PSNR
	#0	none	none	32.948 dB

ACKNOWLEDGMENTS

This work was supported by public funds received through GEOSUD, a project (ANR-10-EQPX-20) of 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.

The mean value after image color correction is linked to the mean value before correction by µ ′ ij = a i µ 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 µ as a function of vector

We develop the expression:

By swaping i and j in the sums, we notice that

and we obtain:

Using the definition of the three n × n matrices H µ , K and L given in Eqn. [START_REF] Victor | Automated image mosaicking[END_REF], we obtain:

Using matrix notation, we get:

When the color correction model is a linear transformation with zero y-intercept, the vector of parameters is limited to a = (a 1 , . . . , a n ). In this restricted case, the cost function has the same expression as C µ with all b i = 0, and we directly obtain its final expression by putting b = 0 in Eqn. [START_REF] Ruderman | Statistics of cone responses to natural images: Implications for visual coding[END_REF]:

No matter which type of linear color correction model is used (Eqn. (4a) or Eqn.(4b)), the standard