Fast and Accurate Multiplicative Decomposition for Fringe Removal in Interferometric Images

Daniel-Chen Soncco, Clara Barbanson, Mila Nikolova, Senior Member, IEEE, Andrés Almansa, Member, IEEE, Y. Ferrec

Abstract-Airborne Hyperspectral images can be efficiently obtained with imaging static Fourier transform spectrometers. However, to be effective on any location, this technology requires to know the relief of the scene. This is not a straightforward process, as the horizontal interference fringes on the images, which are necessary for spectrum construction, prevent efficient stereoscopic processing. We present a novel variational model for multiplicative image decomposition to separate the fringes from the panchromatic image of the scene. This multiplicative model is much more physically accurate than previous additive decomposition models inspired by cartoon-texture decomposition. It combines fully smoothed total variation operators and 1D Fourier transform. Smoothed total variation is adopted to avoid staircasing artifacts caused by traditional total variation regularisation. The use of a 1D Fourier transform is suggested by the geometry of the fringes, in order to circumvent the lack of horizontal periodicity in the interferometric pattern. We also present an optimization algorithm. Finally, a second algorithm is introduced, whose convergence is not mathematically guaranteed. However it systematically approaches the solution of the first one in much less computation time. Our experimental evaluation on real and simulated images shows that the proposed model separates fringes from the panchromatic image very accurately and that this accuracy significantly improves subpixel stereo matching results.

Index Terms-computational hyperspectral imaging, interferometric imaging, inverse problems, multiplicative image decomposition, image restoration, nonconvex variational models, bi-convex optimisation, proximal algorithms, multi-view stereo

I. INTRODUCTION

H YPERSPECTRAL imaging is a much sought-after tech- nology, as it provides both detailed spatial and spectral information from the object of interest, with a spectral range that may go, depending on the instrument, from ultraviolet to thermal infrared. Such technology is widely used in Earth sciences; for instance, geology employs it on airborne or spaceborne vehicles for remote sensing [START_REF] Van Der Meer | Multi-and hyperspectral geologic remote sensing: A review[END_REF], where infrared wavelengths are especially useful for prospection [START_REF] Salles | Hyperspectral remote sensing applied to uranium exploration: A case study at the Mary Kathleen metamorphichydrothermal U-REE deposit, NW, Queensland, Australia[END_REF], [START_REF] Reath | Exploration of geothermal systems using hyperspectral thermal infrared remote sensing[END_REF], soil composition determination [START_REF] Eisele | Advantages using the thermal infrared (TIR) to detect and quantify semi-arid soil properties[END_REF], [START_REF] Tavin | Comparison of Metrics for the Classification of Soils Under Variable Geometrical Conditions Using Hyperspectral Data[END_REF] and even gas plumes analysis in vulcanology [START_REF] Spinetti | Carbon dioxide of Pu'u'O'o volcanic plume at Kilauea retrieved by AVIRIS hyperspectral data[END_REF]. Other gases such as water vapor for atmospheric research can be identified [START_REF] Burke | Satellite hyperspectral IR sensing: investigating water vapor variability using AIRS data[END_REF]. Visible and near-infrared are used in hydrology to assess water contamination [START_REF] Marino | Imaging spectrometry and GIS techniques application for the study of seawater quality and pollution hazard[END_REF] and, in relation to agriculture, to identify different D. Soncco and M. Nikolova are with CMLA, CNRS and Ecole Normale Supérieure de Cachan, Université Paris-Saclay, 94235 Cachan, France. (email chen.soncco@gmail.com, nikolova@cmla.ens-cachan.fr ) C. Barbanson crop fields and monitor water use [START_REF] Fais | Airborne hyperspectral scanner and laser altimeter data application to water reservoirs and water use calculation: first results on two Italian case study[END_REF]. Biology is also wellknown for using infrared bands for vegetation classification [START_REF] Dalponte | Fusion of Hyperspectral and LIDAR Remote Sensing Data for Classification of Complex Forest Areas[END_REF]. Outside Earth sciences, applications also include urban planning [START_REF] Clark | Mapping of land cover in northern California with simulated hyperspectral satellite imagery[END_REF] and military intelligence [START_REF] Briottet | Military applications of hyperspectral imagery[END_REF]. Finally, elevation data are sometimes coupled with hyperspectral imagery to improve detection and analysis [START_REF] Gleeson | Characterization of a sulfur-rich Arctic spring site and field analog to Europa using hyperspectral data[END_REF], [START_REF] Feng | The topographic normalization of hyperspectral data: implications for the selection of spectral end members and lithologic mapping[END_REF] Imaging static Fourier transform spectrometers (ISFTS) are a family of imagers which deliver hyperspectral data. They offer a series of interesting advantages, including large photon collection capability, that allow for a small ground sampling distance, and sturdy optical systems. These instruments can be operated from ground [START_REF] Yang | Research on ground-based lwir hyperspectral imaging remote gas detection[END_REF], [START_REF] Fossi | Miniature and cooled hyperspectral camera for outdoor surveillance applications in the mid-infrared[END_REF] to space [START_REF] Wright | Tircis: thermal infrared compact imaging spectrometer for small satellite applications[END_REF], but they are more often airborne systems, like Lasis, developed by the Xi'an Institute of Optics and Precision Mechanics [START_REF] Wang | Research on lasis interferogram processing[END_REF], that operates in the visible domain, Sieleters, developed by Onera [START_REF] Coudrain | SIELETERS, an airborne infrared dual-band spectro-imaging system for measurement of scene spectral signatures[END_REF], in the midwave infrared (MWIR) and the longwave infrared (LWIR) domain or THI, developed by the University of Hawai'i [START_REF] Lucey | A LWIR hyperspectral imager using a sagnac interferometer and cooled HgCdTe detector array[END_REF], in LWIR.

The optical design of ISFTS is based on the association of an imaging system with a lateral shearing interferometer (apart from a few instruments employing a wedge interferometer). The incident beam emitted from the observed scene therefore proceeds through the interferometer, where it is split into two paths. Two identical images are produced, but they later interfere at a varying optical path difference (OPD) thanks to the interferometer. The resulting image is that of the overflown scene on which interference fringes are superimposed, as shown in Fig. 1.

On the figure, the OPD is different on each row of the image: the fringes are oriented perpendicularly to the linear flight path. Furthermore, the established OPD values do not change over time, which is why such an interferometer is qualified as static. The fringes are thus still. For this reason, spectrum acquisition is not snapshot, but temporal: it is the motion of the aircraft which ensures that any element in the observed scene is seen through different interference states (at different vertical positions in the image). A dense sequence of images, where each one is translated from the former one by about 1 pixel, is required for this purpose. Then, matching each scene element through the sequence of images constructs an interferogram, where each value corresponds to the quantity of detected light for each optical path difference. Computing the Fourier transform of the interferogram gives access to the spectrum for the corresponding element. This process is summarised in Fig. 2.

However, in presence of non-flat scenes (containing for example buildings or natural slopes), artifacts may appear as a consequence of changing the point of view during acquisition.

For instance, due to parallax, elevated objects appear to move at a faster speed than the elements of lesser height. Without a precise estimate of the scene relief, the position of an elevated object across the sequence of images is miscalculated and matching errors occur. The resulting interferograms may then be erroneous, as shown in Fig. 3.

As a result, an accurate estimation of ground elevation is needed in order to identify and correct flawed interferograms. A Digital Elevation Model (DEM) can be achieved by stereoscopy techniques. However, the fringes are a hindrance to stereoscopic methods. The fact that they remain at the same position on all images interferes with relief retrieval algorithms, which exploit the changes between images taken from different positions. Moreover, the fringes induce variations in radiometry that are uncorrelated to relief. As a consequence, before proceeding with relief estimation, the fringes need to be removed from the images. This problem is the topic of this article: how can one remove the interference fringes while best preserving the details of the underlying panchromatic image?

A similar problem has been studied in [START_REF] Wen | Optics and Lasers Fig. 16. On the top line, the two images used for stereoscopy, as they were measured. The area designated by[END_REF]: in this paper, based on images from Lasis, the authors also aim at separating fringes from background, but with the view of image compression. To do this, they use an additive model, with a further one-dimensional total variation (TV) regularization to suppress the residual fringes from the background image. However, the images presented in [START_REF] Wen | Optics and Lasers Fig. 16. On the top line, the two images used for stereoscopy, as they were measured. The area designated by[END_REF] to illustrate their algorithm have a very smooth texture. In this particular setting an additive image model is sufficient, but not in the more general setting we are interested in, which includes irregular textures and contrasted edges. The problem of fringe separation has also been studied in a more general context by [START_REF] Maurel | Locally parallel texture modeling[END_REF], there again with an additive model, and using an adaptive norm to deal with spatial variations of the fringes orientation or frequency. Applications are various, for instance fingerprint analysis or 3D-imaging by strucured-light projection [START_REF] Zhou | Morphological operation-based bi-dimensional empirical mode decomposition for automatic background removal of fringe patterns[END_REF].

The main contribution of our paper consists in the use of a multiplicative model for image decomposition. This model is a better approximation of images obtained by an ISFTS than an additive model. It provides outstanding results for separating the fringes from the image of the scene. The paper is arranged in the following manner. Section II presents the model for interferometric image formation and related work on the texture removal problem. Our proposed solution is presented in Section III. It is declined in two different algorithms using a multiplicative image model. The advantages and drawbacks of the different solutions are also discussed. Section IV presents experimental results and gives a comparison of the various solutions in terms of stereoscopic performance. Finally, conclusion and future work are exposed in Section V.

