Pablo Arias
email: pablo.arias@cmla.ens-cachan.fr

Jean-Michel Morel
email: jean-michel.morel@cmla.ens-cachan.fr

Video denoising via empirical Bayesian estimation of space-time patches

In this paper we present a new patch-based empirical Bayesian video denoising algorithm. The method builds a Bayesian model for each group of similar space-time patches. These patches are not motion compensated, and therefore avoid the risk of inaccuracies caused by motion estimation errors. The high dimensionality of spatio-temporal patches together with a limited number of available samples, poses challenges when estimating the statistics needed for an empirical Bayesian method. We therefore assume that groups of similar patches have a low intrinsic dimensionality, leading to a spiked covariance model. Based on theoretical results about the estimation of spiked covariance matrices, we propose estimators of the eigenvalues of the a priori covariance in high dimensional spaces as simple corrections of the eigenvalues of the sample covariance matrix. We demonstrate empirically that these estimators lead to better empirical Wiener filters. A comparison on classic benchmark videos demonstrates improved visual quality and an increased PSNR with respect to state-of-the-art video denoising methods.

Introduction

Advances in video sensor hardware have steadily improved the acquisition quality. However, due to the reduction in the price of cameras and data storage, video cameras are being used increasingly often in unfavorable conditions, such as low lighting. This results in high levels of noise, which negatively affects the visual quality of the video and hinders its use for any application.

The denoising of videos has received less attention than still image denoising in the literature. In principle, denoising videos should be easier than still images, since the strong temporal redundancy along motion trajectories favours the denoising task. However, this additional form of redundancy creates challenges. On the one hand, the amount of data is much larger than for still images, and efficient ways of navigating through this data are needed. On the other hand, the output of the denoising algorithm should be temporally consistent with the unknown motion of the original sequence, which is an essential quality criterion [START_REF] Liu | A high-quality video denoising algorithm based on reliable motion estimation[END_REF][START_REF] Seshadrinathan | Motion tuned spatio-temporal quality assessment of natural videos[END_REF].

Early works proposed temporal or spatio-temporal filters which either compensate for motion using estimated motion trajectories or used some kind of mechanism to adapt to the temporal changes in the video signal at a fixed location. We refer the reader to [START_REF] Brailean | Noise Reduction Filters for Dynamic Image Sequences: A Review[END_REF] and references therein for more details. Other video denoising works apply temporal filtering to the wavelet coefficients of the frames [START_REF] Jin | Wavelet video denoising with regularized multiresolution motion estimation[END_REF][START_REF] Zlokolica | Wavelet-domain video denoising based on reliability measures[END_REF]. Some methods do not distinguish between the temporal and spatial dimensions, and treat the video as a volume, which is denoised in a transformed domain [START_REF] Rajpoot | Adaptive wavelet restoration of noisy video sequences[END_REF][START_REF] Selesnick | Video denoising using 2D and 3D dual-tree complex wavelet transforms[END_REF].

Until the advent of patch-based approaches, methods without motion compensation failed in the presence of moderate motion, and their main interest resided in their simplicity and low computational cost. The first state-of-the-art results obtained without motion estimation were reported in [START_REF] Buades | Denoising image sequences does not require motion estimation[END_REF], with the non-local means algorithm. The authors argued that for their method, based on averaging similar patches, motion estimation was not only unnecessary but counterproductive: while self-similar regions are problematic for motion estimation, they are a source of a large number of similar patches across different frames. Even if the motion trajectory could be reliably estimated, it makes more sense to use all similar patches, rather than using only the ones on the trajectory. A similar approach was followed by the authors of [START_REF] Dabov | Video denoising by sparse 3D transform-domain collaborative filtering[END_REF]. This method is based on collaborative filtering of similar patches in a spatio-temporal neighborhood. The methods in [START_REF] Buades | Denoising image sequences does not require motion estimation[END_REF] and [START_REF] Dabov | Video denoising by sparse 3D transform-domain collaborative filtering[END_REF] are extensions to video of still-image patch-based denoising algorithms (non-local means [START_REF] Buades | A review of image denoising algorithms , with a new one[END_REF] and BM3D [START_REF] Dabov | Image denoising by sparse 3d transform-domain collaborative filtering[END_REF] respectively). They are based on filtering similar 2D patches searched for in a 3D spatio-temporal neighborhood. These methods exploit both spatial and temporal redundancy. However, since each 2D patch is processed independently, there is no mechanism to impose coherence along trajectories. This results in flickering artifacts which become noticeable for high levels of noise.

An interesting alternative was proposed in [START_REF] Liu | A high-quality video denoising algorithm based on reliable motion estimation[END_REF], based on non-local means. To denoise a target 2D patch, its forward and backward trajectories are estimated using an optical flow algorithm. The patch is denoised as a weighted average of the patches along the trajectory and their nearest neighbors. The weight considers the patch similarity like in the non-local means method, multiplied by a coefficient that decreases smoothly with the temporal distance to the target patch. In this way, 2D patches in the same trajectory are denoised using a smoothly varying non-local average of similar 2D patches, resulting in a video with a high temporal consistency.

A natural next step is to consider 3D spatio-temporal patches. This has two main advantages. First it provides a mechanism to impose some temporal consistency on the result, since then, the temporally contiguous 2D temporal slices of a spatio-temporal patch are filtered jointly. Second, a 3D patch has a higher number of pixels resulting in a reduction of the noise in the distance between patches. To get the same reduction in the distance noise with 2D patches one needs to increase their spatial size which typically leads to a smaller number of similar patches. Two types of 3D patches have been used in the literature: rectangular and motion-compensated. Motion-compensated 3D patches follow a motion trajectory previously estimated.

Motion-compensated patches have the benefit that their temporal slices are highly correlated. As a consequence, temporal filtering can be easily performed by filtering inside the patch in the temporal direction. Because of this, methods that use motion-compensated 3D patches only look for similar motion-compensated 3D patches starting in the same frame.

In [START_REF] Maggioni | Video Denoising Using Separable 4D Nonlocal Spatiotemporal Transforms[END_REF][START_REF] Maggioni | Video denoising, deblocking, and enhancement through separable 4-D nonlocal spatiotemporal transforms[END_REF] the authors introduced V-BM4D, an extension of BM3D to video by collaborative filtering of similar motion compensated 3D patches. Trajectories are computed with a block matching strategy based on the sum of squared differences (SSD) and a temporal regularization term, which favours trajectories with small velocity and low acceleration. Recently, Buades et al. [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF] proposed to denoise each frame by warping the neighboring frames to match the target frame via an optical flow method. This defines a 3D volume with almost identical frames (in the ideal case). A patch-based denoising method is then used to denoise the central frame in the warped volume. Patch similarity is determined using 3D patches in the warped volume, which is equivalent to consider motion compensated patches. The denoising of each 2D patch in the target frame results from a PCA model learnt from the 2D frames of the similar 3D patches.

The main difficulty in using motion-compensated patches is obtaining an accurate motion estimation in the presence of noise. An alternative which does not require motion estimation is to consider rectangular 3D patches. Each rectangular 3D patch results from the combination of a 2D texture and a motion pattern. For arbitrary motion it would seem that each combination of texture and motion is unlikely to repeat across the video. Thus rectangular 3D patches should have less repeatability in the video compared to motion compensated patches. This objection is mitigated in two ways. First, if the motion is uniform, rectangular 3D patches will be similar along the trajectory. Second, if the motion pattern is spatially homogeneous in a region of the image, then the image self-similarity still extends to space-time patches starting at the same frame.

Indeed, many authors have considered rectangular 3D patches. An improved version of nonlocal means, which adaptively controls the shape and size of the spatio-temporal search region is proposed in [START_REF] Boulanger | Space-time adaptation for patch-based image sequence restoration[END_REF]. The authors compare the performance obtained with 2D and 3D patches. Sutour et al. [START_REF] Sutour | Adaptive regularization of the nl-means: Application to image and video denoising[END_REF] proposed to add an adaptive spatial regularization to non-local means, and to apply it to images and videos, also with 3D patches. In [START_REF] Protter | Sparse and Redundant Representations and Motion-Estimation-Free Algorithm for Video Denoising[END_REF][START_REF] Protter | Image sequence denoising via sparse and redundant representations[END_REF] the authors extend the K-SVD [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF] image denoising method to video by learning a dictionary of spatio-temporal rectangular 3D patches. A multi-scale version was proposed in [START_REF] Mairal | Learning multiscale sparse representations for image and video restoration[END_REF]. Ghoniem et al. [START_REF] Ghoniem | Nonlocal video denoising, simplification and inpainting using discrete regularization on graphs[END_REF] model non-local interactions in a video as a weighted graph where nodes are associated to 3D patches and the weights on the edges are given by patch similarity. A graph regularization framework is then proposed and applied to several inverse problems. In [START_REF] Maggioni | A nonlocal transform-domain filter for volumetric data denoising and reconstruction[END_REF] a version of V-BM4D, called BM4D, using rectangular 3D patches is applied to the restoration of volumetric images. The authors report results on video denoising as well, with a performance inferior to that of V-BM4D, particularly in sequences with motion. A similar approach, using shape adaptive patches was proposed in [START_REF] Li | Video denoising using shape-adaptive sparse representation over similar spatiotemporal patches[END_REF] for grayscale sequences.

