A Bayesian Hyperprior Approach for Joint Image Denoising and Interpolation, with an Application to HDR Imaging

TL;DR

This paper introduces a Bayesian hyperprior method for joint image denoising and interpolation, enhancing stability and adaptability for inverse problems like HDR imaging from single, sensor-modified images.

Contribution

It proposes a hyperprior-based Bayesian approach that stabilizes patch modeling, especially for diagonal degradation and signal-dependent noise, improving HDR imaging applications.

Findings

Effective inpainting and zooming for missing data

Handles signal-dependent noise in digital cameras

Demonstrates improved HDR imaging results

Abstract

Recently, impressive denoising results have been achieved by Bayesian approaches which assume Gaussian models for the image patches. This improvement in performance can be attributed to the use of per-patch models. Unfortunately such an approach is particularly unstable for most inverse problems beyond denoising. In this work, we propose the use of a hyperprior to model image patches, in order to stabilize the estimation procedure. There are two main advantages to the proposed restoration scheme: Firstly it is adapted to diagonal degradation matrices, and in particular to missing data problems (e.g. inpainting of missing pixels or zooming). Secondly it can deal with signal dependent noise models, particularly suited to digital cameras. As such, the scheme is especially adapted to computational photography. In order to illustrate this point, we provide an application to high dynamic range imaging from a single image taken with a modified sensor, which shows the effectiveness of the proposed scheme.

Tables

Table 1. TABLE I: Results of the interpolation, combined interpolation and denoising, denoising and zooming tests described in Section  IV-A . Patch size of 8 × 8 8 8 8\times 8 for all methods in all tests. Parameter setting for interpolation, combined interpolation and denoising, and zooming, HBE: α H = 1 subscript 𝛼 𝐻 1 \alpha_{H}=1 , α L = 0.5 subscript 𝛼 𝐿 0.5 \alpha_{L}=0.5 , PLE: σ = 3 𝜎 3 \sigma=3 , ε = 30 𝜀 30 \varepsilon=30 , K = 19 𝐾 19 K=19   [ 10 ] , EPLL: default parameters  [ 24 ] , E-PLE: parameters set as specified in  [ 25 ] . Parameter setting for denoising, HBE: α H = α L = 100 subscript 𝛼 𝐻 subscript 𝛼 𝐿 100 \alpha_{H}=\alpha_{L}=100 , NLB: code provided by the authors  [ 26 ] automatically sets parameters from input σ 2 superscript 𝜎 2 \sigma^{2} , EPLL: default parameters for the denoising example  [ 24 ]

	Interpolation - PSNR (dB)												Interpolation & Denoising - PSNR (dB)
	50%				70%				90%				70%				90%
	HBE	PLE	EPLL	E-PLE	HBE	PLE	EPLL	E-PLE	HBE	PLE	EPLL	E-PLE	HBE	PLE	EPLL	E-PLE	HBE	PLE	EPLL	E-PLE
barbara	39.11	36.93	32.99	35.43	34.69	32.50	27.96	28.77	24.86	23.62	23.30	23.26	33.34	31.99	27.63	27.75	24.57	23.53	23.27	23.20
boat	34.92	34.32	34.21	33.59	31.37	30.74	30.38	30.26	25.96	25.35	24.72	25.43	30.61	30.41	30.15	29.54	25.78	25.45	24.71	25.47
traffic	30.17	30.12	30.19	28.86	27.27	27.12	27.13	26.64	22.84	22.34	21.85	22.27	26.99	26.98	27.05	26.35	22.87	22.43	22.21	22.35
	Denoising - PSNR (dB)													Zooming - PSNR (dB)
${\sigma }^2$	10			30			50			80				$\times 2$
	HBE	NLB	EPLL	HBE	NLB	EPLL	HBE	NLB	EPLL	HBE	NLB	EPLL		HBE	PLE	EPLL	E-PLE	Lanczos
barbara	41.26	41.20	40.56	38.40	38.26	37.32	37.13	36.94	35.84	35.96	35.73	34.51		38.17	37.11	31.34	36.51	28.01
boat	40.05	39.99	39.47	36.71	36.76	36.34	35.41	35.46	35.13	34.30	34.33	34.12		32.35	31.96	31.95	32.08	29.60
traffic	40.73	40.74	40.55	37.03	36.99	36.86	35.32	35.26	35.20	33.78	33.70	33.72		25.05	24.78	25.17	24.91	21.89

Table 2. TABLE II: PSNR values for the extracts shown in Figure  7 .

	PSNR (dB)
	1 (Fig. 7)	2 (Fig. 7)	3	4	5	6
HBE	33.08	33.87	22.95	35.10	36.80	35.66
PLEV	29.65	30.82	22.77	33.99	36.42	34.73
Schöberl et al.	30.38	31.16	21.39	30.04	32.84	31.02
Nayar and Mitsunaga	29.39	30.10	23.24	25.83	30.26	26.90

Equations

$$\mathbf{Z}=\mathbf{D}\mathbf{C}+\mathbf{N},$$

$${\mathbf{Z}}_i={\mathbf{D}}_i{\mathbf{C}}_i+{\mathbf{N}}_i,$$

$$p({\mathbf{Z}}_i|{\mathbf{C}}_i)\sim \mathcal{N}({\mathbf{D}}_i{\mathbf{C}}_i,{\mathbf{\Sigma }}_{N_i})$$

$$\propto |{\mathbf{\Sigma }}_{N_i}^{-1}{|}^{\frac{1}{2}}\exp \left(-\frac{1}{2}({\mathbf{Z}}_i-{\mathbf{D}}_i{\mathbf{C}}_i{)}^T{\mathbf{\Sigma }}_{N_i}^{-1}({\mathbf{Z}}_i-{\mathbf{D}}_i{\mathbf{C}}_i)\right).$$

$$p(\mathbf{\mu },\mathbf{\Lambda })$$

$$\propto |\mathbf{\Lambda }{|}^{1/2}\exp \left(-\frac{\kappa }{2}(\mathbf{\mu }-{\mathbf{\mu }}_0)\mathbf{\Lambda }(\mathbf{\mu }-{\mathbf{\mu }}_0{)}^T\right)$$

$$|\mathbf{\Lambda }{|}^{(\nu -n-1)/2}\exp \left(-\frac{1}{2}\text{tr}(\nu {\mathbf{\Sigma }}_0\mathbf{\Lambda })\right),$$

$$\underset{\{{\mathbf{C}}_i\},\mathbf{\mu },\mathbf{\Lambda }}{\mathrm{argmax}}p(\{{\mathbf{C}}_i\},\mathbf{\mu },\mathbf{\Lambda }|\{{\mathbf{Z}}_i\})=$$

$$=$$

$$f(\{{\mathbf{C}}_i\},\mathbf{\mu },\mathbf{\Lambda })$$

$$\hat{\mathbf{\mu }}$$

$$\widehat{{\mathbf{C}}_i}$$

$${\hat{\mathbf{\Lambda }}}^{-1}=\frac{\nu {\mathbf{\Sigma }}_0+\kappa (\mathbf{\mu }-{\mathbf{\mu }}_0)(\mathbf{\mu }-{\mathbf{\mu }}_0{)}^T+{\sum }_{i=1}^M({\mathbf{C}}_i-\mathbf{\mu })({\mathbf{C}}_i-\mathbf{\mu }{)}^T}{\nu +M-n}.$$

$$d(p,q)=\frac{{\sum }_{j=1}^N({\mathbf{c}}_p^j-{\mathbf{c}}_q^j{)}^2{\omega }_{p,q}^j}{{\sum }_{j=1}^N{\omega }_{p,q}^j},$$