II. PRELIMINARIES

A. Observation model

The primary objective of our hyperspectral interferometric system is to measure the spectral radiance L(x, y, σ) coming from a point (x, y) in image or ground coordinates in a domain D and at a particular wavenumber σ. Color image sensors measure L(x, y, σ) for three wavenumbers σ corresponding to show the interference fringes overlaid on the scene. In red, an area with highly contrasted fringes. In blue, an area with blurred fringes. In green, an area with slightly contrasted fringes. The area with high contrast fringes is limited to the area around the zero optical path difference, due to the broadband observed spectra. Nonetheless, fringes are actually present over the whole image. Image (b) shows the Fourier spectrum of picture (a). The frequencies of the fringes are clearly seen as two very bright vertical lines. These frequencies constitute an area we call Ω 0 . the red, green and blue wavelengths by means of optical filters that only let the desired wavelength intervals through. In hyperspectral imaging, where we do not only want to distinguish three, but hundreds of different wavelengths, filtering methods can be substituted by interferometers.

The geometry of the lateral shearing interferometer [START_REF] Ferrec | Spectro-imagerie aéroportée par transformation de fourier avec un interféromètre statique à décalage latéral : Réalisation et mise en oeuvre[END_REF] shows that the acquired image w(x, y) is related to the apparent spectral radiance 1 L by:

w(x, y) = L(x, y, σ) 1 + F (x, y, σ) dσ ( 1
)
where σ is the wavenumber and F introduces the fringe oscillations in the acquired image w(x, y). More precisely,

F (x, y, σ) = C(σ) cos 2πδ(x, y)σ + φ(σ) (2)
where C is a contrast term, δ(x, y) is the optical path difference (OPD) and φ accounts for phase perturbations. In ISTFS, φ has a little influence and the OPD is nearly planar, δ(x, y) = ax + by + ε(x, y)

1 By apparent we mean that the sensitivity of the detector is included in this term.

where the variations of ε are known to be very small. Without loss of generality, we consider that a b. Consequently, the fringes are nearly horizontal straight lines as seen in Fig. 1. From (2), the wavenumber σ creates fringes at the spatial frequencies -σ(a, b) ∪ σ(a, b) (assuming a flat scene). The sensitivity of the imaging device is given by an interval [σ min , σ max ]. Then the fringes are contained in the spatial frequency domain Ω 0 :

Ω 0 = Ω 0x × Ω 0y = -σ(a, b) ∪ σ(a, b) : σ ∈ [σ min , σ max ]
The goal is to extract the fringe pattern introduced by F . Based on (1), two tractable simplifications can be derived. where u is the sought-after panchromatic image and v is the fringe pattern image.

If we now assume that the scene u only contains smooth objects, then, in the spatial frequency domain, the Fourier transform v of the fringe pattern will be limited to a small area Ω around Ω 0 .

w(x, y) = I(x, y)B tot 1 + 1 B tot B(σ)F (x, y, σ)dσ (6)
We then get the following image model :

w(x, y) = u(x, y) 1 + v(x, y) (7)
where the panchromatic image u and the fringe image v are given by :

u(x, y) = I(x, y)B tot (8) v(x, y) = 1 B tot B(σ)F (x, y, σ)dσ (9)
Since B does not depend on (x, y), and F only has frequencies contained in Ω, we can deduce that the spectral support of v should be exactly limited to Ω 0 . Nevertheless, with real images, our hypothesis that all points have exactly the same spectrum is obviously imprecise, and that is why we still have to enlarge Ω 0 . The finite size of the images also leads to a broadening of Ω 0 . c) Additive versus Multiplicative model: Even though it is still an approximation, the multiplicative model matches much better the actual image formation than the additive one. Indeed, the hypothesis for the latter is quite strong, although it does not appear at first glance when looking at Fig. 1b. From this figure, it is clear that the low spatial frequencies are much more significant than high spatial frequencies. However, on the one hand, depending on the observed scene, or on the area of the observed zone, the content in medium or high spatial frequencies can be much higher. This is for instance the case on Fig. 1b, where we can note that some directions are favoured. For this image, these directions are far enough from the one of the interference fringes, but for other scenes, there could be an overlap with the Ω domain assigned to the interference fringes. On the other hand, the image details at a fine scale are often very useful. This is especially our case, since, as mentioned in the introduction, our further goal is to perform stereoscopy, for which the information content at the pixel scale is quite useful. Thus, the additive model may lead to inadequate results, as it will be shown in the next subsection. Regarding the multiplicative model, we first note from equation ( 1) that, at each wavenumber, there is indeed a multiplication between the scene radiance and the interference pattern. The assumption that all points have the same spectrum is then obviously a simplification, but it is verified at first order for remote sensing images. The first reason, which is quite basic, is that the apparent spectrum of all the points is limited by the spectral sensitivity of the instrument. Besides, the atmospheric transmission and radiance is also shared by all the observed points. Lastly, even if there are spectral variations in natural scenes, we never have very sharp differences between the spectra. Therefore, the multiplicative model is more appropriate for our framework. However, it leads to non-convex restoration models, which are numerically difficult to deal with. This is why the common approach in the literature is to use the simpler additive model, which has the flaw of requiring a much larger dilation of Ω 0 than the multiplicative one. Furthermore, in regions where the image contains sharp and highly contrasted edges, the spectral support assumption for v is no longer valid and additive models introduce artifacts, as it will be seen in the next subsection.

B. Related work

Decomposing an image into an oscillatory and a structure components has been intensely explored through the cartoontexture model [START_REF] Osher | Image Decomposition and Restoration Using Total Variation Minimization and the H 1[END_REF]. Such a decomposition aims at separating an image in the following manner: w = u + v, where u is smooth and v contains the oscillating patterns of the image. This is done by representing u as a smooth function of bounded total variation and v by an object of a Hilbert space. The resulting minimisation problem is of the form:

inf u |∇u| + λ||w -u|| 2

H

Other representations for v have been proposed, depending on the properties that are sought. [START_REF] Aujol | Structure-texture image decomposition-modeling, algorithms, and parameter selection[END_REF] suggests that the texture and structure are not correlated, and minimises the oscillating term using a correlation criterion. [START_REF] Buades | Fast Cartoon + Texture Image Filters[END_REF] determines the structure and texture by identifying them with low-pass and high-pass components respectively.

In the case of interference fringe removal, the desired result is not a smooth image, but an image without fringes, keeping the texture and the edges in the scene. Inspired from cartoontexture decomposition, the idea in [START_REF] Barbanson | Correction des effets de relief en spectro-imagerie aéropotée[END_REF] is to apply an additive model to the measured images w, where this time u is the image of the scene not perturbed by fringes, and v contains solely the interference fringes. It provides an approximation that allows for a simple convex optimisation scheme.

