#  _Cross-boosting of WNNM Image Denoising method by Directional Wavelet
Packets_

Amir Averbuch1 Pekka Neittaanmäki2 Valery Zheludev1 Moshe Salhov1 Jonathan
Hauser3
1School of Computer Science, 3School of Electrical Engineering
Tel Aviv University, Israel
2Faculty of Mathematical Information Technology
University of Jyväskylä, Finland

###### Abstract

The paper presents an image denoising scheme by combining a method that is
based on directional quasi-analytic wavelet packets (qWPs) with the state-of-
the-art Weighted Nuclear Norm Minimization (WNNM) denoising algorithm. The
qWP-based denoising method (qWPdn) consists of multiscale qWP transform of the
degraded image, application of adaptive localized soft thresholding to the
transform coefficients using the _Bivariate Shrinkage_ methodology, and
restoration of the image from the thresholded coefficients from several
decomposition levels. The combined method consists of several iterations of
qWPdn and WNNM algorithms in a way that at each iteration the output from one
algorithm boosts the input to the other. The proposed methodology couples the
qWPdn capabilities to capture edges and fine texture patterns even in the
severely corrupted images with utilizing the non-local self-similarity in real
images that is inherent in the WNNM algorithm. Multiple experiments, which
compared the proposed methodology with six advanced denoising algorithms,
including WNNM, confirmed that the combined cross-boosting algorithm
outperformed most of them in terms of both quantitative measure and visual
perception quality.

##  1 Introduction

High quality denoising is one of the main challenges in image processing. It
tries to achieve suppression of noise while capturing and preserving edges and
fine structures in the image. A huge number of publications related to a
variety of denoising methods (see, for example the reviews [21, 41, 26])
exist.

Currently, two main groups of image denoising methods exist:

  1. 1.

“Classical” schemes, which operate on single images;

  2. 2.

Methodologies based on Deep Learning.

We briefly discuss the relation between these groups of methods in Section
3.4.

Most up to date “classical” schemes, where the proposed algorithm belongs to,
are based on one of two approaches.

Utilization of non-local self-similarity (NSS) in images:

Starting from the introduction of the Non-local mean (NLM) filter in [7],
which is based on the similarity between pixels in different parts of the
image, the exploitation of various forms of the NSS in images has resulted in
a remarkable progress in image denoising. It is reflected in multiple
publications ([15, 16, 22, 44, 18, 33, 46], to name a very few). The NSS is
explored in some denoising schemes based on Deep Learning ([11, 45], for
example).

A kind of benchmark in image denoising remains the BM3D algorithm ( [15]),
which was presented as far as 2007. The algorithm exploits the self-similarity
of patches and sparsity of the image in a transform domain. It collects
similar patches in the image into a 3D array, which is subjected to a
decorrelating 3D transform followed by either hard thresholding or Wiener
filtering. After the inverse transforms, the processed patches are returned to
their original locations with corresponding weights. This method is highly
efficient in restoration of moderately noised images. However, the BM3D tends
to over-smooth and smear the image fine structure and edges when noise is
strong.

Some improvement of the original BM3D algorithm was achieved by using shape
adaptive neighborhoods and the inclusion of the Principal Component Analysis
(PCA) into the 3D transform (BM3D-SAPCA algorithm, [16]). Even better results
compared to BM3D and BM3D-SAPCA are demonstrated by the so-called Weighted
Nuclear Norm Minimization (WNNM) method ([22]), which is based on the
assumption that, by stacking the nonlocal similar patch vectors into a matrix,
this matrix should be a low rank matrix and, as such, must have sparse
singular values. The low rank matrix approximation in [22] is achieved by an
adaptive weighted thresholding SVD values of such matrices. Many denoising
algorithms presented in recent years, which are based on the NSS concept,
report results close to the results produced by BM3D-SAPCA and WNNM. At the
same time, they share, to some extent, the shortcomings of the BM3D algorithm,
especially blurring the fine structure of images restored from the severely
degraded inputs.

Transform domain filtering using directional filters:

A way to capture lines, edges and texture pattern while restoring degraded
images is to use directionional filters and, respectively, dictionaries of
waveforms oriented in multiple directions and having an oscillatory structure.
A number of dictionaries are reported in the literature and applied to image
processing. We mention contourlets [17], curvelets [9, 8], pseudo-polar
Fourier transforms [1, 2] and related to them shearlets [31, 20]. However,
while these transforms successfully capture edges in images, these
dictionaries did not demonstrate a satisfactory texture restoration due to the
shortage of oscillating waveforms in the dictionaries.

A number of publications [30, 37, 27, 6, 24, 23, 25, 28, 29], to name a few,
derive directional dictionaries by the tensor multiplication of complex
wavelets, wavelet frames and wavelet packets (WPs). The tight tensor-product
complex wavelet frames (TP_$\mathbb{C}$TFn)111The index $n$ refers to the
number of filters in the underlying one-dimensional complex tight framelet
filter bank. with different numbers of directions, are designed in [24, 25,
23] and some of them, in particular cptTP_$\mathbb{C}$TF6, TP_$\mathbb{C}$TF6
and TP_$\mathbb{C}$TF${}^{\downarrow}_{6}$, demonstrate impressive performance
for image denoising and inpainting. The waveforms in these frames are oriented
in 14 directions and, due to the 2-layer structure of their spectra, they
possess certain, although limited, oscillatory properties.

In [10] (algorithm _Digital Affine Shear Filter Transform with 2-Layer
Structure (DAS-2)_) the two-layer structure, which is inherent in the
TP_$\mathbb{C}$TF6 frames, is incorporated into shearlet-based directional
filter banks introduced in [47]. This improves the performance of DAS-2 in
comparison to TP_$\mathbb{C}$TF6 on texture-rich images, which is not the case
for smoother images.

Recently, we designed a family of complex WPs ([3], brief outlook of the
design is in [5]), which are referred to as quasi-analytic WPs (qWPs). As a
base for the design, the family of WPs originated from periodic splines of
different orders, which are described in [4] (Chapter 4), is used. The two-
dimensional (2D) qWPs are derived by a standard tensor products of 1D qWPs.
The real parts of the 2D qWPs possess a combination of properties valuable for
image processing: They are oriented in multiple directions, (see Table 2.1);
The waveforms are close to directional cosine waves with multiple frequencies
modulated by localized low-frequency 2D signals; Their DFT spectra form a
refined split of the frequency domain; Both one- and two-dimensional qWP
transforms are implemented in a very fast ways by using the Fast Fourier
transform (FFT). The directional qWPs are successfully applied to image
inpainting ([5]).

Due to the above properties, a qWP-based denoising algorithm (qWPdn), which
utilizes an adapted version of the Bivariate Shrinkage algorithm (BSA [38,
14]), proved to be efficient for image denoising. Experiments with the qWPdn
demonstrate its ability to restore edges and texture details even from
severely degraded images. In most experiments, the qWPdn to be described in
Section 3.1 provides better resolution of edges and fine structures compared
to the cptTP-$\mathbb{C}$TF6, DAS-2 and NSS-based algorithms, which is
reflected in getting higher Structural Similarity Index (SSIM)222[43], ssim.m
Matlab 2020b function. values. On the other hand, the NSS-based algorithms,
especially WNNM, proved to be superior in the noise suppression, especially in
smooth regions of images, thus producing the highest PSNR values in almost all
the experiments. However, some over-smoothing effect on the edges and fine
texture persisted under the BM3D, BM3D-SAPCA, NCSR ([18]) and WNNM algorithms.
Especially, this is the case for severely degraded images.

Therefore, we propose to combine the qWPdn and WNNM algorithms in order to
retain strong features of both algorithms and to get rid of their drawbacks.
The combined qWPdn– WNNM algorithms presented in the paper consist of the
iterated execution of the qWPdn and WNNM algorithms in a way that at each
iteration, the output from one algorithm updates (boosts) the input to the
other. Typically, 2–3 (rarely more than 4) iterations are needed to get an
excellent result.

In multiple experiments, part of which is reported in Section 3.3, the
qWPdn–WNNM algorithms performance is compared with the performance of the NSS-
based BM3D, BM3D-SAPCA, NCSR and WNNM algorithms and the algorithms
cptTP-$\mathbb{C}$TF6 and DAS-2 using directional filters. The hybrid
qWPdn–WNNM algorithms demonstrated noise suppression efficiency that is quite
competitive with all the above methods. It produces PSNR values higher than
BM3D produces and either very close to or higher than the values produced by
the BM3D-SAPCA, NCSR and WNNM algorithms. On the other hand, its performance
related to the edge and fine structures resolution is much better than the
performance of all the participated algorithms, thus, producing significantly
higher SSIM values.

This observation is illustrated by diagrams in Fig. 1.1. The left frame in
Fig. 1.1 shows PSNR values for the restoration of images degraded by Gaussian
noise with STD $\sigma=$ 5, 10, 25, 40 50, 80 and 100 dB, by the above six
methods and two combined qWPdn–WNNM algorithms designated by cbWNNM and
hybrid. The PSNR values are averaged over ten images participated in the
experiments (see Fig. 3.1). The right frame in Fig. 1.1 does the same for the
SSIM values. We can observe that the averaged PSNR values for our methods
cbWNNM and hybrid practically coincide with each other and are very close to
the values produced by BM3D-SAPCA and WNNM. A different situation we see in
the right frame, which displays averaged SSIM values. Again, the averaged
values for our methods cbWNNM and hybrid practically coincide with each other
but they strongly override the values produced by all other methods.

Figure 1.1: Diagrams of PSNR (left frame) and SSIM (right frame) values
averaged over ten images for restoration of the images by eight different
methods. Labels on x-axes indicate STD=$\sigma$ values of the additive noise
(logarithmic scale)

##### Contribution of the paper:

  * •

Development of fast image denoising algorithm qWPdn based on recently designed
directional quasi-analytic wavelet packets.

  * •

Design of the iterative cross-boosting qWPdn–WNNM algorithms, which are highly
efficient in noise suppression and capturing edges and fine structures even in
severely degraded images.

  * •

Experimental comparison of the qWPdn–WNNM algorithms performance with the
performance of multiple state-of-the-art algorithms, which demonstrated a
decisive advantage of thesealgorithms in a SSIM sense and visual perception.

The paper is organised as follows: Section 2 briefly outlines properties of
the qWP transforms in one and two dimensions. Section 3.1 describes the qWPdn
algorithm. Section 3.2 presents the combined qWPdn–WNNM algorithms. In the
multiple experiments in Section 3.3, the performance of these algorithms is
compared with the performance of the cptTP-$\mathbb{C}$TF6, DAS-2, BM3D,
BM3D-SAPCA, NCSR and WNNM algorithms. Section 3.4 briefly discusses the
relation of the proposed methodology to the recently published Deep Learning
denoising methods. Section 4 provides an overview of the results.

##### Notation and abbreviations:

$N=2^{j}$, $\omega\stackrel{{\scriptstyle\mathrm{def}}}{{=}}e^{2\pi\,i/N}$ and
$\Pi[N]$ is a space of real-valued $N$-periodic signals. $\Pi[N,N]$ is the
space of two-dimensional $N$-periodic in both vertical and horizontal
directions arrays. DFT(FFT) means Discrete(Fast) Fourier transform.

The abbreviations WP, dWP and qWP mean wavelet packet, orthonormal spline-
based wavelet packet $\psi^{p}_{[m],l}$ and quasi-analytic wavelet packets
$\Psi^{p}_{\pm[m],l}$, respectively, in a 1D case, and orthonormal WPs
$\psi^{p}_{[m],j,l}$ and quasi-analytic wavelet packets
$\Psi^{p}_{+\pm[m],l,j}$, respectively, in a 2D case.

qWPdn designates the qWP-based image denoising algorithm. qWPdn–WNNM means a
cross-boosting image denoising algorithm combining the qWPdn with the BM3D.

PSNR means Peak Signal-to-Noise ratio in decibels (dB). SSIM means Structural
Similarity Index ([43]) computed by the Matlab 2020b function ssim.m. BSA
stands for Bivariate Shrinkage algorithm ([14, 38]). NSS means non-local self-
similarity.

BM3D stands for _Block-matching and 3D filtering_ ([15]), SAPCA means Shape-
Adaptive Principal Component Analysis ([16]), NCSR means Nonlocally
Centralized Sparse Representation ([18]), WNNM means Weighted Nuclear Norm
Minimization ([22]), cptTP-$\mathbb{C}$TF stands for _Compactly Supported
Tensor Product Complex Tight Framelets with Directionality_ ([48]) and DAS-2
stands for _Digital Affine Shear Filter Transform with 2-Layer Structure_
([10]).

##  2 Preliminaries: Quasi-analytic directional wavelet packets

Recently we designed a family of quasi-analytic wavelet packets (qWPs), which
possess a collection of properties indispensable for image processing. A brief
outline of the qWPs design and the implementation of corresponding transforms
is provided the paper [5], which describes successful application of qWPs to
image inpainting. A detailed description of the design and implementation is
given in [3]. In this section we list properties of qWPs and present some
illustrations.

###  2.1 Properties of qWPs

##### One-dimentional qWPs

