BM3D Frames and Variational Image Deblurring

I Introduction

We consider image restoration from a blurry and noisy observation. Assuming a circular shift-invariant blur operator and additive zero-mean white Gaussian noise the conventional observation model is expressed as

where $||\mathbf{\cdot||}_{2}$ stands for the Euclidean norm, $pen(\cdot)$ is a penalty functional and $\tau>0$ is a regularization parameter.

Image modeling lies at the core of image reconstruction problems. Recent trends are concentrated on sparse representation techniques, where the image is assumed to be defined as a combination of few atomic functions taken from a certain dictionary. It follows that the image can be parameterized and approximated locally or nonlocally by these functions. To enable sparse approximations, the dictionary should be rich enough to grasp all variety of the images. Clearly, bases are too limited for this task and one needs to consider overcomplete systems with a number of elements essentially larger than the dimensionality of the approximated images. Frames are generalization of the concept of basis to the case when the atomic functions are linearly dependent and form an overcomplete system . There is a vast amount of literature devoted to the sparsity based models and methods for imaging. An excellent introduction and overview of this area can be found in the recent book .

The contribution of this paper concerns three main aspects of image deblurring: image modeling, variational problem formulation, and algorithmic reconstruction.

First, the BM3D image modeling developed in is formalized in terms of the overcomplete sparse frame representation. We construct analysis and synthesis BM3D-frames and study their properties. The analysis and synthesis developed in BM3D are interpreted as a general sparse image modeling applicable to variational formulations of various image processing problems.

Second, we consider two different formulations of the image deblurring problem: one given by minimization of the objective function and another based on the Nash equilibrium. The latter approach results in an algorithm where the denoising and the deblurring operations are decoupled.

Third, it is shown by simulation experiments that the best image reconstruction both visually and numerically is obtained by the algorithm based on decoupling of blur inverse and noise filtering. To the best of our knowledge, this algorithm provides results which are the state-of-art in the field.

Here we extend and develop our preliminary ideas sketched in . The BM3D frames are now constructed explicitly, taking into account the particular form of the 3D transform. Proofs of the frame properties are presented. We develop algorithms for the analysis and synthesis-based problem formulations introduced in and provide their convergence analysis. The problem formulation based on the Nash equilibrium and the corresponding decoupled deblurring algorithm are novel developments.

The paper is organized as follows. We start from a presentation of the BM3D image modeling and introduce BM3D-frames (Section II). The variational image reconstruction is a subject of Section III. The algorithms based on the analysis and synthesis formulations are derived in this section. The algorithm based on the Nash equilibrium is presented in Section IV. Convergence results for the proposed algorithms are given in Section V. Implementation of the algorithms is discussed in Section VI. The experiments and comparison of the algorithms are given in Section VII. In Section VIII we discuss the principal differences of the decoupled formulation compared to the analysis and synthesis formulations. Concluding remarks are done in the last section. Proofs of mathematical statements are given in Appendix.

II Overcomplete BM3D image modeling

BM3D is a nonlocal image modelling technique based on adaptive, high order groupwise models. Its detailed discussion can be found in . Below, using the example of the denoising algorithm , we recall the concept of the BM3D modeling. The denoising algorithm can be split into three steps.

Analysis. Similar image blocks are collected in groups. Blocks in each group are stacked together to form 3-D data arrays, which are decorrelated using an invertible 3D transform.

Processing. The obtained 3-D group spectra are filtered by hard thresholding.

Synthesis. The filtered spectra are inverted, providing estimates for each block in the group. These blockwise estimates are returned to their original positions and the final image reconstruction is calculated as a weighted average of all the obtained blockwise estimates.

The blocking imposes a localization of the image on small pieces where simpler models may fit the observations. It has been demonstrated that a higher sparsity of the signal representation and a lower complexity of the model can be achieved using joint 3D groupwise instead of 2D blockwise transforms. This joint 3D transform dramatically improves the effectiveness of image spectrum approximation.

The total number of groupwise spectrum elements is much larger than the image size, and we arrive to an overcomplete or redundant data approximation. This redundancy is important for effectiveness of the BM3D modeling.

Our target is to give a strict frame interpretation of the analysis and synthesis operations in BM3D.

The particular form of the 3-D decorrelating transform constitutes an important part of the BM3D modeling. It is constructed as a separable combination of 2-D intrablock and 1-D interblock transforms. The 2-D transform, in turn, is typically implemented as a separable combination of 1-D transforms. Let $\mathbf{D}_{2}$ and $\mathbf{D}_{1}$ be $\sqrt{N_{bl}}\times\sqrt{N_{bl}}$ and $K\times K$ size matrices representing respectively 1-D interblock and 1-D intrablock transforms. Then the separable 2-D transform for the block $\mathbf{Y}_{j}$ is given by the formula

The vectorization of this formula using the Kronecker matrix product $\mathbf{\otimes}$ gives

Performing vectorization again, we express the 3-D group spectrum coefficients in a compact form:

where $\mathbf{\omega}_{r}$ is the columnwise vectorized matrix $\mathbf{\Omega}_{r}$ and $\mathbf{d}_{j}$ is the $j$-th column of $\mathbf{D}_{1}$. Finally, denoting

