F Soulez

É Thiébaut

J Kluska

M Tallon

A Schutz

A Ferrari

Direct temperature map estimation in optical long baseline interferometry

We present a method to recontruct temperature map directly from spectro-interferometric data. It uses a sparse coding method to describe each pixel as a blackbody. Results are shown on the Herbig Be HD98922. Aside from recovering the mean environment temperature of 1654 K in agreement with the photometry, this technique allows us to study the temperature distribution in the first astronomical units around the star.

INTRODUCTION

The arrival on 2 nd generation beam combiners GRAVITY and MATISSE at the VLTI sheds light on the need of multispectral image reconstruction algorithm. Indeed, as they combine four telescopes with medium and high spectral resolution, it is now conceivable to use them as high spatial resolution integral field spectrograph by mean of image reconstruction. Image reconstruction is a very challenging problem. Recently, as part of the effort of the community to develop the research on multispectral image reconstruction the POLCA project gaves birth to three methods : MiRA3D, 1 PAINTER 2 and SPARCO. [START_REF] Kluska | Sparco: a semi-parametric approach for image reconstruction of chromatic objectsapplication to young stellar objects[END_REF] Spatio-spectral data cubes give very valuable astrophysical information such as velocity, abundance or temperature map. It is important to notice that this astrophysical information represents often few estimated parameters per pixel (e.g. dispersion and mean velocity, temperature, relative abundance of few chemical elements,. . . ) and lies only on a complex but small subset of the spectrum. Using this fact as priors to better constraint the reconstruction, we propose a new algorithm that provides directly a 2D map of parameters of astrophysical interest. In interferometry, such an idea was already implemented using CLEAN in the OYSTER package. [START_REF] Hummel | A novel imaging algorithm for broadband aperture synthesis data[END_REF] It is also behind the principle SPARCO 3 method that only estimates from interferometric data a 2D image and a parametric model with different spectral indexes. In the presented work, we describe a method to reconstruct the temperature map from the polychromatic interferometric measurements.

DICTIONARY BASED POLYCHROMATIC IMAGE RECONSTRUCTION

The interferometric image reconstruction consists in recovering the polychromatic object x given the interferometric measurements (i.e. in most cases: squared visibilities and phase closures). This object is discretized as a vector of K pixels and L wavelengths: x = (w 1,1 , . . . , w K,L ). Most algorithm solves this problem as a constrained minimization:

x

+ = arg min x≥0 (L(x) + µ R(x)) , s.t. ∀ , k x k, = 1. ( 1
)
where L is a likelihood term ensuring that the model is in agreement with the data, R is a regularization function that enforces some priors and µ ≥ 0 an hyper-parameter. The constraints k x k, = 1 ensure the spatial normalization of x as the visibilities are intrinsically normalized in each spectral channel . In our approach, we assume that the spectrum x k of each pixel k can be described by a linear combination of few spectra given by the dictionary D such that x k = D α. In that case, the Equation 1 becomes a sparse coding problem:

5 x + = arg min x L(x) + γ k α k 1 s.t.    ∀ , k x k, = 1. x ≥ 0 ∀k, s × x k = D α k , (2)
The 1 norm of α k acts as the regularization function R and is used to enforce its sparsity (i.e. most elements of α k are null). The level of sparsity is balanced by the hyper-parameter γ ≥ 0. As all spectral channels are normalized, the information about the relative flux between each channel is lost. Hence, each spectrum has to be rescaled by the spectrum s of the whole object before its decomposition on the dictionary. Following the framework developed for the MiRA3D 1 and PAINTER 2 algorithm, we solve Equation 2 using the Alternating Direction Method of Multipliers (ADMM).

RESULTS ON HD98922

Optical interferometry is the only mean to observe the central astronomical unit of the young objects. The study of the phenomena at this scale was the aim of a large observational program on 55 Herbig AeBe with VLTI/PIONIER conducted by J.-P. Berger in 2012/2013 (ID: 090.C-0963). 8 Among these stars, the observation of HD98922 has one of the best (u, v) coverage (see Figure 1). HD98922 is an interesting group II Herbig Be star having an SED suggesting a full disc without evidence for a gap. Its type is however not precisely known as indicated on Table 1. It was observed in the first two months of 2013 with the four 1.8 m Auxiliary Telescopes (ATs) of the Very Large Telescope using PIONIER 9 (see table 3), with using the low spectral resolution in H band with 3 spectral channels ()R ≈ 15).

Building up the dictionary

We assume that the HD98922 SED is dominated by thermal emission and that its circumstellar disk can be modeled as a blackbody with spatially varying temperature. The dictionary D is thus built from a set of 500 spectra of black bodies ranging from 500 K to 20000 K.

No simultaneous spectrum of HD98922 was available to set the scaling vector s. To overcome that issue, we estimated the spectrum s such as the star (i.e. central pixel) had the theoretical spectrum given by its spectral type. Although the type of HD98922 is not precisely known, both types indicated in Table 1 (B9V and A5III) have a quasi identical spectra in H band. Consequently, we have made a reconstruction using the Kurucz spectrum 10 for a B9V of 10500 K with [F e/H] = 0 and log g = 4.

Temperature map reconstruction

A reconstruction (pixel size: 0.25 mas) was performed with γ = 10 -6 . After 200 iterations, the final reduced χ 2 was ≈ 1.2. The three reconstructed channels are shown Figure 2 to Figure 4.

Due to the sparsity constraint, for each pixel k the parameter α k have only very few non-zero elements. In practice, we noticed that for each pixel k, only one element was significativelly non-zero. As each column of the dictionary D are spectra of black bodies, the index of this non-zero element of α k gives directly a good estimate of the temperature of the object at pixel k. We are thus able to recover the temperature map in the first astronomical units around the star directly from the reconstruction without further processing. This temperature map is given on Figure 5. The disk parameters retrieved from this reconstruction are shown on Table 3. On the reconstructed image (Figure 2 to Figure 4) the disk is well resolved and its inner rim is clearly visible. From the temperature map on Figure 5 we estimated the relative flux of the disk of 0.66 and its mean temperature of 1654 K. That is consistent with a relative flux of 0.66 and a temperature of 1493 K estimated from photometric measurements. [START_REF] Lazareff | Few puffed-up rims found in large interferometric survey of herbig aebe stars[END_REF] In addition, some flux seems to be emitted from inside the cavity especially in the bluer part of the spectrum.

A reconstruction at such scale raises several new questions. Taking into account the geometry, one of the sides should be less bright. For the shortest wavelengths, it seems also that the flux emitted from inside the cavity may come from a gaseous internal disk.

CONCLUSION

This work illustrates perfectly the interest of the image reconstruction in interferometry for observing precisely the close environments of YSOs (≈ 1 UA). It provides a new way to estimate quantities of astrophysical interest (here the temperature map) directly from the data. It can be viewed as a generalization of the SPARCO method 3 with as much black-bodies as pixels. It is simpler and more robust than a classical image reconstruction followed by temperature estimation that can introduce additional errors and biases.

More generally it illustrates the potential of approaches based on dictionaries for all the hyperspectral instruments like integral-field spectrographs. Indeed, introducing into the dictionary spectra coming from different models will allow to rapidly discriminate the most probable model.

  Figure 1. (u, v) coverage of the observations. The three colors correspond to the PIONIER spectral channels.

Figure 2 . 1 .

 21 Figure 2. 1.59 µm channel of HD98922 reconstruction Figure 3. 1.68 µm channel of HD98922 reconstruction

Figure 4 . 1 .

 41 Figure 4. 1.77 µm channel of HD98922 reconstruction

  diameter (major axis) inclination log(T env ) (K) T env disk relative flux 4

Figure 5 .

 5 Figure 5. Recovered temperature map of HD98922

Table 1 .

 1 Possible types of HD98922

	Type Distance	T eff	ref.
	B9Ve 1150 930 360 A2III 507	10500 ± 500 Alecian et al. 6 9000 Hales et al. 7

Table 2 .

 2 HD98922 observation log.

Table 3 .

 3 Estimated parameters of the disk

Acknowledgements

This work is supported by the Sinergia project "Euclid: precision cosmology in the dark sector" from the Swiss National Science Foundation and by the French ANR POLCA project (Processing of pOLychromatic interferometriC data for Astrophysics, ANR-10-BLAN-0511).


Matthew Ozon

Éric Thiébaut
email: eric.thiebaut@univ-lyon1.fr

Ferréol Soulez

Jonathan Léger

Global optimization for image restoration in optical/IR interferometry

Keywords: Stochastic global optimization, interferometric imaging, image reconstruction, phase retrieval

Many image reconstruction algorithms have been proposed to cope with optical/IR interferometric data. A first difficulty is to deal with the sparsity of the data and the uneven u-v coverage. This is usually overcome by imposing a priori constraints to the sought image. The reconstruction then amounts to solving an optimization problem where fidelity to the data and to the priors must both be satisfied to a chosen level. Due to the type of interferometric measurements, the image restored by the current algorithms depends on some initial guess and on the optimization strategy. Even though the resulting image is often satisfactory, there are no guarantees that the best image is obtained given the data and the constraints. This is a strong defect which may badly impact the interpretation of interferometric observations by astronomers. In this paper, we consider using stochastic methods such as simulated annealing to solve for the problem of image reconstruction with sparse Fourier amplitudes and unknown Fourier phases. This preliminary study is a first step toward image reconstruction with partial Fourier phase information.

IMAGE RECONSTRUCTION WITH UNKNOWN FOURIER PHASES