The qWPs are derived from the periodic WPs originating from orthonormal
discretized polynomial splines of different orders (dWPs), which are described
in Chapter 4 in [4] (a brief outline is given in [3]). The dWPs are denoted by
$\psi^{p}_{[m],l}$, where $p$ means the generating spline’s order, $m$ is the
decomposition level and $l=0,...2^{m}-1,$ is the index of an $m$-level wavelet
packets. The $2^{m}$-sample shifts
$\left\\{\psi^{p}_{[m],l}(\cdot-2^{m}\,k)\right\\},\;l=0,...,2^{m}-1,\;k=0,...,N/2^{m}-1,$
of the $m$-level dWPs form an orthonormal basis of the space $\Pi[N]$ of
$N$-periodic discrete-time signals. Surely, other orthonormal bases are
possible, for example, wavelet and Best bases ([12]).

The waveforms $\psi^{p}_{[m],l}[k]$ are symmetric, well localized in the
spatial domain and have oscillatory structure. Their DFT spectra form a
refined split of the frequency domain. The shapes of magnitude spectra tend to
rectangular as the spline’s order $p$ grows. A common way to extend 1D WP
transforms to multiple dimensions is by the tensor-product extension. The 2D
dWPs from the level $m$ are:
$\psi_{[m],j,l}^{p}[k,n]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\psi_{[m],j}^{p}[k]\,\psi_{[m],l}^{p}[n]$.
Their $2^{m}$-sample shifts along vertical and horizontal directions form
orthonormal bases of the space $\Pi[N,N]$ of 2D signals $N-$periodic in both
directions. The drawback for image processing is the lack of directionality.
The directionality can be achieved by switching to complex wavelet packets.

For this, we start with application of the Hilbert transform (HT) to the dWPs
$\psi^{p}_{[m],l}$, thus getting the signals
$\tau^{p}_{[m],l}=H(\psi^{p}_{[m],l}),\;m=1,...,M,\;l=0,...,2^{m}-1$. A slight
correction of those signals spectra:

$\hat{\phi}^{p}_{[m],l}[n]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\hat{\psi}^{p}_{[m],l}[0]+\hat{\psi}^{p}_{[m],l}[N/2]+\hat{\tau}^{p}_{[m],l}[n]$
(2.1)

provides us with a set of signals from the space $\Pi[N]$, whose properties
are similar to the properties of the dWPs $\psi^{p}_{[m],l}$. In particular,
their shifts form orthonormal bases in $\Pi[N]$, their magnitude spectra
coincide with the magnitude spectra of the dWPs $\psi^{p}_{[m],l}$. However,
unlike the symmetric dWPs $\psi^{p}_{[m],l}$, the signals $\phi^{p}_{[m],l}$
are antisymmetric for all $l$ except for $l_{0}=0$ and $l_{m}=2^{m}-1$. We
refer to the signals $\phi^{p}_{[m],l}$ as the complementary orthonormal WPs
(cWPs).

The sets of complex-valued WPs, which we refer to as the quasi-analytic
wavelet packets (qWP), are defined as
$\Psi^{p}_{\pm[m],l}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\psi^{p}_{[m],l}\pm
i\phi^{p}_{[m],l},\quad m=1,...,M,\;l=0,...,2^{m}-1$, where $\phi^{p}_{[m],l}$
are the cWPs defined in Eq. (2.1). The qWPs $\Psi^{p}_{\pm[m],l}$ differ from
the analytic WPs by adding two values $\pm i\,\hat{\psi}^{p}_{[m],l}[0]$ and
$\pm i\,\hat{\psi}^{p}_{[m],l}[N/2]$ into their DFT spectra, respectively. The
DFT spectra of the qWPs $\Psi^{p}_{+[m],l}$ are located within positive half-
band of the frequency domain and vice versa for the qWPs $\Psi^{p}_{-[m],l}$.

Figure 2.1 displays the signals ${\psi}^{9}_{[3],l}$ and
${\phi}^{9}_{[3],l},\;l=0,...,7$, from the third decomposition level and their
magnitude spectra (right half-band), that coincide with each other. Addition
of $\hat{\psi}^{9}_{[3],l}[0]$ and $\hat{\psi}^{9}_{[3],l}[N/2]$ to the
spectra of ${\phi}^{9}_{[3],l},\;l=0,7,$ results in an antisymmetry
distortion. These WPs provide a collection of diverse symmetric and
antisymmetric well localized waveforms, which range from smooth wavelets for
$l=0,1$ to fast oscillating transients for $l=5,6,7$. Thus, this collection is
well suited to catching smooth as well as oscillating local patterns in
signals. In the 2D case, these valuable properties of the spline-based wavelet
packets are completed by the directionality of the tensor-product waveforms.

Figure 2.1: Top to bottom: signals ${\psi}^{9}_{[3],l}$; signals
${\phi}^{9}_{[3],l},\;l=0,...,7$: their magnitude DFT spectrum (right half-
band); magnitude DFT spectra of complex qWPs ${\Psi}^{9}_{+[3],l}$; same for
${\Psi}^{9}_{-[3],l},\;l=0,...,7$

##### Two-dimensional qWPs

Similarly to the 2D dWPs $\psi_{[m],j,l}^{p}[k,n]$, the 2D cWPs
$\phi_{[m],j,l}^{p}[k,n]$ are defined as the tensor products of 1D WPs such
that
$\phi_{[m],j,l}^{p}[k,n]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\phi_{[m],j}^{p}[k]\,\phi_{[m],l}^{p}[n].$
The $2^{m}$-sample shifts of the WPs
$\left\\{\phi_{[m],j,l}^{p}\right\\},\;j,l=0,...,2^{m}-1,$ in both directions
form an orthonormal basis for the space $\Pi[N,N]$ of arrays that are
$N$-periodic in both directions.

##### 2D qWPs and their spectra

The 2D dWPs $\left\\{\psi_{[m],j,l}^{p}\right\\}$ as well as the cWPs
$\left\\{\phi_{[m],j,l}^{p}\right\\}$ lack the directionality property which
is needed in many applications that process 2D data. However, real-valued 2D
wavelet packets oriented in multiple directions can be derived from tensor
products of complex quasi-analytic qWPs $\Psi_{\pm[m],\rho}^{p}$. The complex
2D qWPs are defined as follows:

$\displaystyle\Psi_{++[m],j,l}^{p}[k,n]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\Psi_{+[m],j}^{p}[k]\,\Psi_{+[m],l}^{p}[n],\quad\Psi_{+-[m],j,l}^{p}[k,n]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\Psi_{+[m],j}^{p}[k]\,\Psi_{-[m],l}^{p}[n],$
(2.2)

where $m=1,...,M,\;j,l=0,...,2^{m}-1,$ and $k,n=0,...,N-1$. The real parts of
these 2D qWPs are

$\begin{array}[]{lll}\theta_{+[m],j,l}^{p}[k,n]&\stackrel{{\scriptstyle\mathrm{def}}}{{=}}&\mathfrak{Re}(\Psi_{++[m],j,l}^{p}[k,n])=\psi_{[m],j,l}^{p}[k,n]-\phi_{[m],j,l}^{p}[k,n],\\\
\theta_{-[m],j,l}^{p}[k,n]&\stackrel{{\scriptstyle\mathrm{def}}}{{=}}&\mathfrak{Re}(\Psi_{+-[m],j,l}^{p}[k,n])=\psi_{[m],j,l}^{p}[k,n]+\phi_{[m],j,l}^{p}[k,n].\\\
\end{array}$ (2.3)

The block-scheme in Fig. 2.2 illustrates the design of qWPs.

Figure 2.2: Block-scheme of the qWP design (left) and quadrants of frequency
domain (right)

The DFT spectra of the 2D qWPs $\Psi_{++[m],j,l}^{p},\;j,l=0,...,2^{m}-1,$ are
tensor products of the one-sided spectra of the qWPs
$\hat{\Psi}_{++[m],j,l}^{p}[p,q]=\hat{\Psi}_{+[m],j}^{p}[p]\,\hat{\Psi}_{+[m],l}^{p}[q]$
and, as such, they fill the quadrant $\mathbf{q}_{0}$ of the frequency domain,
while the spectra of $\Psi_{+-[m],j,l}^{p},\;j,l=0,...,2^{m}-1,$ fill the
quadrant $\mathbf{q}_{1}$ (see Fig. 2.2). Figure 2.3 displays magnitude
spectra of the ninth-order 2D qWPs $\Psi_{++[2],j,l}^{9}$ and
$\Psi_{+-[2],j,l}^{9}$ from the second decomposition level.

Figure 2.3: Magnitude spectra of 2D qWPs $\Psi_{++[2],j,l}^{9}$ (left) and
$\Psi_{+-[2],j,l}^{9}$ (right) from the second decomposition level

Figure 2.3 shows that the DFT spectra of the qWPs $\Psi_{+\pm[m],j,l}^{9}$
effectively occupy relatively small squares in the frequency domain. For
deeper decomposition levels, sizes of the corresponding squares decrease as
geometric progression. Such configuration of the spectra leads to the
directionality of the real-valued 2D WPs $\theta_{\pm[m],j,l}^{p}$. The
directionality of the WPs $\theta_{\pm[m],j,l}^{p}$ is discussed in [3]. It is
established that if the spectrum of a WP $\Psi_{+\pm[m],j,l}^{p}$ occupies a
square whose center lies in the point $[\kappa_{0},\nu_{0}]$, then the
respective real-valued WP $\theta_{\pm[m],j,l}^{p}$ defined in Eq. (2.3) is
represented by
$\theta_{\pm[m],j,l}^{p}[k,n]\approx{\cos\frac{2\pi(\kappa_{0}k+\nu_{0}n)}{N}}\,\underline{\theta}[k,n],$
where $\underline{\theta}[k,n]$ is a spatially localized low-frequency
waveform which does not have a directionality. But the 2D signal
$\cos\frac{2\pi(\kappa_{0}k+\nu_{0}n)}{N}$ is oscillating in the direction D,
which is orthogonal to the vector $\vec{V}=\kappa_{0}\vec{i}+\nu_{0}\vec{j}$.
Therefore, WP $\theta_{\pm[m],j,l}^{p}$ can be regarded as the directional
cosine wave modulated by the localized low-frequency signal
$\underline{\theta}$. The cosine frequencies in the vertical and horizontal
directions are determined by the indices $j$ and $l$, respectively, of the WP
$\theta_{\pm[m],j,l}^{p}$. The bigger is the index, the higher is frequency in
the respective direction. The situation is illustrated in Fig. 2.4. The
imaginary parts of the qWPs $\Psi_{+\pm[m],j,l}^{p}$ have a similar structure.

Figure 2.4: Left: magnitude spectrum of 2D qWP $\Psi_{++[3],2,5}^{p}[k,n]$.
Right WP $\theta_{++[3],2,5}^{p}=\mathfrak{Re}(\Psi_{++[3],2,5}^{p})$

Figure 2.5 displays WPs $\theta_{+[2],j,l}^{9},\;j,l=0,1,2,3,$ from the second
decomposition level and their magnitude spectra.

Figure 2.5: WPs $\theta_{+[2],j,l}^{9}$ from the second decomposition level
and their magnitude spectra

Figure 2.6 displays WPs $\theta_{-[2],j,l}^{9},\;j,l=0,1,2,3,$ from the second
decomposition level and their magnitude spectra.

Figure 2.6: WPs $\theta_{-[2],j,l}^{9}$ from the second decomposition level
and their magnitude spectra

Figure 2.7 displays WPs $\theta_{+[3],j,l}^{9}$ and
$\theta_{-[3],j,l}^{9},\;j,l=0,1,2,3,$ from the third decomposition level.

Figure 2.7: WPs $\theta_{+[3],j,l}^{9}$ (left) and $\theta_{-[3],j,l}^{9}$
(right) from the third decomposition level

###### Remark 2.1

Note that all the WPs $\theta_{+[m],j,l}^{p}$ whose spectra are located along
the vector $\vec{V}$ have approximately the same orientation. It is seen in
Figs. 2.6, 2.5 and 2.7. Consequently, the number of orientations of the $m$-th
level WPs is less than the number of WPs, which is $2\cdot 4^{m}$. For
example, all the “diagonal” qWPs
$\left\\{\theta_{\pm[m],j,j}^{p}\right\\},\;j=0,...,2^{m}-1,$ are oscillating
with different frequencies in the directions of either $135^{o}$ (for
$\theta_{+}$) or $45^{o}$ (for $\theta_{-}$). Orientation numbers are given in
Table 2.1.

Table 2.1:  Numbers of different orientations of qWPs $\left\\{\theta_{\pm[m],j,l}^{p}\right\\},\;j,l=0,...,2^{m}-1$, for different decomposition levels level $m$ | 1 | 2 | 3 | 4 | 5 | 6 | …
---|---|---|---|---|---|---|---
# of directions | 6 | 22 | 86 | 318 | 1290 | 5030 | …

###  2.2 Outline of the implementation scheme for 2D qWP transforms