The matrix $\mathbf{\Phi}$ defined by the formulas (3)-(4) gives an explicit representation of the BM3D analysis operation.

The synthesis matrix is derived similarly. First, the inverse 3-D transform is applied to each group spectrum $\mathbf{\omega}_{r}$ and then obtained block estimates are returned to their original positions by $\mathbf{P}_{j}^{T},j\in J_{r}$. The estimate obtained from the $r$-th group spectrum is expressed as $\mathbf{\Psi}_{r}\mathbf{\omega}_{r}$, where

The final image estimate is defined as the weighted mean of the groupwise estimates using weights $g_{r}>0$. Hence the synthesis operation has the form

normalizes the weighted mean. $\mathbf{W}$ is a diagonal matrix, since all products $\mathbf{P}_{j}^{T}\mathbf{P}_{j}$ are diagonal matrices. The $m$-th diagonal element of $\mathbf{P}_{j}^{T}\mathbf{P}_{j}$ is $1$ if the $m$-th pixel of $\mathbf{y}$ belongs to the $j$-th block, otherwise it is . Thus, the $m$-th diagonal elements of the matrix-sum $\sum_{j\in I_{r}}\mathbf{P}_{j}^{T}\mathbf{P}_{j}$ indicates the number of blocks in the $r$-th group containing $m$-th pixel.

The matrix $\mathbf{\Psi}$ defined by the formulas (5)-(7) gives the matrix representation of the BM3D synthesis operation.

II-B Frame interpretation

The following equations hold for the matrices $\mathbf{\Phi}$ and $\mathbf{\Psi}$ defined by (4) and (6):

The frame $\left\{\mathbf{\phi}_{n}\right\}$ is not tight because $a\neq b$. This follows from the fact that the elements on the diagonal of matrix $\sum_{r}\sum_{j\in I_{r}}\mathbf{P}_{j}^{T}\mathbf{P}_{j}$ count the number of blocks containing a given pixel. These values are different for different pixels, since pixels from the blocks possessing higher similarity to other blocks participate in a larger number of groups.

Similarly, using (9) we can show that columns of $\mathbf{\Psi}$ constitute a non-tight frame $\left\{\mathbf{\psi}_{n}\right\}$. From equation (10) it follows that $\left\{\mathbf{\phi}_{n}\right\}$ is dual to $\left\{\mathbf{\psi}_{n}\right\}$. In general $\left\{\mathbf{\phi}_{n}\right\}$ is an alternative dual and becomes canonical dual only when all weights $g_{r}$ are equal.

We would like to emphasize that since groups and weights are selected data adaptively, the constructed frames are also data adaptive.

The presented frame interpretation allows to extend the scope of the BM3D modeling to the modern variational image reconstruction techniques.

III Variational image deblurring

The frame based variational image reconstruction problem allows two different formulations depending on what kind of image modeling, analysis or synthesis is used . In the analysis formulation the relation between the image and spectrum variables is given by the analysis equation $\mathbf{\omega}=\mathbf{\Phi y}$. The problem is formalized as a constrained optimization:

where $\left\|\cdot\right\|_{p}$ is the standard notation of the $l_{p}$-norm.

In the synthesis formulation the relation is given by the synthesis equation $\mathbf{y=\Psi\omega},$ leading to the constrained optimization:

These problems have equivalent unconstrained forms in which they usually encounter in literature. To obtain them it is enough to eliminate $\mathbf{\omega}$ and $\mathbf{y}$ respectively from (13) and (14). The analysis problem is then formulated as the minimization in the image domain

Similarly, the synthesis problem is formulated as the minimization in the spectrum domain

Despite of the algebraic similarity, the analysis and synthesis formulations generally lead to different solutions. A detailed discussion of the nontrivial connections between the analysis and synthesis formulations can be found in .

The problems (13)-(16) and the corresponding solution techniques recently become a subject of an intensive study. In particular, several algorithms have been suggested for the convex $l_{1}$-norm penalty. These algorithms sharing many common ideas are known under different names such as split Bregman iterations , iterative shrinkage algorithms , alternating direction method of multipliers , majorization-minimization algorithms . In this paper similar to we confine ourself to the Augmented Langrangian (AL) technique, using it as a simple and efficient tool for an explicit derivation of the reconstruction algorithms. This AL technique, introduced independently by Hestenes and Powell is now widely used for minimization of convex functionals under linear equality constraints.

The AL criterion for the analysis formulation (13) takes the form:

Finding the saddle point requires minimization of $L_{\text{a}}$ with respect to the variables $\mathbf{y},\mathbf{\omega}$ and maximization with respect to $\mathbf{\lambda}$. A common practical approach is to find the saddle point by performing alternating optimization. Applied to (17) it results in the following iterative scheme:

Here maximization with respect to $\mathbf{\lambda}$ is produced as a step (20) in the direction of the gradient $\nabla_{\mathbf{\lambda}}L_{\text{a}}$, with a step-size $\beta>0$. The convergence of the scheme (18)-(20) is studied in .

