A new approach for atmospheric turbulence removal using low-rank matrix factorization

LRMF

Consider $X=\left[x_1,\dots ,x_n\right]\in R^{d\times n}$ representing a sequence of turbulent images. Here, n stands for the number of images and d signifies the dimension of each turbulent image, with each image expressed as a vector. The input dataset can be decomposed into $X=U{V}^T$, where U is the basis matrix and V is the coefficient matrix (Yong et al., 2017). Consequently, each element within the input matrix X can be modeled as follows:

(2) $x_{ij}=u_i{v}_j^T+{\epsilon }_{ij}$where $u_i$ and $v_j$ represent the $i^{th}\mathrm{a}\mathrm{n}\mathrm{d}{j}^{th}$ row vectors, respectively, of matrices U and V. The term ${\epsilon }_{ij}$ corresponds to the noise within the $j^{th}$ pixel of the $i^{th}$ turbulent image. This noise is the spatiotemporal varying blur inherent in turbulence images that is modeled as a Mixture of Gaussians (MOG).

Hence, the noise can be formulated as follows:

(3) $${\epsilon }_{ij}\sim {\sum }_{k=1}^K{\pi }_{k}N\left(0,{\sigma }_k^2\right)$$where $N\left(0,{\sigma }_k^2\right)$ denotes a Gaussian distribution with zero mean and ${\sigma }^2$ variance and ${\pi }_k$ signifies the mixing proportion of the MOG. By merging Eqs. (2) and (3), the following relationship is established:

(4) $$x_{ij}\sim {\sum }_{k=1}^K{\pi }_{k}N(x_{ij}|u_i{v}_j^T,{\sigma }_k^2)$$

Subsequently, by considering (Eq. (5)) as a cost function and finding its maximum value, the estimations for the Gaussian noise parameters (Π, Σ) and the low-rank matrices U and V are determined. To find the maximum value of Eq. (5), the EM method is employed, as elaborated upon in the subsequent section.