The spectra of 2D qWPs
$\left\\{\Psi_{++[m],j,l}^{p}\right\\},\;j,l=0,...,2^{m}-1$ fill the quadrant
$\mathbf{q}_{0}$ of the frequency domain (see Fig. 2.2), while the spectra of
2D qWPs $\left\\{\Psi_{+-[m],j,l}^{p}\right\\}$ fill the quadrant
$\mathbf{q}_{1}$. Consequently, the spectra of the real-valued 2D WPs
$\left\\{\theta_{+[m],j,l}^{p}\right\\},\;j,l=0,...,2^{m}-1$, and
$\left\\{\theta_{-[m],j,l}^{p}\right\\}$ fill the pairs of quadrant
$\mathbf{q}_{+}=\mathbf{q}_{0}\bigcup\mathbf{q}_{2}$ and
$\mathbf{q}_{-}=\mathbf{q}_{1}\bigcup\mathbf{q}_{3}$, respectively.

By this reason, none linear combination of the WPs
$\left\\{\theta_{+[m],j,l}^{p}\right\\}$ and their shifts can serve as a basis
for the signal space $\Pi[N,N]$. The same is true for WPs
$\left\\{\theta_{-[m],j,l}^{p}\right\\}$. However, combinations of the WPs
$\left\\{\theta_{+[m],j,l}^{p}\right\\}$ and
$\left\\{\theta_{-[m],j,l}^{p}\right\\}$ provide frames of the space
$\Pi[N,N]$.

The transforms are implemented in the frequency domain using modulation
matrices of the filter banks, which are built from the corresponding wavelet
packets. It is important to mention that the structure of the filter banks
$\mathbf{Q}_{+}$ and $\mathbf{Q}_{-}$ for the first decomposition level is
different for the transforms with the “positive” $\Psi_{+[m],l}^{p}$ and
“negative” $\Psi_{-[m],l}^{p}$ qWPs, respectively. However, the transforms
from the first to the second and further decomposition levels are executed
using the same filter bank $\mathbf{H}_{m}$ for the “positive” and “negative”
qWPs. This fact makes it possible a parallel implementation of the transforms.

The one-level 2D qWP transforms of a signal
$\mathbf{X}=\left\\{X[k,n]\right\\}\in\Pi[N,N]$ are implemented by a tensor-
product scheme. To be specific, for the transform with $\Psi^{p}_{++[1]}$, the
1D transform of rows from the signal $\mathbf{X}$ is executed using the filter
bank $\mathbf{Q}_{+}$, which is followed by the 1D transform of columns of the
produced coefficient arrays using the same filter bank $\mathbf{Q}_{+}$. These
operations result in the transform coefficient array
$\mathbf{Z}_{+[1]}=\bigcup_{\j,l=0}^{1}\mathbf{Z}_{+[1]}^{j,l}$ comprising of
four blocks of size $N/2\times N/2$. The transform with $\Psi^{p}_{+-[1]}$ is
implemented by the subsequent application of the filter banks $\mathbf{Q}_{+}$
and $\mathbf{Q}_{-}$ to rows from the signal $\mathbf{X}$ and columns of the
produced coefficient arrays, respectively. This results in the coefficient
array $\mathbf{Z}_{-[1]}=\bigcup_{\j,l=0}^{1}\mathbf{Z}_{-[1]}^{j,l}$,

The further transforms starting from the arrays $\mathbf{Z}_{+[1]}$ and
$\mathbf{Z}_{-[1]}$ produce two sets of the coefficients
$\left\\{\mathbf{Z}_{+[m]}=\bigcup_{\j,l=0}^{2^{m}-1}\mathbf{Z}_{+[m]}^{j,l}\right\\}$
and
$\left\\{\mathbf{Z}_{-[m]}=\bigcup_{\j,l=0}^{2^{m}-1}\mathbf{Z}_{-[m]}^{j,l}\right\\},\;m=2,...,M$.
The transforms are implemented by the application of the same filter banks
$\mathbf{H}_{m},\;m=2,...,M$ to rows and columns of the “positive” and
“negative” coefficient arrays. The coefficients from a level $m$ comprise of
$4^{m}$ “positive” blocks of coefficients
$\left\\{\mathbf{Z}_{+[m]}^{j,l}\right\\},\;l,j=0,...,^{2^{m}-1},$ and the
same number of “negative” blocks $\left\\{\mathbf{Z}_{-[m]}^{j,l}\right\\}$.

The coefficients from a block are inner products of the signal
$\mathbf{X}=\left\\{X[k,n]\right\\}\in\Pi[N,N]$ with the shifts of the
corresponding wavelet packet:

$\begin{array}[]{cc}Z_{\pm[m]}^{j,l}[k,n]=&\sum_{\lambda,\mu=0}^{N-1}X[\lambda,\mu]\,{\Psi}^{p}_{+\pm[m],j,l}[\lambda-2^{m}k,\mu-2^{m}n],\\\
Y_{\pm[m]}^{j,l}[k,n]=&\mathfrak{Re}(Z_{\pm[m]}^{j,l}[k,n])=\sum_{\lambda,\mu=0}^{N-1}X[\lambda,\mu]\,{\theta}^{p}_{\pm[m],j,l}[\lambda-2^{m}k,\mu-2^{m}n].\end{array}$
(2.4)

The inverse transforms are implemented accordingly. Prior to the
reconstruction, some structures, possibly different, are defined in the sets
$\left\\{\mathbf{Z}_{+[m]}^{j,l}\right\\}$ and
$\left\\{\mathbf{Z}_{-[m]}^{j,l}\right\\},\;m=1,...M,$ (for example, 2D
wavelet, Best Basis or single-level structures) and some manipulations on the
coefficients, (for example, thresholding, shrinkage, $l_{1}$ minimization) are
executed. The reconstruction produces two complex arrays $\mathbf{X}_{+}$ and
$\mathbf{X}_{-}$. The signal X is restored by
$\mathbf{\tilde{X}}=\mathfrak{Re}(\mathbf{X}_{+}+\mathbf{X}_{-})/8$.

Figure 2.8 illustrates the image “Fingerprint” restoration by the 2D signals
$\mathfrak{Re}(\mathbf{X}_{\pm})$. The signal $\mathfrak{Re}(\mathbf{X}_{-})$
captures oscillations oriented to _north-east_ , while
$\mathfrak{Re}(\mathbf{X}_{+})$ captures oscillations oriented to _north-
west_. The signal
$\tilde{\mathbf{X}}=\mathfrak{Re}(\mathbf{X}_{+}+\mathbf{X}_{-})/8$ perfectly
restores the image achieving PSNR=312.3538 dB.

Figure 2.8: Left to right: 1.Image $\mathfrak{Re}(\mathbf{X}_{+})$. 2. Its
magnitude DFT spectrum. 3.Image $\mathfrak{Re}(\mathbf{X}_{-})$. 4. Its
magnitude DFT spectrum

##  3 Image denoising

In this section, we describe application of the directional qWP transforms
presented in Section 2 to the restoration of an image $\mathbf{X}$ from the
data $\mathbf{\check{X}}=\mathbf{X}+\mathbf{E}$, where $\mathbf{E}$ is the
Gaussian zero-mean noise whose STD=$\sigma$.

###  3.1 Denoising scheme for 2D qWPs

The degraded image $\mathbf{\check{X}}$ is decomposed into two sets
$\left\\{\mathbf{\check{Z}}_{+[m]}^{j,l}\right\\}$ and
$\left\\{\mathbf{\check{Z}}_{-[m]}^{j,l}\right\\},\;m=1,...,M,\;j,l=0,2^{m}-1,$
of the qWP transform coefficients, then a version of the Bivariate Shrinkage
algorithm (BSA)[38, 14] is implemented and the image
$\tilde{\mathbf{X}}\approx\mathbf{X}$ is restored from the shrunken
coefficients. The restoration is executed separately from the sets of
coefficients belonging to several decomposition levels and the results are
averaged with some weights.

####  3.1.1 Image restoration from a single-level transform coefficients

Consider the reconstruction of an image $\mathbf{X}\in\Pi[N,N]$ from the
fourth-level transform coefficients of the degraded array $\mathbf{\check{X}}$
of size $N\times N$.

The denoising algorithm, which we refer to as qWPdn, is implemented by the
following steps:

  1. 1.

In order to eliminate boundary effects, the degraded image
$\mathbf{\check{X}}$ is symmetrically extended to the image
$\mathbf{\check{X}}_{T}$ of size $N_{T}\times N_{T}$, where $N_{T}=N+2T$.
Typically, either $T=N/4\mbox{ or }T=N/8$.

  2. 2.

The Bivariate Shrinkage (BSA) utilizes the interscale dependency of the
transform coefficients. Therefore, the direct 2D transforms of the image
$\mathbf{\check{X}}_{T}$ with the complex qWPs $\Psi^{p}_{++}$ and
$\Psi^{p}_{+-}$ are executed down to the fifth decomposition level. As a
result, two sets
$\mathbf{\check{Z}}_{+[m]}^{j,l}=\left\\{\mathbf{\check{Z}}_{+[m]}^{j,l}[k,n]\right\\}$
and
$\mathbf{\check{Z}}_{-[m]}^{j,l}=\left\\{\mathbf{\check{Z}}_{+[m]}^{j,l}[k,n]\right\\},\;m=1,...,5,\;j,l=0,...,2^{m}-1,\;k,n=0,...,N_{T}/2^{m}-1,$
of the qWP transform coefficients are produced.

  3. 3.

The noise variance is estimated by
$\tilde{\sigma}_{e}^{2}=\frac{\mathrm{median}(|{\check{Z}}_{+[1]}^{1,1}[k,n]|)}{0.6745}.$

  4. 4.

Let
$\check{c}_{4}[k,n]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\check{Z}_{+[4]}^{j,l}[k,n]$
denote a coefficient from a block $\mathbf{\check{Z}}_{+[4]}^{j,l}$ at the
fourth decomposition level. The following operations are applied to the
coefficient $\check{c}_{4}[k,n]$:

    1. (a)

The averaged variance
$\bar{\sigma}_{c}[k,n]^{2}=\frac{1}{W_{4}^{2}}\sum_{\kappa,\nu=-W_{4}/2}^{W_{4}/2-1}\check{c}_{4}[k+\kappa,n+\nu]^{2}$
is calculated. The integer $W_{4}$ determines the neighborhood of
$\check{c}_{4}[k,n]$ size.

    2. (b)

The marginal variance for $\check{c}_{4}[k,n]$ is estimated by
$\tilde{\sigma}[k,n]^{2}=(\bar{\sigma}_{c}[k,n]^{2}-\tilde{\sigma}_{e}^{2})_{+}.$333$s_{+}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\max\\{s,0\\}$.

    3. (c)

In order to estimate the clean transform coefficients from the fourth
decomposition level, coefficients from the fifth level should be utilized. The
size of the coefficient block $\mathbf{\check{Z}}_{+[4]}^{j,l}$ is
$N_{T}/16\times N_{T}/16$. The coefficients from that block are related to the
qWP $\Psi^{p}_{++[4],j,l}$, whose spectrum occupies, approximately, the square
$\mathbf{S}_{+[4]}^{j,l}$ of size $N_{T}/32\times N_{T}/32$ within the
quadrant $\mathbf{q}_{0}$ (see Fig. 2.2). The spectrum’s location determines
the directionality of the waveform $\Psi^{p}_{++[4],j,l}$. On the other hand,
four coefficient blocks
$\left\\{\mathbf{\check{Z}}_{+[5]}^{2j+\iota,2l+\lambda}\right\\},\;\iota,\lambda=0,1,$
of size $N_{T}/32\times N_{T}/32$ are derived by filtering the block
$\mathbf{\check{Z}}_{+[4]}^{j,l}$ coefficients followed by downsampling. The
coefficients from those blocks are related to the qWPs
$\Psi^{p}_{++[5],2j+\iota,2l+\lambda}$, whose spectra occupy, approximately,
the squares $\mathbf{S}_{+[5]}^{2j+\iota,2l+\lambda}$ of size $N_{T}/64\times
N_{T}/64$, which fill the square $\mathbf{S}_{+[4]}^{j,l}$. Therefore, the
orientations of the waveforms $\Psi^{p}_{++[5],2j+\iota,2l+\lambda}$ are close
to the orientation of $\Psi^{p}_{++[4],j,l}$. Keeping this in mind, we form
the joint fifth-level array $\mathbf{c}_{5}^{j,l}$ of size $N_{T}/16\times
N_{T}/16$ by interleaving the coefficients from the arrays
$\left\\{\mathbf{\check{Z}}_{+[5]}^{2j+\iota,2l+\lambda}\right\\}$. To be
specific, the joint array $\mathbf{c}_{5}^{j,l}$ consists of the quadruples:

$\mathbf{c}_{5}^{j,l}=\left\\{\left[\begin{array}[]{cc}\check{Z}_{+[5]}^{2j,2l}[\kappa,\nu]&\check{Z}_{+[5]}^{2j,2l+1}[\kappa,\nu]\\\
\check{Z}_{+[5]}^{2j+1,2l}[\kappa,\nu]&\check{Z}_{+[5]}^{2j+1,2l+1}[\kappa,\nu]\\\
\end{array}\right]\right\\},\quad\kappa,\nu=0,...N_{T}/32-1.$

    4. (d)

Let $\breve{c}_{5}[k,n]$ denote a coefficient from the joint array
$\mathbf{c}_{5}^{j,l}$. Then, the transform coefficient $Z^{j,l}_{+[4]}[k,n]$
from the fourth decomposition level is estimated by the bivariate shrinkage of
the coefficients $\check{Z}^{j,l}_{+[4]}[k,n]$:

$Z^{j,l}_{+[4]}[k,n]\approx\tilde{Z}^{j,l}_{+[4]}[k,n]=\frac{\left(\sqrt{\check{c}_{4}[k,n]^{2}+\breve{c}_{5}[k,n]^{2}}-\frac{\sqrt{3}\,\tilde{\sigma}_{e}^{2}}{\tilde{\sigma}[k,n]}\right)_{+}}{\sqrt{\check{c}_{4}[k,n]^{2}+\breve{c}_{5}[k,n]^{2}}}\,\check{c}_{4}[k,n].$

  5. 5.

As a result of the above operations, the fourth-level coefficient array
$\mathbf{\tilde{Z}}_{+[4]}=\left\\{\mathbf{\tilde{Z}}_{+[4],j,l}\right\\},\;j,l=0,...,15$,
is estimated, where
$\mathbf{\tilde{Z}}_{+[4],j,l}=\left\\{\tilde{Z}^{j,l}_{+[4]}[k,n]\right\\},\;k,n=0,...N_{T}/16-1$.

  6. 6.

The inverse qWP transform is applied to the coefficient array
$\mathbf{\tilde{Z}}_{+[4]}$ and the result shrinks to the original image size
$N\times N$. Thus, the sub-image $\mathbf{\tilde{X}}_{+}^{4}$ is obtained.

  7. 7.

The same operations are applied to $\check{\mathbf{Z}}^{j,l}_{-[m]},\;m=4,5$,
thus resulting in the sub-image $\mathbf{\tilde{X}}_{-}^{4}$.

  8. 8.

The clean image is estimated by
$\mathbf{X}\approx\mathbf{\tilde{X}}^{4}=\frac{\mathfrak{Re}(\mathbf{\tilde{X}}_{+}^{4}+\mathbf{\tilde{X}}_{-}^{4})}{8}.$

####  3.1.2 Image restoration from several decomposition levels

More stable estimation of the image $\mathbf{X}$ is derived by the weighted
average of several single-level estimations
$\left\\{\mathbf{\tilde{X}}^{m}\right\\}$. In most cases, the estimations from
the second, third and fourth levels are combined, so that $m=2,3,4$.

The approximated image $\mathbf{\tilde{X}}^{3}$ is derived from the third-
level coefficients $\mathbf{\check{Z}}_{\pm[3]}^{j,l}$. The fourth-level
coefficients that are needed for the Bivariate Shrinkage of the coefficients
$\mathbf{\check{Z}}_{\pm[3]}^{j,l}$ are taken from the “cleaned” arrays
$\mathbf{\tilde{Z}}_{\pm[4],j,l}$ rather than from the “raw” ones
$\mathbf{\check{Z}}_{\pm[4],j,l}$.

Similarly, the image $\mathbf{\tilde{X}}^{2},$ is derived from the coefficient
arrays $\mathbf{\check{Z}}_{\pm[2]}^{j,l}$ and
$\mathbf{\tilde{Z}}_{\pm[3]}^{j,l}$. The final operation is the weighted
averaging such as

$\mathbf{\tilde{X}}=\frac{\alpha_{2}\mathbf{\tilde{X}}^{2}+\alpha_{3}\mathbf{\tilde{X}}^{3}+\alpha_{4}\mathbf{\tilde{X}}^{4}}{\alpha_{2}+\alpha_{3}+\alpha_{4}}.$
(3.1)

###### Remark 3.1

Matlab implementation of all the operations needed to transform the degraded
array $\mathbf{\check{X}}$ of size $512\times 512$ into the estimation
$\mathbf{\tilde{X}}$ given by Eq. (3.1) takes 1 second. Note that the noise
STD is not a part of the input. It is evaluated as indicated in Item 3.

###### Remark 3.2

In some cases, restoration from third, fourth and fifth levels is preferable.
Then, the degraded array $\mathbf{\check{X}}_{T}$ is decomposed down to the
sixth level.

###### Remark 3.3

The algorithm comprises a number of free parameters which enable a flexible
adaptation to the processed class of objects. These parameters are the order
“p” of the generating spline, integers $W_{4}$, $W_{3}$ and $W_{2}$, which
determine the sizes of neighborhoods for the averaged variances calculation,
and the weights $\alpha_{2}$, $\alpha_{3}$ and $\alpha_{4}$.

###### Remark 3.4

Fragments of the Matlab functions denoising_dwt.m and bishrink.m from the
websites http://eeweb.poly.edu/iselesni/WaveletSoftware/denoising_dwt.html and
http://eeweb.poly.edu/iselesni/WaveletSoftware/denoise2.html, respectively,
were used as patterns while compiling our own denoising software.

###  3.2 qWPdn–WNNM: Hybrid algorithm

Equations (2.4) imply that the coefficients of the qWP transforms reflect the
correlation of the image under processing with the collection of waveforms,
which are well localized in the spatial domain and are oscillating in multiple
directions with multiple frequencies. By this reason, these transforms are
well suited for capturing edges and texture patterns oriented in different
directions. Experiments with the qWPdn image denoising algorithm demonstrate
its ability to restore edges and texture details even from severely degraded
images. In most conducted experiments, the qWPdn provides better resolution of
edges and fine structures compared to the most state-of-the-art algorithms
based on non-local self-similarity (NSS), including the BM3D, NCSR and WNNM
algorithms, which is reflected in higher SSIM values. On the other hand, the
NSS-based algorithms proved to be superior in noise suppression especially in
smooth regions in images, thus producing the high PSNR values.

One of the best existing denoising algorithms is WNNM ([22]) introduced in
2014. It produces high PSNR and SSIM values and, in most cases, a better
visual perception of restored images compared to other algorithms. WNNM is an
iterative image denoising algorithm, whose main components are stacking non-
local patches that are similar to a given patch into a low rank matrix,
computing the singular value decomposition of this matrix and minimization of
the the nuclear norm of the matrix by soft thresholding with adaptive weights
the singular values of this matrix. Even the recently designed methods (for
example, [46], which exploits the similarity of patches from external images
in addition to the inner patches from the given image, and Deep Learning-based
methods such as reported in [11, 45]) produce marginal, if any, improvements
compared to WNNM. Nevertheless, some drawbacks common to NSS-based methods are
inherent in WNNM. Namely, some over-smoothing effect on the edges and fine
texture persists when restoration of severely degraded images.

We propose to combine qWPdn with WNNM algorithms to benefit from the strong
features of both algorithms.

Denote by $\mathbf{Q}$ and $\mathbf{W}$ the operators of application of the
qWPdn , which is described in Section 2, and WNNM denoising algorithms,
respectively, to a degraded array $\mathbf{A}$:
$\mathbf{Q}\,\mathbf{A}=\mathbf{D}_{Q}$ and
$\mathbf{W}\,\mathbf{A}=\mathbf{D}_{W}$.

Assume that we have an array $\check{\mathbf{X}}^{0}=\mathbf{X}+\mathbf{E}$,
which represents an image $\mathbf{X}$ degraded by additive Gaussian noise
$\mathbf{E}$ whose STD is $\sigma.$ The denoising processing is implemented
along the following cross-boosting scheme.

First step:

Apply the operators $\mathbf{Q}$ and $\mathbf{W}$ to the input array
$\check{\mathbf{X}}^{0}$:
$\mathbf{Y}_{Q}^{1}=\mathbf{Q}\,\check{\mathbf{X}}^{0}$ and
$\mathbf{Y}_{W}^{1}=\mathbf{W}\,\check{\mathbf{X}}^{0}$.

Iterations:

$i=1,...,I-1$

  1. 1.

Form new input arrays
$\check{\mathbf{X}}_{Q}^{i}=\frac{\check{\mathbf{X}}^{0}+\mathbf{Y}_{Q}^{i}}{2},\quad\check{\mathbf{X}}_{W}^{i}=\frac{\check{\mathbf{X}}^{0}+\mathbf{Y}_{W}^{i}}{2}.$

  2. 2.

Apply the operators $\mathbf{Q}$ and $\mathbf{W}$ to the input arrays:
$\mathbf{Y}_{Q}^{i+1}=\mathbf{Q}\,\check{\mathbf{X}}_{W}^{i},\quad\mathbf{Y}_{W}^{i+1}=\mathbf{W}\,\check{\mathbf{X}}_{Q}^{i}.$

Estimations of the clean image:

Three estimations are used:

  1. 1.

The cross-boosted WNNM estimation
$\tilde{\mathbf{X}}_{uW}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{Y}_{W}^{I}$
(cbWNNM ).

  2. 2.

The cross-boosted qWPdn estimation
$\tilde{\mathbf{X}}_{uQ}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{Y}_{Q}^{I}$
(cbqWP).

  3. 3.

The hybrid estimation
$\tilde{\mathbf{X}}_{H}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(\mathbf{Y}_{W}^{I}+\mathbf{Y}_{Q}^{I})/2$
(hybrid).

###  3.3 Experimental results

In this section, we compare the performance of our denoising schemes
designated as cbWNNM , cbqWP and hybrid on the restoration of degraded images
with the performances of the state-of-the-art algorithms such as BM3D ([15]),
BM3D-SAPCA ([16]), WNNM ([22]), NCSR ([18]), cptTP-$\mathbb{C}$TF6 ([48]) and
DAS-2 ([10]). To produce results for the comparison, we used the software
available at the websites http://www.cs.tut.fi/ foi/GCF-
BM3D/index.html#ref_software (BM3D and BM3D-SAPCA),
http://staffweb1.cityu.edu.hk/xzhuang7/softs/index.html#bdTPCTF
(cptTP-$\mathbb{C}TF_{6}$ and DAS-2), https://github.com/csjunxu/WNNM_CVPR2014
(WNNM), and
https://www4.comp.polyu.edu.hk/ cslzhang/NCSR.htm (NCSR).

The restored images were evaluated by the visual perception, by Peak Signal-
to-Noise ratio (PSNR) (see Eq. (3.2))444
$PSNR(\mathbf{x},\mathbf{\tilde{x}})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}10\log_{10}\left(\frac{K\,255^{2}}{\sum_{k=1}^{K}(x_{k}-\tilde{x}_{k})^{2}}\right)\;dB.$
(3.2) and by the Structural Similarity Index (SSIM) ([43], it is computed by
the ssim.m Matlab 2020b function). The SSIM measures the structural similarity
of small moving windows in two images. It varies from 1 for fully identical
windows to -1 for completely dissimilar ones. The global index is the average
of local indices. Currently, SSIM is regarded as more informative
characteristics of the image quality compared to PSNR and Mean Square Error
(MSE) (see discussion in [42]).

For the experiments, we used a standard set of benchmark images: “Lena”,
“Boat”, “Hill”, “Barbara”, “Mandrill”, “Bridge”, “Man”, “Fabric” and
“Fingerprint”. One image that represents a stacked seismic section is
designated as “Seismic”. The ”clean” images are displayed in Fig. 3.1.

Figure 3.1: Clean images: “Lena”, “Boat”, “Hill, “Barbara”, “Mandrill”,
“Bridge”, “Man”, “Fabric”, “Fingerprint” and “Seismic”

The images were corrupted by Gaussian zero-mean noise whose STD was
$\sigma=$5, 10, 25, 40, 50, 80 and 100 dB. Then, the BM3D, BM3D-SAPCA, NCSR,
WNNM, cbWNNM , cbqWP, hybrid, cptTP-$\mathbb{C}TF_{6}$ and DAS-2 denoising
algorithms were applied to restore the images. In most experiments the
algorithm cbWNNM performed better than cbqWP. However, this was not the case
in experiments with the “Seismic” image. Therefore, in the “Seismic” block in
Table 3.2 and pictures in Fig. 3.7 we provide results from experiments with
the cbqWP rather than with the cbWNNM algorithm.

Table 3.1 summarizes experimental results from the restoration of the
“Barbara”, “Boat”, “Fingerprint”, “Lena” and “Mandrill” images corrupted by
additive Gaussian noise. PSNR and SSIM values for each experiment are given.

Table 3.1: PSNR/SSIM values from restoration of “Barbara”, “Boat”,
“Fingerprint”, “Lena” and “Mandrill” images. Boldface highlights the best
results. Noise STD=$\sigma$