Ideally, interferometric measurements consist in samples of the Fourier spectrum of the object brightness distribution. The k-th datum then writes:

d k = I λ (ν k ) + n k ( 1
)
where I λ (ν k ) is the Fourier transform at spatial frequency ν k of the specific brightness distribution I λ (θ) of the object in the angular direction θ (at wavelength λ) and n k is a nuisance term which accounts for the noise. In what follows, we will assume that there are M sampled spatial frequencies.

In order to perform image reconstruction, we must choose a parametric representation of the object brightness distribution. For instance, [START_REF] Thiébaut | Image reconstruction with optical interferometers[END_REF] we assume that:

I λ (θ) = N j=1 b(θ -θ j ) x j , ( 2
)
where b : Θ → R is some basis function over Θ the synthesized field of view, x = (x 1 , . . . , x N ) t ∈ R N are the image parameters and {θ j ∈ Θ} j=1,...,N is a rectangular grid of angular directions over Θ. If the b j are cardinal interpolation functions, x is a conventional pixelized image whose components are x j ≈ I λ (θ j ). The Fourier transform of the image model writes:

I λ (ν) = N j=1 b(ν) exp(-2 i π ν t θ j ) x j , (3)
with b(ν) the Fourier transform of b(θ). Then the model of the data in Eq. (1) becomes:

d k = N j=1 H k,j x j + n k (∀k = 1, . . . , M ) ⇐⇒ d = H x + n ( 4
)
where the right hand side expression makes use of the matrix notation, M is the number of measurements and the coefficients of the operator H ∈ C M ×N are given by:

H k,j = b(ν k ) exp(-2 i π θ t j ν k ) . ( 5
)
When the complex Fourier spectrum is fully measured, image reconstruction amounts to solving a constrained optimization problem, 2 typically:

x = arg min x≥0 H x -d 2 W + µ g(x) (6)
where H x -d 2 W is a data fidelity term, g(x) is a regularization term which implements prior knowledge about the sought image x and µ > 0 is a tuning parameter which controls the trade-off between agreement to the data and to the priors. The constraint x ≥ 0 holds componentwise and means that the sought image must be nonnegative everywhere. The notation H x -d 2 W suggests that the noise n follows a Gaussian distribution and thus maximum likelihood leads to minimize some weighted quadratic norm of the difference between the data and their model. Provided g(x) is convex, the solution to the Problem ( 6) is unique and can be efficiently computed by existing algorithms.

With optical interferometers, the situation is not so simple because not all Fourier information can be measured for the sampled frequencies {ν k } k=1,...,M . In particular, due to the turbulence, Fourier phase can only be partially measured by, for instance, so-called phase closure. The image reconstruction problem is then nonconvex and different strategies have been devised to cope with these issues. [START_REF] Meimon | Reconstruction method for weak-phase optical interferometry[END_REF][START_REF] Ireland | Monte-Carlo imaging for optical interferometry[END_REF] In this paper we consider various stochastic methods to find the global optimal solution of the problem and compare these methods to other approaches implemented by existing image reconstruction algorithms. For the sake of generality, although this is not the easiest problem to solve, we assume that only the Fourier amplitudes have been measured and the image reconstruction is closely related to phase restoration problems commonly encountered in other domains as in crystallography. Assuming that measures at different frequencies are independent, the easy image reconstruction Problem (6) becomes a joint image and phase restoration problem:

x, φ ∈ arg min

(x,φ)∈X×Φ f X×Φ (x, φ) = k w k |(H x) k -ρ k e i φ k | 2 + µ g(x) , ( 7
)
where ρ k are the measured amplitudes (we assume that ρ k ≥ 0 whatever k), w k ≥ 0 are given statistical weights and φ = (φ 1 , . . . , φ M ) t ∈ Φ = (-π, +π] M are unknown Fourier phases and X = {x ∈ R N | x ≥ 0} is the set of nonnegative images. In the following we consider different approaches to solve Problem (7). Note that the solution may not be unique.

Self Calibration Approach

Observing that f (x, φ) given in Eq. ( 7) is convex with respect to x knowing φ and that finding the best phase given x has a simple solution, a first possibility is to solve the joint problem by the following alternating method: 0. Choose an initial image x (0) , then for t = 1, 2, . . . and until convergence repeat:

1. Phase restoration step (self-calibration):

φ (t) = arg min φ f X×Φ x (t-1) , φ = arg min φ k w k H x (t-1) k -ρ k e i φ k 2 , ( 8
)
the problem is separable and, ∀k, the trivial solution is φ

(t) k = ∠ H x (t-1)
k where ∠z yields the phase of z ∈ C and provided ρ k > 0.

2. Image restoration step:

x (t) = arg min x≥0 f X×Φ x, φ (t) = arg min x≥0 k w k (H x) k -ρ k e i φ (t) k 2 + µ g(x) , (9)
which amounts to solving a strictly convex problem and has a unique solution.

The above algorithm is similar to the strategy implemented by so-called self-calibration methods 5 for radiointerferometry or more recently by Wisard [START_REF] Meimon | Reconstruction method for weak-phase optical interferometry[END_REF][START_REF] Meimon | Convex approximation to the likelihood criterion for aperture synthesis imaging[END_REF] for optical interferometry. The main difference with these latter approaches is that rigorously the same criterion f (x, φ) is minimized by the alternating method. It is possible to generalize this method to the other kinds of data that can be provided by optical interferometers (not just the Fourier amplitudes). 7

Hierarchical Optimization

As seen in the previous section, finding the best unknown Fourier phases φ given the image x has a closed form solution which is easy to compute. This suggests to get rid of these phases by replacing them by their optimal values given the image parameters x. This yields a criterion f X (x) = min φ f X×Φ (x, φ) which only depends on the image parameters. The image restoration problem then writes:

x = arg min x∈X f X (x) = min φ∈Φ f X×Φ (x, φ) = k w k (|H x| k -ρ k ) 2 + µ g(x) . ( 10
)
The criterion f X (x) is however non-convex and global optimization is de facto required to find its global minimum.

Nevertheless, monotonic optimization algorithms may be able to find a good local minimum. This is the approach implemented by Bsmem 8, 9 and Mira 10 which, starting from an initial image, perform a local nonlinear constrained optimization with algorithm MemSys [START_REF] Skilling | Classic maximum entropy[END_REF][START_REF] Gull | Developments in maximum entropy data analysis[END_REF] for Bsmem and Vmlmb 13 for Mira. The result then depends on the initial image which may be randomly chosen.

Stochastic Methods

The size of an image reconstruction problem forbids to apply a systematic search in the parameter space. Among the algorithms proposed to solve global optimization problems of large dimensions, stochastic methods like simulated annealing (SA) have some desirable characteristics. For short, the main advantage is that a theoretical framework proves the convergence and gives hint on how to tune the algorithm so that the convergence is almost sure. Furthermore, while the brute force algorithms grow exponentially with the dimensionality, SA only grows linearly. This motivate the use of SA in the presented study. Note that stochastic methods are already exploited by existing algorithms. For example, algorithms Macim 4 or Squeeze 14 applied to the problem in Eq. ( 10) would attempt to find the global optimum by means of a simulated annealing method or Monte Carlo Markov Chain (MCMC) method.

Applying a stochastic method to the variables x and find the global minimum of f X (x) = min φ∈Φ f X×Φ (x, φ) may be prohibitively slow because there is a lot of variables (as many as the number N of pixels). As f X×Φ (x, φ) in x, the function:

f Φ (φ) = min x∈X f X×Φ (x, φ) = f X×Φ ( x(φ), φ) , (11)
where:

x(φ) = arg min x∈X f X×Φ (x, φ) = k w k |(H x) k -ρ k e i φ k | 2 + µ g(x) , (12)
is properly defined and only depend on φ. The idea is then to apply the simulated annealing method to the M unknown phases which, in optical interferometry, is usually much smaller than the number N of pixels. For each proposed phase, an image reconstruction problem has to be solved and is likely to be the largest consumer of computational time. So the computational burden scales as ∼ t img (N ) t glob (M ) where t img (N ) represents the time spent to reconstruct a N -pixel image while t glob (M ) represents the time spent by a simulated annealing strategy to solve a global optimization problem with M variables. As the cost of the image reconstruction may be too important, a last possibility is to stochastically sample the whole set of parameters, x and φ, with the simulated annealing method. Then the time spent should scale as ∼ t glob (M + N ).

SIMULATED ANNEALING

We consider the problem of finding a global minimum of a real valued multi-variate function:

ω ∈ arg min ω∈Ω f Ω (ω) , (13)
where, depending on the chosen problem, the parameters ω and the corresponding space Ω can be ω = x and Ω = X for f Ω = f X defined in Eq. ( 10), ω = φ and Ω = Φ for f Ω = f Φ defined in Eq. ( 11) or ω = (x, φ) and Ω = X × Φ. Whatever the choice of the parameter space Ω, the criterion f Ω (ω) is non-convex in ω and monotonic optimization algorithms do not apply and we propose to Simulated Annealing (SA) to sample the global minima of f Ω (ω). However, theoretical results ensuring the convergence of SA have only been derived for problems defined over finite parameter spaces. We will therefore assume that the parameters ω ∈ Ω take a finite number of possible values, possibly after proper discretization.

