1 Introduction
Multispectral imaging systems have emerged as a new technology that is able to deal with various problems encountered with broadband systems.
The development of new cameras and filters enables us to see beyond the visible spectrum such as in the Infrared [700nm1mm], Ultraviolet [10nm380nm], Xray [0.01nm10nm] spectrums. A multispectral image is a set of several monochrome images of the same scene, each of which is taken at a specific wavelength. Each monochrome image is referred to as a band or channel. A multispectral image may be seen as a three dimensional image cube with two spatial dimensions consisting of the vertical and horizontal axes, and a spectral third dimension. A monochrome image has one spectral band where each pixel is represented by one scalar value. A multispectral image consists of at least two bands. A pixel is represented by a vector of components where is the number of spectral bands. In that sense, a color image may be seen as a multispectral image with three spectral bands. However, the term ”multispectral” is commonly used for images with more than three spectral bands. Images with more than a hundred bands are commonly called hyperspectral images.
Using multispectral images is justified by two main reasons. First, narrow spectral bands exhibit more relevant information compared with conventional broadband color and black and white images. Indeed, we obtain a unique spectral signature of the objects being captured. Such information can be used to enhance the accuracy of image processing applications. Second, by using multispectral images, we are able to separate the illumination information from object reflectance, in contrast to broadband images where it is almost impossible to do so. This separated information can now be used to normalize images. For instance, in face recognition applications, nearinfrared spectral band can be combined with the visible image. This approach has been widely used to construct more effective biometric systems
(Meraoumia et al., 2014; Pan et al., 2003, 2004). Thermal infrared images have also been widely used. Thermal infrared sensors detect the heat energy radiated from the face which is independent from the illumination as in the case of reflectance (Kong et al., 2007, 2005). Furthermore, thermal infrared is less sensitive to scattering and absorption by smoke or dust and invariant in case of illumination change (Socolinsky et al., 2003). It also allows to reveal anatomical information which is very useful in detecting disguised faces (Pavlidis and Symosek, 2000).The quality of a multispectral image has great implications on the efficiency of image processing applications. Environmental corruption such as noise and blur is a common phenomena for any captured images due to many factors. In particular, a multispectral image can be subject to quality degradation due to the imperfectness of sensors (Corner et al., 2003b). Particularly, noise is inevitable in all real broadband and multichannel images. Thus, it is essential to have techniques that ensure noise removal to adequate levels in order to increase performance in many image processing applications such as classification and segmentation (Abrams and Cain, 2002; Corner et al., 2003a).
In our work, we propose an improved version of the NLM filter called Optimized Vector NLM (OVNLM), where we take into account the spectral dimension of data. In fact, the NLM denoising filter exploits the redundancy of information in image. In this paper, we propose an optimization approach to tune the filter parameters. These parameters are optimized using Stein’s Unbiased Risk Estimator (SURE) rather than using ad hoc means. Furthermore, we propose a modification to the NLM filter in order to improve its performance. Indeed, the OVNLM is proposed to attenuate the computational complexity of the NLM filter by considering only pixels which are most similar to each other in computing a restored pixel. This attenuation is achieved using a similarity measure based on a probabilistic approach.
In Section 2, we present an overview of the previous works related to multispectral image denoising and give more details about our contribution. Section 3 summarizes the NLM filter and describes the proposed OVNLM filter. In Section 4, we present experimental results of our framework and compare them with similar schemes performed on multispectral images. Conclusions are drawn in Section 5.
2 Related works
Several techniques have been proposed to tackle the problem of multispectral image denoising.
The work of Luisier et al. (Luisier and Blu, 2008) represents the state of the art in multispectral image denoising. Authors proposed a denoising algorithm parameterized as a linear expansion of thresholds (Blu and Luisier, 2007). Optimization is carried out using Stein’s Unbiased Risk Estimator (SURE) (Stein, 1981; 2015, ; 2014, ). The thresholding function is point wise and wavelet based. A nonredundant orthonormal wavelet transform is applied on the noisy input image. Next, a vectorvalued thresholding of individual multichannel wavelet coefficients is performed. Finally, an inverse wavelet transform is applied to obtain the denoised image. The application of an orthonormal wavelet transform is justified by two main properties. First, assuming a white Gaussian noise in the image domain , its wavelet coefficients remain also Gaussian and are independent between subbands. Second, the Mean Square Error (MSE) in is equal to the sum of subbands’ MSEs.
Another SURE based approach was proposed in (Chaux et al., 2008). Authors used a generalized form of shrinkage estimate. The optimal choice of parameters is based on the minimization of the quadratic risk or MSE that depends on the original data which is unknown. Parameters are chosen so as to minimize the obtained risk. Note that the proposed denoising framework was built around a waveletbased approach. Two decomposing schemes were proposed: a decimated Mband wavelet transform and an Mband dualtree wavelet decomposition. For each case, the associated estimator is obtained.
Another scheme was proposed in (Duijster et al., 2009). The algorithm which jointly removes noise and blur from images is based on the Expectation and Maximization (EM) algorithm (Hogg and Craig, 1978)
. The noisy signal is decomposed into two independent parts: the first one represents the blurring problem while the second represents the denoising one. The latter is performed in wavelet domain. A Gaussian scale mixture is used to model the probability density of the wavelet coefficients. Besides that, a coregistered auxiliary noisefree image of the same scene is included in the framework to improve the restoration process. In fact, this inclusion provides an extra prior information to the model. In
(Liu et al., 2012), authors proposed a partial differential equation denoising system based on the Total Variation (TV) denoising method used in
(Chambolle, 2004) which proposes to minimize an objective function. For this, authors used the time marching method (Rudin et al., 1992). The denoising task is then reduced to a partial differential equation (PDE) problem. Authors injected in this PDE problem an auxiliary image as a prior. This approach is justified by the fact that edge directions and texture information of the auxiliary image are similar to those of the noisy image. Thus, a smoothing term that takes into account the contribution of this prior information is added. Although this approach offered better noise smoothing and details conservation, the availability of a reference image as a prior is not straightforward.A nonlocal multidimensional TV model has been recently proposed in (Li et al., 2015)
. Authors presented the denoising problem as a minimization of a mean square cost function that depends on a regularization term. The nonlocal property is not restricted to patches from one band but also to other bands with high correlation. Thus, for a given pixel, the similarity between patches from other bands is considered in the computation of the weight. The multichannel image is first divided into many groups. For a given band, bands with high correlation are grouped together. In addition, the regularization parameters are computed adaptively for each band and derived from the estimated noise standard deviation using the coefficient of the highest frequency wavelet subband. The obtained minimization problem is solved using Bregmanized operator splitting
(Zhang et al., 2010; sb1, ; sb2, ) which introduces an auxiliary variable. The unconstrained problem is treated using Bregman iteration method which leads to an update algorithm where GaussSeidel and shrinkage methods are used. The proposed framework was jointly used for multichannel image denoising and inpainting. Although the nonlocal approach offers good denoising performance, it is still computationally expensive and memory space consuming.In (Zhao and Yang, 2015)
, Zhao et al proposed a denoising framework based on sparse presentation and low rank constraint. Authors analyzed the difference in rank between a clean and a noisy image and concluded that the rank of the clean image is far smaller than the size of the multichannel image. However, this is not true for the noisy image. Thus, an assumption is made: a low rank is a characteristic of a noise free multichannel image. This information is incorporated in the cost function. Furthermore, the cost function requires patch extraction. To avoid the problem of curse of dimensionality and large error, authors suggested to reshape the 3D spectral cube into a 2D matrix by converting each band into a vector, then patches are extracted. The optimization with respect to some variables is carried out by fixing some other variables. The overall complexity of this approach is
where and are the height and length of the spatial dimension. Although good denoising performance was obtained, this approach doesn’t perform well in the presence of high level of noise since it is based on dictionary learning for sparsity representationYuan et al. studied in (Yuan et al., 2014) the noise in multichannel images, and concluded that there are two types of noise distributions: one distribution in spatial domain and one in spectral domain. Thus, two TV models are used: one applied for multichannel image denoising in spatial domain and the other one is applied in the spectral domain. The two models are both optimized with the split Bregman method where the regularization parameter is selected as the one with highest mean Peak SignaltoNoise Ratio and Structural Similarity index. Authors studied also the complementary nature of both schemes and concluded that both denoising results can complement each other and that a fusing process can bring additional improvements. By using the metric proposed in (Zhu and Milanfar, 2010), a fusion scheme between bands from each denoising result is proposed and the final denoised multichannel image is obtained. This approach exhibited good denoising performance but can be improved by adaptively adjusting the regularization parameter on which the denoising performance is highly dependent.
Yuan et al. in (Yuan et al., 2015) proposed also another denoising method where the regularization term in the cost function is often approached by a kernel model. However, this approach has three main drawbacks when applied for multichannel image denoising. First, the spectral information is not considered. Second, since the spatial resolution is lower than the spectral resolution, this approach is inefficient. Finally, noise differs from one band to another. This fact is not considered. Given these challenges, authors suggested two strategies. In the first one, a spectralspatial kernel model is considered where the spatial and spectral information are simultaneously used. In the second one, noise distributions in spectral bands are considered different and a local kernel is used to balance the contributions between bands. However, the regularization parameter which balances the contribution between the regularization term and fidelity term is not adaptively estimated.
In (Xu et al., 2015), authors proposed a denoising framework based on the Bayesian least squares optimization problem. This framework requires the computation of the posterior distribution based on Monte Carlo sampling (Lemieux, 2009)
. Given the noisy pixel, the procedure starts by choosing some neighbor pixels. Then, the acceptance probability of the sampled pixel given the noisy one is used to decide whether the sampled pixel is to be considered or not. This decision is based on a comparison between the acceptance probability and the random variable drawn from the uniform distribution. After selecting sample pixels, the importanceweighted Monte Carlo posterior estimate is computed using the weighted histogram approach proposed in
(Wong et al., 2011), then finally the denoised pixel is obtained.Peng et al. proposed in (Peng and Rao, 2009)
a vector version of the bilateral filter. The basic assumption behind this filter is that pixels which have influence on the restored pixels are not just neighbor pixels but neighbor pixels with similar values. Typically, in a similar way to Gaussian filter, bilateral filter is defined as a weighted average of neighbor pixels. However, in order to preserve edges, the difference in value with the neighbor pixels is taken into account. In their work, Peng et al. extended the bilateral filter to the vector form. The dissimilarity measure is now expressed as a multivariate Gaussian function. For simplification purposes and to avoid the computation of the noise covariance matrix, data are projected into subspace using Principle Component Analysis (PCA), and noise variance of individual channels is computed using the median absolute deviation method
(Rousseeuw and Leroy, 1987). However, this ad hoc method makes the results enormously dependent on the choice of filter parameters.Authors in (Peng et al., 2012, 2014a) proposed an optimization framework for the vector bilateral filter using SURE. They proved that within a neighborhood of a given edge pixel, a high Signal to Error (SER) measure is obtained by maximizing the weight attributed to neighbor pixels with similar values and minimizing the weight given to pixels with significant different values. Authors have also demonstrated that the SER of the vector version of the bilateral filter is always greater than the component wise 2D bilateral filter. The optimization scheme is based on the minimization of the MSE. However, the underlying difficulty of this measure is that it involves the original image which is unknown. MSE is seen as a random variable that depends on the noise. Its expected value is called the risk. To overcome this issue, filter parameters are obtained by minimizing the unbiased risk expression of the SURE estimator. The obtained minimization problem is nonlinear and is solved numerically using Sequential Quadratic Programming (SQP). Experiments on color and multispectral images have been conducted and comparison using the Peak Signal to Noise Ratio (PSNR) is presented.
Maggioni et al. (Maggioni et al., 2013) presented an extension to the BM3D denoising algorithm (Dabov et al., 2007) called BM4D. Based on the paradigm of grouping and collaborative filtering, cubes of voxels are stacked and processed in the transform domain which exploits correlation within cubes and the nonlocal correlation between the corresponding voxels of different cubes. This approach leads to an effective separation between signal and noise through coefficient shrinkage.
Peng et al. (Peng et al., 2014b)
proposed the TDL algorithm. Authors focused on the spatial nonlocal similarity and the spectral correlation of multispectral images. A nonlocal tensor dictionary learning model is developed. This model is constrained by groupblock sparsity. In addition, the proposed model is decomposed into a series of lowrank tensor approximation problems. These problems are approached using higherorder statistics.
Manjon et al. (Manjon et al., 2007) have recently proposed a new algorithm for multispectral image denoising based on the NonLocal Mean (NLM) filter (Buades et al., 2005). NLM filter is designed so that it takes advantage of the redundancy exhibited in the image. This redundancy is no longer pixel based but window based. In other words, every small window centered on a pixel is supposed to have many similar windows. These windows can be located anywhere in the image domain and are no longer restricted to the neighborhood. In the multispectral framework, information from various bands are combined and a new weight is proposed. Also NLM filter is very effective, it is highly dependent on the choice of three parameters: the radius of the search window, the radius of the neighborhood window and a smoothing parameter that controls the degree of the smoothing. The latter is very important. Indeed, with a small value, little noise will be removed. On the other hand, with a high value, the image will be blurred. Authors have set these parameters manually.
Motivated by the successful applications of NLM filter in image denoising and details preservation compared to other filters (e.g. bilateral), we propose a modified version of the NLM filter called Optimized Vector NLM (OVNLM) filter. Our contribution consists in:

