NON RIGID GEOMETRIC DISTORTIONS CORRECTION - APPLICATION TO ATMOSPHERIC TURBULENCE STABILIZATION

A BSTRACT . A novel approach is presented to recover image degraded by atmospheric turbulence. Given a sequence of frames affected by turbulence, we construct a variational model to characterize the static image. The optimization problem is solved by the Bregman iteration and operator splitting method. Our algorithm is simple and efﬁcient, and can be easily generalized for different scenarios.


INTRODUCTION
In the last decade, long range imaging systems have been developed to improve target identification.One of the main visual effect is distortions due to atmospheric turbulence (known in the literature as "image dancing").It may occur in many other scenarios.For example, underwater imaging systems are subject to scattering effects and video shooting in summer suffers from hot air near the ground, and so on.Weak turbulence does not really affect human observers, but it can cause problems for an automatic target recognition algorithm because the shape of the object may be very different from those learned by the algorithm.Fig 1 shows some examples obtained by a camera in real scenarios.For each video we arbitrarily choose three frames to display here.
Previous methods have been developed to deal with the turbulence effect in astronomical images.In [19], local filters (Wiener filter, Laplacian regularization and so on) were utilized and local properties were obtained by block partitioning of the image.As a result, some block artifacts appear on the restored images.
An interesting work about turbulence modelization for mitigation algorithms was made by Frakes [11,12].The authors modeled the turbulence phenomenon by using two operators: (1) where u is the static original scene we want to retrieve, f i is the observed image at time i, H is a blurring kernel, and D i is an operator which represents the geometric distortions caused by the turbulence at time i.Based on this model, the authors of [14] proposed a scheme to evaluate the H −1 and D −1 operators.The H −1 operator is obtained by blind deconvolution, while the correction of the geometrical distortions D −1 is computed by an elastic registration algorithm based on diffeomorphic mappings.This approach gives nice results but it has two main drawbacks.First, it is time consuming to perform the calculations due to the two iterative processes involved in the algorithm.Secondly, the performance is sensitive to the choice of the parameters.Another kind of approach for this problem is to utilize the Kalman filter, which is a statistical tool that recovers a static object from a time series of observations.In [27], the authors successfully use this filter in the turbulence reconstruction problem.However, this method requires a strong time dependence of the frames, therefore the frame rate has to be sufficiently high.It treats the warped frames ordered in time as governed by fluid dynamics, and thus can be characterized by time-dependent differential equations, which is not a practical assumption in some applications.
More recently, some efforts were made to propose new mitigation algorithms.In [20], assuming long exposure video capture, the authors propose to use a Principal Component Analysis to find the statistically best restored image from a sequence of acquired frames.In [2], the authors use the assumption that for a fixed location in the image, its neighborhood has some high probability to appear with better quality through the time.Then the restored region is a fusion of the best ones.Some spatially variant deblurring was proposed in [18] True scene The model of deformation used in this paper.
but the algorithm doesn't specially address the problem of geometrical distortions.Close to the modelization of Frakes [11], in [30] the authors propose to inverse the geometric distortions and the blur.They use some B-Spline registration algorithm embedded in a Bayesian framework with bilateral total variation regularization.
The goal of this paper is to propose new investigations on this problem, principally on the correction of geometrical distortions.We develop a unified framework which uses both an optical flow scheme to estimate the geometrical distortion, and a nonlocal-TV based regularization process to recover the original observed scene.The paper is organized as follows: in section 2, we describe the basic model used throughout the whole paper.Section 3 deals with Bregman iteration and the operator splitting operator used in the optimization process.Section 4 provides the whole algorithm and implementation details.Section 5 present many numerical results obtained by the proposed method on real data.We conclude in section 6.