The family of SA algorithms is well suited for large combinatorial optimization and guarantees to reach a state that belongs to the set of global optimum given a tolerance margin. Contrary to monotonic descent algorithms, SA allows unfortunate configurations to be temporarily accepted in order to escape local optimal configurations (similar to metastable states during annealing in metallurgy [START_REF] Silliman | Principles of Physics, Or Natural Philosophy[END_REF][START_REF] Ewing | The crystalline structure of metals. second paper[END_REF] ) and eventually reach a global minimum. The controlled acceptance of uphill configurations -those that increase the energy -is the relevant distinction from deterministic algorithms or other stochastic samplers.

General Implementation

The SA is an iterative procedure that generates a sequence of states ω (t) ∈ Ω. A typical implementation is given by Algorithm 1 whose output is a sequence (ω (t) ) 0≤t≤N SA |Ω| which is a Markov chain (MC). This chain is defined by three elements: the energy function f Ω , a temperature schedule (T n ) 1≤n≤N SA , and a communication kernel q(ω | ω , T ) that may depend on the temperature T and which gives the probability to generate a state ω given a previous state ω and a temperature T . More specifically, SA draws every candidate ω (but the first one) at random as ω ∼ q(ω | ω , T ) with ω the previously accepted candidate and T the current temperature. The kernel can depend on the temperature so that the state space is largely sampled at the beginning of the process and narrow down with the iterations. In other words, at the beginning of the SA process, the kernel may generate large steps in order to visit a large regions of the landscape while at the end only small steps can be proposed. This way, the convergence speeds up in terms of iterations.

Pick an initial state ω (0) ∈ Ω at random. Set the iterator t = 0. for n = 1, . . . , N SA do Update temperature: set T = T n . for j = 1, 2, . . . , |Ω| do Draw a candidate ω from the probability distribution q(ω | ω (t) , T ).

Set ∆ = f Ω (ω) -f Ω ω (t) .
Set ω (t+1) = ω with probability min 1, e -∆/T , otherwise set ω (t+1) = ω (t) . Set t = t + 1.

Algorithm 1: Simulated Annealing minimizing the energy function f Ω (ω).

The final state ω = ω (N SA |Ω|) of the sequence produced by Algorithm 1 belongs to the set of global minimizers of f Ω under weak assumptions that are discussed in Section 2.2. To put things clear, the convergence of SA is not defined in terms of distance to the solution, which would be a non-sense for the function having more than one global minimizer. Rather, the convergence properties arise from the probability of belonging to the set of minimizers. Hence, considering a margin of error > 0, the probability that the last element ω = ω (N SA |Ω|) of the Markov chain does not belong to the set of global optimizers is given by:

Pr f Ω (ω ) ≥ min ω∈Ω f Ω (ω) + ≤ (N SA |Ω|) -α , ( 14
)
for some α > 0. Clearly the convergence speed depends on α.

It is interesting to notice that a global minimum may not be unique and the Algorithm 1 is a uniform sampler over the set of global minima. Having several minima is not an aberration, but merely the representation of the complexity of the reconstruction problem that may accept several good solutions.

Assumptions

Fortunately, the assumptions ensuring the convergence * of an SA process in a finite number of iterations are weak. The requirements concerning the temperature, the communication kernels and the energy function are few and not too restrictive. First, the energy function f Ω must be defined over a finite state space, yet no other constraints add up. Then, the family of communication kernels q(ω | ω , T ) and its limit while T → 0 must satisfy the three following assumptions:

Irreducibility. The irreducibility expresses the idea that at a fixed temperature of the cooling schedule, or at least before the end of the annealing, all states are reachable from another state in a finite number of steps. In other words, there exists a path (ω (t) ) 0≤t≤n starting at any state ω (0) ∈ Ω and reaching any other ω (n) ∈ Ω in a finite number of steps n ∈ N * < ∞ with a non null probability:

∀t ∈ [[1, n]], q ω (t) ω (t-1) , T > 0 . ( 15
)
Symmetry of the support. This property means that it is always possible to turn back. If the transition from ω 1 ∈ Ω to ω 2 ∈ Ω is possible, that is to say q(ω 2 | ω 1 , T ) is not null, then the transition from ω 2 to ω 1 must also be possible which implies that q(ω 1 | ω 2 , T ) is not null either. The two probabilities of transition are however not necessary the same. Formally, the following must stand:

q(ω 2 | ω 1 , T ) = 0 =⇒ q(ω 1 | ω 2 , T ) = 0 , (16)
for any pair (ω 1 , ω 2 ) ∈ Ω 2 . It implies that the support ∆q of the communication kernel q:

∆q = (ω 1 , ω 2 ) ∈ Ω 2 q(ω 1 | ω 2 , T ) > 0 , (17)
is symmetric.

Equality of the supports. During the cooling process, the support of q may depends on the temperature T , however, starting from a temperature T > T end , the support of the communication kernels must be the same. Formally, the following must stand:

q(ω 2 | ω 1 , T ) = 0 =⇒ q(ω 2 | ω 1 , T ) = 0 , ( 18
)
for any pair (ω 1 , ω 2 ) ∈ Ω 2 and any temperatures T ≤ T and T ≤ T .

Temperature estimation

Since the early 90's, the SA process has been proven to converge in a finite number of iterations when using exponential piecewise constant decreasing temperature sequences of the form:

17 ∀n ∈ [[0, N SA ]], T n = T init (T end /T init ) n/N SA , ( 19
)
in Algorithm 1. The estimation of the temperature parameters T init and T end is based on the average acceptance rate χ T of uphill moves for a given temperature T and computed for two chains of states that increase the energy:

χ T = 1 K K t=1 e -∆ (t) T /T , ( 20
)
where the sequence (∆

(t)
T ) 0≤t≤K is the energy variation induced by the uphill chain, and K is its length. In order to have a good exploration of the energy landscape, an SA process should be characterized by the uphill acceptance rates in [0.6, 0.9] at the beginning and in [10 -5 , 10 -3 ] at the end. Furthermore, it is desirable that χ T decrease linearly with the iterations. We noticed that it is interesting to finish the process with a few extra steps ( 10) at null temperature so that the algorithm does not finish after uphill moves.

Definition of the communication kernel

The idea of the communication kernel is to sample the state space in a clever way. Indeed, as it is mentioned before, the SA process should be able to visit a large panel of states so that it does not get trapped by local optima attraction basins, while at the end of the process, only small perturbations are required in order to overcome the frictions due to the ruggedness of the energy landscape. Hence a good communication kernel must change during the process. Therefore, we propose two elementary kernels that depend on T and form the base of communication kernels used in this study. One, q Φ , explores the phase space and the other, q X , acts over the pixel space. Three kernels are used in our study: q Φ is used for SA on the phases, q X is used for SA on the pixels and a composite kernel q X×Φ is used for Full SA Reconstruction (FSAR), that is to say SA on both the pixels and the phases:

q X×Φ (ω | ω , T ) = τ q Φ (φ | φ , T ) + (1 -τ ) q X (x | x , T ) . ( 21
)
with ω = (x, φ) and ω = (x , φ ). The scalar τ ∈ [0, 1] is the probability of choosing the phase space over the pixel space. The dependence on the temperature in ( 21) is only supported by the two elementary kernels, though τ could vary with the iterations as a way of speeding up the convergence. For convenience, and because the proposed kernels depend linearly on the iteration, we define the normalized index ζ(T ) = log(T /T init )/ log(T end /T init ) = n/N SA that linearly varies in the range [0, 1] with the iterations.

Pixel kernel

This kernel aims at modifying the scene x ∈ X that is being reconstructed. For the sake of simplicity, our pixel kernel q X (x | x , T ) implements the most elementary moves. Considering any given state x ∈ X, the proposed state x is generated by the following procedure:

1) draw a pixel index j 0 ∼ U([[1, N ]]); 2) define the maximum amplitude δ x (T ) = δ max x -(δ max x -1) ζ(T ) ∈ [[1, δ max x ]]; 3) draw the level of variation ∼ U([[-δ x (T ), +δ x (T )]]); 4) ∀j ∈ [[1, N ]], set x j = max{0, min{1, x j + /(L x -1)}} if j = j 0 , x j otherwise;
where L x ∈ N * is the number of discrete values in the range [0, 1] of allowed pixel intensities and δ max x is the maximum amplitude of variation. The amplitude δ x (T ) decreases linearly with the iterations from δ max x at T init to δ min x = 1 at T end . In practice, we take δ max x = (L x -1)/2. Such a kernel takes large steps at the beginning and only small ones at the end so that the energy landscape is vastly explored at the beginning and refined at the end. Note that q X enforces the constraints of belonging to the closed set X, i.e. the positivity and the bounded nature of the brightness.

Phase kernel

The proposed phase kernel implements simple atomic modifications in the phase space. One of the phase is picked at random over all the frequencies but the null one (which is assumed to be at index k = 1), then an additive variation is applied. The frequencies are sampled following the distribution L(ρ) that linearly depends on the amplitude of the complex visibilities. The idea is to visit more frequently the phases of larger influence than those inducing small energy variation. Therefore, the main components are well reconstructed while the details are neglected. Our phase kernel is implemented as follows:

1) draw a phase index k 0 ∼ L(ρ);

2) define the maximum amplitude δ φ (T ) ∈ [[1, δ max φ ]]; 3) draw the level of variation ∼ U([[-δ φ (T ), +δ φ (T )]]); 4) ∀k ∈ [[1, M ]], set φ k = φ k + (2 π/ρ k ) ( /L φ ) mod 2 π if k = k 0 , φ k otherwise;
where L φ ∈ N * is the number of discrete values in the range [-π/ρ k , +π/ρ k ] of phase variation and δ max φ is the maximum amplitude of variation.