As explained in II-A, the frequencies of the fringes are mainly contained in a small area, called Ω. Only the frequencies inside Ω are considered to be altered by fringes, so the rest of the spectrum is kept intact. To proceed, the total variation operator [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is used. It allows to smooth the image by penalising oscillations and is also useful to extrapolate the frequencies of the spectrum inside Ω, [START_REF] Malgouyres | Edge Direction Preserving Image Zooming: A Mathematical and Numerical Analysis[END_REF], [START_REF] Almansa | Measuring and Improving Image Resolution by Adaptation of the Reciprocal Cell[END_REF]. This procedure is expressed in the following convex minimisation problem, as presented in [START_REF] Barbanson | Relief Computation from Images of a Fourier Transform Spectrometer for Interferogram Correction[END_REF]:

inf u T V (u) subject to û| Ω c = ŵ| Ω c ( 10
)
where the symbol "hat" stands for Fourier transform. Although the results are good for many images, residuals can remain in areas containing sharp contrasted edges, see, e.g., Figs. 4a and4b. This image was acquired when a building with a periodic pattern was in the area of the highly contrasted fringes. In the Fourier domain, such a building creates a DC term, a fundamental and multiple harmonics, all aligned along a roughly horizontal direction. The interference fringes v then act as a modulation term, with the carrier frequency inside the Ω domain (in white on Fig. 4d). Thus, each frequency of the building pattern is replicated (see Fig. 4c or Fig. 12b).

When applying the additive model, only the Ω domain (the white rectangle on Fig. 4d) is modified. Consequently, only the replica of the DC term is regularized, but not the replicas of the fundamental or of the harmonics, which creates artifacts in the final image (Fig. 4b). Moreover, since these replicas are retained, the TV regularization creates spurious information inside Ω. This is visible on Fig. 4c, or on Fig. 12c. The fact that the scene is modulated by the frequency of the fringes corresponds, in the spectral domain, to a convolution of the interference term by the frequencies of the scene. Therefore, if the observed scene contains elements with high gradients, frequencies of the fringes will be found not only on the lower frequencies of the spectrum where Ω is set (Fig. 4d), but also on higher frequencies, as is shown in Fig. 4c.

One possibility is to use a larger domain Ω, but this is at the expense of the loss of some high frequency details in u. A better solution is to adopt a physically more accurate multiplicative image model as explained in Section III. Multiplicative data-fitting along with TV regularization has been considered in the context of Poisson denoising (see [START_REF] Aubert | A Variational Approach to Removing Multiplicative Noise[END_REF]); in the quite different context of fringe removal the solution proposed here is one-of-a-kind.

Another issue concerns the discrete total variation operator used and the 2D-Fourier transform, which are bound to produce artifacts such as staircasing, [START_REF] Moisan | How to discretize the Total Variation of an image?[END_REF]. All in all, this scheme is not accurate enough but provides an example of mainstream algorithms for comparison in Section IV.

C. Main contributions

In this work we propose a new variational approach for removing fringes from the images acquired by an ISFTS. Our approach introduces key ingredients that have a radical implication on the restored image. A summary is given next.

• Our method is robust with respect to sharp edges. In comparison with usual additive models, our multiplicative model has important advantages: it is much more accurate when restoring regions near sharp edges -a scenario where additive models cannot clearly separate the fringes v from the panchromatic image u. • Fringe frequencies estimation. Since the quality of the image restoration critically depends on the accuracy of the fringe pattern frequencies, we provide a simple and robust algorithm for their estimation based only on the image content. • Oracle estimation. The spectral constraint on v is enforced only in the direction orthogonal to the fringes. Compared to two-dimensional constraints, the onedimensional constraint alleviates some possible artifacts arising from the DFT approximation. Furthermore, using the fringe pattern frequencies for an observed scene w, we reconstruct quite a faithful oracle image which is used as initialization in the subsequent algorithms. • Frugal regularization. Our variational model is composed out of a multiplicative data term, a spectral penalty on the fringes along with a smoothed TV term on the rows of the observed image and on the columns of the fringe pattern which is more efficient in our setting and less expensive to compute than a standard two-dimensional TV regularization. • Absence of staircase artifacts. The smoothed TV avoids the well-known staircase artifacts (commonly associated to non-differentiable TV minimization) that are particularly damaging for subpixel stereo correspondence of the resulting restored images. • Convergent algorithm with explicit proximal iterations.

Our model is biconvex and we alternate iterations with proximal functions which, explicitly computed, become simple.

• Simple filtering algorithm. We provide a simple and fast algorithm for which no convergence proof is available. An extensive experimentation on a large dataset shows that the result of the fast algorithm is always extremely close to the solution of the bi-convex multiplicative model.

III. DESCRIPTION OF OUR METHOD

Assuming the multiplicative model, the relationship between the measured image w and both the panchromatic image u and the fringe image v can be summarized as

w = u • (1 + v) + η (11)
where • stands for the Hadamard (componentwise) product and η accounts for noise perturbations.

The method presented in this work was optimized on a set of nine ground-truth images where the decomposition into fringes component v and panchromatic component u is known exactly.

A. Nine ground-truth panchromatic images

To obtain noise-free simulated data, we proceeded in the following way. We used two co-registered high resolution infrared images, to deduce, for each pixel, a temperature and an emissivity. Note that these values are neither the actual temperature nor the true emissivity of the pixels, however they are realistic values, which is quite enough. From these temperature and emissivity maps, and taking into account the Planck's law and a standard atmosphere transmission and selfemission, we computed a spectral cube, corresponding to the spectral radiance at entrance pupil level. This spectral cube was then converted to a high resolution interferometric image, based on instrumental characteristics of Sieleters-MWIR (optical path difference map, spectral sensitivity, contrast and phase) and on equations 1-2. Lastly, this image was spatially filtered and downsampled, to obtain an image with the same size and resolution as Sieleters-MWIR. We thus obtained a set of images, based on the physical model. We later refer to these images as measured images, for they simulate the images obtained with an ISFTS. We also computed a panchromatic image, by suppressing the interference term in equations 1-2, and a fringe image, from the average spectrum of all objects in the scene. The panchromatic image will later be called ground truth, as it is true enough to be used as such. One example is given in Fig. 5. These images are given on a discrete grid of size m × n where m = 424 and n = 1000. In these images, the noise perturbation can be neglected, i.e., η ≈ 0 in [START_REF] Clark | Mapping of land cover in northern California with simulated hyperspectral satellite imagery[END_REF].

B. Notation

For a m × n matrix u, its rows are denoted by u i, * and its columns by u * ,i . We use the vertical mirror function M : R m×n → R 3m×n given by M u = (ū, u, ū)

T where T stands for transpose and ū is the matrix consisting of rows m, m -1, . . . , 1 of u. Its left-inverse, denoted by M † , extracts u from z = M u by taking rows m + 1, . . . , 2m of z. I m is the identity operator on R m . A vector or an image with ones of arbitrary dimension will be denoted by 1. The notation F will signify the Hamming windowed 1D Discrete Transform Fourier applied to the columns an image.

C. On the Fourier domain constraint

Our focus is on methods that constrain the Fourier transform of the image to vanish on the frequencies interval produced by the fringe pattern.

The fringe pattern v is nearly horizontal, typical of interferometric images. The 1D Fourier transform of the rows of the measured image w are not expected to change a lot except for some low frequencies. Thus the 1D Fourier transform of the rows of the measured image w and the rows of the ground truth u should be nearly the same; see Fig. 6 left and middle. The main difference is in the 1D Fourier transforms of the columns of the ground truth u and the measured image w; see Fig. 6 right. Therefore, we focus on using only the 1D Fourier transforms of the columns. a) Estimation of the spectral support of the fringe pattern: From the observation model in ( 1) and ( 2) it is known that the fringes are concentrated in a region delimited by f min and f max , the frequencies of the spectral support of v. These can be determined analytically from the optical path difference map from Eq.( 3). However, they usually shift during the acquisition. Since the restoration quality critically depends on the accuracy of f min and f max we propose a robust method for their estimation for each image w. The method is motivated by Fig. 6 and the algorithm is described next. The constraint set Ω is given by a binary m × n mask where the only null rows correspond to the frequencies in

[-f max , -f min ] ∪ [f min , f max ]
b) An oracle image: Using the measured data w and the estimated constraint set Ω, one can build a initial oracle image u (0) . It amounts to a filtering of the columns of w:

u (0) = M † F -1 Ω • F(M w) (12)
and v (0) = (w/u (0) ) -1, with "/" denoting the pointwise division.

D. Variational formulation

One notices the good reconstruction done in the oracle image u (0) as seen in Table I. An image showing the oracle u (0) is given in Fig. 15. The restoration approach is based on an inspection of the oracles u (0) and v (0) for all 9 simulated images. The main facts we notice are listed next. (a) The columns of u (0) have noticeable oscillations especially in the region of the fringes; the rows of u (0) are very close to those of the ground truth and do not need particular smoothing. (b) The 1D DFT of the fringe columns oracle are still perturbed in the frequency range -[f max , f min ] ∪ [f min , f max ] hence we need to penalize Ω • F(M w). (c) The rows of v (0) are quite more perturbed than its columns; this together with (b) leads us to apply a smoothing on the rows of v. (d) Finally, we want to keep close to the measured image w since the noise in the multiplicative model, see [START_REF] Clark | Mapping of land cover in northern California with simulated hyperspectral satellite imagery[END_REF], is very small.

Combining all these observations leads to a variational model of the form J(u, v) given by

λΦ α1 (u)+Ψ α2 (v)+β ||T (v)|| 2 2 2 +γ ||w -u • (1 + v)|| 2 2 2 (13)
where we set

T (v) := Ω • F(M • v) and Φ α1 (u) := ϕ α1 (u i,j -u i-1,j ), (14)
Ψ α2 (v) := ϕ α2 (v i,j -v i,j-1 ). (15)
The choice of ϕ is important. We focus on a convex edge preserving ϕ [START_REF] Aubert | Mathematical problems in image processing[END_REF]. A popular choice is ϕ(t) = |t| which amounts to median prior [START_REF] Besag | Digital image processing[END_REF] and to a particular discretization of the TV model. It is well known that such a ϕ introduces constant segments in the column and the row estimates of u and v, respectively, which is not desired. Therefore, we have to choose a smooth approximation of |t|. A frequent choice is ϕ α (t) = √ t 2 + α, but its derivative is involved for numerical Another choice that leads basically to the same solution is

ϕ α (t) = |t| -α log 1 + |t| α . (16)
Its derivative ϕ α (t) = t α+|t| is simple to compute which saves up to 20% running time compared to the square-root function [START_REF] Nikolova | Fast ordering algorithm for exact histogram specification[END_REF].

The parameters in (13) were fixed using the simulated images. It is quite obvious that β and γ should be large numbers by (b) and (d) because they represent "constraints". The others (quite simple from the simulated images) are given in a subsequent section.

E. Minimization Algorithm

The objective J in (13) is biconvex. A common practice for solving biconvex problems is the Alternate Convex Search where u and v are updated alternatively by fixing one of them [START_REF] Gorski | Biconvex sets and optimization with biconvex functions: A survey and extensions[END_REF]. More precisely, our algorithm follows a proximal regularization of the Gauss-Seidel scheme, see also [START_REF] Combettes | Signal Recovery by Proximal Forward-Backward Splitting[END_REF]. The objective J in ( 13) is of the form

F (u) + G(v) + H(u, v): J(u, v) = F (u) + G(v) + H(u, v), (17)
where

H(u, v) = γ 2 ||u • (v + 1) -w|| 2 2 , F (u) = λΦ α1 (u), G(v) = β ||T (v)|| 2 2 2 + Ψ α2 (v).
Let D h and D v denote the horizontal and vertical difference operators on m × n images. Using vector notation, one has

∇F (u) = λD T v ϕ α1 (D v u), ∇G(u) = βT T T v + D T h ϕ α2 (D h v)
, which gradients are used in the forward step in Algorithm 2 given next.

