# Distortion-driven Turbulence Effect Removal using Variational Model

Yuan Xie, Wensheng Zhang, Dacheng Tao, , Wenrui Hu, Yanyun Qu, and Hanzi Wang,
Y. Xie, W. Zhang and W. Hu are with the State Key Lab. of Complex Systems and
Intelligence Science, Institute of Automation, Chinese Academy of Sciences,
Beijing, 100190, China; E-mail: {yuan.xie, wensheng.zhang<EMAIL_ADDRESS>D. Tao is with the Center for Quantum Computation & Intelligent Systems and
the Faculty of Engineering & Information Technology, University of Technology,
Sydney, Australia; E-mail<EMAIL_ADDRESS>Y. Qu and H. Wang are with
School of Information Science and Technology, Xiamen University, Fujian,
361005, China; E-mail: {yyqu<EMAIL_ADDRESS>
###### Abstract

It remains a challenge to simultaneously remove geometric distortion and
space-time-varying blur in frames captured through a turbulent atmospheric
medium. To solve, or at least reduce these effects, we propose a new scheme to
recover a latent image from observed frames by integrating a new variational
model and distortion-driven spatial-temporal kernel regression. The proposed
scheme first constructs a high-quality reference image from the observed
frames using low-rank decomposition. Then, to generate an improved registered
sequence, the reference image is iteratively optimized using a variational
model containing a new spatial-temporal regularization. The proposed fast
algorithm efficiently solves this model without the use of partial
differential equations (PDEs). Next, to reduce blur variation, distortion-
driven spatial-temporal kernel regression is carried out to fuse the
registered sequence into one image by introducing the concept of the near-
stationary patch. Applying a blind deconvolution algorithm to the fused image
produces the final output. Extensive experimental testing shows, both
qualitatively and quantitatively, that the proposed method can effectively
alleviate distortion and blur and recover details of the original scene
compared to state-of-the-art methods.

###### Index Terms:

Image restoration, atmospheric turbulence, variational model, distortion-
driven kernel

##  I Introduction

Figure 1: Illustration of the visual effects caused by atmospheric turbulence.
$d_{1}$ to $d_{n}$ are local deformation fields, which means that a distortion
effect exists. Blur kernels $k_{1}$ to $k_{n}$ and $k_{n}^{\prime}$ indicate
the space-time-varying blur effect.

Atmospheric turbulence can severely degrade the quality of images produced by
long range observation systems, rendering the images unsuitable for vision
applications such as surveillance or scene inspection. The main visual effects
caused by atmospheric turbulence are geometric distortion and space-time-
varying blur (see examples in Fig. 1). The distortion is primarily generated
by (1) optical turbulence and (2) scattering and absorption by particulates;
aerosols, for example, diffuse light will also cause blur [1]. Several
approaches have been used to restore images, including adaptive optics
techniques (e.g. [2, 3]) and pure image processing-based methods (such as [4,
5, 6, 7, 8]). Due to random fluctuations in turbulence, calculating a
reasonable estimation of atmospheric modulation transfer function (MTF) is
extremely difficult. However, this function is of critical importance for
optics-based restoration methods. Therefore, in this article, we only focus on
using image processing to handle the degradation caused by turbulence.
Supposing that the scene and the image sensor are both static, we adopt the
mathematical model used in [7, 9, 10] to interpret imaging processing through
the turbulence:

$f_{i}(\mathbf{x})=D_{i,blur(\mathbf{x})}\Big{(}H(u(\mathbf{x}))\Big{)}+\varepsilon_{i},\quad\forall
i,\quad i\in[1,\ldots,N]$ (1)

where $u$ is the static original scene needed to be retrieved, $f_{i}$ is the
observed frame at time $i$, $N$ represents the number of observed frames, the
vector $\mathbf{x}=(x,y)^{T}$ is a $2$-D spatial location, and
$\varepsilon_{i}$ denotes the sensor noise. $H$ is a blurring operator, which
is caused by sensor optics, and corresponds to a space-invariant diffraction-
limited point spread function (PSF) $h$. $D$ is a deformation operator, and
its subscript indicates that both the local deformation and the space-varying
blur (the PSF is characterized by $h_{i,x}$) exist concurrently.

Simultaneously removing space-time-varying blur and distortion is a nontrivial
problem. Li et al. [11] explicitly formulated multichannel image deconvolution
as a principal component analysis (PCA) problem, and applied it to the
restoration of atmospheric turbulence-degraded images. However, this spectral
method does not fully correct the deformation. Moreover, due to the fact that
high-frequency information is discarded, the local texture of the true scene
is also poorly recovered. Hirsch et al. [9] utilized space-varying blind
deconvolution to alleviate turbulence distortion. While reasonably effective,
it does not take sensor noise into account when estimating the local PSF,
which results in deblurring artifacts.

Image selection and fusion methodologies have been deployed to produce a high-
quality latent image. The “lucky frame” methods [12, 13, 14, 15] select the
relatively high quality frame from a degraded sequence by using sharpness as
the image quality measurement. However, the so-called “lucky frame” is
unlikely to exist in a short exposure video stream. To alleviate this problem,
Aubailly et al. [16] proposed a local version of the “lucky frame” method
referred to as the “lucky region” method. In this method, a high-quality image
from the output is fused with many small lucky regions detected using a local
sharpness metric. However, even though the space-time-varying blur can be
removed during the fusion process, the final output is still susceptible to
the blur caused by the diffraction-limited PSF [7].

The success of some recently proposed turbulence removal methods stems from
the use of diffeomorphic warping and image sharpening techniques [8, 17, 18,
19]. Shimizu et al. [8] applied a temporal median filter to build a reference
image, and then fixed the geometric distortion using B-spline non-rigid
registration associated with an additional stabilization term. A high-
resolution latent image is then obtained by employing a super-resolution
method. Mao and Gilles [19] combined optical flow-based geometric correction
with a non-local total variance-based regularization process to recover the
original scene. However, this method only focuses on correcting distortion and
ignores the rich information in multiple frames when recovering the detail of
the image.

Zhu et al. [7] proposed the diffeomorphic warping and image sharpening
approach by combining a symmetric constraint-based B-spline registration
associated with near-diffraction-limited image reconstruction to reduce the
space and time-varying problem to a shift invariant one. However, this method
has two main limitations: (1) the method uses the temporal mean of the
observed frame to calculate the reference image, which leads to a poor
registration result (especially in the case of strong turbulence); and (2) the
method constructs a single image from the near-diffraction-limited detection
based on the assumption that the distortion can be effectively removed,
meaning that the noise introduced by registration error cannot be reduced.

In this paper, we propose a new scheme for restoring a single sharp image of
the original scene from a sequence of observed frames degraded by atmospheric
turbulence. The proposed method can effectively remove both local deformation
and spatially-varying blur and recover the image details, even in the case of
strong turbulence. The scheme consists of the following four steps. First, we
construct a reference image by _low-rank decomposition_ , resulting in a
sharper and less noisy image than the traditional method. Second, the
reference image is iteratively optimized and enhanced using a variational
model involving a new _spatial-temporal regularization_ , which helps generate
an improved registered sequence. Moreover, we design _a simple but highly
efficient algorithm_ without the use of PDEs in order to solve the variational
model, which is supported by a rigorous proof of convergence. Third, by
introducing the concept of the near-stationary patch, a single image with
reduced space-varying blur is produced using a fusion process. In this step,
to avoid noise effects, _distortion-driven spatial-temporal kernel regression_
is employed to eliminate the noise both in the image and the temporal domains.
To the best of our knowledge, using the information from the deformation field
to guide the construction of the local kernel has not previously been
described. Finally, a blind space-invariant deconvolution algorithm can be
used to generate the final output.

Figure 2: Block diagram of the proposed restoration scheme. (1) Step 1:
construct a high-quality reference image from observed sequence using low-rank
decomposition; this reference image will be used as the initial value of the
variational model in the next step. (2) Step 2: enhance the reference image
iteratively using a new variational model (as the green arrow shown). After
the convergence of optimization, the observed sequence is registered to the
enhanced reference image to produce a registered sequence with little local
deformation (as the dashed arrow line indicated). (3) Step 3: a distortion-
driven fusion process is employed to fuse the registered sequence to one image
with only space-invariant blur effect existing. (4) Step 4: final restored
image is produced by a commonly space-invariant deconvolution method.

This paper is structured as follows. Section II provides an overview of the
proposed restoration algorithm. In Section III, we detail the method used to
optimize the reference image. Image fusion from the registered sequence is
presented in Section IV, while single image deconvolution is presented in
Section V. Experimental results are presented in Section VI and we discuss the
methods and conclude in Section VII.

##  II Restoration Scheme Overview

The proposed restoration framework has four main steps (see Fig. 2). Given an
observed sequence $\\{f_{i}\\}$, step $1$ applies the low-rank decomposition
to $\\{f_{i}\\}$ to generate a high-quality reference image. In multi-frame
registration, the reference image is usually obtained by computing the
temporal mean of the observed frames [6, 7, 8, 17, 18, 19, 20]. An averaged
image such as this is always blurry and noisy, especially when turbulence is
strong. Inspired by [21], which utilizes sparse decomposition for background
modeling, we use matrix decomposition to obtain the low-rank part and
construct the reference image for registration. Our rationale for using this
approach is that the low-rank part is a stable component of a ”dancing image”,
and therefore corresponds to the original scene to some extent. The
decomposition can be defined as follows:

$\text{minimize }||L||_{\ast}+\lambda||S||_{1}\qquad\text{subject to }L+S=G$
(2)

where $G\in R^{m\times n}$ is the distorted sequence matrix with each column
being a distorted frame vector $f_{i}$, $m$ denotes the total pixels in each
frame, and $n$ denotes the number of frames in the sequence. $L\in R^{m\times
n}$ is the low-rank component of $G$, $S\in R^{m\times n}$ is the sparse
component of $G$, $||\cdot||_{\ast}$ is the nuclear norm defined by the sum of
all singular values, $||\cdot||_{1}$ is the $l_{1}$-norm defined by the
component-wise sum of absolute values of all the entries, and $\lambda$ is a
constant providing a trade-off between the sparse and low-rank components.
Some decomposed results are illustrated in Fig. 3. We exploit robust principal
component analysis (RPCA) to estimate the reference image, which acts as the
initial value for the next step.

Step $2$ enhances the reference image using a variational model employing
spatial-temporal regularization. This optimization is the inverse process of
the following problem

$f_{i}=\Phi_{i}u+\varepsilon,\quad\forall i\in[1,\ldots,N]$ (3)

where $\Phi_{i}$ is the linear operator corresponding to the deformation field
$\overleftarrow{F_{i}}$ obtained by B-spline registration [22, 23] and $u$
represents the reference image to be enhanced. The variational model optimizes
with only a few rounds of non-rigid registration and subproblem optimization.
After the convergence of the optimization procedure, the observed frames
$\\{f_{i}\\}$ are registered to the enhanced reference image to achieve the
registered sequence $\\{R_{i}\\}$. This step effectively removes local
deformation.

Step $3$ restores a single image $Z$ from the registered sequence
$\\{R_{i}\\}$ by distortion-driven fusion in order to eliminate the space-
time-varying blur. A near-stationary patch can be detected for each local
region from the patch sequence through the temporal domain. Spatial-temporal
kernel regression is then carried out to reduce the noise caused by optics and
registration error. Fusing all the denoised near-stationary patches generates
an image $Z$. While the output $Z$ is still a blurred image, it can be
approximately restored using a common spatially invariant deblurring
technique.

In the final step, based on the statistical prior of natural images, a single
image blind deconvolution algorithm is implemented on $Z$ to further remove
blur and enhance image quality. Details of steps $3$, $4$ and $5$ will be
given in the following sections.

Figure 3: Low-rank decomposition of the distorted frames in city_strong and
chimney sequences

##  III Optimize the Reference Image Using a Variational Model

The obtained reference image from RPCA is not completely ready for non-rigid
registration because the inner structure of the reference image cannot be
easily recovered. In other words, the edges representing the object’s profile
tend to be well defined, while the edges that characterize the inner structure
of the object tend to be intermittent and noisy. Therefore, the reference
image needs to be refined. Assuming that the observed image sequence is
$\\{f_{i}\\}_{i=1}^{N}$ and the reference image that we want to optimize is
$u$, then directly reconstructing $u$ from Eq. (3) is an ill-posed problem.
There is a need to introduce a regularization technique. If we denote the
regularization term of the image as $J(u)$, then we can get an unconstrained
optimization problem

$\displaystyle\min_{u}E(u)+J(u)$ (4)
$\displaystyle=\min_{u}\underbrace{\sum\nolimits_{i}\|\Phi_{i}u-f_{i}\|_{2}^{2}}_{E(u)}+\underbrace{\Big{(}\mu_{1}J_{s}(u)+\mu_{2}J_{t}(u)\Big{)}}_{J(u)}$
(5)
$\displaystyle=\min_{u}\underbrace{\sum\nolimits_{i}\|\Phi_{i}u-f_{i}\|_{2}^{2}}_{E(u)}+\underbrace{\Big{(}\mu_{1}\overbrace{|u|_{NLTV}}^{J_{s}(u)}+\mu_{2}\overbrace{|(u-u_{p})|_{TV}}^{J_{t}(u)}\Big{)}}_{J(u)}$
(6)

where $E(u)$ measures the fidelity of the observed data, and $J(u)$ includes
the spatial regularizer $J_{s}(u)$ and the temporal regularizer $J_{t}(u)$ to
constraint the solution of $u$. $J_{s}(u)$ is the non-local total variation
which can explore repetitive structures to preserve important details of an
image while effectively removing artifact. $J_{t}(u)$ is the standard total
variation of the difference between two consecutively optimized results
($u_{p}$ denotes the result obtained by the previous step), which does not
only restrict the smoothness between the iterations but also forces the local
energy of the optimized result to converge to a local structure (e.g. the edge
of the object). The parameters $\mu_{1}$ and $\mu_{2}$ are chosen as a trade-
off between the two regularizers. Because $J(u)$ applies the sparsifying
transform in both the spatial and temporal domains, we refer to it as the
spatial-temporal regularization. In the following subsections, we develop a
fast algorithm to solve the optimization problem (6).

###  III-A Optimization by Bregman Iteration

Bregman iteration is a concept that originated from functional analysis and is
commonly used to find the extrema of convex functions [24]. It was introduced
in [25] for total variation-based image processing, before being extended to
other applications such as wavelet-based denoising [26] and magnetic resonance
imaging [27]. In this subsection, we briefly describe how to restore the
optimization problem (6) using this technique.

The Bregman distance [24] associated with a convex function $J$ between points
$u$ and $v$ is

$D_{J}^{p}(u,v)=J(u)-J(v)-\langle p,u-v\rangle\geq 0,\text{where }p\in\partial
J(v)$ (7)

where $p$ is the subgradient of $J$ at $v$. Because $D_{J}^{p}(u,v)\neq
D_{J}^{p}(v,u)$, $D_{J}^{p}(u,v)$ is not a distance in the usual sense.
However, it measures the closeness between $u$ and $v$ in the sense that
$D_{J}^{p}(u,v)>0$, and $D_{J}^{p}(u,v)\geq D_{J}^{p}(w,v)$ for any point $w$
being a convex combination of $u$ and $v$. According to [25], and using the
Bregman distance, the optimization problem (6) can be solved by the following
iteration

$\displaystyle u^{k}=\operatorname*{argmin}_{u}\mu
J(u)+\sum_{i=1}^{N}\|\Phi_{i}u-f_{i}^{k}\|^{2}$ (8a) $\displaystyle
f_{i}^{k+1}=f_{i}^{k}+f_{i}-\Phi_{i}u^{k}$ (8b)

where $f_{i}^{0}=f_{i}$. It has been proved in [25] that the above iteration
converges to the solution of (6). Moreover, as shown in [28, 29], Bregman
iteration is actually equivalent to alternatively decreasing the primal
variable and increasing the dual variable of the Lagrangian of the problem
(6).

In subproblem (8a), the data-fidelity term $E(u)$ co-occurs with the
regularization term $J(u)$, which leads to a sophisticated solution. However,
this subproblem can be solved efficiently using the forward-backward operator
splitting method [30, 31, 32]:

$\displaystyle
v^{k+1}=u^{k}-\delta\sum_{i=1}^{N}\Phi_{i}^{T}(\Phi_{i}u^{k}-f_{i}^{k})$ (9a)
$\displaystyle u^{k+1}=\operatorname*{argmin}_{u}\bigg{(}\mu
J(u)+\frac{1}{2\delta}\|u-v^{k+1}\|^{2}\bigg{)}$ (9b)

where $\Phi_{i}^{T}$ denotes the adjoint operator of $\Phi_{i}$. The parameter
$\delta$ is set to $1$ in experiments. The advantage of the above method is
the separation of $\Phi_{i}$ and $J(u)$. The first step is called the forward
step, which is actually the gradient descent of the data-fidelity term $E(u)$,
and the second step is called the backward step, which can be solved
efficiently for many choices of regularizer $J(u)$. For example, if we choose
the total variation as the regularizer, then the subproblem is actually a
standard Rudin-Osher-Fatemi (ROF) model [33], which can be efficiently solved
via the graph-cut [34] or the split Bregman method [35]. In our method, $J(u)$
is a mixed regularizer that contains non-local total variation (TV) and TV, so
we refer to the subproblem (9b) as the mixed-ROF model.

###  III-B Finding an Efficient Solution for the Mixed-ROF Model

No matter which ROF-like model we choose, current approaches always involve
solving PDEs in each iteration, including the split Bregman method [28, 35].
Here, we extend the method proposed in [36], which only handles the standard
ROF model, to the general mixed ROF case, and provide a very simple but highly
efficient and effective optimization algorithm without using PDEs. By
replacing $J(u)$ with $J_{s}(u)$ and $J_{t}(u)$, the subproblem (9b) can be
expanded as follows:

$\min_{u}\mu_{1}|\nabla_{w}u|_{1}+\mu_{2}\Big{(}|\nabla_{x}(u-u_{p})|_{1}+|\nabla_{y}(u-u_{p})|_{1}\Big{)}+\frac{1}{2}\|u-v\|_{2}^{2}$
(10)

The following new iteration schema can find the unique solution $u^{\ast}$
efficiently for the minimization problem (10).

Let $b_{w}^{0}=0,b_{x}^{0}=0,b_{y}^{0}=0$ and $u^{1}=v$, for $k=1,2,\ldots$,
the iteration is as follows:

$\displaystyle
b_{w}^{k}=cut(\nabla_{w}u^{k}+b_{w}^{k-1},\frac{\mu_{1}}{\lambda_{1}})$ (11)
$\displaystyle
b_{x}^{k}=cut(\nabla_{x}u^{k}-\nabla_{x}u_{p}+b_{x}^{k-1},\frac{\mu_{2}}{\lambda_{2}})$
(12) $\displaystyle
b_{y}^{k}=cut(\nabla_{y}u^{k}-\nabla_{y}u_{p}+b_{y}^{k-1},\frac{\mu_{2}}{\lambda_{2}})$
(13) $\displaystyle
u^{k+1}=v+\frac{\lambda_{1}}{\mu_{1}}div_{w}b_{w}^{k}-\frac{\lambda_{2}}{\mu_{2}}(\nabla_{x}^{T}b_{x}^{k}+\nabla_{y}^{T}b_{y}^{k})$
(14)

where $shrink(x,\gamma)=x-\frac{x}{|x|}\ast\max(|x|-\gamma,0)$.

In the above equations, $\nabla_{x}$ and $\nabla_{y}$ are the difference
operators along the $x$ and $y$ directions, respectively, $\nabla_{w}$ denotes
the non-local gradient operator, and $\text{div}_{w}$ is the non-local graph
divergence operator (see the detailed mathematical definitions of all these
operators in the Appendix). We prove the convergence of the proposed fast
algorithm in the following theorem:

######  Theorem 1.

For $k=1,2,\ldots$, let $b_{w}^{k},b_{x}^{k},b_{y}^{k}$ and $u^{k+1}$ be given
by the iteration (11) to (14). If $0<20\lambda_{1}+4\lambda_{2}<1$, then
$\lim_{k\rightarrow\infty}u^{k}=u^{\ast}$.

###### Proof.

See the proof in Section VIII-B. ∎

##  IV Distortion-Driven Image Fusion

Figure 4: Near-stationary patch vs. near-diffraction-limited patch in the
airport sequence. (a) Ground truth; (b) Zoomed ground truth patch from (a);
(c) Frame contains near-stationary patch; (d) Zoomed near-stationary patch
from (c); (e) Frame contains near-diffraction-limited patch; (f) Zoomed near-
diffraction-limited patch from (e); (g) Movement energy of the local patch
over 60 frames; (h) Intensity variance of local patch over 60 frames.

1Given a registered sequence $\\{R_{k}\\}$ and the push-forward motion field
sequence $\\{\overleftarrow{F_{k}}\\}$ by inversing
$\\{\overrightarrow{F_{k}}\\}$, divide each frame into $L\times L$ overlapping
patches centered at each pixel, and calculate the intensity variance and the
movement energy of each patch as a local stationary measure.

2

3For a patch sequence $\\{r_{k}\\}$ centered at location $\mathbf{x}$, detect
the most stationary one $r_{k^{*}}$ by local variance and movement energy.

4

5Set $r_{k^{*}}$ as a reference patch, and restore its center pixel value
using spatial-temporal kernel regression. Assign this value to the pixel
$Z[\mathbf{x}]$.

6

Go to the next pixel and return to step 2.

Algorithm 1 Procedure for reconstructing a stationary image from registered
frames

In this section, by introducing the concept of the near-stationary patch, we
fuse the registered sequence $\\{R_{k}\\}$ into a single image $Z$, which can
be deblurred using a space-invariant deconvolution method. The fusion steps
are summarized in Algorithm 1. Different from our approach, [7] builds the
image $Z$ by fusing the diffraction-limited patches, which can be detected
using local sharpness measures. However, due to the registration error, some
detected diffraction-limited patches still contain local distortion. Examples
are shown in Fig. 4 where the patches (d) and (f) are obtained by near-
stationary detection and near-diffraction-limited detection, respectively.
Compared with (f), the near-stationary patch (d) shows less local deformation
and more closely resembles the latent true patch (b). Therefore, employing a
local sharpness measure alone is insufficient to recover the image details. In
the following subsections, we will first describe how to detect near-
stationary patches, and then detail the proposed distortion-driven spatial-
temporal kernel regression used for reducing the noise caused by sensor and
registration errors during the fusion process.

###  IV-A Near-Stationary Patch Detection

In turbulence free or near turbulence free conditions, the light reflected
from the scene through the atmosphere is perpendicular to the camera.
Therefore, any scene inside local regions should be clearly observed by the
camera without distortion or blur. Such local regions are referred to as
stationary patches. Suppose the $k^{\ast}$th patch is found to be near-
stationary: we propose to detect this near-stationary patch by taking both
local sharpness and local movement energy into account. According to [7], the
local sharpness of a $L\times L$ patch is determined by its variance in
intensity, which is

$s_{\mathbf{x}}^{k}=\frac{1}{L^{2}-1}\sum_{\mathbf{x}}(r_{k}[\mathbf{x}]-\bar{r}_{k})^{2}$
(15)

where $\bar{r}_{k}$ represents the mean value of patch $r_{k}$. For further
details on the $s_{\mathbf{x}}^{k}$ calculation, please refer to [7].

The local movement energy of an $L\times L$ patch $r_{k}$ can be calculated
as:

$e_{\mathbf{x}}^{k}=\sum_{\mathbf{x}\in\Omega(r_{k})}\mid\overleftarrow{f_{k}}[\mathbf{x}]\mid^{2}$
(16)

where $\overleftarrow{f_{k}}[\mathbf{x}]$ denotes the deformation vector at
position $\mathbf{x}$ from the $k$-th deformation field in sequence
$\\{\overleftarrow{F_{i}}\\}$ (the deformation field $\overrightarrow{F_{i}}$,
which warps the distorted frame $f_{i}$ to the reference image, and its
inverse field $\overleftarrow{F_{i}}$ are calculated from B-spline
registration), and $\Omega(r_{k})$ is the support region of the local patch
$r_{k}$. Consequently, given the stationary measurement, we select the energy
patch $r_{k^{\ast}}$ with the lowest movement as the near-stationary patch
from the ten sharpest ones in the patch sequence $\\{r_{k}\\}$. Then,
$r_{k^{\ast}}$ is used as a reference patch to restore its center pixel value;
this is described in the next subsection.

###  IV-B Fusion by Spatial-Temporal Kernel Regression

This subsection describes how one single image $Z$ is generated by fusing the
detected near-stationary patches. To avoid possible artifacts that may appear
during subsequent deblurring, noise in the selected near-stationary patches
needs to be suppressed. Zhu and Milanfar [7] employed zero-order kernel
regression to estimate the value of a pixel at $\mathbf{x}$ in its reference
frame $\delta$ by

$\widehat{q}_{\delta}[\mathbf{x}]=\frac{\sum_{k}U(\mathbf{x};k,\delta)r_{k}[\mathbf{x}]}{\sum_{k}U(\mathbf{x};k,\delta)}$
(17)

where $U(\cdot)$ denotes a patch-wise photometric distance which can be
calculated by a Gaussian kernel function

$U(\mathbf{x};k,\delta)=\exp\bigg{(}\frac{2\sigma_{n}^{2}}{\mu^{2}}-\frac{\parallel
r_{k}-r_{\delta}\parallel^{2}}{L^{2}\mu^{2}}\bigg{)},$ (18)

where $\sigma_{n}^{2}$ is the noise variance (we set $\sigma_{n}^{2}=2$ in our
experiments), $L^{2}$ denotes the total number of pixels in the patch, and the
scalar $\mu$ is the smoothing parameter [37].

However, this denoising method only takes temporal information into account
and ignores the registration error in the image domain. From Fig. 4 we can see
that the more the movement energy of a patch, the less likely it is to have
its deformation thoroughly removed. Therefore, $r_{k}[\mathbf{x}]$ in Eq. (17)
may not represent the true value of the pixel $\mathbf{x}$. Hence, to
reconstruct a single image $Z$, a spatial-temporal kernel regression is used
to restore each pixel. In the spatial domain (image domain), distortion-driven
asymmetric steering kernel regression is used to remove the deviation of each
central pixel of the patch of $r_{k}$. In the temporal domain, patch-wise
temporal regression is used to correct the central pixel and further reduce
the noise level [7]. The details of the spatial-temporal kernel regression are
described in the following two subsections.

####  IV-B1 Point-Wise Spatial Asymmetric Steering Kernel Regression

Figure 5: Illustration of distortion-driven asymmetric steering kernel.

Directly finding a true pixel value is a nontrivial problem. We resort to
calculating $\widetilde{r}_{k}[\mathbf{x}]$, which denotes the corrected value
of a pixel at position $\mathbf{x}$ in the registered frame $R_{k}$ to replace
$r_{k}[\mathbf{x}]$ in Eq. (17). Then, $\widetilde{r}_{k}[\mathbf{x}]$ can be
estimated using the zero-order kernel regression (Nadaraya-Watson estimator
[38])

$\widetilde{r}_{k}[\mathbf{x}]=\frac{\sum\limits_{x_{j}\in
r_{k}}K_{dskr}(\mathbf{x}_{j}-\mathbf{x})r_{k}[\mathbf{x}]}{\sum\limits_{x_{j}\in
r_{k}}K_{dskr}(\mathbf{x}_{j}-\mathbf{x})}$ (19)

where $K_{dskr}(\cdot)$ represents a distortion-driven kernel function that
depends on the local deformation field. Such a kernel can indicate a pixel’s
direction of deviation; in other words, the weights should be larger in the
deviated direction (see Fig. 5 $s1$ to $s5$). Inspired by the steering kernel
[39], we propose a two-step approach to construct a distortion-driven kernel
function.

First, measure the dominant orientation of the local deformation field and
build the symmetric kernel function. The orientation is the singular vector
corresponding to the largest (nonzero) singular value of the local deformation
matrix

$M_{k}=\begin{bmatrix}\vdots&\vdots\\\
\overleftarrow{f_{x}}(\mathbf{x}_{j})&\overleftarrow{f_{y}}(\mathbf{x}_{j})\\\
\vdots&\vdots\end{bmatrix}=U_{k}S_{k}V_{k}^{T},~{}~{}~{}\mathbf{x}_{j}\in\Omega(r_{k}),$
(20)

where $\overleftarrow{f_{x}}$ and $\overleftarrow{f_{y}}$ are the displacement
vectors along the $x$ and $y$ direction, and $\Omega(r_{k})$ is the local
support region of $r_{k}$. $U_{k}S_{k}V_{k}^{T}$ is the truncated singular
value decomposition of $M_{k}$, and $S_{k}$ is a diagonal $2\times 2$ matrix
representing the energy in the dominant direction. Then, in the first column
of the orthogonal matrix $V_{k}$, $v_{1}=[\nu_{1},\nu_{2}]^{T}$ defines the
dominant orientation angle $\theta_{k}=\arctan(\frac{\nu_{2}}{\nu_{1}})$. This
means that the singular vector corresponding to the largest singular value of
$M_{k}$ represents the dominant orientation of the local deformation field.
The elongation parameter $\sigma_{i}$ and scaling parameter $\gamma_{i}$ are
calculated by

$\sigma_{i}=\frac{s_{1}+\lambda^{{}^{\prime}}}{s_{2}+\lambda^{{}^{\prime}}},\quad\lambda^{{}^{\prime}}\geq
0,\quad\gamma_{i}=\bigg{(}\frac{s_{1}s_{2}+\lambda^{{}^{\prime\prime}}}{M}\bigg{)}^{\frac{1}{2}},$
(21)

where $s_{1}$ and $s_{2}$ are the diagonal elements of the matrix $S_{k}$, and
$\lambda^{{}^{\prime}}$ and $\lambda^{{}^{\prime\prime}}$ denote the
regularization parameters. Given the parameters mentioned above, we can
calculate the steering kernel function as in [39] (hereafter, the symmetric
steering kernel is referred to as $K^{s}_{dskr}(\cdot)$). Fig. 5 is a visual
illustration of the steering kernel footprints for different local deformation
regions. Image (a) is one registered frame from $\\{R_{i}\\}$, and the patches
$d1$ to $d5$ (cropped from deformation field (b)) are the local deformation
fields of the pixels marked by red dots in image (a). As shown in Fig. 5, the
size and shape of each symmetric steering kernel ($s1$ to $s5$) are locally
adapted to the corresponding deformation field. For example, $s2$ is blunt due
to the small movement energy of its corresponding deformation field, while
$s3$ is sharp because the movement energy of $d3$ is much larger than that of
$d2$.

Second, construct the asymmetric steering kernel using the dominant
orientation. State-of-the-art non-rigid registration approaches are usually
unable to correct a local deformation when the movement energy is very large.
Nevertheless, due to the smoothness of the B-spline-based registration, the
true pixel value $\widetilde{r}_{k}[\mathbf{x}]$ often appears in the opposite
direction of the dominant orientation. As shown in Fig. 5, patches (g) and (f)
are cropped from one observed frame (d) and its corresponding registered frame
(a), respectively, and patch (e) is cropped from the ground truth image (c)
(the central pixels of (e), (f), and (g) are marked by red dots). Comparing
patch (e) with patch (f), the two central pixels represent different positions
due to the registration error. In fact, the central pixel of patch (e)
corresponds to the pixel marked in green in patch (f). Consequently, to
estimate $\widetilde{r}_{k}[\mathbf{x}]$ more precisely, we need to further
shrink the footprint of the local kernel; the weights along the inverse
dominant orientation should be larger than those along the dominant
orientation. Hence, we construct a local asymmetric kernel function
$K^{a}_{dskr}(\cdot)$ by utilizing the Asymmetric Gaussian proposed in [40],
where the dimension needed to be asymmetric can be calculated by

$\displaystyle\textit{A}(z_{i};\mu^{z}_{i},\rho^{2}_{i},r_{i})=\frac{2}{\sqrt{2\pi}}\frac{1}{\sqrt{\rho^{2}_{i}}(r_{i}+1)}\left\\{\begin{aligned}
&exp\Big{(}-\frac{(z_{i}-\mu^{z}_{i})^{2}}{2\rho^{2}_{i}}\Big{)}\qquad\text{if}\quad
z_{i}>\mu^{z}_{i}\\\
&exp\Big{(}-\frac{(z_{i}-\mu^{z}_{i})^{2}}{2r^{2}_{i}\rho^{2}_{i}}\Big{)}\qquad\text{otherwise}\\\
\end{aligned}\right.$ (22)

where $z_{i}$ is along the dominant orientation of the symmetric kernel
$K^{s}_{dskr}(\cdot)$, $\mu^{z}_{i}$ and $\rho_{i}^{2}$ are the mean and
variance corresponding to $z_{i}$, and $r_{i}$ denotes the asymmetric
coefficient with $r_{i}=1$ being equivalent to the ordinary Gaussian. In our
implementation, $r_{i}$ is usually set to $0.5\sqrt{\sigma_{i}}$, where
$\sigma_{i}$ is the elongation coefficient defined in Eq. (21). An example of
an asymmetric kernel is shown in Fig. 5 $k1$ to $k5$. Due to the computational
complexity, we only use spatial kernel regression for the pixels whose
movement energy is greater than a pre-defined threshold. For the remaining
pixels with relatively low movement energy, we directly assign
$r_{k}[\mathbf{x}]$ to $\widetilde{r}_{k}[\mathbf{x}]$.

####  IV-B2 Patch-Wise Temporal Kernel Regression

Once the corrected pixel value $\widetilde{r_{k}}[\mathbf{x}]$ in each frame
has been calculated, we can further estimate the pixel value at position
$\mathbf{x}$ by means of spatial-temporal kernel regression for a single image
$Z$

$\displaystyle\widehat{q}_{\delta}[\mathbf{x}]=\frac{\sum_{k}U(\mathbf{x};k,\delta)\widetilde{r_{k}}[\mathbf{x}]}{\sum_{k}U(\mathbf{x};k,\delta)}=\frac{\sum_{k}U(\mathbf{x};k,\delta)\frac{\sum\limits_{x_{j}\in
r_{k}}K^{a}_{dskr}(\mathbf{x}_{j}-\mathbf{x})r_{k}[\mathbf{x}]}{\sum\limits_{x_{j}\in
r_{k}}K^{a}_{dskr}(\mathbf{x}_{j}-\mathbf{x})}}{\sum_{k}U(\mathbf{x};k,\delta)}.$
(23)

Spatial-temporal kernel regression is illustrated in Fig. 6: the yellow box
denotes the spatial regression that exploits the local deformation structure
to restore a certain pixel only in the image domain, and the red line
represents the temporal regression, which fuses the information across all the
corrected pixels $\\{\widetilde{r_{k}}[\mathbf{x}]\\}_{k=1}^{N}$ in the
temporal domain. Due to utilization of both spatial and temporal information,
we refer to such regression as the Spatial-Temporal Kernel Regression.

Figure 6: Illustration of distortion-driven spatial-temporal kernel
regression.

##  V Space-invariant Deconvolution

A post-process is needed to correct the diffraction-limited blur that still
exists in $Z$. We exploit the blind deconvolution algorithm [41] to calculate
a final output (corresponding to step $18$ in Algorithm 2). In the degradation
model

$Z=L\otimes h+\varepsilon$ (24)

the deblurring procedure can be described using the following optimization
problem

$\operatorname*{argmin}_{F,h}\|Z-L\otimes
h\|_{2}^{2}+\gamma_{1}R_{f}(L)+\gamma_{2}R_{h}(h)$ (25)

where $\varepsilon$ denotes the error caused by the fusion process in Section
IV, and $R_{f}$ and $R_{h}$ are the regularizers used to restrict the latent
sharp image $L$ and the blur kernel $h$ based on their own prior knowledge. As
suggested in [41], the regularization term for $L$ is defined as

$R_{f}(L)=\|\rho(\partial_{x}L)+\rho(\partial_{y}L)\|_{1}$ (26)

where $\partial_{x}L$ and $\partial_{y}L$ represent the derivatives of $L$ in
horizontal and vertical directions, respectively. The $\rho(\cdot)$ is defined
as

$\rho(x)\left\\{\begin{aligned} &-\theta_{1}|x|\qquad&x\leq l_{t}\\\
&-(\theta_{2}x^{2}+\theta_{3})\qquad&x>l_{t}\\\ \end{aligned}\right.$ (27)

where $l_{t}$, $\theta_{1}$, $\theta_{2}$, and $\theta_{3}$ are all fixed
parameters. The sparse constraint is also imposed on the PSF $h$:
$R_{h}(h)=\|h\|_{1}$. The detail of the algorithm used for solving the
optimization problem (25) can be found in [41]. In our experiments, two types
of parameter setting are used for the simulated data and the real data,
respectively. In simulated experiments, “kernelWidth”, “kernelHeight”,
“noiseStr”, and “deblurStrength” are set to $5$, $5$, $0.03$, and $0.2$,
respectively; all other parameters use the default settings, and the
description of all the parameters can be found in [41]. For real sequences,
setting the same parameters to $(7,7,0.03,0.5)$ can reproduce the displayed
results. Finally, we summarize the proposed atmospheric turbulence removal
algorithm in Algorithm 2.

Input: Distorted image set $V\in\mathbb{R}^{w\times h\times
n}=\\{f_{1},\ldots,f_{n}\\}$

Output: Restored image $L$

1 $u\leftarrow\text{Low\\_Rank\\_Decomposition}(V)$

2 // out_loop

3 while  _$iter <iter\\_max$ _ do

4 // middle_loop

5 while  _$\sum_{i}\|\Phi_{i}u-f_{i}\| >\varepsilon_{1}$_ do

6 Initialize the $\tilde{f}_{i}=f_{i}$

7 Compute each $\Phi_{i}$ which warps $u$ to the observed image $f_{i}$ using
B-spline registration

8 // inner_loop

9 while  _$\sum_{i}\|\Phi_{i}u-\tilde{f}_{i}\| >\varepsilon_{2}$_ do

10 $v=u-\delta\sum_{i}\Phi_{i}^{T}(\Phi_{i}u-\tilde{f}_{i})$

11 Solve $u=\operatorname*{argmin}\nolimits_{u}\mu
J(u)+\frac{1}{2\delta}\|u-v\|^{2}$ by iteration in Section 3.2

12

13 end while

14 $\tilde{f}_{i}\leftarrow\tilde{f}_{i}+f_{i}-\Phi_{i}u$

15

16 end while

17 $\hat{u}=u$

18 $\\{R_{i}\\}\leftarrow$ register each frame $f_{i}$ to $\hat{u}$

19 $V\leftarrow\\{R_{i}\\}$

20 $u\leftarrow\text{Low\\_Rank\\_Decomposition}(V)$

21

22 end while

23

24Construct single frame $Z$ from $\\{R_{i}\\}$ using Spatial-Temporal kernel
regression in Section 4

25 Deblur the image $Z$ to get the latent sharp image $L$ in Section 5

26 Return $L$

Algorithm 2 Atmospheric Turbulence Removal Algorithm

##  VI Experimental Results and Analysis

This section presents extensive experimental validation of the proposed
restoration method. We first show that low-rank decomposition improves the
quality of the reference image compared to traditional methods. Next, we
compare the results optimized using the proposed variational model to those
using a spatial regularizer alone, to illustrate the advantages of our model.
Finally, both qualitative and quantitative methods are used to evaluate the
performance of the proposed method in comparison with several state-of-the-art
methods. For quantitative evaluation, the Peak Signal to Noise Ratio (PSNR)
and Structural Similarity Index (SSIM) are adopted to objectively evaluate the
quality of the restored images.

For all the experiments, the intervals of the control points in the
registration are set to $\varepsilon_{x}=\varepsilon_{y}=16$ pixels, and the
patch size $L$ of $r_{k}$ is set to $13$. For non-local total variance, the
local patch size is set to $5\times 5$ (the support region for a certain
pixel), the search window size is set to $21\times 21$ (the region for
searching for similar patches), and the number of best neighbors is set to
$10$ (the number of accepted similar pixels in the search window). Moreover,
for the number of optimized iterations (Algorithm 2 in Section V),
$out\\_loop=1$ is sufficient for all the test sequences and most of
degradation cases, and the $middle\\_loop$ and $inner\\_loop$ are set to $3$
and $10$, respectively. Furthermore, in the variational model, the parameters
$\lambda_{1}$ and $\lambda_{2}$ should satisfy the condition defined in Lemma
1 in Section VIII-B, so they are chosen to be equal:
$\lambda_{1}=\lambda_{2}=0.02$. With respect to the trade-off between spatial
and temporal regularization, we choose $\mu_{1}=0.5$ and $\mu_{2}=0.25$. The
proposed method is implemented in Matlab with MEX, and all the experiments are
performed on a standard Intel Core i7 2.8GHz computer. The code and data of
the proposed method are available at
https://sites.google.com/site/yuanxiehomepage/.

We compare the proposed method with five representative algorithms: the lucky
region method [16] (Lucky region), principal components analysis for
atmospheric turbulence [11] (PCA), the data-driven two-stage approach for
image restoration [42] (Twostage), Bregman iteration and non-local total
variance for atmospheric turbulence stabilization [19] (BNLTV), and near-
diffraction-limited-based image reconstruction for removing turbulence [7]
(NDL). It is worth noting that the Twostage method was originally designed for
recovering the true image of an underwater scene from a sequence distorted by
water waves. Even though the medium is different, Twostage still achieves
reasonable results on images distorted by air turbulence, and therefore the
comparison and inclusion of this algorithm is valid. The respective authors
provide the Twostage [42] and NDL [7] code, and the parameters remain
unchanged. We have implemented the code for the other algorithms, which
remains faithful to the original papers.

###  VI-A Quality of the Low-rank-based Reference Image

Figure 7: Comparison of the reference images constructed by three methods

In this subsection we compare the visual quality of the reference images
generated by our proposed method, temporal averaging, and the lucky region
approach [16]. The reference images are presented in Fig. 7, where the first
column (a) contains three observed frames from the ’chimney’, ’airport’, and
’city_strong’ sequences; the other three columns are the results of temporal
averaging (b), lucky-region (c), and low-rank decomposition (d), respectively.
In the chimney sequence, the low-rank decomposed reference image is slightly
sharper than the other two images. Nevertheless, due to the image being
affected by only weak turbulence, the observed difference between the three
reference images is small. However, in the case of severe turbulence (airport
sequence), the effects are more noticeable. The magnified area of the region
of interest (marked by the green box) can be seen in the lower right hand
corner of each image. Low-rank decomposition produces sharp edges, while the
images produced by temporal averaging and the lucky region methods remain
blurred; similar results are seen in the city_strong sequence. In addition,
the other methods produce heavy edge artifacts, e.g., in the region of the
bridge, with much sharper and clearer edges produced by our method, even in
the case of relatively indistinct edges. These indistinct edges can
subsequently be enhanced by using the variational model.

###  VI-B Advantages of the Spatial-Temporal Regularization

Figure 8: Optimization of the two regularizers for the city_strong sequence.
Figure 9: Optimization of the two regularizers for the airport sequence.

This subsection illustrates the advantages of spatial-temporal regularization.
We therefore compare the reference image $u$ optimized using the two
regularization terms: (1) the spatial-temporal regularizer as proposed; and
(2) the spatial regularizer (traditional NLTV). Figs. 8 and 9 show the process
of optimization using the two regularizers on the city_strong and airport
sequences, respectively. Image (a) is the initial reference image obtained by
low-rank decomposition, and the other three columns from left to right present
the first, second, and the final optimization steps (the three steps
correspond to the $middle\\_loop$ in Algorithm 2).

As can be seen by comparing (a) with the bottom row of (d) in Fig. 8, the
spatial regularizer enhances ’strong’ edges (we refer to object’s profile as
the strong edges, and consider the edges which are inside an object to
describe local structures as the weak edges), but the weak edges are not
optimally recovered and remain blurry. In contrast, the spatial-temporal
regularizer sharpens both strong and weak edges (top row Fig. 8). This
phenomenon can be explained by the following two factors: 1) NLTV is limited
by its preservation of texture and local structures in an image due to the
smoothness caused by the weighted averaging of non-local self-similar
patterns. Unfortunately, these local structures usually include weak edges,
resulting in different degrees of recovery of strong and weak edges. 2) The
total variation employed by $J_{t}(u)$ can force the energy difference of two
optimized results to converge on the sparse discontinuities, on which the
edges lie in the functional space. Therefore, $J_{t}(u)$ can enhance both
strong and weak edges. In summary, combining the non-local TV-based spatial
regularizer with the local TV-based temporal regularizer does noticeably
improve the quality of the reference image. More illustrative results can be
seen in Fig. 9.

Additionally, we analyze the computational complexity and CPU time of the
proposed fast algorithm and the split Bregman method (the iterative
optimization algorithm is presented in Section VIII-A (from Eqs. (28) to
(31)). For a fixed number of iterations, both the split Bregman method and our
approach are linear in $N$ (the number of pixels), since each step only
contains addition and scalar multiplication operations. Therefore, we compare
them in a different way by accounting for the number of atom operators in the
key steps of each method, which in this case are Eqs. (14) and (28),
respectively, because the complexity difference between $cut(\cdot)$ and
$shrink(\cdot)$ is a constant. Considering one pixel access as an atom
operator (ao), we can compare the complexity of the two methods in detail.
$\nabla_{x}$ and $\nabla_{x}^{T}$ are both $2$ aos, as are $\nabla_{y}$ and
$\nabla_{x}^{T}$. According to
$\Delta=-\nabla_{x}^{T}\nabla_{x}-\nabla_{y}^{T}\nabla_{y}$, $\Delta$ is $8$
aos. $\nabla_{w}$, $\Delta_{w}$ and $\text{div}_{w}$ are $20$ aos if we choose
ten similar patterns in NLTV. Therefore, solving Eq. (28) requires $86N$ aos,
while Eq. (14) only requires $25$ aos. So, for a large image or a large number
of optimized iterations, the proposed algorithm can significantly reduce the
computational time. Table I compares the split Bregman method and the proposed
fast algorithm for solving the spatial-temporal regularizations in terms of
the CPU seconds. The absence of PDEs leads to at least a one-quarter reduction
in running time, and the proposed fast algorithm is therefore highly
efficient.

Figure 10: Three sharp images used in simulated experiments. TABLE I: Comparison of CPU time (average) for the split Bregman and the proposed fast algorithm for differently sized images | $240\times 240$ | $260\times 260$ | $320\times 240$
---|---|---|---
Split Bregman | 0.6496 | 0.6231 | 1.1588
Proposed | 0.4482 | 0.4693 | 0.7805
Reduction | $31.00\%$ | $24.68\%$ | $32.65\%$

###  VI-C Simulated Experiments

Figure 11: Image restoration results on the simulated city_weak sequence. Figure 12: Image restoration results on the simulated city_strong sequence. TABLE II: Comparison of PSNR and SSIM obtained by different methods Sequence | Average | Lucky-Region | PCA Based | Twostage | BNLTV | NDL | Proposed
---|---|---|---|---|---|---|---
Airport | $69.7747$ | $69.7758$ | $66.7340$ | $68.7511$ | $69.3926$ | $68.3036$ | $\mathbf{69.9754}$
$0.2745$ | $0.2731$ | $0.2042$ | $0.2566$ | $0.3830$ | $0.3520$ | $\mathbf{0.4011}$
City_strong | $68.2606$ | $66.9073$ | $63.2749$ | $66.4798$ | $67.0916$ | $\mathbf{68.6392}$ | $67.4820$
$0.2479$ | $0.2405$ | $0.1329$ | $0.2378$ | $0.3413$ | $0.2854$ | $\mathbf{0.3544}$
City_weak | $68.8989$ | $67.3415$ | $64.8207$ | $67.4399$ | $67.3981$ | $\mathbf{68.9629}$ | $68.7218$
$0.3611$ | $0.3531$ | $0.2215$ | $0.2686$ | $0.3631$ | $0.3956$ | $\mathbf{0.4022}$
Architecture | $70.2634$ | $68.1599$ | $64.1939$ | $68.8656$ | $68.3711$ | $67.9243$ | $\mathbf{70.5919}$
$0.4299$ | $0.4282$ | $0.2929$ | $0.3656$ | $0.4349$ | $0.4340$ | $\mathbf{0.5851}$
Chimney | $78.9262$ | $79.1731$ | $64.4229$ | $78.7460$ | $75.0367$ | $79.5622$ | $\mathbf{79.9611}$
$0.1194$ | $0.1137$ | $0.0539$ | $0.1193$ | $0.1120$ | $0.1075$ | $\mathbf{0.1257}$
Building | $72.1991$ | $72.1987$ | $62.4768$ | $73.2032$ | $72.9388$ | $74.3453$ | $\mathbf{74.3843}$
$0.3201$ | $0.3271$ | $0.1858$ | $0.3882$ | $0.3855$ | $0.4857$ | $\mathbf{0.5023}$
Figure 13: Image restoration results on the simulated airport sequence.

To quantitatively evaluate the performance of the proposed method, we
generated a set of simulated degraded sequences with different degrees of
turbulence. The latent sharp images ($260\times 260$) are shown in Fig. 10.
The simulation algorithm is similar to [7] but with a slight difference. The
algorithm includes three key components: the deformation field, spatially
variant PSFs, and a spatially-invariant diffraction-limited PSF. The
deformation field is determined by a set of control points whose random
offsets have a Gaussian distribution with mean $\mu_{d}=0$ and variance
$\sigma_{d}^{2}$, with $\sigma_{d}^{2}$ indicating the turbulence strength.
The space-varying blur is generated by convolution with a set of PSFs, each of
which is also a Gaussian function with variance proportional to the magnitude
of the corresponding local motion energy. In our implementation, the number of
control points is characterized by their interval $d_{g}$, and the number of
PSFs is equal to the number of control points. To efficiently apply a
spatially-variant blur, we use the fast algorithm described in [43], which is
based on overlap-add convolution schemes [44] and linear interpolation of
measured PSFs for spatially-invariant blurs. The spatially-invariant
diffraction-limited blur is produced using a disc function.

To test the performance of the different restoration methods with different
degrees of turbulence, we produced two degraded sequences for the same true
image (Fig. 10 (a)) named city_weak and city_strong, respectively. Considering
limited page length, we only simulate two kinds of turbulence for one image
(city); the other two images are used to produce the airport sequence with
strong turbulence and the arch sequence with weak turbulence, respectively.
The parameter configuration for the two cases are listed as follows: for weak
turbulence, $\sigma_{d}^{2}=4$, $d_{g}=32$, and the variance of the added
Gaussian noise $\sigma_{n}^{2}=3$; for strong turbulence, $\sigma_{d}^{2}=10$,
$d_{g}=16$, and the variance of the added Gaussian noise $\sigma_{n}^{2}=16$.
Table II compars the PSNR and SSIM values for all the outputs of the seven
different restoration algorithms. Each sequence has two sub rows, the top one
denoting the PSNR values, and the bottom one denoting the SSIM values.

The restoration results of the city_weak sequence are illustrated in Fig. 11.
The lucky region method (Fig. 11 (d)) does not seem to offer an improvement
over the mean image (Fig. 11 (c)). The PCA-based method produces unnatural
components due to the loss of high frequency information. BNLTV restores the
structure of the objects well but cannot recover local details. This is due in
large part to the smoothness caused by the weighted averaging of non-local
self-similar patterns. Consequently, the smoothness prevents the restoration
algorithm from recovering the local details. Twostage and NDL achieve similar
results (Figs. 11 (f) and (h)); the distortion has been corrected thoroughly
but some blur still exists. The proposed approach significantly improves
visual quality and recovers many high-frequency details in the image. Even
more noticeable differences are shown in Fig. 12, depicting the results of the
city_strong sequence with severe turbulence (see one of the observed frames in
Fig. 12 (b)). Twostage and BNLTV are superior to NDL, as shown in Fig. 12 (h),
which produces many artifacts. Only the proposed method can remove large
deformations and simultaneously recover detail. The superiority of the
proposed algorithm is due to the high-quality reference image and its
subsequently optimized version facilitating the removal of distortion, and the
spatial-temporal kernel regression recovering sharp local details as well as
reducing the noise introduced by the sensor and the registration error. More
experimental results are shown in Fig. 13 and Fig. 1 in the Appendix.

###  VI-D Real Video Experiments

Figure 14: Image restoration results on the chimney sequence.

We test several video sequences captured through real atmospheric turbulence
to illustrate the performance of the proposed restoration algorithm. In the
first sequence, a chimney is captured through hot air exhausted by a
building’s vent (the original size of the videos ’chimney’ and ’building’ is
$237\times 237\times 100$; we resize them to $240\times 240\times 100$ in
order to facilitate non-rigid registration). The final outputs of all the
methods are shown in Fig. 14: the proposed algorithm provides the best
restoration result and faithfully recovers details of the object. The PSNR and
SSIM values also indicate that the proposed method outperforms the other
methods for the chimney sequence. A noteworthy phenomenon is that all the PSNR
values of the average outputs are relatively high (shown in Table II), because
PSNR is known to sometimes correlate poorly with human perception [45].

Similar restoration results are seen in the experiment on the second test
sequence, ’building’ ($240\times 240\times 100$). Most restoration methods
cannot remove diffraction-limited blur, except for NDL and our algorithm. Both
of these methods can produce sharp restored images, but the NDL output (Fig.
15 (g)) contains halo artifact near the edges, such as around the windows of
the building. Such artifacts can be attributed to the limited accuracy of the
PSF estimation and the existence of noise caused by sensor and registration
error. As mentioned earlier, NDL and the proposed method utilize the same
deconvolution algorithm with identical parameters to produce the final
outputs. Since the major cause of artifact appears to be noise, we can
conclude that the noise can be effectively reduced by spatial-temporal kernel
regression. More real data experimental results are presented in the Appendix.

Figure 15: Image restoration results on the building sequence.

##  VII Discussion and Conclusions

In this paper, we propose a new method to restore a high-quality image from a
given image sequence degraded by atmospheric turbulence. Geometric distortion
and space-time-varying blur are the major challenges that need to be overcome
during restoration. To effectively remove local distortion, the proposed
method constructs a high-quality reference image using low-rank decomposition.
Then, to further improve the registration, the proposed method applies a
variational model with a novel spatial-temporal regularization term to
iteratively optimize the reference image; the proposed fast algorithm can
efficiently solve this variational model. To reduce blur variation, the
registered frames can be fused into a single image with reduced PSF variation
by applying near-stationary patch detection and distortion-driven spatial-
temporal kernel regression. The fused image can be deblurred by a space-
invariant blind deconvolution method in order to produce the final output.
Thorough empirical studies on a set of simulated data and real videos
demonstrate the effectiveness and the efficiency of the new framework for
recovering a single image from a degraded sequence.

Additionally, we test the performance of the proposed restoration method on
images distorted by water rather than air (see the Appendix in supplemental
material), with promising results. This method is therefore suitable for use
with multiple distorting media, although further adaptations of the proposed
method are likely.

##  VIII Proof

###  VIII-A Split Bregman for Solving Mixed ROF Model

The split Bregman method for solving the mixed ROF subproblem (10) is
presented below:

Let $b_{w}^{0}=0,b_{x}^{0}=0,b_{y}^{0}=0$ and $u^{0}=v$, for $k=1,2,\ldots$,
the iteration is as follows

$\displaystyle u^{k+1}=$
$\displaystyle(I-2\lambda_{1}\Delta_{w}-\lambda_{2}\Delta)^{-1}(v+\lambda_{1}\text{div}_{w}(b_{w}^{k}-d_{w}^{k})+\lambda_{2}\nabla_{x}^{T}(d_{x}^{k}+\nabla_{x}u_{p}-b_{x}^{k})$
(28)
$\displaystyle+\lambda_{2}\nabla_{y}^{T}(d_{y}^{k}+\nabla_{y}u_{p}-b_{y}^{k}))$
$\displaystyle
d_{w}^{k+1}=shrink\Big{(}\nabla_{w}u^{k+1}+b_{w}^{k},\frac{\mu_{1}}{\lambda_{1}}\Big{)}$
(29) $\displaystyle
d_{x}^{k+1}=shrink\Big{(}\nabla_{x}u^{k+1}-\nabla_{x}u_{p}+b_{x}^{k},\frac{\mu_{2}}{\lambda_{2}}\Big{)}$
(30) $\displaystyle
d_{y}^{k+1}=shrink\Big{(}\nabla_{y}u^{k+1}-\nabla_{y}u_{p}+b_{y}^{k},\frac{\mu_{2}}{\lambda_{2}}\Big{)}$
(31)

where $\nabla_{x}$ and $\nabla_{y}$ are the difference operators along the $x$
and $y$ directions, respectively,
$\Delta:=-\nabla_{x}^{T}\nabla_{x}-\nabla_{y}^{T}\nabla_{y}$ is the discrete
Laplacian operator, $\nabla_{w}$ denotes the non-local gradient operator and
$\Delta_{w}$ represents the non-local graph Laplacian, and $\text{div}_{w}$ is
the non-local graph divergence operator (see the mathematical definitions of
all those operators in Appendix). Moreover, the function $shrink(\cdot,\cdot)$
is defined as:

$shrink(x,\gamma)=\frac{x}{|x|}\ast\max(|x|-\gamma,0)$ (32)

###  VIII-B Proof of Theorem 1

The detailed proof of the theorem 1 is given below.

######  Lemma 1.

Given $0<20\lambda_{1}+4\lambda_{2}<1$, the real symmetric linear operator
$I+2\lambda_{1}\Delta_{w}+\lambda_{2}\Delta$ is positive definite.

######  Lemma 2.

For $k=1,2,\ldots$, let
$d_{w}^{k},d_{x}^{k},d_{y}^{k},b_{w}^{k},b_{x}^{k},b_{y}^{k}$ and $u^{k+1}$ be
given by the iteration (11) to (14). Then
$\lim_{k\rightarrow\infty}(u^{k+1}-u^{k})=0$, and

$\displaystyle
b_{w}^{k}=cut(\nabla_{w}u^{k}+b_{w}^{k-1},\frac{\mu_{1}}{\lambda_{1}})$
$\displaystyle
b_{x}^{k}=cut(\nabla_{x}u^{k}-\nabla_{x}u_{p}+b_{x}^{k-1},\frac{\mu_{2}}{\lambda_{2}})$
$\displaystyle
b_{y}^{k}=cut(\nabla_{y}u^{k}-\nabla_{y}u_{p}+b_{y}^{k-1},\frac{\mu_{2}}{\lambda_{2}})$

######  Lemma 3.

For $k=1,2,\ldots$, let
$d_{w}^{k},d_{x}^{k},d_{y}^{k},b_{w}^{k},b_{x}^{k},b_{y}^{k}$ and $u^{k+1}$ be
given by the iteration (11) to (14). Then all the sequences
$(d_{w}^{k})_{k=1,2,\ldots},(d_{x}^{k})_{k=1,2,\ldots},(d_{y}^{k})_{k=1,2,\ldots},(b_{w}^{k})_{k=1,2,\ldots},(b_{x}^{k})_{k=1,2,\ldots},(b_{y}^{k})_{k=1,2,\ldots}$
and $(u^{k})_{k=1,2,\ldots}$ are bounded. Moreover,

$\displaystyle
u^{k+1}=v+\frac{\lambda_{1}}{\mu_{1}}div_{w}b_{w}^{k}-\frac{\lambda_{2}}{\mu_{2}}(\nabla_{x}^{T}b_{x}^{k}+\nabla_{y}^{T}b_{y}^{k})$

######  Lemma 4.

For $k=1,2,\ldots$, let
$d_{w}^{k},d_{x}^{k},d_{y}^{k},b_{w}^{k},b_{x}^{k},b_{y}^{k}$ and $u^{k+1}$ be
given by the iteration (11) to (14). Then

$\displaystyle\lim_{k\rightarrow\infty}(d_{w}^{k}-\nabla_{w}u^{k})=0$
$\displaystyle\lim_{k\rightarrow\infty}(d_{x}^{k}-\nabla_{x}(u^{k}-u_{p}))=0$
$\displaystyle\lim_{k\rightarrow\infty}(d_{y}^{k}-\nabla_{y}(u^{k}-u_{p}))=0$

The proofs of the above lemmas are given in Appendix in supplemental material.

######  Theorem 1.

For $k=1,2,\ldots$, let $b_{w}^{k},b_{x}^{k},b_{y}^{k}$ and $u^{k+1}$ be given
by the iteration (11) to (14). If $0<20\lambda_{1}+4\lambda_{2}<1$, then
$\lim_{k\rightarrow\infty}u^{k}=u^{\ast}$.

###### Proof.

Let $F(u):=(1/2)\|u-v\|_{2}^{2}$, then $\partial F(u)=u-v$. For
$w\in\mathbb{R}^{N^{2}}$ we can get

$F(u^{k+1}+w)-F(u^{k+1})-\langle u^{k+1}-v,w\rangle\geq 0$ (33)

From Lemma 3, we have
$-(u^{k+1}-v)=\frac{\lambda_{2}}{\mu_{2}}\nabla_{x}^{T}b_{x}^{k}+\frac{\lambda_{2}}{\mu_{2}}\nabla_{y}^{T}b_{y}^{k}-\frac{\lambda_{1}}{\mu_{1}}\text{div}_{w}b_{w}^{k}$,
and moreover $\langle\text{div}_{w}p,q\rangle=-\langle p,\nabla_{w}q\rangle$.
Then,

$F(u^{k+1}+w)-F(u^{k+1})+\langle\frac{\lambda_{2}}{\mu_{2}}b_{x}^{k},\nabla_{x}w\rangle+\langle\frac{\lambda_{2}}{\mu_{2}}b_{y}^{k},\nabla_{y}w\rangle+\langle\frac{\lambda_{1}}{\mu_{1}}b_{w}^{k},\nabla_{w}w\rangle\geq
0$ (34)

Recall that $G(d)=\|d\|_{1}$, $\frac{\lambda_{1}}{\mu_{1}}b_{w}^{k}\in\partial
G(d_{w}^{k}),\frac{\lambda_{2}}{\mu_{2}}b_{x}^{k}\in\partial G(d_{x}^{k})$ and
$\frac{\lambda_{2}}{\mu_{2}}b_{y}^{k}\in\partial G(d_{y}^{k})$, so,

$\displaystyle\|d_{w}^{k}+\nabla_{w}w\|_{1}-\|d_{w}^{k}\|_{1}-\langle\frac{\lambda_{1}}{\mu_{1}}b_{w}^{k},\nabla_{w}w\rangle\geq
0$ (35)
$\displaystyle\|d_{x}^{k}+\nabla_{x}w\|_{1}-\|d_{x}^{k}\|_{1}-\langle\frac{\lambda_{2}}{\mu_{2}}b_{x}^{k},\nabla_{x}w\rangle\geq
0$ (36)
$\displaystyle\|d_{y}^{k}+\nabla_{y}w\|_{1}-\|d_{y}^{k}\|_{1}-\langle\frac{\lambda_{2}}{\mu_{2}}b_{y}^{k},\nabla_{y}w\rangle\geq
0$ (37)

Adding (34) $\sim$ (37) gives

$\|d_{w}^{k}\|_{1}+\|d_{x}^{k}\|_{1}+\|d_{y}^{k}\|_{1}+F(u^{k+1})\leq\|d_{w}^{k}+\nabla_{w}w\|_{1}+\|d_{x}^{k}+\nabla_{x}w\|_{1}+\|d_{y}^{k}+\nabla_{y}w\|_{1}+F(u^{k+1}+w)$
(38)

Suppose that $(k_{j})_{j=1,2,\ldots}$ is an increasing sequence of positive
integers such that the sequence $(u^{k_{j}})_{j=1,2,\ldots}$ converges to the
limit $\tilde{u}$. By Lemma 2, we have
$\lim_{k\rightarrow\infty}(u^{k+1}-u^{k})=0$. Therefore,
$\lim_{j\rightarrow\infty}u^{k_{j}+1}=\tilde{u}$. Moreover, we have the
following via the Lemma 4:

$\lim_{j\rightarrow\infty}d_{w}^{k}=\lim_{j\rightarrow\infty}[(d_{w}^{k}-\nabla_{w}u^{k})+\nabla_{w}u^{k}]=\nabla_{w}\tilde{u}$
(39)

also have

$\lim_{j\rightarrow\infty}d_{x}^{k}=\lim_{j\rightarrow\infty}[(d_{x}^{k}-\nabla_{x}(u^{k}-u_{p}))+\nabla_{x}(u^{k}-u_{p})]=\nabla_{x}(\tilde{u}-u_{p})$
(40)
$\lim_{j\rightarrow\infty}d_{y}^{k}=\lim_{j\rightarrow\infty}[(d_{y}^{k}-\nabla_{y}(u^{k}-u_{p}))+\nabla_{y}(u^{k}-u_{p})]=\nabla_{y}(\tilde{u}-u_{p})$
(41)

Replacing $k$ by $k_{j}$ in (38) and let $j\rightarrow\infty$, we have:

$\displaystyle\|\nabla_{w}\tilde{u}\|_{1}+\|\nabla_{x}(\tilde{u}-u_{p})\|_{1}+\|\nabla_{y}(\tilde{u}-u_{p})\|_{1}+F(\tilde{u})\leq\|\nabla_{w}(\tilde{u}+w)\|_{1}+$
$\displaystyle\|\nabla_{x}(\tilde{u}-u_{p}+w)\|_{1}+\|\nabla_{y}(\tilde{u}-u_{p}+w)\|_{1}+F(\tilde{u}+w)$

The above equations hold for all the $w\in\mathbb{R}^{N^{2}}$. On the other
hand, $u^{\ast}$ is the unique solution to the minimization problem (10).
Therefore, we must have $\tilde{u}=u^{\ast}$. Since $(u^{k})_{k=1,2,\ldots}$
is a bounded sequence, we have

$\lim_{k\rightarrow\infty}u^{k}=u^{\ast}.$ (42)

This completes the proof of the Main Theorem 1.

∎

## References

  * [1] Y. Yitzhaky, I. Dror, N. S. Kopeika. Restoration of atmospherically blurred images according to weather-predicted atmospheric modulation transfer functions. Optical Engineering, 3064-3072(36), 1997.
  * [2] J. E. Pearson. Atmospheric turbulence compensation using coherent optical adaptive techniques. Applied Optics, 622-631(15), 1976.
  * [3] K. T. Robert. Adaptive optics compensation of atmospheric turbulence: the past, the present, and the promise. Proceedings of SPIE, Atmospheric Propagation and Remote Sensing III, 404-412(2222), 1994.
  * [4] M. A. Vorontsov. Parallel image processing based on an evolution equation with anisotropic gain: integrated optoelectronic architectures. Journal of the Optical Society of America A, 1623-1637(16), 1999.
  * [5] M. A. Vorontsov, G. W. Carhart. Anisoplanatic imaging through turbulent media: image recovery by local information fusion from a set of short-exposure images. Journal of the Optical Society of America A, 1312-1324(18), 2001.
  * [6] X. Zhu, P. Milanfar. Image Reconstruction from videos distorted by atmospheric turbulence. Proceedings of SPIE, 2010.
  * [7] X. Zhu, P. Milanfar. Removing atmospheric turbulence via space-invariant deconvolution. IEEE Trans. on Pattern Analysis and Machine Intelligence, 157-170(35), 2013.
  * [8] M. Shimizu, S. Yoshimura, M. Tanaka, M. Okutomi. Super-resolution from image sequence under influence of hot-air optical turbulence. Computer Vision and Pattern Recognition, 2008.
  * [9] M. Hirsch, S. Sra, B. Sch$\ddot{o}$lkopf, S. Harmeling. Efficient filter flow for space-variant multiframe blind deconvolution. IEEE Conference on Computer Vison and Pattern Recognition, 607-614, 2010.
  * [10] N. M. Law. Lucky imaging: Diffraction-limited astronomy from the ground in the visible. Ph.D. Thesis, Cambridge University, 2003.
  * [11] D. Li, R. M. Mersereau, S. Simske. Atmospheric turbulence-degraded image restoration using principal components analysis. IEEE Geoscience and Remote Sensing Letters, 340-344(4), 2007.
  * [12] M. C. Roggemann, C. A. Stoudt, B. M. Welsh. Image-spectrum signal-to-noise-ratio improvements by statistical frame selection for adaptive-optics imaging through atmospheric turbulence. Optical Engineering, 3254-3264(33), 1994.
  * [13] D. L. Fried. Probability of getting a lucky short-exposure image through turbulence. Journal of the Optical Society of America A, 1651-1658(68), 1978.
  * [14] M. Aubailly, M. A. Vorontsov, G. W. Carhat, M. T. Valley. Image enhancement by local information fusion with pre-processing and composed metric. Proceedings of SPIE, 2008.
  * [15] S. John, M. A. Vorontsov. Multiframe selective information fusion from robust error estimation theory. IEEE Trans. on Image Processing, 577–584(14), 2005.
  * [16] M. Aubailly, M. A. Vorontsov, G. W. Carhat, M. T. Valley. Automated video enhancement from a stream of atmospherically-distorted images: the lucky-region fusion approach. Proceedings of SPIE, 2009.
  * [17] S. Gepshtein, A. Shtainman, B. Fishbain, L. Yaroslavsky. Restoration of atmospheric turbulent video containing real motion using filtering and elastic image registration. Proceeding of the Eusipco, 2004.
  * [18] D. Frakes, J. Monaco, M. Smith. Suppresion of atmospheric turbulence in video using an adaptive control grid interpolation approach. IEEE International Conference on Acoustics, Speech and Signal Processing, 1881-1884(3), 2001.
  * [19] Y. Mao, J. Gilles. Nonrigid geometric distortions correction - Application to atmospheric turbulence stabilization. Inverse Problems and Imaging, 531–546(6), 2012.
  * [20] J. Gilles, T. Dagobert, C. D. Franchis. Atmospheric turbulence restoration by diffeomorphic image registration and blind deconvolution. Proceedings of Advanced Concepts for Intelligent Vison Systems, 400-409, 2008.
  * [21] E. J. Candes, X. Li, Y. Ma, J. Wright. Robust principal compnent analysis? Journal of ACM, 1–37, 58(1), 2009
  * [22] D. Rueckert, L. Sonoda, C. Hayes, D. Hill, M. Leach, D. Hawkes. Nonrigid registration using free-form deformations: application to breast mr images. IEEE Trans. on Medical Imaging, 712-721, 18(8), 1999.
  * [23] A. Myronenko, X. Song. Intensity-based image registration by minimizing residual complexity. IEEE Trans. on Medical Imaging. 1882-1891(29), 2010.
  * [24] L. Bregman. The relaxation method of finding the common points of convex sets and its application to the solution of problems in convex optimization. USSR Computational Mathematics and Mathematical Physics, 200–217(7), 1967.
  * [25] S. Osher, M. Burger, D. Goldfarb, J. Xu, W. Yin. An iterative regularization method for total variation-based image restoration. SIAM Journal on Multiscale Modeling and Simulation, 460–489(4), 2005.
  * [26] J. Xu, S. Osher. Iterative regularization and nonlinear inverse scale space applied to wavelet-based denoising. IEEE Trans. on Image Processing, 534–544(16), 2006.
  * [27] L. He, T. C. Chang, S. Osher, T. Fang, P. Speier. MR image reconstruction by using the iterative refinement method and nonlinear inverse scale space methods. UCLA CAM Report, 06–35, 2006.
  * [28] X. Zhang, M. Burger, X. Bresson, S. Osher. Bregmanized nonlocal regularization for deconvolution and sparse reconstruction. SIAM Journal on Imaging Sciences, 253–276(3), 2010.
  * [29] X. Zhang, M. Burger, S. Osher. A unified primal-dual algorithm framework based on Bregman iteration. Journal of Scientific Computing, 1–27(46), 2010.
  * [30] P. Lions, B. Mercier. Splitting algorithm for the sum of two nonlinear operators. SIAM Journal on Numerical Analysis, 964–979(16), 1979.
  * [31] PL. Combettes, VR. Wajs. Signal recovery by proximal forward-backward splitting. Multiscale Modeling and Simulation, 1168–1200, 2005.
  * [32] G. B. Passty. Ergodic convergence to a zero of the sum of monotone operators in Hilbert space. Journal of Mathematical Analysis and Applications, 383–390(72), 1979.
  * [33] L. I. Rudin, S. Osher, E. Fatemi. Nonlinear total variation based noise removal algorithm. Journal of Physics D: Applied Physics, 259–268(60), 1992.
  * [34] D. Goldfarb, W. Yin. Parametric maximum flow algorithm for fast total variation minimization. CAAM Technical report, 2007.
  * [35] T. Goldstein, S. Osher. The split Bregman method for L1-regularized problems. SIAM Journal on Imaging Sciences, 323–343(2), 2009.
  * [36] R. Jia, H. Zhao. A fast algorithm for the total variation model of image denoising. Advances in Computational Mathematics, 231-241(33), 2010.
  * [37] M. P. Wand, M. C. Jones. Kernel smoothing. Series Monographs on Statistics and Applied Probability, Chapman and Hall, 1995.
  * [38] E. A. Nadaraya. On estimating regression. Theory of Probability and Its Application, 141-142, 1964.
  * [39] H. Takeda, S. Farsiu, P. Milanfar. Kernel regression for image processing and reconstruction. IEEE Trans. Image Processing, 349-366(16), 2007.
  * [40] T. Kato, S. Omachi, H. Aso. Asymmetric gaussian and its application to pattern recognition. Proceedings of the Joint IAPR International Workshop on Structural, Syntactic, and Statistical Pattern Recognition, 405-413, 2002.
  * [41] Q. Shan, J. Jia, A. Agarwala. High-quality motion deblurring from a single image. ACM Trans. Graphics, vol. 27, 2008.
  * [42] O. Oreifej, G. Shu, T. Pace, M. Shah. A two-stage reconstruction approach for seeing through water. Computer Vision and Pattern Recognition, 1153-1160, 2011.
  * [43] J. G. Nagy, D. P. O’leary. Fast iterative image restoration with a space-varying PSF. Advanced Signal Processing Algorithms, Architectures, and Implementations IV, 388-399, 1997.
  * [44] T. Stockham Jr. High-speed convolution and correlation. ACM Spring Joint Computer Conference, 229-233, 1966.
  * [45] W. Xue, L. Zhang, X. Mou, A. C. Bovik. Gradient magnitude similarity deviation: a highly efficient perceptual image quality index. arXiv:1308.3052, 2013.


# Ultrafast 3.5 µm fiber laser

Nathaniel Bawden  1,* Ori Henderson-Sapir  1,2 Stuart D. Jackson  3 and David
J. Ottaway1  1Department of Physics and Institute of Photonics and Advanced
Sensing, The University of Adelaide, Adelaide, SA, Australia
2Mirage Photonics, Oaklands Park, SA, Australia
3MQ Photonics, School of Engineering, Faculty of Science and Engineering,
Macquarie University, North Ryde, NSW, Australia
<EMAIL_ADDRESS>
###### Abstract

We report the first mode-locked fiber laser to operate in the femtosecond
regime well beyond 3 µm. The laser uses dual-wavelength pumping and non-linear
polarisation rotation to produce 3.5 µm wavelength pulses with minimum
duration of 580 fs at a repetition rate of 68 MHz. The pulse energy is 3.2 nJ,
corresponding to a peak power of 5.5 kW.

††journal: osajournal††articletype: Letter

##  1 Introduction

Lasers operating around 3.5 µm are of great interest because they directly
excite and probe C-H and N-O vibrational bonds. This makes them an interesting
source for polymer processing [1], molecular spectroscopy [2] and breath
analysis [3]. Prior to this work, sub-picosecond pulses at this wavelength
have only been achieved by using non-linear techniques to frequency shift
shorter wavelength lasers [4, 5]. Fiber based sources offer several potential
advantages over conventional methods of producing ultrafast pulses in the mid-
infrared wavelength region using OPOs and OPAs [6]. These include greater
average output power, improved beam quality, robustness and affordability.
Ultrafast fiber systems could therefore see broader implementation and use in
potential new applications.

Previous demonstrations of mode-locked fiber lasers at 3.5 µm have produced
picosecond scale pulses. The first demonstration by Qin et al. used a black
phosphorous saturable absorber to mode-lock with a pulse energy of 1.38 nJ at
28.91 MHz [7]. A pulse duration of 3.8 ps was estimated from spectral
broadening assuming a transform limited temporal Gaussian shape. We have
previously used the frequency shifted feedback technique to achieve tunable
pulsed operation, producing pulses as short as 53 ps with energies up to 1.38
nJ over a 215 nm wavelength range [8].

At the time of writing, the majority of ultrafast mid-infrared fibre laser
research has focused on the 3 µm band. Passive mode-locking has been achieved
using several different saturable absorber materials, producing picosecond
scale pulses with peak powers around the 1 kW level [9]. The greatest peak
power of 1.86 kW (25 ps pulse duration) was demonstrated by Tang et al. using
a semiconductor saturable absorber mirror [10]. While these systems have the
advantage of a simpler resonator configuration, shorter pulses with higher
peak power have been produced using the more complex non-linear polarisation
rotation (NPR) technique [6, 9], which utilizes polarisation controlling
optics to create an intensity dependent loss resulting in a preference for a
short pulse rather then CW operation [11]. With the inclusion of additional
amplifier and non-linear stages, NPR mode-locked fiber lasers have shown
themselves to be strong candidates for generating mid-infrared frequency combs
and supercontinuum sources [12, 13, 14, 5].

Fiber laser development generally follows a logical progression for any
particular lasing transition. Continuous wave demonstrations are followed by
Q-switching and mode-locking using saturable absorbers. Ultrashort pulse
generation is then achieved using NPR. The first systems to reach this point
in the 3 µm wavelength range utilised the 2.8 µm transition in erbium [15,
16]. Holmium-praseodymium co-doped systems followed shortly after, with the
2.9 µm operating wavelength avoiding the strong water vapor present around 2.8
µm [17]. Ultrashort pulses have recently been demonstrated using dysprosium
doped fiber [18]. Wang et al. reported pulses as short as 828 fs with a peak
power of 4.2 kW operating near 3.1 µm, the longest wavelength sub-picosecond
fiber laser prior to this work.

In this work, we advance the central operating wavelength of femtosecond fiber
lasers to 3.5 µm, well within the 3 µm to 5 µm atmospheric window. The NPR
mode-locked system uses an erbium-doped zirconium fluoride fiber and the dual-
wavelength pumping method originally developed by Henderson-Sapir et al. [19].
Direct operation at this wavelength range gives a high brightness, mode-locked
spectrum which could be used to directly probe the absorption lines of
specific molecules such as $\text{CH}_{4}$ and HCl, opening up possible
applications in precision spectroscopy and remote sensing [4, 20].

##  2 Experimental setup and results

The experimental setup for NPR mode-locking experiments at 3.5 µm is shown in
Fig. 1. This configuration is based on those used in the first NPR mode-
locking experiments at 3 µm and adapted for 3.5 µm [15, 16].

Figure 1: Schematic diagram of the NPR mode-locked laser. Inset: A simplified
energy level diagram showing the dual-wavelength pumping scheme. HWP - half
wave plate. QWP - quarter wave plate. DM - dichroic miror. OC - output
coupler. $\text{P}_{1}$ \- 977 nm pump. $\text{P}_{2}$ \- 1973 nm pump.

The 977 nm and 1973 nm pump sources are a commercial fiber-coupled laser diode
(LIMO HLU30F200-977) and a thulium-doped silica fiber laser (Mirage Photonics
TFL-2000). The combined pumps are passed through a 45 degree dichroic mirror
and launched into the 1 mol% $\text{Er}^{3+}\text{:ZrF}_{4}$ fiber by a custom
aspheric lens. The 2.3 m long, double-clad fiber (Le Verre Fluoré) has a 16 µm
diameter core with a 0.12 NA at 3.5 µm and a 240/260 µm double-truncated
circular inner cladding. The 977 nm and 1973 nm pumps are launched into the
inner cladding and core, respectively. The output end of the fiber is angle
polished to prevent back reflections and parasitic lasing at 2.8 µm. Light
from this end is collimated by another aspheric lens and incident on a 45
degree output coupler mirror with 77% reflectivity at 3.5 µm. A second 45
degree dichroic mirror then completes the ring laser resonator. The length of
the free space section in the laser is approximately 0.9 m long.

A free space Faraday rotator (Thorlabs) and two high power polarisers
(FastPulse Technology) form an optical isolator that is placed between the two
dichroic mirrors, providing unidirectional operation. The transmission of the
complete isolator system is measured to be 86% at 3.5 µm for linearly
polarised light. Two quarter wave plates and a half wave plate (Thorlabs) are
also added as shown in Fig. 1 to implement NPR mode-locking. The quarter wave
plate at the input end of the isolator is present to reduce the loss at the
following polariser and was necessary to achieve mode-locked operation. The
laser output is passed through a 3 µm longpass filter to remove the residual
pump light when making measurement of the mode-locked system.

For a 977 nm pump power of 6.8 W and a 1973 nm pump power of 5.5 W, mode-
locked operation can be started with careful adjustment of the wave plates and
knocking one of the resonator mirrors which creates an intensity spike that
can preferentially pass through the polarisation controlling optics and the
fiber. The average output power is measured to be 216 mW during mode-locked
operation at these pump powers. A 2.4 µm to 3.2 µm bandpass filter and grating
spectrometer were used to check for parasitic lasing at 2.8 µm. No lasing at
this wavelength was observed during mode-locked operation.

Time domain measurements of the mode-locked output made with a 1 ns rise time
photodetector (VIGO PEM-10.6) show continuous pulsed operation. A peak at 68.0
MHz is displayed in the RF spectrum, equal to the inverse round trip time of
the ring resonator, with harmonics extending up to the 1 GHz limit of the
spectrum analyser (Fig. 2). The fundamental beat frequency has a signal to
noise ratio of nearly 70 dB for a RBW of 3 kHz, comparable to NPR mode-locked
systems around 3 µm [15, 17]. An average pulse energy of 3.2 nJ is calculated
from the output power and repetition rate.

Figure 2: RF spectrum of the mode-locked laser output with a fundamental
repetition rate of 68.0 MHz.

An autocorrelator based on two-photon absorption was constructed to measure
the pulse width (same as in [8]). The fit to the autocorrelation trace in Fig.
3 gives a FWHM of 580 fs for a $\text{sech}^{2}$ shaped pulse after
deconvolution. A peak power of 5.5 kW is obtained based on these values, on
par with the results obtained at 2.8 µm and 3.1 µm [16, 18].

Figure 3: Autocorrelation trace of the mode-locked pulses at an average output
power of 216 mW. The corresponding pulse width is 580 fs FWHM for a
$\text{sech}^{2}$ shaped pulses. The 977 nm and 1973 nm pump powers were 6.8 W
and 5.5 W, respectively.

The optical spectrum of the mode-locked laser was measured using a grating
spectrometer (Fig. 4). The operating wavelength is centred on 3.54 µm and
broadens to 23 nm FWHM when mode-locking compared to the 0.5 nm linewidth
typically measured during CW operation. The time-bandwidth product for these
pulses is 0.32 indicating that they are close to transform limited (0.315 for
$\text{sech}^{2}$).

Figure 4: Optical spectrum of the mode-locked laser at maximum pump power
showing Kelly sidebands. The spectrum broadens to 23 nm FWHM.

The narrow peaks in the optical spectrum either side of the central wavelength
indicate that this system operates in the soliton regime. The spacing of the
sidebands from the center frequency can be used to calculate the resonator
dispersion using the equation [21]:

$(f_{m}-f_{c})^{2}=\frac{m}{\pi\absolutevalue{B_{2}}}-\frac{1}{4\pi^{2}\tau^{2}}$
(1)

where $m$ is the order of the sideband, $f_{m}$ is the corresponding sideband
frequency, $f_{c}$ is the frequency of the central peak, $B_{2}$ is the total
resonator dispersion and $\tau$ is the pulse width (for a pulse given by
$\text{sech}^{2}\quantity(\frac{t}{\tau})$, the FWHM is $1.763\tau$). By
fitting to $(f_{m}-f_{c})^{2}$ as a function of $m$ we obtain a total
resonator dispersion of $\quantity(-0.44\pm 0.02)\times
10^{-24}\,\text{s}^{2}$. If we assume the total resonator dispersion is due
entirely to the fiber, dividing this result by the 2.3 m fiber length gives
the group velocity dispersion $\beta_{2}$ of the fiber as $\quantity(-0.191\pm
0.009)\times 10^{-24}\,\text{s}^{2}\,\text{m}^{-1}$. This agrees with the
value of $-0.181\times 10^{-24}\,\text{s}^{2}\,\text{m}^{-1}$ we calculate
from numerical fiber eigenmode analysis [22]. The refractive indices used in
this calculation were obtained from data provided by fiber manufacturer.

The pulse energy and duration of solition lasers are often estimated using the
soliton area theorem [11]:

$E\tau=\frac{2\absolutevalue{\beta_{2}}}{\gamma},$ (2)

where $\beta_{2}$ is the group velocity dispersion, $\gamma$ is the nonlinear
coefficient and $\tau$ is the pulse duration ($\text{FWHM}=1.763\tau$). We
calculate $\gamma=$ $0.83\times 10^{-4}\,\text{W}^{-1}\text{m}^{-1}$ using a
mode field radius of 12.0 µm and a nonlinear refractive index of $2.1\times
10^{-20}\,\text{m}^{2}\text{W}^{-1}$ [23]. We use the value of $-0.191\times
10^{-24}\,\text{s}^{2}\text{m}^{-1}$ for $\beta_{2}$ that was calculated
earlier from the Kelly sidebands present in the optical spectrum. Using the
measured pulse energy of 3.2 nJ and accounting for the 77% reflectivity output
coupler mirror, the soliton area theorem gives a pulse duration of 590 fs.
Although this agrees well with our measured value of 580 fs, this is still
only an approximation as these are not strictly solitons but are more
accurately described as average solitons [24]. In this operating regime, the
pulse characteristics are determined by the balance between the accumulated
phase in a resonator round trip due to the averaged dispersion and
nonlinearity [25]. The pulse energy varies significantly as it propogates
around the resonator due to gain/loss and the average is markedly lower than
that at the output end of the fiber. The pulse duration and spectral width
however, remain relatively constant throughout a resonator round trip, based
on numerical simulations by Woodward [25]. The pulses therefore disagree with
the soliton area theorem at many points in the resonator such that the
averaged effects of self phase modulation balance anomalous dispersion.

Experiments using 72% and 87% reflectivity output coupler mirrors resulted in
mode-locked pulses with different durations and energies. The results are
summarised in Table 1. For the 87% reflectivity mirror, instabilities in the
mode-locked output limited the incident pump powers. The optical spectrum in
this case contained stronger sidebands, up to ten times the intensity of the
central operating wavelength. The larger output coupling fraction of the 72%
reflectivity mirror resulted in the greatest pulse enrgy of 3.6 nJ.

Table 1: Measured parameters for mode-locked operation with various output coupler mirrors. R [%] |  $\text{P}_{1}$ [W] |  $\text{P}_{2}$ [W] | Power [mW] | Energy [nJ] | Duration [fs] | Peak power [kW]
---|---|---|---|---|---|---
72 | 8.0 | 6.1 | 244 | 3.6 | 650 | 5.5
77 | 6.8 | 5.5 | 216 | 3.2 | 580 | 5.5
87 | 6.2 | 5.3 | 141 | 2.1 | 600 | 3.5

Our NPR mode-locking results presented in this work significantly improve upon
previous mode-locking demonstrations at this wavelength because the shorter,
sub-picosecond pulses result in a greater peak power. We have shown that NPR
is a robust method for short pulse generation well into the mid-infrared. Our
peak power and pulse duration are comparable to those demonstrated in similar
experiments at 2.8 µm and 3.1 µm, respectively, even though the erbium 3.5 µm
transition is a lower gain lasing transition.

Useful information is also obtained about the fluorozirconate fiber. The
measured group velocity dispersion begins to create a useful set of parameter
values for this wavelength range when combined with the measured values from
similar experiments. These are summarised in Table 2.

Table 2: Measured values of the group velocity dispersion from NPR mode-locking experiments in the mid-infrared. Wavelength [µm] |  $\beta_{2}$ [$\text{s}^{2}\text{m}^{-1}$] | Ref.
---|---|---
2.8 | 0.094 | [16]
3.1 | 0.120 | [18]
3.5 | 0.191 | This work

The erbium 3.5 µm transition has previously been shown to have a tunability of
over 450 nm in CW operation from 3.33 µm and almost reaching 3.8 µm [26].
Future experiments with increased pump powers and dispersion management should
therefore be able to achieve further spectral broadening and hence shorter
pulses while still allowing for wavelength tuning, making this system more
appealing for spectroscopy applications.

We have demonstrated mode-locked operation of a fiber laser at 3.5 µm using
NPR. The average output power and repetition rate were 216 mW and 68 MHz,
corresponding to a pulse energy of 3.2 nJ. Autocorrelation measurements gave a
pulse duration of 580 fs, the first sub-picosecond pulses from a fiber laser
operating directly at this wavelength. The resulting peak power for the laser
was 5.5 kW.

## Acknowledgments

The authors thank Matthew R. Majewski and Robert I. Woodward for helpful
discussions. The authors also acknowledge the expertise, equipment, and
support provided by the Australian National Fabrication Facility (ANFF) at The
University of Adelaide. This research was supported by the Australian Research
Council, the Asian Office of Aerospace R&D.

## Funding

Asian Office of Aerospace R&D (AOARD) Grants FA-9550-20-1-0160 and
FA2386-19-0043.

## Diclosures

OHS: Mirage Photonics (I,P)

## References

  * [1] C. Frayssinous, V. Fortin, J.-P. Bérubé, A. Fraser, and R. Vallée, “Resonant polymer ablation using a compact 3.44 µm fiber laser,” Journal of Materials Processing Technology 252, 813–820 (2018).
  * [2] F. Adler, P. Masłowski, A. Foltynowicz, K. C. Cossel, T. C. Briles, I. Hartl, and J. Ye, “Mid-infrared Fourier transform spectroscopy with a broadband frequency comb,” Optics Express 18, 21861 (2010).
  * [3] C. Wang and P. Sahay, “Breath analysis using laser spectroscopic techniques: Breath biomarkers, spectral fingerprints, and detection limits,” Sensors 9, 8230–8262 (2009).
  * [4] A. Schliesser, N. Picqué, and T. W. Hänsch, “Mid-infrared frequency combs,” Nature Photonics 6, 440–449 (2012).
  * [5] S. Duval, J.-C. Gauthier, L.-R. Robichaud, P. Paradis, M. Olivier, V. Fortin, M. Bernier, M. Piché, and R. Vallée, “Watt-level fiber-based femtosecond laser source tunable from 2.8 to 3.6 µm,” Optics Letters 41, 5294 (2016).
  * [6] S. D. Jackson and R. K. Jain, “Fiber-based sources of coherent MIR radiation: Key advances and future prospects (invited),” Optics Express 28, 30964 (2020).
  * [7] Z. Qin, T. Hai, G. Xie, J. Ma, P. Yuan, L. Qian, L. Li, L. Zhao, and D. Shen, “Black phosphorus Q-switched and mode-locked mid-infrared Er:ZBLAN fiber laser at 3.5 µm wavelength,” Optics Express 26, 8224 (2018).
  * [8] O. Henderson-Sapir, N. Bawden, M. R. Majewski, R. I. Woodward, D. J. Ottaway, and S. D. Jackson, “Mode-locked and tunable fiber laser at the 3.5 µm band using frequency-shifted feedback,” Optics Letters 45, 224 (2020).
  * [9] X. Zhu, G. Zhu, C. Wei, L. V. Kotov, J. Wang, M. Tong, R. A. Norwood, and N. Peyghambarian, “Pulsed fluoride fiber lasers at 3 µm [Invited],” Journal of the Optical Society of America B 34, A15 (2017).
  * [10] P. Tang, Z. Qin, J. Liu, C. Zhao, G. Xie, S. Wen, and L. Qian, “Watt-level passively mode-locked Er3+-doped ZBLAN fiber laser at 2.8 µm,” Optics Letters 40, 4855 (2015).
  * [11] H. A. Haus, “Mode-locking of lasers,” IEEE Journal of Selected Topics in Quantum Electronics 6, 1173–1185 (2000).
  * [12] J. Huang, M. Pang, X. Jiang, F. Köttig, D. Schade, W. He, M. Butryn, and P. S. J. Russell, “Sub-two-cycle octave-spanning mid-infrared fiber laser,” Optica 7, 574 (2020).
  * [13] D. D. Hudson, S. Antipov, L. Li, I. Alamgir, T. Hu, M. El Amraoui, Y. Messaddeq, M. Rochette, S. D. Jackson, and A. Fuerbach, “Toward all-fiber supercontinuum spanning the mid-infrared,” Optica 4, 1163 (2017).
  * [14] R. I. Woodward, D. D. Hudson, A. Fuerbach, and S. D. Jackson, “Generation of 70-fs pulses at 2.86 µm from a mid-infrared fiber laser,” Optics Letters 42, 4893 (2017).
  * [15] T. Hu, S. D. Jackson, and D. D. Hudson, “Ultrafast pulses from a mid-infrared fiber laser,” Optics Letters 40, 4226 (2015).
  * [16] S. Duval, M. Bernier, V. Fortin, J. Genest, M. Piché, and R. Vallée, “Femtosecond fiber lasers reach the mid-infrared,” Optica 2, 623 (2015).
  * [17] S. Antipov, D. D. Hudson, A. Fuerbach, and S. D. Jackson, “High-power mid-infrared femtosecond fiber laser in the water vapor transmission window,” Optica 3, 1373 (2016).
  * [18] Y. Wang, F. Jobin, S. Duval, V. Fortin, P. Laporta, M. Bernier, G. Galzerano, and R. Vallée, “Ultrafast Dy3+:fluoride fiber laser beyond 3 µm,” Optics Letters 44, 395 (2019).
  * [19] O. Henderson-Sapir, J. Munch, and D. J. Ottaway, “Mid-infrared fiber lasers at and beyond 3.5 µm using dual-wavelength pumping,” Optics Letters 39, 493 (2014).
  * [20] S. D. Jackson, “Towards high-power mid-infrared emission from a fibre laser,” Nature Photonics 6, 423–431 (2012).
  * [21] N. J. Smith, K. J. Blow, and I. Andonovic, “Sideband generation through perturbations to the average soliton model,” Journal of Lightwave Technology 10, 1329–1333 (1992).
  * [22] D. Marcuse, “Interdependence of waveguide and material dispersion,” Applied Optics 18, 2930 (1979).
  * [23] C. Agger, C. Petersen, S. Dupont, H. Steffensen, J. K. Lyngsø, C. L. Thomsen, J. Thøgersen, S. R. Keiding, and O. Bang, “Supercontinuum generation in ZBLAN fibers—detailed comparison between measurement and simulation,” Journal of the Optical Society of America B 29, 635 (2012).
  * [24] S. M. J. Kelly, K. Smith, K. J. Blow, and N. J. Doran, “Average soliton dynamics of a high-gain erbium fiber laser,” Optics Letters 16, 1337 (1991).
  * [25] R. I. Woodward, “Dispersion engineering of mode-locked fibre lasers,” Journal of Optics 20, 033002 (2018).
  * [26] O. Henderson-Sapir, S. D. Jackson, and D. J. Ottaway, “Versatile and widely tunable mid-infrared erbium doped ZBLAN fiber laser,” Optics Letters 41, 1676 (2016).