$\sigma$ | 5 | 10 | 25 | 40 | 50 | 80 | 100
---|---|---|---|---|---|---|---
Barbara
noised | ${34.16}/{0.77}$ | ${28.14}/{0.6}$ | ${20.18}/{0.35}$ | ${16.09}/{0.23}$ | ${14.16}/{0.18}$ | ${10.07}/{0.09}$ | ${8.14}/{0.06}$
WNNM | $\mathbf{{38.77}}/{0.836}$ | $\mathbf{35.51}/{0.7733}$ | ${\mathbf{31.27}}/{0.68}$ | ${{28.77}/{0.5993}}$ | ${{27.8}}/{0.5628}$ | ${{25.48}/{0.4499}}$ | ${24.36}/{0.3891}$
SAPCA | ${38.37}/{0.8445}$ | ${35.1}/{0.7756}$ | ${{31}}/{0.6797}$ | ${{28.68}/{0.6086}}$ | ${{27.5}}/{0.5605}$ | ${24.3}/{0.402}$ | ${23.09}/{0.339}$
NCSR | ${38.37}/{0.8302}$ | ${35.4}/{0.7674}$ | ${30.62}/{0.6579}$ | ${28.19}/{0.5687}$ | ${27}/{0.5176}$ | ${24.37}/{0.3902}$ | ${23.21}/{0.3296}$
cptTP-$\mathbb{C}$TF6 | ${37.71}/{0.8382}$ | ${34.04}/{0.763}$ | ${29.24}/{0.6358}$ | ${26.75}/{0.5375}$ | ${25.61}/{0.4833}$ | ${23.44}/{0.3633}$ | ${22.55}/{0.3092}$
DAS-2 | ${37.7}/{0.8438}$ | ${33.98}/{0.764}$ | ${29.4}/{0.6422}$ | ${27.09}/{0.5595}$ | ${26.01}/{0.5121}$ | ${23.72}/{0.3987}$ | ${22.63}/{0.3398}$
BM3D | ${38.3}/{0.8387}$ | ${34.97}/{0.7738}$ | ${{30.71}}/{0.6707}$ | ${27.98}/{0.5788}$ | ${27.23}/{0.5389}$ | ${24.79}/{0.4192}$ | ${23.63}/{0.3608}$
cbWNNM | ${38.74}/{{0.853}}$ | ${{35.4}}/{0.7833}$ | ${31.11}/{\mathbf{0.6852}}$ | $\mathbf{{28.9}/{0.6171}}$ | $\mathbf{27.81}/{\mathbf{0.576}}$ | ${\mathbf{25.51}/{0.4681}}$ | $\mathbf{24.47}/{{0.4076}}$
Hybrid | ${{38.59}}/\mathbf{0.8543}$ | ${35.25}/{\mathbf{0.7859}}$ | ${30.87}/{0.6834}$ | ${28.67}/{0.6154}$ | ${27.59}/{{0.5745}}$ | ${{25.35}/\mathbf{0.4688}}$ | ${24.33}/{{\mathbf{0.4114}}}$
Boat
noised | ${34.16}/{0.79}$ | ${28.14}/{0.58}$ | ${20.18}/{0.29}$ | ${16.09}/{0.18}$ | ${14.16}/{0.13}$ | ${10.07}/{0.07}$ | ${8.14}/{0.05}$
WNNM | ${{37.35}}/{0.8042}$ | ${{34.07}}/{0.6803}$ | ${\mathbf{30.05}}/{0.5266}$ | ${\mathbf{27.96}}/{0.4317}$ | ${{26.98}}/{0.3833}$ | ${25.03}/{0.2853}$ | ${24.11}/{0.2389}$
SAPCA | ${\mathbf{37.51}}/\mathbf{0.8288}$ | ${\mathbf{34.1}}/{0.6925}$ | ${{30.03}}/{0.5282}$ | ${{27.92}}/{0.4384}$ | ${{26.89}}/{0.3931}$ | ${24.68}/{0.2993}$ | ${23.69}/{0.2609}$
NCSR | ${37.33}/{0.8101}$ | ${33.92}/{0.6848}$ | ${29.78}/{0.5053}$ | ${27.64}/{0.4003}$ | ${26.4}/{0.3551}$ | ${24.61}/{0.257}$ | ${23.69}/{0.217}$
cptTP-$\mathbb{C}$TF6 | ${36.86}/{0.8044}$ | ${33.34}/{0.6713}$ | ${29.12}/{0.504}$ | ${27.02}/{0.4044}$ | ${26.09}/{0.3567}$ | ${24.27}/{0.2626}$ | ${23.48}/{0.2232}$
DAS-2 | ${36.85}/{0.8294}$ | ${33.15}/{0.688}$ | ${28.88}/{0.5182}$ | ${26.78}/{0.4228}$ | ${25.81}/{0.3763}$ | ${23.84}/{0.2792}$ | ${22.9}/{0.2361}$
BM3D | ${{37.29}}/{0.8065}$ | ${{33.91}}/{0.6805}$ | ${{29.9}}/{0.5296}$ | ${27.74}/{0.4395}$ | ${26.77}/{0.3899}$ | ${{24.87}}/{0.2952}$ | ${23.96}/{0.2533}$
cbWNNM | ${37.42}/{{0.8257}}$ | ${34.03}/{0.702}$ | ${29.99}/{0.5455}$ | ${{27.95}}/{0.4568}$ | ${\mathbf{27}}/{0.4113}$ | ${\mathbf{25.09}}/{0.3116}$ | $\mathbf{24.26}/{{0.2638}}$
Hybrid | ${{37.13}/{0.8194}}$ | ${33.86}/{\mathbf{0.7105}}$ | ${29.72}/{\mathbf{0.5521}}$ | ${27.66}/{\mathbf{0.4616}}$ | ${26.78}/{\mathbf{0.4158}}$ | ${24.92}/{\mathbf{0.3182}}$ | ${24.14}/\mathbf{0.2725}$
Fingerprint
noised | ${34.16}/{0.97}$ | ${28.14}/{0.91}$ | ${20.18}/{0.67}$ | ${16.09}/{0.48}$ | ${14.16}/{0.38}$ | ${10.07}/{0.21}$ | ${8.14}/{0.15}$
WNNM | ${{35.16}/{0.9801}}$ | ${{32.68}/\mathbf{0.9655}}$ | $\mathbf{27.94}/{0.9014}$ | ${25.6}/{0.8388}$ | $\mathbf{24.7}/{0.8093}$ | ${22.77}/{0.7289}$ | ${21.85}/{0.6777}$
SAPCA | ${36.66}/{0.9859}$ | ${32.64}/{0.9648}$ | ${27.8}/{0.8981}$ | ${25.54}/{0.839}$ | ${24.53}/{0.8057}$ | ${22.3}/{0.7068}$ | ${21.18}/{0.6449}$
NCSR | ${{36.8}}/{0.9860}$ | ${\mathbf{32.69}}/{0.9645}$ | ${{27.83}}/{0.8968}$ | ${25.52}/{0.8290}$ | ${24.48}/{0.7894}$ | ${22.37}/{0.6869}$ | ${21.4}/{0.6295}$
cptTP-$\mathbb{C}$TF6 | ${36}/{0.9845}$ | ${32.23}/{0.9619}$ | ${27.33}/{0.889}$ | ${25.07}/{0.8219}$ | ${24.04}/{0.7815}$ | ${21.98}/{0.6746}$ | ${21.03}/{0.6129}$
DAS-2 | ${36.25}/{0.9843}$ | ${32.03}/{0.9601}$ | ${27.07}/{0.886}$ | ${24.86}/{0.8246}$ | ${23.88}/{0.7895}$ | ${21.87}/{0.7013}$ | ${20.96}/{0.652}$
BM3D | ${36.5}/{0.9854}$ | ${32.45}/{0.9634}$ | ${27.71}/{0.8955}$ | ${25.3}/{0.8334}$ | ${24.53}/{0.8019}$ | ${22.2}/{0.7165}$ | ${21.61}/{0.6643}$
cbWNNM | $\mathbf{36.91}/\mathbf{0.9864}$ | ${32.1}/{0.9615}$ | ${27.85}/{0.9039}$ | ${\mathbf{25.68}/{0.851}}$ | ${{24.67}}/\mathbf{0.8215}$ | ${\mathbf{22.79}/{0.7417}}$ | ${\mathbf{21.89}/{{0.7031}}}$
Hybrid | ${36.8}/{{0.9861}}$ | ${32.51}/{{0.9645}}$ | ${27.88}/{\mathbf{0.904}}$ | ${25.63}/\mathbf{{0.8512}}$ | ${24.59}/{{0.8213}}$ | ${22.73}/\mathbf{{0.7425}}$ | ${21.77}/\mathbf{0.7061}$
Lena
noised | ${34.16}/{0.65}$ | ${28.14}/{0.43}$ | ${20.18}/{0.2}$ | ${16.09}/{0.12}$ | ${14.16}/{0.09}$ | ${10.07}/{0.04}$ | ${8.14}/{0.03}$
WNNM | ${\mathbf{38.8}}/{0.7131}$ | ${{36.05}}/{0.6193}$ | ${\mathbf{32.25}}/{0.5036}$ | ${{30.1}}/{0.425}$ | ${{29.24}}/{0.3942}$ | ${27.22}/{0.3208}$ | ${26.19}/{0.278}$
SAPCA | ${{38.67}}/{0.6893}$ | ${\mathbf{36.07}}/{0.6233}$ | ${{32.23}}/{0.5077}$ | ${{30.11}}/{0.4411}$ | ${{29.07}}/{0.408}$ | ${26.5}/{0.314}$ | ${25.1}/{0.2771}$
NCSR | ${38.73}/{0.7123}$ | ${35.85}/{0.6214}$ | ${31.92}/{0.4841}$ | ${29.91}/{0.4155}$ | ${28.9}/{0.382}$ | ${26.72}/{0.3025}$ | ${25.71}/{0.2685}$
cptTP-$\mathbb{C}$TF6 | ${38}/{0.705}$ | ${35.46}/{0.616}$ | ${31.53}/{0.4992}$ | ${29.41}/{0.4313}$ | ${28.41}/{0.3965}$ | ${26.36}/{0.3199}$ | ${25.41}/{0.2833}$
DAS-2 | ${38.18}/{0.7272}$ | ${35.2}/{0.625}$ | ${31.09}/{0.4920}$ | ${28.9}/{0.4148}$ | ${27.84}/{0.3757}$ | ${25.53}/{0.29}$ | ${24.4}/{0.2285}$
BM3D | ${{38.72}}/{0.7078}$ | ${{35.92}}/{0.6233}$ | ${{32.07}}/{0.5026}$ | ${29.87}/{0.4265}$ | ${{29.04}}/{0.3957}$ | ${{26.98}}/{0.3214}$ | ${{25.95}}/{0.2851}$
cbWNNM | ${38.79}/{{0.7095}}$ | ${36}/{0.633}$ | ${32.20}/{{0.5155}}$ | ${\mathbf{30.2}}/{0.4479}$ | ${\mathbf{29.25}/{0.4139}}$ | ${\mathbf{27.24}/\mathbf{0.3359}}$ | ${\mathbf{26.36}/{0.304}}$
Hybrid | ${{38.76}}/\mathbf{0.7218}$ | ${35.89}/{\mathbf{0.6408}}$ | ${32.02}/{\mathbf{0.5202}}$ | ${30}/{\mathbf{0.4528}}$ | ${29.06}/{\mathbf{0.4193}}$ | ${27.07}/{{0.3349}}$ | ${26.21}/{\mathbf{0.312}}$
Mandrill
noised | ${34.16}/{0.7}$ | ${28.14}/{0.52}$ | ${20.18}/{0.51}$ | ${16.09}/{0.34}$ | ${14.16}/{0.26}$ | ${10.07}/{0.14}$ | ${8.14}/{0.1}$
WNNM | ${{34.7}}/{0.8902}$ | ${{30.36}}/{0.7994}$ | ${{25.44}}/{0.605}$ | ${\mathbf{23.43}}/{0.4739}$ | ${\mathbf{22.58}}/{0.4024}$ | ${\mathbf{21.1}}/{0.2645}$ | $\mathbf{20.46}/{0.1949}$
SAPCA | ${\mathbf{35.2}}/{0.9281}$ | ${\mathbf{30.59}}/{0.8278}$ | ${\mathbf{25.54}}/{0.6244}$ | ${{23.4}}/{0.4843}$ | ${{22.54}}/{0.413}$ | ${{20.92}}/{0.2649}$ | ${20.3}/{0.2048}$
NCSR | ${35.07}/{0.9129}$ | ${30.38}/{0.7986}$ | ${25.36}/{0.5908}$ | ${23.18}/{0.4325}$ | ${22.35}/{0.3669}$ | ${20.82}/{0.2215}$ | ${20.23}/{0.162}$
cptTP-$\mathbb{C}$TF6 | ${35.06}/{0.9252}$ | ${30.32}/{0.8198}$ | ${25.3}/{0.6082}$ | ${23.16}/{0.4531}$ | ${22.25}/{0.371}$ | ${20.72}/{0.2148}$ | ${20.2}/{0.1617}$
DAS-2 | ${35.02}/{\mathbf{0.9301}}$ | ${30.24}/{{0.8329}}$ | ${25.24}/{0.64}$ | ${23.2}/{0.51}$ | ${22.34}/{0.4423}$ | ${20.8}/{0.2997}$ | ${20.14}/{0.2376}$
BM3D | ${{{34.98}}}/{0.9209}$ | ${{30.34}}/{0.8135}$ | ${{25.27}}/{0.6095}$ | ${23}/{0.4613}$ | ${{22.27}}/{0.3813}$ | ${{20.9}}/{0.2439}$ | ${{20.39}}/{0.1956}$
cbWNNM | ${34.34}/{{0.9158}}$ | ${{30.39}}/{\mathbf{0.8359}}$ | ${25.53}/{\mathbf{0.6631}}$ | ${{23.35}}/\mathbf{0.5477}$ | ${{{22.39}}/{0.4802}}$ | ${{{20.68}}/{0.3513}}$ | ${{20.03}/{0.2841}}$
Hybrid | ${{33.23}}/{0.9028}$ | ${30.02}/{{0.8303}}$ | ${25.4}/{{0.6625}}$ | ${23.14}/{{0.5471}}$ | ${22.28}/{\mathbf{0.4887}}$ | ${20.31}/{\mathbf{0.3565}}$ | ${19.59}/{\mathbf{0.2945}}$