This algorithm is initialized with the oracle image u (0) and the corresponding v (0) . For practical reasons the algorithm below is applied with w normalized as in steps (a)-(b) below. This means we apply Algorithm 1 with the normalized input w to obtain the corresponding values f min and f max . Further, τ 1 2 Lip∇F and τ 2 2 Lip∇G , see [START_REF] Combettes | Signal Recovery by Proximal Forward-Backward Splitting[END_REF]. The Lipschitz constants here are calculated as the upper bounds of ∇ 2 F and of ∇ 2 G, respectively. Using that

D T h D h = D T v D v 4, they read as Lip∇F = 4λ α1 and Lip∇G = β + 4 α2 . Algorithm 2 INPUT: f min , f max , u (0) , v (0) , τ 1 , τ 2 , n iter , α 1 , α 2 , λ, β, γ (a) Compute c 1 := mean(w) and c 2 := 8 s.t.d.(w) (b) Normalize w ← 1 + w-c1 c2 u (0) ← 1 + u (0) -c1 c2 for k = 0 to n iter -1 do u (k+1) = prox τ1H(•,v k ) (u (k) -τ 1 ∇F (u (k) )) {minimize in u } v (k+1) = prox τ2H(u k+1 ,•) (v (k) -τ 2 ∇G(v (k) )) {minimize in v} end for OUTPUT: recovered image u = c 1 + (u niter -1) c 2
The proximity operators have a simple closed-form expression. Noticing that

∇ u H(u, v) = γ(v + 1) • ((v + 1) • u -w) and that ∇ v H(u, v) = γu • (u • (1 + v) -w), one obtains prox τ1H(.,v) (z) = z + τ 1 γ • (v + 1) • w 1 + τ 1 γ(v + 1) • (v + 1)
,

prox τ2H(u,.) (z) = z + τ 2 γu • (u -w) 1 -τ 2 γu • u ,
where division is componentwise.

The convergence of the algorithm to a critical point of J can be found in [START_REF] Nikolova | Alternating proximal forward-backward minimization for nonconvex regularized problems[END_REF].

Remark 1. Using the normalization in (a)-(b), the new measured image (or sequence of images) has mean 1 and standard deviation 0.125. This is helpful for two reasons: (1) one can process some erroneously measured images w involving negative pixels; and (2) one can fix the "free" parameters α 1 , α 2 , λ, β and γ so that they can be applied to any image recorded with a similar device. At the end of the iterations, this normalization is restituted in the recovered image u.

F. Iterative filtering to solve the multiplicative model

Let us have another look at the observations (a)-(d) given in section III-D. We want to use them in a different way to improve the oracle image u (0) in ( 12). In the light of (a) and (c), we want to minimize both Φ α1 and Ψ α2 for the same values of (α 1 , α 2 ) (assuming thus that our prior modeling on u and v in Algorithm 2 was correct enough). Furthermore, by (d) we want to closely satisfy the data-fidelity equation in [START_REF] Clark | Mapping of land cover in northern California with simulated hyperspectral satellite imagery[END_REF] and by (b) -the knowledge that the 1D columnwise Fourier transform of the fringe pattern v belongs to the interval

-[f max , f min ] ∪ [f min , f max ]. Formally, this reads as minimize Φ α1 (u) + Ψ α2 (v), subject to w ≈ u • (v + 1) and v = M † • F -1 (1 -Ω) • F(M • v) .
The first constraint is inexact (approximated) and nonconvex. . This is summarized in Algorithm 3.

Algorithm 3 INPUT: f min , f max , u (0) , v (0) , α 1 , α 2 , δ 1 , δ 2 , n iter (a) Compute c 1 := mean(w) and c 2 := 8 s.t.d.(w) (b) Normalize w ← 1 + w-c1 c2 u (0) ← 1 + u (0) -c1 c2 for k = 0 to n iter -1 do ū(k) = u (k) -δ 1 ∇Φ α1 (u (k) ) {decrease Φ(u) } v(k) = w/ū (k) -1 {preserve data-fidelity} ṽ(k) = M † • F -1 (1 -Ω) • F(M • vk ) {filtering v to approach T (v) = 0 } v (k+1) = ṽ(k) -δ 2 ∇Ψ α2 (ṽ (k) ) {decrease Ψ(v) } u (k+1) = w/(1 + v (k+1) ) {preserve data-fidelity} end for OUTPUT: recovered image u = c 1 + (u niter -1) c 2
This algorithm amounts to a nonlinear filtering. We do not have guarantees for its convergence and we use it just as a filtering, a limited number of steps. This question is discussed in detail in the next subsection.

G. Comparison of the Algorithms 2 and 3

The performance of Algorithm 2 is illustrated on Fig. 7 for the image 9 (see Table 1). On the left we observe that the PSNR value at iteration 350 equals 59.6107 and at the last iteration 500 its value is 59.6524. Further increase of the iteration number would not yield much improvement. The value of the objective J in ( 13) decreases but the decrease becomes slower. The relative error between two consecutive iterates, u (k) -u (k-1) / u (k-1) is 4.71 × 10 -7 for k = 400 while it is 3.76 × 10 -7 for k = 500. One can say that for n iter = 500 Algorithm 2 has nearly reached a minimum. We tested the filtering Algorithm 3 on all 9 images with ground truth, see subsection III-A, for 100 iterations. The best PSNR value occurred for iterations between 33 and 57, with a median value 41 and mean 41.44. What is more, its decrease after the peak is slow -for 7 cases it was less then 0.5 dB. The evolution of the PSNR for the image 9 (considered in Fig. 7) is depicted on Fig. 8 left and the relative error between two consecutive iterates on the right using a log-log scale. The PSNR value is maximum for k = 41 and the relative error is 9.87 × 10 -6 . Based on these observations, choosing n iter = 20 seems quite a safe decision. It is corroborated also by various numerical tests in section IV, as well as in the subsequent paragraph.

H. Normalisation and tools for comparing the algorithms

In many experiments we observed that starting from the Oracle initialization in [START_REF] Briottet | Military applications of hyperspectral imagery[END_REF], Algorithm 2 and Algorithm 3 yield nearly the same solution (the relative error is always less than 0.3% between the output of both algorithms as we shall see below). This observation raises important theoretical questions (when can we simplify a hard optimisation problem by a simple filtering-thresholding algorithm ?). This theoretical question is out of the scope of this report.

In section IV, we give some numerical evidence that the solution obtained by the algorithm 3 and by Algorithm 2 (see the objective in ( 13)) are nearly the same.

For this purpose, the normalization in (a)-(b) in Algorithms 2 and 3 is important. It guarantees that the images are of mean 1 and standard deviation 0.125.

To compare the different algorithms, the relative error computed in percents is computed. We can express it as follows, with u * denoting the minimizer of the objective in [START_REF] Gleeson | Characterization of a sulfur-rich Arctic spring site and field analog to Europa using hyperspectral data[END_REF] and û the result of the Algorithm 2:

ε = u * -û 2 u * 2 × 100.
The Peak Signal-to-Noise Ratio (PSNR, in dB) is also computed to quantify the difference between the algorithms. It is defined as usual:

PSNR = 10 log 10 max x |u * (x)| 2 u * -û 2 .
For the other datasets where no ground truth is available, we directly compared the outputs of Algorithms 2 and 3.

IV. EXPERIMENTAL RESULTS

A quantitative evaluation of the different algorithms are presented in this section. First, a comparison is made between the results in terms of image quality and running time. Then, some qualitative validations are presented and an assessment of the quality of stereoscopic processing results allows to further appraise the different methods. Using stereoscopy is motivated by the need for relief estimation for the correction of Fourier Transform spectrometers measures, which is the reason behind the process of removing fringes, as explained in Section I. Another reason that was be pointed out is that elevation models can be complementary to hyperspectral data for detection and analysis.

A. Available Data

Two types of data are used to compare results. First, the simulated data, as presented in subsection III-A.

The second type of data are experimental images acquired with Sieleters 2 [START_REF] Coudrain | SIELETERS, an airborne infrared dual-band spectro-imaging system for measurement of scene spectral signatures[END_REF] which is a Fourier transform spectrometer. Sieleters has two channels of acquisition for infrared MWIR and LWIR. The sequences of images were obtained from an airborne platform. The plane flew at 2000m, resulting in a ground sampling distance of 50 cm. The scenes describe rural and small urban areas in France in relatively flat areas.

B. Parameter Settings

The Algorithms 1, 2 and 3 presented in the previous sections require several parameters to be set. Different principles are used to choose the right parameter values for different parameter groups, namely: 2 Developed by ONERA/DGA

• Initialization Parameters like f min and f max are chosen by robust fitting to training data using Algorithm 1. Once these are known the initial (oracle) values u (0) and v (0) are chosen using Equation ( 12). • Numerical Parameters like τ 1 = 1.9

Lip∇F , τ 2 = 1.9

Lip∇G

were chosen in Algorithm 2 to be as large as possible while still ensuring theoretical convergence properties.

For Algorithm 3 we do not have convergence guarantees, but by analogy to the similar Algorithm 2 we chose these parameter values δ 1 = 1.99 Lip∇Ψα 1 , δ 2 = 1.99 Lip∇Ψα 2 in the same manner.

• Model Parameters like α 1 and α 2 , λ, β, γ and n iter should be chosen in such a way that the local minima of the variational model in equation ( 13) fit the ground truth solutions when available. We already know that:

β and γ should have large values in order to enforce the constraints T (v) = 0 and w = u • (v + 1) as strongly as possible without making the optimization problem too stiff (we chose β = 2500 and γ = 10 4 ); α 1 and α 2 should have quite small for better edge preservation in u and v (we chose α 1 = 5 × 10 -5 , α 2 = 5 × 10 -3 ); λ should be relatively small in order to avoid removing textures from u, but not too small otherwise some fringes may go into u (we chose λ = 0.001); n iter should be as large as possible for Algorithm 2 to reach convergence. In practice we have found that after about n iter = 500 iterations very little was improved in terms of PSNR. For the filtering Algorithm 3 the number of iterations can have a stronger impact. We have found that optimal results where obtained for n iter = 20. After the range of possible values for the different parameters was established according to the discussion above, the near-optimality of the chosen parameter values was checked with respect to the only Sieleters dataset, where a ground truth was available (namely the 9 images described in III-A). As it can be seen in Figures 9 and 10 neither Algorithm 2 nor 3 is sensitive to the exact value of these parameters and the PSNR is about 60 dB. Indeed the PSNR between the solution of Algorithm 2 varied by less than 1 dB, and 2 dB for Algorithm 3 when any of the model parameters was modified by a factor of 2.

C. Results

Experimental results and discussion on the Algorithms 2 and 3 are presented here. A comparison between results of the additive and multiplicative image models is provided. a) Validating Algorithm 3: Our first validation is to compare the results of Algorithm 3 (Fast algorithm) with Algorithm 2 (Optimization algorithm), because of the ensured convergence of this one. Table I shows the difference in performance and accuracy among the different algorithms applied to the 9 simulated images. Accuracy is measured in PSNR 4 with respect to ground truth, whereas computational performance is measured in terms of CPU 5 . We observe that for all 9 simulated images the PSNR of Algorithms 2 and 3 are nearly the same up to 0.2 dB, meaning that their solutions cannot be more than 0.4 dB apart in terms of PSNR. These values thus indicated the excellent behaviour of Algorithm 3: the difference between the PSNR obtained with Algorithm 3 and with Algorithm 2 is negligible, but Algorithm 3 is more than 20 times faster. Then, the histograms comparing the relative error ε and the PSNR(û, u * ) (following the notation in subsection III.F) between the minimizer of Algorithm 2 and the result of Algorithm 3 are presented in Fig. 11. They are computed with Sieleters images, and Algorithm 2 is used as "ground truth". The relative error appears negligible as it takes its values between 0.02% and 0.03%. The PSNR takes its lower values at 70, confirming the adecuacy of Algorithm 3 in terms of image quality compared to Algorithm 2. Furthermore, on this test, Algorithm 2 was 23.2 times slower than Algorithm 3.

Finally, we propose a qualitative evaluation. Fig. 14 offers to look at results on different landscapes and textures. Fringes were successfully removed without harming image resolution and contrast. Some slight residues can however be found in the texture in Fig. 14d. This can be explained by the fact that the texture has some patterns parallel to the fringes. In Fig. 15 images with intersecting regions were compared. The intersecting region was covered with fringes on one image, and in an area with blurred fringes on the other. Note that the latter still shows some slight fringes. We can observe that no artifact appeared and that image quality was not decreased.

b) Validating the multiplicative model: As was mentioned in subsection II-B, using an additive image formation model gives good results on many images but leaves fringe residuals on very contrasted areas (Fig. 13b). This is due to the fact that only frequencies within the small area Ω are regularised for fringe removal, when, actually, the image follows a rather multiplicative model where frequencies for the fringes are found on the whole spectrum.

If we compare the results of the additive model with the results of the multiplicative model of Fig. 4, extracted from Fig. 12a, it is clear that the multiplicative model succeeds in removing the impact of fringes on the periodic pattern of the building, as well the DC term as the fundamental and harmonics. Meantime, the additive model fails to remove those frequencies (Fig. 12). This explains why the multiplicative model leads to a significantly better image than the additive model (see Fig. 13).

Regarding quantitave results, Table I shows that the PSNR is indeed significantly better for the multiplicative model. What is more, the oracle gives a better scores also. This is explained by the fact that the additive model uses a 2D Fourier transform which induces artifacts, as the first and last columns do not match well due to the fact that the fringes are not perfectly aligned with the lines of the detector.

D. Impact on Stereoscopic Performance

As explained in the introduction, stereoscopic performance is studied as it is, in the case of Fourier transform imaging spectrometers, a likely second step for interferogram reconstruction. For this assessment we used the block matching algorithm described in [START_REF] Sabater | Meaningful Matches in Stereovision[END_REF] and the corresponding code which is available on [START_REF] Monasse | Quasi-Euclidean Epipolar Rectification[END_REF].

This algorithm provides a statistical test to decide which blocks in the first image could find a statistically significant optimal correspondance in the second image. When no correspondance could be found the algorithm outputs a cyan pixel (meaning undecidable) instead of a gray-level pixel that represents the estimated disparity. As it was shown in [START_REF] Sabater | Meaningful Matches in Stereovision[END_REF], for good stereo pairs less than 30% of the pixels are undecidable, and a larger proportion of undecidable pixels means that at least one of the two images is severely corrupted with some kind of noise or does not correspond to the same scene.

As is shown in Fig. 16, the best result is clearly obtained with the images where fringes are removed using the multiplicative model associated with the 1D Fourier transform. With the measured images, where no attempt is made to erase fringes, no information of relief could be retrieved on the lower part of the image, where fringes are most contrasted. The mainstream additive model presented in II-B gives poor results, especially in homogeneous areas. As mentioned before, this can be explained by the fact that the artifacts generated from using a discrete total variation prevent the subpixel stereo algorithm from performing well.

Numerical values corroborate the visual observations, as 85% of the image could be matched with the multiplicative model, whereas this score falls to 58% and 46% for the measured image and the additive model respectively, as seen in Table II.

From these results, we conclude that the multiplicative model gives the best results, with more coverage. It therefore brings more coherence to the image and improves results on relief estimation. 4 PSNR values are given in dB 5 CPU running times are given in seconds on an Intel QuadCore processor at 1.8 GHz running Matlab R2016b.

V. CONCLUSION

In this paper, we have introduced a novel algorithm to remove fringes from images obtained with an imaging static Fourier transform spectrometer. This algorithm exploits a numerical model close to physical model as well as the spectral properties of interference fringes to devise operators that minimise the emergence of image artifacts during processing.

We show that it is possible to find accurate solutions to the non-convex problem of using a mutiplicative image formation model. In turn, such a model allows for a better description of the interferometric image formation process. Furthermore, a good approximation to this solution can be found in a competitive computational time.

The experimental results have shown that the use of smoothed total variation regularization and adopting a 1D Fourier transform, well-adapted to the vertical periodicity of interference fringes, considerably improves results. In addition, our 1D regularization scheme retains the well-known spectral extrapolation capability of two-dimensional TV minimization, albeit with some limitations for images with horizontal patterns at the exact same frequency as the fringe pattern.

The use of the smoothed TV was also beneficial for subpixel stereo matching when compared to non-differentiable TV. However, it comes at the price of smoothing out small gradients in the measured image, with an arbitrary contrast threshold. A more accurate (even though more computationally expensive) approach that will be investigated in future research would be to reduce staircasing artifacts by using the non-differentiable Shannon Total Variation proposed by Moisan and Abergel [START_REF] Moisan | How to discretize the Total Variation of an image?[END_REF], [START_REF] Abergel | The Shannon Total Variation[END_REF].

VI. ACKNOWLEDGEMENT

This work has been partially funded by the French Research Agency (ANR) under grant No ANR-14-CE27-001 (MIRIAM).

  and A. Almansa are with LTCI, CNRS & Telecom ParisTech. A. Almansa is with MAP5, CNRS UMR 8145, Université Paris-Descartes, Sorbonne Paris Cité, 75720, Paris, France. (email andres.almansa@parisdescartes.fr) C. Barbanson and Y. Ferrec are with ONERA,91123 Palaiseau, France. (email {clara.barbanson,yann.ferrec}@onera.fr)

Fig. 1 .

 1 Fig.1. On top, an image obtained with Sieleters. Subimages are magnified to show the interference fringes overlaid on the scene. In red, an area with highly contrasted fringes. In blue, an area with blurred fringes. In green, an area with slightly contrasted fringes. The area with high contrast fringes is limited to the area around the zero optical path difference, due to the broadband observed spectra. Nonetheless, fringes are actually present over the whole image. Image (b) shows the Fourier spectrum of picture (a). The frequencies of the fringes are clearly seen as two very bright vertical lines. These frequencies constitute an area we call Ω 0 .

Fig. 2 .

 2 Fig. 2. Image processing flow of an airborne static Fourier transform spectrometer. All elements of the scene, like the one pointed out by the red cross on each picture, are matched through the sequence of images to build their corresponding interferogram as a function of optical path difference. Computing the Fourier transform provides the spectrum.

Fig. 3 .

 3 Fig. 3. On Fig. a, two images of a sequence acquired by the ISFTS CaHyD [24] are shown. Values from matching scene elements are used to construct the interferogram, as on Fig. b. The crosses on the images correspond to matches made without taking the height and geometry of the building into account. Hence, in the case of the red cross, part of the rooftop and part of the wall are used to construct a flawed interferogram. The green crosses are also mismatched, but are contained in an homogeneous area, thus providing a correct interferogram.

  a) Additive model: The simplest way to approximate the model in (1)-(2) is given by w(x, y) = u(x, y) + v(x, y) where u(x, y) = L(x, y, σ)dσ (4) v(x, y) = L(x, y, σ)F (x, y, σ)dσ (5)

  b) Multiplicative model: Assuming that all points in D have the same spectral response B(σ), one has : L(x, y, σ) = B(σ)I(x, y) By setting the constant B tot := σ B(σ)dσ, we obtain from (1)

  Fig. 4. Extract from an image with fringe residuals. The fringes correspond to the horizontal stripes. It can be seen that most have been removed in image (b). Yet, some remain on the objects with high radiometric variations. These correspond in (c) to the harmonics of the fundamental frequency which lay outside Ω.