In this work we present a new empirical Bayesian patch-based video denoising method using 3D rectangular patches. The proposed method assumes that similar patches can be modelled as independent, identically distributed (IID) samples from an unknown a priori distribution. Each patch is denoised by computing the linear minimum mean square error estimator (LMMSE) with the estimated first and second order moments of the prior. Equivalently, the proposed approach can be described as a maximum a posteriori (MAP) under the assumption that the prior distribution is Gaussian. 1In recent years several works have proposed Gaussian models (or Gaussian mixture models) as priors for image patches, achieving state-of-the-art results in image denoising and other inverse problems. However, the potential of these techniques for video restoration remains mostly unexplored.

A Gaussian mixture model was proposed in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] as a fixed prior for image patches. The mixture comprising hundreds of Gaussians, is learnt off-line from a large patch database. In [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF] an image dependent mixture is learnt for each image. These works consider a single GMM for the whole image. Extending these approaches to video is challenging because of the high variability of spatio-temporal patches. Thus we shall stick to the idea of using video redundancy to build a local Gaussian model for each patch from its most similar patches as in [START_REF] Aguerrebere | Inverse Problems in Imaging: a Hyperprior Bayesian Approach[END_REF][START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Zhang | Two-stage image denoising by principal component analysis with local pixel grouping[END_REF]. Such an approach is applied in [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF] to video denoising, but restricted to 2D patches. The authors argue that the main limitation for considering 3D patches is that, as the dimensionality of the patch grows, the number of similar patches required to estimate the covariance matrix grows larger than the number of similar patches available.

However, we can expect that similar patches have only a few modes of variation, resulting in a low intrinsic dimensionality, or equivalently a low-rank covariance matrix. Under this assumption we propose estimators of a low-rank a priori covariance matrix in high dimensional spaces, which are motivated by recent theoretical results and supported by an empirical analysis. We apply these estimators to derive empirical Wiener filters adapted to sets of similar 3D patches. We tested with 3D patches of dimensionality d = 100 and d = 200 and different combinations of spatial and temporal sizes. The proposed approach produces results which compare favourably to the current state-of-the-art, both quantitatively and qualitatively, specifically in terms of temporal consistency.

A preliminary version of this work using 2D patches was presented in [START_REF] Arias | Towards a Bayesian video denoising method[END_REF]. Here we generalize it by considering 3D patches together with an improved estimation of the eigenvalues of the low-rank covariance matrices.

The rest of the paper is organized as follows: in §2 we describe the proposed algorithm. In §3 we address the problem of estimating a low-rank covariance matrix in high dimensional spaces and propose two simple estimators derived from the asymptotic results of [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF]. In §4 we explain the selection of the parameters. Some results, including a comparison with the state-of-the-art methods, are shown in §5. Concluding remarks are given in §6.

Bayesian video denoising

Our main assumption, following [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF], is that video patches similar to a given patch are independent random samples drawn from some a priori distribution. The problem of denoising each patch can then be formulated as a problem of optimal Bayesian inference, where the parameters of the prior are learnt from the noisy data.

A non-local Bayesian principle

Consider a grayscale video u :

Ω × {1, • • • , T } → R,
where Ω is the spatial domain (a rectangular discrete grid). We assume that v is a noisy version of u, contaminated with additive, Gaussian, white noise (AGWN) z of known variance σ 2 ,

v = u + z.
We denote a rectangular 3D patch of size s x × s x × s t of the noisy video v as q, and the corresponding patch in the clean video u as p. We consider the patches as vectors in R d , with d = s 2 x s t . We have q = p + z, where z ∼ N (0, σ 2 I d ) is a patch of AGWN (I d denotes the d × d identity matrix).

For a noisy patch q, we select the n closest 3D patches q 1 = q, q 2 , . . . , q n according to the Euclidean distance in R d . The search is conducted in a spatio-temporal region centered at q. The spatial section of the search region is of size w x × w x and it extends for w t frames. The main assumption justifying a Bayesian method is that the corresponding patches in the clean video, p 1 , p 2 , . . . , p n , are independent random samples from an a priori distribution with mean E{p} = µ and covariance

E{(p -µ)(p -µ) T } = C.
We can then estimate the clean patches p i using the linear estimator minimizing the expected mean square error (MSE), or LMMSE:

p i = µ + C(C + σ 2 I d ) -1 (q i -µ) (1)
If we assume that the a priori distribution is Gaussian, then this estimator minimizes the MSE across all estimators, and coincides also with the maximum-a-posteriori (MAP) estimator.

The matrix C(C + σ 2 I d ) -1 in the above equation is the known Wiener filter and can be diagonalized in the basis of principal directions of the covariance matrix C. Let C = U ΛU T denote the spectral decomposition of C, where U is the orthogonal matrix of eigenvectors and Λ = Diag(λ 1 , λ 2 , . . . , λ d ) is a diagonal matrix with the eigenvalues sorted in non-increasing order. The eigenvalue λ i is the variance of the a priori distribution along the ith principal direction. Then we have that

U T ( p -µ) = Λ(Λ + σ 2 I d ) -1 U T (q -µ).
The diagonal operator S = Λ(Λ + σ 2 I d ) -1 has the coefficients

s ii = λ i λ i + σ 2 = (1 + snr -1 i ) -1 (2)
which apply a linear shrinkage of the principal components depending on the signal-to-noise ratio snr i = λ i /σ 2 of each component.

The LMMSE estimator requires the knowledge of the mean and covariance of the a priori distribution, which have to be estimated from the data. The resulting data-dependent filter is sometimes referred to as empirical Wiener filter.

The mean and covariance from the a posteriori distribution can be estimated from the data. Since the noise has zero mean and is independent from u, we have that

E{q} = E{p} = µ, E{(q -µ)(q -µ) T } = C q = C + σ 2 I d .
These can be estimated as the sample mean and sample covariance matrix

µ = 1 n n i=1 q i and C q = 1 n n i=1 (q i -µ)(q i -µ) T . ( 3
)
These correspond to the maximum likelihood estimators (MLE) if the a priori is Gaussian (in which case the a posteriori distribution is Gaussian as well, since the noise is also Gaussian). 2 The above method is used in [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] for still image denoising, where the a priori covariance matrix is estimated as C = C q -σ 2 I d . This may lead to a matrix which is not positive semi-definite (and therefore not a valid covariance matrix). We will address this issue later in Section 3 where we propose a more principled estimator, valid in cases in which the space dimensionality is high. Fig. 2: For the two reference patches displayed in Fig. 1 we show the 10 nearest neighbors (top four rows), the corresponding MAP estimates (middle four rows) and the sample mean and first 9 principal directions (bottom four rows). The four rows in each case correspond to the four frames of the 9 × 9 × 4 patches displayed vertically.

Figures 1 and2 show two examples of sets of similar 9 × 9 × 4 patches. Figure 1 displays the positions of the similar 200 patches found in a local search region and Figure 2 shows the 10 nearest neighbors in more detail. Each patch is displayed by vertically stacking its four 9×9 temporal slices. For the two examples the figure shows the noisy patches, the result of the LMMSE estimation 2 In [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] the unbiased estimator for the covariance matrix is used, which differs from the MLE in that the sum is divided by n -1 instead of n. In practice, except for small values of n, the choice between these two estimators has little effect on the results. Here we prefer to use the MLE because it is the one studied by Paul in [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF]. Later, in §3.2 we will exploit Paul's results to derive estimators for the eigenvalues of the covariance matrix.

and the mean patch together with the leading 9 principal directions. The Gaussian models were learnt from the n = 200 similar patches shown in Figure 1. The temporal slices of the principal directions show oscillatory patterns, with increasing frequency, oriented in the direction of the dominant edges in the patches. In some cases these patterns are localized in a smaller region. Note that even if the temporal slices of the patch may seem unrelated (particularly for the example on the left), the model captures the correlations between them, and in this sense it differs from a decoupled model in which the frames are independent. Taking into consideration these temporal correlations is crucial for the temporal consistency of the result.

Description of the algorithm

The pseudo-code of the proposed video denoising method is described in Algorithms 1 and 2. As in [START_REF] Dabov | Video denoising by sparse 3D transform-domain collaborative filtering[END_REF][START_REF] Dabov | Image denoising by sparse 3d transform-domain collaborative filtering[END_REF][START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Maggioni | Video denoising, deblocking, and enhancement through separable 4-D nonlocal spatiotemporal transforms[END_REF] we perform two stages. The first stage computes a basic estimate u (1) . The second stage computes the final estimate u (2) using the basic estimate as an oracle: patches from the basic estimate are used to compute the patch distances and to estimate the mean and covariance matrix of the a priori distribution. For the latter, we consider the residue of the basic estimate as white Gaussian noise with variance σ (1) . This of course is not true, but it is a rough approximation useful to set the parameters of the 2nd step. A similar approximation is done in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] also to control a parameter in an iterative denoising method.

Note that all the patches of each group are denoised, and these estimates are aggregated via averaging to form the denoised image. There are two aggregation levels. On the one hand, the value of a pixel in the denoised image is computed as the average of all denoised patches that overlap the pixel. In addition, if a patch belongs to several groups, it will be denoised multiple times and these estimates are averaged as well. The aggregation is a key element of patch-based denoising methods [START_REF] Dabov | Image denoising with block-matching and 3D filtering[END_REF][START_REF] Lebrun | Secrets of image denoising cuisine[END_REF].

Algorithm 1 Video NL-Bayes

Require: Noisy video v, noise σ Ensure: Final estimate of noiseless video u (2) 1: (1) , σ (1) )

