1. Introduction

The non-uniform distribution of temperature yields random atmosphere refractive-index fluctuations, resulting in optical turbulence that distorts the optical wave front. This phenomenon causes spatiotemporal distortions and spatial image blur in long-distance imaging. Techniques to mitigate optical turbulence have resulted in image degradation, such as the hardware-based adaptive optics approach and the software-based numerical processing approach. The adaptive optics approach [1] is applied in large telescopes and reduces the turbulence effects in real time. However, its working mode makes this approach unsuitable for horizontal imaging in atmospheric boundary layers.

The numerical processing approach can be divided into two classes. The first is video-to-image type methods. These methods reconstruct a single high-quality image from atmospherically distorted video, such as lucky imaging [2], speckle imaging [3, 4], image restoration methods [5, 6], non-rigid registration methods [7], registration and fusion hybrid methods [8]. The other class consists of video-to-video type methods, namely video stabilization methods, such as the super-resolution method [9], the low-rank matrix decomposition based method [10], the Sobolev sharpening flow based method [11] and the phase based method [12]. In this paper, we focus on investigating video stabilization problems. The super-resolution method first generates a reference frame using a temporal pixel-wise median over a block of frames. Video stabilization is obtained by computing a map of pixel displacements between every frame and its reference frame. To some extent, the effects of stabilization depend on elastic local registration results. The low-rank matrix decomposition method decomposes the turbulence sequence into three components: the background, the turbulence and the object. However, matrix decomposition is time consuming and complicated for high-resolution imaging. The Sobolev sharpening flow-based method proposes a PDE model for video stabilization, and numerical acceleration techniques are used to improve efficiency. However, the target properties in infrared images and the contrasts in panchromatic images are affected using the Sobolev sharpening flow, and details cannot be restored in the event of destructive turbulent degradation. The phase-based method alleviates turbulence distortion through local phase temporal filtering within a complex-valued steerable pyramid. However, the results with this method can be affected by noisy disturbances because the method did not process the amplitude variations.

Wave-front angle-of-arrival (AOA) fluctuations are the phase related parameters of optical wave propagating through atmospheric turbulence, which have two primary degrading effects on recorded video sequences: spatiotemporal distortions and long-exposure imaging blur. The phase of a signal implies the structure and motion information [13]; therefore, phase information has been widely applied in signal processing and image processing fields [8, 12, 14–17]. In one-dimensional (1-D) signal processing, the Hilbert transform based analytic signal is used as a complex-valued representation of a 1-D signal for local frequency analysis. For 2-D signal processing problems, such as image processing, the Riesz transform based monogenic signal [18] has been proposed for analyzing the local amplitude and the local phase.

In this paper, we present an image complex pyramid representation, the Laplacian-Riesz pyramid, which is used for video stabilization in turbulence. The Laplacian pyramid is characterized by its compact structure and perfect reconstruction [19], and it is completely built in the spatial domain; consequently, it possesses high computational efficiency. The Riesz transform is a steerable Hilbert transformer and can compute phase-shift and translate image features in the direction of the dominant orientation at every pixel [17]. Therefore, in this work we combine the Laplacian pyramid with the Riesz transform to form a Laplacian-Riesz pyramid (LRP) for image representation. Based on this complex pyramid representation, local phase and amplitude variations of the coefficients of a Laplacian-Riesz pyramid over time can be extracted. Local phase variations correspond to turbulent-induced wiggles. Additionally, AOA effects also result in local amplitude fluctuations. Therefore, low-pass filters are used to temporally attenuate these local phase and amplitude variations to remove distracting changes. The high spatial frequencies are undesirable for the temporal filtering due to camera sensor and quantization noise. Therefore, we decompose each frame into the pyramid structure first for the purpose of temporal processing.

The paper is organized as follows: Section 2 describes the Laplacian-Riesz pyramid-based video stabilization framework in detail. Experimental results are given in Section 3 to show the performance compared with other methods; finally, we conclude and discuss the direction of future research in Section 4.

2. Laplacian-Riesz pyramid based video stabilization method

The proposed video stabilization framework contains three steps (see the diagram in Fig. 1), (1) the Laplacian-Riesz pyramid representation of each frame of a video; (2) the temporal filtering of local amplitude and local phase; (3) the reconstruction of the corrected frame.

 

Fig. 1 Block diagram for the proposed video stabilization framework.