Fig. 5 .

 5 Fig. 5. Simulated data. Top: measured image. Middle: fringe pattern. Bottom: ground truth.

3 )

 3 column mean log |F(M • w)| {where ⊗ is the Kronecker product} w = a is the nearest integer to a } c = robust fitting( w) {using a cubic curve with Cauchy weighting for parameter 1 } [f min , f max ] = the largest interval such that c > w OUTPUT:[f min , f max ]

Fig. 6 .

 6 Fig. 6. Log of 1D Fourier transforms of the ground truth (blue) and the measured image w (in red). The plot on the left shows the mean spectrum of rows. The one in the middle is a zoom of the first around the zero frequency. The plot on the right shows the mean spectrum of columns.

  For the second constraint, we use the values f min and f max estimated by Algorithm 1 with normalized input w as in (a)-(b) below. We adopt an inexact gradient projection descent approach with step-sizes satisfying δ

Fig. 7 .

 7 Fig. 7. Algorithm 2. Left: Increase of PSNR values with the iteration number. Right: the value of the objective J in (13) in a log-log scale.

Fig. 8 .

 8 Fig. 8. Algorithm 3. Left: Evolution of the PSNR values with the iteration number. Right: The relative error between two consecutive iterates in a log-log scale.

Fig. 9 .

 9 Fig.9. Sensitivity analysis of parameters for Algorithm 2: median of psnr on 9 simulated images. The plots show the variations of PSNR when the parameters are moved around their optimal value on a range from half their value to twice their value. PSNR* corresponds to the PSNR of the optimum image found by the algorithms. Globally, variations do not exceed 1dB, thus showing that the final result is robust to variations in parameter values. Parameters γ, β and the number of iterations show an optimum. This is not the case for λ, α 1 and α 2 , for which barely any difference is registered.

Fig. 10 .

 10 Fig. 10. Sensitivity analysis of parameters for Algorithm 3: median of psnr on 9 simulated images. The plots show the variations of PSNR when the parameters are moved around their optimal value on a range from half their value to twice their value. PSNR* corresponds to the PSNR of the optimum image found by the algorithms. The variations go up to 2dB, but vary very slowly after n iter = 20, our chosen value.

Fig. 11 .

 11 Fig. 11. Comparison between the results of Algorithm 3 with respect to the desired minimizer u * obtained with Algorithm 2. The first column shows histograms of the relative error as described in III.F. The second column computes the PSNR. Top line : Results are based on 100 evenly spaced images from Sieleters-MWIR. The images overlap, but it is a necessary compromise to have a large number of samples. Top line : Results are obtained with 53 disjoint (non-overlapping) images using the same instrument.

  (a) Measured image from Sieleters (b) Fourier transform from the measured image (c) Fourier transform of image obtained by additive scheme (d) Fourier transform of image obtained by multiplicative scheme

Fig. 12 .

 12 Fig. 12. Comparison between the spectra of the different schemes. On Fig. (b), the frequencies of the fringes can be seen as a very bright vertical line. Due to the image formation model, replicates of this line can be seen. The additive model, on Fig (b) manages to remove the main frequencies, but the replicas, which are outside the designated area Ω, remain. The multiplicative model removes them completely, as seen on Fig. (c).

Fig. 13 .

 13 Fig. 13. Extract of the image shown in Fig. 12. Subimage (a) comes from the measured image. The interference fringes are the horizontal lines. The image in the background corresponds to a building. Subimage (b) comes from the image obtained using the additive scheme with a 2D Fourier transform. Residues of the fringes can be observed. Subimage (c) is the result of the multiplicative algorithm, which successfully removes the residus in (b).

Fig. 14 .

 14 Fig. 14. Results for various types of landscapes. On the left column, measured images from Sieleters-MWIR, with interference fringes. On the right, the fringes have been removed using Algorithm 3.

  (a) First image: measured (b) First image: oracle (c) First image: result (d) Second image: measured

Fig. 15 .

 15 Fig. 15. Image (a) shows the measured image (from Sieleters-MWIR), (b) the oracle and (c) the image resulting from Algorithm 3 for image (a). (d) is extracted from another image which did not contain contrasted fringes (but light fringes can be observed at the top and bottom of the image) on the same area, for visual comparison.

TABLE I RESULTS

 I FROM SIMULATED DATASET WHERE THE GROUND TRUTH IS KNOWN.

	Measured image		Oracle	Fast Method	Variational Method	Additive Model
			Equation (12)	(Algorithm 3)	(Algorithm 2)	Equation (10)
		PSNR 4	PSNR 4	PSNR 4	CPU 5	PSNR 4	CPU 5	PSNR 4	CPU 5
	Image 1	21.27	53.20	59.09	1.61	59.10	40.81	38.07	1828
	Image 2	20.88	56.36	60.10	1.63	60.04	40.48	38.76	2225
	Image 3	21.12	54.33	61.38	2.40	61.32	40.58	36.94	1850
	Image 4	20.58	57.00	64.03	1.66	63.91	40.39	39.91	1684
	Image 5	24.22	59.88	67.16	1.63	66.91	40.70	35.18	1585
	Image 6	32.78	52.95	61.47	1.66	61.49	40.86	35.26	1749
	Image 7	21.31	53.29	61.13	1.70	61.15	38.95	35.73	1822
	Image 8	21.40	51.89	58.30	1.66	58.34	40.08	40.08	1990
	Image 9	20.04	50.42	59.41	1.71	59.40	38.94	36.35	2203

TABLE II STEREOSCOPIC

 II PERFORMANCE FOR THE ADDITIVE MODEL, THE MEASURED IMAGES AND THE MULTIPLICATIVE MODEL.

	Method	Percent of matched points
	Additive TV+2DFT	46%
	Measured image	58%
	Multiplicative (Algorithm 2)	85%


Paul Riot

Andrés Almansa

Yann Gousseau

Florence Tupin

A Correlation-based Dissimilarity Measure for Noisy Patches

published or not. The documents may come

A Correlation-based Dissimilarity Measure for Noisy Patches

Paul Riot, Andrés Almansa, Yann Gousseau, Florence

Introduction

In this article, we define a new dissimilarity measure between patches from natural images that are corrupted by an additive white Gaussian noise. The common whiteness assumption, although inaccurate on processed images, is realistic on RAW images on which the denoising should preferably be performed [START_REF] Colom | Nonparametric noise estimation method for raw images[END_REF]. However, most of the available algorithms do not use the white noise hypothesis to its full potential. Indeed, it is almost never enforced nor evaluated. In particular, non-local denoising methods often use the L 2 -norm to compare patches and the SNR or MSE are very widely used to evaluate denoising performances. Here, we propose to evaluate patches whiteness using their autocorrelation as an alternative to the classical L 2 -norm to make it more robust to noise and more discriminative. We demonstrate the potential of the proposed dissimilarity measure on three applications, namely: (a) similar patch detection, (b) non-local denoising, and (c) no-reference (as without oracle contrarily to SSIM or the SNR that are full-reference) denoising quality assessment.

The search for optimal dissimilarity measures between patches is not new and has already emerged in all three application areas. In a recent work [START_REF] Deledalle | How to compare noisy patches? patch similarity beyond gaussian noise[END_REF] Deledalle et al. use probabilistic arguments to derive optimal dissimilarity measures between noisy image patches for different kinds of noise. This leads to state of the art denoising algorithms in radar imaging, but in the case of Gaussian white noise their methodology leads to a disappointingly simple optimal measure, namely the L 2 norm, due to the agnostic statistical model that was adopted as a natural image prior. On the other hand, based on a relatively simple prior for natural image patches, Sabater et al. [START_REF] Sabater | Meaningful Matches in Stereovision[END_REF] propose a metric to discern similar patches from dissimilar ones that works much better for stereo vision than the L 2 norm. However this approach does not take into account the noise distribution and does not work well for noisy images.

In this work we aim to develop a dissimilarity measure that combines the advantages of both previous works. Our approach uses the autocorrelation as a measure to discriminate unstructured white noise from structured natural image features. To the best of our knowledge, only a few works proposed to take advantage of the noise statistics in this manner [START_REF] Brunet | The use of residuals in image denoising[END_REF], [START_REF] Fehrenbach | Bilevel image denoising using gaussianity tests[END_REF], [START_REF] Lanza | Variational image restoration with constraints on noise whiteness[END_REF], [START_REF] Lanza | Variational image denoising based on autocorrelation whiteness[END_REF], [START_REF] Lanza | Variational image denoising while constraining the distribution of the residual[END_REF], [START_REF] Riot | Penalizing local correlations in the residual improves image denoising performance[END_REF], [START_REF] Teuber | Denoising by second order statistics[END_REF]. However all these works focus on a variational formulation where a whiteness term is used to complement the data fitting and regularization terms in order to further constrain the solution to one that more closely complies with the first and second order statistics of the noise we want to estimate and remove.