(5) $$\underset{U,V,\mathrm{\Pi },\mathrm{\Sigma }}{max}L\left(U,V,\mathrm{\Pi },\mathrm{\Sigma }\right)=max{\sum }_{i,j}\log {\sum }_{k=1}^K{\pi }_{k}N(x_{ij}|u_i{v}_j^T,{\sigma }_k^2)$$where $\mathrm{\Pi }=\left\{{\pi }_1,{\pi }_2,\dots ,{\pi }_K\right\},$ and $\mathrm{\Sigma }=\left\{{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_K\right\}$.

Estimation of LRMF model parameters using the EM method

Equation (5) can be solved through the utilization of the EM method. This iterative technique computes the parameters of Eq. (5) in two steps: the Expectation (E) step and the Maximization (M) step.

E step: Let $z_{ijk}\in \left\{0,1\right\}$ $\left({\sum }_k{z}_{ijk}=1\right)$: $z_{ijk}$ takes the value of one ( $z_{ijk}=1$) if the ${\epsilon }_{ij}$noise is generated by the $k^{th}$ Gaussian distribution. If not, it equals zero ( $z_{ijk}=0$). So, the probability of the partial membership of the noise of data $x_{ij}$( ${\epsilon }_{ij}$) to the $k^{th}$ Gaussian model can be formulated as follows:

(6) $$E\left(z_{ijk}\right)={\gamma }_{ijk}={\displaystyle \frac{{\pi }_{k}N(x_{ij}|u_i{v}_j^T,{\sigma }_k^2)}{{\sum }_{k=1}^K{\pi }_{k}N(x_{ij}|u_i{v}_j^T,{\sigma }_k^2)}}$$

M step 1 (updating Π, Σ): To estimate the MOG parameters, including Π and Σ, Eq. (5) is reformulated as follows:

(7) $$p(X,Z|U,V,\mathrm{\Pi },\mathrm{\Sigma })={\prod }_{i,j}{\sum }_{k=1}^K{\left[{\pi }_{k}N(x_{ij}|u_i{v}_j^T,{\sigma }_k^2)\right]}^{z_{ijk}}$$

Recall that $z_{ijk}$ is equal to one if $x_{ij}$ is generated by the $k^{th}$ Gaussian model; otherwise, it is zero.

Now, as shown in Eqs. (8) and (9), through maximizing the logarithm of Eq. (7), the parameters Π and Σ can be estimated, as presented in Eq. (10).

(8) $$L\left(X,Z|U,V,\mathrm{\Pi },\mathrm{\Sigma }\right)=\sum _{i,j}{\sum }_{k=1}^K{z}_{ijk}(\mathit{log}{\pi }_k-\mathit{log}\sqrt{2\pi }{\sigma }_k-{\displaystyle \frac{{(x_{ij}-u_i{v}_j^T)}^2}{2{\sigma }_k^2})}$$

$$E_{z}L\left(X,Z|U,V,\mathrm{\Pi },\mathrm{\Sigma }\right)=\sum _{i,j}{\sum }_{k=1}^{K}E(z_{ijk})(\mathit{log}{\pi }_k-\mathit{log}\sqrt{2\pi }{\sigma }_k-{\displaystyle \frac{{(x_{ij}-u_i{v}_j^T)}^2}{2{\sigma }_k^2}}$$

(9) $$=\sum _{i,j}{\sum }_{k=1}^K{\gamma }_{ijk}(\mathit{log}{\pi }_k-\mathit{log}\sqrt{2\pi }{\sigma }_k-{\displaystyle \frac{{(x_{ij}-u_i{v}_j^T)}^2}{2{\sigma }_k^2})}$$

(10) $${\pi }_k={\displaystyle \frac{N_k}{{\sum }_{k=1}^K{N}_k},{\sigma }_k^2={\displaystyle \frac{1}{N_k}{\sum }_{i,j}{\gamma }_{ijk}{\left(x_{ij}-u_i{v}_j^T\right)}^2}}$$where $N_k={\sum }_{i,j}{\gamma }_{ijk}$.

M step 2 (updating U, V): Employing the acquired estimations for the parameters Π and Σ, the values of U and V are determined through the employment of the Weighted Alternating Least Squares (WALS) method (Pan et al., 2008), which is described in the subsequent subsection.

Hence, the general procedure for estimating the low-rank matrices U and V involves the initial random initialization of the Π, Σ, and V matrices. Subsequently, by employing the EM steps and updating the U and V matrices using Eqs. (15) and (16) respectively, the process continues until the convergence criterion is satisfied (Pan et al., 2008).

Weighted alternating least squares (WALS)

Consider $X=\left[x_1,x_2,\dots ,x_n\right]\in R^{d\times n}$ as the data matrix, and $W=\left[w_1,w_2,\dots ,w_n\right]\in R^{d\times n}$ as the matrix of non-negative weights. These weights are computed based on the estimated noise values for each pixel, and they are formulated as $W=\sqrt{{\sum }_{k=1}^K{\displaystyle \frac{{\gamma }_{ijk}}{2\pi {\sigma }_k^2}}}$. A higher weight $w_{ij}$ indicates that the data $x_{ij}$ has a more influence on the generation of the restored image.

As demonstrated in the aforementioned equation, a greater value of ${\gamma }_{ijk}$ implies an increased likelihood of partial membership of the data $x_{ij}$ within the $k^{th}$ Gaussian model. Consequently, the weight assigned to this data also rises. The weighted low-rank matrix approximation of matrix X is computed by minimizing the subsequent function:

(11) $$L\left(Y\right)=\left|\right|W\odot \left(X-Y\right)|{|}_F^2$$

where $\left|\right|.|{|}_F^2$ signifies the Frobenius norm and $\odot$ represents the Hadamard product. Let’s assume that the matrix Y, comprising turbulence-free images, can be decomposed into two lower-rank U and V matrices ( $r\ll min\left(n,d\right)$): $Y=U{V}^T$. Consequently, Eq. (11) can be reformulated as follows:

(12) $$L\left(U,V\right)={\displaystyle {\sum }_{ij}{\left|\left|w_{ij}\odot \left(X_{ij}-U_i{V}_j^T\right)\right|\right|}_F^2}$$

To avoid overfitting, by adding new terms, Eq. (12) becomes the following:

(13) $$L\left(U,V\right)={\displaystyle {\sum }_{ij}{\left|\left|w_{ij}\odot \left(X_{ij}-U_i{V}_j^T\right)\right|\right|}_F^2}+\lambda \left({\left|\left|w_{ij}\odot U_i\right|\right|}_F^2+{\left|\left|w_{ij}\odot V_j\right|\right|}_F^2\right)$$

Here, λ represent regularization parameter. By assuming V as a constant and solving the equation ${\displaystyle \frac{1}{2}{\displaystyle \frac{\mathrm{\partial }L\left(U,V\right)}{\mathrm{\partial }{U}_i}=0}}$, the estimation of U is accomplished, as illustrated in Eq. (15):

(14) $${\displaystyle \frac{1}{2}{\displaystyle \frac{\mathrm{\partial }L\left(U,V\right)}{\mathrm{\partial }{U}_i}=U_i\left(V^{T}diag\left(w_i\right)V+\lambda \left({\sum }_j{w}_{ij}\right)I\right)-x_{i}diag\left(w_i\right)V}}$$

(15) $$u_i=x_{i}diag\left(w_i\right)V{\left(V^{T}diag\left(w_i\right)V+\lambda \left({\sum }_j{w}_{ij}\right)I\right)}^{-1},1\le i\le d$$

Here $diag\left(w_i\right)\in R^{n\times n}$ denotes a diagonal matrix with elements $w_i$ and I represents an identity matrix with dimensions $r\times r$.

By performing the same procedure, the value of V can also be calculated as follows:

(16) $$v_j=x_j^{T}diag\left(w_j\right)U{\left(U^{T}diag\left(w_j\right)U+\lambda \left({\sum }_i{w}_{ij}\right)I\right)}^{-1},1\le j\le n$$

Algorithm 1 outlines the pseudo-code for parameter estimation using the LRMF method.

Since the count of Gaussian models has a significant impact on both the convergence speed and the accuracy of the algorithm, the determination of the count for these models was conducted in an adaptive manner.

As stated in Meng & De La Torre (2013), the process commences with the initialization of the count of Gaussian models, set at a value suitable for modeling the noise distribution (denoted as $Ma{x}_k$). Subsequent to each iteration of the EM algorithm, the distance between each pair of Gaussian distributions is computed employing the Kullback-Leibler (KL) divergence technique. If the distance between the two Gaussian models falls below the pre-established threshold, these models are regarded as similar, prompting their amalgamation. Given that the mean value of the Gaussian distributions utilized for noise modeling is consistently taken as zero, the amalgamation of Gaussian models necessitates the computation solely of the standard deviation and the mixing proportion for the resultant hybrid Gaussian model.

Distortion model

Alongside inducing spatiotemporal varying blur, turbulence also introduces image distortion. Distortion encompasses abrupt shifts in the signal that lead to its degradation. Distortions in turbulent images emerge when the pixel displacement within the local vicinity can be modeled as a geometric transformation. Determining the optimal alignment for each pixel reduces distortion and strives to minimize the rank of the transformed image. With regard to the transformed image, the optimization equation for LRMF is presented as follows:

(17) $$\underset{U,V}{min}{\left|\left|W\odot \left(x\circ \tau -U{V}^T\right)\right|\right|}_F$$

As $x\circ \tau$ involves a nonlinear operation, it is very difficult to optimize it within the same formulation. However, considering the slight changes in τ, by linearizing (Eq. (17)), the iterative estimation of τ becomes feasible, as delineated in Eq. (18).

(18) $$X\circ \left(\tau +\mathrm{\Delta }\tau \right)\simeq X\circ \tau +{\sum }_{i=1}^n{J}_i\mathrm{\Delta }{\tau }_i$$where $J_i$ represents the Jacobian matrix corresponding to the $i^{th}$ image with respect to the transformation parameters ${\tau }_i$:

(19) $${\mathrm{J}}_{\mathrm{i}}={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\epsilon }({\mathrm{x}}_{\mathrm{i}}\circ \epsilon ){/}_{\epsilon ={\tau }_{\mathrm{i}}}}$$

The Robust Alignment by Sparse and Low-rank decomposition (RASL) method (Peng et al., 2012) has demonstrated the robustness of this approximation in cases of linearly correlated images, even in the presence of drastic variations in brightness and occlusion. Considering Eq. (19), the modification of the optimization function (Eq. (17)) for the $i^{th}$ image transpires as follows:

(20) $$\underset{U,V}{min}{\left|\left|W_i\odot \left(x_i\circ {\tau }_i+J_i\Delta {\tau }_i-U{V}_i^T\right)\right|\right|}_F$$where $V_i,\mathrm{a}\mathrm{n}\mathrm{d}{W}_i$ denote the $i^{th}$ columns of matrices V and W, respectively, pertaining to the $i^{th}$ image ( $x_i$). With the known values of U and V, the cost function is established in Eq. (21) to estimate the transformation matrix:

(21) $$\underset{\Delta \tau }{min}{\left|\left|e\right|\right|}_1,e=x_i\circ {\tau }_i+J_i\Delta {\tau }_i-U{V}_i^T$$

As Eq. (21) is a variation of the least absolute deviation problem, it can be effectively solved using the Alternating Direction Method of Multipliers (ADMM) method (Boyd, Parikh & Chu, 2011).

Subsequently, by formulating the Lagrangian function as Eq. (22) and incorporating the input values of $U,V_i,$ ${\tau }_{\mathrm{i}},$ and $J_i$, the optimal value of $\mathrm{\Delta }{\tau }^{P+1}$ is computed utilizing the ADMM approach by Eq. (23):

(22) $$L\left(U,V,e,\Delta \tau ,\lambda \right)={\left|\left|e\right|\right|}_1+{\lambda }^{T}h\left(e,\Delta \tau \right)+\frac{\mu }{2}{\left|\left|h\left(e,\Delta \tau \right)\right|\right|}_2^2$$where $h\left(e,\mathrm{\Delta }\tau \right)=U{V}_i^T+e-x_i\circ {\tau }_i-J_i\mathrm{\Delta }{\tau }_i$.

(23) $$\mathrm{\Delta }{\tau }^{P+1}={\left(J_i{J}_i^T\right)}^{-1}{J}_i^T\left(U{V}_i^T+e^P-x_i\circ {\tau }_i+{\displaystyle \frac{1}{{\mu }^P}{\lambda }^P}\right)$$where $e^{P+1}=S_{{\displaystyle \frac{1}{{\mu }^P}}}\left(x_i\circ {\tau }_i+J_i\mathrm{\Delta }{\tau }_i^{P+1}-U{V}_i^T-{\displaystyle \frac{1}{{\mu }^P}{\lambda }^P}\right)$, ${\lambda }^{P+1}={\lambda }^P+{\mu }^{P}h\left(e^{P+1},\mathrm{\Delta }{\tau }_i^{P+1}\right)$, and ${\mu }^{P+1}=\rho {\mu }^P$.

In these equations, $S_{{\displaystyle \frac{1}{{\mu }^P}}}$ stands for the elementwise soft thresholding operator (as described in Boyd, Boyd & Vandenberghe (2004)). Here, $\rho >1$ signifies the penalty value that increases the μ value in the ascending order, and P denotes the step within the ADMM algorithm. The overall procedure of this estimation is depicted in Algorithm 2.

As previously mentioned, the task entails determining the locally optimal transformation matrix for every pixel. Thus, for each image pixel, a 3 × 3 transformation matrix needs to be identified to effectively model the distortions. Rather than focusing on the transformation of individual pixels, it is feasible to model the transformation of a window encompassing the pixel at its center. After the application of the transformation matrix to each window, the value of each transformed pixel ( ${\mathrm{x}}_{\mathrm{i}}\circ {\tau }_{\mathrm{i}}$) within the $i^{th}$ image corresponds to the median value across all windows containing that pixel.

To achieve this goal, the sliding windows approach is employed. The dimensions of each window is adaptive depending on the entropy of the window. Initially, the window’s size is set to the minimum size $W_{min}$ and it expands until its entropy becomes lower than $\alpha$ percent ( $0<\alpha <1$) of the image’s entropy. The maximum size allowed for windows is $W_{max}$. $W_{min},$ $W_{max}$ and $\alpha$ constitute the free parameters of the proposed method. Since higher entropy in the image indicates more intricate details, lower entropy (signifying image smoothness) results in greater uncertainty when calculating the transformation matrix. Thus, if the window’s size is expanded to the maximum permissible extent without surpassing α percent of the overall image entropy, the calculation of the transformation matrix for that window is skipped. Instead, its value is approximated using the transformation matrix from neighboring windows. Despite distortion implying abrupt shifts in pixel positions, these changes cannot be very drastic. Consequently, a mechanism must be considered for removing transformation matrices that yield substantial pixel displacements. To address this, the mean and variance of all estimated $\mathrm{\Delta }{\tau }_i$ matrices are calculated. Those matrices exceeding a threshold of β percent, in terms of their distance from the mean value, relative to the total variance, are designated as outliers. Subsequently, pixels with the $\mathrm{\Delta }{\tau }_i$ matrix that are identified as outliers (pixels with a window size equal to $W_{max}$ and those with a $\mathrm{\Delta }{\tau }_i$ matrix with a drastic value) are assigned the mean value of $\mathrm{\Delta }{\tau }_i$ matrices from neighboring pixels. Ultimately, the transformation matrix for each window is determined by accumulating the estimated $\mathrm{\Delta }{\tau }_i$ values from different iterations. ( ${\tau }_i^{k+1}=\mathrm{\Delta }{\tau }_i+{\tau }_i^k$).

Given that the distortion affecting pixels arises from a limited range of distinct transformations, a MOG model is computed using the estimated transformation matrices. This MOG model serves as a representation of all the distortions affecting image pixels. Employing the computed MOG, the final transformation matrix for each window is established as the mean value of the Gaussian model, with the transformation matrix of that window being a member of said Gaussian model.

As a reminder, similar to the adaptive determination of the number of Gaussian models for noise (as outlined in Algorithm 1), the number of Gaussian models for distortions is also computed adaptively using the same approach. The comprehensive process of distortion removal is presented in Algorithm 3.

Selecting a subset of images using the RANdom SAmple Consensus (RANSAC) method

Turbulent images can exhibit varying degrees of distortion and blurring due to different conditions and intensities. Consequently, identifying and excluding images that exhibit distinct turbulence patterns can lead to improved outcomes. The RANSAC method is employed within the proposed approach to achieve this objective (refer to Algorithm 4). To accomplish this, a subset of m images $\left(m\ll n\right)$) is randomly selected, and the proposed method (as described in Algorithm 5) is applied to this subset, yielding a restored image (denoted as r1). Subsequently, each of the remaining images is processed using Algorithm 5 along with the initial set of m images. This leads to the generation of an output matrix $R=U{V}^T,$ containing $m+1$ images. Instead of averaging the images, the ${\left(m+1\right)}^{th}$ image is selected as the final image for the process.

