Charles-Alban Deledalle

Jérôme Gilles
email: jgilles@sdsu.edu

BATUD: Blind atmospheric turbulence deconvolution BATUD: Blind Atmospheric TUrbulence Deconvolution

Keywords:

come

Introduction

Restoration of long range images acquired through the atmospheric turbulence has gain a lot of interest these last two decades. In astronomy, a preferred solution is to use adaptive deformable mirrors which try to compensate distortions created by the atmosphere. If such systems work reasonably well for astronomical images, they are very expansive, fragile and of large dimension making them unusable for ground-to-ground or air-to-ground applications, and hard to embed in terrestrial or aerial vehicles. An alternative approach is to use image processing techniques to mitigate the atmospheric turbulence. Such methods were widely investigated these last two decades. If a fluid mechanics point of view seems a natural way to consider the turbulence, the corresponding equations are too complicated to be used within an image processing framework. Instead, Frakes et al. [START_REF] Frakes | Suppression of atmospheric turbulence in video using an adaptive control grid interpolation approach[END_REF] proposed to use a general operator point of view of the problem (usually not based on physics models). If we denote x the static original scene we want to recover, v i the observed image at time i, D i some geometric distortions due to the scattering effects of the turbulent atmosphere at time i and M a blurring kernel, a general atmospheric turbulence model impacting the observations can be written as

v i = D i (M (x)) + noise. ( 1
)
The blur M is commonly considered stationary (i.e it does not change on a short period of time) compared to the geometric distortions. The authors of [START_REF] Frakes | Suppression of atmospheric turbulence in video using an adaptive control grid interpolation approach[END_REF] then use some adaptive grid interpolation technique to reconstruct a clean image from the input sequence (see also [START_REF] Patel | An adaptive image registration technique to remove atmospheric turbulence[END_REF] for a more efficient approach). In [START_REF] Gepshtein | Restoration of atmospheric turbulent video containing real motion using rank filtering and elastic image registration[END_REF], the authors followed the same idea but replace the interpolation process with a combination of elastic registration and low-rank filtering. In [START_REF] Lemaitre | Etude de la turbulence atmosphérique en vision horizontale lointaine et restauration de séquences dégradées dans le visible et l'infrarouge[END_REF], the author proposed to use some local filtering process by block partitioning the image, however some block artifacts appear in the restored image. A system point of view is adopted in [START_REF] Micheli | A linear systems approach to imaging through turbulence[END_REF] where a Kalman filter is used to stabilize the images, followed by a Nonlocal Total Variation [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] deconvolution step to remove the blur. Several statistical arguments were use in the literature to mitigate atmospheric turbulence. In [START_REF] Micheli | The centroid method for imaging through turbulence[END_REF] and [START_REF] Meinhardt | Implementation of the Centroid Method for the Correction of Turbulence[END_REF], the authors develop a centroid based approach, while (assuming long exposure video capture), Principal Component Analysis is used in [START_REF] Li | Atmospheric turbulence-degraded image restoration using principal components analysis[END_REF]; an objective function combining temporal and spatial information is optimized by the SURE (Stein's Unbiased Risk Estimator) method in [START_REF] Song | Diffraction-limited image reconstruction with SURE for atmospheric turbulence removal[END_REF]. A method called "luckyregion/imaging" approach was proposed in [START_REF] Aubailly | Automated video enhancement from a stream of atmospherically-distorted images: the lucky-region fusion approach[END_REF], it is based on the assumption that, at specific locations and time, a neighborhood of good quality will appear with high probability. The method then aims at fusing the highest quality neighborhoods to recreate a clean image. If this method is originally motivated by the optical behavior of the atmosphere, its major drawback is the question of how to measure the quality of a neighborhood. In [START_REF] Yang | Distorted image restoration via non-rigid registration and lucky-region fusion approach[END_REF], Yang et al. proposed to combine some non-rigid registration and the above lucky-region method. Gilles et al. in [13], seek to estimate the inverse operators M -1 and D -1 i by using some blind-deconvolution algorithm after performing a diffeomorphic registration step. Based on the same strategy, Zhu et al. in [14] used a B-Spline registration within a Bayesian framework, involving a bilateral total variation regularization, to estimate the inverse operators. A two steps method using first a multiscale optical flow estimation and then the FRTAAS algorithm (First Register Then Average And Subtract) was proposed in [15,16] to obtained a restored image. The authors of [START_REF] Halder | Model-free prediction of atmospheric warp based on artificial neural network[END_REF] and [START_REF] Bai | Restoration of turbulencedegraded images based on deep convolutional network[END_REF] respectively use a generalized regression neural network and a convolutional neural network to learn turbulence induced deformations. These neural networks are then used to predict and compensate the turbulence impact. A fusion technique in a dual-tree complex wavelet domain combined with a non-rigid registration algorithm is used in [START_REF] Anantrasirichai | Atmospheric turbulence mitigation using complex wavelet-based fusion[END_REF]. This approach has two main drawbacks: 1) it requires a human intervention to select informative regions of interest, and 2) a segmentation step, whose efficacy varies, is required in the fusion strategy. In [START_REF] Xie | Removing turbulence effect via hybrid total variation and deformation-guided kernel regression[END_REF], the authors designed a multi-step algorithm which firstly estimates a reference frame via low-rank approximation; secondly, a variational model is used to enhance the quality of that reference frame. Next, this enhanced reference is used to estimate the deformation fields in the sequence which are then compensated to obtain a registered sequence. Finally a deformation guided fusion combined with a semi-blind deconvolution provides the restored image. In [START_REF] Lau | Variational models for joint subsampling and reconstruction of turbulence-degraded images[END_REF], Lau et al. proposed to minimize a functional, using different regularization terms, to perform specific subsampling of frames without any registration. A procedure based on optimal transport has been proposed in [START_REF] Nichols | Transport-based model for turbulence-corrupted imagery[END_REF]. Several similar methods also use the strategy of combining a registration step with a deconvolution step, see for instance [START_REF] Gal | Progress in the restoration of image sequences degraded by atmospheric turbulence[END_REF][START_REF] Nieuwenhuizen | Deep learning for software-based turbulence mitigation in long-range imaging[END_REF][START_REF] Arora | Key challenges and resolutions for atmospheric turbulence effect removal[END_REF][START_REF] Patel | A roadmap to mitigation techniques: Bedrock for atmospheric turbulence[END_REF][START_REF] Lau | Restoration of atmospheric turbulence-distorted images via RPCA and quasiconformal maps[END_REF].