Table 3.2 summarizes experimental results from the restoration of the “Hill”,
“Seismic”, “Fabric”, “Bridge” and“Man” images corrupted by additive Gaussian
noise. The PSNR and SSIM values for each experiment are given.

Table 3.2: PSNR/SSIM values from restoration of “Hill”, “Seismic”, “Fabric”,
“Bridge” and“Man” images. Boldface highlights the best results. Noise
STD=$\sigma$.

$\sigma$ | 5 | 10 | 25 | 40 | 50 | 80 | 100
---|---|---|---|---|---|---|---
Hill
noised | ${34.16}/{0.8}$ | ${28.14}/{0.59}$ | ${20.18}/{0.27}$ | ${16.09}/{0.15}$ | ${14.16}/{0.11}$ | ${10.07}/{0.05}$ | ${8.14}/{0.03}$
WNNM | ${{37.02}}/{0.8239}$ | ${{33.7}}/{0.7121}$ | ${{29.96}}/{0.5225}$ | ${\mathbf{28.18}}/{0.4203}$ | ${27.34}/{0.373}$ | ${25.63}/{0.2822}$ | ${24.75}/{0.236}$
SAPCA | ${\mathbf{37.3}}/{0.8457}$ | ${\mathbf{33.84}}/{0.7316}$ | ${{29.95}}/{0.5335}$ | ${{28.09}}/{0.4281}$ | ${27.2}/{0.3815}$ | ${25.18}/{0.2867}$ | ${24.29}/{0.2491}$
NCSR | ${37.17}/{0.8350}$ | ${33.7}/{0.7246}$ | ${29.72}/{0.5138}$ | ${27.84}/{0.3936}$ | ${26.99}/{0.3463}$ | ${25.16}/{0.2480}$ | ${24.35}/{0.2089}$
cptTP-$\mathbb{C}$TF6 | ${36.83}/{0.8388}$ | ${33.17}/{0.7157}$ | ${29.13}/{0.4992}$ | ${27.37}/{0.3903}$ | ${26.59}/{0.3425}$ | ${25}/{0.2495}$ | ${24.29}/{0.2108}$
DAS-2 | ${36.68}/{{0.8447}}$ | ${33.02}/{0.7292}$ | ${29.15}/{0.5327}$ | ${27.31}/{0.4245}$ | ${26.47}/{0.3747}$ | ${24.62}/{0.2733}$ | ${23.7}/{0.2272}$
BM3D | ${37.14}/{0.8384}$ | ${{33.61}}/{0.7193}$ | ${{29.85}}/{0.5311}$ | ${27.99}/{0.4305}$ | ${{27.19}}/{0.3825}$ | ${{25.42}}/{0.2878}$ | ${{24.59}}/{0.2448}$
cbWNNM | $37.2/{{0.8515}}$ | ${33.72}/0.7461$ | $\mathbf{30.01}/{0.5545}$ | ${{28.13}}/{0.4584}$ | ${{\mathbf{27.4}}/{0.4033}}$ | ${\mathbf{25.64}}/{0.3148}$ | ${{\mathbf{24.87}}/{0.2718}}$
Hybrid | ${{37.09}}/\mathbf{0.8516}$ | ${33.64}/{\mathbf{0.7511}}$ | ${29.86}/{\mathbf{0.5655}}$ | ${28.01}/{\mathbf{0.4642}}$ | ${27.29}/{\mathbf{0.4118}}$ | ${{25.62}}/{\mathbf{0.3215}}$ | ${24.84}/{\mathbf{0.278}}$
Seismic
noised | ${34.16}/{0.8}$ | ${28.14}/{0.55}$ | ${20.18}/{0.22}$ | ${16.09}/{0.12}$ | ${14.16}/{0.08}$ | ${10.07}/{0.04}$ | ${8.14}/{0.03}$
WNNM | ${39.26}/{0.9205}$ | ${34.95}/{0.765}$ | ${30.65}/{0.4556}$ | ${28.93}/{0.3092}$ | ${28.34}/{0.2524}$ | $\mathbf{27.11}/{0.1621}$ | $\mathbf{26.56}/{0.1305}$
SAPCA | ${39.09}/{0.9167}$ | ${35.04}/{0.7871}$ | ${30.8}/{0.5016}$ | ${29.04}/{0.3654}$ | ${28.2}/{0.3124}$ | ${26.34}/{0.23}$ | ${25.4}/{0.1934}$
NCSR | ${38.97}/{0.9131}$ | ${34.86}/{0.78}$ | ${30.62}/{0.4643}$ | ${29}/{0.2948}$ | ${28.27}/{0.2374}$ | ${26.81}/{0.139}$ | ${26.25}/{0.1089}$
cptTP-$\mathbb{C}$TF6 | ${38.81}/{0.9107}$ | ${34.83}/{0.7696}$ | ${30.72}/{0.4616}$ | ${29.1}/{0.3188}$ | ${28.4}/{0.2618}$ | ${26.91}/{0.1616}$ | ${26.16}/{0.127}$
DAS-2 | ${38.8}/{0.9141}$ | ${34.83}/{0.795}$ | ${30.67}/{0.5383}$ | ${28.87}/{0.3938}$ | ${28.05}/{0.3318}$ | ${26.22}/{0.2222}$ | ${25.29}/{0.1785}$
BM3D | ${39}/{0.9126}$ | ${{34.9}}/{0.77}$ | ${{30.8}}/{0.4984}$ | ${29.08}/{0.3606}$ | ${\mathbf{28.45}}/{0.2853}$ | ${{26.97}}/{0.19}$ | ${{26.31}}/{0.1547}$
cbqWP | $\mathbf{{39.57}/{0.9318}}$ | $\mathbf{{35.41}/{0.8321}}$ | ${30.89}/{\mathbf{0.5849}}$ | ${{28.98}}/{\mathbf{0.445}}$ | ${{27.97}}/{\mathbf{0.3836}}$ | ${25.42}/{\mathbf{0.2755}}$ | ${23.36}/{\mathbf{0.2299}}$
Hybrid | ${{{39.13}}/{0.9241}}$ | ${{35.41}/{0.8266}}$ | ${\mathbf{30.99}}/{{0.5745}}$ | ${\mathbf{29.26}}/{{0.4298}}$ | ${{28.4}}/{{0.369}}$ | ${26.28}/{{0.2657}}$ | ${24.75}/{{{0.2227}}}$
Fabric
noised | ${34.16}/{0.8}$ | ${28.14}/{0.55}$ | ${20.18}/{0.22}$ | ${16.09}/{0.12}$ | ${14.16}/{0.08}$ | ${10.07}/{0.04}$ | ${8.14}/{0.03}$
WNNM | ${\mathbf{38.75}/{0.833}}$ | ${{34.92}/{0.7318}}$ | ${{30.76}/{0.5488}}$ | ${{28.83}/{0.4485}}$ | ${\mathbf{28.01}}/{0.4113}$ | $\mathbf{26.19}/{0.3361}$ | $\mathbf{25.26}/{0.2992}$
SAPCA | ${38.72}/{0.8351}$ | $\mathbf{{35}/{0.7432}}$ | $\mathbf{{30.8}/{0.5678}}$ | ${{28.82}/{\mathbf{0.4726}}}$ | ${{27.89}}/{0.3905}$ | ${25.54}/{0.3275}$ | ${24.51}/{0.2915}$
NCSR | ${38.58}/{0.8263}$ | ${34.78}/{0.7298}$ | ${30.53}/{0.5408}$ | ${28.54}/{0.4384}$ | ${27.69}/{0.3988}$ | ${25.77}/{0.3185}$ | ${24.8}/{0.2838}$
cptTP-$\mathbb{C}$TF6 | ${37.99}/{0.81}$ | ${34.07}/{0.6978}$ | ${29.73}/{0.5127}$ | ${27.83}/{0.4266}$ | ${26.94}/{0.3905}$ | ${25.12}/{0.3185}$ | ${24.23}/{0.2856}$
DAS-2 | ${37.74}/{0.8017}$ | ${33.78}/{0.6985}$ | ${29.48}/{0.5223}$ | ${27.48}/{0.4283}$ | ${26.52}/{0.3858}$ | ${24.45}/{0.3017}$ | ${23.42}/{0.2641}$
BM3D | ${38.47}/{0.8303}$ | ${{34.63}}/{0.7262}$ | ${{30.57}}/{0.5476}$ | ${{28.65}}/{0.4548}$ | ${{27.85}}/{0.4182}$ | ${{26}}/{0.3402}$ | ${{25.06}}/{0.3032}$
cbWNNM | ${{38.69}/\mathbf{0.8354}}$ | ${{34.92}/{0.7426}}$ | ${30.74}/{{0.5649}}$ | ${\mathbf{28.84}}/{\mathbf{0.4726}}$ | ${27.96}/{{0.4279}}$ | ${26.15}/{{0.3489}}$ | $\mathbf{25.26}/{\mathbf{0.3128}}$
Hybrid | ${{38.52}}/{0.8294}$ | ${34.76}/{{0.7349}}$ | ${{30.52}}/{0.561}$ | ${{28.61}}/{0.4718}$ | ${{27.73}}/\mathbf{0.4283}$ | ${25.9}/{\mathbf{0.352}}$ | ${25.01}/{{0.3096}}$
Bridge
noised | ${34.16}/{0.7}$ | ${28.14}/{0.52}$ | ${20.18}/{0.51}$ | ${16.09}/{0.34}$ | ${14.16}/{0.26}$ | ${10.07}/{0.14}$ | ${8.14}/{0.1}$
WNNM | ${35.81}/{0.9387}$ | ${31.17}/{0.854}$ | ${26.29}/{0.635}$ | ${{24.45}}/{0.4992}$ | ${23.67}/{0.43}$ | $\mathbf{22.24}/{0.3095}$ | ${21.6}/{0.2568}$
SAPCA | $\mathbf{{35.89}/{0.9439}}$ | ${\mathbf{31.35}}/{0.8671}$ | ${26.44}/{0.6587}$ | ${\mathbf{24.5}}/{0.5152}$ | $\mathbf{23.7}/{0.4465}$ | ${22.17}/{0.3273}$ | ${21.48}/{0.2783}$
NCSR | ${35.78}/{0.9375}$ | ${31.2}/{0.8551}$ | ${26.29}/{0.6389}$ | ${24.31}/{0.4814}$ | ${23.54}/{0.4169}$ | ${22.02}/{0.2832}$ | ${21.36}/{0.2344}$
cptTP-$\mathbb{C}$TF6 | ${35.56}/{0.9388}$ | ${30.94}/{0.8575}$ | ${26}/{0.6362}$ | ${24.01}/{0.4803}$ | ${23.22}/{0.4067}$ | ${21.78}/{0.2744}$ | ${21.17}/{0.2263}$
DAS-2 | ${35.44}/{{0.9373}}$ | ${30.76}/{{0.8573}}$ | ${25.91}/{0.6609}$ | ${24.03}/{0.5287}$ | ${23.22}/{0.4625}$ | ${21.66}/{0.3316}$ | ${20.94}/{0.2768}$
BM3D | ${{{35.78}}}/{0.9415}$ | ${{31.17}}/{0.8597}$ | ${{26.22}}/{0.641}$ | ${24.3}/{0.5007}$ | ${{23.58}}/{0.4286}$ | ${{22.22}}/{0.3168}$ | ${\mathbf{21.61}}/{0.2716}$
cbWNNM | ${35.77}/{0.9412}$ | ${31.28}/{\mathbf{0.869}}$ | $\mathbf{26.45}/{{0.6838}}$ | ${{24.49}}/{0.5587}$ | ${{{23.68}}/{0.4963}}$ | ${22.05}/{0.3717}$ | ${{21.37}}/{0.3158}$
Hybrid | ${{35.45}}/{0.9362}$ | ${31.12}/{{0.8643}}$ | ${26.23}/{\mathbf{0.684}}$ | ${24.24}/{\mathbf{0.5659}}$ | ${23.32}/{\mathbf{0.5033}}$ | ${21.57}/{\mathbf{0.3822}}$ | ${20.77}/{\mathbf{0.3292}}$
Man
noised | ${34.16}/{0.75}$ | ${28.14}/{0.55}$ | ${20.18}/{0.27}$ | ${16.09}/{0.16}$ | ${14.16}/{0.12}$ | ${10.07}/{0.06}$ | ${8.14}/{0.04}$
WNNM | ${{37.99}/{0.8267}}$ | ${{34.19}}/{0.7234}$ | ${{29.79}}/{0.5403}$ | ${{27.8}}/{0.4326}$ | ${{26.95}}/{0.3842}$ | $\mathbf{25.18}/{0.2911}$ | $\mathbf{24.34}/{0.2455}$
SAPCA | $\mathbf{38.09}/{0.8387}$ | ${\mathbf{34.28}}/{0.7325}$ | ${\mathbf{29.84}}/{0.5488}$ | ${\mathbf{27.86}}/{0.445}$ | ${\mathbf{26.97}}/{0.3997}$ | ${24.92}/{0.3078}$ | ${23.9}/{0.2677}$
NCSR | ${37.88}/{0.8216}$ | ${34.08}/{0.7212}$ | ${29.62}/{0.5289}$ | ${27.58}/{0.4107}$ | ${26.7}/{0.3632}$ | ${24.89}/{0.2671}$ | ${24.05}/{0.2276}$
cptTP-$\mathbb{C}$TF6 | ${37.44}/{0.8302}$ | ${33.59}/{0.7235}$ | ${29.14}/{0.5341}$ | ${27.14}/{0.4241}$ | ${26.25}/{0.3725}$ | ${23.99}/{0.2713}$ | ${23.73}/{0.229}$
DAS-2 | ${37.1}/{{0.7525}}$ | ${33.16}/{{0.7184}}$ | ${28.81}/{0.5374}$ | ${26.84}/{0.4326}$ | ${25.93}/{0.3815}$ | ${23.99}/{0.2783}$ | ${23.06}/{0.2336}$
BM3D | ${{37.85}}/{{0.8306}}$ | ${{34.02}}/{0.7225}$ | ${{29.65}}/{0.5429}$ | ${27.68}/{0.442}$ | ${{26.84}}/{0.3919}$ | ${{25.09}}/{0.3009}$ | ${{22.26}}/{0.2598}$
cbWNNM | ${37.97}/{{0.8422}}$ | ${{34.19}}/{{0.7428}}$ | ${29.81}/{0.5646}$ | ${{27.81}}/{0.4622}$ | ${{{26.82}}/{0.4173}}$ | ${{{25.04}}/{0.3216}}$ | ${{{24.11}}/{0.2806}}$
Hybrid | ${{37.75}}/\mathbf{0.8426}$ | ${34.02}/{\mathbf{0.7452}}$ | ${29.65}/{{\mathbf{0.5674}}}$ | ${27.61}/{\mathbf{0.4664}}$ | ${26.63}/{\mathbf{0.4195}}$ | ${24.84}/{\mathbf{0.3257}}$ | ${23.88}/{\mathbf{0.285}}$