Minimization with respect to $\mathbf{y}$. Since $L_{\text{a}}$ is quadratic with respect to $\mathbf{y}$ the optimal solution is defined by the linear equation

We denote by $\hat{Y}_{\text{a}}\left(\mathbf{\omega,\mathbf{\lambda}}\right)$ the operator giving the solution of (21).

Minimization with respect to $\mathbf{\omega}$. Regrouping the terms in $L_{\text{a}}$ we arrive to the following formula

Since the first and the last terms do not depend on $\mathbf{\omega}$, the problem is reduced to the optimization

For $p\leq 1,$ the $l_{p}$-norm is non-differentiable which makes optimization on $\mathbf{\omega}$ non-trivial. Nevertheless, for $p=0$ and $p=1$ there are well known analytical solutions.

Let us denote $\mathbf{\mathbf{b}}=\mathbf{\Phi y-\mathbf{\lambda}}$, then (22) takes the form

Depending on the used norm the solution of (23) is given either by the hard or soft thresholding according to the formula:

Here all vector operations are elementwise, and ’$\circ$’ stands for the elementwise product of two vectors. We use $\mathfrak{Th}_{\tau}\left(\mathbf{b}\right)$ as a generic notation for the thresholding operator. Note, that for a given $\tau$ the thresholding levels for the hard and soft thresholdings are calculated differently.

Applying the general formula (26) to (22) we obtain the solution in the form

Following (18)-(20) and using (21) and (27) we define the analysis-based iterative algorithm which is presented in Figure 1. In each iteration it first updates the image estimate using the linear filtering (21). Then, the difference between the spectrum $\mathbf{\Phi y}_{t}$ and $\mathbf{\lambda}_{t}$ is thresholded, what corresponds to the optimization with respect to $\mathbf{\omega}$. Finally, the Lagrange multipliers are updated in the direction of the gradient $\mathbf{\omega}_{t+1}-\mathbf{\Phi y}_{t+1}$. Process is iterated until some convergence criteria is satisfied. Particularly, the iterations can be stopped as soon as the difference between consecutive estimates becomes small enough.

III-B Synthesis-based reconstruction

The AL criterion for the synthesis formulation (14) takes form:

In $L_{\text{s}}$, as opposed to $L_{\text{a}}$, the spectrum variable $\mathbf{\omega}$ enters the quadratic term with a matrix factor $\mathbf{\Psi}$. It makes the thresholding formula (26) inapplicable for minimizing $L_{\text{s}}\left(\mathbf{y},\mathbf{\omega},\mathbf{\lambda}\right)$ with respect to $\mathbf{\omega}$. One option is to apply one of the iterative shrinkage methods , but we prefer to follow a different approach which leads to a simpler solution. We modify (28) by introducing a splitting variable $\mathbf{u}\in R^{M},$ used as an auxiliary estimate of the spectrum $\omega$. The modified AL takes the form:

The corresponding saddle point problem is

where optimization with respect to the splitting variable $\mathbf{u}$ is required.

Minimization with respect to $\mathbf{y}$ is given by the solution of the linear equation

Minimization with respect to $\mathbf{u}$ satisfies the linear equation

Minimization with respect to $\mathbf{\omega}$, thanks to the splitting variable $\mathbf{u}$, can be obtained by the thresholding (26) with the parameter $\tau\xi$:

We denote by $\hat{Y}_{\text{s}}\left(\mathbf{u,\omega,\lambda}\right)$ and $\hat{U}_{\text{s}}\mathbf{(y,}\mathbf{\omega},\mathbf{\mathbf{\lambda})}$ the operators giving the solutions of (31) and (32).

Using (31)-(33) we define the synthesis-based iterative deblurring algorithm which is presented in Figure 2. At the first two steps the estimates for the image $\mathbf{y}_{t}$ and the splitting variable $\mathbf{u}_{t}$ are updated by solving (31) and (32). Then, the splitting variable $\mathbf{u}_{t+1}$ is thresholded reducing the complexity of the spectrum estimate $\mathbf{\omega}$. Finally, the Lagrange multipliers are updated in the direction of the gradient $\mathbf{y}_{t+1}\mathbf{-\Psi u}_{t+1}$. Process is iterated until some convergence criteria is satisfied.

IV Decoupling of blur inversion and denoising

Above we considered algorithms based on the minimization of a single objective function. In this section we present an alternative approach based on formulation of the deblurring as a Nash equlibrium problem for two objective functions. This approach allows to split the deblurring problem into two subproblems: a blur inversion and denoising, which are then solved sequentially. Such a decoupling has several advantages:

The decoupled algorithms are simpler in design and parameter selection;

The blur inversion can be implemented efficiently using Fast Fourier Transform (FFT);

Various denoising algorithms can be used in this scheme selected independently with respect to deblurring;

In many cases decoupled algorithms demonstrate better performance than the algorithms where deblurring and denoising are performed jointly.

Examples of the decoupled deblurring can be found in works , , and , where the regularized inverse is followed by different types of filtering (wavelet, shape-adaptive DCT, BM3D, pyramidal). An interesting development of this technique is demonstrated in where an iterative algorithm is derived by alternating optimization of multiple objective functions.

