1

Accelerating GMM-based patch priors for image
restoration: Three ingredients for a 100× speed-up

arXiv:1710.08124v1 [] 23 Oct 2017

Shibin Parameswaran,Charles-Alban Deledalle, Loïc Denis and Truong Q. Nguyen

Abstract—Image restoration methods aim to recover the underlying clean image from corrupted observations. The Expected
Patch Log-likelihood (EPLL) algorithm is a powerful image
restoration method that uses a Gaussian mixture model (GMM)
prior on the patches of natural images. Although it is very
effective for restoring images, its high runtime complexity makes
EPLL ill-suited for most practical applications. In this paper, we
propose three approximations to the original EPLL algorithm.
The resulting algorithm, which we call the fast-EPLL (FEPLL),
attains a dramatic speed-up of two orders of magnitude over
EPLL while incurring a negligible drop in the restored image
quality (less than 0.5 dB). We demonstrate the efficacy and
versatility of our algorithm on a number of inverse problems
such as denoising, deblurring, super-resolution, inpainting and
devignetting. To the best of our knowledge, FEPLL is the first
algorithm that can competitively restore a 512×512 pixel image
in under 0.5s for all the degradations mentioned above without
specialized code optimizations such as CPU parallelization or
GPU implementation.
Index Terms—Image restoration, image patch, Gaussian mixture model, efficient algorithms

I. I NTRODUCTION
Patch-based methods form a very popular and successful
class of image restoration techniques. These methods process
an image on a patch-by-patch basis where a patch is a small
sub-image (e.g., of 8×8 pixels) that captures both geometric
and textural information. Patch-based algorithms have been at
the core of many state-of-the-art results obtained on various
image restoration problems such as denoising, deblurring, superresolution, defogging, or compression artifact removal to name
a few. In image denoising, patch-based processing became
popular after the success of the Non-Local Means algorithm [2].
Subsequently, continued research efforts have led to significant
algorithmic advancements in this area [1], [8], [39], [11],
[27], [37], [24]. Other inverse problems such as image superresolution and image deblurring have also benefited from patchbased models [9], [16], [30], [34], [14], [25].
Among these various patch-based methods, the Expected
Patch Log-Likelihood algorithm (EPLL) [39] deserves a special
mention due to its restoration performance and versatility.
The EPLL introduced an innovative application of Gaussian
S. Parameswaran, C.-A. Deledalle and T. Nguyen are with the Department
of Electrical and Computer Engineering, University of California, San
Diego, 9500 Gilman Dr, La Jolla, CA 92093, e-mail: sparames@ucsd.edu,
cdeledalle@eng.ucsd.edu and tqn001@eng.ucsd.edu.
C. Deledalle is also with IMB, CNRS, Univ. Bordeaux, Bordeaux INP,
F-33405 Talence, France, e-mail: charles-alban.deledalle@math.u-bordeaux.fr
L. Denis is with Univ Lyon, UJM-Saint-Etienne, CNRS, Institut d Optique
Graduate School, Laboratoire Hubert Curien UMR 5516, F-42023, SaintEtienne, France, e-mail: loic.denis@univ-st-etienne.fr.

Mixture Models (GMMs) to capture the prior distribution
of patches in natural images. Note that a similar idea was
introduced concurrently in [37]. The success of this method is
evident from the large number of recent works that extend the
original EPLL formulation [35], [26], [4], [29], [31], [18].
However, a persistent problem of EPLL-based algorithms
is their high runtime complexity. For instance, it is orders
of magnitude slower than the well-engineered BM3D image
denoising algorithm [8]. However, extensions of BM3D that
perform super-resolution [10] and other inverse problems [20]
require fundamental algorithmic changes, making BM3D far
less adaptable than EPLL. Other approaches that are as versatile
as EPLL [33], [5], [19] either lack the algorithmic efficiency
of BM3D or the restoration efficacy of EPLL.
Another class of techniques that arguably offers better runtime performance than EPLL-based methods (but not BM3D)
are those based on deep learning. With the advancements in
computational resources, researchers have attempted to solve
some classical inverse problems using multi-layer perceptrons
[3] and deep networks [6], [13], [21]. These methods achieve
very good restoration performance, but are heavily dependent
on the amount of training data available for each degradation
scenario. Most of these methods learn filters that are suited to
restore a specific noise level (denoising), blur (deblurring) or
upsampling factor (super-resolution), which makes them less
attractive to serve as generic image restoration solutions. More
recently, Zhang et al. [38] demonstrated the use of deep residual
networks for general denoising problems, single-image superresolution and compression artifact removal. Unlike earlier deep
learning efforts, their approach can restore images with different
noise levels using a single model which is learned by training on
image patches containing a range of degradations. Even in this
case, the underlying deep learning model requires retraining
whenever a new degradation scenario different from those
considered during the learning stage is encountered. Moreover,
it is much harder to gain insight into the actual model learned by
a deep architecture compared to a GMM. For this reason, even
with the advent of deep learning methods, flexible algorithms
like EPLL that have a transparent formulation remain relevant
for image restoration.
Recently, researchers have tried to improve the speed of
EPLL by replacing the most time-consuming operation in
the EPLL algorithm with a machine learning-based technique
of their choice [36], [32]. These methods were successful in
accelerating EPLL to an extent but did not consider tackling
all of its bottlenecks. In contrast, this paper focuses on
accelerating EPLL by proposing algorithmic approximations to
all the prospective bottlenecks present in the original algorithm

2