By utilizing cross-correlation, the similarity between the restored image (r1) and the ${\left(m+1\right)}^{th}$ image is calculated, denoted as C1. Similarly, the distances between two noise distributions of the m images and the $m+1$ images and two transformation matrix distributions $\left(\tau \right)$of the m images and the $m+1$ images are computed using the Kullback-Leibler (KL) divergence. These are respectively termed C2 and C3.

If the weighted average of these three criteria (C1, C2, and C3) falls below a predetermined threshold, this image is regarded as a member of the support set of the initial m images. The procedure for calculating the KL divergence between two MOG distributions is elaborated in Appendix A.

This process is iterated multiple times, and in the end, the m images with the largest support set are selected. The selected $m$ images and their associated support set are combined to create the new subset. The final restored image is obtained by applying the proposed method to this refined subset. In this process, the excluded images are the ones that have different turbulence characteristics compared to the others. It is important to note that the excluded images may not necessarily have lower turbulence than those that were not selected. Nonetheless, due to the fact that the selected images have the same conditions compared to each other, empirical evidence indicates that the proposed method efficiently estimates the blur and distortion characteristics within this subset, leading to an improved output.

Overall procedure

As stated, the turbulence present in the images is modeled as being influenced by two distinct factors: distortion, characterized by geometric transformations that alter pixel positions, and spatiotemporal varying blur, which leads to variations in pixel brightness values.