OUR BASIC MODEL
We denote the observed image sequence as {f i } i=1,...,N and the true image that needs to be reconstructed as u.We assume (2) where φ i corresponds to the geometric deformation on the i-th frame (let remark that the φ i are the deformations between the true image and the observed frame i and not the continuous movement flow from frame to frame, see Fig. 2).
If we fix φ i , u(•) → u(φ i (•)) can be treated as a linear operator on u, so we can write the right hand side of (2) as where Φ i is the linear operator corresponding deformation φ i , then (2) becomes (4) which gives our fidelity equations in the model.On the other hand, it is reasonable to assume that our image has certain regularization features.For natural image the total variation has been proved to be satisfactory for avoiding noise while preserving sharp edges in image [25].There are many other modern models such as nonlocal total variation that has been developed and investigated throughly as well [5,13].If we denote the regularization term of the image as J(u), then we can formulate our problem as (5) min We want to remark that rather than investigating the time-dependent fluid dynamics model behind the warping effect, which is often extremely complicated, we simply treat the frames as arbitrary samples of the image after random morphing without using any sequential information of the video frames.This can greatly simplify the model and make the model available even if the sequential relevance is not strong enough in practical data.
The regularization term of u, J(u), has many different choices.Most notably, total variation based optimization [25] have been successful for edge preserving regularization.However, the use of pure total variation models for realistic images have been shown to produce artificial patches; this is due to the choice of the bounded-variation space and the corresponding total variation norm.More recently, the nonlocal means regularization model [5,13] has been introduced which modifies the intensity of a pixel by considering the nearby pixel values with similar patterns.The basic assumption behind is that a natural image contains repeating structures instead of repeating pixels.This method has been proved to be successful to remove artifacts while keeping the regular pattern and texture contained in the image and has been extended to include variational method using functionals with nonlocal regularization, and proven superior to many other image regularization methods as it considers the large-scale structure of the image beside the local differences between pixels, which makes it capable of preserving important detailed features in an image while removing artifacts effectively.For these reasons, in this paper we utilize the nonlocal regularization.Thereafter, for the reader's convenience, we recall the expression of J(u) in the nonlocal regularization case.A detailed introduction of this regularization method and its numerical implementation can be found in [13].The nonlocal Total Variation (NLTV) is defined by ( 6) where the weight w(x, y) corresponds to the similarity between patches centered on pixel x and y.The more the patches are similar, the more they are taken into account into the regularized image.
There are many problems which have similar forms with (5) that have been studied in many recent references.For example, in [8,9,17], the blind deconvolution problem was modeled as (7) min where H(k) is another regularization term on the unknown convolution kernel k.This kind of model can be solved by the alternative optimization method, i.e. optimizing over different variables alternatively.We want to remark that ( 7) is a convex problem for u and k respectively, but not convex for the joint variables, and model (5) has the same feature.The nature of non-convexity makes it hard to analyze the convergence behaviour for the optimization procedure, but the alternative optimization method at least guarantees that the functional is always decreasing over the iterations.In our model ( 5), if we have a good guess on u, then the optimal φ i can be estimated by (2) via certain optical flow algorithms (e.g. the methods developed in [3,4,26]).On the other hand, given fixed {φ i } the model ( 5) can be solved as a constrained problem, which we will discuss in the next section.
3. ALGORITHM 3.1.BREGMAN ITERATION.The Bregman iterative method was originally introduced to the image processing by [23].It solves the following constrained optimization problem (8) min by solving the following series of problems ( 9) with f 0 = f and A a linear operator (the deformations in our case).It has been shown that this iteration converges to the solution of (8).The Bregman iterative method ( 9) is actually equivalent to alternatively descent the primal variable and ascent the dual variable of the Lagrangian of ( 8), as pointed out by many papers, e.g.[28,29].9) is an unconstrained problem.This formulation has appeared in many practical imaging or signal processing problems, see [17,22,28] for examples.It can be solved by the forward-backward operating splitting method, which was first proposed by Lions and Mercier [21] and Passty [24] and generalized by Combettes and Wajs [10].The scheme can be described as follows: to solve the unconstrained problem (10) min