Let us formulate the deblurring problem as the following constrained optimization:

where $\varepsilon_{1},\varepsilon_{2}>0$. This problem can be replaced by the equivalent unconstrained one:

and $\gamma,\xi$ are constants selected correspondingly to the values of $\varepsilon_{1},\varepsilon_{2}.$

In terms of the game theory the problem (35) can be interpreted as a game of two players identified, respectively, with two variables $\mathbf{y}$ and $\mathbf{\omega}$ ,. An interaction between the players is noncooperative because minimization of $L_{\text{inv}}(\mathbf{y,\omega})$ with respect to $\mathbf{y}$ in general results in increase of $L_{\text{den}}(\mathbf{y,\omega})$ and minimization of $L_{\text{den}}(\mathbf{y,\omega})$ with respect to $\mathbf{\omega}$ increases $L_{\text{inv}}(\mathbf{y,\omega})$. The equilibrium of this game called Nash equilibrium defines the fixed point $\left(\mathbf{y}^{\ast},\mathbf{\omega}^{\ast}\right)$ of the optimization. For $p=1$, problem (35) is convex.

The objective functions $L_{\text{inv}}$ and $L_{\text{den}}$ allow the following interpretation. In $L_{\text{inv}}$ the fidelity term $\dfrac{1}{2\sigma^{2}}\left\|\mathbf{z}-\mathbf{Ay}\right\|_{2}^{2}$ evaluates the divergency between the observation $\mathbf{z}$ and its prediction $\mathbf{Ay}$. This fidelity is penalized by the norm $\left\|\mathbf{y-\Psi\omega}\right\|_{2}^{2}$ defining a difference between $\mathbf{y}$ and its prediction $\mathbf{\Psi\omega}$ through $\mathbf{\omega}$. The term $\dfrac{1}{2\xi}\left\|\mathbf{\omega}-\mathbf{\Phi y}\right\|_{2}^{2}$ in $L_{\text{den}}$ evaluates a difference between the spectrum $\mathbf{\omega}$ and the spectrum prediction $\mathbf{\Phi y}$ obtained from $\mathbf{y}$. The error between $\mathbf{\omega}$ and $\mathbf{\Phi y}$ is penalized by the norm $\left\|\mathbf{\omega}\right\|_{p}$.

Hence the Nash equilibrium provides a balance between the fit of the reconstruction $\mathbf{y}$ to the observation $\mathbf{z}$ and the complexity of the model $\left\|\mathbf{\omega}\right\|_{p}$. This can be contrasted with the analysis and synthesis-based problem formulations where the balance is provided within a single criterion. As we demonstrate later the form of the balance plays an essential role in the reconstructions with non-tight frames.

IV-B IDD-BM3D algorithm

To solve (35) we consider the following iterative procedure:

The iterative algorithm (38) models the selfish behavior, where each variable minimizes only its own objective function. These iterations converge to the fixed point $\left(\mathbf{y}^{\ast},\mathbf{\omega}^{\ast}\right)$ of (35), the corresponding result is formulated in Section V.

Minimization of $L_{\text{inv}}$ with respect to $\mathbf{y}$ is given by the solution of the linear equation

This step performs regularized inversion of the blur operator.

The minimization of $L_{\text{den}}$ with respect to $\mathbf{\omega}$ is obtained by thresholding with the threshold parameter $\tau\xi$:

Thus, in (38) the blur inversion and the denoising steps are fully decoupled.

The algorithm based on (38) is presented in Figure 3. We call this algorithm Iterative Decoupled Deblurring BM3D (IDD-BM3D).We wish to note that IDD-BM3D is similar but not identical to our Augmented Lagrangian BM3D deblurring (AL-BM3D-DEB) algorithm presented earlier in . The AL-BM3D-DEB algorithm is derived from the analysis-based formulation (17). The regularized inverse step (21) in AL-BM3D-DEB is replaced by the inverse (31) obtained from the synthesis-based formulation (28). In this replacement is treated as an approximation and is not mathematically rigorous. The presence of the Lagrange multipliers discriminates the AL-BM3D-DEB algorithm from the IDD-BM3D.

V Convergence

The main motivation of the AL technique is to replace a constrained optimization with a simpler saddle-point problem. The equivalence of these two problems is not a given fact. The classical results stating equivalence are formulated for the convex and differentiable functions . Since $l_{p}$-norms with $p\leq 1$ are non-differentiable these results are inapplicable. Nevertheless, for the $l_{1}$-norm the equivalence can be shown, provided that the constraints in the problem are linear. In the recent paper the equivalence statement is proved for the total variation penalty. This proof remains valid for any convex and non-differentiable penalties, in particularly for the $l_{1}$-norm based penalties. The equivalence result is formulated as following:

$(\mathbf{\hat{y},\hat{\omega}})$ is a solution of the analysis or synthesis problems if and only if there exist a saddle-point of the corresponding ALs.

Practically it means that the saddle-point of the AL optimization can be used in order to obtain the solutions of the considered optimization problems.

The convergence properties for the analysis and synthesis-based algorithms are formulated in the following proposition.