Some articles in the literature focus more specifically on either the geometric distortion compensation or the deblurring. For instance, in [START_REF] Mao | Non rigid geometric distortions correction -application to atmospheric turbulence stabilization[END_REF][START_REF] Gilles | Mao-Gilles Stabilization Algorithm[END_REF], the authors proposed a variational model combining optical flow and a Nonlocal Total Variation regularization term while, in [START_REF] Zwart | Improved motion estimation for restoring turbulence distorted video[END_REF], the authors used a rapid control grid optical flow estimation to get rid of the geometric distortions. A comparison of several registration techniques applied to turbulence distortions can be found in [START_REF] Gezer | An investigation of registration methods on the images degraded by atmospheric turbulence[END_REF]. The deblurring aspect has been considered by Hirsch et al. in [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF] where the authors used a spatially variant deconvolution technique. Li et al. in [START_REF] Li | Atmospheric turbulence degraded image restoration by kurtosis minimization[END_REF] used a Wiener filtering approach with a blur model whose parameters are estimated by minimizing the kurtosis of the restored image. A combination of Sobolev gradients and Laplacian in a unified framework applied on a previously stabilized image was proposed in [START_REF] Lou | Video stabilization of atmospheric turbulence distortion[END_REF]. Inspired by the lucky-imaging principle, a wavelet burst technique was proposed in [START_REF] Gilles | Wavelet burst accumulation for turbulence mitigation[END_REF] to deblur a sequence of stabilized images.

All these techniques have their specific pros and cons where generally most efficient methods, in terms of restoration quality, are the most computationally expansive. It is then a tradeoff between quality and computation which drive the user to choose one method or the other given the specific constraints for his application. As mentioned earlier, none of the methods listed above use any knowledge coming from a physical model. However, in 1966, Fried [START_REF] Fried | Optical resolution through a randomly inhomogeneous medium for very long and very short exposures[END_REF] investigated the physics involved on how the atmospheric turbulence impacts the optical resolution of images. His work was revisited by Tofsted [START_REF] Tofsted | Reanalysis of turbulence effects on short-exposure passive imaging[END_REF] who provided a modern version of this model. The main result of this work is an analytical expression of the modulation transfer function (MTF), i.e the convolution kernel, corresponding to the impact of the atmospheric turbulence. In [START_REF] Gilles | Fried deconvolution[END_REF], Gilles et al. investigated the use of this MTF in a framelet based deconvolution algorithm. If the deblurring results are superior than the ones obtained by "generic" deconvolution algorithms, the proposed method has a major drawback: the expression of this kernel depends on four physical parameters. If it is reasonable to assume that three of them can be known/measured, the fourth one, called C 2 n and representing how strong the turbulence is, is less obvious to know. In [START_REF] Gilles | Fried deconvolution[END_REF], the authors proposed some sort of "brute force" approach to estimate this parameter but it is computationally intensive since many deconvolutions must be processed. In this paper, we revisit this deconvolution approach and propose a simplified version of the Fried kernel depending only into two parameters which can be estimated numerically in an efficient way using the proposed BATUD (Blind Atmospheric TUrbulence Deconvolution) algorithm.

The rest of this paper is organized as follows. In Section 2, we recall the original analytical form of the Fried kernel. We proposed our simplified expression of the Fried kernel as well as the BATUD algorithm in Section 3. In Section 4, we present numerous experiments on both simulated and real images, and show the superiority of the proposed approach against several standard blind deconvolution algorithms. We finally conclude this work in Section 5.