Proposing a vector version of the NLM filter where we take into consideration the additional information brought by the spectral imaging system.

Automatic tuning of the filter parameters. Unlike the ad hoc method proposed in (Manjon et al., 2007), we use an optimization approach to properly choose these parameters in a way that guarantees better denoising performance.

Reducing the computation complexity. The main advantage of NLM filter is its nonlocal property which means each pixel is influenced by all pixels in the image. This comes unfortunately with more computation burden. We alleviate this burden by proposing a similarity measure used to decide whether we should take the pixel contribution or not during the pixel restoration process. Our experimental results demonstrate that this approach not only reduces the computation time but also ensures good denoising results.
We prove through quantitative evaluation the advantages of the proposed method compared to other denoising algorithms derived from the classic NLM filter as well as from other theories. Indeed, our method achieves better denoising performance compared to other algorithms. Furthermore, we show how OVNLM is capable to preserve image details while conserving its nonlocal property and ensuring acceptable computational efficiency.
3 Optimized vector NonLocal Mean filter for multispectral image denoising: OVNLM
We consider the following additive noise model:
(1) 
Where and are the noisy and original pixels respectively, is the Gaussian noise and is the pixel coordinates in the spatial domain.
3.1 NonLocal Mean filter
The basic assumption behind the definition of the NLM filter is that we need to take advantage of the high degree of redundancy in the image: the neighborhood of a pixel is any set of pixels in the image domain such that a local window surrounding is similar to the local window surrounding (Buades et al., 2005). The general case of NLM filter is given by:
(2) 
is the weight calculated for each pixel. It is computed based on a similarity measure between pixels in position and . satisfies the following constraints:

(3) 
The similarity between two pixels and is measured as a decreasing function of the Gaussian weighted Euclidean distance , where is the standard deviation of the Gaussian kernel. Let and be the pixel vectors of the gray level intensity within a squared neighborhood centered at positions and respectively.
(4) 
acts as a smoothing parameter. is a normalization constant which ensures that .
(5) 
The Gaussian weighted Euclidean distance is given by:
(6) 
where is a local window and is defined as:
(7) 
Thus, we can distinguish two main characteristics: the restored pixel is obtained by taking into account the contribution of pixels in the whole image and the weight computation is based on the similarity between local windows. Such characteristics have triggered researchers to design various novel methods (Buades et al., 2005).
3.2 Vector NLM filter
To take advantage of the additional information brought by the spectral dimension, we extend the NLM filter to the vector case. In the multispectral context, we have the reflectance intensity at a given position in different spectral bands. Thus, we are operating on a set of pixel vectors . We define the vector NLM (VNLM) filter as:
(8) 
where the new formulation of the weight between two pixels at position and is defined as:

(9) 
If , where
is the identity matrix, we get the classical Euclidean distance.

(10) 
3.3 Optimization framework for vector NLM
In our framework design, we target two main objectives: optimize the parameters of the filter and reduce its computation complexity. First we use both the classical Euclidean distance and Mahalanobis distance where is a covariance matrix. In addition, we preselect for each pixel a subset of the most similar pixels based on a probabilistic similarity measure.
The filter depends on two parameters: the smoothing parameter and the covariance matrix . Thus, we have:

(11) 
where is a nonlinear estimator and is the filter parameter.
Our aim is to optimize so that we can ensure the choice of the optimal parameters in order to obtain the best denoising result. The performance of the estimator is generally evaluated using the mean square error (MSE):
(12) 
However, the problem of such estimator is that the ground truth image is unknown. MSE can be seen as a random variable of the noise. Its expected value is designated as the Risk and expressed as:
(13) 
The problem of estimating the risk without the need to have the underlying image is approached by Stein’s Unbiased Risk Estimator (SURE) (Chaux et al., 2008; Luisier and Blu, 2008). Thus, we have (Stein, 1981):

(14) 
and:

(15) 
is the transpose operator. If we consider zero mean multivariate Gaussian noise, we get (Luisier and Blu, 2008):
(16) 
where is the noise covariance matrix.
By combining eq. 14 and eq. 15, we end up with an expression without :

(17) 
Therefore, the risk is the unbiased risk estimator of MSE in eq. 12 and is given by:

(18) 
where is the Jacobian matrix with respect to . is given by (Peng et al., 2014a):

(19) 
where is the delta function and is defined as:

(20) 
With the derivation of (see appendix), we formulate the problem of vector NLM filter as a constrained optimization problem:
(21) 
Note that in case of using the Euclidean distance, the only parameter to be optimized is .
3.4 Relevant pixel selection
If we go back to eq. 8, we can clearly see that in order to restore every pixel, we need to go through every other pixel in the domain . This is obviously a very time consuming process.
To attenuate the computation complexity of the proposed VNLM filter, we propose to preselect for each processed pixel, a set of most relevant pixels based on the similarity measure proposed in (Matsushita and Lin, 2007). This measure is based on a probabilistic approach to compute the similarity between two pixels based on the noise distribution.
In the grayscale case, the similarity measure is defined as:

(22) 
where is the maximum value of the true intensity, is a constant, and is the error function defined as:
(23) 
Fig. 1 illustrates the form of the similarity measure for =100.
For a given , eq. 22 illustrates a Gaussian function. We consider that all values beyond the width of 1/ of maximum () are zeros. In the case of RGB color images, the similarity between two pixels and is defined as:
(24) 
We generalize this similarity measure for the multispectral case, such that the similarity measure between and is defined as:
(25) 
The proposed Optimized VNLM (OVNLM) filter becomes:
(26) 
3.5 OVLNM algorithm
The proposed approach is detailed in what follows. We solve the constrained nonlinear optimization problem using Sequential Quadratic Programming. Given a noisy image, the noise covariance matrix is estimated with the standard median absolute deviation method (Huber et al., 1981). The diagonal elements of are calculated as follows:
(27) 
. The offdiagonal elements are defined as:

(28) 
where: and and . We minimize the risk value based on an optimal choice of parameters until we reach the maximum number of iteration or the risk value decreases below a preset threshold . We implemented this approach in Matlab (R2015a). The minimization is conducted using the function with the risk as an objective function to minimize and SQP as the optimization approach. We use a neighborhood window of and we set up .
4 Experiments
4.1 Data sets
To assess the performance of our approach, we conducted experiments on real world multispectral images. In one of the experiments, we used the Salinas scene collected using the Airborne Visible InfraRed Imaging Spectrometer (AVIRIS) at NASA’s Jet Propulsion Laboratory^{1}^{1}1http://aviris.jpl.nasa.gov/. Sample bands are shown in Fig. 2. The Salinas Valley image is a high spatial resolution image consisting of a collection of 224 spectral band images taken over Salinas Valley California in the range from 0.4m to 2.5m at a resolution of 3.7 meters per pixel. Spectral Bands [108112], [154167] and 220 are discarded due to water absorption and noise. Before processing, images are resized to pixels. In another experiment we used multispectral face images from the IRIS Lab database at the University of Tennessee Chang et al. (2006). The IRIS Lab database was built between August 2005 and March 2006 and consists of 2624 multispectral face images taken along the visible spectrum in addition to thermal images with a resolution of pixels. RGB images are also generated with a resolution of pixels. These images are taken in different lighting conditions: Halogen light, daylight and fluorescent light. The total size of the database is 8.91GB. A total of 82 participants were involved from different genders (76% male, 24% female), ethnicities as depicted in Table 1, ages, facial expression, genders and hair characteristics. We conducted experiments on 8 multispectral images referred to as . Figures 3 and 4 illustrate the multispectral image samples taken with Halogen light and used in our experiments.
Caucasian  Asian  Asian Indian  African descent  
%  57%  23%  12%  8% 
4.2 Evaluation results
We applied the proposed OVLNM approach and compared its denoising performance with several state of the art multispectral image denoising algorithms: MNLM (Manjon et al., 2007), the multichannel SURELET (MSURE) (Luisier and Blu, 2008), BM4D (Maggioni et al., 2013) and TDL (Peng et al., 2014b). Note that MNLM is also inspired from the NLM filter and adapted for multispectral image denoising with choice of parameters conducted using ad hoc means. Experiments on multispectral image denoising are conducted by contaminating original images with an additive Gaussian noise at different levels, then denoising algorithms are applied on the noisy images. We use the Peak Signal to Noise ratio (PSNR) expressed in dB and the Structure Similarity Index Measure (SSIM):
(29) 
(30) 
where and are pixel values at position in the original and output images respectively.
PSNR is the ratio of the maximum possible value of the signal in term of its power and the power of the distortion caused by the noise. PSNR is expressed in the logarithmic decibel scale. The higher the PSNR is, the better is the result. SSIM is an index that measures the similarity between two images and with , , and are respectively the mean of image , the mean of image , the standard deviation of image and the standard deviation of image . and are two stabilization constants and is the covariance of and . The best result is highlighted by the highest value of SSIM.
We show in Fig. 5 the denoising results of the Salinas image with PSNR=19 dB. We notice that oversmoothness is witnessed with MSURE, BM4D and more with TDL. Some patterns in the images are better preserved particularly with MNLM and OVNLM. Indeed, one of the intrinsic property of an NLM based method is its ability to preserve image details beside its good performance even with high noise level. These characteristics are now strengthened with the fine tuning of the parameters.
In our experiments, we noticed that more detail preservation can be obtained with TDL by reducing the voxel block size but this comes with a loss in term of PSNR. Fig. 6 presents the PSNR of denoising results at different level of input PSNR. We can clearly see that OVNLM exhibits the best performance with high PSNR (low noise level) and also with low PSNR (high noise level). This performance is confirmed with the SSIM results presented in Tab. 2.
Input PSNR  MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.82  0.87  0.88  0.85  0.83 
19  0.77  0.80  0.83  0.80  0.77 
17  0.75  0.74  0.78  0.77  0.75 
15  0.70  0.66  0.72  0.74  0.73 
13.5  0.66  0.58  0.66  0.68  0.70 
12.5  0.62  0.53  0.61  0.66  0.67 
12  0.60  0.49  0.59  0.64  0.66 
11  0.56  0.43  0.53  0.58  0.62 
10.5  0.53  0.40  0.51  0.56  0.60 
Denoising results for image from the IRIS Lab multispectral face image are illustrated in Fig. 7 with input PSNR=19 dB. Performance evaluation is shown in Fig. 8. It demonstrates the improvement of OVNLM over the other algorithms with different noise levels and with distinctive performance in low PSNR (11 dB and 10.5 dB). SSIM results presented in Tab. 3 also demonstrate the outperformance of OVNLM. With low PSNR, MSURE exhibits the lowest performance with SSIM=0.47 while with OVNLM we have SSIM=0.63. Fig. 9 shows the outputs of the denoising algorithms with input PSNR=19 dB for . OVNLM presents the best performance with 33.07 dB compared to 28.4 dB for MSURE, 27.7 dB for MNLM, 27.5 dB for BM4D and 25 dB for TDL. This performance is confirmed by Fig. 10 which illustrates the PSNR variations with respect to the input PSNR. SSIM results in Tab. 4 confirm the outperformance of OVNLM particularly in presence of high noise level. For example, with input PSNR=10.5 dB, we have SSIM=0.67 for OVNLM which is better than the rest of the algorithms.
Input
PSNR 
MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.88  0.90  0.87  0.85  0.86 
19  0.82  0.85  0.82  0.79  0.82 
17  0.75  0.80  0.77  0.74  0.78 
15  0.75  0.73  0.73  0.68  0.76 
13.5  0.69  0.66  0.69  0.65  0.73 
12.5  0.64  0.61  0.66  0.61  0.70 
12  0.62  0.58  0.64  0.58  0.69 
11  0.56  0.51  0.60  0.55  0.65 
10.5  0.53  0.47  0.58  0.51  0.63 
Input PSNR  MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.86  0.88  0.87  0.87  0.88 
19  0.82  0.82  0.82  0.83  0.84 
17  0.78  0.77  0.80  0.77  0.81 
15  0.72  0.73  0.76  0.72  0.79 
13.5  0.67  0.66  0.68  0.72  0.76 
12.5  0.63  0.61  0.62  0.70  0.74 
12  0.60  0.58  0.61  0.68  0.72 
11  0.56  0.51  0.57  0.62  0.69 
10.5  0.53  0.47  0.54  0.60  0.67 
Fig. 12 shows PSNR variations of from IRIS Lab multispectral face image database. These results confirm the outperformance of OVNLM algorithm over the other algorithms. This performance is emphasized by the SSIM results for presented in Tab. LABEL:tab_sub3,tab_sub4,tab_sub5,tab_sub6,tab_sub7,tab_sub8. According to these SSIM values, denoised images obtained with OVNLM are the most similar ones to the original images.
To emphasize the importance of a good parametrization of the proposed filter, we analyzed its performance by varying parameter which controls the degree of smoothing for a fixed covariance matrix with different noise levels as shown in Fig. 11 for multispectral image. We noticed that the resulting PSNR depends enormously on the choice of parameter which also implies a good choice of parameter .
The overall experimental results confirm the following findings. First, NLM based algorithm (MNLM) without a proper choice of parameters cannot deal with high level noise. Second, dictionary methods are known by their lack of robustness which is confirmed by TDL results. MSURE, a waveletbased method and BM4D, a grouping and collaborative filtering method, do not exhibit stable denoising performance. Furthermore, SSIM results of MSURE showed that this algorithm does not preserve details. Indeed, for and with PSNR=10.5, MSURE exhibited the lowest value which is too much smaller than BM4D, TDL and OVNLM values. For BM4D, some details such as the ones of Salinas image in Fig. 5 are lost. The experimental findings reflect also the intrinsic characteristics of the proposed method: (i) OVNLM performs better in term of details preservation. This is quantitatively confirmed by the SSIM results. (ii) The fine tuning of the parameters allows OVNLM to achieve better denoising performance particularly with high level noise. Indeed, by comparing OVNLM results to MNLM, we deduce the importance of setting up a proper choice for these parameters. (iii) Our relevant pixel selection strategy conserves the nonlocal property but chooses the most influential pixels for each pixel restoration. PSNR results have demonstrated the effectiveness of our approach. Next, we focus on the computation complexity of our method and compare it with other NLM based algorithms.
4.3 Computation Complexity
To evaluate the relevance of choosing a particular subset for each pixel in term of computation complexity, we evaluate the computation time of the OVNLM algorithm and compare it against the computation time of several variations of the NLM filter: classic NLM filter applied on each spectral band separately, the MNLM algorithm, the proposed OVNLM algorithm applied with and without selecting subsets for each pixel and referred as OVNLM v1 and OVNLM v2 respectively.
We use the multispectral image of with input PSNR=22.5 dB. We test the algorithms on a Windows machine, Intel Core i7, 3.1 GHz with 8GB RAM. Results are illustrated in Table 5.
Method  Computational time (s)  PSNR 