Taking into account the methodologies discussed in the LRMF and distortion model sections, the comprehensive procedure of the proposed approach for removing turbulence in a collection of images captured from a single scene can be outlined as follows:

To begin, employ the RANSAC method as detailed in the selecting a subset of images using the RANSAC method section to identify the most suitable subset of input images. Then calculate the R matrix using optimization described by Eq. (17), which is the main optimization function of the proposed method, where $R=U{V}^T$ is the set of final images. Since the proposed procedure effectively eliminates turbulence, the R matrix includes n images by minimal turbulence. Consequently, a straightforward averaging of these n images yields the turbulence-free image (reference image). The estimation of both the U and V matrices entails an iterative methodology based on the process elucidated in the LRMF section. Nonetheless, due to the inherently nonlinear character of calculating the $x\circ \tau$ transformation, as expounded in the distortion model section, the transformation is executed iteratively by estimating ∆τ values. Consequently, to integrate these two iterative procedures, the previously outlined steps are executed concurrently. At the initiation of the procedure, all pixel transformation matrices are initially set as identity matrices ( $\mathrm{\forall }i{\tau }_i=\mathrm{I}$). Subsequently, the algorithm commences with an iteration of Algorithm 1 (as detailed in the LRMF section) to estimate the values of $V^k\mathrm{a}\mathrm{n}\mathrm{d}{U}^k$ (with k representing the current iteration step of the algorithm). To achieve this, the initial step involves estimating the parameters of the MOG. Following this, the weights needed for the WALS method are calculated based on these parameters. Subsequently, utilizing the calculated $V^k\mathrm{a}\mathrm{n}\mathrm{d}{U}^k$ values obtained through the WALS method (as described in the LRMF section), the second step is executed. This involves computing ∆τ for all images in accordance with Algorithm 2 (elaborated upon in the distortion model section), and subsequently generating $x_i\circ {\tau }_i^k$ (the transformed image in the $k^{th}$step). The algorithm’s detailed steps are presented in Algorithm 3. After this, the previously described procedure is repeated, and the U and V matrices are updated using the calculated $x_i\circ {\tau }_i^k$ values. This iterative process continues until the predetermined stopping condition is met. The complete procedure is visually depicted in Fig. 2. The main part of the proposed method is the calculation of the low rank matrix, which is shown by the green box in the mentioned figure. A detailed algorithmic representation of this part is provided in Algorithm 5.