Fried kernel

The Fried kernel was first proposed in [START_REF] Fried | Optical resolution through a randomly inhomogeneous medium for very long and very short exposures[END_REF] and then revisited in [START_REF] Tofsted | Reanalysis of turbulence effects on short-exposure passive imaging[END_REF].

In this paper, we use the formalism of the latter article. The Fried kernel can be viewed as a combination of two terms. The first one, denoted M 0 , corresponds to a combination of the system plus atmosphere MTFs when the turbulence is negligible. The second term, denoted M SA and also called the short-term exposure MTF, models the impact (in term of blur) of phase-tilt due to the turbulence.

Denoting ω the spatial frequency (in 2D we consider an isotropic kernel and ω is the frequency modulus), M 0 (ω) can be expressed by

M 0 (ω) = 2 π arccos(|ω|) -|ω| 1 -|ω| 2 |ω| < 1 0 |ω| > 1, (2)
and M SA (ω) is given by

M SA (ω) = exp -(2.1X|ω|) 5/3 (1 -V (Q, X)|ω| 1/3 ) . (3)

If we denote by

• D: the system entrance pupil diameter (we recall from the geometrical optics that D = f /N where f is the focal length and N the optics F-number),

• L: the path length (distance from the sensor to the observed scene),

• λ: the wavelength, • C 2 n : the refractive index structure representing the turbulence intensity of the atmosphere [START_REF] Tunick | Characterization of optical turbulence (C 2 n ) data measured at the ARL A_LOT facility[END_REF], we can define the following quantities: k = 2π λ , the coherence diameter r 0 defined by r 0 = 2.1ρ 0 where ρ 0 = 1.437(k 2 LC 2 n ) -3/5 (the coherence length) and P = √ λL. Then, the quantities Q, X and V (Q, X) in eq. ( 3) can be defined by

Q = D P , X = D r 0 , (4)
V (Q, X) = A + B 10 exp - (γ + 1) 3 3.5 , (5)
where γ = log 10 (X), q = log 2 (Q) and

A =          0.840 + 0.116Σqa with qa = 1.35(q + 1.50) if q > -1.50 0.840 + 0.280Σqc with qc = 0.51(q + 1.50) if q -1.50 , (6)
and Σq = e q -1 e q +1 . The coefficient B is defined by

B = 0.805 + 0.265Σ qb with qb = 1.45(q -0.15). (7)
Finally, in the Fourier domain, Fried's MTF M F (ω) is defined by the product of M 0 (ω) and M SA (ω):

M F (ω) = M 0 (ω)M SA (ω). (8)
Practically, this kernel depends on four parameters: D, L, λ and C 2 n . The first three clearly depend on the acquisition system and the imaging scene and are eventually known. The last parameter, C 2 n , represents the turbulence intensity and, accordingly to experimental measurements [START_REF] Tunick | Characterization of optical turbulence (C 2 n ) data measured at the ARL A_LOT facility[END_REF], is generally in the range [10 -16 m -2/3 , 10 -12 m -2/3 ] corresponding respectively to weak and strong turbulence. It requires some very specific complicated equipment to measure this parameter and unfortunately such equipment are impossible to use in real practical scenarios.

Atmospheric deconvolution

A simplified formulation of the Fried kernel

In this section, we propose a general kernel model of which the Fried kernel is a particular case. Our general formulation is given by

Ma(ω) = A(ω) exp N i=1 a i f i (ω) , (9)
where N ∈ N, A : R + → R and f i : R + → R are fixed and known while a = (a 1 , . . . , a N ) T is a set of parameters to be estimated. We can notice that the Fried kernel corresponds to (the notation a b mean "a plays the role of b")

A(ω) = M 0 (ω), N = 2, ( 10
)
a 1 (2.1X) 5/3 , a 2 (2.1X) 5/3 V (Q, X), f 1 (.) = -(•) 5/3 and f 2 (.) = (•) 2 .
Note that this model also encompasses the traditional Gaussian kernel of bandwidth γ by choosing

A(ω) = 1, N = 1, a 1 = γ 2 /2 and f 1 (•) = -(•) 2 .

A blind deconvolution model

In this paper, we model the spectrum y of a blurry image as

y(ω) = Ma(ω) • x(ω) + ε(ω), (11)
where Ma is the unknown underlying Fried kernel, x the spectrum of the underlying clean image and ε is the spectrum of an additive white Gaussian noise component with standard deviation σ > 0. We consider the blind deconvolution problem of estimating jointly the image x and the vector of coefficients a solely from y. This blind deblurring problem is ill-posed such that regularization based on prior knowledge on the unknowns must be employed. To this end, we propose estimating x and a by looking for the image x , the vector a and an auxiliary variable λ * > 0, minimizing the following regularized least square problem, defined for α > 0, as