(a) If there exists a saddle point $\mathbf{(y}^{\ast}\mathbf{,\omega}^{\ast}\mathbf{,\lambda}^{\ast}\mathbf{)}$ of $L_{\text{a}}\left(\mathbf{y,\omega,\lambda}\right)$ (17), then $\mathbf{y}_{t}\mathbf{\rightarrow y}^{\ast}\mathbf{,\omega}_{t}\mathbf{\rightarrow\omega}^{\ast}\mathbf{,\lambda}_{t}\mathbf{\rightarrow\lambda}^{\ast}$.

(b) If there exists a saddle point $\mathbf{(y}^{\ast}\mathbf{,\omega}^{\ast}\mathbf{,u}^{\ast}\mathbf{,\lambda}^{\ast}\mathbf{)}$ of $L_{\text{s}}\left(\mathbf{y,\omega,u},\mathbf{\lambda}\right)$ (29), then $\mathbf{y}_{t}\mathbf{\rightarrow y}^{\ast}\mathbf{,\omega}_{t}\mathbf{\rightarrow\omega}^{\ast}\mathbf{,u}_{t}\mathbf{\rightarrow u}^{\ast}\mathbf{,\lambda}_{t}\mathbf{\rightarrow\lambda}^{\ast}$.

On the other hand, if no such saddle point exists, then at least one of the sequences \mathbf{\{y}_{t}\mathbf{\}}\or $\mathbf{\{\lambda}_{t}\mathbf{\}}$ must be unbounded.

V-B IDD-BM3D algorithm

For any set of parameters $\sigma,\tau,\gamma,\xi$ the sequence $\left(\mathbf{y}_{t},\mathbf{\omega}_{t}\right)$ generated by the IDD-BM3D algorithm with equal group weights $g_{r}$, converges to the fixed point $\left(\mathbf{y}^{\ast},\mathbf{\omega}^{\ast}\right)$ defined by the equations (35), if the fixed point exists.

The proof of the proposition is given in Appendix B. It is not required that the fixed point is unique. Depending on a starting point $(\mathbf{y}_{0},\mathbf{\omega}_{0})$ the limit point of the algorithm can be different but should satisfy the fixed point equations.

VI Implementation

Grouping and frame operators. To build the groups, we use the block-matching procedure from and apply it to the image reconstructed by the BM3DDEB deblurring algorithm . The found locations of the similar blocks constitute the set $J$ that is necessary to construct the analysis and synthesis frames. Multiplications against the matrices $\mathbf{\Phi,\Phi}^{T}\mathbf{,\Psi}$ and $\mathbf{\Psi}^{T}$ are calculated efficiently since all of them involve only groupwise separable 3-D transformations of the data (possibly with some averaging of the estimates). In our experiments the 3-D transform is performed by first applying the 2-D discrete sine transform (DST) to each block in the group followed by the 1-D Haar transform applied along the third dimension of the group. The image block size is $4\times 4$, and the number of blocks in the group is $8$.

Choice of the group weights. Since image blocks are overlapping, for each pixel we obtain several estimates. The weighted averaging can be used to improve the final aggregated estimate. For the one-step (non-iterative) algorithms the weights can be adaptively selected so to minimize the variance of the final aggregated estimate, based on the variance of each of the estimates (e.g. , , ). In the considered iterative algorithms the influence of the weights on the final estimate is complex, and deriving a formula for the optimal weights is rather involved. Instead, following the idea of the sparse representations, we suggest giving the preference to the estimates obtained from the sparser groups. In our implementations we use weights inversely proportional to the number of significant spectrum coefficients of the groups $g_{r}=1/\left\|\mathfrak{Th}_{\epsilon}\left(\mathbf{\omega}_{r}\right)\right\|_{0}$, where significant coefficients are found by the hard thresholding of the group spectra using a small threshold $\epsilon$.

The grouping and the adaptive group weights are calculated only once, using the initial image estimate $\mathbf{y}_{\text{init}}$ and remain unchanged through the subsequent iterations.

Choice of the regularization parameters. The parameters $\tau,\gamma,\xi$ are optimized to provide best reconstruction quality. Optimization has been performed separately for each algorithm and each deblurring scenario. The parameter $\beta$ is always set to 1.

Initialization. We experimentally confirmed the convergence to an asymptotic solution that is independent of the initialization $\mathbf{y}_{0}$ and $\mathbf{\omega}_{0}$. Nevertheless, initialization with a better estimate, for example with the reconstruction obtained by BM3DDEB (which we also use to define grouping) results in a much faster convergence.

Solution of the large-scale linear equations. All proposed algorithms contain steps involving solution of large-scale linear equations. For a circular shift-invariant blur operator, the solution of the equations (31) and (39) can be calculated in the Fourier domain using the FFT. The more complex equations (21) and (32) are solved using the conjugate gradient method. The conjugate gradient method allows avoiding explicit calculations of the matrices $\Phi^{T}\Phi$ and $\Psi^{T}\Psi$, since it requires only evaluating products of these matrices against vectors.

Practical considerations. The two steps of the IDD-BM3D algorithm can be merged into a single one

