Denoising images corrupted with Gaussian noise is a classical image processing application, on which tremendous amount of research has been performed. Generalised K-Means clustering (K-SVD) [1] , autoencoders [2] and deep neural networks [3] are worth mentioning works in this field. For the past two decades, non-local patch based methods [4, 5, 6, 7] have been among the state-of-art methods for this application. Especially, 3D Block Matching (BM3D) in combination with some non-linear transformations and inverse transformations is a widely used method for denoising images which are corrupted with a variety of noise distributions including Gaussian, Poisson and Poisson-Gaussian mixture types of noise [8, 9, 10, 11] . All the above mentioned research work, however, has been performed on denoising single images. Denoising multi-frame images is a less explored field. Electron microscopy, CT imaging, multi-spectral imaging are some of the fields where we encounter the problem of obtaining one image from multiple noisy images of the same sample/scene. Notable works in the field of multi-frame image denoising include non-local patch-based denoising methods [12, 13] , lowrank tensor approximation framework with Laplacian scale mixture [14, 15] , Stein's unbiased risk estimator and linear expansion of threshold combining methods [16, 17] , multiscale sparsity denoising of spectral domain optical coherence tomography images [18] , non-local energy functional minimisation approach involving spatio-temporal image patches [19] , low-rank tensor modelling based denoising approach [20] , combining blind source separation and block matching for denoising CT images [21] , and tensor based modeling for denoising multi-frame images corrupted with speckle noise [22] .
Three different types of baseline strategies can exist that are relevant in extending patch-based methods for denoising multi-frame images: First, average the noisy images and then denoise the averaged image [12] . The second is a trivial strategy introduced by us: Denoise the individual images and then average the denoised images. Finally, a non-trivial strategy [13] : In the initial stage of patch-based denoising methods, the group of patches that are most similar to a patch are acquired from images other than just the reference image also. Once this group of patches are acquired, a traditional single image denoising approach is performed.
Our Contribution In this work, our main contribution is a novel strategy that differs from all the above mentioned strategies at three stages. The first and second differences are at the 3D patch grouping stage and the third difference is at the aggregation stage of acquiring the final denoised image. These ideas enhance the denoising quality of the algorithm significantly. We perform a detailed and comprehensive evaluation of all the above mentioned extensions, specifically for BM3D. Such a detailed study of extensions for multi-frame image denoising is missing for any of the patch-based methods in the existing works.
Paper Structure In Section 2, we introduce the various extensions of BM3D for multi-frame images. In Section 3, we present the denoising results of the experiments performed using various extensions of BM3D. We also perform an analysis on why our proposed extension outperforms the other existing approaches. Finally, in Section 4, we conclude the paper with a summary of our work and some ideas for future research. [5] is a two step method for removing additive white Gaussian noise (AWGN). Each step consists of three substeps. In the first sub-step of the first main step, a 3D group is created by finding the most similar patches to a reference patch using L 2 distance. In the second sub-step, a 2D biorthogonal spline wavelet transform is applied to this 3D group, followed by a 1D Walsh-Hadamard transform in the third dimension. After a hard thresholding, the 3D group undergoes corresponding 1D and 2D inverse transformations. In the third sub-step, a weighted aggregation of all the denoised pixels at a particular position in the 2D image domain is performed to obtain the final denoised pixel at that particular position. Thus, we have a denoised image after the first step. In the first sub-step of the second step, the grouping is done using the denoised image obtained after the first step. In the second sub-step, a 2D discrete cosine transform followed by a 1D Walsh-Hadamard transform is applied on both 3D groups obtained from the noisy image and the denoised image from the first step. Then a Wiener filter is applied on a combination of both of the above transformed 3D groups. The resulting 3D group is then back transformed. In the third sub-step, the same corresponding strategy as in the first main step gives the final denoised image. Thus, the two main steps primarily only differ at the second sub-step. We refer the reader to the original work [5] for more details on the method. The final denoised image using the weighted aggregation process [6] is described as
Here, x is a position in the image domain Q and u final is the final denoised image. The symbol u wien Q,P (x) denotes the estimation of the value at pixel position x belonging to the patch Q obtained after the Wiener filtering of the reference patch P . Also, w wien P is obtained from the emperical Wiener coefficients of patch P , P(P ) is the set of most similar patches to the reference patch P , and finally χ Q (x) = 1 if x ∈ Q and 0 otherwise.
The first extension is a trivial existing one [12] : Average all the noisy frames after aligning them and then denoise the averaged image (BM3D-1). Another trivial extension, proposed by us, is to denoise every single image and then average the aligned denoised images (BM3D-2). A third non-trival existing [13] extension (BM3D-3) would be to first consider a reference image. Then while obtaining a 3D group of similar patches in the first step, an L 2 distance evaluation is performed with patches from all the other images also instead of just those from the reference image. The rest of the algorithm remains the same.
Let us now describe the final extension which is the main extension proposed by us: For every image, we carry out a grouping using L 2 distance evaluation from patches in all the other images. Thus, we use the same grouping strategy as the one proposed in BM3D-3, but we do this for every image. Then we perform the same filtering process as the original BM3D method, for the 3D groups in every image as a part of the second sub-step of the first step. In the third sub-step, the same aggregation process as in the original BM3D method is carried out on every image. Thus, we have as many denoised images as we have input images. In the first sub-step of the second step, we perform the same grouping strategy as we have in the first sub-step of the first step, but using the denoised images of the first step. In the second sub-step, we implement the same filtering process used in the original BM3D method on 3D groups in every image. Finally, in the third sub-step, we carry out an aggregation at a particular position in image domain, using a weighted aggregation of the denoised pixels obtained in "all the images" at this particular position.
As mentioned in Section 1, our method differs from the other three strategies in three aspects. First, we do not have a reference image like in BM3D-3, but follow the denoising procedure for all the images. This is done in order to maximise the use of available information. Second, in the detailed study for selection for BM3D parameters [6] and in the original work [5] , two parameters control the number of patches in a 3D group: The threshold parameter for L 2 distance and the maximum number of patches. In this work we do not use the former parameter as we have the risk of losing some similar patches in the case of highly noisy images. Finally, to obtain the final denoised image, all the other existing extensions perform the weighted aggregation of pixels from 3D groups in just one image. However, since we denoise every image unlike the other three strategies, we perform the weighted aggregation from denoised 3D groups in all the images. The equation that governs the weighted aggregation to obtain the final denoised image is modified as follows:
In the above equation, we can see the extra summation over all the images denoted by index . Reference patches P must be considered from every image . The rest of the symbols have the same meaning as described in (1).
Generally, the various noisy images of the same scene in a multi-frame image denoising problem are first registered before denoising them to obtain the final image. For the registration, one can use a motion estimation method. However, in this work, our aim is to carry out a comprehensive performance study of the extensions for denoising methods. Thus, we assume that all the images are pre-registered so that we can straight away compare the denoising qualities of various algorithms. It is not difficult to combine the denoising strategies proposed by us with state-of-art motion estimation methods. Hence, to satisfy this experimental setting, we have corrupted Lena, House, Peppers and Bridge 1 images with AWGN of standard deviations σ noise = 80, 100, 120. Two datasets for each image for each noise standard deviation have been created, one with five realisations of noise and the other with ten realisations of noise. The threshold parameters for L 2 distance while forming the 3D groups have been excluded in both the main steps of BM3D for all the four extensions to keep them on an equal footing. This particular step is advantageous to all the extensions for high standard deviation of noise. One can introduce back these parameters in applications where there is less amount of noise (σ noise < 80). All the other parameters have been chosen according to the detailed study of parameter selection for BM3D in [6] . It also has to be mentioned that, in the results we showcase shortly, we present the mean squared error (MSE) value for BM3D-3 which is the best of the 5/10 denoised images obtained when every image is selected as the reference frame.
From Table 1 and Figures 1 and 2 , we can observe that BM3D-M outperforms all other methods significantly, in terms of both MSE and visually. The images obtained using BM3D-M in particular, appear sharper. Moreover, for the Lena 10-image dataset and σ noise = 80, the next best method other than BM3D-M is worse by 28%. Thus, BM3D-M has the capability of producing significantly better results than the other algorithms. We attribute these results to the crucial differences in the modelling: inter-frame connectivity strategy in both main steps of BM3D to acquire 3D groups and a weighted aggregation that connects all the frames instead of just one frame.
It is not clear if the BM3D-1 or the BM3D-2 method is better. This is because BM3D-1 does not retain the same explicit noise model that was introduced in the original image after averaging the noisy images and BM3D-2 does not have enough signal in each of the input images in each dataset. Also, we would generally expect the non-trivial BM3D-3 method to perform better than the BM3D-1 and BM3D-2 methods. The three algorithms other than BM3D-3 perform some sort of averaging while BM3D-3 does not. It simply does not make full use of the information available in all the images. This explains its surprisingly inferior performance.
Experiments on an Intel(R) Core(TM) i7-6700 CPU @3.4
GHz using OpenMP indicates that BM3D-M takes 11 times more time than single channel BM3D for the 5-image 256 × 256 sized House dataset. Also, an ANSI C code combined with a CUDA implementation of BM3D-M on an NVIDIA Quadro P5000 architecture takes 2.7 seconds.
We have performed the first comprehensive study of extensions for non-local patch-based methods like BM3D for denoising multi-frame image datasets. Our proposed extension which uses novel inter-frame patch and pixel connectivity strategies, gives significantly better results than all the other existing trivial and non-trivial extensions. This improvement of the denoising performance is achieved without using any additional parameters. The new ideas we introduce can also be applied to other similar patch-based methods like BM3D.
In the future, we plan to combine BM3D-M with stateof-art motion registration algorithms and also study the error incurred due to false registration of pixels. We will also study the performance of BM3D-M when the images are corrupted with noise of different standard deviations. Finally, we would apply the algorithm to multi-frame bio-physical image datasets obtained from electron microscopes.