OPERATOR SPLITTING METHOD. The first step in (
we want to find the u such that 0 ∈ ∂J(u) + λA ⊤ (Au − f ), where ∂J(u) denotes the sub-derivative of J(u).This leads to the following fixed point algorithm: This method is practically efficient.The first line is called the forward step which is nothing but the gradient descent of Au − f 2 .The second line is called the backward step and can be solved efficiently for many choice of J(u).For example, if we choose the total variation as the regularization term (13) J(u) = |∇u(x)|dx, then ( 12) is nothing but the standard ROF model and can be solved very efficiently via graph-cut method [15].In our method J(u) is the nonlocal total variation, then (12) can be solved as a typical nonlocal denoising problem as discussed in many references, e.g.[13].
Apply this method to (9), we get the following Algorithm 1 to solve (8).

Algorithm 1 The Bregman Iteration
Initialize: Start from some initial guess u.Let f = f .while Au − f 2 not small enough do while Au − end while 3.3.ALGORITHM.We now apply algorithm 1 to our problem.Besides putting Au − and J(u) to be the nonlocal regularization, we also need to notice that in our problem the Φ i is unknown and should be updated as well.
In our implementation we will combine the updating step for the Φ i into the Bregman updating loop, similar as in [17].We will choose the optical flow method described in [3] to update Φ i .The overall algorithm 2 is as follows:

Algorithm 2 The Alternative Optimization Algorithm
Initialize: Start from some initial guess u.Let fi = f i .while i Φ i u − f i 2 not small enough do Estimate Φ i which maps u onto f i from (2) via optical flow scheme.The initial value of u is chosen as the temporal average of the frames.In Fig. 3 an example is shown.We can see that the average of the frames is very blur but gives good initial guess of the rough shape of the object.4.2.NUMBER OF FRAMES.The quantity of frames determines the quality of the reconstruction.On the other hand, the more frames are used, the larger computational cost is consumed due to the fact that the registration between the frames are most time costing.In our numerical experiments we use less than 100 frames and generally we can always get satisfactory results with only 10 frames.In the synthetic example shown in Fig. 8 we only use 20 frames.

COMPUTATION OF Φ i AND Φ ⊤
i .As we described above, the turbulence warping φ i is estimated by optical flow scheme.In our implementation we use the method developed in [26], while other schemes can be applied as well.Given φ i , Φ i u = u(φ i (•)) can be evaluate by interpolation.
Φ ⊤ i is the adjoint operator of Φ i , that is to say, Mathematically, this can be written as ( 15) which is sometimes denoted as ( 16) where φ i #u(x) is called the push-forward action of φ i onto u(x).Given φ i , Φ ⊤ i u(x) can be numerically evaluated as follows.Let (17) v y (x) = 1, y = x 0, y = x be the 'single-pixel spike' function at pixel y, then the value of Φ ⊤ i u at pixel y can be calculated as (18) ( This computation is very efficient because Φ i v y is a very simple function and can be directly evaluated from φ i .
4.4.REGULARIZATION.As discussed in section 3 the regularization step is a well-studied optimization problem and can be solved by many existing routine.We use the method developed in [13] to implement this step.

PARAMETERS.
There are two parameters in our algorithm 2. δ is the step size for the gradient descent of the fidelity term.As shows in [10], the step size should be chosen such that u → u−δ i Φ ⊤ i (Φ i u− f ) is contractive to guarantee the convergence.The parameter λ is not a crucial factor and numerical results also indicate that it is not sensitive.In practice we suggest a small initial λ to make the image regularized enough at the very beginning, then increase the λ gradually.This method is also used in other image reconstruction method, e.g.[1].

EXPERIMENTATIONS
5.1.TESTS RESULTS.In Fig. 4, 5, 6 and 7 some frames of test sequences from real videos are shown, as well as the magnified details.Our reconstruction results are shown in the last row for Fig. 4, 5, 6 and the last column for Fig. 7.Only 5 iterations are implemented in our algorithm to get this satisfactory result (especially in Fig. 7 where the letters on the board are much readable on processed image than in the original frames).Fig. 8 shows a synthetic example, where the wave filter with unknown parameter is imposed on the static image.Only 20 frames are used in this example.

INFLUENCE OF THE OPTICAL FLOW SCHEME CHOICE.
As mentioned before, the optical flow scheme used was one developped in [26].This algorithm gives good optical flow estimation but is time consuming.In Fig. 9, we present some results both obtained by the above mentioned algorithm and the classical Lucas-Kanade [4] optical flow algorithm which is faster.The different experiments seems to show that the results are very similar and then our restoration algorithm is not sensitive to the choice of the optical flow scheme.5.3.THE CASE OF NONSTATIC SCENE.In Fig. 10, we show some results when the proposed algorithm is applied on a sequence where a pedestrian is moving.A restored frame N is processed with only the ten initial frames {N − 10, N }.The number of frames used by the algorithm may clearly depends on the velocity difference between the movement of the pedestrian and movements dues to turbulence.The results show an improved sequence where the turbulence deformations are mostly compensated without altering the pedestrian movement which will be easier to detect.5.4.COMPARISON WITH EXISTING METHODS.Fig. 11 shows a comparison between results obtained by our method and two state of the art other methods: the algorithm based on PCA [20] and the algorithm using the Lucky-Region Fusion [2].In these tests, we deal with short exposure sequences; then it is not surprising that the PCA method failed as in this case the assumption of a gaussian kernel for the blur is not relevant.The Lucky-Region Fusion approach gives some good results but our method shows a better geometry reconstruction (see for example near the high spatial frequencies in the second row, the separation between each bar is more clear in our results).The other advantage of our algorithm is that it need few frames to get a good result, compared to the 100 frames used in the Lucky-Region Fusion.

CONCLUSION
We proposed a novel approach to restore an image from a video sequence contaminated by geometric deformations like the ones observed in atmospheric turbulence.Our method is based on a variational model solved by Bregman iteration and operator splitting method.Both state-of-the-art and classical optical flow methods were utilized to estimate the warping effect.No restriction on the frame rate or sequence order is needed in our method.The algorithm is simple, concise and computationally efficient.
In order to deal with the complete turbulence problem and following the model of Frakes [11], we try some deblurring at the end of the process.We test both nonblind deconvolution (by assuming a gaussian blur kernel, where its size is choosen experimentally) based on total variation regularization [16] or framelet decomposition [6,7], and bayesian blind deconvolution (the deconv function available in Matlab).
The results, shown in Fig. 12, obtained by adding these deblurring effects doesn't give real good improvements, for the two nonblind technics the differences between the inputs and outputs are rather small.Another point, is that different experiments seems to show that the gaussian assumption for the kernel is not relevant for this kind of blur.The Bayesian blind algorithm gives images with sharper edges, but as previously, the improvements are not really significants.Some investigations are currently done to find the best way to deblur.The main questions are: which kind of deblurring method?And when is the best moment to apply deblurring (at the end of the process of during the Bregman iteration)?Finally, we look forward to exploring the broader capabilities of this method in different applications like underwater imaging.

FIGURE 1 .
FIGURE 1. Sample images.Each row contains three arbitrary frames from different testing turbulence videos.

FIGURE 3 .
FIGURE 3. The average of the frames of one example video and its magnification of the top right corner.

FIGURE 4 .
FIGURE 4. The first three rows are example frames and the magnifications of the right part of the frames.The last row shows our reconstructed result.

FIGURE 5 .FIGURE 6 .
FIGURE 5.The first three rows are example frames and the magnification of the top right part of the frames.The last row shows our reconstructed result.

FIGURE 7 .
FIGURE 7. The first two figures are example frames.The last figure is our reconstructed result.Only 10 frames are used in this example.

FIGURE 8 .
FIGURE 8.The first two figures are the distorted Lena.The last figure is our reconstructed result.Only 20 frames are used in this example.

FIGURE 9 .
FIGURE 9. Results obtained by using different optical flow schemes.Black-Anandan scheme on the left column and Lukas-Kanade scheme on the right column.

FIGURE 10 .
FIGURE 10. Results obtained on a sequence with a moving pedestrian (original frames on top and restored frames on bottom).

FIGURE 11 .
FIGURE 11. Results obtained with different algorithms: our method on left, the PCA based approach in the middle and the Lucky-Region Fusion on the right.

FIGURE 12 .
FIGURE 12. Outputs of three different deblurring algorithms: Nonblind TV-deblurring on left, a Nonblind Framelet based approach in the middle and the Matlab Bayesian Blind deconvolution on right.