$$\kappa=M\alpha,\quad\alpha=\left\{\begin{array}[]{rl}\alpha_{L}&\text{if $P$ and $M$}>\text{threshold}\\ \alpha_{H}&\text{otherwise}.\end{array}\right.$$

$$\nu =M\alpha +n$$

$${\mathbf{\mu }}_0=\frac{1}{M}\sum _{j=1}^M{\tilde{\mathbf{c}}}_j,{\mathbf{\Sigma }}_0=\frac{1}{M-1}\sum _{j=1}^M({\tilde{\mathbf{c}}}_j-{\mathbf{\mu }}_0)({\tilde{\mathbf{c}}}_j-{\mathbf{\mu }}_0{)}^T.$$

$$\mathrm{Z}(p)\sim \mathcal{N}(g{a}_p\tau \mathrm{C}(p)+{\mathbf{\mu }}_R,g^2{a}_p\tau \mathrm{C}(p)+{\sigma }_R^2),$$

$$\mathrm{Z}(p)\sim \mathcal{N}(g{o}_p{a}_p\tau \mathrm{C}(p)+{\mathbf{\mu }}_R,g^2{o}_p{a}_p\tau \mathrm{C}(p)+{\sigma }_R^2).$$

$$\mathrm{Y}(p)=\frac{\mathrm{Z}(p)-{\mathbf{\mu }}_R}{g{o}_p{a}_p\tau }.$$

$$(\boldsymbol{D})_{p}=\left\{\begin{array}[]{rl}1&\text{if }\boldsymbol{\mu}_{R}<\mathrm{Z}(p)<z_{sat},\\ 0&\text{otherwise},\end{array}\right.$$

$$\mathrm{Y}(p)|(\mathbf{D}{)}_p\sim \mathcal{N}\left((\mathbf{D}{)}_p\mathrm{C}(p),\frac{g^2{o}_p{a}_p\tau (\mathbf{D}{)}_p\mathrm{C}(p)+{\sigma }_R^2}{(g{o}_p{a}_p\tau {)}^2}\right).$$

$$\mathrm{Y}=\mathbf{D}\mathbf{C}+\mathbf{N},$$

$$({\mathbf{\Sigma }}_{\mathbf{N}}{)}_j=\frac{g^2{o}_p{a}_p\tau (\mathbf{D}{)}_p\mathrm{C}(p)+{\sigma }_R^2}{(g{o}_p{a}_p\tau {)}^2}.$$

$$\underset{\Vert (\{{\mathbf{C}}_i\},\mathbf{\mu },\mathbf{\Lambda })\Vert \to +\infty }{\lim }f(\{{\mathbf{C}}_i\},\mathbf{\mu },\mathbf{\Lambda })$$

$$f(\{{\mathbf{C}}_i\},\mathbf{\mu },\mathbf{\Lambda })$$

$$+\frac{\kappa }{2}(\mathbf{\mu }-{\mathbf{\mu }}_0{)}^T\mathbf{\Lambda }(\mathbf{\mu }-{\mathbf{\mu }}_0)$$

$$f(\{{\mathbf{C}}^{l_k}\},{\mathbf{\mu }}^{l_k},{\mathbf{\Lambda }}^{l_k})\le f(\{{\mathbf{C}}^{l_k}\},{\mathbf{\mu }}^{l_k},\mathbf{\Lambda }),\forall \mathbf{\Lambda }\in S_n^{++}(\mathbb{R})$$

$$f(\{{\mathbf{C}}^{\ast }\},{\mathbf{\mu }}^{\ast },{\mathbf{\Lambda }}^{\ast })=\underset{\mathbf{\Lambda }\in S_n^{++}(\mathbb{R})}{\mathrm{argmin}}f(\{{\mathbf{C}}^{\ast }\},{\mathbf{\mu }}^{\ast },\mathbf{\Lambda }).$$

$$(\{{\mathbf{C}}^{\prime }\},{\mathbf{\mu }}^{\prime })$$

$${\mathbf{\Lambda }}^{\prime }$$

$$f(G(\{{\mathbf{C}}^{l_{k-1}}\},{\mathbf{\mu }}^{l_{k-1}},{\mathbf{\Lambda }}^{l_{k-1}}))$$

Full text

I Introduction

Digital images are subject to a wide variety of degradations, which in most cases can be modeled as

[TABLE]

where $\mathbf{Z}$ is the observation, $\boldsymbol{D}$ is the degradation operator, $\mathbf{C}$ is the underlying ground-truth image and $\mathbf{N}$ is additive noise. Different settings of the degradation matrix $\boldsymbol{D}$ model different problems such as zooming, deblurring or missing pixels. Different versions of the noise term $\mathbf{N}$ include the classical additive Gaussian noise with constant variance or more complicated and realistic models such as signal dependent noise.

Due to the inherent ill-posedness of such inverse problems, standard approaches impose some prior on the image, in either variational or Bayesian approaches. Popular image models have been proposed through the total variation [1], wavelet decompositions [2] or the sparsity of image patches [3]. Buades et al. [4] introduced the use of patches and the self-similarity hypothesis to the denoising problem leading to a new era of patch-based image restoration techniques. A major step forward in fully exploiting the potential of patches was achieved by several state-of-the-art restoration methods with the introduction of patch prior models, in a Bayesian framework. Some methods are devoted to the denoising problem [5, 6, 7, 8], while others propose a more general framework for the solution of image inverse problems [9, 10], including inpainting, deblurring and zooming. The work by Lebrun et al. [11, 7] presents a thorough analysis of several recent restoration methods, revealing their common roots and their relationship with the Bayesian approach.

Among the state-of-the-art restoration methods, two noticeable approaches are the patch-based Bayesian approach by Yu et al. [10], namely the piece-wise linear estimators (PLE), and the non-local Bayes (NLB) algorithm by Lebrun et al. [7]. PLE is a general framework for the solution of image inverse problems under Model (1), while NLB is a denoising method ($\boldsymbol{D}=Id$). Both methods use a Gaussian patch prior learnt from image patches through iterative procedures. In the case of PLE, patches are modeled according to a Gaussian Mixture Model (GMM), with a relatively small number of classes (19 in all their experiments), whose parameters are learnt from all image patches111Actually, the authors report the use of $128\times 128$ image sub-regions in their experiments, so we may consider PLE as a semi-local approach.. In the case of NLB, each patch is associated with a single Gaussian model, whose parameters are computed from similar patches chosen from a local neighbourhood, i.e., a search window centered at the patch. We refer hereafter to this kind of per-patch modelling as local.

Zoran and Weiss [9] (EPLL) follow a similar approach, but instead of iteratively updating the GMM from image patches, they use a larger number of classes that are fixed and learnt from a large database of natural image patches. Wang and Morel [8] claim that, in the case of denoising, it is better to have fewer models that are updated with the image patches (as in PLE) than having a large number of fixed models (as in EPLL).

All of the previous restoration approaches share a common Bayesian framework based on Gaussian patch priors. Relying on local priors [7, 8] has proven more accurate for the task of image denoising than relying on a mixture of a limited number of Gaussian models [9, 10]. In particular, NLB outperforms PLE for this task [12], mostly due to its local model estimation. On the other hand, PLE yields state-of-the-art results in other applications such as interpolation of missing pixels (especially with high masking rates), deblurring and zooming.

As a consequence we are interested in taking advantage of a local patch modelling for more general inverse problems than denoising. The main difficulty lies in the estimation of the models, especially when the image degradations involve a high rate of missing pixels, in which case the estimation is seriously ill-posed.

In this work we propose to model image patches according to a Gaussian prior, whose parameters, the mean $\boldsymbol{\mu}$ and the covariance matrix $\boldsymbol{\Sigma}$, will be estimated locally from similar patches. In order to tackle this problem, we include prior knowledge on the model parameters making use of a hyperprior, i.e. a probability distribution on the parameters of the prior. In Bayesian statistics, $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ are known as hyperparameters, while the prior on them is called a hyperprior. Such a framework is often called hierarchical Bayesian modelling [13]. Its application to inverse problems in imaging is not new. In particular, in the field of image restoration, this methodology was proposed by Molina et al. [14, 15], and was more recently applied to image unmixing problems [16] and to image deconvolution and the estimation of the point spread function of a camera [17]. However, to our knowledge, this is the first time that such a hierarchical Bayesian methodology is used to reduce ill-posedness in patch-based image restoration. In this context, the use of a hyperprior compensates for the patches missing information.

There are two main contributions of this work: First, as described above, we propose a robust framework enabling the use of Gaussian local priors on image patches for solving a useful family of restoration problems by drawing on a hierarchical Bayesian approach. The second advantage of the proposed framework is its ability to deal with signal dependent noise, therefore making it adapted to realistic digital photography applications.

Experiments on both synthetic and real data show that the approach is well suited to various problems involving a diagonal degradation operator. First, we show state-of-the-art results in image restoration problems such as denoising, zooming and interpolation of missing pixels. Then we consider the generation of high dynamic range (HDR) images from a single snapshot using spatially varying pixel exposures [18] and demonstrate that our approach significantly outperforms existing methods to deal with this inverse problem. It is worth mentioning that modified sensors enabling such approaches have been recently made available by Sony but are not yet fully exploited by available smartphones and digital cameras.

The article is organized as follows. Section II introduces the proposed approach while Section III presents the main implementation aspects. Supportive experiments are presented in Section IV. Section V is devoted to the application of the proposed framework to the HDR imaging problem. Last, conclusions are summarized in Section VI.

II Hyperprior Bayesian Estimator

The proposed restoration method, called Hyperprior Bayesian Estimator (HBE), assumes a Gaussian prior for image patches, with parameters $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$. A joint maximum a posteriori formulation is used to estimate both the image patches and the parameters $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$, thanks to a Bayesian hyperprior model on these parameters, stabilizing the local estimation of the Gaussian statistics. As a consequence, we can exploit the accuracy of local model estimation for general restoration problems, in particular with missing values (e.g. for interpolation or zooming). Figure 1 illustrates the proposed approach which is described in detail in the following.

II-A Patch degradation model

The observed image $\mathbf{z}$ is decomposed into $I$ overlapping patches $\{\mathbf{z}_{i}\}_{i=1,\dots,I}$ of size $\sqrt{n}\times\sqrt{n}$. Each patch $\mathbf{z}_{i}\in\mathbb{R}^{n\times 1}$ is considered to be a realization of the random variable $\mathbf{Z}_{i}$ given by

[TABLE]

where $\boldsymbol{D}_{i}\in\mathbb{R}^{n\times n}$ is a degradation operator, $\mathbf{C}_{i}\in\mathbb{R}^{n\times 1}$ is the original patch we seek to estimate and $\mathbf{N}_{i}\in\mathbb{R}^{n\times 1}$ is an additive noise term, modeled by a Gaussian distribution $\mathbf{N}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma}_{N_{i}})$. Therefore, the distribution of $\mathbf{Z}_{i}$ given $\mathbf{C}_{i}$ can be written as