Datasets

The evaluation of the proposed method encompasses both simulated and real datasets. The real datasets consist of the Chimney, Building, Moon-surface, and Water-tower sets, each of which includes corresponding ground truth images (Zhu & Milanfar, 2012). Meanwhile, the simulated datasets encompass the Car-front (Lau, Lai & Lui, 2019a) and Road (Lau, Lai & Lui, 2019a) sets. The simulated dataset is generated using two distinct approaches. The first method involves warping through a Gaussian deformation vector and introducing blurring through Gaussian noise. This dataset is referred to as the Gaussian simulated dataset. The second approach employs inter-modal and spatially correlated Zernike coefficients (Chimitt & Chan, 2020). This dataset is referred to as the Zernike simulated dataset. Further comprehensive information about these datasets is available in Appendix B.

Criterion for comparison of algorithms

In this article, the two criteria of PSNR (peak signal-to-noise ratio) and SSIM (structural similarity index) are used to evaluate the proposed method. The PSNR criterion uses the root mean square error between two images. Larger PSNR values indicate greater similarity between images. Unlike PSNR, which uses error estimation, SSIM is a perception-based approach that incorporates important perceptual information such as structure and contrast.

The value of SSIM is in the range of [0, 1] and equals one if the two images are exactly equal.

Free parameters and implementation details