Download Full Size | PDF

2.1 Laplacian-Riesz pyramid

2.1.1 Laplacian pyramid

A Laplacian pyramid [19] is the subtraction between the lower level’s up-sampling and the higher level in a Gaussian pyramid; this means the Laplacian pyramid decreases the correlation between each level. The Laplacian pyramid forms a series of band-pass components and a low-pass component of the original image. Without losing generality, we used a 5 × 5 Gaussian kernel w(.) as the low-pass filter to build the pyramid. We set the standard deviation of w(.) to 2 because it provides adequate filtering with low reconstruction errors. When the turbulence effects are particularly severe, the standard deviation should be larger to avoid grid noise.

The original frame I convolves with w to get the first level $g_1$ of the Gaussian pyramid, which is treated as the highest level here. The function REDUCE denotes a down-sample and low-pass process, which turns a higher sub-band $g_{k-1}$ with a size of $2C_k\times 2R_k$ into a lower sub-band $g_k$ with a size of $C_k\times R_k$

Inversely, the function EXPAND denotes the up-sample and interpolation operation, which expands the lower sub-band $g_k$ with a size of $C_k\times R_k$ to a higher sub-band $g_{k-1}$ with a size of $2C_k\times 2R_k$ by interpolating new values between the given values:

Suppose each sub-band of the Laplacian pyramid is $L^{\left\{k\right\}}$_{, $k=1,2,\dots ,N$}. N indicates the maximum level. The Laplacian pyramid can then be built as follows:

Equations (3)-(5) show that the Laplacian pyramid can be reconstructed exactly to recover the original image. $g_{N-1}$ can be computed by adding $L^{\left\{N-1\right\}}$ to the EXPAND of $L^{\left\{N\right\}}$. Then, $g_{N-2}$ is yielded by adding the EXPAND of $g_{N-1}$ to $L^{\left\{N-2\right\}}$. The rest can be done in the same manner:

Finally, a reconstructed image $I_{out}$ is recovered from the pyramid:

2.1.2 Riesz transform

The Riesz transform is the multi-dimensional extension of the Hilbert transform. The Riesz operator’s frequency response is $-j\omega /\Vert \omega \Vert$. That is, there are N Riesz filters for an N-dimension signal, and when N = 1 the Riesz transform is equal to the Hilbert transform. As for a 2D signal $I(x)$, $x=\left[x,y\right]$, the Riesz transform is expressed as:

The frequency response of the Riesz operator indicates that the Riesz transform is a 90-degree phase-shift system. The phase-shift property enables the Riesz transform to compute the local amplitude and local phase along the dominant orientation $\theta$ at every point of an image I.

The input I and its Riesz transform (R₁,R₂) together form a triple that can output the local amplitude A and local phase $\phi$:

The local amplitude A implies the local intensity and energy, and it is defined as:

The local phase $\phi$ contains the structural and motion information, $\phi$ is defined as:

2.1.3 Laplacian-Riesz pyramid

After building the Laplacian pyramid, we performed the Riesz transform to each sub-band $L^{\left\{k\right\}}$. The LRP can then be created: $\left\{R_1{}^{\left\{k\right\}},R_2{}^{\left\{k\right\}}\right\}$. There are two components on each level of the LRP. Figures 2(a) and 2(b) illustrate the LRP structure. With Eqs. (10) and (11), the LRP provides a multi-resolution method to analyze the local amplitude $A^{\left\{k\right\}}$ and local phase ${\phi }^{\left\{k\right\}}$, as is shown in Fig. 2(c).

 

Fig. 2 LRP and Local information (a) Input image, (b) LRP, (c) Local amplitude and local phase.

Download Full Size | PDF

2.2 Temporal filtering of the local amplitude and local phase

The impact of turbulence on imaging appears as pixel intensity blinking and pixel wiggles, which can be represented and analyzed by the local amplitude and the local phase of the Riesz transform. In this paper, the Riesz transform is applied to reveal the turbulence distortion effects.

Assume we are given an image sequence without turbulence effects: $I\left(x,y,t\right)$, where t denotes the sequence frame number. In Eq. (9), $I\left(x,y,t\right)=A\left(x,y,t\right)cos\phi \left(x,y,t\right)$, $A\left(x,y,t\right)$ and $\phi \left(x,y,t\right)$ are the amplitude and the phase of the Riesz transform of non-turbulence videos. $A\left(x,y,t\right)$ and $\phi \left(x,y,t\right)$ are the temporal constants in areas that are not influenced by turbulence. Therefore, we use $A\left(x,y\right)$ and $\phi \left(x,y\right)$ to replace $A\left(x,y,t\right)$ and $\phi \left(x,y,t\right)$.