[TABLE]

In this noise model, the matrix $\boldsymbol{\Sigma}_{N_{i}}$ is only assumed to be diagonal (the noise is uncorrelated). It can represent a constant variance, spatially variable variances or even variances dependent on the pixel value (to approximate Poisson noise).

This degradation model is deliberately generic. We will see in Section IV that keeping a broad noise model is essential to properly tackle the problem of HDR imaging from a single image. The model also includes the special case of multiplicative noise.

II-B *Joint Maximum A Posteriori *

We assume a Gaussian prior for each patch, with unknown mean $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$, $p(\mathbf{C}_{i}\;|\;\boldsymbol{\mu},\boldsymbol{\Sigma})\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}).$ To simplify calculations we work with the precision matrix $\boldsymbol{\Lambda}=\boldsymbol{\Sigma}^{-1}$. As it is usual when considering hyperpriors, we assume that the parameters $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ follow a conjugate distribution. In our case, that boils down to assuming a Normal-Wishart222The Normal-Wishart distribution is the conjugate prior of a multivariate normal distribution with unknown mean and covariance matrix. $\mathcal{W}$ denotes the Wishart distribution [19]. prior for the couple ($\boldsymbol{\mu},\boldsymbol{\Lambda}$),

[TABLE]

with parameters $\boldsymbol{\mu}_{0}$, $\boldsymbol{\Sigma}_{0}$, scale parameter $\kappa>0$ and $\nu>n-1$ degrees of freedom.

Now, assume that we observe a group $\{\mathbf{Z}_{i}\}_{i=1,\dots,M}$ of similar patches and that we want to recover the restored patches $\{\mathbf{C}_{i}\}_{i=1,\dots,M}$. If these unknown $\{\mathbf{C}_{i}\}$ are independent333We rely on the classical independence assumption made in the patch-based literature, even if it is wrong in case patches overlap. and follow the same Gaussian model, we can compute the joint maximum a posteriori

[TABLE]

In this product, the first term is given by the noise model (3), the second one is the Gaussian prior on the set of patches $\{\mathbf{C}_{i}\}$ and the third one is the hyperprior (4). In the last equality we omit the explicit dependence on $\boldsymbol{\mu}$ and $\boldsymbol{\Lambda}$ in $p(\{\mathbf{Z}_{i}\}\;|\;\{\mathbf{C}_{i}\},\boldsymbol{\mu},\boldsymbol{\Lambda})$, since these parameters are completely determined by the set $\{\mathbf{C}_{i}\}$.

II-C Optimality conditions

Computing the joint maximum a posteriori amounts to minimizing

[TABLE]

over the set $\mathbb{R}^{nM}\times\mathbb{R}^{n}\times S_{n}^{++}(\mathbb{R})$, with $S_{n}^{++}(\mathbb{R})$ the set of real symmetric positive definite matrices of size $n$.

The function $f$ is biconvex respectively in the variables $(\{\mathbf{C}_{i}\},\boldsymbol{\mu})$ and $\mathbf{\boldsymbol{\Lambda}}$. To minimize this energy for a given set of hyper-parameters $(\boldsymbol{\mu}_{0},\boldsymbol{\Lambda}_{0})$, we will use an alternating convex minimization scheme. At each iteration, $f$ is first minimized with respect to $(\{\mathbf{C}_{i}\},\boldsymbol{\mu})$ with $\boldsymbol{\Lambda}$ fixed, then viceversa.

Differentiating $f$ with respect to each variable, we get explicit optimality equations for the minimization scheme. The proofs of the following propositions are straightforward and available in the supplementary material.

Proposition 1**.**

Assume that $\boldsymbol{\Lambda}$ is fixed and that the covariance $\boldsymbol{\Sigma}_{N_{i}}$ does not depend on the $\{\mathbf{C}_{i}\}$. The function $(\{{\mathbf{C}_{i}}\},\boldsymbol{\mu})\mapsto f(\{{\mathbf{C}_{i}}\},\boldsymbol{\mu},\boldsymbol{\Lambda})$ is convex on $\mathbb{R}^{n(M+1)}$ and attains its minimum at $(\{\hat{\mathbf{C}_{i}}\},\hat{\boldsymbol{\mu}})$, given by

[TABLE]

[TABLE]

with $\mathbf{A}_{i}=\boldsymbol{\Lambda}^{-1}\boldsymbol{D}_{i}^{T}(\boldsymbol{D}_{i}\boldsymbol{\Lambda}^{-1}\boldsymbol{D}_{i}^{T}+\boldsymbol{\Sigma}_{N_{i}})^{-1}$.