The proposed method involves thirteen free parameters, which can be fine-tuned by experts and users who utilize the method. The specific values of these parameters, employed in all experiments within this section, are presented in Table 1. The free parameters can be classified into six distinct groups. The first group encompasses the convergence criteria of the EM method. This parameter involves a trade-off between runtime and algorithm accuracy. Given that turbulence removal is not a real-time computation but rather an offline problem, the chosen value strikes a balance between achieving satisfactory accuracy and ensuring a reasonable runtime. The intention is to avoid an excessively stringent threshold that could lead to impractical computation times while still maintaining attainable levels of accuracy. The second group of parameters pertains to the calculation of $\mathrm{\Delta }\tau$ matrices and includes $\alpha$, $W_{min},$ and $W_{max},$, which are utilized in determining the size of the sliding window. Additionally, β is employed to identify outlier ∆τ. These parameter values were chosen through a process of trial and error, involving thorough examination of the datasets. It is important to note that the value of $W_{min}$ should be chosen small enough to encompass local transformation information, while $W_{max}$ should not be too large, as it might include a substantial portion of the image, causing the estimated transformation to become out of locality. The third group of free parameters is used to determine the optimal value of ∆τ through the ADMM method. This group comprises parameters $\rho \mathrm{a}\mathrm{n}\mathrm{d}{\epsilon }^{tolerance}$ with their values set to the ones according to the recommendations of He et al. (2014). The fourth group of free parameters pertains to the number of MOG models used in noise and distortion estimation. This group encompasses the maximum number of models and the threshold for model integration, determined through a trial-and-error process. When selecting the value of the integration threshold, ensure that it is not too small, which would preserve similar models, and not excessively large, as it might lead to the merging of dissimilar models. The fifth group of free parameters encompasses the RANSAC parameters: the parameter m represents the minimum number of images required for the algorithm to function accurately. Through the examination of the datasets, it has been determined that the proposed method can effectively remove turbulence with a set of 25 images. Therefore, the parameter m has been set to 25. Within this group, a threshold value is defined for the purpose of identifying the support set. If the weighted mean of the three criteria introduced in selecting a subset of images using the RANSAC method section becomes lower than this threshold, the current image is regarded as a member of the support set of the initial set of m images. The value of this threshold is also determined through trial and error. The stopping criteria for the EM method (Algorithm 1), RANSAC method (Algorithm 4), and calculating low rank matrix (Algorithm 5) are the number of iterations, which are set to 40, 10, and 100 iterations, respectively. The sixth group of free parameters involves the λ parameter applied in the WALS method. The determination of this parameter value was achieved through an iterative process of trial and error, conducted specifically on the applied datasets.