NLM  195  30.21 
MNLM  35  30.29 
OVNLM v1  29  32.03 
OVNLM v2  175  33.1 
First of all, we notice that MNLM and OVNLM v1 exhibit the best performance in term of computational time. Classic NLM filter and OVNLM v2 need to consider the contribution of every pixel in the image domain which justifies their high computational time. Results also demonstrate that the computational reduction that characterizes OVLM v2 comes with a slight decrease of PSNR compared to OVNLM v2. SURE procedure itself does not greatly affect the computation time which is also confirmed in (Peng et al., 2014a; Van De Ville and Kocher, 2009).
We analyze the effect of the choice of the parameter which controls the amount of points to be considered when restoring each pixel. We experiment on multispectral image. The findings in Tab. 6 show that at a certain level ( 100), we achieve a considerable improvement in term of computational time. Meanwhile, the improvement in term of PSNR is not very relevant (around 1dB). We conclude that a proper adjustment of this parameter can ensure a good denoising performance while keeping a reasonable computational time.
1  10  50  100  500  

PSNR  27.38  28.52  29.76  32.88  33.19  33.34  33.83 
Time (s)  20  33  43  47  57  60  68 
Finally, we can state that the proposed denoising filter demonstrated its effectiveness compared to other algorithms. Experiments on real multispectral images have shown good results in terms of PSNR. While denoising, we tried to take advantage of the spectral information by extending the NLM filter to the vector case, optimized the choice of its parameter and reduced the computation complexity.
5 Conclusion
In this paper, we have proposed a novel multispectral image denoising algorithm. The proposed scheme is an improvement of NonLocal Mean filter. This improvement is obtained by extending NLM filter to the vector case.
However, the choice of the filter tuning parameters still poses a challenge. In this work,
filter parameters are tuned automatically using an optimization technique based on SURE.
Indeed, unbiased risk estimator is applied and filter parameters are obtained by minimizing the expression of SURE. An easily sequential quadratic programming is used to solve the nonlinear minimization problem.
Experiments performed on color and multispectral face images demonstrate the superiority of the proposed framework compared with two other wellknown similar algorithms.
Good performance in terms of PSNR and SSIM is obtained. Nevertheless, the computation burden is still an important challenge with NLM in general. Further techniques to speed up the computation should be investigated.
Acknowledgment
This publication was made possible by NPRP grant # 41165 2453 from the Qatar National Research Fund (a member of Qatar Foundation). The statements made herein are solely the responsibility of the authors.
Input
PSNR 
MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.86  0.87  0.88  0.85  0.88 
19  0.81  0.79  0.83  0.80  0.84 
17  0.76  0.73  0.80  0.76  0.82 
15  0.70  0.66  0.77  0.72  0.79 
13.5  0.64  0.60  0.72  0.64  0.77 
12.5  0.59  0.54  0.69  0.63  0.74 
12  0.56  0.51  0.68  0.59  0.72 
11  0.52  0.45  0.63  0.58  0.69 
10.5  0.49  0.41  0.60  0.55  0.67 
Input PSNR  MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.85  0.86  0.87  0.85  0.87 
19  0.80  0.79  0.83  0.79  0.83 
17  0.75  0.65  0.79  0.76  0.81 
15  0.68  0.57  0.75  0.69  0.77 
13.5  0.62  0.56  0.70  0.64  0.73 
12.5  0.57  0.50  0.67  0.62  0.69 
12  0.55  0.47  0.65  0.58  0.67 
11  0.50  0.40  0.60  0.57  0.61 
10.5  0.47  0.37  0.57  0.56  0.59 
Input PSNR  MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.89  0.89  0.89  0.88  0.90 
19  0.85  0.84  0.86  0.83  0.88 
17  0.81  0.79  0.82  0.79  0.86 
15  0.76  0.72  0.78  0.75  0.81 
13.5  0.70  0.66  0.75  0.70  0.80 
12.5  0.66  0.61  0.72  0.65  0.78 
12  0.64  0.58  0.69  0.64  0.76 
11  0.58  0.51  0.65  0.61  0.71 
10.5  0.56  0.48  0.64  0.60  0.69 
Input PSNR  MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.87  0.89  0.89  0.88  0.89 
19  0.83  0.83  0.85  0.83  0.86 
17  0.79  0.78  0.82  0.79  0.83 
15  0.74  0.72  0.77  0.74  0.80 
13.5  0.69  0.66  0.74  0.70  0.77 
12.5  0.65  0.61  0.72  0.66  0.74 
12  0.63  0.59  0.69  0.62  0.72 
11  0.58  0.53  0.65  0.58  0.68 
10.5  0.56  0.50  0.62  0.57  0.66 
Input PSNR  MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.89  0.90  0.89  0.89  0.90 
19  0.85  0.84  0.85  0.83  0.87 
17  0.82  0.80  0.83  0.78  0.85 
15  0.77  0.74  0.78  0.75  0.82 
13.5  0.72  0.66  0.74  0.70  0.78 
12.5  0.68  0.64  0.72  0.67  0.77 
12  0.66  0.62  0.69  0.64  0.75 
11  0.61  0.56  0.64  0.58  0.71 
10.5  0.58  0.52  0.63  0.58  0.68 
Input PSNR  MNLM  MSURE  BM4D  TDL  OVNLM 