where the analysis-thresholding-synthesis operation $\Psi\mathfrak{Th}_{\tau\xi}\left(\mathbf{\Phi}\mathbf{y}_{t}\right)$ can be calculated groupwise without need to obtain the whole spectrum $\omega_{t}$ explicitly. Here $\mathbf{h}$ denotes the vectorized blurring kernel corresponding to the blur operator $\mathbf{A}$, and ’$\circ$’ stands for the elementwise product of two vectors. The operator $\mathcal{F}\left(\cdot\right)$ reshapes the input vector into a 2-D array, performs 2-D FFT and vectorizes the obtained result. $\mathcal{F}^{-1}\left(\cdot\right)$ works analogously, performing inverse FFT.

Complexity. Application of the frame operators is the most computationally expensive part of the proposed algorithms. However, due to their specific structure, the complexity of the frame operators $\mathbf{\Phi}$ and $\mathbf{\Psi}$ is growing only linearly with respect to the number of the pixels in the image. To give an estimate of the complexity of the IDD-BM3D algorithm, we mention that, on a $256\times 256$ image, one iteration takes about 0.35 seconds, and about 50 iterations are typically sufficient. This timing has been done on dual core 2.6 GHz processor for an implementation where the computationally most intensive parts have been written in C++.

VII Experiments

We consider six deblurring scenarios used as the benchmarks in many publications (e.g., and ). The blur point spread function (PSF) $h\left(x_{1},x_{2}\right)$ and the variance of the noise $\sigma^{2}$ for each scenario are summarized in Table I. PSFs are normalized so that $\sum h=1$. Each of the scenarios was tested with the four standard images: Cameraman, Lena, House and Barbara.

All three proposed algorithms, namely: analysis-based, synthesis-based and IDD-BM3D are evaluated in the scheme with the soft thresholding and unit group weights ($g_{r}=1$). Additionally, the IDD-BM3D algorithm is tested with the adaptive group weights ($g_{r}=1/\left\|\mathfrak{Th}_{\epsilon}\left(\mathbf{\omega}_{r}\right)\right\|_{0}$) using the soft and hard thresholdings.

In Table II we present improvement of signal-to-noise ratio (ISNR) values achieved by each algorithm for the Cameraman image. From these values we can conclude that the synthesis-based algorithm performs essentially worse than the IDD-BM3D algorithm, with the analysis-based algorithm being in-between. We can also see that the adaptive weights indeed provide a noticeable restoration improvement. Finally, comparing the last two rows, we conclude that hard thresholding enables better results than the soft thresholding, and combined with the adaptive weights it provides the best results among the considered algorithms.

Convergence properties of the IDD-BM3D algorithm are demonstrated in Figure 4.

The experiments with the IDD-BM3D algorithm can be reproduced using the Matlab program available as a part of the BM3D packagehttp://www.cs.tut.fi/~foi/GCF-BM3D.

VII-B Experiment 2 - comparison with the state of the art

Table III presents a comparison of the IDD-BM3D algorithm versus a number of algorithms including the current state of the art. The ISNR values for ForWaRD , SV-GSM , SA-DCT and BM3DDEB are taken from our previous paper , while the results for L0-AbS , TVMM , CGMK are obtained by the software available online. We use the default parameters suggested by the authors of the algorithms. The IDD-BM3D algorithm in this comparison employs the hard thresholding and the adaptive weights.

The proposed IDD-BM3D algorithm provides the best results with significant advantage over closest competitors. Particularly interesting is the comparison against the BM3DDEB algorithm. BM3DDEB is a two-stage non-iterative algorithm. On the first stage it utilizes the BM3D image modeling to obtain the initial estimate, which is then used on the second stage for an empirical Wiener filtering. Better performance of the IDD-BM3D algorithm demonstrates that considered decoupled formulation (35) enables more effective exploiting of the BM3D-modeling than the two-stage approach of BM3DDEB.

The visual quality of some of the restored images can be evaluated from Figures 5 and 6, where for a comparison we show results by the closest competitors , and . One can see that the proposed algorithm is able to suppress the ringing artifacts better than BM3DDEB and provides sharper image edges. This latter effect is achieved in particular due to the smaller block size used in IDD-BM3D compared to BM3DDEB.

VIII Discussion

In the experiments of the previous section we observed a clear advantage of the IDD-BM3D algorithm over the analysis-based one. This result is rather surprising, since in the case of the tight frames the IDD-BM3D and the analysis-based algorithms are almost identical.

Indeed, if we assume that $\left\{\mathbf{\phi}_{n}\right\}$ is a tight frame and require that all group weights will be equal, then $\mathbf{\Phi}^{T}\mathbf{\Phi}=\alpha\mathbf{I}$ and $\mathbf{\Psi}=\left(\mathbf{\Phi}^{T}\mathbf{\Phi}\right)^{-1}\mathbf{\Phi}^{T}=\alpha^{-1}\mathbf{\Phi}^{T}$. Substituting these expressions into equation (21) of the analysis-based algorithm we obtain

Comparing it with the equation (39) we see that up to the presence of the Lagrange multipliers the analysis-based algorithm is identical to the IDD-BM3D algorithm. This observation rises a question: what makes the algorithms behave differently when the frame is not tight?