Compare the results of the proposed method with the existing methods

As mentioned, the results of the proposed method are compared with the results of the other eleven methods. The results of the Gaussian simulated dataset are compared with the Two-Stage, NDL, SGL, Centroid, RPCA, and JSUB methods, and the results of the Chimney and Building image sequences are compared with the PCA-Based, Lucky Region, Two-Stage, BNLTV, NDL, HVDKR, SGL, Centroid, and UCSPF methods. These comparative results are presented in Tables 2 and 3, respectively. The results of the proposed method on the sequences of the Moon-surface and Water-tower images are provided in Appendix C. As none of the previous methods reported results on these two sequences and their code was unavailable, only the results of the proposed method are presented for these datasets. Similarly, Table 4 contains results only for the Zernike simulated dataset. According to the results presented in Tables 2 and 3, it can be observed that the proposed method outperformed the previous methods in terms of the SSIM criterion for all datasets. However, in terms of the PSNR criterion, the proposed method exhibited lower performance compared to the previous methods for the Road and Chimney datasets.

Figures 3–6 provide visual representations of the restored images obtained from the proposed method and the results of the previous methods. It is evident from these images that the Centroid, RPCA, JSUB, Two-Stage, NDL, HVDKR, and UCSPF methods yield higher visual quality in terms of the final restored image compared to the other methods.

The Centroid method exhibits a restored image with lower quality compared to other methods. Despite the restoration of the geometric structure, the resulting image appears blurry and overly brightened due to the effects of temporal averaging.

The RPCA method, lacking a mechanism to effectively remove distortion, primarily focuses on selecting a subset of turbulent images with minimal distortion and maximal sharpness during the subsampling stage. As a result, the restored image obtained using the RPCA method continues to exhibit some degree of distortion. Nonetheless, due to its emphasis on the selection of a subset of sharp images and using the appropriate energy function to merge these subsampled images, coupled with the application of the blind deconvolution method, the restored image has little blur.

The JSUB method closely resembles the RPCA method except that its energy function uses TV regularization, the $L_2$ norm, and adaptive Gaussian noise for joint subsampling. This method also lacks a blind deconvolution step, resulting in a mildly distorted but blurred image.