Proposition 2**.**

Assume that the variables $(\{{\mathbf{C}_{i}}\},\boldsymbol{\mu})$ are fixed. The function $\boldsymbol{\Lambda}\rightarrow f(\{{\mathbf{C}_{i}}\},\boldsymbol{\mu},\boldsymbol{\Lambda})$ is convex on $S_{n}^{++}(\mathbb{R})$ and attains its minimum at $\hat{\boldsymbol{\Lambda}}$ such that

[TABLE]

The expression of $\widehat{\mathbf{C}}_{i}$ in (7) is obtained under the hypothesis that the noise covariance matrix $\mathbf{\Sigma}_{N_{i}}$ does not depend on $\mathbf{C}_{i}$. Under the somewhat weaker hypothesis that the noise $N_{i}$ and the signal $\mathbf{C}_{i}$ are uncorrelated, this estimator is also the affine estimator $\tilde{\mathbf{C}}_{i}$ that minimizes the Bayes risk $\mathbb{E}[(\tilde{\mathbf{C}}_{i}-\mathbf{C}_{i})^{2}]$ (c.f. supplementary material). The uncorrelatedness of $N_{i}$ and $C_{i}$ is a reasonable hypothesis in practice. This includes various noise models, such as $N_{i}=f(C_{i})\varepsilon_{i}$ with $\varepsilon_{i}$ independent of $C_{i}$, which approximates CMOS and CCD raw data noise [20].

From (6), we find that the MAP estimator of $\boldsymbol{\mu}$ is a weighted average of two terms: the mean estimated from the similar restored patches and the prior $\boldsymbol{\mu}_{0}$. The parameter $\kappa$ controls the confidence level we have on the prior $\boldsymbol{\mu}_{0}$. With the same idea, we observe that the MAP estimator for $\boldsymbol{\Lambda}$ is a combination of the prior $\boldsymbol{\Lambda}_{0}$ on $\boldsymbol{\Lambda}$, the covariance imposed by ${\boldsymbol{\mu}}$ and the covariance matrix estimated from the patches ${\mathbf{C}}_{i}$.

II-D Alternating convex minimization of $f$

The previous propositions imply that we can derive an explicit alternating convex minimization scheme for $f$, presented in Algorithm 1. Starting with a given value of $\boldsymbol{\Lambda}$, at each step, $\boldsymbol{\mu}^{l}$ and $\mathbf{C}^{l}$ are computed according to Equations (6) and  (7), then $\boldsymbol{\Lambda}^{l}$ is updated according to (8).

We show in Appendix A the following convergence result for the previous algorithm. The proof adapts the arguments in [21] to our particular case.

Proposition 3**.**

The sequence $f(\{\mathbf{C}_{i}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})$ converges monotonically when $l\rightarrow+\infty$. The sequence $\{\{\mathbf{C}_{i}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l}\}$ generated by the alternate minimization scheme has at least one accumulation point. The set of its accumulation points forms a connected and compact set of partial optima and stationary points of $f$, and they all have the same function value.

In practice, we observe in our experiments that the algorithm always converges after a few iterations.

II-E Full restoration algorithm

The full restoration algorithm used in our experiments is summarized in Algorithm 2 and illustrated by Figure 1. It alternates between two stages: the minimization of $f$ using Algorithm 1, and the estimation of the hyper-parameters $\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0}$. In order to estimate these parameters, we rely on an oracle image computed by aggregation of all the patches estimated on the first stage (details are provided in Section III-D).

III Implementation details

III-A Search for similar patches

The similar patches are all patches within a search window centered at the current patch, whose $L_{2}$ distance to the central patch is less than a given threshold. This threshold is given by a tolerance parameter $\varepsilon$ times the distance to the nearest neighbour (the most similar one). In all our experiments, the search window was set to size $25\times 25$ (with a patch size of $8\times 8$) and $\varepsilon=1.5$. The patch comparison is performed in an oracle image (i.e. the result of the previous iteration), so all pixels are known. However, it may be useful to assign different confidence levels to the known pixels and to those originally missing and then restored. For all the experimental results presented in Section IV, the distance between patches $\mathbf{c}_{p}$ and $\mathbf{c}_{q}$ in the oracle image $\mathbf{C}_{\mathrm{oracle}}$ is computed as

[TABLE]

where $j$ indexes the pixels in the patch, $\omega^{j}_{p,q}=1$ if $\boldsymbol{D}^{j}_{p}=\boldsymbol{D}^{j}_{q}=1$ (known pixel) and $\omega^{j}_{p,q}=0.01$ otherwise (originally missing then restored pixel) [22]. With this formulation, known pixels are assigned a much higher priority than unknown ones. Variations on these weights could be explored.

III-B Optional speed-up by Collaborative Filtering

The proposed method computes one Gaussian model per image patch according to Equations (6) and (8). In order to reduce the computational cost, we can rely on the collaborative filtering idea previously introduced for patch-based denoising techniques [7, 23]. Based on the hypothesis that similar patches share the same model, we assign the same model to all patches in the set of similar patches (as defined in Section III-A).

The number of similar patches jointly restored depends on the image and the tolerance parameter $\varepsilon$, but it is often much smaller than what would result from the patch clustering performed by methods that use global GMMs such as PLE or EPLL. Performance degradation is observed in practice when using a very large tolerance parameter ($\varepsilon>3$), showing that mixing more patches than needed is detrimental. The collaborative filtering strategy helps accelerating the algorithm up to a certain point, but a trade-off with performance needs to be considered.

III-C Parameter choices

The four parameters of the Normal-Wishart distribution: $\kappa$, $\nu$, the prior mean $\boldsymbol{\mu}_{0}$ and the prior covariance matrix $\boldsymbol{\Sigma}_{0}$, must be set in order to compute $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$.

Choice of $\kappa$ and $\nu$

The computation of $\boldsymbol{\mu}$ according to (6) combines the mean $\sum_{i=1}^{M}\mathbf{A}_{i}\mathbf{Z}_{i}$ estimated from the similar patches and the prior mean $\boldsymbol{\mu}_{0}$. The parameter $\kappa$ is related to the degree of confidence we have on the prior $\boldsymbol{\mu}_{0}$. Hence, its value should be a trade-off between the confidence we have in the prior accuracy vs. the one we have in the information provided by the similar patches. The latter improves when both $M$ (i.e. the number of similar patches) and $P=\operatorname{trace}(\boldsymbol{D}_{i})$ (i.e. the number of known pixels in the current patch) increase. These intuitive insights suggest the following rule to set $\kappa$:

[TABLE]

A similar reasoning leads to the same rule for $\nu$,

[TABLE]

where the addition of $n$ ensures the condition $\nu>n-1$ required by the Normal-Wishart prior to be verified.

This rule is used to obtain the experimental results presented in Section IV, and proved to be a consistently good choice despite its simplicity. However, setting these parameters in a more general setting is not a trivial task and should be the subject of further study. In particular we could explore a more continuous dependence of $\alpha$ on $P$, $M$, and possibly a third term $Q=\sum_{i=1}^{n}S_{ii}$ where $S=\sum_{j=1}^{M}\boldsymbol{\Lambda}^{-1}\boldsymbol{D}_{j}\boldsymbol{\Lambda}^{*}_{j}\boldsymbol{D}_{j}$. This term estimates to what an extent similar patches cover the missing pixels in the current patch.

Setting of $\boldsymbol{\mu}_{0}$ and $\boldsymbol{\Sigma}_{0}$