To find an answer, let us look again at the equation (21). Its solution requires inversion of the matrix $\frac{1}{\sigma^{2}}\mathbf{A}^{T}\mathbf{A}+\frac{1}{\gamma}\mathbf{\Phi}^{T}\mathbf{\Phi}$, whoes condition number depends not only on the properties of the blur operator but also on the properties of the frame. In the case of the non-tight analysis BM3D-frame, $\mathbf{\Phi}^{T}\mathbf{\Phi}$ is a diagonal matrix, its entries are defined by the data grouping and count number of times each pixel appears in different groups. Experiments demonstrate that the variation of these entries can be very large (up to hundreds times). The large differences in magnitude of the diagonal elements of $\mathbf{\Phi}^{T}\mathbf{\Phi}$ make the matrix $\frac{1}{\sigma^{2}}\mathbf{A}^{T}\mathbf{A}+\frac{1}{\gamma}\mathbf{\Phi}^{T}\mathbf{\Phi}$ ill-conditioned and result in degradation of image reconstruction compared to IDD-BM3D.

Presence of the matrix $\mathbf{\Phi}^{T}\mathbf{\Phi}$ in the reconstruction formulas is inevitable as long as one uses criterion containing norms both for the image and spectrum domain. Formulation based on the Nash equilibrium allows to overcome this problem and have norms only from one domain in each criterion.

IX Conclusions

The frame based formulation opens new perspectives for the use of BM3D modeling within the variational reconstruction techniques. The developed deblurring algorithm demonstrates state-of-the-art performance, confirming a valuable potential of BM3D-frames as an advanced image modeling tool. For non-tight frames, we argue the validity of image reconstruction by minimizing a single objective function and propose an alternative formulation, based on Nash equilibrium problem.

Appendix A

The proof is based on use of the following Kronecker matrix product formulas.

If $\mathbf{A}$ is an $m\times n$ matrix and $\mathbf{B}$ is a $p\times q$ matrix, then the Kronecker product $\mathbf{A}\otimes\mathbf{B}$ is the $mp\times nq$ block matrix and

Also, matrix equation $\mathbf{AXB}=\mathbf{C}$ can be vectorized columnwise with respect to $\mathbf{X}$ and $\mathbf{C}$ as following

To simplify notation we denote $\mathbf{G}=\left(\mathbf{D}_{1}\otimes\mathbf{D}_{1}\right)$. Then the formula (8) from Proposition 1 is proved as following

The last identity holds since $\sum_{r}g_{r}^{2}\sum_{j\in J_{r}}\mathbf{P}_{j}^{T}\mathbf{P}_{j}$ and $\mathbf{W}^{-1}$ are diagonal matrices.

The formula (10) in Proposition 1 is valid since

Appendix B

Let us consider constrained optimization problem given in the following general form

The link between the main variable $\mathbf{u}$ and the auxiliary splitting variable $\mathbf{v}$ is given by the linear equation $\mathbf{Cv}+\mathbf{Du=b}$. If $\mathbf{C}$ is the identity matrix, then $\mathbf{v=b-Du}$ and the convergence of the corresponding iterative algorithm can be obtained from the Eckstein-Bertsekas’s theorem (, Theorem 8). However, if $\mathbf{Cv}+\mathbf{Du=b}$ is not resolved with respect to $\mathbf{v}$ then the theorem is not applicable in its original form. The techniques exploited in our paper leads to the relations between the variables which cannot be resolved with respect to $\mathbf{v}$. In order to analyze the convergence of the proposed algorithm we use a novel formulation of the Eckstein-Bertsekas’s theorem adapted to the general linear link between the variables $\mathbf{v}$ and $\mathbf{u}$. This new Eckstein-Bertsekas’s theorem is given in the following form .

If there exists a saddle point $\left(\mathbf{v}^{\ast}\mathbf{,u}^{\ast}\mathbf{,\lambda}^{\ast}\right)$ for $L\left(\mathbf{u,v,\lambda}\right)$ (43), then $v_{t}\rightarrow v^{\ast},u_{t}\rightarrow u^{\ast},\lambda_{t}\rightarrow\lambda^{\ast}$. On the other hand, if no such a saddle point exists, then at least one of the sequences $\left\{\mathbf{u}_{t}\right\}$ or $\left\{\mathbf{\lambda}_{t}\right\}$ must be unbounded.

This formulation of the convergence concerns approximate solutions on each optimization step, where the parameters $\sigma_{t}^{2}$ and $\nu_{t}$ controls the accuracy at each step. The finite sums $\sum_{t}\sigma_{t}^{2}<\infty,\sum_{t}\nu_{t}<\infty$ mean that $\sigma_{t}^{2}\mathbf{,}\nu_{t}\rightarrow 0$, i.e. the accuracy should asymptotically improve.

Armed with this theorem we can proceed to the proof of Proposition 2.