Our approach is different and in a certain way more general. We propose a novel dissimilarity measure between patches whose interest in the denoising framework goes beyond the variational framework proposed in [START_REF] Riot | Penalizing local correlations in the residual improves image denoising performance[END_REF]. We use this measure as a substitute of the L 2 -norm to more robustly identify similar patches in non-local denoising algorithms. The new dissimilarity measure improves the performances of the Non Local Means (NLM) algorithm [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] by a large margin. It can also be applied to detect similar patches in a very noisy environment. Once again, performances are improved using the autocorrelation instead of the L 2 -norm.

Finally, we tackle the problem of denoising quality assessment. It is well known that the widely used SNR does not provide sufficient structural information. The SSIM presented in [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] was designed to offer a more perceptually pertinent dissimilarity measure. However, doubts can be raised regarding its relevance [START_REF] Dosselmann | A comprehensive assessment of the structural similarity index[END_REF] in the context of image denoising. We propose an alternative (no-reference) way to assess the quality of image denoising, relying on the autocorrelation of the difference between patches and we show on several examples that it conveys some interesting structural information.

In the first part of this article, we explain how the autocorrelation can be used to define a new patch similarity criterion. We also explain how it can be used to perform denoising evaluation. In the second part of this work, we present a few experiments demonstrating the potential of this tool for the previously defined tasks.

Proposed Similarity Criterion

In the sequel we denote the noisy image by g = f + n where f ∈ R N is the noise-free image (N is the number of pixels), n is a centered white Gaussian noise of variance σ 2 . The denoising task consists in computing an estimation u of the noise-free image. We call residual R = u -g what was removed from the noisy image. The no-reference quality assessment task consists in determining the quality of u from a statistical analysis of the residual R.

Evaluate Whiteness with the Autocorrelation

Our goal is to enrich the L 2 -norm of a noisy patch to include higher order statistics that convey more structural information. As previously stated, there exist many different tools to evaluate the whiteness [START_REF] Brunet | The use of residuals in image denoising[END_REF]. Here, we chose to use the autocorrelation that allows to detect structures and dependencies.

We can define the classical autocorrelation estimator r P on a noisy patch P of size s × s for all possible lags (l, m) ∈ {-s 2 , ..., s 2 } 2 :

r P (l, m) = (x,y) P (x, y)P (x + l, y + m) (1)
where the indices are circular. Computed on a white noise, the r P (l, m) for (l, m) = (0, 0) should hover around zero. On the other hand, computed on a structured image, the modulus of the autocorrelation coefficients would reach higher values showing dependencies. We note that r P (0, 0) is equal to the L 2norm of the patch. Choosing a circular convolution has several advantages. First, it speeds up the computation as the circular convolution can easily be computed in the Fourier domain. Second, it simplifies the calculation as every lag involves the same number of samples. Finally, it improves the estimation of the higher lags since more samples are used. Note that, in the perfect case of two identical patches corrupted with white noise, the circular convolution is harmless.

In the following, we will see how this tool can be used to enrich the L 2norm and design a new patch dissimilarity measure as well as a new denoising evaluation method.

The Autocorrelation to Define a New Patch Dissimilarity Measure

In the non-local denoising methods, the L 2 -norm of the patch difference P 1 -P 2 is a natural way to evaluate the distance. Indeed, supposing that the underlying noise free patches are equal, we know that P 1 -P 2 is a Gaussian white noise of variance 2σ 2 . Thus, the L 2 -norm can be mathematically justified as shown in [START_REF] Deledalle | How to compare noisy patches? patch similarity beyond gaussian noise[END_REF] as it can be directly interpreted as a generalized likelihood ratio test. However, it is not robust for strong noise levels.

When using the L 2 -norm for the patch distance, it has been shown for the NLM algorithm [START_REF] Buades | Non-Local Means Denoising[END_REF] that a small patch size is optimal for low noise levels. However, when σ 2 increases the distance estimation becomes unstable on small patches and s must be increased to obtain a reliable estimate. The trade-off is that there are now less similar patches. When the signal to noise ratio is low, it has been shown in [START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF] that the L 2 -norm often matches noise realizations over image content. This phenomenon shows that the L 2 distance is quite sensitive to the noise. Indeed, it only evaluates the energy of the patch difference P 1 -P 2 and does not take into account any structural information.

We show that computing the autocorrelation matrix of D is more reliable than the L 2 norm and appears to contain more structural content. This idea is applied to design a new patch dissimilarity measure using the autocorrelation estimator r P1-P2 defined previously. One can note that r P1-P2 (0, 0) = ||P 1 -P 2 || 2 gives the energy of P 1 -P 2 which is equal to the square of the L 2 norm. The other lags give information about the structural content of P 1 -P 2 . The higher the |r P1-P2 (l, m)| for (l, m) = (0, 0) are, the more structural content there is in P 1 -P 2 , meaning that the underlying noise free patches are not similar. From there, we can design a new patch dissimilarity measure that we will refer to as Whiteness Dissimilarity Measure (WDM) :

d 2 w (P 1 , P 2 ) = ||r P1-P2 || 2 2 . ( 2
)
This measure is experimentally more robust to high noise levels as it contains more structural information. It is also more reliable on big patch sizes as the autocorrelation estimation becomes more robust when more samples are available. However, we will show in the experiments that the results are also improved considering small patch sizes.

It is important to note that WDM is not a distance. Indeed, it is possible to show that it verifies the properties of non-negativity, identity of indiscernibles and symmetry. But, we can show that it does not satisfy the triangle inequality.

Denoising Evaluation Beyond the SNR

The SNR, and all its derivatives like the MSE or the PSNR, are the most commonly used metrics to evaluate and compare denoising performances. It is defined as :

SN R(u) = 10 log 10 ( ||u|| 2 ||f -u|| 2 ) (3)
Fundamentally, this measures the magnitude of the error by computing the error energy ||f -u|| 2 . It is perfectly relevant from a mathematical point of view. However, it is not always informative for a visual evaluation. Indeed, it is well known not to correlate well with visual evaluation. We can also point the fact that many state of the art denoising methods are reaching very close SNR. However, some present much more visible artifacts than others. This can't be well evaluated using the SNR. This is the reason why new metrics were developed such as SSIM [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]. However, severe criticism has emerged regarding its structural properties as demonstrated in [START_REF] Dosselmann | A comprehensive assessment of the structural similarity index[END_REF]. Both SSIM and the SNR use an oracle to estimate the denoising quality.

In this work, we aim at showing that the residual autocorrelation can be used as a complement to the SNR to evaluate the amount of structural content removed from the noisy image. Also, we claim that this can be linked to visual quality as it is sensitive to damaged textures and colored artifacts. On top of that, our metric is designed to work without any oracle by evaluating the residual whiteness. It also provides a local map to evaluate the denoising quality.

We propose to estimate a normalized residual autocorrelation on every patch R k of the residual :

A(R k ) = ||r R k || 2 2 -r R k (0, 0) 2 σ 4 k s 2 (s 2 -1) (4)
where

σ 2 k = (i,j) R k (i, j) 2
is the energy of the residual estimated on the patch R k . The normalization corresponds to the expected value of the estimator in the numerator under the hypothesis that R k is a white Gaussian noise with zero mean and variance σ 2 k . Then, we can aggregate the A(k) coefficients computed over each patch to obtain a quality map at each pixel (i, j) :

W-map(i, j) = k χ k (i, j)A(R k ) k χ k (i, j) (5)
where χ k (i, j) = 1 if and only if (i, j) belongs to the support of R k . This map gives a local estimate of how much structural content was removed from the image. A global score can be obtained with a simple normalized sum :

W-score = 1 N (i,j) W-map(i, j) (6)
with N the total number of pixels. This total score gives an estimate of the quality of denoising. When it is very close to 1, the denoising algorithm did not damage the underlying image. When it increases, it means that the residual is structured and thus some content from the image was removed. We claim that this is an interesting visual property as it is usually preferred to preserve content and only slightly denoise images, especially in the textured areas, rather than overly denoise and blur out some structural content. Also, it can be useful for later image processing tasks such as detection or recognition as we can estimate if the image content was damaged. We note that the W-score should not be evaluated on its own as it does not contain any information on the residual energy because we removed the term r R k (0, 0) 2 . Indeed, adding a white noise to an already noisy image would lead to a perfect W-score of 1 as the residual, equal to the added noise, would be white. Thus, it should be coupled with the SNR to be properly interpreted.

We chose to remove r R k (0, 0) 2 that contains the residual energy information because it is not clear what should be done with it. Furthermore, it is not simple to evaluate it at the same time as the other coefficients that contain a different information.

Experimental Results

In this section, we propose three experiments demonstrating the relevance of the proposed patch dissimilarity measure, WDM. Only a small number of visual results are shown due to length constraints, but more results are visible on the website https://wdm-ssvm2017.telecom-paristech.fr.

First, we evaluate the ability of our measure to retrieve similar patches corrupted by a strong additive white noise. Second, we show that simply replacing the L 2 norm by our measure in the Non Local Means (NLM) algorithm yields a significant performance gain. Last, the no-reference quality measure introduced in paragraph 2.3 is shown on some experiments to usefully complement the SNR.

Detecting Similar Patches with High Noise Levels