RESULTS

We apply the proposed SA variants to simulated data.

Simulated data

We consider the object of the 2014 Interferometric Imaging Beauty Contest. [START_REF] Lawson | An interferometry imaging beauty contest[END_REF] We simulated noisy data with independent Gaussian noise with given signal to noise ratio (SNR) as follows:

d k = N j=1 H k,j x true j + n k (∀k = 1, . . . , M ), (22)
with x true j the true object displayed in Fig. 1 and independent Gaussian noise such that Re(n k ) ∼ N (0, σ 2 k ) and Im(n k ) ∼ N (0, σ 2 k ). To ensure a given SNR, the standard deviation σ k of the noise at k-th frequency is given by (see Appendix A):

σ k = max(|z k |, ρ min ) 2 SNR 2 + 2 SNR 1 + SNR 2 ≈ SNR 1 max(|z k |, ρ min ) 2 SNR , ( 23
)
with z k = N j=1 H k,j x true j and ρ min > 0. For our simulations, we consider very good data and take SNR = 100 and ρ min = 10 -2 . The number of sampled frequencies (see Fig. 1b) is M = 196 including the zero-th frequency.

Common Settings

For the regularization we use the loose support constraint [START_REF] Le Besnerais | Advanced imaging methods for long-baseline optical interferometry[END_REF] and tune the regularization level µ so that the best image (by visual assessment) can be reconstructed from the complete data (Fourier amplitudes and phases) by solving Problem (6). The corresponding image is shown by Fig. 3a.

In order to minimize the number of parameters, we take the largest pixel size as possible (because it yields small images and hence fast computations in all respect) while keeping enough degrees of freedom. To determine the largest pixel size, we proceed as follows. We perform a convex image restoration with the complete data (Fourier amplitudes and phases) by solving Problem (6) with a given coarse pixel size. We then perform a second reconstruction, again solving Problem (6), with a much smaller pixel from the same simulated Fourier amplitudes but Fourier phases computed from the restored coarse image and check that this second image compares well with the best one shown in Fig. 3a. The largest pixel size selected in this manner leads us to apply the SA methods to 14 × 14 = 196 pixel images.

SA on the Phases

For reconstructions based on SA on the phases, the phase kernel q Φ is used with L φ = 180 and δ φ (T ) = 2 if ζ(T ) < 1/2 and δ φ (T ) = 1/2 for cooler temperatures, so that the angular precision is ≈ 2 degrees.The total number of SA iterations is set to N SA = 5000 which is roughly 5000 times the number of iterations for the algorithm used for image reconstructions.

For the image reconstructions, we use VMLMB, [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF] a quasi-Newton method with bound constraints, to carry the constrained optimization of the convex criterion defined in (12). VMLMB is used with a maximum of 20 iterations and enforcing the pixel values of being in the set X = [0, 1] N . To speed up the computations, each reconstruction for a new candidate is started with the image obtained for the previous accepted candidate.

SA on the Pixels

For reconstructions based on SA on the pixels, the pixel kernel q X is used with L φ = 50/ max(x true ) with x true the true image. The total number of SA iterations is set to N SA = 200000.

SA on the Pixels and on the Phases (FSAR)

For Full SA Reconstruction (FSAR) which applies SA on both the pixels and the phases, the composite kernel q X×Φ is used with τ = 1/2. The settings for the pixel kernel q X and the phase kernel q Φ are the same as above. The total number of SA iterations is set to N SA = 50000.

Analysis of the Convergence

Figure 2 shows the evolution of the acceptance rate and the energy during the SA iterations of the considered problems. A few remarks can be made. (i) The uphill move acceptance rate varies from 0.8 ∼ 0.9 to ∼ 10 -3 which indicates that at the beginning, the space is widely visited, then the criterion of minimizing the energy is more and more enforced, and, at the end, only few uphill moves are accepted so that the minimization terminates without being stuck because of the roughness of the energy surface. (ii) The energy tends to decrease fast at the beginning, then it varies slowly. The convergence is clear, but it is not obviously in the absolute minimum unfortunately. (iii) Two interesting features of the figures 2b, 2c, 2e and 2f: the discontinuity indicates the change in precision of the phase kernel (which is not visible for the pixel kernel because the precision varies slowly). As a consequence, there are almost two SA process in one, such as a multiscale approach. First a coarse convergence then the refinement.

Analysis of the Estimated Parameters

The estimated phases and restored images are shown in Figure 3. Sub-figures 3a and 3e shows the best results we can hope for, considering the reconstruction model, the regularization parameter and the optimization algorithm. The phases are organized by increasing visibility magnitude. It is clear that they are organized as groups (which should be integrated in the communication kernel in order to improve the quality of the result). Sub-figures 3b and 3f show the result when applying SA on the pixels: only the very first phase are well organized, as a consequence the reconstruction is poor. Clearly the method was not successful with the allowed number of SA iterations. By contrast, the phase restored by applying SA on the phases (Fig. 3c) or on both the pixels and the phases (Fig. 3d) are closer to the actual phases. As a result the restored images (Fig. 3g and Fig. 3h) are better approximations of the truth. As we expected, the phase and the image obtained by applying SA on the smaller number of parameters, that is on the phases with Ω = Φ, yields the best results. However applying SA on the phases and on the pixels is not so bad.

DISCUSSION AND PERSPECTIVES

The paper presents a preliminary study of exploiting Simulated Annealing (SA) to solve the difficult problem of image restoration from sparse measurements of the Fourier amplitudes. This also amounts to solving a Fourier phase restoration problem. Our results show that applying SA on the unknown phase coupled with a convex image reconstruction is the most effective of the considered stochastic variants. In particular, applying SA to the pixels rather than to the phases is easier (because there are no needs to implement a convex image reconstruction algorithm) but yields much worse results. This preliminary study is a first step toward image reconstruction with partial Fourier phase information which will be considered in a following work. (28)

To generate artificial complex-value data with a given SNR (for the powerspectrum), the variance must be:

σ 2 = |E{ a}| 2 SNR( a) 2 + 2 SNR( a) 1 + SNR( a) 2 . ( 29
)
As lim |E{z}|→0 σ = 0, we use a minimal value for σ.

Figure 1 :

 1 Figure 1: Simulated data for the experiments: a) the true object, 18 b) Fourier plan coverage and c) dirty beam.

  (a) Acceptance rate for Ω = X. (b) Acceptance rate for Ω = Φ. (c) Acceptance rate for Ω = X × Φ. (d) Evolution of f Ω for Ω = X. (e) Evolution of f Ω for Ω = Φ. (f) Evolution of f Ω for Ω = X × Φ.

Figure 2 :

 2 Figure 2: Evolution of the acceptance rate and the energy with the iterations during the SA process for the considered variants.

Figure 3 :

 3 Figure 3: Phase estimations (top row) and image reconstructions (bottom row). The phases are sorted with respect to the amplitude of the true complex visibilities. (e) is the best image that can be reconstructed from the complete data (Fourier amplitudes and phases) by solving Problem (6); (f) is the final image with SA on the pixels, i.e. with Ω = X; (g) is the final image with SA on the phases, i.e. with Ω = Φ; (h) is the final image with SA on the pixels and on the phases, i.e. with Ω = X × Φ.

κ = 2 2 4 σ 4 + 4 |E{z}| 2 σ 2 .

 222 degrees of freedom and parameter λ = E{Re(z/σ)} 2 + E{Im(z/σ)} 2 = |E{z/σ}| 2 . Expectation and variance of a/σ 2 are given by:E{a/σ 2 } = κ + λ , (24) Var{a/σ 2 } = 2 κ + 4 λ ,(25)Hence:E{a} = (κ + λ) σ 2 = |E{z}| 2 + 2 σ 2 , (26) Var{a} = (2 κ + 4 λ) σ 4 = 4 σ 4 + 4 |E{z}| 2 σ 2 , (27)Therefore a = a -2 σ 2 is an unbiased estimator of the power spectrum and its signal to noise ratio (SNR) is:

Proc. of SPIE Vol. 9907 99071E-2 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

Proc. of SPIE Vol. 9907 99071E-3 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

Proc. of SPIE Vol. 9907 99071E-4 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

* in terms of probability of belonging to the set of minimizers Proc. of SPIE Vol. 9907 99071E-5 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

Proc. of SPIE Vol. 9907 99071E-6 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

Proc. of SPIE Vol. 9907 99071E-7 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

Proc. of SPIE Vol. 9907 99071E-8 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

Proc. of SPIE Vol. 9907 99071E-9 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

Proc. of SPIE Vol. 9907 99071E-10 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

Proc. of SPIE Vol. 9907 99071E-11 Downloaded From: http://proceedings.spiedigitallibrary.org/ on 08/30/2016 Terms of Use: http://spiedigitallibrary.org/ss/termsofuse.aspx

ACKNOWLEDGMENTS

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2013-2016) under grant agreement number 312430 (OPTICON).

APPENDIX A. COMPUTATION OF THE SIGNAL TO NOISE RATIO

Consider the powerspectrum a = |z| 2 of z ∈ C with i.i.d. noise in the real and imaginary parts, i.e. Re(z) ∼ N (Re(E{z}), σ 2 ) and Im(z) ∼ N (Im(E{z}), σ 2 ). Then a/σ 2 follows a non-central χ-squared distribution with


F Soulez

F Courbin

M Unser

Back-propagating the light of field stars to probe telescope mirrors aberrations