22.5  0.88  0.89  0.89  0.88  0.88 
19  0.83  0.83  0.85  0.83  0.83 
17  0.79  0.78  0.82  0.79  0.82 
15  0.72  0.70  0.78  0.73  0.79 
13.5  0.66  0.63  0.74  0.68  0.75 
12.5  0.61  0.57  0.71  0.66  0.73 
12  0.59  0.54  0.70  0.62  0.71 
11  0.53  0.47  0.63  0.58  0.67 
10.5  0.50  0.43  0.62  0.55  0.64 
Appendix
Let: .


References
References
 Abrams and Cain (2002) Abrams, M.C., Cain, S.C., 2002. Sampling, radiometry, and image reconstruction for polar and geostationary meteorological remote sensing systems. Proc. SPIE 4792, 207–215.
 Blu and Luisier (2007) Blu, T., Luisier, F., 2007. The SURELET Approach to Image Denoising. IEEE Transactions on Image Processing 16, 2778–2786.
 Buades et al. (2005) Buades, A., Coll, B., Morel, J., 2005. A review of image denoising algorithms, with a new one. Multiscale Modeling & Simulation 4, 490–530.
 Chambolle (2004) Chambolle, A., 2004. An algorithm for total variation minimization and applications. Journal of Mathematical Imaging and Vision 20, 89–97.

Chang et al. (2006)
Chang, H., Harishwaran, H.,
Yi, M., Koschan, A.,
Abidi, B., Abidi, M.,
2006.
An indoor and outdoor, multimodal, multispectral and multiilluminant database for face recognition, in: Computer Vision and Pattern Recognition Workshop, 2006. CVPRW ’06. Conference on, pp. 54–54.
 Chaux et al. (2008) Chaux, C., Duval, L., BenazzaBenyahia, A., Pesquet, J., 2008. A nonlinear steinbased estimator for multichannel image denoising. Signal Processing, IEEE Transactions on 56, 3855–3870.
 Corner et al. (2003a) Corner, B., Narayanan, R., Reichenbach, S., 2003a. Noise estimation in remote sensing imagery using data masking. International Journal of Remote Sensing 24, 689–702.
 Corner et al. (2003b) Corner, B.R., Narayanan, R.M., Reichenbach, S.E., 2003b. Noise estimation in remote sensing imagery using data masking. International Journal of Remote Sensing 24, 689–702.
 Dabov et al. (2007) Dabov, K., Foi, A., Katkovnik, V., Egiazarian, K., 2007. Image denoising by sparse 3d transformdomain collaborative filtering. Image Processing, IEEE Transactions on 16, 2080–2095.
 Duijster et al. (2009) Duijster, A., Scheunders, P., De Backer, S., 2009. Waveletbased em algorithm for multispectralimage restoration. Geoscience and Remote Sensing, IEEE Transactions on 47, 3892–3898.
 Hogg and Craig (1978) Hogg, R.V., Craig, A.T., 1978. Introduction to mathematical statistics; 4th ed. Macmillan, New York, NY.
 Huber et al. (1981) Huber, P., Wiley, J., InterScience, W., 1981. Robust statistics. Wiley New York.
 Kong et al. (2007) Kong, S., Heo, J., Boughorbel, F., Zheng, Y., Abidi, B., Koschan, A., Yi, M., Abidi, M., 2007. Multiscale fusion of visible and thermal ir images for illuminationinvariant face recognition. International Journal of Computer Vision 71, 215–233.
 Kong et al. (2005) Kong, S.G., Heo, J., Abidi, B.R., Paik, J., Abidi, M.A., 2005. Recent advances in visual and infrared face recognition—a review. Computer Vision and Image Understanding 97, 103 – 135.
 Lemieux (2009) Lemieux, C., 2009. Monte Carlo and QuasiMonte Carlo Sampling. SpringerVerlag New York.
 Li et al. (2015) Li, J., Yuan, Q., Shen, H., Zhang, L., 2015. Hyperspectral image recovery employing a multidimensional nonlocal total variation model. Signal Processing 111, 230 – 248.
 Liu et al. (2012) Liu, P., Huang, F., Li, G., Liu, Z., 2012. Remotesensing image denoising using partial differential equations and auxiliary images as priors. Geoscience and Remote Sensing Letters, IEEE 9, 358–362.
 Luisier and Blu (2008) Luisier, F., Blu, T., 2008. Surelet multichannel image denoising: Interscale orthonormal wavelet thresholding. Image Processing, IEEE Transactions on 17, 482–492.
 Maggioni et al. (2013) Maggioni, M., Katkovnik, V., Egiazarian, K., Foi, A., 2013. Nonlocal transformdomain filter for volumetric data denoising and reconstruction. Image Processing, IEEE Transactions on 22, 119–133.
 Manjon et al. (2007) Manjon, J.V., Robles, M., Thacker, N.A., 2007. Multispectral mri denoising using nonlocal means. Medical Image Understanding and Analysis , 41–46.
 Matsushita and Lin (2007) Matsushita, Y., Lin, S., 2007. A probabilistic intensity similarity measure based on noise distributions, in: Computer Vision and Pattern Recognition, 2007. CVPR ’07. IEEE Conference on, pp. 1–8.
 Meraoumia et al. (2014) Meraoumia, A., Chitroub, S., Bouridane, A., 2014. Biometric recognition systems using multispectral imaging, in: Hassanien, A.E., Kim, T.H., Kacprzyk, J., Awad, A.I. (Eds.), Bioinspiring Cyber Security and Cloud Services: Trends and Innovations. Springer Berlin Heidelberg. volume 70, pp. 321–347.
 Pan et al. (2003) Pan, Z., Healey, G., Prasad, M., Tromberg, B., 2003. Face recognition in hyperspectral images. Pattern Analysis and Machine Intelligence, IEEE Transactions on 25, 1552–1560.
 Pan et al. (2004) Pan, Z., Healey, G.E., Prasad, M., Tromberg, B.J., 2004. Hyperspectral face recognition under variable outdoor illumination. Proc. SPIE 5425, 520–529.
 Pavlidis and Symosek (2000) Pavlidis, I., Symosek, P., 2000. The imaging issue in an automatic face/disguise detection system, in: Computer Vision Beyond the Visible Spectrum: Methods and Applications, 2000. Proceedings. IEEE Workshop on, pp. 15–24.
 Peng and Rao (2009) Peng, H., Rao, R., 2009. Hyperspectral image enhancement with vector bilateral filtering, in: Image Processing (ICIP), 2009 16th IEEE International Conference on, pp. 3713–3716.
 Peng et al. (2012) Peng, H., Rao, R., Dianat, S., 2012. Optimized vector bilateral filter for multispectral image denoising, in: Image Processing (ICIP), 2012 19th IEEE International Conference on, pp. 2141–2144.
 Peng et al. (2014a) Peng, H., Rao, R., Dianat, S., 2014a. Multispectral image denoising with optimized vector bilateral filter. Image Processing, IEEE Transactions on 23, 264–273.
 Peng et al. (2014b) Peng, Y., Meng, D., Xu, Z., Gao, C., Yang, Y., Zhang, B., 2014b. Decomposable nonlocal tensor dictionary learning for multispectral image denoising, in: Computer Vision and Pattern Recognition (CVPR), 2014 IEEE Conference on, pp. 2949–2956.