Assuming an oracle image $\mathbf{C}_{\mathrm{oracle}}$ is available (see details in Section II-E), $\boldsymbol{\mu}_{0}$ and $\boldsymbol{\Sigma}_{0}$ can be computed using the classical MLE estimators from a set of similar patches $(\tilde{\mathbf{c}}_{1},\dots,\tilde{\mathbf{c}}_{M})$ taken from $\mathbf{C}_{\mathrm{oracle}}$

[TABLE]

This is the same approach followed by Lebrun et al. [7] to estimate the patch model parameters in the case of denoising.

III-D Initialization

A good initialization is crucial since we aim at solving a non-convex problem through an iterative procedure. To initialize the proposed algorithm we follow the approach proposed by Yu et al. [10] (described in detail in Appendix A in the supplementary material). They propose to initialize the PLE algorithm by learning the $K$ GMM covariance matrices from synthetic images of edges with different orientations as well as from the DCT basis to represent isotropic patterns. As they state, in dictionary learning, the most prominent atoms represent local edges which are useful to represent and restore contours. Hence, this initialization helps to correctly restore corrupted patches even in quite extreme cases. The oracle of the first iteration of the proposed approach is the output of the first iteration of the PLE algorithm.

III-E Computational complexity

With the original per-patch strategy, the complexity of the algorithm is given by step 3 in Algorithm 1: $[(4n^{3}+n^{3}/3)M+n^{3}/3]$, so the total complexity is $[(4n^{3}+n^{3}/3)M+n^{3}/3]\times maxIts\times maxIts_{2}\times T$ (where $T=$ total number of patches to be restored and assuming the Cholesky factorization is used for matrix inversion). The collaborative filtering strategy reduces this value by a factor that depends on the number of groups of similar patches, which depends on the image contents and the distance tolerance parameter $\varepsilon$. The main difference with the PLE algorithm complexity ($(3n^{3}+n^{3}/3)\times its_{PLE}\times T$) is a factor given by the number of groups defined by the collaborative filtering approach and the ratio between $its_{PLE}$ and $maxIts\times maxIts_{2}$. As mentioned by Yu et al. [10], computational complexity can be further reduced in the case of binary masks by removing the zero rows and inverting a matrix of size $n^{2}/S\times n^{2}/S$ instead of $n^{2}\times n^{2}$ where $S$ is the masking ratio. Moreover, the proposed algorithm can be run in parallel in different image subregions thus allowing for even further acceleration in multiple-core architectures. The complexity comparison to NLB needs to be made in the case where the degradation is additive noise with constant variance (translation invariant degradation), which is the task performed by NLB. In that case, the complexity of the proposed approach (without considering collaborative filtering nor parallelization, which are both done also in NLB), is $11n^{3}/3\times maxIts\times maxIts_{2}\times T$ whereas that of NLB is $2\times(4n^{3}/3)$.

IV Image Restoration Experiments

In this section we illustrate the ability of the proposed method to solve several image inverse problems. Both synthetic (i.e., where we have added the degradation artificially) and real data (i.e., issued from a real acquisition process) are used. The considered problems are: interpolation, combined interpolation and denoising, denoising, and zooming. The reported values of peak signal-to-noise ratio ($PSNR{}=20\log_{10}(255/\sqrt{MSE})$) are averaged over 10 realizations for each experiment (variance is below 0.1 for interpolation and below 0.05 for combined interpolation and denoising and for denoising only). Similar results are obtained with the structural similarity index (SSIM) which is included in the supplementary material (Appendix B).

IV-A Synthetic degradation

Interpolation

Random masks with 50%, 70% and 90% missing pixels are applied to the tested ground-truth images. The interpolation performance of the proposed method is compared to that of PLE [10], EPLL [9] and E-PLE [25] using a patch size of $8\times 8$ for all methods. PLE parameters are set as indicated in [10] ($\sigma=3$, $\varepsilon=30$, $K=19$). We used the EPLL code provided by the authors [24] with default parameters and the E-PLE code available in [25] with the parameters set as specified in the companion demo. The parameters for the proposed method are set to $\alpha_{H}=1$, $\alpha_{L}=0.5$ ($\alpha_{H}$ and $\alpha_{L}$ define the values for $\kappa$ and $\nu$, see Section III-C). The PSNR results are shown in Table I. Figure 2 shows some extracts of the obtained results, the PSNR values for the extracts and the corresponding difference images with respect to the ground-truth. The proposed method gives sharper results than the other considered methods. This is specially noticeable on the reconstruction of the texture of the fabric of Barbara’s trousers shown in the first row of Figure 2 or on the strips that appear through the car’s window shown in the second row of the same figure.

Combined interpolation and denoising

For this experiment, the ground-truth images are corrupted with additive Gaussian noise with variance 10, and a random mask with 70% and 90% of missing pixels. The parameters for all methods are set as in the previous interpolation-only experiment. Table I summarizes the PSNR values obtained by each method. Figure 3 shows some extracts of the obtained results, the PSNR values for the extracts and the corresponding difference images with respect to the ground-truth. Once again, the results show that the proposed approach outperforms the others. Fine structures, such as the mast and the ropes of the ship, as well as textures, as in Barbara’s headscarf, are much better preserved.

Denoising

For the denoising task the proposed approach should perform very similarly to the state-of-the-art denoising algorithm NLB [7]. The following experiments are conducted in order to verify this.

The ground-truth images are corrupted with additive Gaussian noise with variance $\sigma^{2}=10,30,50,80$. The code provided by the authors [26] automatically sets the NLB parameters from the input $\sigma^{2}$ and the patch size, in this case $8\times 8$. For this experiment, there are no unknown pixels to interpolate (the mask $\boldsymbol{D}$ is the identity matrix).

The results of both methods are very similar if HBE is initialized with the output of the first step of NLB [7] (instead of using the initialization described in Section II-E) and the parameters $\kappa$ and $\nu$ are large enough. In this case, $\boldsymbol{\mu}_{0}$ and $\boldsymbol{\Sigma}_{0}$ are prioritized in equations (6) and (8) and both algorithms are almost the same. That is what we observe in practice with $\alpha_{H}=\alpha_{L}=100$, as demonstrated in the results summarized in Table I. The denoising performance of HBE is degraded for small $\kappa$ and $\nu$ values. This is due to the fact that $\boldsymbol{\mu}_{0}$ and $\boldsymbol{\Sigma}_{0}$, as well as $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ in NLB, are computed from an oracle image resulting from the first restoration step. This restoration includes not only the denoising of each patch, but also an aggregation step that greatly improves the final result. Therefore, the contribution of the first term of (6) to the computation of $\hat{\boldsymbol{\mu}}$ degrades the result compared to that of using $\boldsymbol{\mu}_{0}$ only (i.e. using a large $\kappa$).

Zooming

In order to evaluate the zooming capacity of the proposed approach, ground-truth images are downsampled by a factor 2 (no anti-aliasing filter is used) and the zooming is compared to the ground-truth. The results are compared with PLE, EPLL, E-PLE and Lanczos interpolation. Table I summarizes the obtained PSNR values. Figure 4 shows extracts from the obtained results, the PSNR values for the extracts and the corresponding difference images with respect to the ground-truth. Again, HBE yields a sharper reconstruction than the other methods.

IV-B Real data

For this experiment, we use raw images captured with a Canon 400D camera (set to ISO 400 and exposure time 1/160 seconds). The main noise sources for CMOS sensors are: the Poisson photon shot noise, which can be approximated by a Gaussian distribution with equal mean and variance; the thermally generated readout noise, which is modeled as an additive Gaussian distributed noise and the spatially varying gain given by the photo response non uniformity (prnu) [20, 27]. We thus consider the following noise model for the non saturated raw pixel value $\mathrm{Z}(p)$ at position $p$