We propose a wavefront-based method to estimate the PSF over the whole field of view. This method estimate the aberrations of all the mirrors of the telescope using only field stars. In this proof of concept paper, we described the method and present some qualitative results.

MOTIVATION

The Euclid 1 and WFIRST missions 2 will probe dark matter distribution using weak gravitational lensing. The precision needed on galaxy shape measurements required by weak lensing imposes stringent requirements on the PSF knowledge. The anisoplanatism of such wide-field telescope can not be neglected and the PSF have to be estimated for every position in the field of view. Field stars can give PSF measurements at random positions across the field of view. However, for weak lensing, PSF must be computed at each galaxy position (i.e. between field stars). The problem is thus twofold:

• PSF estimation at the position of each field star from its noisy observations,

• PSF interpolation at each galaxy position.

There are mainly two approaches to solve the PSF estimation problem: (i) image domain methods that parameterize PSF with pixels 3 and (ii) pupil domain methods. In the latter case, the PSF is described as a function of aberrations in the entrance pupil of the telescope. These pupil estimation methods relies on phase retrieval algorithms and most of it were conceived to estimate Hubble Space Telescope aberrations at the beginning of the 90's. [START_REF] Fienup | Hubble space telescope characterized by using phase-retrieval algorithms[END_REF][START_REF] Fienup | Phase retrieval for undersampled broadband images[END_REF][START_REF] Roddier | Combined approach to the hubble space telescope wave-front distortion analysis[END_REF][START_REF] Redding | Hubble space telescope prescription retrieval[END_REF][START_REF] Lyon | Hubble space telescope faint object camera calculated point-spread functions[END_REF][START_REF] Krist | Phase-retrieval analysis of pre-and post-repair hubble space telescope images[END_REF] The interpolation problem is then solved using a model of pupil aberration variation across the field of view.

In this paper, we propose to solve both problems jointly using a wavefront based method to estimate the PSF over the whole field of view. Indeed, the PSFs at every position of the field of view are fully characterized by aberrations of each optical surface of the telescope and can be computed using Fourier optics propagation.Although these aberrations can be calibrated on ground, it is probable that they will not remain stable enough after launch. One possible way to measure the wavefront on orbit would be to strongly defocus and refocus the telescope, an operation that is risky and therefore highly unlikely to be implemented by space agencies.

In this proof of concept paper, we propose a method to use scientific observations to estimate wavefront aberrations on the few optical surfaces of a space telescope. It uses each observed bright star as a source of a coherent plane wave to probe these aberrations as done for diffraction tomography. [START_REF] Kamilov | Learning approach to optical tomography[END_REF] This method can monitor the surface of every telescope mirrors bringing a new access to all its optical component status without any need to move optical elements. In addition, as it use stars present in the scientific channel, it does not require any additional calibration time. Finally, the knowledge of these optical surfaces will give the mean to estimate the PSF at all wavelengths and in each point of the field of view solving the problem of PSF interpolation on positions of lensed galaxies.

Determining mirrors aberrations using many images of stars is solved in an inverse problem framework. For each star, the forward model consists of free space propagation of a plane wave (whose angle is given by the star position) across the telescope optics ended by intensity recording in the detector plane. This model is non linear, however, as propagation between each mirror is a linear operation, modeling errors can be back-propagated and used to update the estimated aberrations of each mirrors. These back-propagated errors for many stars across the field of view are used by a continuous optimization algorithm (VMLMB, [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF] ) to probe precisely these aberrations. In this algorithm, the phase retrieval problem given the measured intensity is solved by the mean of an adapted proximity operator. [START_REF] Schutz | Painter: a spatiospectral image reconstruction algorithm for optical interferometry[END_REF]

IMAGE FORMATION MODEL

The forward model links the incoming wave w 1 (x, y) arriving on the telescope and the image recorded on the detector given the telescope parameters and its aberrations α. This model has two main parts:

• the propagation of the incoming wave w 1 (x, y) through the telescope to the detector plane,

• the measurement by the detector which records only the intensity (i.e. the squared modulus of the complex amplitude of the light in the detector plane) and is plagued by measurement noises such as both photon noise and read out noise.

The telescope model

The incoming wavefront emitted at wavelength λ by a single star at angular position (θ 1 , θ 2 ) relatively to the telescope optical axis, can be modeled as a plane wave. Its complex amplitude in the first mirror (M1) plane is given by

w 1 (x, y) = exp (ı (x sin(θ 1 )/λ + y sin(θ 2 )/λ)) . (1)
To define our forward model, this wave is adequately sampled on N pixels and we adopt a vector representation: w 1 = (w 1,1 , . . . , w 1,N ). The propagation of this wave through the telescope can be decomposed as a sequence of K similar operations, where K is the number of optical interfaces (mirrors, lenses and the detector). For each interface k, the incoming wave (i.e. the wave right before the interaction with the interface) can be itself modeled as a sequence of linear operations

w k (α) = H k-1 M k-1 A k-1 (α k ) w k-1 (α) (2)
where M k-1 and A k-1 (α k ) are two diagonal operators accounting for the effect of the (k -1) th interface and its aberration respectively and H k-1 is a propagation operator from the interface k -1 to the interface k. All these operators are square matrices in C N ×N . The aberration operator A k-1 (α k ) is a function of the unknown aberration parameters α k that will be estimated in our methods. Whereas the model described in Equation ( 2) is linear in w, it is highly non-linear in α.

Mirrors

A mirror M k modifies the incoming wave in two ways; (i) it cuts the light outside of its support p k (x, y) and (ii) it adds a space-varying phase term:

m k (x, y) = p k (x, y) exp (ı 2 φ k (x, y)) (3)
where φ k (x, y) is the sagitta of the mirror defined by

φ k (x, y) = x 2 + y 2 R k + R 2 k -(1 -2 k ) (x 2 + y 2 ) , (4)
where R k and k are the radius of curvature and the eccentricity respectively.

The mirror support p k is defined by

p k (x, y) = 1 , if (x, y) is inside the support 0 , otherwise. (5)
With an adequate sampling of m k (x, y), the discrete operator M k is diagonal and writes

M k = diag (m k ) . (6)

Aberrations

The aberrations are due to errors in the polishing of mirrors and optical misalignement. They are described by additional phase terms in the plane of each mirror. As the support of each mirror is usually a disk, the Zernike polynomials Z provide a suitable basis to express these aberrations. The aberrations of the k th mirror are described in the zernike basis Z k with the parameters α k . The aberration operator is

A k (α k ) = diag (exp (ı Z k α k )) . (7)

Propagation

The propagation operator H k from the interface k to the interface k + 1 is modeled using the paraxial approximation. Given the size of most telescopes, the Fresnel number is in general very high (∝ 10 6 ) and we define the propagation operator H k using the angular spectrum method.

Measurements

The detector measures only the intensity of the light wave. The forward model that links the complex amplitude in the detector plane w K ∈ C N to the measured image intensities d ∈ R N + is then

d n = |w K,n | 2 + e n , (8)
where e n is some measurement noise with spatially varying variance σ 2 n and |w K,n | 2 denotes the squared modulus of w K,n .

ALGORITHM

The goal of our algorithm is to estimate the vector of aberration parameters α using observations of S stars randomly distributed across the field of view. Assuming Gaussian measurement noise e, the estimated aberration parameters α + is the solution of the minimization problem:

α + = arg min α S s=1 N n=1 1 σ 2 n w s K,n (α) 2 -d s n 2 . ( 9
)
w s K,n is the complex amplitude at n th pixel of the detector of the light emitted by the s th star. It is modeled using the Equation 2 with w s 1 given by the Equation 1. This problem can be reformulated as a constrained problem:

α + arg min α S s=1 N n=1 1 σ 2 n |y s n | 2 -d s n 2 subject to y s n = w s K,n (α) , (10)
The Augmented Lagrangian formulation of this constrained problem is:

L(α, t, u) = S s=1 N n=1 1 σ 2 n |t s n | 2 -d s n 2 + ρ 2 S s=1 w s K (α) -t s -u s 2 2 , ( 11
)
where u s are the scaled Lagrange multipliers and ρ > 0 is the augmented penalty parameter.

Following Mourya et al., [START_REF] Mourya | Augmented lagrangian without alternating directions: Practical algorithms for inverse problems in imaging[END_REF] we solve this problem in a hierarchical way:

α = arg min α S s=1 N n=1 w s K,n (α) -t s n (α) -u s n 2 (12)
with t s n (α) = arg min The Equation 12 is solved using a continuous iterative optimization method (e.g. quasi-Newton method). The inner Equation 13is separable and consists on solving S × N small 1D problems that can be easily parallelized. At the end of each iteration k, we update the Lagrangian parameters:

t∈C 1 σ 2 n |t| 2 -d s n 2 + ρ 2 t -w s K,n (α) + u s n 2 (13)
u (k+1) = u (k) + w (k) -t (k) (14)

The phase retrieval problem

The inner minimization problem in Equation 13 is a phase retrieval problem. It is separable and can be defined as the solution of the proximity operator of the function f :

f (x) = 1 σ 2 (|x| 2 -d) 2 , ( 15
) prox 1/ρ f (t) = arg min x∈C 1 ρ f (x) + 1 2 |x -t| 2 . ( 16
)
This proximity operator has a closed form solution described in Schutz et al. [START_REF] Schutz | Painter: a spatiospectral image reconstruction algorithm for optical interferometry[END_REF]

The tomography problem

The outer minimization in Equation 12 is a tomography problem. It can be rewritten as:

α = arg min α S s=1 w K (α) -t -u 2 2 (17)
We solve this non linear problem using VMLMB, 11 a continuous optimization routine. The needed derivatives are computed using the recursive back-propagation algorithm described by Kamilov et al. [START_REF] Kamilov | Learning approach to optical tomography[END_REF]

RESULTS

We have tested our algorithm on simulations. We have simulated a Richtey-Chrétien telescope similar to the Hubble Space Telescope. Its characteristics are given by the Table 1. We introduce aberrations by drawing random coefficients of Zernike basis in Equation 7. We used 56 and 10 coefficients for the aberration of the first and the second mirrors respectively.

The dataset was generated with the telescope model described in Section 2. 50 stars distributed randomly across the field of view were generated. Their positions are shown on Figure 1. Their fluxes were adjusted such that 26000 photons on average were recorded per star ; that corresponds to a maximum intensity of 256 photons in the brightest pixel of the 1500 × 1500 pixels PSF. To simulate the detector, we add background noise of 5 e - and generate the data d s using the Poisson distribution P:

d s n = P( w s K,n 2 + 5) . (18)
The 300 × 300 pixels central part of the observation of the star indicated by an A on Figure 1 is shown on Figure 2.

We minimized Equation 12 using the unaberrated telescope as a starting point α = 0. To assess the performance of our method, we simulate observations of stars that were not in the data-set (denoted as B and C on Figure 1) using our aberrations estimate, the true aberrations and without any aberrations. These stars B and C and one of the star used in aberrations estimation (A) shown on Figure 3. On this figure, we can see that, beginning from an aberration free model, our algorithm successfully converges toward a PSF very similar to the ground truth.

CONCLUSION AND FUTURE WORKS

In this proof of concept paper, we show the validity of our approach to estimate PSFs at every positions across the field of view without any need to carry out in-flight telescope defocusing to measure directly the wavefront. It achieves to provide qualitatively good PSF even in noisy conditions and we are working on quantitative results in term of ellipticiy and size of the estimated PSF. In addition a lot of works has to be done to be able to process real data, particularly it has to handle undersampled and broadband PSFs.

Figure 1 .

 1 Figure 1. Distribution of the observed stars (blue stars) and test stars (red square) across the 0.6 °× 0.6 °field of view.

Figure 2 .

 2 Figure 2. 300 × 300 pixels central part of the recorded intensity for the star indicated by the letter A on Figure 1

Table 1 .

 1 Telescope simulation parameters

	Mirror 1
	Diameter (m)	2.4
	Curvature radius (m)	11.040
	Conic constant	-1.0022985
	Mirror 2
	Diameter (m)	0.281
	Curvature radius (m)	-1.358
	Conic constant	-1.496
	Distance (m)
	M1 to M2	4.9069
	M2 to detector	6.4062
	Detector
	Field of view (°)	0.6 × 0.6
	Pixel size (µm)	5
	Wavelength	500 nm

Acknowledgements

This work is supported by the Sinergia project "Euclid: precision cosmology in the dark sector" from the Swiss National Science Foundation


Member, IEEE Ihsen Hedhli

Senior Member, IEEE Gabriele Moser

Fellow, IEEE Sebastiano B Serpico

Fellow, IEEE Josiane Zerubia

Classification of Multisensor and Multiresolution Remote Sensing Images through Hierarchical Markov Random Fields

Keywords: Multisensor, multiresolution, multifrequency fusion, quad-tree, hierarchical Markov random field (MRF), marginal

whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

G IVEN the current amount and variety of data from high or very-high resolution (HR/VHR) satellite missions, a major challenge is to develop classifiers that can benefit from multiresolution, multisensor, and multifrequency input imagery [START_REF] Pohl | Remote sensing image fusion: an update in the context of digital Earth[END_REF]- [3]. The use of multiband and multiresolution data has been shown to favor accuracy and spatial precision of the classification maps [START_REF] Gómez-Chova | Multimodal classification of remote sensing images: a review and future directions[END_REF]. Furthermore, the use of multisensor data generally allows complementary properties, such as those of SAR and optical imagery, to be exploited [START_REF] Pohl | Remote sensing image fusion: an update in the context of digital Earth[END_REF]. Hence, there is a definite need for classification methods that can fuse images taken on the same area both by different sensors and at different resolutions. One way to address multisensor data classification is to compute joint class-conditional distributions given the marginal probability density function (PDF) of the data collected by each sensor. Meta-Gaussian [START_REF] Storvik | On the combination of multisensor data using meta-Gaussian distributions[END_REF], copula [START_REF] Voisin | Supervised classification of multi-sensor and multi-resolution remote sensing images with a hierarchical copula-based approach[END_REF], or non-parametric models [START_REF] Landgrebe | Signal theory methods in multispectral remote sensing[END_REF] can be used for this purpose. However, finding such joint statistical models is generally complex, time-demanding, and possibly prone to over-fitting.

In this paper, we build on the theory of hierarchical Markov random fields (MRFs) [START_REF] Li | Markov random field modeling in image analysis[END_REF] to address the challenging problem of the supervised classification of data that are both multiresolution and multisensor (and possibly multifrequency). The rationale is to benefit from the fusion capabilities of I. Hedhli and J. Zerubia are with UCA (Université Côte d'Azur), Inria, France (e-mail: ihsen.hedhli@inria.fr, josiane.zerubia@inria.fr). I. Hedhli, G. Moser, and S. B. Serpico are with the University of Genoa, Italy (e-mail: gabriele.moser@unige.it). The authors would like to thank the French Space Agency (CNES), the Italian Space Agency (ASI). and the Canadian Space Agency (CSA) for providing Pléiades, COSMO-SkyMed, and RADARSAT-2 images, respectively. Part of the data was procured within the SOAR-ASI 5245 project framed within the joint ASI-CSA announcement of opportunity. The support is gratefully acknowledged.

hierarchical MRFs and avoid the computation of joint statistics of multisensor data. Inference is accomplished by extending the approach that was previously developed in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF] for multitemporal fusion and is based on the marginal posterior mode (MPM) criterion. PDF estimation is addressed using mixtures of Gaussian and generalized Gamma distributions.

Two methods are proposed in this framework. In the case of the first one, the image acquired by each sensor is associated with a separate quad-tree according to its resolution. In regard to the second method, the focus is on the specific case of the fusion of multifrequency X-band COSMO-SkyMed (CSK) and C-band RADARSAT-2 (RS2) SAR data together with optical visible and near-infrared (VNIR) Pléiades data. Distinct quadtrees are used again, but here, optical and SAR data are both included in each quad-tree to benefit from the finest resolution available from each sensor. While the first technique addresses the fusion of data from generally arbitrary SAR and optical sensors, the second method focuses on the synergy among three specific and topical satellite data sources. The main novel contribution of the proposed framework is the formalization of the joint problem of multisensor, multiresolution, and possibly multifrequency fusion in terms of a multiple quadtree topology and of hierarchical MRFs.

The letter is organized as follows. Section II-A describes the proposed framework and the quad-tree structures of the two methods. Sections II-B and II-C discuss the MPM inference and PDF estimation issues. Experimental results are presented in Section III and conclusions are drawn in Section IV.

II. PROPOSED MULTISENSOR HIERARCHICAL APPROACH A. Overview of the proposed methodology

To address multisensor and possibly multifrequency fusion, the proposed framework extends and adapts the approach developed in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF] for multitemporal fusion and based on multiple quad-trees in cascade and on hierarchical MRFs (Fig. 1). This hierarchical formulation also allows multiresolution data to be naturally incorporated. In particular, a pair of quadtrees in cascade is constructed using the input multisource imagery. Given a training set for M thematic classes, inference can be efficiently formulated on this topology using MPM (see Section II-B), while the statistics of the data from each sensor are modeled using finite mixtures and the stochastic expectation-maximization (SEM) algorithm [START_REF] Celeux | Stochastic versions of the EM algorithm: an experimental study in the mixture case[END_REF] (see Section II-C). In this framework, two methods, which differ in how data are incorporated in the quad-trees, are defined. 1) Quad-trees of the first proposed method: Let an optical and a SAR image, acquired at the same time over the same area, be given. In the first method, these images are included in the finest-scale levels (the leaves) of two quad-trees. Coarserscale levels of each quad-tree are filled in using wavelet transforms of the finest-scale level [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The roots of the quad-trees are assumed to have the same spatial resolution. Hence, there generally are different numbers of levels in the quad-trees, from the roots at a common resolution to the input resolutions of the optical and SAR data. The quad-tree topology also implies that the resolutions of the two input images are in a power-of-2 relation. With current HR satellite sensors, this condition is easily met up to minor resampling.

2) Quad-trees of the second proposed method: The second method is aimed at classifying multisensor, multifrequency (Xband radar, C-band radar, VNIR optical), and multiresolution data from CSK, RS2, and Pléiades. These prominent sensors support multiple resolutions, up to 0.5 m for Pléiades and approximately 1 m for CSK and 1-2 m for RS2. Given input CSK, RS2, and Pléiades images acquired on the same area at the same time, case-specific quad-trees are constructed to classify at the finest of the input resolutions. Each SAR image is inserted in a separate quad-tree according to its own resolution. These resolutions are expected to be coarser than the finest resolution achieved using Pléiades. Therefore, the input Pléiades image is embedded in the finest-resolution levels of both quad-trees. Empty levels are filled in using wavelet transforms of the Pléiades data. The comment on the power-of-2 relation among the input resolutions holds again.

B. Multisensor hierarchical formulation of MPM