(a) Comparing the AL (17) with (42) we note that $f\left(\mathbf{u}\right)=\dfrac{1}{2\sigma^{2}}\left\|\mathbf{z}-\mathbf{Ay}\right\|_{2}^{2}$ and the equality $\mathbf{Cv}+\mathbf{Du=b}$ takes the form $\mathbf{\omega}-\mathbf{\Phi y}=0$, where $\mathbf{\omega}$ corresponds to $\mathbf{v}$ and $\mathbf{u}$ corresponds to $\mathbf{y}$. Thus, $\mathbf{C=I}_{M\times M}$ and $\mathbf{D=-\Phi}$.

We have two conditions of the theorem to be tested: $\mathbf{C}$ has full column rank and $f\left(\mathbf{u}\right)\mathbf{+}\left\|\mathbf{Du}\right\|_{2}^{2}$ is strictly convex. In our case, $\mathbf{C=I}_{M\times M}$ has full column rank, $\left\|\mathbf{Du}\right\|_{2}^{2}=\left\langle\mathbf{\Phi}^{T}\mathbf{\Phi u},\mathbf{u}\right\rangle$. Due to (8) $\mathbf{\Phi}^{T}\mathbf{\Phi=W>0}$, thus $\left\|\mathbf{Du}\right\|_{2}^{2}$ is strongly convex and the same holds for $\dfrac{1}{2\sigma^{2}}\left\|\mathbf{z}-\mathbf{Ay}\right\|_{2}^{2}+\left\|\mathbf{Du}\right\|_{2}^{2}$. Thus, all conditions of the theorem are satisfied and the analysis-based algorithm converges to the saddle-point of the AL (17), if it exists. It proves the first part of the proposition.

(b) Comparing the formulation (28) with (42) we note that $f\left(\mathbf{u}\right)=\dfrac{1}{2\sigma^{2}}\left\|\mathbf{z}-\mathbf{Ay}\right\|_{2}^{2}$ and the equality $\mathbf{Cv}+\mathbf{Du=b}$ takes the form $\mathbf{y-\Psi u=0}$ and $\mathbf{\omega-u=0}$. Assuming $\mathbf{v}\rightarrow\dbinom{\mathbf{y}}{\mathbf{u}},\mathbf{u}\rightarrow\mathbf{\omega}$ these equations give

The matrix $\mathbf{C}$ is square triangular with elements of the main diagonal equal to $1$. It has full column rank. For $\left\|\mathbf{\mathbf{Du}}\right\|_{2}^{2}$ we have $\left\|\mathbf{Du}\right\|_{2}^{2}\rightarrow\left\|\mathbf{\omega}\right\|_{2}^{2}$. Thus $\left\|\mathbf{Du}\right\|_{2}^{2}$ is strongly convex and the both conditions of the theorem are fulfilled. It proves the second part of the proposition.

B-B Proof of Proposition 3

Each iteration of the IDD-BM3D algorithm consists of two steps

where $\mathbf{M}=\frac{\gamma}{\sigma^{2}}\mathbf{A}^{T}\mathbf{A}+\mathbf{I}>0$.

Introducing the operator $O_{\text{d}}\left(\mathbf{\omega}\right)=\mathbf{\Phi M}^{-1}[\frac{\gamma}{\sigma^{2}}\mathbf{\mathbf{A}}^{T}\mathbf{z}+\mathbf{\Psi\omega}]$ and denoting $\mathbf{q}_{t}=\mathbf{\Phi y}_{t}$ we rewrite (44) in a compact form

It is shown in (Proposition 3.1) that the soft thresholding is a nonexpansive operator

Hence the operator $\mathfrak{Th}_{\tau\xi}\left(\mathbf{\cdot}\right)$ in (45) is nonexpansive.

To prove that the operator $O_{\text{d}}$ in (45) is also nonexpansive, we first notice that

To find the norm of the matrix $\mathbf{\Phi M}^{-1}\mathbf{\Psi}$ we evaluate its eigenvalues. For the matrix $\mathbf{\Phi M}^{-1}\mathbf{\Psi}$, the corresponding characteristic equation is defined as a determinant of the equation

Multiplication by $\mathbf{W}^{-1}\mathbf{\Phi}^{T}$ in (48) is legitimate because it preserves the rank of this system of the linear equations. Since $\mathbf{W}^{-1}\mathbf{\Phi}^{T}\mathbf{\Phi=I}$, (48) takes the form

Here $\lambda$ and $\mathbf{v}$ become the eigenvector and eigenvalue for the matrix $\mathbf{M}^{-1}$. The eigenvalues of the matrix $\mathbf{M}^{-1}=\left[\frac{\gamma}{\sigma^{2}}\mathbf{A}^{T}\mathbf{A}+\mathbf{I}\right]^{-1}$ are positive and take values less than or equal to $1$.

The passage from (47) to (49) proves that nonzero eigenvalues of the matrix $\mathbf{\Phi M}^{-1}\mathbf{\Psi}$ are equal to the eigenvalues of the matrix $\mathbf{M}^{-1}.$Thus all eigenvalues of the matrix $\mathbf{\Phi M}^{-1}\mathbf{\Psi}$ are nonnegative and take values less than or equal to 1. Hence, the matrix norm $\rho\left(\mathbf{\Phi M}^{-1}\mathbf{\Psi}\right)$ is less than or equal to one, and the operator $O_{\text{d}}$ is nonexpansive due to the inequality

References