Rousseeuw and Leroy (1987)
Rousseeuw, P.J., Leroy, A.M.,
1987.
Robust Regression and Outlier Detection.
John Wiley & Sons, Inc., New York, NY, USA.  Rudin et al. (1992) Rudin, L.I., Osher, S., Fatemi, E., 1992. Nonlinear total variation based noise removal algorithms. Phys. D 60, 259–268.
 Socolinsky et al. (2003) Socolinsky, D.A., Selinger, A., Neuheisel, J.D., 2003. Face recognition with visible and thermal infrared imagery. Computer Vision and Image Understanding 91, 72 – 114.

Stein (1981)
Stein, C.M., 1981.
Estimation of the mean of a multivariate normal distribution.
The Annals of Statistics 9, pp. 1135–1151.  (34) Ulfarsson, M.O., Solo, O.,2015. Selecting the Number of Principal Components with SURE. IEEE Signal Processing Letters 22, pp. 239–243.
 (35) Wue, Y., Tracey, B.H., Natarajan, P., Noonan, J.P.,2014. Fast blockwise SURE shrinkage for image denoising. Signal Processing 103, pp. 45–59.
 Van De Ville and Kocher (2009) Van De Ville, D., Kocher, M., 2009. Surebased nonlocal means. Signal Processing Letters, IEEE 16, 973–976.

Wong et al. (2011)
Wong, A., Mishra, A.,
Zhang, W., Fieguth, P.,
Clausi, D.A., 2011.
Stochastic image denoising based on markovchain monte carlo sampling.
Signal Processing 91, 2112 – 2120.  Xu et al. (2015) Xu, L., Li, F., Wong, A., Clausi, D., 2015. Hyperspectral image denoising using a spatialspectral monte carlo sampling approach. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of PP, 1–14.
 Yuan et al. (2014) Yuan, Q., Zhang, L., Shen, H., 2014. Hyperspectral image denoising with a spatialspectral view fusion strategy. Geoscience and Remote Sensing, IEEE Transactions on 52, 2314–2325.
 Yuan et al. (2015) Yuan, Y., Zheng, X., Lu, X., 2015. Spectralspatial kernel regularized for hyperspectral image denoising. Geoscience and Remote Sensing, IEEE Transactions on 53, 3815–3832.
 Zhang et al. (2010) Zhang, X., Burger, M., Bresson, X., Osher, S., 2010. Bregmanized nonlocal regularization for deconvolution and sparse reconstruction. SIAM Journal on Imaging Sciences 3, 253–276.
 (42) Hao, Y., Xu, J., Li, S., Zhang, X., 2015. A variational model based on split Bregman method for multiplicative noise removal. International Journal of Electronics and Communications 69, 1291–1296.

(43)
He, L., Wang, Y.,
2015.
Iterative support detectionbased split bregman method for wavelet framebased image inpainting.
IEEE Transactions on Image Processing 69, 1291–1296.  Zhao and Yang (2015) Zhao, Y.Q., Yang, J., 2015. Hyperspectral image denoising via sparse representation and lowrank constraint. Geoscience and Remote Sensing, IEEE Transactions on 53, 296–308.
 Zhu and Milanfar (2010) Zhu, X., Milanfar, P., 2010. Automatic parameter selection for denoising algorithms using a noreference measure of image content. Image Processing, IEEE Transactions on 19, 3116–3132.
Comments
There are no comments yet.