proposed by Zoran et al. [39]. To this end, we first provide a
complete computational and runtime analysis of EPLL, present
a new and efficient implementation of original EPLL algorithm
and then finally propose innovative approximations that lead to
a novel algorithm that is more than 100× faster compared to
the efficiently implemented EPLL (and 350× faster than the
runtime obtained by using the original implementation [39]).
Contributions: The main contributions of this work are
the following. We introduce three strategies to accelerate patchbased image restoration algorithms that use a GMM prior. We
show that, when used jointly, they lead to a speed-up of the
EPLL algorithm by two orders of magnitude. Compared to the
popular BM3D algorithm, which represents the current state-ofthe-art in terms of speed among CPU-based implementations,
the proposed algorithm is almost an order of magnitude faster.
The three strategies introduced in this work are general enough
to be applied individually or in any combination to accelerate
other related algorithms. For example, the random subsampling
strategy is a general technique that could be reused in any
algorithm that considers overlapping patches to process images;
the flat tail spectrum approximation can accelerate any method
that needs Gaussian log-likelihood or multiple Mahalanobis
metric calculations; finally, the binary search tree for Gaussian
matching can be included in any algorithm based on a GMM
prior model and can be easily adapted for vector quantization
techniques that use a dictionary.
For reproducibility purposes, we release our software on
GitHub along with a few usage demonstrations (available at
https://goo.gl/xjqKUA).
II. E XPECTED PATCH L OG -L IKELIHOOD (EPLL)
We consider the problem of estimating an image x ∈ RN
(N is the number of pixels) from noisy linear observations
y = Ax + w, where A : RN → RM is a linear operator
and w ∈ RM is a noise component assumed to be white and
Gaussian with variance σ 2 . In a standard denoising problem
A is the identity matrix, but in more general settings, it can
account for loss of information or blurring. Typical examples
for operator A are: a low pass filter (for deconvolution),
a masking operator (for inpainting), or a projection on a
random subspace (for compressive sensing). To reduce noise
and stabilize the inversion of A, some prior information is used
for the estimation of x. The EPLL introduced by Zoran and
Weiss [39] includes this prior information as a model for the
distribution of patches found in natural images. Specifically, the
EPLL defines the restored image as the maximum a posteriori
estimate, corresponding to the following minimization problem:
X
P
2
argmin
||Ax
−
y||
−
log p (Pi x)
(1)
2σ 2
x

Algorithm 1 The five steps of an EPLL iteration [39]
for all i ∈ I
z̃ i ← Pi x̂
ki?

← argmin
16ki 6K

(Patch extraction)
log wk−2
i

+ log

Σki + β1 IdP

−1

z̃ i
z̃ ti Σki + β1 IdP

−1
ẑ i ← Σki? + β1 IdP
Σki? z̃ i
X t −1 X t
Pi ẑ i
x̃ ←
Pi Pi
i∈I

+
(Gaussian selection)
(Patch estimation)
(Patch reprojection)

i∈I

x̂ ← At A + βσ 2 IdN

−1

At y + βσ 2 x̃



(Image estimation)

under the white Gaussian noise assumption), the second term
regularizes the solution x by favoring an image such that all its
patches fit the prior model of patches in natural images. The
authors of [39] showed that this prior can be well approximated
(upon removal of the DC component of each patch) using a
zero-mean Gaussian Mixture Model (GMM) with K = 200
components, that reads for any patch z ∈ RP , as


K
X
1
1 t −1
exp
−
z
Σ
z
, (2)
p (z) =
wk
k
2
(2π)P/2 |Σk |1/2
k=1
P
where the weights wk (such that wk > 0 and
k wk = 1)
P ×P
and the covariance matrices Σk ∈ R
are estimated using
the Expectation-Maximization algorithm [12] on a dataset
consisting of 2 million “clean” patches extracted from the
training set of the Berkeley Segmentation (BSDS) dataset [28].
Half-quadratic splitting: Problem (1) is a large nonconvex problem where A couples all unknown pixel values
x and the patch prior is highly non-convex. A classical
workaround, known as half-quadratic splitting [15], [23], is to
introduce N auxiliary unknown vectors z i ∈ RP , and consider
instead the penalized optimization problem that reads, for
β > 0, as
argmin
x,z 1 ,...,z N

P
||Ax − y||2
2σ 2
X
βX
||Pi x − z i ||2 −
log p (z i ) . (3)
+
2
i∈I

i∈I

When β → ∞, the problem (3) becomes equivalent to the
original problem (1). In practice, an increasing sequence of β
is considered, and an alternating optimization scheme is used:


β
ẑ i ← argmin ||Pi x̂ − z i ||2 − log p (z i )
(4)
2
zi
i=1..N
P
βX
x̂ ← argmin
||Ax − y||2 +
||Pi x − ẑ i ||2 . (5)
2
2σ
2
x
i∈I

i∈I