argmin x,a,λ 1 2σ 2 y -Ma • x 2 2 + R(x) + λ Ma 2 2 + α 2 λ . (12)
The first term in [START_REF] Yang | Distorted image restoration via non-rigid registration and lucky-region fusion approach[END_REF] is the standard 2 2 data fidelity term that typically results from a Bayesian perspective when assuming that the noise component is additive white Gaussian with standard deviation σ. Though this term is convex with respect to Ma and x, it is not jointly convex with respect to both variables.

The second term R(x) ≥ 0 is a regularization term enforcing x to fit some image prior knowledge. In this paper, we consider an image regularization term based on a zero-mean Gaussian Mixture Model (GMM) prior of natural image patches, as introduced in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], and defined as

R(x) = - n i=1 log K i=1 w k N (P i F -1 x; 0, Σ k ) , (13)
where w k and Σ k are learned weights and covariance matrices, and P i is the operator extracting in the spatial domain a discrete 8 × 8 zero-mean patch at pixel i. The Fourier transform and its inverse are denoted by F and F -1 , respectively. Numerically, the Fourier transform is computed via the unitary discrete Fourier transform. This regularization term enforces all spatial patches of x (with overlap) to fit the GMM prior model. In practice, we use the GMM obtained in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], composed of K = 200 components, trained by the Expectation-Maximization algorithm on a collection of 2 million patches extracted from the Berkeley Segmentation Data Set (BSDS) [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. Due to the multi-modality of GMMs, this regularization term is highly non-convex. The third and fourth terms in [START_REF] Yang | Distorted image restoration via non-rigid registration and lucky-region fusion approach[END_REF] act jointly on the structure of M a . Since R(x) always get smaller as the images become smoother, it subsequently favors the Fried kernel that over-amplifies high frequencies. To counter-balance this effect, we prevent the kernel from taking arbitrarily large values by controlling for its total energy M a 2 2 . The variable λ > 0 plays the role of a regularization parameter: the larger λ, the smaller M a 2 2 is expected to be. We observed that this parameter cannot be blindly chosen as a fixed constant for all turbulence scenarios. For this reason, we let this variable be optimized in a way that it cannot get too small as controlled by the last term in [START_REF] Yang | Distorted image restoration via non-rigid registration and lucky-region fusion approach[END_REF] and parameterized by α > 0. Though convex in λ and Ma, this problem is neither convex in a, nor jointly in λ and a.

A blind deconvolution optimization procedure

Due to the non-convexity of the joint optimization problem, we adopt in this paper a greedy approach that we coin Blind Atmospheric TUrbulence Deconvolution (BATUD). BATUD attempts to minimize (12) by using the following alternate minimization scheme

x t+1 ∈ argmin x 1 2σ 2 y -M a t • x 2 2 + R(x), (14)
λ t+1 ∈ argmin λ λ M a t 2 2 + α 2 λ , (15)
a t+1 ∈ argmin a 1 2σ 2 y -Ma • x t+1 2 2 + λ t+1 Ma 2 2 , ( 16
)
where t ≥ 0 is the time step and the initialization for a 0 will be discussed in Section 4.1 (note that we do not need to initialize x 0 and λ 0 ). Problem ( 14) can be seen as a non-blind deconvolution problem for the fixed convolutional kernel M a t . This deconvolution problem, with R(•) the GMM patch based prior given in (13), can be solved with the Expected Patch Log-Likelihood (EPLL) algorithm [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. EPLL is an iterative algorithm that repeats three steps: (i) extracting all image patches, (ii) filtering all patches assuming a GMM patch prior, and (iii) performing a Wiener-like deconvolution step that combines all previously estimated patches. In this paper, we use Fast EPLL (FEPLL), a fast approximation of EPLL introduced in [START_REF] Parameswaran | Accelerating GMM-Based Patch Priors for Image Restoration: Three Ingredients for a 100× Speed-Up[END_REF]. The corresponding code, provided by the authors, is freely available on GitHub * .

Problem (15) is convex and has a trivial closed form solution given by

λ t+1 = α M a t 2 . ( 17
)
Equation [START_REF] Halder | Model-free prediction of atmospheric warp based on artificial neural network[END_REF] shows that if M a t 2 gets larger with t, λ t decreases until a trade-off is achieved. Problem ( 16) can be recast as the minimization of

E(a) = 1 2 |Ma(ω)x(ω)| 2 -x(ω) * Ma(ω) * y(ω) + λσ 2 |Ma(ω)| 2 dω, (18)
where * denotes the complex conjugate, and we removed the dependency of x and λ with t for the sake of notation simplicity. Equation ( 18) can be minimized by Newton's descent consisting in iterating the following update

   a 1 . . . a N    ←    a 1 . . . a N    -     ∂ 2 E(a) ∂a1∂a1
. . .

∂ 2 E(a) ∂a1∂a N . . . . . . ∂ 2 E(a) ∂a N ∂a1
. . .

∂ 2 E(a) ∂a N ∂a N     -1     ∂E(a) ∂a1 . . . ∂E(a) ∂a N     , (19)
where

∂E(a) ∂a i = f i (ω)x(ω) * Ma(ω) * [Ma(ω)x(ω) -y(ω)] + λσ 2 f i (ω)|Ma(ω)| 2 dω, (20)
and

∂ 2 E(a) ∂a i ∂a j = f i (ω)f j (ω)x(ω) * Ma(ω) * [2Ma(ω)x(ω) -y(ω)] + 2λσ 2 f i (ω)f j (ω)|Ma(ω)| 2 dω . (21)
In practice, at every time step t > 1, we use a warm start strategy to initialize Newton's descent by using the solution a t-1 obtained at the previous time step. By doing so, we observed a very fast convergence of this numerical scheme. Though the objective function might have regions of non-convexity where its Hessian matrix is not non-negative definite, we always observed a stable behavior of this algorithm by using such an initialization. Nevertheless, as a safeguard, if any numerical issues were arising from the resolution of ( 19), the optimization of a would be instead carried out by the simplex search method of Lagarias et al.

[43] † .

Stopping criteria

Regarding the overall optimization scheme, at each time step t ≥ 1 we compute the relative errors between two consecutive estimates for the three quantities to be optimized as

δ t x = x t -x t-1 2 x t-1 2 , δ t λ = |λ t -λ t-1 | |λ t-1 | , (22)
and

δ t a = M a t -M a t-1 2 M a t-1 2 . ( 23
)
We consider that the algorithm has converged as soon as these relative errors are under a threshold τ during at least ∆ consecutive iterations. In other words, we stop the algorithm at time step t when

max t-∆≤t ≤t max(δ t x , δ t λ , δ t a ) < τ. (24)
Similarly, for Problem (16), we run Newton's descent by iteratively applying eq. ( 19) until two consecutive estimates for a achieve a relative error smaller than τ .

Experiments

In this section, we conducted several studies and experiments to understand the behavior and demonstrate the performance of our proposed blind atmospheric turbulence deconvolution algorithm. For the sake of reproducible research, our Matlab implementation of BATUD has been made open-source and is freely available online at:

https://www.charles-deledalle.fr/batud

Parameter setting and initialization

In practice, we choose α 2 = γn/σ 2 with a constant parameter γ = 8.5 • 10 -3 that is kept fixed in all our experiments, and n being the number of pixels. The initial blur parameter a 0 has been chosen empirically to correspond to the Fried kernel with X = 2.75 and V = 1.06. FEPLL [START_REF] Parameswaran | Accelerating GMM-Based Patch Priors for Image Restoration: Three Ingredients for a 100× Speed-Up[END_REF] was used with default parameters which consists in performing five iterations of patch filtering and aggregation.

On experiments involving real data, we estimate the noise standard deviation σ using a classical procedure based on the Median Absolute Deviation (MAD) of detail coefficients, as explained in [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF].

Otherwise specified, all experiments have been carried out by performing at most 150 iterations of the proposed alternate optimization scheme with the stopping criterion parameterized by τ = 2 • 10 -4 and ∆ = 3.

Experimental study on the simulated scenarios

We start to assess the proposed approach on simulated blurs applied to the Kodak dataset * composed of 23 images of size 768 × 512 (or 512 × 768) and encoded on 8-bits and converted to grayscale. We apply three types of blur, by using the original Fried kernel model given in Section 2, corresponding to different levels of turbulence: Weak (X = 2.25, V = 1.29), Medium (X = 3.00, V = 1.00) and † using the Matlab function fminsearch. * http://r0k.us/graphics/kodak/ Strong (X = 5.00, V = 0.80). In addition we add white Gaussian noise with three different levels: σ 1 = 1.3, σ 2 = 1.9 and σ 3 = 2.5.

Note that this choice of levels of turbulence and noise captures the range of scenarios we observed in practice. Our deconvolution procedure is then applied to each simulated image. Quantitative performance measurements are computed by calculating the Peak Signal to Noise Ratio (PSNR) criterion defined in decibels (dB) by

PSNR(x, x) = 10 log 10 255 2 1 n x -x 2 2 , (25)
where x denotes the restored image and x the reference image (both images are assumed to be encoded on 8-bits). We also consider Structural Similarity Index (SSIM) introduced in [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]] that provides a quantitative measure of quality between 0 and 1, and supposedly being more faithful to visual perception than the PSNR is. The larger these deblurring criteria are, the better the deblurring quality is considered to be. We also assess the performance of the algorithm to estimate the original kernel used in the simulations. To quantitatively measure this performance, we compute the relative error between the original and estimated kernels (we will call it the Kernel Relative Error, i.e KRE) defined by

KRE(â, a) = M â -Ma 2 2 Ma 2 2 , (26)

4.2.1

Behavioral study of the proposed algorithm: In order to assess the quality of BATUD (our blind deconvolution approach), we computed the oracle images, i.e the deconvolved images using the real kernel which was used in the simulation (this is basically equivalent to perform a non-blind deconvolution). In addition, we also computed oracle kernels, i.e the estimated kernels using the original clean image used in the simulation. The performance of these oracle estimators provides in some sense an ideal baseline for the expected performance of our blind deconvolution algorithm.

Comparisons in terms of PSNR, SSIM and KRE (averaged over the 23 images) are given in Table 1 for the three levels of turbulence and the three levels of noise we considered. Visual comparisons on two images of the dataset are provided in Figs. 1 and2 for the three considered levels of turbulence, and the two levels of noise σ 1 and σ 3 , respectively. Additional results obtained from the BATUD algorithm on different simulated scenarios are also provided in Fig. 3.

For Medium and Strong turbulence levels, we can see that the estimated kernels have, in average, relative errors of no more than 4%. Given that the Fried kernels are well estimated in this case, the deblurring results are close to the one provided by the non-blind oracle, with a drop of PSNR smaller than 1dB and a drop of SSIM smaller than 0.03. Visual inspections of the results of blind and non-blind deconvolution, in Medium and Strong turbulence levels (second and third column in the figures), show that these drops of quality measures are indeed visually imperceptible.

In the Weak turbulence level, the estimation of the Fried kernel is less accurate than in the other two scenarios. The relative errors of BATUD increase around 16%. Even the oracle kernel estimator, that was providing relative errors smaller than 1% in the other two scenarios, now provides relative errors around 10%. Inspecting our estimation of the Fried kernel in Fig. 1 (first column) shows that our approach tends to slightly under-estimate the high-frequency responses of the filter (leading to a KRE of about 23% in this case). As a result the deblurring procedure over-amplifies these high frequencies leading to structures appearing sharper than they are supposed to be. Small Gibbs-like artifacts can also be observed around high-contrasted objects. These effects can be seen in Fig. 1 (first column) by comparing the ability of the BATUD and the non-blind deconvolution at recovering the waves patterns, the folds and the numbers on the sail. Though barely visible, Table 1 shows that these over-sharpening effects and Gibbs artifacts contribute to a loss of PSNR of about 7dB compared to the non-blind oracle. Nevertheless, compared to the oracle, these artifacts do not influence much the SSIM, the drop being smaller than 0.01, hence reflecting our visual inspections.

As an additional study, Figs. 1 and 2 (last rows) provide the evolution of the relative errors (in percentage) between two consecutive estimates with respect to the iteration index t of BATUD. These quantities correspond to δ t x , δ t λ and δ t Ma defined in [START_REF] Nichols | Transport-based model for turbulence-corrupted imagery[END_REF]. The curves are computed until the algorithm stops when the stopping criterion in [START_REF] Nieuwenhuizen | Deep learning for software-based turbulence mitigation in long-range imaging[END_REF] is satisfied. We first observe that under any turbulence and noise scenario, these relative errors converges to the required level of tolerance 0.02% in at most 150 iterations. Under Weak and Medium levels of turbulence, the procedure takes at most 70 iterations and the relative errors for λ and Ma appears to be monotonically non-increasing. Under a Strong level of turbulence, BATUD behaves differently. It converges in about 150 iterations and presents two different phases. The first ∼65 iterations can be seen as an exploration phase where λ and Ma keep changing (two consecutive estimates differ more than 25%) until they reach a probing configuration. Afterwards, the procedure enters an exploitation phase in which this promising configuration is iteratively refined until no significant changes arise.

Given that one iteration of BATUD takes about 1.4s to process a Kodak image * , the computation time ranges from 56s to 3min 30s per image depending on the underlying level of turbulence.

4.2.2

Comparisons with state-of-the-art blind deconvolution algorithms: In this section, we compare our algorithm against several standard and state-of-the-art algorithms. The first algorithm, directly available in Matlab via the deconvblind function, is based on a Maximum Likelihood estimator to recover both the convolution kernel and the restored image [START_REF] Holmes | Light microscopic images reconstructed by maximum likelihood deconvolution[END_REF][START_REF] Hanisch | Deconvolutions of Hubble space telescope images and spectra[END_REF] (we will denote it MLE). The second algorithm (denoted OBD, available at † ) described in [START_REF] Harmeling | Online blind deconvolution for astronomical imaging[END_REF] (this algorithm was initially designed to process a set of images but the authors mentioned that it also works using a single image), minimizes a loss function based on 2 assuming a Poisson model on the data and the positivity of both the data and the kernel. In [START_REF] Kotera | Blind deconvolution with model discrepancies[END_REF], an estimation method (denoted KOT, whose code is available at ‡ ) based on a modified of the Automatic Relevance Determination model is proposed. The authors include a Student's t-prior as well as some Gamma distribution within a variational Bayesian framework to solve this model. Cai et al.,[START_REF] Cai | Framelet-based blind motion deblurring from a single image[END_REF], propose a blind deconvolution algorithm using a variational approach (we will denote this algorithm CLS and use the code available at § ). The authors use a sparsity assumption in a framelet domain for both the reconstructed image and the kernel, an 2 constraint is also added to the kernel. In [START_REF] Shan | High-quality motion deblurring from a single image[END_REF], Shan et al. propose a unified probabilistic model to estimate both the kernel and the latent image (this algorithm will be denoted SJA and is available at ¶ ). They use a Gaussian model for the noise, an exponential prior for the kernel and a combination of a local and global priors for the image. This algorithm has been used in actual turbulence mitigation techniques in [14,[START_REF] Xie | Removing turbulence effect via hybrid total variation and deformation-guided kernel regression[END_REF] and [START_REF] Lau | Restoration of atmospheric turbulence-distorted images via RPCA and quasiconformal maps[END_REF]. Finally, the last algorithm (denoted ZWZ and available at ) we consider combines probabilistic and variational models (an 2 and p, with 0 < p < 1, * Computation time evaluated on an Intel(R) Xeon(R) CPU E5-2640 v2 @ 2.00GHz, using one single core only.

† http://pixel.kyb.tuebingen.mpg. constraints are used to enforce a sparsity prior). A Maximum A Posteriori approach is then used to numerically solve this model. For all these algorithms, we used the default parameters recommended by the respective authors.

We run these algorithms on the Kodak dataset used in Section 4.1 and re-utilize the PSNR and SSIM metrics (unfortunately some of the available codes do not return the estimated kernels thus we were not able to use the KRE metric). Table 2 provides the averaged PSNR and SSIM for each algorithm. It is clear that our method outperforms the state-of-the-art algorithms.

We observe that CLS and ZWZ locally modify the image dynamics and introduce some ringing artifacts (especially in the case of weak blur). Moreover, CLS also performs poorly (i.e limited deblurring) on the weak turbulence case, this is probably due to the fact the algorithm is made to remove motion blur hence it is not capable of estimating small isotropic type kernels. The SJA algorithm is very sensitive to the presence of noise that leads to the introduction of a lot of artifacts in such case. Kotera is not capable to deal with strong blurs and sometimes it tends to do some image sharpening than actual deblurring. Some of the observed artifacts for the CLS, KOT, SJA and ZWZ algorithms are illustrated in Fig. 4. The MLE and OBD algorithms are inefficient at estimating the underlying kernel and therefore almost do not perform any deblurring at all. If these algorithms perform reasonably well on some images, in average, they perform worse than keeping the actual blurry images! These results clearly show that a deconvolution model based on a physical model is preferable to a blind algorithm.

Moreover, note that the CLS and ZWZ algorithms take about 8min and 5min to process a Kodak image * , respectively, while BATUD takes only about 2min.

Experiments on real images

In this section, we experiment the BATUD algorithm on real images coming from several datasets: the OTIS dataset [START_REF] Gilles | Open Turbulent Image Set (OTIS)[END_REF], the NATO-SET152 New Mexico field trials, the U.S. Army RDECOM CERDEC Night Vision & Electronic Sensors Directorate data and the test images available from the IPOL demonstration webpage † . The sequences from these datasets are first stabilized using the Mao-Gilles algorithm [START_REF] Mao | Non rigid geometric distortions correction -application to atmospheric turbulence stabilization[END_REF][START_REF] Gilles | Mao-Gilles Stabilization Algorithm[END_REF] to obtain the input image y. These images and their deconvolved versions, as well as the estimated kernel are provided in Fig. 5. We can observe that the algorithm is capable of removing the blur and retrieve more details in the images (the texts are easier to read, most textures are recovered and thin structures * We provide only approximate values since the CLS and ZWZ implementations perform the reading and writing of images within their code. † http://demo.ipol.im/demo/46/ are sharpened). As expected, we can also notice that the stronger the turbulence, the narrower are the estimated kernels.

In Fig. 6, we compare our algorithm against the CLS and ZWZ blind deconvolution algorithms mentioned in Section 4.2.2 (we choose only these two algorithms since they are the ones giving the best results in medium or strong turbulence). On the first image, we can observe that the CLS algorithm provides a slightly better result (especially in the textured regions) than ZWZ and BATUD while for the two other images BATUD gives superior results. In particular, on the second image, the word "A BEAR" is easier to read from the BATUD output. The same observation can be made on the third image where BATUD is capable of better restoring smaller text than the other algorithms.

Conclusion

We introduced BATUD a new algorithm for blind atmospheric turbulence deconvolution. Unlike concurrent algorithms, ours relies on a physical of blur induced by atmospheric turbulence, namely the Fried kernel. Our approach is fast and fully automatic as it estimates jointly the deconvolved image and the Fried kernel while auto-adjusting the regularization parameter. Experiments on simulated and real data have demonstrated the improved performance of BATUD in several atmospheric turbulence scenarios against state-ofthe-art approaches. In future works, we plan to develop algorithmic recipes to make our implementation even faster and to port the implementation on Graphics Processing Units (GPUs). We also aim at adapting the approach to low-light vision systems and hyperspectral imagery. Extensions of this work to submarine imagery by using appropriate physical models of underwater turbulence are also under investigation.

Fig. 1 :

 1 Fig. 1: (a) Blurry image y under medium noise level σ = 0.005. (b) Oracle deconvolved image with BATUD. (c) Deconvolved image BATUD. (d) True and estimated Fried kernel with the proposed approach (the zero frequency is at the center and amplitudes are displayed in log-scale). (e) Relative errors between consecutive estimates with respect to the time step t. (a-c) PSNR and SSIM values are displayed on the bottom-left. (d) KRE is displayed on the bottom-left.

Fig. 2 :

 2 Fig. 2: (a) Blurry image y under strong noise level σ = 0.01. (b) Oracle deconvolved image with FEPLL. (c) Deconvolved image with BATUD. (d) True and estimated Fried kernel with the proposed approach (the zero frequency is at the center and amplitudes are displayed in log-scale). (e) Relative errors between consecutive estimates with respect to the time step t. (a-c) PSNR and SSIM values are displayed on the bottom-left. (d) KRE is displayed on the bottom-left.

Fig. 3 :

 3 Fig. 3: BATUD results on different simulated scenarios of blur and noise from the Kodak dataset. PSNR and SSIM values are displayed on the bottom-left.

Fig. 4 :

 4 Fig. 4: Illustration of artifacts appearing with the CLS, KOT, SJA and ZWZ blind deblurring algorithms.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Results obtained on height images with real atmospheric turbulence. (a) Blurry images y with different levels of turbulence. (b) Deconvolved images. (c) Estimated Fried kernels (the zero frequency is at the center and amplitudes are displayed in log-scale).

Table 1

 1 Simulated atmospheric turbulence. PSNR, SSIM and Kernel relative

	errors (KRE) averaged on the Kodak dataset
			Weak		Medium			Strong
		σ1	σ2	σ3	σ1	σ2	σ3	σ1	σ2	σ3
	PSNR	Blurry BATUD 34.53 34.63 34.12 32.37 31.27 30.52 27.44 27.13 26.88 30.45 30.28 30.05 26.75 26.67 26.57 24.39 24.35 24.29 Oracle 41.21 38.97 37.36 33.29 31.99 31.15 27.71 27.38 27.14
	SSIM	Blurry BATUD 0.925 0.902 0.879 0.863 0.822 0.789 0.694 0.658 0.629 0.858 0.834 0.810 0.705 0.686 0.667 0.482 0.469 0.454 Oracle 0.932 0.905 0.882 0.868 0.831 0.800 0.709 0.674 0.646
	KRE	BATUD 0.165 0.162 0.164 0.023 0.032 0.040 0.020 0.024 0.028 Oracle 0.064 0.107 0.147 0.003 0.005 0.008 0.000 0.000 0.001

Table 2

 2 Comparison with state-of-the-art blind deconvolution algorithms.PSNR and SSIM averaged on the Kodak dataset.

				Weak			Medium			Strong
			σ1	σ2	σ3	σ1	σ2	σ3	σ1	σ2	σ3
		Blurry	30.45 30.28 30.05 26.75 26.67 26.57 24.39 24.35 24.29
		BATUD	34.53 34.63 34.12	32.37 31.27 30.52	27.44 27.13 26.88
	PSNR	CLS KOT MLE	26.26 26.15 25.82 30.14 29.12 28.29 26.62 26.43 26.14 27.70 27.57 27.57 28.18 28.02 27.97 25.89 25.79 25.68 30.45 30.28 30.05 26.75 26.67 26.57 24.39 24.35 24.29
		OBD	30.45 30.28 30.05 26.75 26.67 26.57 24.39 24.35 24.29
		SJA	26.52 25.49 23.86 26.47 25.93 24.08 24.89 24.09 22.28
		ZWZ	31.26 31.03 30.19 28.51 28.13 27.78 25.83 25.76 25.63
		Blurry	0.943 0.940 0.935 0.863 0.859 0.854 0.734 0.731 0.726
		BATUD	0.987 0.984 0.978	0.966 0.949 0.934	0.869 0.853 0.839
	SSIM	CLS KOT MLE	0.864 0.857 0.847 0.949 0.931 0.912 0.851 0.839 0.824 0.923 0.915 0.907 0.913 0.904 0.896 0.820 0.809 0.798 0.943 0.940 0.935 0.863 0.859 0.854 0.734 0.731 0.726
		OBD	0.943 0.940 0.935 0.863 0.859 0.854 0.734 0.731 0.726
		SJA	0.909 0.877 0.840 0.844 0.818 0.724 0.780 0.718 0.562
		ZWZ	0.910 0.904 0.894 0.914 0.900 0.887 0.807 0.800 0.790

Acknowledgments

This work was partially supported by the Air Force Office of Scientific Research under the grant number FA9550-15-1-0065. The authors want to thank the members of the NATO SET156 (ex-SET072) Task Group and the U.S. Army RDECOM CERDEC Night Vision & Electronic Sensors Directorate for the opportunity of using their respective data. We also want to thanks the different authors who made their code publicly available.