[TABLE]

where $g$ is the camera gain, $a_{p}$ models the prnu factor, $\tau$ is the exposure time, $\mathrm{C}(p)$ is the irradiance reaching pixel $p$, $\boldsymbol{\mu}_{R}$ and $\sigma^{2}_{R}$ are the readout noise mean and variance. The camera parameters have to be estimated by a calibration procedure [20]. The noise covariance matrix $\boldsymbol{\Sigma}_{N}$ is thus diagonal with entries that depend on the pixel value $(\boldsymbol{\Sigma}_{N})_{p}=g^{2}a_{p}\tau\mathrm{C}(p)+\sigma^{2}_{R}$.

In order to evaluate the interpolation capacity of the proposed approach, we consider the pixels of the green channel only (i.e. 50% of the pixels in the RGGB Bayer pattern) and interpolate the missing values. We compare the results to those obtained using an adaptation of PLE to images degraded with noise with variable variance (PLEV) [28]. The results for the EPLL and E-PLE methods are not presented here since these methods are not suited for this kind of noise. Figure 6 shows extracts of the obtained results (see Figure 5 for a JPEG version of the raw image showing the location of the extracts). As it was already observed in the synthetic data experiments, fine details and edges are better preserved. Compare for example the reconstruction of the balcony edges and the wall structure in the first row of Figure 6, as well as the structure of the roof and the railing in the second row of the same image.

IV-C Discussion

In all the previous experiments, the results obtained with HBE outperform or are very close to those obtained by the other evaluated methods. Visually, details are better reconstructed and images are sharper.

We interpret this as the inability of a fixed set of patch classes (19 for PLE) to accurately represent patches that seldom appear in the image, such as edges or textures (as in Barbara’s trouser). The fact that methods such as PLE are actually semi-local (classes are estimated on $128\times 128$ regions [10]) does not solve this issue. On the contrary, a local model estimation as the one performed by HBE correctly handles those cases.

The performance difference is much more remarkable for the higher masking rates. In such cases, the robustness of the estimation is crucial. Indeed the proposed method performs the estimation from similar patches in a local window. The hypothesis of self-similarity being reinforced by considering local neighbourhoods, such a strategy restricts the possible models, therefore making the estimation more robust. Furthermore, the local model estimation, previously shown to be successful at describing patches [7], gives a better reconstruction even when a very large part of the patch is missing.

EPLL uses more mixture components (200 components are learnt from $2\times 10^{6}$ patches of natural images [9]) in its GMM model than PLE. It was observed in [8] that this strategy is less efficient than PLE for the denoising task. In this work, we also observe that the proposed approach outperforms EPLL, not only in denoising, but also in inpainting and zooming. However, here it is harder to tell if the improvement is due to the local model estimation performed from similar patches or to the different restoration strategies followed by these methods.

V High dynamic range imaging from a single snapshot

In this section, we apply the proposed framework to generate HDR images from a single shot. HDR imaging aims at reproducing an extended dynamic range of luminosity compared to what can be captured using a standard digital camera, which is often not enough to produce an accurate representation of real scenes. In the case of a static scene and a static camera, the combination of multiple images with different exposure levels is a simple and efficient solution [29, 30, 27]. However, several problems arise when either the camera or the elements in the scene move [31, 32].

An alternative to the HDR from multiple frames is to use a single image with spatially varying pixel exposures (SVE) [18]. An optical mask with spatially varying transmittance is placed adjacent to a conventional image sensor, thus controlling the amount of light that reaches each pixel (see Figure 7) [18, 33, 34].

The greatest advantage of this acquisition method is that it avoids the need for image alignment and motion estimation. Another advantage is that the saturated pixels are not organized in large regions. Indeed, some recent multi-image methods tackle motion problems by taking a reference image and then by estimating motion or reconstructing the image relative to this reference [35, 31]. A problem encountered by these approaches is the need to inpaint very large saturated and underexposed regions in the reference frame. The SVE acquisition strategy avoids this problem since, in general, all scene regions are sampled by at least one of the exposures.

Taking advantage of the ability of the proposed framework to simultaneously estimate missing pixels and denoise well-exposed ones, we propose a novel approach to generate HDR images from a single shot acquired with spatially varying pixel exposures. The proposed approach shows significant improvements over existing methods.

V-A Spatially varying exposure acquisition model

An optical mask with spatially varying transmittance is placed adjacent to a conventional image sensor to give different exposure levels to the pixels. This optical mask does not change the acquisition process of the sensor. Hence, the noise model (13) can be adapted to the SVE acquisition by including the per-pixel SVE gain $o_{p}$444Some noise sources not modeled here, such as blooming, might have an impact in the SVE acquisition strategy and should be considered in a more accurate image modeling.:

[TABLE]

In the approach proposed by Nayar and Mitsunaga [18], the varying exposures follow a regular pattern. Motivated by the aliasing problems of regular sampling patterns, Schöberl et al. [36] propose to use spatially varying exposures on a non-regular pattern. Figure 7 shows examples of both acquisition patterns. This observation led us to choose the non-regular pattern in the proposed approach.

V-B Hyperprior Bayesian Estimator for Single Shot High Dynamic Range Imaging

V-B1 Problem statement

In order to reconstruct the dynamic range of the scene we need to solve an inverse problem. We want to estimate the irradiance image $C$ from the SVE image $\mathbf{Z}$, knowing the exposure levels of the optical mask and the camera parameters. For this purpose we map the raw pixel values to the irradiance domain $\mathrm{Y}$ with

[TABLE]

We take into account the effect of saturation and under-exposure by introducing the exposure degradation matrix $\boldsymbol{D}$, a diagonal matrix given by

[TABLE]

with $z_{sat}$ equal to the pixel saturation value. From (14) and (16), $\mathrm{Y}(p)$ can be modeled as

[TABLE]

Notice that (17) is the distribution of $\mathrm{Y}(p)$ for a given exposure degradation factor $(\boldsymbol{D})_{p}$, since $(\boldsymbol{D})_{p}$ is itself a random variable that depends on $\mathrm{Z}(p)$. The exposure degradation factor must be included in (17) since the variance of the over or under exposed pixels no longer depends on the irradiance $\mathrm{C}(p)$ but is only due to the readout noise $\sigma^{2}_{R}$. From (17) we have

[TABLE]

where $\mathbf{N}$ is zero-mean Gaussian noise with diagonal covariance matrix $\boldsymbol{\Sigma}_{\mathbf{N}}$ given by

[TABLE]

Then the problem of irradiance estimation can be stated as retrieving $C$ from $\mathrm{Y}$, which implies denoising the well-exposed pixel values ($(\boldsymbol{D})_{p}=1$) and estimating the unknown ones ($(\boldsymbol{D})_{p}=0$).

V-B2 Proposed solution

From (18), image $\mathrm{Y}$ is under the hypothesis of the HBE framework introduced in Section II. The proposed HDR imaging algorithm consists of the following steps: 1) generate $\boldsymbol{D}$ from $\mathrm{Z}$ according to (16), 2) obtain $\mathrm{Y}$ from $\mathrm{Z}$ according to (15), 3) apply the HBE approach to $\mathrm{Y}$ with the given $\boldsymbol{D}$ and $\boldsymbol{\Sigma}_{\mathbf{N}}$.

V-C Experiments

The proposed reconstruction method was thoroughly tested in several synthetic and real data examples. A brief summary of the results is presented in this section.

V-C1 Synthetic data

Sample images are generated according to Model (18) using the HDR image in Figure 7 as the ground-truth. Both a random and a regular pattern with four equiprobable exposure levels $o=\{1,8,64,512\}$ are simulated. The exposure time is set to $\tau=1/200$ seconds and the camera parameters are those of a Canon 7D camera set to ISO 200 ($g=0.87$, $\sigma^{2}_{R}=30$, $\boldsymbol{\mu}_{R}=2048$, $\text{v}_{\text{sat}}=15000$) [27].