where I = {1, . . . , N } is the set of pixel indices, Pi : RN →
RP is the linear operator extracting a patch with P pixels
centered at the pixel with location i (typically, P = 8 × 8),
and p(.) is the a priori probability density function (i.e.,
the statistical model of noiseless patches in natural images).
While the first term in eq. (1) ensures that Ax is close to
the observations y (this term is the negative log-likelihood

Algorithm: Subproblem (5) corresponds to solving a linear
inverse problem with a Tikhonov regularization, and has an
explicit solution often referred to as Wiener filtering:
!−1
!
2 X
2 X
βσ
βσ
x̂ = At A +
Pit Pi
At y +
Pit ẑ i ,
P
P
i∈I

i∈I

(6)

3

TABLE I
C OMPARISON OF THE EXECUTION TIME OF OUR IMPLEMENTATION OF EPLL WITH AND WITHOUT PROPOSED ACCELERATIONS . E XPERIMENT CONDUCTED
ON A 481 × 321 IMAGE DENOISING PROBLEM . P ROFILING WAS CARRIED OUT USING MATLAB (R2014 B ) ON A PC WITH I NTEL (R) C ORE (TM) I 7-4790K
CPU @4.00GH Z AND 16 GB RAM. E XECUTION TIMES ARE REPORTED IN SECONDS ( S ) AND AS A PERCENTAGE OF THE TOTAL TIME (%).
Step
(Gaussian selection)
(Patch estimation)
(Patch extraction)
(Patch reprojection)
Others
Total

Without accelerations
43.53s
0.95s
2%
0.46s
1%
1%
0.23s
0.52s
1%
45.69s

95 %

where Pit Pi is a diagonal matrix whose i-th diagonal element
corresponds to the number of patches overlapping the pixel
of index i. This number is a constant equal to P (assuming
proper boundary conditions are used), which allows to split
the computation into two steps Patch reprojection and Image
estimation as shown in Alg. 1. Note that the step Patch
reprojection is simply the average of all overlapping patches.
In contrast, subproblem (4) cannot be obtained in closed
form as it involves a term with the logarithm of a sum of
exponentials. A practical solution proposed in [39] is to keep
only the component ki? maximizing the likelihood for the given
patch assuming it is a zero-mean Gaussian random vector
with covariance matrix Σki + β1 IdP . With this approximation,
the solution of (4) is also given by Wiener filtering, and the
resulting algorithm iterates the steps described in Alg. 1. The
authors of [39] found that using T = 5 iterations, with the
sequence β = σ12 {1, 4, 8, 16, 32}, for the initialization x̂ = y,
provides relevant solutions in denoising contexts for a wide
range of noise level values σ 2 .

With the proposed accelerations
0.23s
66 %
0.05s
13 %
0.03s
7%
0.01s
4%
0.03s
10 %
0.35s

number of operations per iteration of the alternate optimization
scheme. The steps Patch extraction and Patch reprojection
share a complexity of O(N P ). Finally, the complexity of
step Image estimation depends on the transform A. In many
scenarios of interest, At A can be diagonalized using a fast
transform and the inversion of At A+βσ 2 IdN can be performed
efficiently in the transformed domain (since IdN is diagonal
in any orthonormal basis). For instance, it leads to O(N )
operations for denoising or inpainting, and O(N log N ) for
periodical deconvolutions or super-resolution problems, thanks
to the fast Fourier transform. If A cannot be easily diagonalized,
this step can be performed using conjugate gradient (CG)
method, as done in [39], at a computational cost that depends
on the number of CG iterations (i.e., on the conditioning of
At A + βσ 2 IdN ). In any case, as shown in the next section,
this step has a complexity independent of P and K and is
one of the faster operations in the image restoration problems
considered in this paper.
B. Computation time analysis

A. Complexity via eigenspace implementation
The algorithm summarized in the previous section may
reveal cumbersome computations as it requires performing
numerous matrix multiplications and inversions. Nevertheless,
as the matrices Σk are known prior to any calculation, their
eigendecomposition can be computed offline to improve the
runtime. If we denote the eigendecomposition (obtained offline)
of Σk = U k S k U tk , such that U k ∈ RP ×P is unitary and S k is
diagonal with positive diagonal elements ordered in decreasing
order, steps Gaussian selection and Patch estimation can be
expressed in the space of coefficients c as


t
k
c̃i ← U k z̃ i k=1..K
O(N KP 2 ) (7)
i=1..N



ki? ← argmin ιk +
16k6K



k?

k?

[ĉi ]j ← γj i [c̃i i ]j


?
ẑ i ← U ki ĉi


P 
X
[c̃ki ]2j
log νjk + k
νj
i=1..N
j=1
O(N KP )

(8)

O(N P )

(9)

O(N P 2 )

(10)


j=1..P
i=1..N

i=1..N

where [c̃]j denotes the j-th entry of vector c̃, ιk = −2 log wk ,
νjk = [S k ]jj + β1 , and γjk = [S k ]jj /νjk with [S k ]jj the jth entry on the diagonal of matrix S k . The complexity of
each operation is indicated on its right and corresponds to the

In order to uncover the practical computational bottlenecks
of EPLL, we have performed the following computational
analysis. To identify clearly which part is time consuming, it is
important to make the algorithm implementation as optimal as
possible. Therefore, we refrain from using the MATLAB code
provided by the original authors [39] for speed comparisons.
Instead, we use a MATLAB/C version of EPLL based on the
eigenspace implementation described above, where some steps
are written in C language and interfaced using mex functions.
This version, which we refer to as EPLLc, provides results
identical to the original implementation while being 2-3 times
faster. The execution time of each step for a single run of
EPLLc is reported in the second column of Table I. Reported
times fit our complexity analysis and clearly indicate that the
Step Gaussian selection causes significant bottleneck due to
O(N P 2 K) complexity.
In the next section, we propose three independent modifications leading to an algorithm with a complexity of
O(N P r̄ log K/s2 ) with two constants 1 6 s2 6 P and
1 6 r̄ 6 P that control the accuracy of the approximations
introduced. The algorithm, in practice, is more than 100 times
faster as shown by its runtimes reported in the third column.
III. FAST EPLL: THE THREE KEY INGREDIENTS
We propose three accelerations based on (i) scanning only a
(random) subset of the N patches, (ii) reducing the number of

4

95%

Original spectrum
Flat tail approx.
0.95-Threshold rk

5%

Eigenvalue

0.4

0.1
10

20

30

40

50

60

Index of eigen dimension

Eigenvalue

10−1

95%

Fig. 2. (left) Our obtained search tree (with numbers of nodes for each level).
(right) Four patches well represented by each of the eight nodes along the
branch highlighted in red in the tree. These patches are randomly sampled
from the generative model encoded at each node of the branch (patches on
the same column are generated from the same random seed).

5%

10−3

where νPk = λk + β1 , γPk = λk /νPk . As ||z̃ i ||2 can be computed
once for all k, thePcomplexity of each step is divided by
10
20
30
40
50
60
K
1
Index of eigen dimension
P/r̄, where r̄ = K
k=1 rk is the average rank after which
eigenvalues are considered constant.
Fig. 1. Flat tail approximation: (a) with eigenvalues display on linear and (b)
In practice, covariance matrices Σk are not flat-tail but can
logarithmic scale. The rk = 15 first eigen components explains ρ = 95% of
be approximated by a flat-tail matrix by replacing the lowest
the variability.
eigenvalues by a constant λk . To obtain a small value of r̄
(hence a large speed-up), we preserve a fixed proportion ρ ∈
mixture components matched, and (iii) projecting on a smaller (0, 1] of the total variability and replace the smallest eigenvalues
subspace of the covariance eigenspace. We begin by describing accounting for the remaining 1−ρ fraction of the variability
this latter acceleration strategy in the following paragraph.
by their average (see Fig. 1): rk is the smallest integer such
that Tr(S̄ k ) > ρ Tr(S k ). Choosing ρ = 0.95 means that 5% of
the variability, in the eigendirections associated to the smallest
A. Speed-up via flat tail spectrum approximation
eigenvalues, is assumed to be evenly spread in these directions.
k
To avoid computing the P coefficients of the vector c̃i in
In practice, the choice of ρ = 0.95 leads to an average rank of
eq. (7), we rely on a flat-tail approximation. The k-th Gaussian
r̄ = 19.6 (for P = 8 × 8 ) for a small drop of PSNR as shown
model is said to have a flat tail if there exists a rank rk such
in Fig. 4. Among several other covariance approximations that
that for any j > rk , the eigenvalues are constant: [S k ]j,j = λk .
c
c
we tested, for instance, the one consisting in keeping only the
P ×rk
P ×rk
Denoting by Ū k ∈ R
(resp. Ū k ∈ R
) the matrix
rk first directions, the flat tail approximation provided the best
formed by the rk first (resp. rkc = P −rk last) columns of U k ,
trade-off in terms of acceleration and restoration quality.
c
c t
t
we have Ū k (Ū k ) = IdP − Ū k Ū k . It follows
t

(Σk + β1 IdP )−1 = Ū k (S̄ k + β1 Idrk )−1 Ū k
t

+ (λk + β1 )−1 (IdP − Ū k Ū k ),

(11)

t

(Σk + β1 IdP )−1 Σk = Ū k (S̄ k + β1 Idrk )−1 S̄ k Ū k
t

+ λk (λk + β1 )−1 (IdP − Ū k Ū k ), (12)
where S̄ k ∈ Rrk ×rk is the diagonal matrix formed by the rk
first rows and columns of S k . Steps Gaussian selection and
Patch estimation can thus be rewritten as


t
k
c̃i ← Ū k z̃ i k=1..K
O(N KP r̄) (13)
i=1..N

||z̃ i ||2
ki? ← argmin ιk + rkc log νPk + k
νP
16k6K


r
k
k
k
2
2
X
[c̃i ]j [c̃i ]j
+
log νjk + k − k
O(N K r̄) (14)
νj
νP
i=1..N
j=1


k?
k?  k? 
[ĉi ]j ← (γj i −γPi ) c̃i i j j=1..rki?
O(N P r̄) (15)
i=1..N


k?
ẑ i ← Ū ki? ĉi +γPi z̃ i
O(N P r̄) (16)
i=1..N

B. Speed-up via a balanced search tree
As shown in Table I, the step Gaussian selection has a
complexity of O(N P 2 K), reduced to O(N P K r̄) using the
flat tail spectrum approximation. This step remains the biggest
bottleneck since each query patch has to be compared to all
the K components of the GMM. To make this step even more
efficient, we reduce its complexity using a balanced search
tree. As described below, such a tree can be built offline by
repeatedly collapsing the original GMM to models with fewer
components, until the entire model is reduced to a single
Gaussian model.
We progressively combine the GMM components from one
level to form the level above, by clustering the K components
into L < K clusters of similar ones, until the entire model is
reduced to a single component. The similarity between two
zero-mean Gaussian models with covariance Σ1 and Σ2 is
measured by the symmetric Kullback-Leibler (KL) divergence
KL(Σ1 , Σ2 ) =

1
2

−1
Tr(Σ−1
2 Σ1 + Σ1 Σ2 − 2IdP ).

(17)

5

Based on this divergence, at each level n, we look for a partition
Ωn of the K Gaussian components into L clusters (with about
equal sizes) minimizing the following optimization problem
argmin
Ωn

such that

L
X

X

KL(Σk1 , Σk1 ),

(18)

l=1 k1 ,k2 ∈Ωn
l

(a) Regular patch subsampling at (i0 , j0 )

SL

n
n
n
n
l=1 Ωl = [K] and Ωl1 ∩ Ωl2 = ∅, where Ωl is
l1 6=l2

the l-th set of Gaussian components for the GMM at level n.
This clustering problem can be approximately solved using the
genetic algorithm of [22] for the Multiple Traveling Salesmen
Problem (MTSP). MTSP is a variation of the classical Traveling
Salesman Problem where several salesmen visit a unique set
(b) Stochastic patch subsampling at (i, j)
of cities and return to their origins, and each city is visited by
exactly one salesman. This attempts to minimize the total Fig. 3. Illustration of patch subsampling. Instead of extracting all patches,
distance traveled by all salesmen. Hence, it is similar to only a subset of patches is extracted either (a) regularly or (b) with some
our original problem given in eq. (18) where the Gaussian randomizations. Patches are represented by 8×8 squares and the red intensity
components and the clusters correspond to K cities and L represents the number of patches overlapping the corresponding pixel.
salesmen, respectively. Given the clustering at level n, the new
GMM at level n − 1 is obtained by combining the zero-mean the complexity of the four most time-consuming steps, see
Table I. One approach, followed by BM3D [8], consists of
Gaussian components such that, for all 1 6 l 6 L:
restricting
√ the set I to locations on a regular grid with spacing
X
1 X n n
n−1
n−1
n
s
∈
[1,
P ] pixels in both directions, leading to a reduction
wl
=
wk and Σl
= n−1
wk Σk , (19)
wl
of complexity by a factor s2 . We refer to this approach as
n
k∈Ω
k∈Ωn
l
l
the regular patch subsampling. A direct consequence is that
n
n
where Σk and wk are the corresponding covariance matrix and |I| = N/s2 and the complexity is divided by s2 . However, we
weight of the k-th Gaussian component at level n. Following observed that this strategy consistently creates blocky artifacts
this scheme, the original GMM of K = 200 components is revealing the regularity of the extraction pattern. A random
collapsed into increasingly more compact GMMs with K = 64, sampling approach, called "jittering", used in the computer
32, 16, 8, 4, 2 and 1 components. The main advantage of using graphics community [7] is preferable to limit this effect. This
MTSP compared to classical clustering approaches, is that this procedure ensures that each pixel is covered by at least one
procedure can be adapted easily to enforce approximately equal patch. The location (i , j ) of a point of the grid undergoes a
0 0
sized clusters, simply by enforcing that each salesman visits random perturbation, giving a new location (i, j) such that
at least 3 cities for the last level and 2 for the other ones.
j√
k
j√
k
P −s
P −s
We also experimented with other clustering strategies such as
i0 −
6
i
6
i
+
0
2
2
j√
k
j√
k
the hierarchical kmeans-like clustering in [17] and hierarchical
P −s
P −s
and
j
−
6
j
6
j
+
,
(20)
0
0
agglomerative clustering. With no principled way to enforce
2
2
even-sized clusters, these approaches, in general, lead to
where b·c denotes the flooring operation. We found experimenunbalanced trees (with comb structured branches) which
tally that independent and uniform perturbations offered the
result in large variations in computation times from one
best performance against all other tested strategies. In addition,
image to another. Although they all lead to similar denoising
we also resample these positions at each of the T iterations
performances, we opted for MTSP based clustering to build our
and add a (random) global shift to ensure that all pixels have
Gaussian tree in favor of obtaining a stable speed-up profile
the same expected number of patches covering them.
for our resulting algorithm.
Figure 3 illustrates the difference between a regular grid and
In Fig. 2 we show that the tree obtained using MTSP-based
a jittered grid of period s = 6 for patches of size P = 8×8. In
clustering is almost a binary tree (left) and also display the
both cases, all pixels are covered by at least one patch, but the
types of patches it encodes along a given path (right). Such a
stochastic version reveals an irregular pattern.
balanced tree structure lets one avoid testing each patch against
Nevertheless, when using random subsampling, a major
all K components. Instead, a patch is first compared to the
bottleneck occurs when At A is not diagonal because the
two first nodes in level 1 of the tree, then the branch providing
inversion involved in eq. (6) cannot be simplified as in Alg. 1.
the smallest cost is followed and the operation is repeated at
Using a conjugate gradient is a practical solution but will negate
higher levels until a leaf has been reached. Using this balanced
the reduction of complexity gained by using subsampling. To
search tree reduces the complexity of step Gaussian selection
the best of our knowledge, this is the main reason why patch
to O(N P r̄ log K).
subsampling has not been utilized to speed up EPLL. Here, we
follow a different path. We opt for approximating the solution
C. Speed-up via the restriction to a random subset of patches of the original problem (involving all patches) rather than
The simplest and most effective proposed acceleration evaluating the solution of an approximate problem (involving
consists of subsampling the set I of N patches to improve random subsample of patches). More precisely, we speed up

6

one of the steps of EPLL namely the Gaussian selection step.
Both use machine learning techniques to reduce its runtime.
In [36], the authors use a binary decision tree to approximate
the mapping z̃ i 7→ ki∗ performed in step Gaussian selection. At
each node k and level n of the tree, the patch z̃ i is confronted
with a linear separator in order to decide if the recursion should
continue on the left or right child given by

1

10

l3
Al

ub
+S
ee

T+
Fl

Tr

at

at
Fl

bs
Su

S

bS
Su

ee
T+

pl
am

Fl

O

Tr

in

g

ee
Tr

ta
il
at

rig

in
a

l

100

29.5

29

ub

3

ee

+S

Al
l

S

S
Tr

Fl

at

T+

Su
b

ee
T+
Tr
at

pl
am
bs

Su

Fl

in

g

ee
Tr

il
ta
at
Fl

O

rig

in

al

28.5

Fig. 4. (left) Average speed-up and (right) average PSNR for our three
accelerations, and all possible combinations of them, on the 40 images of the
BSDS validation set.

hank , z̃ i i + bnk > 0

(21)

where (ank , bnk ) are the parameters of the hyperplane for the k-th
node at level n. These separators are trained offline on all pairs
of (z̃ i , ki∗ ) obtained after the first iteration of EPLL for a given
β and noise level σ. Once a leaf has been reached, its index
provides a first estimate for the index ki∗ . To reduce errors due
to large variations among the neighboring pixels, this method
further employs a Markov random fields on the resulting map
of Gaussian components which runs in O(N K) complexity.
Hence, their overall approach reduces the complexity of step
Gaussian selection from O(N KP 2 ) to O(N (P D+K)), where
D = 12 is the depth of the learned decision tree.
In [32], the authors approximate the Gaussian selection step,
by using a gating (feed-forward) network with one hidden layer

Q 
X
z̃ i 7→ ιk +
log νjk +ωjk (d̃i )2j
with d̃i = V t z̃ i (22)

j=1
Alg. 1 by replacing the complete set of indices by the random
subset of patches. In this case, step Patch reprojection consists where Q is the size of the hidden layer. The matrix V ∈ RP ×Q
of averaging only this subset of overlapping restored patches. encodes the weights of the first layer, ω k corresponds to the
This novel and nuanced idea avoids additional overhead and weights of the hidden layer and they are learned discriminatively
attains dramatic complexity improvements compared to the to approximate the exact posterior probability:
standard approach.

P 
X
(c̃ki )2j
Experiments conducted on our validation dataset show that
k
k
z̃
→
7
ι
+
log
ν
+
with c̃ki = U tk z̃ i (23)
i
j
k
this strategy used with s = 6 leads to an acceleration of about
ν
j
j=1
36× with less than a 0.2dB drop in PSNR. In comparison, for
a similar drop of PSNR, the regular patch subsampling can that we encounter in eq. (7) and (8). Theoretically, a new
network will need to be trained for each type of degradations,
only achieve a 9× acceleration with s = 3.
noise levels and choices of β (recall that νjk = (S k )jj + β1 ).
However, the findings of [32] indicate that applying a network
D. Performance analysis
learned on clean patches and with β1 = 0 is effective regardless
Figure 4 shows the image restoration performance and speed- of the type of degradation or the value of β. Their main
up obtained when the three ingredients are applied separately advantage can be highlighted by comparing eq. (22) and (23)
or in combination. The results are averaged over 40 images where complexity is reduced from O(N KP 2 ) to O(N Q(K +
from the test set of BSDS dataset [28] that is set aside for P )). The authors utilize this benefit by choosing Q = 100.
Unlike these two approaches, our method does not try to learn
validation purposes. The speed-up is calculated with respect to
the EPLLc implementation which is labeled “original” in Fig. 4. the Gaussian selection rule directly (which depends on both
Among the three ingredients, random subsampling or jittering the noise level through 1/β and the prior model through the
(labeled “subsampling”) leads to the largest speed-up (32×), GMM). Instead, we simply define a hierarchical organization
while the usage of the search tree provides more than 7× faster of the covariance matrices Σk . In other words, while the two
processing. The average speed-up obtained when combining other approaches try to infer the posterior probabilities (or
all three ingredients is around 179× on our validation set directly the maximum a posteriori), our approach provides
consisting of images of size 481×321, for an average drop of an approximation to the prior model. During runtime, this
approximation of the prior is used in the posterior for the
PSNR less than 0.5dB.
Gaussian selection task. Please note that the value of β does
not play a role in determining the prior. This allows us to
IV. R ELATED METHODS
use the same search tree independently of the noise level,
To the best of our knowledge, there are only two other degradations, etc. Given that the main advantage of EPLL is
approaches [36], [32] that have attempted to accelerate EPLL. that the same model can be used for any type of degradations, it
Unlike our approach, these methods focus on accelerating only is important that this property remains true for the accelerated

7

TABLE II
PSNR, SSIM AND EXECUTION TIME ON THE BSDS TEST SET ( AVERAGE ON 60 IMAGES OF SIZE 481×321), AND ON SIX STANDARD IMAGES ( EACH OF
0
SIZE 512×512) FOR THE PROPOSED FEPLL AND FEPLL , EPLL ( WITH TIMING GIVEN FOR BOTH EPLL M [39] AND OUR EPLL C ), BM3D [8], CSF [33]
AND RO G [32] WITH 3 DIFFERENT LEVELS OF NOISE .
σ

5

20

60

Algo.
FEPLL
FEPLL0
EPLL
RoG
BM3D
CSF3×3
FEPLL
FEPLL0
EPLL
RoG
BM3D
CSF3×3
FEPLL
FEPLL0
EPLL
RoG
BM3D
CSF3×3
FEPLL
FEPLL0
EPLLc
EPLLm
RoG
BM3D
CSF3×3
gpu
CSF3×3

Berkeley
36.8
37.1
37.3
37.1
37.3
36.8
29.1
29.3
29.5
29.3
29.4
29.0
24.5
24.5
24.8
24.6
24.8
22.0

/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/

.959
.962
.963
.959
.962
.952
.812
.831
.836
.828
.824
.805
.614
.620
.631
.622
.637
.489

0.27
0.96
44.86
82.68
1.19
1.64
0.87
0.41

Barbara
36.9
37.4
37.6
36.7
38.3
37.0
29.1
29.7
29.9
28.4
31.8
28.3
23.5
23.8
24.0
23.3
26.3
21.4

/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/

.958
.961
.962
.954
.964
.955
.853
.869
.874
.838
.905
.822
.634
.646
.658
.624
.756
.487

0.42
1.91
78.66
145.15
1.99
2.79
1.09
0.41

Boat

Couple

Fingerprint

PSNR/SSIM
36.6 / .930
36.7 / .944
35.4 / .984
36.8 / .933
37.0 / .948
35.9 / .985
36.8 / .933
37.3 / .949
36.5 / .987
36.5 / .920
37.2 / .946
36.4 / .987
37.3 / .939
37.4 / .949
36.5 / .987
36.7 / .929
37.1 / .945
36.2 / .986
30.1 / .802
29.9 / .812
27.7 / .909
30.4 / .815
30.2 / .827
28.1 / .920
30.6 / .821
30.4 / .833
28.3 / .924
30.4 / .815
30.2 / .826
28.3 / .922
30.8 / .825
30.7 / .842
28.8 / .929
30.2 / .802
29.9 / .812
28.0 / .917
25.5 / .644
25.0 / .626
22.1 / .721
25.6 / .646
25.1 / .634
22.4 / .746
25.8 / .658
25.3 / .646
22.6 / .754
25.6 / .653
25.1 / .638
22.4 / .739
25.9 / .671
25.5 / .662
23.7 / .800
22.7 / .504
22.5 / .503
21.2 / .727
Time (in seconds)
0.40
0.41
0.38
1.90
1.90
1.90
78.71
78.78
78.17
143.68
144.30
144.13
1.97
1.96
1.94
2.91
2.83
2.45
1.09
1.13
1.09
0.41
0.41
0.41

version. Last but not least, the training of our search tree takes
a few minutes while the training steps for the above mentioned
approach take from several hours to a few days [36].
In the next section, we show that our proposed accelerations produce restoration results with comparable quality to
competing methods while requiring a smaller amount of time.

Lena
38.1
38.4
38.6
38.3
38.7
38.2
32.3
32.4
32.6
32.5
33.0
31.8
27.4
27.2
27.6
27.3
28.2
23.4

/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/

.941
.943
.944
.938
.944
.938
.863
.865
.869
.866
.877
.838
.742
.729
.746
.744
.775
.509

0.43
1.90
78.13
144.15
1.98
2.87
1.09
0.40

Mandrill
35.1
35.2
35.2
35.0
35.2
34.8
26.4
26.6
26.7
26.5
26.6
26.1
21.4
21.6
21.7
21.5
21.7
20.3

/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/
/

.959
.960
.960
.956
.959
.953
.786
.805
.808
.796
.796
.779
.468
.501
.506
.480
.501
.467

0.42
1.90
78.34
143.86
1.96
2.84
1.08
0.41

[36] and got confirmation that the implementation of their
algorithm (briefly described in IV) is not publicly available.
Due to certain missing technical details, we were unable to
reimplement it faithfully. However, the results reported in [36]
indicate that their algorithm performs in par with BM3D in
terms of both PSNR and time. Hence, BM3D results can be
used as a faithful proxy for the expected performance of Wang
et al.’s algorithm [36].
V. N UMERICAL EXPERIMENTS
To explicitly illustrate the quality vs. runtime tradeoff of
In this section, we present the results obtained on various
FEPLL,
we include results obtained using a slightly slower
image restoration tasks. Our experiments were conducted on
version
of
FEPLL referred to as FEPLL0 , that does not use the
standard images of size 512 × 512 such as Barbara, Boat,
Couple, Fingerprint, Lena, Mandrill and on 60 test images of balanced search tree and uses a flat tail spectrum approximation
0
size 481×321 from the Berkeley Segmentation Dataset (BSDS) with ρ = 0.98. Please note that FEPLL is not meant to be
[28] (the original BSDS test set contains 100 images, the other better or worse than FEPLL, it is just another version running
40 was used for validation purposes while setting parameters at a different PSNR/time tradeoff which allows us to compare
ρ and s). For denoising, we compare the performance of our our algorithm to others operating in different playing fields.
Finally, to illustrate the versatility of FEPLL, we also include
fast EPLL (FEPLL) to the original EPLL algorithm [39] and
BM3D [8]. For the original EPLL, we have included timing results for other inverse problems such as devignetting, superresults given by our own MATLAB/C implementation (EPLLc) resolution, and inpainting.
and the MATLAB implementation provided by the authors
Parameter settings: In our experiments, we use patches
(EPLLm). We also compare our restoration performance and of size P = 8 × 8, and the GMM provided by Zoran et
runtime against other fast restoration methods introduced to al. [39] with K = 200 components. The 200-components
achieve competitive trade-off between runtime efficiency and GMM is progressively collapsed into smaller GMMs with
image quality. These methods include RoG [32] (a method K = 64, 32, 16, 8, 4, 2 and 1, and then all Gaussians of the tree
accelerating EPLL based on feedforward networks described are modified offline by flat-tail approximations with ρ = 0.95.
in Sec. IV), and CSF [33] (a fast restoration technique using The final estimate for the restored image is obtained after 5
random field-based architecture).
iterations of our algorithm with β set to λσ −2 {1, 4, 8, 16, 32}
2
For deblurring experiments, we additionally compare with where λ = min{N −1 ||At A||2F ||A||−2
2 , 250σ }. For denoising,
field-of-experts (FoE)-based non-blind deconvolution [5] de- where A is identity, λ = 1 which boils down to the setting
noted as iPiano. We contacted the corresponding author of used by Zoran et al. [39]. For inverse problems, we found

8

30.6 / .872 (1.68s)

22.1 / .368

(a) Reference x

(b) Noisy image y

(c) BM3D result x̂

30.6 / .873 (44.88s)

(d) EPLLc result x̂

30.2 / .862 (0.36s)

(e) FEPLL result x̂

Fig. 5. Illustration of a denoising problem with noise standard deviation σ = 20. Part of: (a) the original image (b) its noisy version (c-e) denoised results of
competitive methods with PSNR/SSIM and time (inset).
29.6

29.4

29.2

29

28.8
10 -1

10 0

(a)

10 1

10 2

10 1

10 2

0.84

0.83

0.82

0.81

0.8

0.79
10 -1

10 0

(b)
Fig. 6. (a) PSNR and (b) SSIM versus time for different restoration methods
in a denoising problem with noise standard deviation σ = 20. PSNR, SSIM
and time are averaged on the 60 BSDS test images, each of size 481 × 321.
Optimal methods tend to be in the top-left corner.

that the initialization x̂ = (At A + 0.2σ 2 /λ∇)−1 At y, with ∇
the image Laplacian, provides relevant solutions whatever the
linear operator A and the noise level σ 2 . While the authors of
[39] do not provide any further direction for setting β and the
initialization in general inverse problems, our proposed setting
leads to competitive solutions irrespective of A and σ 2 . For
BM3D [8], EPLLm [39], RoG [32], CSF [33] and iPiano [5]
we use the implementations provided by the original authors
and use the default parameters prescribed by them.
Denoising: Table II shows the quantitative performances
of FEPLL on the denoising task compared to EPLLm [39],
EPLLc (our own MATLAB/C implementation), RoG [32]
BM3D [8] and CSF [33]. We evaluate the algorithms under

low-, mid- and high-noise settings by using Gaussian noise
of variance 52 , 202 and 602 , respectively. The result labeled
“Berkeley” is an average over 60 images from the BSDS testing
set [28]. Figures 6 provide graphical representations of these
performances in terms of PSNR/SSIM versus computation time
for the BSDS images for the noise variance setting σ 2 = 202 .
On average, FEPLL results are 0.5dB below regular EPLL
and BM3D; however, FEPLL is approximately 7 times faster
than BM3D, 170-200 times faster than EPLLc and over 350
times faster than EPLLm. FEPLL outperforms the faster CSF
algorithm in terms of both PSNR and time. In this case, FEPLL
is even faster than the GPU accelerated version of CSF (CSFgpu ).
Our approach is 4 times faster than RoG with a PSNR drop of
0.1-0.3dB. Nevertheless, if we slow down FEPLL to FEPLL0 ,
we can easily neutralize this quality deficit while still being
faster than RoG. Note that these accelerations are obtained
purely based on the approximations and no parallel processing
is used. Also, in most cases, a loss of 0.5dB may not affect
the visual quality of the image. To illustrate this, we show a
sample image denoised by BM3D, EPLL and FEPLL in Fig. 5.
Deblurring: Table III shows the performance of FEPLL
when used for deblurring as compared to RoG [32], iPiano
[5] and CSF [33]. For these experiments, we use the blur
kernel provided by Chen et al. [5] along with their algorithm
implementation. The results under the label “Berkeley” are
averaged over 60 images from the BSDS test dataset [28].
The results labeled “Classic” is averaged over the six standard
images (Barbara, Boat, Couple, Fingerprint, Lena and Mandrill).
FEPLL consistently outperforms its efficient competitors both
in terms of quality and runtime. Although the GPU version of
CSF is faster, the restoration quality obtained by CSF is 2-3dB
lower than FEPLL. The proposed algorithm outperforms RoG
by 1-1.8dB while running 3 and 5 times faster on “Berkeley”
and “Classic” datasets, respectively.
The superior qualitative performance of FEPLL is demonstrated in Fig. 7. For brevity, we only include the deblurring
results obtained from the top competitors of FEPLL algorithm
in terms of both quality and runtime. As observed, FEPLL
provides the best quality vs. runtime efficiency trade-off.
In contrast, a deblurring procedure using the regular EPLL
is around 350 times slower than FEPLL with the original

9

31.4 / .891 (0.17∗s)

24.9 / .624

(a) Reference x / Blur kernel

(b) Blurry image y

(c) CSFpw result x̂

32.2 / .910 (1.17s)

32.7 / .924 (0.46s)

(d) RoG result x̂

(e) FEPLL result x̂

Fig. 7. Illustration of a deblurring problem with noise standard deviation σ = 0.5. Part of: (a) the original image and the blur kernel (inset), (b) blurry version
(c-e) deblurred results of competitive methods with PSNR/SSIM and time (inset). The ’*’ indicates runtime on GPU while others are CPU times.
TABLE III
PSNR, SSIM AND EXECUTION TIME ON THE BSDS TEST SET ( AVERAGE OF
60 IMAGES OF SIZE 481×321), AND ON STANDARD IMAGES ( AVERAGE OF 6
IMAGES OF SIZE 512×512) FOR THE PROPOSED FEPLL AND FEPLL 0 , CSF
[33], RO G [32] AND C HEN ET AL .’ S [5] METHOD THAT IS CALLED iPiano IN
THEIR IMPLEMENTATION . T HE BLUR KERNEL USED IS THE ONE PROVIDED
IN ALONG WITH I P IANO IMPLEMENTATION AND NOISE IS SET TO σ = 0.5.
Algo.
iPiano
CSFpw
RoG
FEPLL
FEPLL’

Berkeley
PSNR/SSIM
Time (s)
29.5 / .824
29.53
30.2 / .875
0.50 (0.14*)
31.3 / .897
1.19
33.1 / .928
0.40
33.2 / .930
1.01

Classic
PSNR/SSIM
Time (s)
29.9 / .848
59.10
30.5 / 0.870
0.47 (0.14*)
31.8 / .915
2.07
32.8 / .931
0.46
33.0 / .933
1.82

implementation [39]. Specifically, on the sample image shown
in Fig. 7, EPLL provides a qualitatively similar result (not
shown in the figure) with a PSNR of 32.7 dB and SSIM of
0.922 in 142 seconds.
Other inverse problems: Unlike BM3D, EPLL and FEPLL
are more versatile and handle a wide range of inverse problems
without any change in formulation. In Fig. 8, we show the
results obtained by FEPLL on problems such as (a) devignetting,
which involves a progressive loss of intensity, (b) superresolution and (c) inpainting. To show the robustness of our
method, the input images of size 481×321 were degraded with
zero-mean Gaussian noise with σ = 2. All of the restoration
results were obtained within/under 0.4 seconds and with the
same set of parameters explained above (cf. Parameter settings).
VI. C ONCLUSION
In this paper, we accelerate EPLL by a factor greater than
100 with negligible loss of image quality (less than 0.5dB).
This is achieved by combining three independent strategies: a
flat tail approximation, matching via a balanced search tree,
and stochastic patch sampling. We show that the proposed
accelerations are effective in denoising and deblurring problems,
as well as in other inverse problems such as super-resolution
and devignetting. An important distinction of the proposed
accelerations is their genericity: the accelerated EPLL prior
can be applied to many restoration tasks and various signalto-noise ratios, in contrast to existing accelerations based on
learning techniques applied to specific conditions (such as

image size, noise level, blur kernel, etc.) and that require an
expensive re-training to address a different problem.
Since the speed-up is achieved solely by reducing the
algorithmic complexity, we believe that further inclusion of
accelerations based on parallelization and/or GPU implementations will allow for real-time video processing. Moreover,
the acceleration techniques introduced in this work are general
strategies that can be used to speed up other image restoration
and/or related machine learning algorithms. For reproducibility
purposes, the code of FEPLL is made available on GitHub1 .
R EFERENCES
[1] M. Aharon, M. Elad, and A. Bruckstein. k -svd: An algorithm for
designing overcomplete dictionaries for sparse representation. IEEE
Transactions on Signal Processing, 54(11):4311–4322, Nov 2006.
[2] A. Buades, B. Coll, and J.-M. Morel. A review of image denoising
algorithms, with a new one. Multiscale Modeling & Simulation, 4(2):490–
530, 2005.
[3] H. Burger, C. Schuler, and S. Harmeling. Image denoising: Can plain
neural networks compete with BM3D? In International Conference on
Computer Vision and Pattern Recognition (CVPR), 2012.
[4] N. Cai, Y. Zhou, S. Wang, B. W.-K. Ling, and S. Weng. Image denoising
via patch-based adaptive gaussian mixture prior method. Signal, Image
and Video Processing, 10(6):993–999, 2016.
[5] Y. Chen, R. Ranftl, and T. Pock. Insights into analysis operator learning:
From patch-based sparse models to higher order mrfs. IEEE Transactions
on Image Processing, 23(3):1060–1072, March 2014.
[6] Y. Chen, W. Yu, and T. Pock. On learning optimized reaction diffusion
processes for effective image restoration. In Proceedings of the IEEE
Conference on Computer Vision and Pattern Recognition, pages 5261–
5269, 2015.
[7] R. L. Cook. Stochastic sampling in computer graphics. ACM Transactions
on Graphics (TOG), 5(1):51–72, 1986.
[8] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian. Image denoising by
sparse 3-D transform-domain collaborative filtering. IEEE Transactions
on Image Processing, 16(8):2080–2095, August 2007.
[9] A. Danielyan, A. Foi, V. Katkovnik, and K. Egiazarian. Image upsampling
via spatially adaptive block-matching filtering. In Signal Processing
Conference, 2008 16th European, pages 1–5. IEEE, 2008.
[10] A. Danielyan, A. Foi, V. Katkovnik, and K. Egiazarian. Spatially adaptive
filtering as regularization in inverse imaging: Compressive sensing superresolution and upsampling. Super-Resolution Imaging, pages 123–154,
2010.
[11] C.-A. Deledalle, J. Salmon, and A. S. Dalalyan. Image denoising with
patch based PCA: local versus global. In BMVC 2011-22nd British
Machine Vision Conference, pages 25–1. BMVA Press, 2011.
[12] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from
incomplete data via the EM algorithm. Journal of the royal statistical
society. Series B (methodological), pages 1–38, 1977.
1 https://goo.gl/xjqKUA

10

11.1 / .662

20.8 / .598

8.31 / .112

36.8 / .972 (0.38s)

23.3 / .738 (0.29s)

27.0 / .905 (0.36s)

(a) devignetting

(b) ×3 super-resolution

(c) 50% inpainting

Fig. 8. FEPLL on various inverse problems. All inputs contain Gaussian noise with σ = 2. Top row: (a) the observation in a devignetting problem, (b) the
bi-cubic interpolation and the actual low-resolution size image (inset) in a ×3 super-resolution problem and (c) the observation in an inpainting problem with
50% of missing pixels shown in red. Bottom row: respective FEPLL results all obtained in less than 0.4s.
[13] C. Dong, C. C. Loy, K. He, and X. Tang. Image super-resolution using
deep convolutional networks. IEEE transactions on pattern analysis and
machine intelligence, 38(2):295–307, 2016.
[14] K. Egiazarian and V. Katkovnik. Single image super-resolution via
BM3D sparse coding. In Signal Processing Conference (EUSIPCO),
2015 23rd European, pages 2849–2853. IEEE, 2015.
[15] D. Geman and C. Yang. Nonlinear image recovery with half-quadratic
regularization. IEEE Transactions on Image Processing, 4(7):932–946,
1995.
[16] D. Glasner, S. Bagon, and M. Irani. Super-resolution from a single image.
In IEEE 12th International Conference on Computer Vision, 2009.
[17] J. Goldberger and S. T. Roweis. Hierarchical clustering of a mixture
model. In NIPS, volume 17, pages 505–512, 2004.
[18] A. Houdard, C. Bouveyron, and J. Delon. High-Dimensional Mixture
Models For Unsupervised Image Denoising (HDMI). Preprint hal01544249, Aug. 2017.
[19] J. Jancsary, S. Nowozin, and C. Rother. Regression tree fields – an
efficient, non-parametric approach to image labeling problems. In 25th
IEEE Conference on Computer Vision and Pattern Recognition (CVPR),
2012.
[20] V. Katkovnik and K. Egiazarian. Nonlocal image deblurring: Variational
formulation with nonlocal collaborative L0-norm prior. In Local and
Non-Local Approximation in Image Processing, 2009. LNLA 2009.
International Workshop on, pages 46–53. IEEE, 2009.
[21] J. Kim, J. Kwon Lee, and K. Mu Lee. Accurate image super-resolution
using very deep convolutional networks. In Proceedings of the IEEE
Conference on Computer Vision and Pattern Recognition, pages 1646–
1654, 2016.
[22] J. Kirk. Source code for "multiple traveling salesmen problem genetic algorithm", https://fr.mathworks.com/matlabcentral/fileexchange/
19049-multiple-traveling-salesmen-problem-genetic-algorithm. MATLAB Central File Exchange. Retrieved May 6, 2014, 2008.
[23] D. Krishnan and R. Fergus. Fast image deconvolution using hyperlaplacian priors. In Advances in Neural Information Processing Systems,
pages 1033–1041, 2009.
[24] M. Lebrun, A. Buades, and J.-M. Morel. A nonlocal bayesian image
denoising algorithm. SIAM Journal on Imaging Sciences, 6(3):1665–1688,
2013.
[25] E. López-Rubio. Superresolution from a single noisy image by the median
filter transform. SIAM Journal on Imaging Sciences, 9(1):82–115, 2016.
[26] E. Luo, S. H. Chan, and T. Q. Nguyen. Adaptive image denoising by

[27]
[28]

[29]
[30]
[31]
[32]

[33]
[34]
[35]
[36]
[37]
[38]
[39]

mixture adaptation. IEEE Transactions on Image Processing, 25(10):4489–
4503, 2016.
J. Mairal, F. Bach, and J. Ponce. Task-driven dictionary learning. IEEE
Transactions on Pattern Analysis and Machine Intelligence, 34(4):791–
804, 2012.
D. Martin, C. Fowlkes, D. Tal, and J. Malik. A database of human
segmented natural images and its application to evaluating segmentation
algorithms and measuring ecological statistics. In Computer Vision, 2001.
ICCV 2001. Proceedings. Eighth IEEE International Conference on,
volume 2, pages 416–423. IEEE, 2001.
V. Papyan and M. Elad. Multi-scale patch-based image restoration. IEEE
Transactions on image processing, 25(1):249–261, 2016.
M. Protter, M. Elad, H. Takeda, and P. Milanfar. Generalizing the
nonlocal-means to super-resolution reconstruction. IEEE Transactions
on image processing, 18(1):36–51, 2009.
Y. Ren, Y. Romano, and M. Elad. Example-based image synthesis via
randomized patch-matching. arXiv preprint arXiv:1609.07370, 2016.
D. Rosenbaum and Y. Weiss. The return of the gating network: combining
generative models and discriminative training in natural image priors. In
Advances in Neural Information Processing Systems, pages 2683–2691,
2015.
U. Schmidt and S. Roth. Shrinkage fields for effective image restoration.
In Proceedings of the IEEE Conference on Computer Vision and Pattern
Recognition, pages 2774–2781, 2014.
A. Singh, F. Porikli, and N. Ahuja. Super-resolving noisy images. In
Proceedings of the IEEE Conference on Computer Vision and Pattern
Recognition, pages 2846–2853, 2014.
J. Sulam and M. Elad. Expected patch log likelihood with a sparse prior.
In International Workshop on Energy Minimization Methods in Computer
Vision and Pattern Recognition, pages 99–111. Springer, 2015.
Y. Wang, S. Cho, J. Wang, and S.-F. Chang. Discriminative indexing for
probabilistic image patch priors. In European Conference on Computer
Vision, pages 200–214. Springer, 2014.
G. Yu, G. Sapiro, and S. Mallat. Solving inverse problems with piecewise
linear estimators: From gaussian mixture models to structured sparsity.
IEEE Transactions on Image Processing, 21(5):2481–2499, 2012.
K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang. Beyond a Gaussian
Denoiser: Residual Learning of Deep CNN for Image Denoising. arXiv
preprint arXiv:1608.03981, 2016.
D. Zoran and Y. Weiss. From learning models of natural image patches
to whole image restoration. In International Conference on Computer
Vision, pages 479–486. IEEE, November 2011.