Owing to the quad-tree structure, a hierarchical MRF model, which consists of a series of random fields linked, at various scales, by suitable transition relations, can be defined (Fig. 1). The formulation of MPM proposed in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF] for multitemporal classification is extended here to multisensor fusion. MPM is used instead of the maximum a-posteriori (MAP) rule because MAP formulations on quad-trees, although they are feasible [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF], [START_REF] Dawid | Applications of a general propagation algorithm for probabilistic expert systems[END_REF], are affected by underflow and penalize errors regardless of the scale at which they occur -an undesired behavior for multiscale classification [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF], [START_REF] Dawid | Applications of a general propagation algorithm for probabilistic expert systems[END_REF], [START_REF] Pérez | Noniterative manipulation of discrete energy-based models for image analysis[END_REF].

Collecting all the (optical, SAR, and wavelet) features of both quad-trees in the vector y, the posterior marginal p(x s |y) of the label x s of each site s in the second quadtree is expressed, under suitable conditional independence assumptions, as a function of the posterior marginals p(x s -|y) and p(x s = |y) of the parent node s -in the same quad-tree and the parent node s = in the first quad-tree, i.e. [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF]:

p(x s |y) = x s -,x s = p(x s , x s -, x s = |y d(s) )p(x s -|y)p(x s = |y) xs p(x s , x s -, x s = |y d(s) ) , (1)
where it is proven in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF] that, under appropriate assumptions:

p(x s , x s -, x s = |y d(s) )p(x s ) = (2) = p(x s |x s -, x s = )p(x s -|x s = )p(x s = )p(x s |y d(s) ),
and where y d(s) collects the features of all descendants of s, p(x s ) is the prior probability, p(x s |x s -, x s = ) is the parentchild transition probability, p(x s -|x s = ) is the transition probability between sites at the same scale, and p(x s |y d(s) ) is the partial posterior marginal. While these quantities are made available, ( 1) and ( 2) recursively calculate p(x s |y) on each site s of the second quad-tree. This is obtained using one bottomup and two top-down recursive passes, which are described in the next subsections. Note that, by ( 1), the quad-trees are ordered. In the first proposed method, SAR and optical data are related to the first and second quad-trees, respectively. In the second method, two options are possible, i.e., CSK in the first quad-tree and RS2 in the second one or vice versa.

1) Initialization on the first quad-tree: First, classification is performed only on the first quad-tree using a traditional MPM on a single quad-tree [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF]. As in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF] (to which we refer for details), the prior distribution at the root is modeled using a Potts MRF to favor spatial consistency. After the initialization stage, p(x f |y) is known for each site f of the first quad-tree and p(x f |y d(f ) ) is derived as a by-product.

2) First top-down pass on the second quad-tree: The first top-down pass proceeds from the root to the leaves of the second quad-tree to recursively compute p(x s ). As in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF], this prior is initialized in each site s of the root of the second quadtree as p(x s ) = p(x f |y d(f ) ), where f is the site of the root of the first quad-tree with the same location as s. Then, the recursion proceeds along the other levels down to the leaves:

p(x s ) = x s - p(x s |x s -)p(x s -). (3)
This favors the same parent-child labeling and characterizes the statistical interactions between consecutive levels. We model the parent-child transition probability p(x s |x s -) through the parametric form in [START_REF] Bouman | A multiscale image model for Bayesian image segmentation[END_REF]. As a result of this top-down pass, p(x s ) is derived on each site s of the second quad-tree.

3) Bottom-up pass on the second quad-tree: A bottomup recursion is performed to compute p(x s , x s -, x s = |y d(s) ) starting from the leaves of the second quad-tree and using (2) until the root is reached. In addition to p(x s ), three sets of probabilities are required: (i) the transition probabilities at the same scale p(x s -|x s = ); (ii) the parent-child transition probabilities p(x s |x s -, x s = ); and (iii) the partial posterior marginals p(x s |y d(s) ). Details of the modeling of (i) and (ii), which is performed parametrically, can be found in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF]. Concerning (iii), it has been proved that [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF]:

p(x s |y d(s) ) ∝ p(y s |x s )p(x s ) u∈s + xu p(x u |y d(u) )p(x u |x s ) p(x u ) (4)
where s + and y s are the set of four children and the feature vector of site s, respectively. First, p(x s |y d(s) ) is initialized on the leaves of the second quad-tree with p(x s |y d(s) ) = p(x s |y s ) ∝ p(y s |x s )p(x s ). Then, p(x s |y d(s) ) is recursively estimated with (4) up to the root. This recursion involves the pixelwise class-conditional PDF p(y s |x s ) (see Section II-C).

After the bottom-up pass, p(x s , x s -, x s = |y d(s) ) is computed on each site s of the second quad-tree. 4) Second top-down pass on the second quad-tree and output labeling: First, the posterior marginal is initialized at the root of the second quad-tree as p(x s |y) = p(x s |y d(s) ). Then, given the quantities computed by the previous stages, p(x s |y) is derived on all sites s of all other levels of the second quad-tree in a top-down pass using the recursion (1). Then, in principle, the class label x s that maximizes p(x s |y) over the set of M classes could be assigned to s directly. This approach is often avoided in the literature of hierarchical MRFs because of its computational burden and of possible blocky artifacts [START_REF] Voisin | Supervised classification of multi-sensor and multi-resolution remote sensing images with a hierarchical copula-based approach[END_REF], [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF]. As an alternate approach, the casespecific formulation of the modified Metropolis dynamics algorithm [START_REF] Berthod | Bayesian image classification using Markov random fields[END_REF] that was combined in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF] with MPM is extended here as well. The labeling result on the leaves of the second quad-tree represents the output classification map.

C. PDF estimation through finite mixtures

For each class m ∈ {1, 2, . . . , M } and each level n between the leaves and the root of each quad-tree q ∈ {1, 2}, we use a finite-mixture model for the pixelwise PDF p(y s |x s = m). This choice is explained by the flexibility of finite mixtures in characterizing heterogeneous statistics such as those of HR imagery [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF]. If the data at level n of quad-tree q are optical, then a multivariate Gaussian mixture is used [START_REF] Landgrebe | Signal theory methods in multispectral remote sensing[END_REF]. The same model is used for wavelet transforms of optical data, because of the linearity of the wavelets. If the data at level n of quad-tree q are SAR, a mixture of generalized Gamma distributions is used. The generalized Gamma model has been found accurate for HR SAR [START_REF] Li | On the empirical-statistical modeling of SAR images with generalized Gamma distribution[END_REF], [START_REF] Krylov | On the method of logarithmic cumulants for parametric probability density function estimation[END_REF]; here, it is also extended as an empirical model for wavelet transforms of SAR data.

The mixture parameters are estimated using SEM, which is an iterative stochastic algorithm developed for problems characterized by data incompleteness. SEM generates a random sequence of estimates, which does not converge almost surely but, under suitable assumptions, forms an ergodic homogeneous Markov chain converging to a unique stationary distribution, which is expected to concentrate around the global loglikelihood maxima [START_REF] Celeux | Stochastic versions of the EM algorithm: an experimental study in the mixture case[END_REF]. Hence, initialization is not critical and is usually performed randomly with uniform distribution [START_REF] Celeux | Stochastic versions of the EM algorithm: an experimental study in the mixture case[END_REF], [START_REF] Moser | Dictionary-based stochastic expectation-maximization for sar amplitude probability density function estimation[END_REF]. The number of mixture components is also automatically estimated using the specific SEM formulation in [START_REF] Krylov | Enhanced dictionary-based SAR amplitude distribution estimation and its validation with very high-resolution data[END_REF]. For each scale n of each quad-tree q, SEM is separately applied to the training samples of each class m. For generalized Gamma mixtures, it is also combined with the method of log-cumulants (details can be found in [START_REF] Moser | Dictionary-based stochastic expectation-maximization for sar amplitude probability density function estimation[END_REF], [START_REF] Krylov | Enhanced dictionary-based SAR amplitude distribution estimation and its validation with very high-resolution data[END_REF]).

III. EXPERIMENTAL RESULTS

A. Results of the first proposed method

An HR dataset collected in 2010 over Port-au-Prince, Haiti, and composed of an HH-polarized single-look CSK stripmap image with 2.5-m pixel spacing (325 × 400 pixels; Fig. 2(a)) and of a GeoEye-1 RGB image resampled at 2.5-m resolution (Fig. 2(b)), was used to experimentally validate the first proposed method. 2D Daubechies wavelets of order 10 were used. Five land cover classes were considered (see Table I and Fig. 2(c)). We compared the result of the proposed method (Fig. 2(i)) with those of: (i) the single-sensor multiscale classification of the image of each sensor, obtained by applying the original MPM to a single quad-tree [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF] with Gaussian or generalized Gamma mixtures for the class-conditional PDFs of optical and SAR data, respectively (Figs. 2(d) and(e)); (ii) the multisensor multiresolution method in [START_REF] Voisin | Supervised classification of multi-sensor and multi-resolution remote sensing images with a hierarchical copula-based approach[END_REF], in which a hierarchical MRF on a single quad-tree is used and multisensor fusion is addressed through copulas (Fig. 2(f)); (iii) the multisensor single-resolution approach in [START_REF] Storvik | On the combination of multisensor data using meta-Gaussian distributions[END_REF], in which the joint class-conditional PDFs of multisensor data are modeled with meta-Gaussian distributions and combined with the maximum likelihood rule (Fig. 2(g)); and (iv) a support vector machine (SVM) applied to the stacked optical-SAR data as a nonparametric pixelwise benchmark (Fig. 2(h)). The Gaussian radial basis function kernel was used and the SVM parameters were optimized using the method in [START_REF] Moser | Combining support vector machines and Markov random fields in an integrated framework for contextual image classification[END_REF].