Table 3.3 provides the PSNR and SSIM values from Tables 3.1 and 3.2, which are
averaged over ten images participating in the experiments. Respective diagrams
are drawn in Fig. 1.1.

Table 3.3: PSNR/SSIM values averaged over 10 images. Boldface highlights the
highest values. Noise STD=$\sigma$.

$\sigma$ | 5 | 10 | 25 | 40 | 50 | 80 | 100
---|---|---|---|---|---|---|---
WNNM | ${\mathbf{37.63}}/{0.8649}$ | ${\mathbf{33.81}}/{0.768}$ | ${\mathbf{29.46}}/{0.5949}$ | ${{27.34}}/{0.489}$ | $\mathbf{26.56}/{0.4405}$ | $\mathbf{24.8}/{0.3405}$ | $\mathbf{23.95}/{0.2916}$
SAPCA | ${{37.57}}/{0.8685}$ | ${{33.8}}/{0.7745}$ | ${{29.44}}/{0.6048}$ | ${{27.4}}/{0.5038}$ | ${26.45}/{0.4511}$ | ${24.29}/{0.3466}$ | ${23.29}/{0.3}$
NCSR | ${37.44}/{0.8586}$ | ${33.65}/{0.765}$ | ${29.2}/{0.5819}$ | ${27.13}/{0.4664}$ | ${26.2}/{0.4163}$ | ${24.33}/{0.3105}$ | ${23.48}/{0.2662}$
cptTP-$\mathbb{C}$TF6 | ${37.03}/{0.8586}$ | ${33.2}/{0.7596}$ | ${28.72}/{0.578}$ | ${26.69}/{0.4688}$ | ${25.78}/{0.4163}$ | ${23.96}/{0.311}$ | ${23.23}/{0.2669}$
DAS-2 | ${36.98}/{{0.8565}}$ | ${33.02}/{0.7668}$ | ${28.56}/{0.597}$ | ${26.54}/{0.494}$ | ${25.61}/{0.4432}$ | ${23.67}/{0.3376}$ | ${22.74}/{0.2874}$
BM3D | ${37.4}/{0.8613}$ | ${{33.59}}/{0.7652}$ | ${{29.28}}/{0.5969}$ | ${27.15}/{0.4928}$ | ${{26.36}}/{0.4414}$ | ${{24.54}}/{0.3432}$ | ${{23.54}}/{0.2993}$
cbWNNM | $37.54/{\mathbf{0.8692}}$ | ${33.74}/0.7848$ | $\mathbf{29.46}/{0.6266}$ | ${\mathbf{27.43}}/{0.5317}$ | ${{{26.5}}/{0.4831}}$ | ${{24.56}}/{0.3841}$ | ${{{23.6}}/{0.3373}}$
Hybrid | ${{37.25}}/{0.8668}$ | ${33.65}/{\mathbf{0.7854}}$ | ${29.31}/{\mathbf{0.6275}}$ | ${27.28}/{\mathbf{0.5326}}$ | ${26.37}/{\mathbf{0.4852}}$ | ${{24.46}}/{\mathbf{0.3868}}$ | ${23.53}/{\mathbf{0.3421}}$

It is seen from Tables 3.1, 3.2 and 3.3 that our methods cbWNNM and hybrid
produce PSNR values very close to (sometimes higher than) those produced by
the WNNM and BM3D-SAPCA methods and higher than those from BM3D, NCSR,
cptTP-$\mathbb{C}TF_{6}$ and DAS-2 denoising algorithms. On the other hand,
the SSIM values for the images restored by cbWNNM and hybrid significantly
exceed the SSIM values from all other methods. Especially it is true for the
restoration of images in presence of a strong noise. This observation reflects
the fact that directional qWPs have exceptional capabilities to capture fine
structures even in severely degraded images. This fact is illustrated by
Figures 3.2, 3.3, 3.4, 3.5, 3.6 and 3.7, which display results of restoration
of several images corrupted by strong Gaussian noise. In all those cases the
SSIM values from the cbWNNM and hybrid algorithms significantly exceed values
from all other algorithms. Respectively, the restoration of the images’ fine
structure by the cbWNNM and hybrid algorithms is much better compared to the
restoration by other algorithms.

Each of the mentioned figures comprises 12 frames, which are arranged in a
$4\times 3$ order: $\left(\begin{array}[]{ccc}f_{11}&f_{12}&f_{13}\\\
f_{21}&f_{22}&f_{23}\\\ f_{31}&f_{32}&f_{33}\\\ f_{41}&f_{42}&f_{43}\\\
\end{array}\right).$ Here frame $f_{11}$ displays noised image; frame $f_{21}$
– image restored by BM3D; $f_{12}$ – image restored by BM3D-SAPCA; $f_{22}$ –
image restored by WNNM; $f_{13}$ – image restored by cbWNNM555For the
“Seismic” image cbqWP instead of cbWNNM. ; $f_{23}$ – image restored by
hybrid. Frame $f_{31}$ displays a fragment of the original image. The
remaining frames $\left\\{f_{32},f_{33},f_{41},f_{42},f_{43}\right\\}$ display
the fragments of the restored images shown in frames
$\left\\{f_{12},f_{13},f_{21},f_{22},f_{23}\right\\}$, which are arranged in
the same order.

Figure 3.2: Restoration of the “Hill” image corrupted by Gaussian noise with
STD $\sigma=50$ dB. PSNR/SSIM for BM3D– 27.19/0.3825; for WNNM– 27.34/0.373;
for BM3D-SAPCA–27.2/0.3815; for cbWNNM–27.4/0.4033; for hybrid–27.29/0.4118

Figure 3.3: Restoration of “Barbara” image corrupted by Gaussian noise with
STD $\sigma=80$ dB. PSNR/SSIM for BM3D– 24.79/0.4192; for WNNM– 25.48/0.4499;
for BM3D-SAPCA–24.33/0.402; for cbWMMN–25.51/0.4681; for hybrid–25.35/0.4688

Figure 3.4: Restoration of “Fingerprint” image corrupted by Gaussian noise
with STD $\sigma=80$ dB. PSNR/SSIM for BM3D– 22.2/0.7165; for WNNM–
22.77/0.7289; for BM3D-SAPCA–22.3/0.7068; for upWMMN–22.79/0.7417; for
hybrid–22.73/0.7425

Figure 3.5: Restoration of “Mandrill” image corrupted by Gaussian noise with
STD $\sigma=80$ dB. PSNR/SSIM for BM3D– 20.9/0.2439; for WNNM– 21.1/0.2645;
for BM3D-SAPCA–20.92/0.2649; for cbWNNM–20.68/0.3513; for hybrid–20.31/0.3565

Figure 3.6: Restoration of “Bridge” image corrupted by Gaussian noise with STD
$\sigma=40$ dB. PSNR/SSIM for BM3D– 24.3/0.5007; for WNNM– 24.45/0.4992; for
BM3D-SAPCA–24.5/0.5152; for cbWNNM–24.49/0.5587; for hybrid–24.24/0.5659

Figure 3.7: Restoration of “Seismic” image corrupted by Gaussian noise with
STD $\sigma=25$ dB. PSNR/SSIM for BM3D– 30.8/0.4984; for WNNM– 30.65/0.4556;
for BM3D-SAPCA–30.8/0.5016; for cbqWP–30.89/0.5849; for hybrid–30.99/0.5745

###  3.4 Relation of the proposed algorithms to the Deep Learning methods

In recent years, the focus of the image denoising research shifted to the Deep
Learning methods, which resulted in a huge amount of publications. We mention
a few of them: [11, 45, 13, 19, 39, 32, 35, 40]. Much more references can be
found in the reviews [41, 26]. One of advantages of the Deep Learning methods
is that, once the Neural Net is trained (which can involve extended datasets
and take several days), its application to test images is very fast.
Therefore, the experimental results in most related publications are presented
via the PSNR and, sometimes, the SSIM values averaged over some test datasets
such as, for example, Set12 introduced in [45].

Set12 partially overlap with the set of 10 (Set10) images that we used in our
experiments. Namely the images “Barbara”, Boat”, “Fingerprint”, “Hill”, “Lena”
and “Man” participate in both sets. The structure of the remaining images
“Seismic”, “Fabric”, “Mandrill” and “Bridge” from Set10 is more complicated
compared to the images “Camera”, “Couple”, “House”, “Monarch”, “Pepper” and
“Straw” from Set12. Therefore, the averaged results comparison from these two
datasets is quite justified. For even better compatibility, we compare the
gains of results from different methods over the corresponding results from
BM3D:
$P_{method}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{PSNR_{method}}{PSNR_{BM3D}}$
and
$S_{method}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{SSIM_{method}}{SSIM_{BM3D}}$.
Recall that for the calculation of SSIM we use the function ssim from Matlab
2020b, whereas in most publications SSIM is computed by some other schemes.

We compare results from the recent state-of-the-art algorithms Cola_Net ([34],
(2022)), CDNet ([35], (2021)), FLCNN ([40], (2021)), DRCNN ([32], (2020)), and
the non-Deep Learning algorithm presented in [46] (2021), which we mark as
SRENSS, averaged over Set12 with the results from WNNM [22] and the proposed
hybrid algorithm averaged over Set10666Averaged results from upBM3D are almost
identical to those from hybrid algorithm. The PSNR and SSIM values for all
methods except for WNNM and hybrid are taken from the corresponding
publications. Table 3.4 shows the results of this comparison.

Table 3.4: $P_{method}/S_{method}$, where the PSNR and SSIM values are
averaged over either Set12 or Set10. Noise STD=$\sigma$.

$\sigma$ | 10 | 15 | 25 | 50 | 70 | 80 | 100
---|---|---|---|---|---|---|---
Cola_Net | - | 1.028/1.016 | 1.031/1.025 | 1.04/1.046 | 1.041/1.06 | - | -
CDNet | - | 1.019/1.016 | 1.016/1.01 | 1.025/1.074 | 1.027/1.057 | 1.027/1.065 | -
FLCNN | - | 1.022/1.02 | 1.017/1.011 | 1.028/1.039 | - | - | -
DRCNN | - | 1.016/1.036 | 1.02/1.033 | 1.02/1.039 | - | - | -
SRENSS | 1.044/1.011 | - | - | 1.02/1.028 | 1.011/1.028 | - | 1.006/1.0038
WNNM | 1.007/1.004 | - | 1.006/0.997 | 1.007/0.998 | - | 1.01/0.992 | 1.018/0.974
hybrid | 1.002/1.026 | - | 1.008/1.051 | 1.001/1.099 | - | 0.997/1.127 | 0.999/1.143

We can observe from the table that all participated up-to-date schemes,
including the non-Deep Learning algorithm SRENSS, demonstrate a moderate gain
over BM3D in both PSNR and SSIM values averaged over Set12 (far from being a
breakthrough). The values from WNNM averaged over Set10 are very close to
those from BM3D averaged over the same set. The same can be said for the PSNR
values from the hybrid algorithm. However, the SSIM values from the hybrid
algorithm demonstrate a strong gain over BM3D, which is clearly seen in Fig.
1.1 and Table 3.3. Especially it is true for a strong Gaussian noise with
$\sigma=80,\,100$ dB. This fact highlights the ability of the qWP-based
algorithms to restore edges and fine structures even in severely damaged
images.