u (0) = v, σ (0) = σ 2: u (1) = video_nlbayes_step(v, σ, u (0) , σ (0) ) 3: σ (1) = estimate_noise( u (1) ) 4: u (2) = video_nlbayes_step(v, σ, u

Implementation details

Handling of color. We express the video in a luminance-chrominance color space, by applying the opponent color transform (see [START_REF] Dabov | Image denoising by sparse 3d transform-domain collaborative filtering[END_REF]). In the first step, patch distances are computed using the luminance only, in the second step the three channels are used. Using the n most similar patches, a group is built for each channel. These groups are filtered independently (a Gaussian model is learnt for each of them).

This differs from [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF], where for the second step a joint a priori Gaussian model is used for the color patches. In principle this should allow for a better treatment of color by capturing the Algorithm 2 Video NL-Bayes step Require: Noisy video v, noise σ, basic estimate u (1) , basic noise σ (1) Ensure: Final estimate of noiseless video u (2) 1: Set P = {q : q patch of v} // eligible reference patches 2: while P = ∅ do 3:

Get a patch q from P // reference patch 4:

Retrieve the n nearest neighbors to q in a spatio-temporal volume around q. The distance is computed between the basic estimates p (1) .

5:

Estimate µ (2) and C (2) from the basic estimates p (1) i (see §3).

6:

for all n neighbors q i of q do 7:

Compute LMMSE estimate p (2)
i , Eq. ( 1).

8:

Aggregate p

(2) i on u (2) .

9:

Remove q i from P.

// for speed-up 10:

end for 11: end while correlations between the color channels. We found that this option leads to slightly better results (an improvement of 0.2dB up to 0.4dB), but at a much higher computational cost. This is due to the fact that if the channels are treated independently the patch dimensionality is s 2

x s t , whereas a joint treatment requires a single Gaussian model of tripled dimensionality. As we show in the next section, the processing of a group of patches scales with the patch dimensionality d as O(d 2 ) .

Patch distance threshold. When building the group of similar patches during the second stage, more than n patches are allowed if their distances with respect to the reference patch are smaller than a threshold ρ. As in [START_REF] Lebrun | Implementation of the "Non-Local Bayes" (NL-Bayes) Image Denoising Algorithm[END_REF], we set ρ = 4. The distance between patches is normalized by the number of elements in the patch, the number of channels, and the values of each channel are in the range [0, 255]. 3 Thus ρ can be interpreted as the maximum average pixel difference between two similar patches. It is a small value to guarantee that only very similar patches are considered.

In the same way the distance threshold ρ is used to include a larger number of similar patches, some methods have mechanisms to exclude patches if they are too far away from the reference patch. In [START_REF] Niknejad | Image restoration using Gaussian Mixture Models with spatially constrained patch clustering[END_REF] a Gaussian is estimated using a fixed number of nearest neighbors. Then a weighted aggregation is performed which gives small weights to the patches that are not likely according to the learnt model. In this work we did not explore such options.

Sample mean in the second step. In the second step if the sample mean µ (2) is computed using the patches from the basic estimate, the result tends to be too smooth. Textures and details are better preserved if the noisy patches are used instead. In doing so, however, some noise is kept in flat areas, particularly if the number of similar patches n is small and the noise is high. To balance these two phenomena we perform a simple to determine if a region is flat. Specifically, we compute the sample variance of all pixels of all the patches in the group (we use for this the noisy patches).

If this variance is lower than σ 2 , then the area is considered flat. This test is essentially a χ 2 test for the variance. For flat patches we compute the sample mean using the basic estimate, whereas for non-flat ones we use the noisy patches. In [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] flat patches are also given a special treatment, and are estimated as the average of all the pixel values of all patches in the set.

Computational complexity and speed-ups

Processing each group of patches requires (1) searching for the nearest neighbors, (2) estimating the a priori covariance matrix and (3) computing the LMMSE estimates. The search for the nearest neighbors is of asymptotic order O(w 2

x w t d). To estimate the covariance matrix and the LMMSE estimates it is convenient to work on the basis of principal directions. As we will see in Section 3 it is usually enough to consider only the r < d leading principal directions. This requires nd 2 operations to compute the sample covariance matrix, rd 2 to compute its r leading eigenvectors and eigenvalues, and nrd operations to apply the Wiener filter on the r principal components. Therefore the computational complexity of processing each group of similar patches is of asymptotic order O(w 2

x w t d + nd 2 + rd 2 + nrd). The most expensive operations are the computation of the covariance matrix and of its principal directions. There are methods for approximating the r principal directions of a set of points that do not require computing the covariance matrix [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. Instead these methods efficiently approximate a truncated SVD of the data matrix, and would allow to improve the computational complexity to O(w 2

x w t d + drn). The computation can be substantially accelerated by reducing the number of groups of patches that are jointly processed. Each group of patches is associated to the reference patch from which the nearest neighbors are computed. We apply two tricks considered in [START_REF] Dabov | Image denoising by sparse 3d transform-domain collaborative filtering[END_REF][START_REF] Lebrun | Implementation of the "Non-Local Bayes" (NL-Bayes) Image Denoising Algorithm[END_REF] in order to reduce the number of reference patches. The first one consists in visiting only those patches on a coarser spatial grid with vertical and horizontal separation of s x /2. This reduces the number of reference patches by a factor of s 2

x /4. The second acceleration trick is to give up taking as reference patch the patches that have been already once filtered in another group (step 9 in Algorithm 2). The actual number of aggregated estimates for each pixel remains anyway high.

Due to these speed-ups, the actual number of groups of similar patches processed depends on the patch size s x , s t and n. In practice, for the found optimal parameters (see §4 it ranges roughly from 1% to 5% of the total number of pixels in the sequence.

The running time depends on the patch size. The following table shows the running times of our non-optimized C++ implementation for different patch sizes. The patch size is indicated as s 2

x × s t and the time is expressed in seconds per 10 6 pixels (s/MPx): patch size These times were computed by averaging the running times obtained for 5 color test sequences for several noise levels. For comparison we also show the running time of V-BM4D [START_REF] Maggioni | Video denoising, deblocking, and enhancement through separable 4-D nonlocal spatiotemporal transforms[END_REF] using the implementation available online [START_REF] Maggioni | Implementation of V-BM4D[END_REF] with the "modified profile", which is the one with the highest complexity.

3 2 × 2 4 2 × 3 7 2 × 2 10 2 × 1 V-BM4D-

Related works

Using a Gaussian patch model for image denoising has several antecedents. The BLS-GSM method [START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF] (Bayes Least Square estimate of Gaussian Scale Mixture) models noiseless "wavelet coefficient neighborhoods" with a Gaussian scale mixture defined as a random scaling of a zero-mean Gaussian density. The wavelet coefficient neighborhood turns out to be a patch of an oriented channel of the image at a given scale. In [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF], the authors introduce a framework for solving inverse problems, based on the assumption that the patches of an image are distributed according to a Gaussian Mixture Model (GMM) which is learnt from the image with an EM-like algorithm. A related approach is proposed in [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF], where the patches of the image are clustered according to their "structural similarity" and a Gaussian model is fit to each cluster. A GMM is also used in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], but it is trained from a database of 2 • 10 6 patches randomly sampled from the Berkeley database.

A hierarchical covariance data structure termed the covariance tree was introduced in [START_REF] Guillemot | Covariance trees for 2d and 3d processing[END_REF]. The structure captures the covariance of image patches at several scales. The covariance tree can be learnt from a noisy image in the case of denoising, or from a database of images, in which case it can be applied to other inverse problems.

Instead of learning a GMM (or covariance tree) for the whole image in [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] a Gaussian model is fit to each patch and its nearest neighbors. Aguerrebere et al. [START_REF] Aguerrebere | Inverse Problems in Imaging: a Hyperprior Bayesian Approach[END_REF] extended this approach to an hyper-Bayesian model for similar image patches by considering priors on the mean and the covariance matrix. This allows handling more complex inverse problems, where the mean and covariance matrix cannot be directly estimated from the observed data.

Other related approaches are those that use PCA on groups of patches. In [START_REF] Muresan | Adaptive principal components and image denoising[END_REF] a PCA model is learnt for the patches contained on a window sliding over the image. A version of the BM3D algorithm using shape adaptive patches and PCA as a transformed domain was proposed in [START_REF] Dabov | BM3D Image Denoising with Shape-Adaptive Principal Component Analysis[END_REF]. In [START_REF] Zhang | Two-stage image denoising by principal component analysis with local pixel grouping[END_REF] an algorithm that applies PCA on groups of similar patches is introduced. PCA models of patches are is also studied in [START_REF] Deledalle | Image denoising with patch based PCA: local versus global[END_REF], where the authors compare the performance of different estimators (linear and non-linear) in the domain of the principal components.

Although most of these works focus on white Gaussian noise, PCA models have also been extented to Poisson noise in [START_REF] Salmon | Poisson noise reduction with non-local pca[END_REF]. This requires a specific PCA model suitable for that type of noise. A multiscale version of [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] was proposed in [START_REF] Lebrun | The Noise Clinic: a Blind Image Denoising Algorithm[END_REF] to handle structured noise, focusing on JPEG compression artifacts. The authors in [START_REF] Laus | A nonlocal denoising algorithm for manifold-valued images using second order statistics[END_REF] propose an extension of [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] for the denoising of images taking values on Riemannian manifolds.

Gaussian/PCA models for patches lead to several successful image restoration methods. Yet, their application to video denoising has been limited. In [START_REF] Li | Video denoising using shape-adaptive sparse representation over similar spatiotemporal patches[END_REF], the PCA model of [START_REF] Dabov | BM3D Image Denoising with Shape-Adaptive Principal Component Analysis[END_REF] was adapted to video denoising by considering 3D shape-adaptive patches. Buades et al. [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF] denoise 2D patches by using a PCA model (in a first step) or a Gaussian prior (in a second step) learnt from similar patches. The similar patches are taken as the 2D slices of 3D motion compensated patches.

All these works assume that the patch point density function can be locally approximated with Gaussians (locally in the patch space). Related ideas of local Gaussian approximations can also be found applied to manifold learning, where a manifold of dimensionality r is approximated as a mixture of Gaussians of rank r [START_REF] Kambhatla | Dimension reduction by local principal component analysis[END_REF][START_REF] Tipping | Mixtures of probabilistic principal component analysis[END_REF]. Their success does not necessarily mean that similar patches (or nearby points) do follow Gaussian distributions. Nevertherless, regardless of the exact shape of the distribution of a group of similar patches, the best linear estimator in the MSE sense is still the Wiener filter.

This being said, other models for local patches have been proposed recently. For example the authors in [START_REF] Gu | Weighted nuclear norm minimization with application to image denoising[END_REF] propose to penalize a weighted nuclear norm of the matrices formed by stacking a fixed number of similar patches. The resulting estimator is similar to the Wiener filter, in the sense that it is a shrinkage operator on the principal directions.

Empirical Wiener filters in high dimensions

Wiener filters have been used extensively for image and video restoration, and denoising in particular. Here we consider the case of a random vector x ∈ R d contaminated with additive white Gaussian noise y = x + z with z ∼ N (0, σ 2 I d ). The signal of interest x can be a patch of an image or the image itself. In practice, the mean µ and covariance matrix C x of x are unknown and have to be estimated as well. Suppose a set of n noisy observations y 1 , . . . , y n is available.

In this section we discuss principled ways to estimate the eivengalues of the covariance matrix C x , based on results by Paul [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF] about the asymptotic behaviour of the spectral decomposition of the sample covariance matrix when both d and n tend to infinity at the same rate, for Gaussian distributions with spiked covariances. These results fit well the video denoising setting due to the high dimensionality d of spatio-temporal patches.

The resulting asymptotic estimator defines a threshold below which the a priori eigenvalue and the corresponding eigenvector cannot be recovered from the noise. In the asymptotic regime, the value of the threshold can be computed from the noise and the ratio between dimensionality and the number of samples. However, for a finite sample size, the resulting threshold turns out to be too conservative. This will motivate the use of a simple thresholded estimator, where the threshold now is manual parameter which is set during a training stage.

Empirical Wiener filters for denoising

Often in the denoising literature, there is a single observation available (n = 1), which makes the estimation of the a priori parameters µ and C x highly ill-posed. The available observation can be the whole noisy image (as in some Bayesian wavelet shrinkage methods [START_REF] Ghael | Improved wavelet denoising via empirical wiener filtering[END_REF]), a noisy patch (as in block DCT methods [START_REF] Oktem | Image filtering based on discrete cosine transform[END_REF][START_REF] Yu | DCT image denoising: a simple and effective image denoising algorithm[END_REF]), or a 3D stack of noisy patches, as in BM3D [START_REF] Dabov | Image denoising by sparse 3d transform-domain collaborative filtering[END_REF]. These works restrict µ and C x , by assuming that µ = 0 and that the eigenvectors u 1 , . . . , u d of C x are given by a fixed orthonormal basis (such as wavelets, or DCT). In doing so, only the variances λ i have to be estimated. A common strategy is to estimate them from the unique observation available as λ i = max{ y, u i 2 -σ 2 , 0}. In addition, some works use a two-step approach. The first step consists in obtaining a rough estimate of x with non-Bayesian methods, such as universal wavelet thresholding. This so-called pilot or basic estimate is used in the second step to learn the parameters of the a priori distribution [START_REF] Dabov | Image denoising by sparse 3d transform-domain collaborative filtering[END_REF][START_REF] Ghael | Improved wavelet denoising via empirical wiener filtering[END_REF].

More interesting is the case in which several samples are available. Some works fit a Gaussian model (or a mixture of Gaussian models) to a set of n patches [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF][START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. The estimation problem is now "better posed", and estimates of the mean and the covariance matrix (both eigenvectors and eigenvalues) can be provided, at least when n is sufficiently large compared to the dimension d. In our setting, similar to [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF], we estimate the mean and covariance from n patches given by a reference patch and its n -1 nearest neighbors.

Since C y = C x + σ 2 I d the authors in [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] proposed to estimate C x as C x = C y -σ 2 I d , where C y is the sample covariance matrix of the set of observed noisy patches. While the eigenvalues of C y are larger or equal than σ 2 , the estimation errors on C y due to the finite sample size cause small eigenvalues to appear below σ 2 .

Denoting by ξ i the eigenvalues of C q , the eigenvalues of C are λ i = ξ i -σ 2 . The resulting filter coefficient s ii is given by

s ii = λ i λ i + σ 2 = 1 -σ 2 / ξ i ,
and can be negative and of arbitrarily large magnitude as ξ i approaches 0. These large negative weights will amplify the components on the directions of smaller variance, which mostly capture noise, resulting in a poor estimate. To avoid this, it seems reasonable to clip the negative a priori variances as λ i = ( ξ i -σ 2 ) + . Such estimate of the covariance matrix has been also used on the context of patch-based denoising in [START_REF] Chatterjee | Patch-based near-optimal image denoising[END_REF].

Under the assumption than the a priori distribution is Gaussian a natural estimator is the maximum likelihood (ML) estimator. Assuming that the mean is known, the ML estimator of the inverse covariance matrix Q y = C -1 y results from the following convex semidefinite program (for a derivation see [START_REF] Boyd | Convex Optimization[END_REF], §7.1.1):

max log det Q y -tr(Q y C y ) subject to 0 ≺ Q y σ -2 I d .
The inequalities in the constraint force the eigenvalues of S y to be on the interval (0, σ -2 ]. The ML estimator for C x , results from inverting the solution of the above SDP and subtracting σ 2 I d , and is guaranteed to be positive semidefinite. Although there are efficient algorithms for solving such an SDP, they are still prohibitive for the present application.

The covariance matrix estimation has of course many applications other than image restoration, and there is a vast literature in statistics about the topic. Most methods focus on classes of matrices with some sparsity constraint, either on the covariance matrix itself or on its inverse, the precision matrix. See [START_REF] Cai | Estimating structured high-dimensional covariance and precision matrices: Optimal rates and adaptive estimation[END_REF] for a recent survey. Among these matrix classes, the one that fits more naturally our local patch models is the spiked covariance, which expresses the covariance matrix as a low-rank signal component plus noise (AGWN). Several authors [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF][START_REF] Johnstone | On consistency and sparsity for principal components analysis in high dimensions [with comments[END_REF][START_REF] Nadler | Finite sample approximation results for principal component analysis: A matrix perturbation approach[END_REF] studied the relation between the sample covariance matrix and the population one, showing that the sample eigenvectors and eigenvalues are inconsistent estimators. Cai et al. [START_REF] Cai | Optimal estimation and rank detection for sparse spiked covariance matrices[END_REF], impose additionally a group-sparsity constraint on the signal eigenvectors, and proposed estimators with an optimal minimax risk convergence rate. The proposed estimators are mainly of theoretical interest, since require a global search for the support of the eigenvectors.

Here we focus on the simpler problem of estimating the eigenvalues, and propose two estimators derived from the asymptotic results of [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF] that leads to improved reconstruction results in comparison to other estimators used previously in the denoising literature, with the same computational cost.

In the remainder of the section we focus our attention on the estimation of the covariance matrix and assume that the mean µ is zero and it is known or can be well estimated.

Convergence of the sample covariance matrix

We start by reviewing the results of Paul [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF] on the asymptotic behaviour of the spectral decomposition of the sample covariance matrix for spiked covariance models. We consider that samples x 1 , . . . , x n are Gaussian random variables with 0 mean and covariance

C x = U Diag(λ 1 , . . . , λ m , 0, . . . , 0) U T ,
with the eigenvalues sorted in non-increasing order. The noisy observations y 1 , . . . , y n are contaminated with white Gaussian noise of standard deviation σ 2 . The following result can be obtained from Theorems 1 and 2 in [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF] via a change of variables. Theorem 1 (Almost sure limit of the sample eigenvalues) Suppose that d/n → γ ∈ (0, 1) as n → ∞. Let ξ i be the ith largest eigenvector of the sample covariance matrix C y , with 1 i m. Then ξ i converges almost surely to

ξ i →    σ 2 (1 + √ γ) 2 if λ i √ γσ 2 , f (λ i ) := (λ i + σ 2 ) 1 + γσ 2 λ i if λ i > √ γσ 2 .
Furthermore, if λ i has multiplicity one, the scalar product between the corresponding sample eigenvector u i and the population eigenvector u i converges almost surely to

u i , u i 2 →    0 if λ i √ γσ 2 , 1 -γσ 4 /λ 2 i 1 + γσ 2 /λ i if λ i > √ γσ 2 .
The interesting fact here is that there is a threshold √ γσ 2 below which the eigenvalueeigenvector pair cannot be recovered from the noise. The squared root of the ratio γ between the dimensionality and the number of samples acts as a noise amplifier in the threshold. As intuition suggests, the lower the number of samples in relation to the dimensionality, the higher the threshold for recoverability. Additionally, it is shown in [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF] that the recoverable sample eigenvalues and eigenvectors are asymptotically normal.

Estimators for the eigenvalues of C x

We can derive an estimator for λ i by assuming that the asymptotic regime holds. The population eigenvalue λ i and the corresponding eigenvector are recoverable if and only if ξ > σ

2 (1 + √ γ) 2 .
This gives us a recoverability threshold on the sample eigenvalue. If a sample eigenvalue is below that threshold, it is reasonable to set the population eigenvalue to zero, since the corresponding eigenvector is essentially noise. The function f is invertible when restricted to [ √ γσ 2 , ∞). Therefore, if ξ i is above the recoverability threshold, we can estimate it as λ i = f -1 ( ξ i ). We thus propose the following estimator λ S :

λ S ( ξ) = 0 if ξ σ 2 (1 + √ γ) 2 , f -1 ( ξ i ) if ξ > σ 2 (1 + √ γ) 2 .
where

f -1 ( ξ i ) = 1 2 1 + 1 - 4γσ 4 ( ξ -σ 2 (1 + γ)) 2 .
Figure 3 shows a plot of λ S (ξ) for σ = 5 and γ = 1. When ξ σ 2 (1 + √ γ) 2 , the estimator approaches ξ -σ 2 (1 + γ). As the estimator approximates a linear function of ξ, the bias reduces, since

E{ λ S ( ξ)} ≈ E{f -1 ( ξ i )} ≈ f -1 (E{ ξ i }) = f -1 (f (λ i )) = λ i .
The first approximation is because when λ i is large enough, the probability that ξ falls below the recoverability threshold is negligible. And the second approximation is due to the approximate linearity of the estimator. The proposed estimator is valid if the asymptotic regime holds, but we have no optimality guarantees for a finite size sample. Nadler [START_REF] Nadler | Finite sample approximation results for principal component analysis: A matrix perturbation approach[END_REF] provides finite sample size bounds but only for the case m = 1. In practice for our video denoising application, we use patches with dimensionality between d = 100 and d = 200 (for example for a patch of 5 × 5 × 4 or 7 × 7 × 4). Typical sample sizes are of the same order.

In Figure 4 we show results of experiments with simulated data. We consider m = 20 eigenvalues λ 1 = 20, λ 2 = 19, . . . , λ m = 1, and noise with σ = 2. We show the estimated eigenvalues (averaged over 400 simulations) for different values of the dimensionality d and the ratio γ. As expected, both estimators perform better for smaller γ, since the number of samples is higher. The

= 5. If ξ σ 2 (1 + √ γ) 2 , then λ S approximates ξ -σ 2 (1 + γ).
estimator ( ξ -σ 2 ) + has a positive bias. When γ increases the bias increases for the first eigenvalues. Although it drops faster to zero, it still largely overestimates many zero eigenvalues, resulting in a filter which allows many noise components to pass. When d increases, the bias on the m first eigenvectors reduces, but on the other hand, a larger number of zero eigenvalues are estimated as positive. The estimator λ S is much more accurate. The results in the first row confirm that when the asymptotic regime applies, the estimator behaves well. For larger values of γ and lower dimensionality, the estimation accuracy degrades, and many positive eigenvalues are estimated as close to zero, even if they are above the recoverability threshold.

In the light of this, we propose an additional estimator by hard thresholding the difference ξ -σ 2 . The value of the threshold is given by a parameter τ :

λH i = H τ ( ξi -σ 2 ) = ξ -σ 2 if ξ τ σ 2 0 if ξ < τ σ 2 .
For τ = 0 we recover ( ξ -σ 2 ) + . In this way, varying τ we can balance between overestimating zero eigenvalues and underestimating signal components. Determining the optimal value for τ for a finite sample size is a difficult task. Instead, we shall include it as an additional parameter of the denoising algorithm and tune it over a training dataset. Our final objective is not estimating the eigenvalues of the a priori covariance matrix, but the clean data vectors from the noisy ones. Therefore, we evaluate the performance of each estimator by measuring the denoising MSE obtained by the resulting empirical Wiener filter. Figure 5 shows the denoising MSE obtained by different estimators of the spectrum of C x : ( ξ -σ 2 ) + , λ S and λH τ , varying the value of τ > 0. The MSE is estimated as an average of the denoising errors obtained for 400 realizations of the noisy dataset. The obtained MSE is normalized by the minimum MSE and is plotted as a function of τ . Each plot corresponds to a different value of the dimensionality d and the ratio γ as in Figure 4. Additionally, we also show the MSE of an oracular empirical Wiener filter built from the true population variances λ but using the eigenvectors u i of the sample covariance matrix.

The best MSE attained by the empirical Wiener filters depends mainly on the ratio γ. As expected, when n decreases in relation with the dimensionality, the denoising MSE increases. In most cases, the eigenvalues λ S yield the best results, close to the ones obtained with the oracular Fig. 4: Estimated variances when varying the dimensionality d and the ratio γ. Each plot shows the true variances, and the ones estimated from the sample covariance matrix of the noisy data using two estimation methods. The horizontal axis corresponds to the eigenvalue index i = 1, . . . , 60, and the vertical axis to the value of the estimated eigenvalue using the different estimators. From left to right, results with d = 100, 200 and 400. From top to bottom, γ = 0.4, 0.8, 1.6 and 3.2. Note that diagonals have the same number of samples n = d/γ. λ, and even slightly better for d = 400 and γ = 0.2. 4 As expected its denoising performance drops when d is small and γ large. In these situations, better results are obtained with ( ξ -σ 2 ) + . The results obtain with the estimator λ H depend on τ , but the optimal value is always close to the result attained by the best among the other filters.

Parameter selection

In view of the results of the previous section, we shall consider two variants of the video NL-Bayes algorithm described in Algorithms 1 and 2 which differ on the estimator used for the eigenvalues of the a priori covariance matrix. We denote by VNLB-S the version of the algorithm that uses the aymptotic estimator λ S , and by VNLB-H the one corresponding to the thresholded estimator λ H .

For each stage of the algorithm we need to specify the following parameters: the spatial and temporal components of the patch size s x and s t , the spatial and temporal components of the search region w x and w t , and the number of similar patches n. In addition, for VNLB-H we also need to specify the value of τ . When needed, we will add a subscript indicating the stage, e.g. n 1 and n 2 . We consider different parameter settings for color and grayscale videos. 5The parameters of the search region are chosen manually. The rest of the parameters are determined by sampling the parameter space and picking the ones which maximize the average PSNR over a training set of sequences. We consider a grayscale training set containing six sequences of twenty frames, and a color training set of four sequences of ten frames. We selected training sequences having different characteristics on their dominant motion (static, translational, nonrigid) and on their image content (random textures or geometric structures).

Search region

Due to the temporal consistency of natural image sequences, the patches similar to the reference patch are expected to be located around its motion trajectory in neighboring frames. Thus we define the spatio-temporal search region as a sequence of w t square windows of size w x ×w x , whose centers follow the motion trajectory of the reference patch. To estimate the trajectory we compute the forward and backward optical flows of the sequence. The forward half of the trajectory ϕ x,t passing throught (x, t) is computed by integrating the forward optical flow v f as follows:

ϕ x,t (h) = v f ([ϕ x,t (h -1)], h -1) + ϕ x,t (h -1), h = t + 1, . . . , t + w t /2 ,
where [ • ] and • denote the round and floor operators. The backward half of the trajectory is defined analogously using the backward optical flow v b .

In all our experiments we set w x = 27, and w t = 13, yielding a search volume of 9477 patches. Due to the size of the spatial search window, the optical flow does not need to be accurate. We have tested with two optical flow algorithms: the TV-L1 method [START_REF] Zach | A Duality Based Approach for Realtime TV-L1 Optical Flow[END_REF] (we used the implementation of [START_REF] Sánchez | TV-L1 Optical Flow Estimation[END_REF]), and its robust variant described [START_REF] Monzón | Robust Discontinuity Preserving Optical Flow Methods[END_REF], and have found virtually no difference in the denoising result. There is also no significant difference between computing the optical flow on the noisy data, or using an oracular optical flow computed on the ground truth. For our experiments, we use the TV-L1 optical flow since it is fast to compute.

= λ λ = ξ -σ 2 λ = λS ( ξ) λ = H τ (ξ -σ 2 ) 0 1 2 3 4 5 6 7 8 1 1.5 2 2.5 3 3.5 4 τ MSE/MMSE λ = λ λ = ξ -σ 2 λ = λS ( ξ) λ = H τ (ξ -σ 2 ) 0 1 2 3 4 5 6 7 8 1 1.5 2 2.5 3 3.5 4 τ MSE/MMSE λ = λ λ = ξ -σ 2 λ = λS ( ξ) λ = H τ (ξ -σ 2 ) 0 1 2 3 4 5 6 7 8 1 1.5 2 2.5 3 3.5 4 τ MSE/MMSE λ = λ λ = ξ -σ 2 λ = λS ( ξ) λ = H τ (ξ -σ 2 ) 0 1 2 3 4 5 6 7 8 1 1.5 2 2.5 3 3.5 4 τ MSE/MMSE λ = λ λ = ξ -σ 2 λ = λS ( ξ) λ = H τ (ξ -σ 2 ) 0 1 2 3 4 5 6 7 8 1 1.5 2 2.5 3 3.5 4 τ MSE/MMSE λ = λ λ = ξ -σ 2 λ = λS ( ξ) λ = H τ (ξ -σ 2 )
= λ λ = ξ -σ 2 λ = λS ( ξ) λ = H τ (ξ -σ 2 )
Fig. 5: MSE obtained by estimating the n data samples from their noisy observations using empirical Wiener filters with different estimators for the a priori variances λ. The plots correspond to ones in Figure 4. The MSE is normalized by the minimum MSE (attained by the Wiener filter) and it is shown as a function of the parameter τ .

Parameters of Bayesian estimation

The empirical Bayesian estimation is controlled by the number of patches n and the parameters of the variance estimator. We determine these parameters as functions of the noise level σ for each of the patch sizes considered.

Based on the results of §3 on simulated data, using a higher number of samples to estimate the a priori covariance matrix yields a lower MSE. In the context of our denoising algorithm, as we increase the number n of similar patches, we include patches that are not as similar to the reference patch and the assumption that the selected patches are IID samples becomes less valid. Therefore there is trade-off between the error in the sample covariance matrix and the accuracy of the model.

The thresholded eigenvalue estimator used in the VNLB-H variant requires the specification of the additional parameter τ . For a fixed patch size, there is an interplay between n and τ . This can be intuitively understood in the light of the asymptotic theory, where the recoverability threshold depends on the ratio γ between the number of samples and the dimensionality. There is also another interplay between the parameters of both stages n 1 , τ 1 and n 2 , τ 2 . Due to this coupling we sample this four dimensional parameter space and search exhaustively for the parameters that optimize the mean PSNR on the training set, for different values of σ. To simplify (and avoid over-fitting) we restrict the relationship between the parameters and σ to be linear.

The asymptotic eigenvalue estimator λ S used in the VNLB-S variant does not require the additional parameter τ , but depends on the noise level σ. This estimator is based on Paul's asymptotic theory [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF] which applies only for low-rank covariance matrices when the observations are contaminated with AGWN. This suits well to the situation in the first step of the algorithm, where the covariance is estimated using the noisy patches. This is not the case in the second step of the algorithm, where patches of the basic estimate are used. The application of Paul's theory in this case is rather forced, since the residue in the basic estimate is not AGWN. In spite of this, a soft-thresholding of the eigenvalues in our estimator seems reasonable under the assumption that the similar patches should lie close to a lower dimensional space. In this case we can think of λ S as a parametric estimator, such as a λ H τ . In the first step, since the noise AGWN, the theory in [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF] gives the appropriate value of the parameter as the noise level σ. For the second step, we found empirically that an effective way to set the noise parameter is to assign an artificial noise level to the basic estimate, denoted as σ (1) in Algorithms 1 and 2. We estimate this artificial noise level with the noise estimation algorithm [START_REF] Colom | Analysis and Extension of the Percentile Method, Estimating a Noise Curve from a Single Image[END_REF]. 6 The remaining parameters n 1 and n 2 , are determined by a 2D parameter search on the training set. As before, we allow for a linear dependency with the noise level σ.

Patch size

The patch size is a key parameter of patch-based denoising algorithms [START_REF] Buades | A review of image denoising algorithms , with a new one[END_REF][START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF]. Increasing the patch size reduces the variance on the patch distance and therefore on the set of similar patches. However, larger patches are less sensitive to finer scale details and less similar among them. We solved this trade-off empirically, by testing different patch sizes on a set of training sequences.

We consider patches sizes with temporal extent ranging from 1 to 4 frames, and dimensionality d ≈ 100 and d ≈ 200: 1: Parameters used for all results shown in §5, for the two variants of the method and two patch sizes. The parameters w x and w t are the same for both steps. For RGB sequences we only consider the VNLB-S method, since its parameter space is smaller.

d ≈ 100 : 10 × 10 × 1 7 × 7 × 2 6 × 6 × 3 5 × 5 × 4, d ≈ 200 : 14 × 14 × 1 10 × 10 × 2 8 × 8 × 3 7 × 7 × 4.
In Figure 6 we plot for each patch size the best attainable PSNR on each sequence of the grayscale training set. The larger patches with d ≈ 200 show better performance, specially as the noise level increases. We can also see that for sequences with stronger motion (such as football and foreman) it is better to use patches with a shorter time span, of two and even one frame. On the other hand on more static sequences, e.g. mobile and salesman, it is better to use a larger time span. For both dimensionalities, patches with two and three frames perform well for the tested noise levels. For low noise levels 2D patches perform well, but their performance decays drastically when σ increases.

Based on these results, in the following we will consider two patch sizes: 7×7×2 and 10×10×2. The corresponding optimal parameters are shown in Table 1.

Results

In this section we present results obtained with the proposed methods on some classical test sequences with white Gaussian noise added. We consider separately the cases of grayscale and color videos. For grayscale sequences we show results obtained with both variants of our method, VNLB-H and VNLB-S, with patch sizes 7 × 7 × 2 and 10 × 10 × 2. These attained the best average performance on the training sequences among the other patch sizes of the same dimensionality. We used the parameters of Table 1. For color sequence we only show the results of VLNB-S. Some sample sequences can be downloaded from the project website7 in uncompressed format. In Table 2 we compare our results against the state-of-the-art method V-BM4D [START_REF] Maggioni | Video denoising, deblocking, and enhancement through separable 4-D nonlocal spatiotemporal transforms[END_REF] with two different choices of parameters. The method labeled V-BM4D-tip corresponds to the results reported in [START_REF] Maggioni | Video denoising, deblocking, and enhancement through separable 4-D nonlocal spatiotemporal transforms[END_REF]. The results labeled V-BM4D-mp were obtained using the implementation available online [START_REF] Maggioni | Implementation of V-BM4D[END_REF]. This implementation provides three parameter profiles of increasing quality and complexity: "low complexity profile", "normal profile" and "modified profile". We show the results obtained with the modified profile.

Each sequence was contaminated with AGWN of σ = 10, 20, 30, 40. Table 2 shows the PSNR values obtained after denoising. We can observe that using a larger patch size gives better results, specially for higher noise levels, but at a higher computational cost. For our methods and V-BM4Dmp we show the average PSNR of the basic estimate as well. This allows to see the improvement achieved by the second iteration.

In Table 3 we also include a comparison with SPTWO [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF].To denoise each frame, SPTWO uses the optical flow to register the neighboring frames (6 past and future frames) to the target frame.

For each 5 × 5 reference patch in the target frame, 13k similar patches are selected by searching for the k most similar 5 × 5 × 13 patch trajectories in the volume defined by the warped frames. Because patches are warped according to the motion, SPTWO is in principle able to better handle zooms, rotations and non-rigid transformations. This comes at the price of a higher dependency on the optical flow. To compare our method with SPTWO, we computed the result of our algorithm on some of the sequences used in [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF]. The grayscale scale sequences have only 30 frames, and the color sequences have 8. 8 Both variants of our algorithm achieve comparable results. The performance of VNLB-H is slighly better (at most 0.3dB on the average). This is not surprinsing given that VNLB-H uses a hand tuned estimator (on both steps of the algorithm) optimized over the test sequences. With two parameters less, the results of VNLB-S are only slightly inferior. This supports the choice of the basic noise level σ (1) to set the noise parameter of the soft-threshold estimator λ S in the second step of the algorithm.

Our results outperform the ones of V-BM4D for all sequences and noise levels, and with the two variants VNLB-S and VNLB-H and patch sizes. On average, the difference with V-BM4D-mp is 1.8dB. The proposed methods also outperform SPTWO on the average (with a margin varying between 1.2dB and 1.6dB depending on the noise level), but the difference is much smaller for bus and tennis. In these sequences the result of SPTWO is slighly better that the ones we obtain with the patch size of 7 × 7 × 2, except for the highest noise levels.

For a qualitative comparison, we show in Figures 7, and8, details of the results obtained with the proposed VNLB-H method, V-BM4D-mp and SPTWO, 9 with noise levels of σ = 10, 40. For VNLB-H we show results obtained with a patch size of 10 × 10 × 2. We omit the results obtained with the VNLB-S since they are visually hard to distinguish.

The differences between the methods become more apparent for higher noise levels. The VNLB-H preserves more details than V-BM4D. This can be noticed for example for mobile and bicycle. SPTWO achieves very good visual results in moving objects for which the motion can be well approximated. This can be seen in the face of salesman. Its results are temporally consistent with the estimated motion. This creates artifacts when the motion is not well estimated. As an example, the wallpaper texture in tennis adopts the motion of the arm. This dragging effect is common close to the boundaries of moving objects, particularly if the background texture is weaker than the noise level. In the bicycle sequence, for noise σ = 40 the rays of the wheel are not well reconstructed and the sticker depicting a train shows deformations. This is probably also due to problems with the motion estimation. Notice also that the result of the proposed method is smoother for constant regions, such as the dark background in bicycle.

The results shown in Figures 7 and8 show the spatial characteristics of the proposed method. As a way of assessing the temporal consistency of the result we show a slice of the video volume at a fixed row (or column). Figure 10 shows a horizontal slice of the mobile sequence. Each row in the images shown corresponds to a different frame. The gray tube corresponds to a ball which is hit by a black train entering from the left. The result of V-BM4D shows a high frequency variation in the temporal (vertical) dimension, a consequence of the flickering on the video sequence. On the other hand, the result of the video non-local Bayes method, shows a much slower variation in time.

To the best of our knowledge, the best denoising results in grayscale sequences have been reported in [START_REF] Li | Video denoising using shape-adaptive sparse representation over similar spatiotemporal patches[END_REF]. Some of these results are surprising compared to the rest of the state of the 8

In [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF] the error is measured as the root mean square error (RMSE) of the central frame. For color sequences the average channel RMSE is used. To allow a direct comparison, we adopt in Table 3 these error measurements, expressing them as a PSNR in decibels. For color sequences, this amounts to PSNR = 20 log 10 255/ 1 3 (RMSE R + RMSE G + RMSE B ) . 9 We are thankful to the authors of [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF] for kindly sharing with us their results.

art. For example, for salesman with σ = 40 the result of [START_REF] Li | Video denoising using shape-adaptive sparse representation over similar spatiotemporal patches[END_REF] is 4dB higher than the best result among all public state-of-the-art methods, including ours. We have not included these results in the comparison table since the code is not available and we have not been able to reproduce the results. This method builds upon BM4D, which in [START_REF] Maggioni | A nonlocal transform-domain filter for volumetric data denoising and reconstruction[END_REF] is applied to volumetric imagery and video, showing an inferior performance than V-BM4D. The authors in [START_REF] Li | Video denoising using shape-adaptive sparse representation over similar spatiotemporal patches[END_REF] consider shape adaptive patches as in [START_REF] Ercole | Spatio-temporal pointwise adaptive denoising of video: 3D non-parametric regression approach[END_REF] (the patches are quite large 8 × 8 × 8), and perform several iterations. In each iteration, the previous iterate is used as an oracle to compute patch distances and the adaptive shape of the patches. Our method can be iterated as in [START_REF] Li | Video denoising using shape-adaptive sparse representation over similar spatiotemporal patches[END_REF], and some improvement can be obtained in a third iteration, but not enough to compensate for the increase computational cost.

Finally in Table 4 we compare the PSNR obtained on the grayscale sequences of Table 2 without using a motion compensated search region (or equivalently, assuming zero motion). The resulting search region is a 27 × 27 × 13 3D rectangle centered at the reference patch. The aim of this experiment is to determine the impact of motion compensation on the proposed method. The results shown correspond to the method VNLB-H with a patch size of 7×7×2. Some sequences are mainly static, or have a slow motion (salesman, mobile, tennis). They show almost no difference in PSNR when no optical flow is used. The sequences garden and bicycle have more motion, yet the drop in PSNR caused by assuming zero motion is small, between 0.1dB and 0.3dB. As expected, the sequence which is most affected is stefan, which has a very fast camera motion. In this case the difference is quite important: ≈ 0.8dB. Nevertheless the result obtained with no motion compensation still outperforms V-BM4D. This shows that the proposed method does not have a critical dependence on the quality of the optical flow, and in fact can still be competitive without any motion estimation.

Results on color videos

We considered five color test sequences, tennis, coastguard, foreman, bus and football. The first four were chosen because they were used in [START_REF] Maggioni | Video denoising, deblocking, and enhancement through separable 4-D nonlocal spatiotemporal transforms[END_REF]. The fifth one, football, was added to test the performance of our method on a sequence with complex motions.

The results are shown in Table 5. As for grayscale videos, we compared our method with V-BM4D using the two parameter profiles: V-BM4D-tip and V-BM4D-mp. We also include results of CIFIC [START_REF] Dai | Color video denoising based on combined interframe and intercolor prediction[END_REF], a recent method. The proposed method, achieves the highest PSNR values on all the considered sequences and levels of noise, with a significant difference. The difference increases for higher levels of noise. A possible reason for this is that both CIFIC and V-BM4D rely on motion estimation by 2D block matching, which might not be reliable for high noise levels.

To compare with SPTWO, we show in the right half of Table 3 our results on four 8 frame sequences used by [START_REF] Buades | Patch-based video denoising with optical flow estimation[END_REF] (see footnote 8). SPTWO achieves excellent results for the sequence army possible due to the fact that its optical flow can be well estimated. On the other sequences, the proposed methods attain a higher PSNR.

Figure 9 shows results obtained by our method and V-BM4D for the sequence bus, for noise of σ = 10, 40. As for grayscale sequences, it can be noticed that some details are better preserved by the proposed method. V-BM4D also shows some color artifacts for high noise values (noticeable for instance on the street bus).

Conclusions

We presented a Bayesian video denoising algorithm assuming that similar spatio-temporal patches are IID samples from an a priori distribution. The proposed method uses motion estimation only Fig. 7: Details of the results obtained for the grayscale versions of mobile and tennis. For each sequence, the left column corresponds to σ = 10 and the right one to σ = 40. At the top we show the original frame, and the second row shows the noisy frames. In the third row we show the result of V-BM4D-mp for mobile and SPTWO for tennis. On the bottom row we show our result (VNLB-H 10 2 × 2). to guide the search for similar patches, thus it does not rely on an accurate motion field. It provides state-of-the-art results, both in terms of PSNR and qualitatively.

The patches in each group are estimated using an empirical Wiener filter, which requires estimates of the mean and the covariance matrix of the a priori distribution. We studied this estimation problem in the light of the asymptotic results of [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF], assuming a low-rank Gaussian prior. This led us to propose two estimators for the eigenvalues of the a priori covariance matrix. One of them is derived by inverting the asymptotic bias of the eigenvalues of the sample covarianace matrix predicted by [START_REF] Debashis | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF]. This results in a soft-thresholding of such eigenvalues. The resulting estimator performs well when the asymptotic regimes applies, but degrades as the number of patches is reduced. This led us to propose a second estimator consisting on a hard threshold, where the threshold value is controlled by a hand-tuned parameter. Two variants of the video denoising algorithm were proposed: VNLB-S and VNLB-H, based respectively on the soft and hard estimators. VNLB-S is considerably easier to tune, since it has one less parameter per step, and achieves results which are only slightly inferior to those of VNLB-H. The significant improvement over the state-of-the-art suggests that the Gaussian (or PCA) models that have been successfully applied to still images processing, can also contribute to video restoration. The high-dimensionality of spatio-temporal patches poses challenges in the estimation of the moments of the prior distribution. In this setting, better Bayesian estimators can be obtained by a theoretically motivated optimization of the empirical Wiener filter, taking into account the uncertainty on the sample covariance matrix.

Another conclusion is that very good results can be obtained without motion estimation. We do not believe that video processing/restoration should be performed in general without taking motion into account, but it is interesting to see how much can be done without compensating motion.

The proposed approach could be extended to tackle more general inverse problems by considering an hyper-Bayesian framework as in [START_REF] Aguerrebere | Inverse Problems in Imaging: a Hyperprior Bayesian Approach[END_REF]. The main obstacle is the computation time. Another logical extension would contemplate including motion deblurring in the Bayesian formulation.

Fig. 1 :

 1 Fig. 1: Examples of sets of 200 nearest neighbors of two reference spatio-temporal patches, of size 9 × 9 × 4. The images in each row show the five frames of a 37 × 37 × 5 search region. In red we show the positions of the 5 nearest neighbors (one of them, in the center of the middle frame is the reference patch), in green the next 40 and in blue the rest. The points shown correspond to the top-left pixel of the first spatial slice of each patch. To highlight the position of the patches, the color of the images has been attenuated.

Fig. 3 :

 3 Fig. 3: Proposed variance estimator λ S as a function of the sample eigenvalue ξ, for γ = 1 and σ = 5. If ξ σ 2 (1 + √ γ) 2 , then λ S approximates ξ -σ 2 (1 + γ).

  λ = ξσ 2 λ = λS ( ξ) λ = H τ (ξσ 2 ) λ = ξσ 2 λ = λS ( ξ) λ = H τ (ξσ 2 ) λ = ξσ 2 λ = λS ( ξ) λ = H τ (ξσ 2 )

  λ = ξσ 2 λ = λS ( ξ) λ = H τ (ξσ 2 ) λ = ξσ 2 λ = λS ( ξ) λ = H τ (ξσ 2 )

Fig. 6 :

 6 Fig. 6: Highest attainable PSNR for each patch size over the grayscale training set. For each patch size and sequence, we plot on the vertical axis the highest PSNR obtained with that patch size for that sequence minus the highest PSNR value attained for the sequence among all considered patch sizes. The rightmost points on each plot correspond to the average PSNR over the training set. The columns correspond to the noise levels σ = 6, 12, 48. The top row shows the results for patches sizes with d ≈ 100 and those with d ≈ 200 are shown in the bottom row.

Fig. 8 :

 8 Fig. 8: Details of the results obtained for the grayscale versions of bicycle and salesman. For each sequence, the left column corresponds to σ = 10 and the right one to σ = 40. From top to bottom we show the (1) original frame, (2) noisy frames, (3) the result of V-BM4D-mp, (4) the result of SPTWO and (5) our result (VNLB-H 10 2 × 2).

Fig. 9 :

 9 Fig. 9: Details of the results obtained for the color version of bus. The left column corresponds to σ = 10 and the right one to σ = 40. From top to bottom we show the (1) original frame, (2) noisy frames, (3) the result of V-BM4D-mp, and (4) our result (VNLB-S 10 2 × 2).

Fig. 10 :

 10 Fig. 10: Visualizing the temporal coherency of the denoising. The images show a horizontal slice of the mobile video, at row 220, showing the ball (gray tube) as it gets hit by the toy train (black region in the left). The vertical axis is time, and the horizontal axis is the x axis in the video. On the first row: the original sequence and noisy version with noise 10, 20 and 40. The following rows show the results obtained with V-BM4D, and the results for VNLB-H with patch size 10 × 10 × 2. The result of VNLB varies smoothly in time, reducing the flickering effect. The contrast of the images has been enhanced for better visualization.

Table 2 :

 2 Quantitative denoising results for some classic grayscale test sequences. The values in italics correspond to the average of the PSNR of the basic estimates. See text for details.5.1 Results on grayscale videosWe present denoising results on grayscale test sequences used commonly in the literature. The sequences have CIF or SIF resolution and between 150 and 300 frames.

		Grayscale sequences			Color sequences
	σ Method	Bus Tennis Sales. Bike	Ave.	Army Cooper Dog Truck	Ave.
	10 SPTWO	36.07 34.69 36.38 36.74	35.97	39.44	36.09	35.90 37.10	37.13
	VNLB-S 7 2 × 2	36.32 34.67 38.13 39.25	37.09	39.30 37.31 36.66 37.29 37.64
	VNLB-S 10 2 × 2 36.68 34.86 38.51 39.22	37.32 39.62 37.31 36.83 37.31 37.77
	VNLB-H 7 2 × 2	36.68 35.04 38.31 39.45 37.37
	VNLB-H 10 2 × 2 36.98 35.18 38.77 39.38 37.58
	20 SPTWO	32.24 30.59 32.95 33.01	32.20 36.47 32.34	32.94 33.65	33.85
	VNLB-S 7 2 × 2	32.31 30.35 34.44 36.09	33.30	36.06	33.44	33.80 33.75 34.26
	VNLB-S 10 2 × 2 32.77 30.70 34.78 36.11	33.59	36.28 33.49 33.92 33.76 34.36
	VNLB-H 7 2 × 2	32.52 30.29 34.61 36.28	33.42
	VNLB-H 10 2 × 2 33.08 30.91 35.06 36.32 33.84
	30 SPTWO	30.05 27.48 30.95 31.62	30.03 34.67 30.24 31.29 31.61	31.95
	VNLB-S 7 2 × 2	29.98 27.69 32.15 33.85	30.92	34.26	31.28	32.13 31.74 32.35
	VNLB-S 10 2 × 2 30.50 27.91 32.57 34.04	31.25	34.52 31.38 32.32 31.79 32.50
	VNLB-H 7 2 × 2	30.11 27.64 32.33 34.05	31.03
	VNLB-H 10 2 × 2 30.76 27.84 32.83 34.30 31.43
	40 SPTWO	28.51 25.53 29.60 29.77	28.35 33.24 28.85	30.15 30.19	30.61
	VNLB-S 7 2 × 2	28.21 26.10 30.40 32.12	29.21	32.90	29.78	30.88 30.35	30.98
	VNLB-S 10 2 × 2 28.76 26.32 30.98 32.42	29.62 33.29 29.94 31.19 30.44 31.22
	VNLB-H 7 2 × 2	28.28 26.05 30.57 32.30	29.30
	VNLB-H 10 2 × 2 28.97 26.22 31.19 32.71 29.77
	50 SPTWO	27.19 24.96 28.34 28.66	27.29	31.99	27.78	29.09 29.13	29.50
	VNLB-S 7 2 × 2	26.95 25.07 29.03 30.59	27.91	31.82	28.70	29.91 29.27	29.92
	VNLB-S 10 2 × 2 27.56 25.20 29.67 31.00	28.36 32.32 28.93 30.30 29.43 30.24
	VNLB-H 7 2 × 2	27.05 25.00 29.20 30.82	28.02
	VNLB-H 10 2 × 2 27.71 25.13 29.86 31.25 28.49

Table 3 :

 3 Comparison with the method SPTWO. Shown PSNRs correspond to the central frame for each test sequences. See text for details.

	σ Method	Tennis Sales. Garden Mobile Bicycle Stefan Average
	10 VNLB-H 7 2 × 2	35.73	38.23	34.47	36.63	39.62	35.78	36.74
	VNLB-H 7 2 × 2 (no m.c.) 35.71	38.23	34.29	36.60	39.51	35.06	36.59
	20 VNLB-H 7 2 × 2	32.00	34.55	30.26	32.56	36.42	31.60	32.90
	VNLB-H 7 2 × 2 (no m.c.) 31.98	34.55	30.01	32.53	36.27	30.74	32.71
	30 VNLB-H 7 2 × 2	30.26	32.27	27.66	29.94	34.15	29.05	30.56
	VNLB-H 7 2 × 2 (no m.c.) 30.23	32.26	27.40	29.91	33.96	28.23	30.36
	40 VNLB-H 7 2 × 2	29.02	30.49	25.69	27.92	32.31	27.19	28.77
	VNLB-H 7 2 × 2 (no m.c.) 28.99	30.50	25.44	27.89	32.10	26.46	28.59

Table 4 :

 4 Results without a motion compensated search region (no m.c.). See text for details.

	σ	Method	Tennis	Coast.	Fore.	Bus	Foot.	Average
	5	V-BM3D	39.45	40.18	40.16	39.07			39.71
		V-BM4D-tip	39.98	41.13	41.38	40.21			40.68
		V-BM4D-mp	39.60	40.28	40.78	39.56	40.00	39.98	40.06
		VNLB-S 7 2 × 2	40.29	41.43	41.89	41.27	41.23	41.13	41.22
		VNLB-S 10 2 × 2	39.92	41.26	42.13	41.05	40.89	41.14	41.09
	10	V-BM3D	36.04	36.82	37.52	34.96			36.34
		V-BM4D-tip	36.42	37.27	37.92	36.23			36.96
		V-BM4D-mp	35.90	36.30	37.21	35.38	36.08	35.91	36.20
		VNLB-S 7 2 × 2	36.70	37.78	38.64	37.65	37.50	37.46	37.69
		VNLB-S 10 2 × 2	36.88	38.01	39.05	37.94	37.62	37.71	37.97
	15	CIFIC	32.76			33.47	31.59
		V-BM4D-mp	33.67	34.05	35.17	33.02	33.82	33.49	33.98
		VNLB-S 7 2 × 2	34.62	35.69	36.76	35.54	35.34	35.25	35.65
		VNLB-S 10 2 × 2	34.85	35.95	37.13	35.91	35.51	35.59	35.96
	20	V-BM3D	32.54	33.39	34.49	31.03			32.86
		V-BM4D-tip	32.88	33.61	34.62	32.27			33.35
		V-BM4D-mp	31.98	32.44	33.70	31.34	32.22	31.76	32.37
		VNLB-S 7 2 × 2	33.17	34.26	35.43	34.03	33.82	33.66	34.22
		VNLB-S 10 2 × 2	33.35	34.50	35.78	34.42	33.99	34.06	34.51
	25	CIFIC	29.74			30.48	28.82
		V-BM4D-mp	30.73	31.24	32.59	30.07	30.96	30.44	31.16
		VNLB-S 7 2 × 2	32.07	33.15	34.35	32.81	32.62	32.42	33.10
		VNLB-S 10 2 × 2	32.18	33.40	34.74	33.24	32.80	32.86	33.39
	40	V-BM3D	29.20	29.99	31.17	27.34			29.43
		V-BM4D-tip	29.52	30.00	31.30	28.32			29.78
		V-BM4D-mp	28.14	28.73	30.09	27.44	28.35	27.66	28.60
		VNLB-S 7 2 × 2	29.78	30.89	32.07	30.28	30.15	29.87	30.76
		VNLB-S 10 2 × 2	29.92	31.19	32.60	30.69	30.36	30.34	31.10

Table 5 :

 5 PSNRs obtained for the classic color test sequences. The averages corresponds to the first four sequences. The values in italics correspond to the average of the PSNR of the basic estimates. See text for details.

We refer to our method as an empirical Bayesian approach because although it is based on a Bayesian model for patches, the parameters of the prior are determined from the data using frequentist estimators. The empirical (or data-driven) Bayesian approach is a usual choice when the parameters of the prior are unknown, but it is a departure from a strict Bayesian methodology which would require a hyperprior to estimate the parameters (see[START_REF] Aguerrebere | Inverse Problems in Imaging: a Hyperprior Bayesian Approach[END_REF] for a closely related method considering a Bayesian estimation of the parameters of the prior).

The distance threshold ρ is applied in the second step, where the distance is computed using patches from the basic estimate. For this reason (as in[START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF]) we do not consider a noise dependent threshold.

Pablo Arias, Jean-Michel Morel

It might seem counterintuitive that the filter computed from the estimated eigenvalues λ S , outperforms the oracular one. It is a consequence of the error in the eigenvectors. In particular, eigenvectors corresponding to eigenvalues below the recoverability threshold are better discarded.

In order to reduce the parameter space, we did not optimize the choice of the distance threshold ρ.

We use the default parameters in[START_REF] Colom | Analysis and Extension of the Percentile Method, Estimating a Noise Curve from a Single Image[END_REF], with the exception of the number of bins, which we set to one since we are not interested in a noise curve, but in just an average noise level.

http://dev.ipol.im/ ~pariasm/video_nlbayes/

This work is partly founded by BPIFrance and Rgion Ile de France, in the framework of the FUI 18 Plein Phare project; by the Office of Naval research by grant N00014-17-1-2552; by ANR-DGA project ANR-12-ASTR-0035; and by ANR-DGA project ANR-14-CE27-001 (MIRIAM).