Denoising of Impulsive noise in single/multichannel images
Recently, so called annihilating filer-based low rank Hankel matrix (ALOHA) approach was proposed as a powerful image inpainting method. Based on the observation that smoothness or textures within an image patch corresponds to sparse spectral components in the frequency domain, ALOHA exploits the existence of annihilating filters and the associated rank-deficient Hankel matrices in the image domain to estimate the missing pixels. By extending this idea, here we propose a novel impulse noise removal algorithm using sparse + low rank decomposition of an annihilating filter-based Hankel matrix. The new approach, what we call the robust ALOHA, is motivated by the observation that an image corrupted with impulse noises has intact pixels; so the impulse noises can be modeled as sparse components, whereas the underlying image can be still modeled using a low-rank Hankel structured matrix. To solve the sparse + low rank decomposition problem, we propose an alternating direction method of multiplier (ADMM) method with initial factorized matrices coming from low rank matrix fitting (LMaFit) algorithm. To adapt the local image statistics that have distinct spectral distributions, the robust ALOHA is applied patch by patch. Experimental results from two types of impulse noises - random valued impulse noises and salt/pepper noises - for both single channel and multi-channel color images demonstrate that the robust ALOHA outperforms the existing algorithms up to 8dB in terms of the peak signal to noise ratio (PSNR).READ FULL TEXT VIEW PDF
The goal of this paper is to develop a novel numerical method for effici...
Denoising images contaminated by the mixture of additive white Gaussian ...
We consider the case of inpainting single depth images. Without correspo...
In this paper, we propose a new method to reconstruct a signal corrupted...
Transform Invariant Low-rank Textures (TILT) is a novel and powerful too...
We present a mathematical model to decompose a longitudinal deformation ...
In this paper, we propose a new method to reconstruct a signal corrupted...
Denoising of Impulsive noise in single/multichannel images
Impulse noises occur by malfunctioning of detector pixels in camera or from the missing memory elements in imaging hardware . There are two types of impulse noises. The first type includes salt/pepper noises which have the extremal values of dynamic ranges; therefore, the noise pixels can be easily detected by an adaptive median filter (AMF). The second type noises are random valued impulse noises (RVIN) which have random values within the dynamic range of an image pixel. Unlike the salt/pepper noises, the noisy pixels of RVIN cannot be effectively detected by an adaptive median filter. Instead, the adaptive center-weighted median filter (ACWMF) has been widely used to find the locations of noisy pixels. Even with AMF and ACWMF, when the density of noise increases, the denoising performance of these single step algorithms become severely degraded. To compensate for this weakness, two-phase denoising algorithms with “decision-based filter” or “switching filter” were proposed [4, 5, 6, 7, 1]
. More specifically, these algorithms consist of two main parts: detecting noise pixels by AMF, ACWMF, or other outlier finding algorithms; and then replacing the detected noise pixels with the estimated values using the total variation or edge preserving regularizations [1, 8], while leaving other noiseless pixels unchanged.
On the other hand, impulse noise denoising algorithms using proximal optimizations with non-smooth penalties were recently proposed [9, 10, 11]. In particular, in the TVL1 (total variation ) approach [10, 9], the data fidelity term was measured with the norm to deal with impulse outliers and the total variation regularization was used as a penalty. They can effectively remove the impulse noises at sufficiently fast speed. However, the algorithm often causes edge distortions or texture patterns blurrings due to the TV terms. With the advance of compressed sensing theory [12, 13], the impulse noise denoising methods based on compressed sensing were also proposed [14, 15]. In , the authors encouraged the spatio-spectral domain redundancy to be sparsely represented in Fourier domain by using blind compressed sensing framework. This approach demonstrated outstanding recovery performances; however, the algorithm cannot be performed without highly correlated spectral dataset. In , the sparsity level of target signal was minimized just as in the conventional CS approach using noisy measurements contaminated by impulse noise. However, the performance was inferior to the two-phase methods . While a low-rank matrix completion approach for impulse noise denoising for video sequence was proposed , the algorithm only worked for video sequence denoising, because spatio-temporal redudancy should be exploited in this method.
One of the unique characteristics of impulse noises compared to Gaussian or Poisson noise is that an image corrupted with impulse noises has intact pixels; so the impulse noises can be modeled as sparse components, whereas the underlying intact image still retains the original image characteristics. In fact, this was the main idea utilised in TVL1 approach , and this paper also aims at fully exploiting this observation. However, one of the most important contributions of this paper is to demonstrate that there exists a more natural and powerful image modeling method than the TV approach, and that the resulting impulse noise removal problem becomes a sparse + low rank decomposition problem. This is inspired by our novel image modeling and associated inpainting method using so-called the annihilating filter-based low-rank Hankel matrix (ALOHA) approach .
More specifically, in , we demonstrated that the smoothness or textures within an image patch leads to sparse spectrum in the frequency domain, so the sampling theory of signals with finite rate of innovations [18, 19] tells us that there exists annihilating filters that annihilate the pixel values within the corresponding image patch. Moreover, the existence of annihilating filter enables us to construct a rank-deficient Hankel structured matrix whose rank is determined by the sparsity in the spectral domain . Thanks to this observation, an image patch could be modeled using an annihilating filter-based Hankel structured matrix (ALOHA), and image inpainting problems was solved by a low rank matrix completion algorithm. The idea was extended by our group for compressed sensing MRI [20, 21], image deconvolution 
, and interpolation of scanning microscopy. Ongie et al
independently developed similar approaches for super-resolution MRI[24, 25].
One of the most important consequences of ALOHA in the context of impulse noise removal is an observation saying that the construction of Hankel structured matrix is a kind of linear lifting scheme, so the sparse components in image are also sparse in the lifted Hankel matrix. Therefore, we can use a sparse+low rank decomposition of the Hankel structured matrix to decouple the sparse impulse noise components from the underlying image. The new algorithm, what we call robust ALOHA, is applied patch by patch to adapt the local image statistics that have distinct spectral distribution. We are aware that there have been significant progresses on the decomposition of superposed matrix consisting of low-rank and sparse components [26, 27, 28, 29, 30], which is often called Robust Principal Component Analysis (RPCA) . However, the matrix in RPCA is usually unstructured, whereas the robust ALOHA uses an Hankel structured matrix. As will be shown later, the introduction of Hankel structure matrix significantly improves the denoising performance by exploiting the spectral domain sparsity.
During the writing of this paper, we came across a very interesting sparse and low rank decomposition of Hankel structure matrix for predicting target locations under the occlusion. While the optimization framework and algorithm used in  have similarities with ours, the idea of  was originally inspired by the dynamic system identification rather than image modeling using annihilating filter and we are not aware of its applications for impulse noise removal. Therefore, we believe that the application of robust ALOHA to impulse noise removal is sufficiently novel.
To solve the associated sparse+low rank decomposition problem of Hankel structure matrix, an alternating direction method of multiplier (ADMM) is utilized with initial factorized matrices from low rank matrix fitting algorithm (LMaFit) . Furthermore, the denoising algorithm is also extended to exploit the joint sparsity constraint in colour images by stacking Hankel structured matrix from each channel side by side and applying sparse+low rank decomposition of the concatenated Hankel matrix. Using extensive numerical experiments, we demonstrate that the robust ALOHA significantly outperforms all the existing methods.
This paper is organized as follows. Section II discusses the theory behind the robust ALOHA. In Section III, an optimization method for the associated sparse+low rank decomposition of Hankel structured matrix is described. Extension to multi-channel denoising is discussed in Section IV. Experimental results are provided in Section V, which is followed by discussion and conclusion in Section VI and VII.
Because impulse noises occur by malfunctioning of detector or memory elements , only a subset of image pixels are corrupted by the noises. Therefore, if denotes an image measurement corrupted by impulse noises, it can be modeled as
where is the underlying “clean” image, and denotes a sparse matrix composed of impulse noises. This model is quite often used in the existing impulse noise removal algorithm. For example, in TVL1 [10, 9], is considered as a sparse outlier, whereas the underlying image is modeled using the total variations. This leads to the following minimization problem:
where the norm is the norm corresponding to summation of absolute values of each matrix elements for outlier removal, and the denotes the 2-D TV penalty to model the underlying image. In the following, we will explain how the image model in (1) is modified in the proposed method to give superior denoising performance.
In our recent work , we demonstrate that diffusion [34, 35, 36, 37] and/or Gaussian Markov random field (GMRF) approaches for image modelling [38, 39, 40, 41, 42] are closely related to an annihilating filter relationship from the sampling theory of signals with finite rate of innovations (FRI) [18, 19]. More specifically, as shown in Fig. 1(a), a smoothly varying patch usually has spectrum content in the low frequency regions, and the other frequency regions have very little spectral components. A similar spectral domain sparsity can be observed in a texture patch in Fig. 1(b), where the spectral components are mainly concentrated on the fundamental frequencies of the patterns. For the case of an abrupt transition along the edge as shown in Fig. 1(c), the spectral components are mostly localized along the axis.
Mathematically, when a patch has sparse spectral components, we can show that there exists a corresponding annihilating filter in the image domain. More specifically, if the spectrum of an image patch is described by
where denotes the number of non-zero spectral components, then it is easy to find an annihilating function in the spectrum domain that does not overlap with the support of , i.e.
This implies the existence of the annihilating filter in the image domain:
For example, if an image is sufficiently flat with little variations in pixel values, then its spectral component is most concentrated around the zero frequency component (i.e. ); therefore, becomes an annihilating function in the spectral domain. The associated annihilating filtering in the image domain is then given by
where denotes the Laplacian operator, which corresponds to the diffusion operation that is widely used in image denoising, inpainting, and so on [34, 35, 36, 37, 9, 10, 11]. This example clearly shows why the diffusion based approach is closely related to the annihilating filter approach.
If Fourier measurement data is discretized at an appropriate Nyquist sampling rates, the corresponding discrete counterpart is given by
where the discrete filter is now a discrete annihilating filter. Among the various of choices of annihilating function that satisfies (3), Vetterli et al [18, 19] showed that an annihilating function can be constructed using a finite combination of sinusoidals such that the corresponding discrete annihilating filter has a finite filter length. In this case, the convolution (6) becomes finite length convolution, and we can exploit (6) in a matrix representation.
Specifically, let denote the matrix composed of such that . We also define the discrete filter matrix such that . Then, by removing the boundary data beyond the image patch, we can construct the following matrix equation:
where is the vectorisation of the matrix and the overline is the order reversal operation. Moreover, the 2-D Hankel matrix in (7) is defined by
and the 1-D Hankel matrix for the -th column of the matrix is given by
If the underlying signal has -sparse spectral components, we can further show the following key result :
which implies that as long as the annihilating filter size is bigger than the sparsity level, the rank of the associated Hankel matrix is always equal to the sparsity level. This implies the following duality:
By exploiting this duality, our previous work  derived an annihilating filter-based low rank Hankel matrix approach for image inpainting. In this paper, this approach will be recasted into sparse+low rank decomposition for impulse noise removal.
Unlike the lifting scheme used for phase retrieval problems , our lifting scheme to the Hankel structured matrix is linear, so the sparse impulse noise is also lifted to the sparse outliers in the lifted Hankel structured matrix (see Figure 2). Accordingly, if an underlying image is corrupted with sparse impulse noises, then we have
where denotes the low-rank component and represents the sparse components originated from impulse noises, which are both in the Hankel structure. This is the key property we want to exploit in our impulse noise removal algorithm. In fact, to address this type of sparse + low rank decomposition, the robust principal component analysis (RPCA) was actively investigated [26, 27, 28, 29, 30]. More specifically, for a given measurement matrix , the RPCA solves the following minimization problem:
where is the nuclear norm. To minimize this, alternating directions methods  were employed.
Compared to the standard RPCA approach, our sparse+ low rank decomposition problem using (10) requires additional constraint to impose Hankel structures. Therefore the RPCA algorithm should be modified. The specific optimization algorithm under this constraint will be explained in the following sections. Additionally, because the image statistics change across an image with spatially varying annihilating property, a noisy image should be partitioned into overlapped patches, which are processed patch-by-patch using robust ALOHA and the average values are used as described in the algorithm flowchart in Fig. 3.
Note that the Hankel structure matrix in (8) is determined by the underlying image patch () size and the associated annihilating filter () size. For given image patch and annihilating filter, we now denote the associated spaces for the Hankel matrix as . Then, for a given noisy image patch and annihilating filter size, our impulse noise removal algorithm can be implemented by solving the following sparse + low rank decomposition under the Hankel structure matrix constraint:
Since the sparse components in image patch are also sparse in a lifted Hankel structure, () can be further simplified as
where, with a slight abuse of notation, denotes an appropriately scaled version from in . Note that is now in image patch domain, unlike in the lifted Hankel matrix structured matrix domain in . The advantage of over is an associated simpler optimization method. More specifically, if we apply a factorized form of nuclear norm relaxation , then the final problem formulation of the optimization problem can be expressed as
Then, each subproblem is simply obtained from (16). More specifically, we have
It is easy to show that the first step can be simply reduced to a single soft-thresholding in image patch domain rather than in a lifted Hankel matrix space:
where denotes the pixel by pixel soft-thresholding with the threshold value . The simple thresholding step in (23) is the main motivation which is why we prefer over . Now, the second step becomes
where corresponds to the Penrose-Moore pseduo-inverse mapping from our block Hankel structure to a patch, which is calculated as
Note that the adjoint operator adds multiple elements of and put it back to the patch coordinate, and denotes the division by the number of multiple correspondence; hence, the role of the pseudo-inverse is taking the average value and put it back to the patch coordinate. Next, the subproblems for and can be easily calculated by taking the derivative with respect to each matrix. For example, the derivative of cost function for is given by
and the closed-form solution of the subproblem for is obtained by setting . In the similar way, the derivative with respect to can be obtained. Accordingly, the closed-form update equations for and are given by
Even though the original Hankel matrix has large dimension, it is important to note that our algorithm using (26) and (27) only require the matrix inversion of matrix, where denotes the estimated rank of the Hankel matrix. This significantly reduces the overall computational complexity.
Before we apply ADMM, the initial estimatel and have to be determined with an estimated rank. For this, we employed an SVD-free algorithm called low-rank factorization model (LMaFit) . More specifically, for a low-rank matrix , LMaFit solves the following optimization problem:
and is initialized with . LMaFit solves a linear equation with respect to and to find their updates and relaxes the updates by taking the average between the previous iterations. Moreover, the rank update can be done automatically by detecting the abrupt changes of diagonal elements of QR factorization . Even though the problem (28) is non-convex due to the multiplication of and , the convergence of LMaFit to a stationary point was analyzed in detail . However, the LMaFit alone cannot recover the block Hankel structure, which is the reason we use ADMM later to impose the structure.
In many applications, images are obtained through multiple measurement channels. For example, in a colour image, multiple images are measured throughout R (red), G (green) and B (blue) detectors. In multispectral imaging for remote sensing applications, a scene is measured through many spectral bands. In these applications, the underlying structure is identical so that there exists strong correlation between different channel measurements. Regarding the random impulse noise contamination, we may encounter two different scenarios: 1) noisy pixel locations are independent between the channels, and 2) noisy pixel locations are same across the channels. The first scenario is commonly observed when independent detectors are used for each channel. On the other hand, when a spectrometer is used to split an input into multiple channels, then the noisy pixel locations should be common across the channel. Therefore, in this section, we are interested in extending single-channel robust ALOHA to address these two cases.
Let denote an underlying image patch that are common for all channel measurements, and be its spectrum. Then, as shown in Fig. 4(a), the spectrum of the -th channel measurement can be modeled as:
where denotes a spectral modulation function of the -th channel. The model (29) assumes that each channel measurements still retains the textures of the underlying images with a channel specific modulation, which property was previously extensively used in multichannel deconvolution problems [47, 48, 49] illustrated in Fig. 4(b). As a result, it is easy to derive the following inter-channel annihilating filter relation:
To exploit the inter-channel annihilating property in our robust ALOHA, we construct the following matrix:
where denotes the Hankel structured matrix constructed from the -th channel measurement . Then, it is easy to show that
where is defined recursively as follows:
denotes the reversed ordered, vectored spectral modulation filter for the-th channel. Because we have
when the underlying spectrum is -sparse. Hence, by choosing the annihilating filter size sufficiently large, we can make low-ranked and the aforementioned sparse+ low rank model can be used for impulse noise removal.
When the locations of the impulse noise are independent between channels, then the associate optimization problem is very similar to that of the single channel problem. More specifically, if we define
such that the each channel measurement is give by
then the optimization problem becomes
When the noisy pixel locations are common across the channels, we need a major algorithmic change that comes from a common sparsity inducing matrix norm penalty. More specifically, a common support condition of sparse components should be imposed across channels. In this paper, this constraint is formulated using the group-wise mixed norm:
Then, Eq. (16) can be converted to
where, for in (38), is defined as
We first performed denoising experiments using randomly distributed random valued impulse noise (RVIN) that corrupts 25% and 40% of the whole image pixels. The RVIN is given as follows. Let and be the original pixel value at location and the contaminated pixel with impulse noise at location , respectively. When the dynamic range of pixel value is given as , RVIN is described as
where is a random number between
chosen by the uniform random probability density function andis the proportion of noisy pixels with respect to total pixels.
The test sets were Baboon, Barbara, Boat, Cameraman, House, Lena and Peppers images. All test images were rescaled to have values between 0 and 1. For comparison, a median filter method (MATLAB built-in function ‘medfilt2’, indicated as MF in the figures) was used as the simplest reference algorithm, and the existing algorithms such as ACWMF, and TVL1 were also used. The original codes from the original authors were used and the parameters for the comparative algorithms were optimized to have the best performances. The parameters for the proposed method are given in Table. I. The maximum iteration number of ADMM in Eq. (16) was set to , and the stopping criteria was defined as in  with the tolerance set to . For quantitative evaluation, we used the PSNR (peak signal-to-noise ratio). Specifically, when the reference signal () is given, the PSNR of the reconstructed image () is calculated as
We summarized the PSNR results in Table II for all reconstructed images. The proposed method outperformed all other algorithms in both PSNR and visual quality, and the PSNR improvement was up to 8dB. Fig. 5 is the typical reconstruction result that shows noticeable enhancement in both visual quality and quantitative measures. In order to show that the proposed sparse+ low rank decomposition properly decomposes the impulse noises from the images, the decomposed sparse and low-rank components are illustrated in Fig. 6. We can see that sparse component () looks similar to the additive impulse noises as indicated by .
The salt and pepper noise is a special case of RVIN, because it has impulse noises with the intensity of the minimum and maximum values of pixel dynamic range. Specifically, the salt and pepper noises are given by
where variables were defined in Eq. (42).
Due to the extremal pixel values, the salt and pepper noise can be exceptionally well detected by adaptive median filter (AMF), so the AMF-based denoising algorithms with edge preserving prior (AM-EPR) have been proposed . If the position of the salt and noise locations are well detected by AMF, our robust ALOHA can be modified accordingly. More specifically, the estimation of the sparse component in (14) and (15) is no more necessary, and the resulting algorithm becomes identical to the ALOHA image inpainting algorithm . We denote this modification using AMF as AM-ALOHA. To demonstrate that our algorithm still outperforms the existing algorithms, we compared our method with median filtering (MF), adaptive median filtering (AMF), AMF-based denoising algorithms with edge preserving prior (AM-EPR) using 25 % salt and pepper noise. The results in Fig. 7 clearly demonstrates that the proposed AM-ALOHA outperforms all other algorithms.
To verify that the proposed method can be easily extended to multichannel images, we conducted experiments with colour RGB images. As discussed before, the noisy pixel location can be either same across channels or independent for each channel. So we conducted experiments under the two different scenarios. Fig. 8 showed the reconstruction result when channel independent impulse noises were added, whereas Fig. 9 corresponds to the scenario when of impulse noises were added at the same locations across RGB channels. The proposed method provided better detailed structures (e.g. bundle of peppers and the edges of peppers) than TVL1 method as shown in Figs. 8-9. Also, the cartoon-like artifacts were significantly reduced in the proposed method. In the inset images, the detail structures of peppers are magnified to demonstrate superior performance of the proposed method over TVL1.
One of the interesting observations from these experiments was that the proposed reconstruction provided a better PSNR for the channel independent impulse noises. This was because the noiseless pixel values from other channels could improve the image inpainting performance of noisy pixel values by exploiting the correlation between the channels.
To verify that a lifting to a Hankel matrix is essential for performance improvement, we applied the standard RPCA as an impulse noise removal algorithm and compared the results. Two types of RPCA were implemented: one using whole images, and the other using image patch of the same size as our robust ALOHA. Note that the standard RPCA uses an image or patches as they are, without reformulating them into a Hankel structured matrix. For RPCA, we used the software packages provided by the original authors in . We chose the parameters for the best PSNR results in each reconstruction. As you can see in Fig. 10, two forms of RPCA implementations completely failed in removing impulse noises, and the detailed image structures were distorted. On the other hand, the robust ALOHA provided nearly perfect noise removal. Such a remarkable performance improvement was originated from an image modeling using a low rankness of annihilating filter-based Hankel matrix, which again confirm the robust ALOHA is a superior image model and denoising algorithm for images corrupted with impulse noises.
Recall that the proposed robust ALOHA was performed patch by patch without considering additional similar patches. This is of great importance in contrast to other denoising algorithms that use low-rank approaches [16, 52, 53, 54]. While the authors in [16, 52, 53, 54] used the patch-based low-rankness, all methods required additional redundancies from, for example, multiple dynamic frames [16, 52] or groups of similar spectral patches [53, 54]. Even though such additional redundant information may introduce a low-rankness, those approaches could not properly perform denoising without utilising such additional redundancies. On the other hand, the robust ALOHA exploits the low-rankness originated from intrinsic spectral sparsity of an image patch, so no additional redundancy is necessary. Thus, it is more flexible and powerful.
As briefly discussed in the introduction, a recent research  successfully demonstrated accurate prediction of target locations under occlusion using sparse low rank decomposition of Hankel structured matrix. However, unlike our robust ALOHA, one-dimensional trajectories extracted from video sequences are required as inputs to construct the Hankel structured matrix, because the algorithm was derived based on the assumption that those trajectories follow the linear time invariant state-space models, which was first suggested in . On the other hand, the Hankel structured matrix in robust ALOHA is derived from two dimensional patches by exploiting the spectral domain sparsity; thus, the construction of the Hankel matrix is different from . Moreover, we exploit an SVD-free minimization algorithm  instead of an augmented Lagrangian method (ALM) in [31, 55] for saving computational burdens. Therefore, we believe that there exist significant differences between the two approaches.
From the sampling theory point of view, a recent paper by Chen and Chi  provides a theoretical estimate of fundamental performance of our robust ALOHA. In , the authors showed that the required number of samples, , for the recovery of signal with corrupted with sparse outlier using (P) in (III-A) is given by
when the regularization parameter in (III-A) is given by and the a noise corruption fraction is smaller than . In (44), is a numerical constant depending on the corruption fraction, is incoherence parameter, , and is the rank of single channel Hankel matrix. Because in our robust ALOHA for impulse noise removal, the theoretical result in (44) strongly supports our finding that as long as the spectrum of a patch is sufficiently sparse, the proposed sparse and low-rank method can restore the corrupted signal even with significant impulse noises (in our simulation, up to 40 ).
In this paper, we proposed a sparse + low rank decomposition of annihilating filter-based Hankel matrices for impulse noise removal. The new algorithm, called robust ALOHA, extends the conventional RPCA approaches by exploiting the spectral domain sparsity and the associated rank deficient Hankel matrix. The robust ALOHA was implemented using ADMM iteration with initialization using LMaFit algorithms. In our ADMM formulation, factorization based nuclear norm minimization, was used instead of SVD so that the computational gain is achieved. We demonstrated that the robust ALOHA outperformed the existing impulse noise removal algorithms up to 8dB. Furthermore, we showed that the robust ALOHA can be used for salt and pepper noises by incorporating the estimated noise locations. In addition, the extension to impulse noise removals from colour channels was very straightforward by concatenating the Hankel structure matrix side by side and imposing the low-rankness.
The superior performance of the robust ALOHA as well as ALOHA inpainting  clearly shows that image modeling using annihilating filter based Hankel matrix is a very powerful tool with many image processing applications.
The authors would like to thank Dr. Nikolova for sharing source codes of AM-EPR. This work was supported by Korea Science and Engineering Foundation under Grant NRF-2014R1A2A1A11052491.
J. Wright, A. Y. Yang, A. Ganesh, S. S. Sastry, and Y. Ma, “Robust face recognition via sparse representation,”IEEE Trans. Pattern Anal. Mach. Intell., vol. 31, no. 2, pp. 210–227, 2009.
W. Dong, G. Shi, and X. Li, “Nonlocal image restoration with bilateral variance estimation: a low-rank approach,”IEEE Trans. Image Process., vol. 22, no. 2, pp. 700–711, 2013.