When the image sequence is influenced by turbulence, the degraded effects cause the local brightness variation, which is represented by $\Delta \left(x,y,t\right)$, and also create a local ‘dance’ $\delta \left(x,y,t\right)$. Turbulence effects can be expressed as:

Using the Riesz transform and Eqs. (10) and (11), the local amplitude $A\left(x,y\right)+\Delta \left(x,y,t\right)$ and local phase $\phi \left(x,y\right)+\delta \left(x,y,t\right)$ can be extracted. Figure 3 illustrates a video’s local amplitude and the local phase temporal variance. We found that the local amplitude and local phase in turbulence videos are displayed like high-frequency fluctuations. In this case, low-pass filtering may be reasonable to eliminate the distortions and retain the DC components: $A\left(x,y\right)$ and $\phi \left(x,y\right)$.

 

Fig. 3 The high-frequency fluctuations of turbulence distortion. (a) One frame of Sequence 2. (b) Local amplitude and local phase temporal varying of the green point .

Download Full Size | PDF

In our experiments, the first order Butterworth filter was chosen as the low-pass filter for comparison with the results of [12]. Of course, other low-pass filters are also feasible. In general, the higher the order of a filter, the more accurate the approximation, but the calculation cost also rises. Therefore, a balance must be found between the performance and the cost. For this we chose the first order Butterworth filter.

The first order Butterworth low-pass filter’s frequency response can be expressed as:

The low-pass filtering is carried out by IIR filter, whose difference equation is:

After temporal filtering, the original $\left\{A^{\left\{k\right\}},{\phi }^{\left\{k\right\}}\right\}$ of the distorted frame I is modified to be $\{y_A{}^{\left\{k\right\}},y_{\phi }{}^{\left\{k\right\}}\}$. Then the filtered Laplacian pyramid can be expressed as $L_{out}{}^{\left\{k\right\}}=y_A{}^{\left\{k\right\}}cos{y}_{\phi }{}^{\left\{k\right\}}$. The undistorted frame I_out is rebuilt using Eqs. (6) and (7).

3. Experiments

In this section, we conducted experiments to illustrate the performance of the proposed approach using both simulated and real videos. Results of this video-to-video type approach in [11] and [12] are also shown for comparison (the outputs of method [11] were produced through our own implementation; the outputs of method [12] were generated using the original code [20]). All experiments were performed on a PC with an Intel i5-4460 processor and 8G RAM (the algorithm uses only a single CPU core). We implemented these methods in Matlab 2013b.

3.1 Simulated experiments

To quantitatively evaluate the method performance, we generated a set of turbulence degraded image sequences in horizontal views. This simulation method is a hybrid one [21, 22] that combines the physical models of optical waves under Kolmogorov turbulence with imaging processing algorithms (image convolution and interpolation methods). Simulation uses a latent sharp image (280 × 280) as an input image (as shown in Fig. 4) and produces image sequences with 100 frames degraded by specified turbulence, which includes imaging distance, optics diameter, wavelength, turbulence spectral power law, and turbulence strength. Simulation parameters were designed as: the imaging distance T = 1000m, wavelength λ = 0.5μm, turbulent strength (strong turbulence) C_n² = 2.5 × 10⁻¹³m^-2/3, the aperture diameter D = 0.3m. In the following, we used other parameter settings as described in [21]. We refer interested readers to [21, 22] for details.

 

Fig. 4 Latent sharp image used for simulation.

Download Full Size | PDF

The cutoff frequency ω_c is set to several values in order to show the effects of this parameter. Figure 5 shows the performance of the different ω_c. The yellow line profile in each image is plotted in time on the right (for a fixed x, plot of y vs t). As we can see, the smaller ω_c is, the more stable the corresponding result.

 

Fig. 5 The effects of ω_c (see Visualization 1) (a) one frame from the simulated video, (b) ω_c = 0.1, (c) ω_c = 0.3, (d) ω_c = 0.5, (e) ω_c = 0.8.

Download Full Size | PDF