Note that denoising results presented in the overwhelming majority of the Deep
Learning publications deal with the noise level not exceeding 50 dB.

##  4 Discussion

We presented a denoising scheme that combines the qWPdn algorithm based on the
directional quasi-analytic wavelet packets, which are designed in [3], with
the state-of-the-art denoising algorithm WNNM ([22]) considered to be one of
the best in the field. Either of the two algorithms has their strong features
and shortcomings. The qWPdn algorithm described in Section 3.1 demonstrates
the ability to restore edges and texture details even from severely degraded
images. This ability stems from the fact that the designed 2D qWP transforms
provide a variety of 2D testing waveforms, which are close to windowed cosine
waves with multiple frequencies oriented in multiple directions. In most
separate experiments, the qWPdn method demonstrated better resolution of edges
and fine structures compared to the WNNM algorithm, which were reflected by
the higher SSIM values. In turn, the WNNM algorithm is superior for noise
suppression, especially in smooth regions in images, thus producing the
highest PSNR values in almost all the experiments. However, some over-
smoothing effect on the edges and fine texture persisted with the WNNM
algorithm when noise was strong.

The qWPdn and WNNM methods complement each other. Therefore, the idea to
combine these methods is natural. In the iterative hybrid scheme qWPdn–WNNM,
which is proposed in Section 3.2, the output from one algorithm updates the
input to the other. Such a hybrid method has some distant relation to the SOS
boosting scheme presented in [36]. The main distinction between the qWPdn–WNNM
and the SOS boosting is that each of the qWPdn and WNNM algorithms is
“boosted” by the output from the other algorithm. Such a scheme can be
regarded as a _Cross-Boosting_.

The scheme proved to be highly efficient. It is confirmed by a series of
experiments on restoration of 10 multimedia images of various structure, which
were degraded by Gaussian noise of various intensity. In the experiments, the
performance of two combined qWPdn–WNNM algorithms was compared with the
performance of six advanced denoising algorithms BM3D, BM3D-SAPCA ([16]), WNNM
([22]), NCSR ([18]), cptTP-$\mathbb{C}$TF6 ([48]) and DAS-2 ([10]). In almost
all the experiments reported in Section 3.3, the two combined algorithms
produce PSNR values, which are very close to the values produced by WNNM and
BM3D-SAPCA. Their noise suppression efficiency is competitive with that of
WNNM and BM3D-SAPCA. On the other hand, their results in the resolution of
edges and fine structures are much better than the results from all other
algorithms participating in the experiments. This is seen in the images
presented in Section 3.3. Consequently, the SSIM values produced by the cross-
boosted algorithms qWPdn–WNNM are significantly higher than the values
produced by all other participated algorithms.

Disscussion in Section 3.4 shows that the qWPdn–WNNM algorithms can, in some
aspects, compete with the Deep Learning denoising methods.

##### Acknowledgment

This research was supported by the Israel Science Foundation (ISF, 1556/17,
1873/21), Len Blavatnik and the Blavatnik Family Foundation, Israel Ministry
of Science Technology and Space 3-16414 and 3-14481.

## References

  * [1] A. Averbuch, R. R. Coifman, D. L. Donoho, M. Israeli, and Y. Shkolnisky.  A framework for discrete integral transformations I – the pseudopolar Fourier transform.  SIAM Journal on Scientific Computing, 30(2):764–784, 2008.
  * [2] A. Averbuch, R. R. Coifman, D. L. Donoho, M. Israeli, Y. Shkolnisky, and I. Sedelnikov.  A framework for discrete integral transformations II – the 2d discrete Radon transform.  SIAM Journal on Scientific Computing, 30(2):785–803, 2008.
  * [3] A. Averbuch, P. Neittaanmäki, and V. Zheludev.  Directional wavelet packets originating from polynomial splines.  http://arxiv.org/abs/2008.05364.
  * [4] A. Averbuch, P. Neittaanmäki, and V. Zheludev.  Splines and spline wavelet methods with application to signal and image processing, Volume III: Selected topics.  Springer, 2019.
  * [5] A. Averbuch, P. Neittaanmäki, V. Zheludev, M. Salhov, and J. Hauser.  Image inpainting using directional wavelet packets originating from polynomial splines.  Signal Processing: Image Communication, 97, 2021.  http://arxiv.org/abs/2001.04899.
  * [6] I. Bayram and I. W. Selesnick.  On the dual-tree complex wavelet packet and m-band transforms.  IEEE Trans. Signal Process., 56:2298–2310, 2008.
  * [7] A. Buades, B. Coll, and J.-M. Morel.  A review of image denoising algorithms, with a new one.  Multiscale Model. Simul., 4(2):490–530, 2005.
  * [8] E. Candés, L. Demanet, D. Donoho, and L. X. Ying.  Fast discrete curvelet transforms.  Multiscale Model. Simul., 5:861–899, 2006.
  * [9] E. Candés and D. Donoho.  New tight frames of curvelets and optimal representations of objects with piecewise $c^{2}$ singularities.  Commun. Pure Appl. Math., 57:219–266, 2004.
  * [10] Z. Che and X. Zhuang.  Digital affine shear filter banks with 2-layer structure and their applications in image processing.  IEEE Trans. on Image Processing, 27(8):3931–3941, 2018.
  * [11] Y. Chen and T. Pock.  Trainable nonlinear reaction diffusion: a flexible framework for fast and effective image restoration.  IEEE Trans. on Pattern analysis and Machine Intelligence, 39(6):1256–1272, 2017.
  * [12] R. R. Coifman and V. M. Wickerhauser.  Entropy-based algorithms for best basis selection.  IEEE Trans. Inform. Theory, 38(2):713–718, 1992.
  * [13] Cristóvão Cruz, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian.  Nonlocality-reinforced convolutional neural networks for image denoising.  IEEE Signal Processing Letters, 25(8):1216–1220, 2018.
  * [14] L. Şendur and I. W. Selesnick.  Bivariate shrinkage functions for wavelet-based denoising exploiting interscale dependency.  IEEE Trans. Signal Process., 50:2744–2756, 2002.
  * [15] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian.  Image denoising by sparse 3d transform-domain collaborative filtering.  IEEE Trans. Image Process., 16(8):2080––2095, 2007.
  * [16] K. Dabov, A. Foi, V. Katkovnik, and K.O. Egiazarian.  BM3D image denoising with shape adaptive principal component analysis.  In Proceedings of the Workshop on Signal Processing with Adaptive Sparse Structured Representations (SPARS’09), 2009.
  * [17] M. N. Do and M. Vetterli.  Contourlets.  In Beyond Wavelets, G. V. Welland, ed. Academic Press, San Diego, CA, 2008.
  * [18] W.S. Dong, L. Zhang, G.M. Shi, and X. Li.  Nonlocally centralized sparse representation for image restoration.  IEEE Trans. Image Process., 22(4):1620–1630, 2013.
  * [19] Y. Fang and T. Zeng.  Learning deep edge prior for image denoising.  Computer Vision and Image Understanding, 200:103044, 2020.
  * [20] W.-Q Lim G. Kutyniok and X. Zhuang.  Digital shearlet transforms.  In Shearlets: Multiscale Analysis for Multivariate Data, pages 239–282. Birkhäuser, Boston, 2012.
  * [21] B. Goyal, A. Dogra, S. Agrawal, B.S. Sohi, and A. Sharma.  Image denoising review: From classical to state-of-the-art approaches.  Information Fusion, 55:220–244, 2020.
  * [22] S. Gu, L. Zhang, W. Zuo, and X. Feng.  Weighted nuclear norm minimization with application to image denoising.  In 2014 IEEE Conference on Computer Vision and Pattern Recognition, pages 2862–2869, 2014.
  * [23] B. Han, Q. Mo, Z. Zhao, and X. Zhuang.  Directional compactly supported tensor product complex tight framelets with applications to image denoising and inpainting.  SIAM J. Imaging Sci., 12(4):1739–1771, 2019.
  * [24] B. Han and Z. Zhao.  Tensor product complex tight framelets with increasing directionality.  SIAM J. Imaging Sci., 7(2):997–1034, 2014.
  * [25] B. Han, Z. Zhao, and X. Zhuang.  Directional tensor product complex tight framelets with low redundancy.  Appl. Comput. Harmon. Anal., 41(2):603–637, 2016.
  * [26] A. E. Ilesanmi and T. O Ilesanmi.  Methods for image denoising using convolutional neural network: a review.  Complex Intell. Syst., 7:2179–2198, 2021.
  * [27] A. Jalobeanu, L. Blanc-Féraud, and J. Zerubia.  Satellite image deconvolution using complex wavelet packets.  In Proc. IEEE Int. Conf. Image Process. (ICIP), pages 809–812, 2000\.
  * [28] H. Ji, Z. Shen, and Y. Zhao.  Directional frames for image recovery: Multi-scale discrete Gabor frames.  Journal of Fourier Analysis and Applications, 23(4):729–757, 2017\.
  * [29] H. Ji, Z. Shen, and Y. Zhao.  Digital Gabor filters with MRA structure.  SIAM Journal on Multiscale Modeling and Simulation, 16(1):452–476, 2018.
  * [30] N.G. Kingsbury.  Image processing with complex wavelets.  Philos. Trans. R. Soc. London A, Math. Phys. Sci., 357(1760):2543–2560, 1999.
  * [31] G. Kutyniok and D. Labate.  Shearlets: Multiscale Analysis for Multivariate Data.  Birkhäuser, Boston, 2012.
  * [32] Xiaoxia Li, Juan Xiao, Yingyue Zhou, Yuanzheng Ye, Nianzu Lv, Xueyuan Wang, Shunli Wang, and ShaoBing Gao.  Detail retaining convolutional neural network for image denoising.  Journal of Visual Communication and Image Representation, 71:102774, 2020.
  * [33] J. Liu and S. Osher.  Block matching local SVD operator based sparsity and TV regularization for image denoising.  Journal of Scientific Computing, 78:607–624, 2019.
  * [34] Chong Mou, Jian Zhang, Xiaopeng Fan, Hangfan Liu, and Ronggang Wang.  Cola-net: Collaborative attention network for image restoration.  IEEE Transactions on Multimedia, 24:1366–1377, 2022.
  * [35] Yuhui Quan, Yixin Chen, Yizhen Shao, Huan Teng, Yong Xu, and Hui Ji.  Image denoising using complex-valued deep cnn.  Pattern Recognition, 111:107639, 2021.
  * [36] Y. Romano and M. Elad.  Boosting of image denoising algorithms.  SIAM J. IMAGING SCIENCES, 8(2):1187–1219, 2015.
  * [37] I.W. Selesnick, R.G. Baraniuk, and N.G. Kingsbury.  The dual-tree complex wavelet transform.  IEEE Signal Process. Mag., 22(6):123––151, 2005.
  * [38] L. Şendur and I. Selesnick.  Bivariate shrinkage with local variance estimation.  IEEE Signal Process. Letters, 9(12):438–441, 2002.
  * [39] Wuzhen Shi, Feng Jiang, Shengping Zhang, Rui Wang, Debin Zhao, and Huiyu Zhou.  Hierarchical residual learning for image denoising.  Signal Processing: Image Communication, 76:243–251, 2019.
  * [40] Pham Huu Thanh Binh, Cristóvão Cruz, and Karen Egiazarian.  Flashlight CNN image denoising.  In 2020 28th European Signal Processing Conference (EUSIPCO), pages 670–674, 2021.
  * [41] Chunwei Tian, Lunke Fei, Wenxian Zheng, Wangmeng Zuo Yong Xu and, and Chia-Wen Lin.  Deep learning on image denoising: An overview.  Neural Networks, 131, 2020.
  * [42] Z. Wang and A. C. Bovik.  Mean squared error : Love it or leave it?  IEEE Signal Process. Mag., 26:98–117, January 2009.
  * [43] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli.  Image quality assessment: From error visibility to structural similarity.  IEEE Trans. Image Proc., 13(4):600–612, 2004.
  * [44] Y.Hou and D. Shen.  Image denoising with morphology- and size-adaptive block-matching transform domain filtering.  EURASIP Journal on Image and Video Processing, 59, 2018.
  * [45] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang.  Beyond a Gaussian denoiser: residual learning of deep CNN for image denoising.  IEEE Trans. Image Process., 26(7):3142–3155, 2017.
  * [46] T. Zhou, C. Li, X. Zeng, and Y. Zhao.  Sparse representation with enhanced nonlocal self‑similarity for image denoising.  Machine Vision and Applications, 32(5):1–11, 2021.
  * [47] X. Zhuang.  Digital affine shear transforms: fast realization and applications in image/video processing.  SIAM J. Imag. Sci., 9(3):1437–1466, 2016.
  * [48] X. Zhuang and B. Han.  Compactly supported tensor product complex tight framelets with directionality.  In 2019 International Conference on Sampling Theory and Applications (SampTA), Bordeaux, France.