Figure 7 shows extracts of the results obtained by the proposed method, by PLEV [28] (basically an adaptation of PLE to the same single image framework) and by Schöberl et al. [36] for the random pattern and by Nayar et Mitsunaga [18] using the regular pattern. The percentage of unknown pixels in the considered extracts is 50% (it is nearly the same for both the regular and non-regular pattern). Table II shows the PSNR values obtained in each extract marked in Figure 7. The proposed method manages to correctly reconstruct the irradiance on the unknown pixels. Moreover, its denoising performance is much better than that of Schöberl et al. and Nayar and Mitsunaga, and still sharper than PLEV.

V-C2 Real data

The feasibility of the SVE random pattern has been shown in [34] and that of the SVE regular pattern in [33]. Nevertheless, these acquisition systems are still not available for general usage.555 While writing the last version of this article the authors got aware of Sony’s latest sensor IMX378. This sensor has a special mode called SME-HDR, which is a variation of the SVE acquisition principle. Whereas this sensor was adopted by the Google Pixel smartphone in 2016, the special SME-HDR mode is never activated by the Google Pixel phone, according to experts from the company DxO, and we found no way to activate it and access the raw image.

However, as stated in Section V-A, the only variation between the classical and the SVE acquisition is the optical filter. Hence, the noise at a pixel $p$ captured using SVE with an optical gain factor $o_{p}$ and exposure time $\tau/o_{p}$ and a pixel captured with a classical camera using exposure time $\tau$ should be very close. We take advantage of this fact in order to evaluate the reconstruction performance of the proposed approach using real data.

For this purpose, we generate an SVE image $\mathbf{z}_{sve}$ from four raw images $\{\mathbf{z}_{raw}^{i}\}_{i=1,\dots,4}$ acquired with different exposure times. The four different exposure times simulate four different filters of the SVE optical mask. The value at position $(x,y)$ in $\mathbf{z}_{sve}$ is chosen at random among the four available values at that position $\{\mathbf{z}_{raw}^{i}(x,y)\}_{i=1,\dots,4}$. Notice that the Bayer pattern is kept on $\mathbf{z}_{sve}$ by construction. The images $\{\mathbf{z}_{raw}^{i}\}_{i=1,\dots,4}$ are acquired using a remotely controlled camera and a tripod so as to be perfectly aligned. This protocol does not allow us to take scenes with moving objects. Let us emphasize, however, that using a real SVE device, this, as well as the treatment of moving camera, would be a non-issue.

Figures 8 and 9 show the results obtained from two real scenes, together with the masks of well-exposed (white) and unknown (black) pixels (the SVE raw images are included in Appendix B in the supplementary material). Recall that among the unknown pixels, some of them are saturated and some of them are under exposed. Square patches of size 6 and 8 were used for the examples in Figure 9 and Figure 8 respectively. Demosaicing [37] and tone mapping [38] are used for displaying purposes.

We compare the results to those obtained by PLEV [28]. A comparison against the methods by Nayar and Mitsunaga and Schöberl et al. is not presented since they do not specify how to treat raw images with a Bayer pattern. The proposed method manages to correctly reconstruct the unknown pixels even in extreme conditions where more than $70\%$ of the pixels are missing, as for example the last extract in Figure 9.

These examples show the suitability of the proposed approach to reconstruct the irradiance information in both very dark and bright regions simultaneously. See for instance the example in Figure 9, where the dark interior of the building (which can be seen through the windows) and the highly illuminated part of another building are both correctly reconstructed (see the electronic version of the article for better visualization).

VI Conclusions

In this work we have presented a novel image restoration framework. It has the benefits of local patch characterization (that was key to the success of NLB as a state of the art denoising method), but manages to extend its use to more general restoration problems where the linear degradation operator is diagonal, by combining local estimation with Bayesian restoration based on hyperpriors. This includes problems such as zooming, inpainting and interpolation. In this way, all these restoration problems are set under the same framework. It does not include image deblurring or deconvolution, since the degradation operator is no longer diagonal. Correctly addressing deconvolution with large kernels with patch-based approaches and Gaussian prior models is a major challenge that will be the subject of future work.

We have presented a large series of experiments both on synthetic and real data that confirm the robustness of the proposed strategy based on hyperpriors. These experiments show that for a wide range of image restoration problems HBE outperforms several state-of-the-art restoration methods.

This work opens several perspectives. The first one concerns the relevance of the Gaussian patch model and its relation to the underlying image patches manifold. If this linear approximation has proven successful for image restoration, its full relevance in other areas remains to be explored, especially in all domains requiring to compare image patches. Another important related question is the one of the estimation of the degradation model in images jointly degraded by noise, missing pixels, blur, etc. Restoration approaches generally rely on the precise knowledge of this model and of its parameters. In practice however, we often deal with images for which the acquisition process is unknown, and that have possibly been affected by post-treatments. In such cases, blind restoration remains an unsolved challenge.

Finally, we have presented a novel application of the proposed general framework to the generation of HDR images from a single variable exposure (SVE) snapshot. The SVE acquisition strategy allows the creation of HDR images from a single shot without the drawbacks of multi-image approaches, such as the need for global alignment and motion estimation to avoid ghosting problems. The proposed method manages to simultaneously denoise and reconstruct the missing pixels, even in the presence of (possibly complex) motions, improving the results obtained by existing methods. Examples with real data acquired in very similar conditions to those of the SVE acquisition show the capabilities of the proposed approach.

Appendix A Alternate minimization scheme convergence

We study in the following the convergence of the alternate minimization algorithm 1. To show the main convergence result, we need the following lemma

Lemma 1**.**

The function $f$ is coercive on $\mathbb{R}^{n(M+1)}\times S_{n}^{++}(\mathbb{R})$.

Proof.

We need to show that

[TABLE]

Now, $\|(\{\mathbf{C}_{i}\},\boldsymbol{\mu},\boldsymbol{\Lambda})\|\rightarrow+\infty$ if and only if $\|{\mathbf{C}_{i}}\|\rightarrow+\infty$ or $\|\boldsymbol{\mu}\|\rightarrow+\infty$ or $\|\boldsymbol{\Lambda}\|\rightarrow+\infty$.

The matrix $\boldsymbol{\Lambda}$ being positive-definite, the terms $\frac{1}{2}\sum_{i=1}^{M}(\mathbf{C}_{i}-\boldsymbol{\mu})^{T}\boldsymbol{\Lambda}(\mathbf{C}_{i}-\boldsymbol{\mu})$ and $\frac{\kappa}{2}(\boldsymbol{\mu}-\boldsymbol{\mu}_{0})^{T}\boldsymbol{\Lambda}(\boldsymbol{\mu}-\boldsymbol{\mu}_{0})$ are both positive. Thus

[TABLE]

Now, this function of $\boldsymbol{\Lambda}$ is convex and coercive on $S_{n}^{++}(\mathbb{R})$, which implies that $f(\{\mathbf{C}_{i}\},\boldsymbol{\mu},\boldsymbol{\Lambda})\rightarrow+\infty$ as soon as $\|\boldsymbol{\Lambda}\|\rightarrow+\infty$. It also follows that the previous function of $\boldsymbol{\Lambda}$ has a global minimum that we denote by $m_{\boldsymbol{\Lambda}}$. We can now write

[TABLE]

and this function of $(\{\mathbf{C}_{i}\},\boldsymbol{\mu})$ clearly tends towards $+\infty$ as soon as $\|{\mathbf{C}_{i}}\|\rightarrow+\infty$ or $\|\boldsymbol{\mu}\|\rightarrow+\infty$. ∎