To assess the relevance of WDM to retrieve similar patches, we follow the exact same experimental protocol recently proposed in [START_REF] Deledalle | How to compare noisy patches? patch similarity beyond gaussian noise[END_REF]. We compare the respective performances of WDM, L 1 , L 2 and L 4 on a dictionary composed of 196 noisefree patches of size s = 8. The noise-free patches were obtained using the kmeans algorithm on patches extracted from the classical Barbara image. For this experiment, it is important to have different patches in the dictionary for the probability of false alarm to be well defined. The noisy patches are noisy realizations of the noise-free patches under strong noise conditions. Here we choose a white Gaussian noise with σ = 100. This experiment is only valid for strong noise levels (σ > 50) as both criteria perform perfectly for lower noises. The dissimilarity measures are evaluated for all pairs of noisy patches. The process is repeated 20 times with independent noise realizations. Numerically, the performance of the similarity criterion is given in term of their receiver operating characteristic (ROC) curve, i.e., the curve of the probability of detection as a function of the probability of false alarm. Results are given in Figure 1.

It clearly appears on this curve that WDM performs better than L 1 , L 2 and L 4 . As previously said, this is particularly true at low signal to noise ratios. Also, L 2 starts to match noise realizations over noise-free content after a certain noise level as shown in [START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF].

Improving the NLM Algorithm

The new similarity criterion defined previously can be applied to any patch based algorithm. However, for evaluation purposes we will focus on the original NLM algorithm [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. While there are many more recent and efficient denoising methods, we believe, as it is explained in [START_REF] Deledalle | How to compare noisy patches? patch similarity beyond gaussian noise[END_REF], that this algorithm provides a good way to compare similarity criteria. Note also that there are alternative ways to make this algorithm robust to noise, such as the PCA-based approach proposed in [START_REF] Tasdizen | Principal neighborhood dictionaries for nonlocal means image denoising[END_REF]. In this section, we will compare the use of this algorithm when using WDM and when using L 2 . To limit the number of parameters, we fix a search window of size W = 21. To compare the results on a fair basis the bandwidth parameter h is optimized in each case to obtain the best SNR. The SNR and SSIM results for different images using different patch sizes are shown on table 1. We can see that changing the distance from L 2 to WDM always improves the SNR by around 0.5 dB, sometimes more. The SSIM is also consistently improved.

Some visual results are displayed on figure 2. The images were corrupted using a white Gaussian noise with σ = 25. On the building image we can note that the structural components are better retrieved. Indeed, on the first detail, the lines are more contrasted and smooth. The homogeneous parts of the image are also better denoised. On the second detail, we can see that the dark areas are better denoised and again that the structural content is better preserved. Overall, the contrast preservation is much better for WDM than for L 2 . Indeed, the NLM algorithm tends to average patches that present slight differences in their mean value resulting in important losses on low frequency areas. On the contrary, simple computations can show that the WDM measure is more affected than L 2 by a constant bias, thus resulting in a smaller contrast loss. Improved performances are due to the fact that WDM is more robust to the noise. A good way to evaluate this claim is to test different levels of noise and compare the performances for low and high noise levels. We expect the L 2 norm to perform better on very low level noises while our measure would better suit the stronger noise levels. We show the corresponding SNR curve as a function of σ on Figure 3 (a), indeed supporting the previous claim of robustness. Those results can be further confirmed by another experiment. Using the noise free image to compute the weights and then using these weights to denoise a noisy image shows that, on the Lena image, the PSNR is better for the L 2 norm (34.66dB) than for the WDM measure (34.57 dB). It is nevertheless interesting to note that, even in a noise free environment, the WDM distance is still relevant as it performs only very slightly worse than L 2 .

Another interesting property of WDM is that it outperforms L 2 whatever the choice for s is, as can be seen on Figure 3 (c). While this was expected for large patch sizes, it is somehow surprising that this still occurs at small patch sizes (up to s = 3). Furthermore, the performances of WDM appear more stable with respect to the choice of s than those of L 2 , in particular when s increases. This could allow further improvements for other applications that require bigger patch sizes. Finally, Figure 3 (b) shows that the performances of NLM using WDM are somehow more stable with respect to the choice of the bandwidth than when using L 2 . This is an interesting property since the optimal value of such bandwidth parameters can highly vary from images to images, for NLM as for more recent algorithms.

Residual Evaluation

As a last experiment, we investigate the usefulness of the quality measure introduced in Section 2.3, defined by Formula (6) with s = 15. Recall that this index relies on a local (patch-based) whiteness estimation, performed on the residual image (the difference between the noisy and the denoised image). As such, and unlike the SNR or SSIM quality measures, it may be computed without any knowledge of the noise-free image. The experiment is very simple and is only provided as a proof of concept. It consists in comparing the W-scores, as defined by Formula (6), for images having the same SNR but different visual aspects. A first example is visible on Figure 4, where we compare (non-optimal) denoised versions of the Barbara image corrupted by a Gaussian additive noise with σ = 25. These images correspond to different denoising methods yielding relatively poor SNR (19.45dB±0.05dB). The result 4.(a) is relatively conservative and only a small amount of noise is removed from the noisy image. No information was removed from the image. The textures and the contrast are well preserved. As a result, no structure is detected in the residual and our criterion hovers around a value of one on the whole image.

The W-score is equal to one, which is the expected score for a white residual. On the other hand, the result 4.(b) was overly denoised, most of the textures are washed out and the final denoised image suffers heavy contrast losses. A large amount of structure was removed from the noisy image. Our evaluation criterion appears able to detect the most damaged areas (corresponding to large values of the W-map). One can note that the homogeneous areas present a very low local W-map as they were well denoised. However, the textures on the clothes are completely removed and the residual is clearly correlated, as shown by the W-map. Globally, the W-score for this image is 1.72. The result 4.(c) appears less smoothed and damaged than 4.(b) and gets a W-score of 1.23. However, it is visible on the denoised image, the residual and the W-map that some texture was removed on a few areas leading to a higher W-score. Finally, the result 4.(d) clearly damaged the textures in the image. Accordingly, the W-score is 2.72.

These experiments show that the SSIM is barely able to differentiate the results (a) and (c). However, we can see on the W-map and on the residual that a few textures were locally damaged on result (c) as opposed to (a). Furthermore the SSIM cannot tell the difference between (b) and (d) whereas W-score clearly chooses (d) as the worst result. Also, the proposed quality measure is no-reference, contrarily to the SSIM that is full-reference.

Of course, This simple experiment is in no case sufficient to assess the perceptual interest of the proposed quality measure. This task deserves further work, probably through a user study. Nevertheless, we believe this measure is of strong interest. Indeed, it is now widely admitted that the SNR is in no way sufficient to assess the quality of a restoration task. In particular, it is known to be almost blind to the destruction of textured areas. Now, the preservation of such textured areas appears as a major industrial challenge for imaging applications [START_REF] Cao | Measuring texture sharpness of a digital camera[END_REF]. While the SSIM somehow complements the SNR by paying more attention to the structure of the image, its usefulness is not always clear in the case of image denoising. Last, being no-reference, the proposed quality measure can be used as an objective function, in a way similar to the one followed in [START_REF] Riot | Penalizing local correlations in the residual improves image denoising performance[END_REF].

Conclusion

In this work, we have presented an alternative to the L 2 norm to compare noisy patches. This alternative, called WDM, appears to be more sensitive to the structural information than the L 2 norm. The WDM dissimilarity measure, relying on the autocorrelation operator, was applied to different tasks such as non-local denoising, similar patch detection or denoising evaluation. It showed interesting properties and performances in all those applications. This work opens many perspectives. While we only considered Gaussian white noise, the approach could be extended to other types of white noise. Also the proposed no-reference quality measure, showing promising preliminary results, would definitively deserve an in-depth validation. Last, many other applications could benefit from the introduced measure, such as variational denoising, stereo vision or motion tracking. Also we intend to extend this study in order to characterize a wider family of dissimilarity measures that are able to discriminate white noise from natural image structure.

Fig. 1 :

 1 Fig. 1: (left) Patch dictionary. (right) ROC curve obtained when retrieving patches under a strong white Gaussian noise (σ = 100).

Fig. 2 :

 2 Fig. 2: Building denoising with the NLM algorithm (s = 5 and W = 21)

Fig. 3 :

 3 Fig. 3: SNR performances as a function of different parameters for the NLM algorithm using WDM and L 2 (W=21). h was optimized to obtain the best SNR for (a) and (c). The scale is semi-log for (b) and the functions were aligned for an easier comparison.

  (a) and (b) were obtained with the NLM algorithm, (c) with a variationnal variant of the NLM and (d) with ROF. More examples are visible on the dedicated website.

  (a) W-score= 1.00, SSIM= 0.88 (b) W-score= 1.72, SSIM= 0.82 (c) W-score= 1.23, SSIM= 0.87 (d) W-score= 2.72, SSIM= 0.82

Fig. 4 :

 4 Fig. 4: Quality assessment of 4 denoising methods yielding the same SNR but different visual qualities. From left to right, denoised results, residuals, whiteness criterion maps (W-maps). (a) and (b) were obtained with the NLM algorithm, (c) with a variationnal variant of the NLM and (d) with ROF.

Table 1 :

 1 SNR (in dB) and SSIM performances of the NLM algorithm for L 2 and WDM patch distances with σ = 25

		L 2	WDM
	Lena	22.10 / 0.869 22.65 / 0.878
	Building 17.85 / 0.883 18.37 / 0.897
	Mandrill 19.06 / 0.821 19.29 / 0.836
	House	25.73 / 0.811 26.57 / 0.829
	Fingerprint 21.26 / 0.921 21.95 / 0.943
	Flintstones 22.52 / 0.940 23.28 / 0.951
	Boat	22.59 / 0.852 23.24 / 0.869
	Barbara 21.91 / 0.895 22.68 / 0.907