A visual analysis of the classification maps suggest that the first proposed method led to rather accurate results, especially as compared to the separate hierarchical classifications of the images of the two sensors. The results obtained using only SAR data accurately detected roads and poorly discriminated most other classes, while the map generated using only optical data better discriminated classes that were spatially homogeneous. The proposed method took benefit from both sources and generated a result that visually discriminated quite well most classes. Compared to the multisensor single-scale method in [START_REF] Storvik | On the combination of multisensor data using meta-Gaussian distributions[END_REF], the proposed algorithm provided a spatially more regular map thanks to MRF modeling and wavelets.

Quantitative analysis on the test set confirmed that the proposed method provided quite high overall accuracy (OA) and Cohen's κ, and improved as compared to the method in [START_REF] Storvik | On the combination of multisensor data using meta-Gaussian distributions[END_REF] (Table I). Poor discrimination was obtained by all techniques for "containers," because of its large overlapping with "urban" in the multisensor or multispectral feature space. To improve the accuracy for this class, texture features could be integrated as an extension of the proposed method. The multisensor multiscale method in [START_REF] Voisin | Supervised classification of multi-sensor and multi-resolution remote sensing images with a hierarchical copula-based approach[END_REF] also provided high accuracies, similar to those of the proposed algorithm, with slightly lower OA and higher κ. However, as compared to this technique, the proposed method improved in terms of spatial detail (Figs. 2(f) and (i)) and computation time (Table I). The approach to multisensor fusion based on distinct quad-trees related to individual sensors avoids the difficult and possibly time-expensive task of estimating the joint optical-SAR PDFs, whereas copula and meta-Gaussian models are used for this task in [START_REF] Voisin | Supervised classification of multi-sensor and multi-resolution remote sensing images with a hierarchical copula-based approach[END_REF] and [START_REF] Storvik | On the combination of multisensor data using meta-Gaussian distributions[END_REF], respectively. The SVM provided lower accuracy than the other considered methods, which was expected because of its noncontextual formulation. The values of McNemar's Z statistics of the maps of the benchmark methods with respect to the map of the proposed method are all much larger than the 5% significance threshold (|Z| > 1.96), which confirms that the differences among these results are significant.

B. Results of the second proposed method

A challenging VHR dataset, collected again over Port-au-Prince and composed of a Pléiades pansharpened image at 0.5m resolution (1000 × 1000 pixels; Fig. 3 II). Given the 0.5-m resolution of the Pléiades data, we slightly resampled the RS2 image to obtain 2-m pixel spacing and fit with the power-of-2 constraint implied by the quad-tree. Downsampling from 1.56 to 2 m is expected to have a minor impact on the classification map, because the resampling ratio is quite close to 1. Preliminary experiments (not shown for brevity) indicated that the order in which CSK and RS2 data were inserted in the quad-trees did not significantly affect the results.

The second proposed method obtained rather high OA and κ on the test set (Table II) and remarkable spatial regularity (Fig. 3(h)). Again, poor discrimination was obtained for "containers." To investigate the capability of the method to exploit the synergy among X-band, C-band, and VNIR data, comparisons with the results obtained by using the Pléiades image only by itself or together with either SAR image (Table II and Figs. 3(d), (e), and (f)) were conducted. In all these comparisons, the same quad-tree, MPM, and finitemixture formulations as in the proposed method were used. The results obtained using only SAR data poorly discriminated the classes and are not reported for brevity. A comparison with a nonparametric pixelwise benchmark was performed by resampling the CSK and RS2 images on the 0.5-m grid, stacking them together, and applying an SVM with the same model selection strategy as in Section III-A (Fig. 3(g)). Quantitative analysis on the test set confirmed that the joint use of all three sensors led to substantially higher accuracies as compared to the use of a subset of the input data sources. The use of only the Pléiades image allowed "water," "vegetation," and "bare soil" to be effectively discriminated but led to poor detection of "urban." The joint use of this VNIR image and of X-band and C-band SAR through the second proposed method allowed improving the discrimination of "urban" of nearly 30% and outperformed the results obtained using the Pléiades data together with only one of the two SAR images. These results suggest the potential of the second proposed method in exploiting the complementarity of multisensor, multifrequency, and multiresolution data from current VHR satellite sensors. In this case as well, the noncontextual SVM obtained lower test-set accuracy and a more noisy map than the proposed method. McNemar's test indicates that the differences between the results of the proposed method and of each benchmark are again significant. A shortcoming of the proposed technique was the lower accuracy for "water" than when only Pléiades or Pléiades and CSK data were used. Indeed, "water" appears highly textured in the RS2 image (Fig. 3(c)), whereas the proposed technique uses no texture features. On one hand, the accuracy of the proposed method for "water" is still around 91%. On the other hand, incorporating texture features could be a feasible extension of the second proposed method as well.

IV. CONCLUSION

A hierarchical MRF approach to multisensor and multiresolution optical-SAR image classification has been formulated, extending inference procedures originally introduced for multitemporal fusion and incorporating suitable finite mixture models. Two classification methods have been developed in this framework. Experimental results on challenging HR/VHR data sets have suggested the effectiveness of this approach to data fusion and the improvement of the two proposed methods as compared to single-sensor, single-resolution, and/or previous techniques in terms of accuracy, spatial precision, or computational burden. A further advantage of the proposed approach to fusion is its flexibility. While multitemporal fusion was already formalized in this framework in [START_REF] Hedhli | A new cascade model for the hierarchical joint classification of multitemporal and multiresolution remote sensing data[END_REF], the multiple quad-tree structure could be extended to also incorporate data from further sensors or feature extraction stages (e.g., textures or morphological profiles). This extension will be pursued in the near future. Moreover, alternate wavelet de- I.

compositions [START_REF] Mallat | A wavelet tour of signal processing[END_REF] will also be investigated, including wavelet packets, with which the number of transformed samples in the quad-tree is increased, and undecimated stationary wavelets, which preserve image size along different levels and could be incorporated into a multisource MRF energy as in [START_REF] Solberg | A Markov random field model for classification of multisource satellite imagery[END_REF].

Fig. 1 .

 1 Fig. 1. Multisensor hierarchical structure.

  (a)), an HH-polarized CSK spotlight image at 1-m resolution (Fig. 3(b)), and an HH-polarized RS2 ultrafine image with 1.56-m pixel spacing (Fig. 3(c)), was used for experiments. The same five classes of the previous section were considered (Table

Fig. 2 .

 2 Fig. 2. First proposed method: (a) CSK ( c ASI 2010) and (b) GeoEye-1 ( c GeoEye 2010; R-band display) images; (c) ground truth (training and test areas); results of the hierarchical MRF-based classifications of (d) the GeoEye-1 and (e) the CSK image; and (f) results of the techniques in[START_REF] Voisin | Supervised classification of multi-sensor and multi-resolution remote sensing images with a hierarchical copula-based approach[END_REF] and (g) in[START_REF] Storvik | On the combination of multisensor data using meta-Gaussian distributions[END_REF], (h) SVM, and (i) the proposed method. Color legend: see TableI.

TABLE I FIRST

 I PROPOSED METHOD, CSK-GEOEYE-1 DATASET: TRAINING AND TEST SAMPLE SIZES; CLASSIFICATION ACCURACIES AND κ; MCNEMAR'S Z WITH RESPECT TO THE PROPOSED METHOD; AND COMPUTATION TIMES (INTEL I7 QUAD-CORE, 2.40 GHZ, 8 GB RAM, 64 BIT LINUX SYSTEM).

		water	urban	vegetation	bare soil	containers	OA	κ	Z	computation time [s]
	Proposed method	100%	78.12%	89.46%	98.78%	47.12%	82.69%	0.76	-		254
	Method in [4]	99.95%	97.32%	90.81%	96.22%	37.25%	79.44%	0.68	1843.7	298
	Method in [5]	100%	75.24%	87.16%	98.89%	49.31%	82.12%	0.79	89.6		668
	SVM	96.52%	70.31%	44.28%	78.15%	43.20%	66.49%	0.57	2298.3	372
	Training samples	50043	79492	51312	11667	38936
	Test Samples	49057	75508	50688	10333	31064
					TABLE II
	SECOND PROPOSED METHOD, CSK-RS2-PL ÉIADES DATASET: TRAINING AND TEST SAMPLE SIZES; CLASSIFICATION ACCURACIES AND κ; AND
			MCNEMAR'S Z WITH RESPECT TO THE PROPOSED METHOD.
				water	urban	vegetation	bare soil	containers	OA	κ	Z
	Pléiades only		100%	61.66%	81.69%	82.82%	56.72%	76.57%	0.63	1361.7
	Pléiades + CSK		100%	44.32%	83.54%	74.75%	49.12%	70.34%	0.58	2102.7
	Pléiades + RS2		92.56%	44.85%	79.85%	78.62%	42.15%	67.60%	0.57	2163.7
	Proposed method (Pléiades + RS2 + CSK)	90.79%	91.45%	82.59%	81.02%	54.85%	80.14%	0.78	-
		SVM		92.86%	92.56%	61.86%	69.91%	56.64%	74.76%	0.61	2302.9
	Training samples		12943	15492	10312	8667	6016
	Test samples		12057	14508	9688	6333	5984