Figure 6 shows a comparison between the other two methods, and Figs. 6(a) and 6(b) show the outputs of [12] and [11], respectively. As shown [12], effectively alleviates the distortions, but the local brightness is obviously still blinking. The method in [11] can reduce the distortion to some extent, but the result is not stabilized enough or produces cross-color during restoration. Our proposed method outperforms the others.

 

Fig. 6 Simulated comparison results (see Visualization 2) (a) [12] output, (b) [11] output, (c) the proposed output (with ω_c = 0.1).

Download Full Size | PDF

3.2 Quantitative comparisons

To quantitatively evaluate the performance, we measured the Peak Signal-to-Noise Ratio (PSNR) between each frame and reference frame. For the simulated video, we used the latent image as the reference; the results are listed in Table 1. For the real videos, we used the average frame of the sequence as the reference. Table 2 reports the comparison results. It is clear that our method obtained state-of-art results and outperformed the other two methods.

Table 1. Average PSNR (dB) for differentω_c values

View Table | View all tables in this article

Table 2. Comparison of the Average PSNR (dB)

View Table | View all tables in this article

To evaluate the robustness to noise, we added white Gaussian noise with different SNRs to the simulated video. Table 3 shows the performance. Figure 7 shows the visual results.

Table 3. Performance with white Gaussian noise

View Table | View all tables in this article

 

Fig. 7 Results of additional noises (see Visualization 3, Visualization 4, Visualization 5) (a) one frame of simulated video with noise, (b) one frame of [11] output, (c) one frame of [12] output, (d) one frame of the proposed output.

Download Full Size | PDF

In Table 4, the total consuming times are presented. As shown, the proposed approach is the fastest one.

Table 4. Comparison of running time (s)

View Table | View all tables in this article

3.3 Real video experiments

This work tested three real sequences distorted by atmospheric turbulence to illustrate the performance of the proposed approach. Sequence 1 consists of 101 frames degraded by strong turbulence effects. Sequences 2 and 3 are two infrared sequences [23]. Sequence 3 contains a moving man and two moving dogs. Figures 8, 9, and 10 show three turbulence mitigation experimental results. The cutoff frequency is set to 0.3Hz for Sequence 1 and 2, 0.8Hz for Sequence 3. Visual effect comparisons show that the proposed approach obtained the best turbulence mitigation result.

 

Fig. 8 Sequence1 stabilization results (see Visualization 6) (a) one frame from sequence 1, (b) [12] output, (c) [11] output, (d) the proposed output.

Download Full Size | PDF

 

Fig. 9 Sequence 2 stabilization results (see Visualization 7) (a) one frame from sequence 2, (b) [12] output, (c) [11] output, (d) the proposed output.

Download Full Size | PDF

 

Fig. 10 Sequence 3 stabilization results (see Visualization 8) (a) one frame from sequence 3, (b) [12] output, (c) [11] output, (d) the proposed output.

Download Full Size | PDF

4. Discussion and conclusions

In this work we proposed a fast method to stabilize video distorted by atmospheric turbulence. In this approach, each frame was decomposed into non-oriental sub-bands with a Laplacian pyramid that was built with spatial Gaussian kernel convolution and difference, then local amplitude and phase along the dominant orientation was computed directly with the Riesz transform. Therefore, the efficiency of the proposed approach is much higher than the phase-based method [12], which used eight orientation complex steerable pyramid representations. Additionally, our method used temporal filtering of the local amplitude and local phase to mitigate the temporal oscillation of atmospheric turbulence. As a result, our method obtained excellent visual quality.

We noted that our method retained the moving objects while stabilizing turbulence distortions in Sequence 3. However, the movement was weakened to some extent because low-pass filtering is harmful to movement. Eliminating the turbulence effects without harming moving objects remains a challenge. Maybe in the future, when turbulence fluctuation frequencies are known, the low-pass filter can be replaced by a band-stop filter to hold back the turbulence effects. Consequently, moving objects would be better preserved.

Temporal oscillation and blurry images are two of the main degradation effects of atmospheric turbulence. Results of video stabilization using the proposed method still suffer from turbulence-induced blur. To obtain a latent sharp frame, we need to know the atmospheric modulation transfer function (MTF) [24]. And then, deconvolution techniques should be used to restore each sharp frame [25, 26]. These remain open problems with plenty of room for further investigation.