We now show the main convergence result for our alternate minimization algorithm. The proof adapts the arguments in [21] to our case.

Proposition 3**.**

The sequence $f(\{\mathbf{C}_{i}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})$ converges monotonically when $l\rightarrow+\infty$. The sequence $\{\{\mathbf{C}_{i}^{l}\},\boldsymbol{\mu}_{l},\boldsymbol{\Lambda}^{l}\}$ generated by the alternate minimization scheme has at least one accumulation point. The set of its accumulation points forms a connected and compact set of partial optima and stationary points of $f$, all having the same function value.

Proof.

The sequence $f(\{\mathbf{C}_{i}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})$ obviously decreases at each step by construction. The coercivity and continuity of $f$ imply that this sequence is also bounded from below, and thus converges. The convergence of $f(\{\mathbf{C}_{i}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})$ implies that the sequence $\{\{\mathbf{C}_{i}^{l}\},\boldsymbol{\mu}_{l},\boldsymbol{\Lambda}^{l}\}$ is bounded. It follows that it has at least one accumulation point $(\{\mathbf{C}_{i}^{\star}\},\boldsymbol{\mu}^{\star},\boldsymbol{\Lambda}^{\star})$ and that there exists a strictly increasing sequence $(l_{k})_{k\in\mathbb{N}}$ of integers such that $\{\{\mathbf{C}_{i}^{l_{k}}\},\boldsymbol{\mu}^{l_{k}},\boldsymbol{\Lambda}^{l_{k}}\}_{k\in\mathbb{N}}$ converges to $\{\{\mathbf{C}_{i}^{\star}\},\boldsymbol{\mu}^{\star},\boldsymbol{\Lambda}^{\star}\}$.

Now, we can show that such an accumulation point is a partial optimum of $f$, i.e. that $f(\{\mathbf{C}_{i}^{\star}\},\boldsymbol{\mu}^{\star},.)$ attains its minimum at $\boldsymbol{\Lambda}^{*}$ and $f(.,.,\boldsymbol{\Lambda}^{*})$ attains its minimum at $(\{\mathbf{C}_{i}^{\star}\},\boldsymbol{\mu}^{\star})$. By construction,

[TABLE]

which implies by continuity of $f$ that

[TABLE]

Let us denote $G(\{\mathbf{C}\},\boldsymbol{\mu},\boldsymbol{\Lambda})=(\{\mathbf{C}^{\prime}\},\boldsymbol{\mu}^{\prime},\boldsymbol{\Lambda}^{\prime})$ with

[TABLE]

The alternate minimization scheme consists in updating $G$ at each iteration. From Equations (6), (7), and (8), we see that $G$ is explicit and continuous. Since $\{l_{k}\}_{k\in\mathbb{N}}$ is strictly increasing, for each $k\in{\mathbb{N}}^{*}$, $l_{k}\geq l_{k-1}+1$. The sequence $\{f(\{\mathbf{C}_{i}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})\}_{l\in\mathbb{N}}$ decreases, so

[TABLE]

Therefore, as $k\rightarrow+\infty$, since $G$ is continuous, it follows that

[TABLE]

Now, writing $(\{\mathbf{C}^{**}\},\boldsymbol{\mu}^{**},\boldsymbol{\Lambda}^{**})=G(\{\mathbf{C}^{*}\},\boldsymbol{\mu}^{*},\boldsymbol{\Lambda}^{*})$, we get

[TABLE]

We can conclude that all these terms are equal and in particular

[TABLE]

From (20) and (21) we deduce that the accumulation point $(\{\mathbf{C}^{\star}\},\boldsymbol{\mu}^{\star},\boldsymbol{\Lambda}^{\star})$ is a partial optimum of $f$ and since $f$ is differentiable, it is also a stationary point of $f$. Moreover, since $f(.,.,\boldsymbol{\Lambda}^{*})$ is strictly convex and has a unique minimum, it follows from (21) that $(\{\mathbf{C}^{**}\},\boldsymbol{\mu}^{**})=(\{\mathbf{C}^{*}\},\boldsymbol{\mu}^{*})$. As a consequence, $\boldsymbol{\Lambda}^{**}=\operatorname*{arg\,min}_{\boldsymbol{\Lambda}}f(\{\mathbf{C}^{**}\},\boldsymbol{\mu}^{**},\boldsymbol{\Lambda})=\operatorname*{arg\,min}_{\boldsymbol{\Lambda}}f(\{\mathbf{C}^{*}\},\boldsymbol{\mu}^{*},\boldsymbol{\Lambda})=\boldsymbol{\Lambda}^{*}$. Therefore, the accumulation point $(\{\mathbf{C}^{\star}\},\boldsymbol{\mu}^{\star},\boldsymbol{\Lambda}^{\star})$ (and actually any accumulation point of the sequence) is also a fixed point of function $G$.

We have shown that accumulation points of the sequence $\{\{\mathbf{C}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l}\}$ are partial optima of $f$ and fixed points of the function $G$. The set of accumulation points is obviously compact. Let us show that it is also a connected set. First, observe that whatever the norm $\|\|$, the sequence $\|(\{\mathbf{C}^{l+1}\},\boldsymbol{\mu}^{l+1},\boldsymbol{\Lambda}^{l+1})-(\{\mathbf{C}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})\|$ converges to [math] when $l\rightarrow\infty$. If it was not the case, it would be possible to extract a subsequence $(\{\mathbf{C}^{l_{k}}\},\boldsymbol{\mu}^{l_{k}},\boldsymbol{\Lambda}^{l_{k}})$ converging to an accumulation point $(\{\mathbf{C}^{*}\},\boldsymbol{\mu}^{*},\boldsymbol{\Lambda}^{*})$ while $(\{\mathbf{C}^{l_{k}+1}\},\boldsymbol{\mu}^{l_{k}+1},\boldsymbol{\Lambda}^{l_{k}+1})$ converges to a different accumulation point $(\{\mathbf{C}^{{}^{\prime}}\},\boldsymbol{\mu}^{{}^{\prime}},\boldsymbol{\Lambda}^{{}^{\prime}})$, but we know that it is impossible since $(\{\mathbf{C}^{l_{k}+1}\},\boldsymbol{\mu}^{l_{k}+1},\boldsymbol{\Lambda}^{l_{k}+1})=G(\{\mathbf{C}^{l_{k}}\},\boldsymbol{\mu}^{l_{k}},\boldsymbol{\Lambda}^{l_{k}})$ would also tend toward $(\{\mathbf{C}^{*}\},\boldsymbol{\mu}^{*},\boldsymbol{\Lambda}^{*})$. The sequence $(\{\mathbf{C}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})$ being bounded and such that $\|(\{\mathbf{C}^{l+1}\},\boldsymbol{\mu}^{l+1},\boldsymbol{\Lambda}^{l+1})-(\{\mathbf{C}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})\|$ converges to [math], the set of its accumulation points is connected (see [39]). The fact that all accumulation points have the same function value is obvious since the sequence $\{f(\{\mathbf{C}^{l}\},\boldsymbol{\mu}^{l},\boldsymbol{\Lambda}^{l})\}_{l\in\mathbb{N}}$ decreases. ∎

Acknowledgement

We would like to thank the reviewers for their thorough and fruitful contributions, as well as Mila Nikolova and Alasdair Newson for their insightful comments and the authors of [26, 25, 24] for kindly providing their code. This work has been partially funded by the French Research Agency (ANR) under grant nro ANR-14-CE27-001 (MIRIAM), and by the “Investissement d’avenir” project, reference ANR-11-LABX-0056-LMH.