The Two-Stage method’s deficiency caused by from the generation of an unsuitable reference image through the averaging process, leading to errors during the registration stage and ultimately resulting in a blurred final image.

Similar to the Two-Stage method, the NDL method also suffers from the impact of the erroneous reference image generated through the averaging process. This issue substantially affects the performance of the registration stage, and even though the method employs the blind deconvolution technique, the resulting image lacks clarity in terms of capturing image details. HVDKR uses the RPCA method to generate a reference image with reduced blur and distortion. This approach then further enhances the reference image using a variational model. Employing deformation-guided fusion the method not only models and mitigates the impact of the image but also effectively reduces blurriness in the final restored image. But since it is a step-by-step process, the interaction between the two factors of blur and distortion has not been taken into account. As a result, while the output images may exhibit relatively improved sharpness, the overall restored image tends to retain some level of blurriness and distortion, and also the fine details in the image are not well-defined.

In the UCSPF method, by separating the phase and magnitude components in the frequency domain, structural information and brightness intensity are separated to process the two distortion and blur factors separately. Through multiscale and multidimensional operations applied to these distinct components, the UCSPF method effectively produces an output image characterized by reduced levels of distortion and blur. Moreover, this approach facilitates the preservation of sufficient image details. Nonetheless, the processing conducted in the frequency domain lacks access to spatial information. For this reason, frequency-domain processing, although well-restored in detail in high-energy areas, has caused fluctuations in the brightness intensity in smooth regions. Despite the visually satisfactory quality of the final images, it is noteworthy that the computed PSNR and SSIM values are lower than expected.

As evident in Figs. 3–6, the resulting images from the proposed method exhibit sharper edges and finer details in comparison to those generated by other methods. The smooth regions within the ground truth images are returned as smoothly in the output images. The distortion is removed as much as possible, but the brightness of the output image in the Chimney sequence produced by the proposed method appears slightly darker compared to the corresponding ground truth image. Figure 7 presents the restored images generated by the proposed method on the Zernike turbulence dataset. A visual comparison of both Gaussian and Zernike turbulence images reveals that blurring is more obvious than distortion in the Zernike images. While the proposed method effectively restores Zernike turbulence images with notable accuracy, as evident in Fig. 7F, it is worth noting that the details of the asphalt road are not fully restored to their original quality.

Indeed, conventional images typically exhibit characteristics of low-frequency signals, rendering them low rank. Distortion, characterized by sudden change in the signal changes, introduces high-frequency components that elevate the signal’s rank.

The time-varying blur in an image sequence from a single scene results in images that are distinct from one another. Consequently, the overall diversity of the image sequence increases. This phenomenon elevates their temporal rank. Turbulence images, caused by both distortion and spatiotemporal varying blur, possess high ranks. Consequently, by generating a low-rank matrix, it becomes possible to obtain an image that encompasses all high-rank images derived from it. This image aims to closely resemble the reference (turbulence-free) image. Consequently, the creation of a low-rank matrix can be used to eliminate turbulence. Accordingly, the LRMF method was adopted for this particular objective. Through the iterative integration of these two techniques within the LRMF method, a comprehensive framework has been established. This framework is capable of generating turbulence-free images while minimizing the presence of turbulence. The performance of the proposed method was assessed in terms of the PSNR and SSIM criteria using the t-test statistical method. To accomplish this, the maximum values obtained from the previous methods were compared with those of the proposed method. The outcome of the PSNR-based t-test indicated a 50% similarity, signifying that the proposed method’s accuracy is on par with the best outcomes of the previous methods. Conversely, the t-test conducted based on the SSIM yielded a result of 95%, underscoring the superior accuracy of the proposed method in comparison to the preceding previous methods. Since the PSNR criterion is based on MSE, even slight changes in the image’s brightness can lead to a decrease in the numerical values of this criterion. Therefore, the relative darkening observed in the results of the proposed method might contribute to the decrease in the numerical values of this criterion. On the other hand, SSIM is a perception-based criterion that, like the human visual system, considers structural, luminance, and contrast similarities between the two images. Therefore, higher numerical values of this criterion in the results of the proposed method indicate a high perceptual similarity between the resultant images and the reference images. Further detailed evaluation of the proposed method is provided in Appendix D.