A TECHNIQUE FOR CANCELING IMPULSE NOISE IN IMAGES BASED ON COMPRESSIVE SENSING B.S. Adiga, M. Girish Chandra and Swanand Kadhe Innovation Labs, Tata Consultancy Services, Bangalore, India {bs.adiga, m.gchandra and swanand.kadhe}@tcs.com
ABSTRACT This paper presents a novel technique based on Compressive Sensing (CS) for canceling impulse noise in images. The technique is pivoted around exploiting the strong connection between CS and error correction using the complex (or real) field codes. Even though the usage of real field BoseChaudhuri-Hocquenghem (BCH) codes for impulse noise cancellation in images is rather old, bringing the CS framework to address this classical problem in image processing provides a fresh perspective. Specifically, the paper investigates a CS-based product code based on partial Fourier matrices, with the requisite rows chosen based on a Perfect Difference Set (PDS) or consecutively. The decoding algorithms are based on the CS reconstruction and are rather elegant and computationally efficient compared to those considered earlier for real BCH codes. Extensive simulation studies suggest that the novel PDS-based product code is effective in canceling the impulse noise through iterative decoding. Index Terms— Image, Impulse Noise, Compressive Sensing, Complex Field Codes, Pursuit Algorithms 1. INTRODUCTION The Compressed Sensing or Compressive Sensing (CS) has recently emerged as a very powerful field in signal processing. The methodologies of CS enable the acquisition of signals at a rate much smaller than what is commonly prescribed by Shannon-Nyquist [1], [2], [3] if the signals exhibit the property of sparsity or compressibility. A signal or data represented as a vector (say of length N) is referred to as s-sparse if the representation of that vector in a suitable basis or dictionary has at most s nonzero coefficient values. Compressibility is a weaker notion than sparsity, where the coefficients of representation, sorted in decreasing order of magnitude decay with a power law, and thus allowing the signal to be well approximated by, say, s coefficients. In either case, we typically have s << N . An interesting aspect about CS is its close connection to the areas like high-dimensional geometry, sparse approxi-
mation theory, data streaming algorithms and random sampling (see references in [4]). The newer connections are getting established at quick pace and during the process the CS, many a times, is providing a fresh perspective on many of the areas. Also, there are vigorous efforts in exploring the application of CS in many different fields such as data mining, DNA microarrays, astronomy, tomography, digital photography, sensor networks, and A/D converters [5]. Removing or reducing the impulse noise, also known as “salt and pepper” noise, is a classical problem studied in image processing. Various solutions have been proposed over the years, including, applying median filtering, Block Truncation Coding (BTC), “localization and correction strategy” etc (see [6] and the references therein). Another approach suggested in the papers [6] and [7] (also see [8]) is based on the Bose-Chaudhuri-Hocquenghem (BCH) codes in the field of real numbers. That is, the arithmetic in the latter codes is over real (or complex fields, in general) rather than the usual operation of BCH codes over finite fields or Galois fields. The idea is to consider impulses as errors and the “signature” of this error is available in the appropriate syndrome in the frequency domain (see Section 2 for little more elaboration). Using this signature, the decoder suggested in these references first evaluate the number of errors, then locate the errors, and finally find the values of errors towards canceling (subtracting out) the noise. It is well established that the error correction operating in the complex field is strongly connected to CS (see [4], [9] and the references therein). In the light of this link it is logical to think of replacing the real BCH codes of [6], [7] and [8] by encoding and decoding based on CS. Towards this end, in this paper, we propose a technique based on CS for the impulse noise cancellation in images. Following our earlier work [4], the codes considered are based on the Inverse Discrete Fourier Transform (DFT) matrix (of size N × N ). Specifically, the parity check matrices are the partial IDFT matrices with N columns and M rows (with M < N ), where we choose the requisite rows either consecutively or based on a Cyclic Difference Set (CDS) (see Section 2.1.1 for details). The usage of this code can facili-
tate elegant and simple decoding pivoted on sparse-data recovery algorithms of CS, compared to the complicated decoding suggested in [6], [7] and [8]. To the best of our knowledge, the proposed technique is not considered elsewhere. The rest of the paper is organized as follows: Section 2 attempts to capture necessary background information very briefly. The CS and error correction link, partial Fourier codes with CDS based row selection (of [4]) and the real BCH codes (of [6], [7] and [8]) are touched upon in this section. Section 3 captures the proposed technique. Some results and related discussion are provided in Section 4. Section 5 concludes the paper together with some further research issues. 2. JUST ENOUGH INFORMATION FROM THE EXISTNG WORKS 2.1. CS and Error Correction In any coding scheme, we start with an uncoded vector u (dimension K × 1 ) belonging to C K ( C being complex field) and obtain the encoded vector c ∈ C N using the Generator matrix G : c = Gu
(1)
where, G is of dimension N × K with N > K (thus introducing redundancy). Further, G has to be full rank to recover u from c. Now suppose that c is corrupted by an arbitrary vector e of dimension N × 1 : y = c + e = Gu + e
(2)
The task is to reconstruct u from y, for which it is sufficient to reconstruct e , since y − e = Gu , from which u can be recovered due to full rank property of G. By using a matrix H whose kernel is in the range of G, it is possible to recast Eqn.(2) into the form: ~y = H * y = H * (G u + e ) = H * e
(3)
since H ∗G = 0 , and H ∗ is the complex conjugate of the Parity Check matrix H (with dimension M × N , where M = N − K ). ~y is the well known syndrome vector. Since any error correction scheme can tolerate errors up to a certain extent, there is an implicit assumption that the vector e is sparse (say, with number of non-zero entries is ≤ s ) for successful operation. Here, we are essentially reconstructing the sparse vector e* from M linear “measurements” ~y , which are obtained through appropriate projections using the “measurement or sensing” matrix H . Thus, we end up with the same basic problem of CS. As suggested in the CS theory, if the matrix H satisfy certain conditions based on one of these desirable properties: spark, Null-Space Property
(NSP), Restricted Isometry Property (RIP) or Coherence (see [4] and some classical references therein), the recovery of the sparse vector e can be carried out using l1 -norm minimization (also called as basis pursuit) or greedy algorithms such as various matching pursuits.
2.1.1. CDS Based Partial Fourier Parity Check Matrices A set of M distinct integers D = { d 0 , d1 , " , d M −1} is called an (N , M , λ ) Cyclic Difference Set (CDS) if a difference τ ≡ (d i − d j ) mod N takes every integer in 1 ≤ τ ≤ N − 1 exactly λ times. In [4], H is constructed by choosing a subset of rows of the N-point Inverse DFT (IDFT) matrix indexed by a CDS D with the condition that 0 ∉ D . Thus, the elements of H are given by (for 0 ≤ k ≤ M − 1 and 0 ≤ n ≤ N − 1 ):
[H ]k ,n =
1 ⎛ 2π n d k ⎞ exp⎜ j ⎟ N ⎠ M ⎝
(4)
This deterministic and structured construction can facilitate faster reconstruction of e with lower complexity; see [4] and references therein for more details on one such algorithm. Since G is also structured in this case, the construction also enables efficient encoding. This matrix H is guaranteed to recover s-sparse error by l1 - minimization or greedy algorithms as long as 1⎛ M ⎞ s < ⎜⎜ + 1⎟⎟ 2⎝ M −λ ⎠
(5)
It is useful to note the following points: (1) when λ = 1 , a CDS is referred to as a Perfect Difference Set (PDS); (2) from Eqn.5, a larger value of M results in higher capability to correct errors. 2.2. Real BCH Codes for Impulse Noise Cancellation The definition and properties of real BCH codes can be investigated using the DFT ([6],[7] and [8]). The encoding can be carried out by first computing the length K DFT of the real-valued u (the message vector) and the resulting block is then padded with N-K consecutive zeros in such a way that Hermitian symmetry is preserved. Essentially, we are putting zeros in the high-frequency indices of the DFT. The codeword c is then generated by taking the IDFT of the zero-padded block. Using the Maximum Distance Separable property of these codes, we have the relationship N − K = 2t , where t is the number of errors the code can correct. Viewed differently, the DFT of a BCH codeword, i.e., its spectrum has 2t consecutive zeros. If the codeword is corrupted, say, by impulse noise then these 2t consecutive spectral values are not zero. They constitute the syndrome
and can be used as a signature of the impulses to be removed. The operations on this signature towards canceling the impulse noise are rather complex as remarked earlier (see [6], [7] and [8] for more details). Further, towards having more powerful codes, product codes or more explicitly product BCH codes on the real field are suggested in the latter references. The idea is to decode iteratively, first by operating on the rows of the encoded image, and then on the columns. The process is repeated on row and column directions. The rationale being, even when only a few impulses are suitably corrected in the beginning, such correction greatly reduces the task of the following step (in the other direction). 3. CS BASED IMPULSE NOISE CANCELLATION In our proposed technique, we encode the image using CSbased real product code, where each component code in the row and column direction is based on partial Fourier matrix. We have considered parity check matrices obtained through CDS/PDS as mentioned in the previous section as well those obtained by choosing M consecutive rows. The cardinality of the CDS/PDS chosen (same as M) depends on the errorcorrection capability we are looking for each of the component codes, like s = 1 or s = 2 etc. Since we start with parity-check matrix H in our CS framework, for encoding we have to obtain the generator matrix G such that H defines its null space. Due to the great structure of the IDFT matrix, G can be simply obtained by choosing the rows not used by the H matrix (as in [4]). But, this G is not in the systematic form and it is necessary to reduce G to the systematic form since we are adopting the product code. In a product code, we should see the valid codewords in both row and column directions, and this can be easily negotiated by having G in the systematic format and applying G in the two directions one after the other. By using the MATLAB, the G obtained through simple row selection can be put into the systematic form by simply using the reduced row echelon form function “rref”! Essentially, to encode, we take blocks of K pixel values row-wise first and multiply by G matrix to arrive at the encoded blocks. Then, blocks of K values in the column direction (of the row-wise encoded image) are coded by G to arrive at the final product-coded image. For decoding, one can use different CS based algorithmms, like, basis pursuit, Compressive Sampling Matching Pursuit (CoSaMP) [9], or the reduced-complexity algorithm of [4] which exploits the CDS structure. The last two algorithms are computationally efficient, more elegant and straightforward than the decoding algorithm of [6], [7] and [8]. It is useful to note that all these decoding algorithms can indicate the “decoding failure” when there is no convergence (defined differently for different algorithms). This feature can be beneficially exploited towards providing the better image quality, say, in terms of Peak Signal-to-Noise
ratio (PSNR) at the user end, by suitably handling the undecoded blocks. 4. SOME RESULTS AND DISCUSSION Extensive simulations are carried out with the standard 256×256 grayscale Lenna image. Two candidate codes are considered towards the results presented in this paper. The first code corresponds to the 4 × 8 parity check matrix where 4 consecutive rows of the 8 × 8 IDFT matrix are chosen. We refer to this code as Code1 in this paper and the rate of this code is 0.5 (needless to mention the generator matrix of this code is also of dimension 4 × 8 ). The second code, referred to as Code2, corresponds to the 8 × 57 parity check matrix whose rows are chosen using the PDS with parameters (N , M , λ ) = (57, 8,1) . The generator matrix of Code2 is of dimension 49 × 57 and the rate is thus 0.8596. Further, the Code2 is guaranteed to correct 2 errors ( s = 2 from Eqn.5); of course, it can correct some patterns of more than 2 errors. The generator matrices of these two codes are put into systematic format and are applied in rows direction first followed by column directions, resulting in the encoding of the image by CS-based real product codes. The encoded image is corrupted by different percentages of impulse noise like 3%, 4%, 5%, etc. As usual, 3% impulse noise means a pixel is damaged with the probability of 0.03. The pixel locations are chosen randomly and the experimentation is carried out with different realizations of these randomly damaged pixels. For both the product codes, decoding is performed iteratively. The decoding is first carried out over the columns of the corrupted image and then over the rows. This complete cycle of decoding, first over the columns and then over the rows, is referred to as one iteration in the discussion of our results. As mentioned earlier, each of these decodings is carried out using CS-based reconstruction algorithms. For both Code1 and Code2, we have considered the CoSaMP algorithm, which is computationally efficient compared to the l1 -norm minimization. It is to be noted that the decoding of Code2 can also be carried out using the reduced complexity algorithm adopted in [4], and the results are not provided here. The simulation study is also carried out with different number of iterations ( nI ) for each of the specific scenarios. It is observed that when the impulse noise level is up to 3%, both Code1 and Code2 can completely remove the noise. Further, three ( nI = 3 ) iterations are sufficient for both the codes with CoSaMP decoding. For 4% and 5% impulse noise cases, it is observed that the PDS-based technique can completely remove the noise for most of the cases in 5 iterations (different cases corresponds to different locations of damaged pixels, as mentioned earlier). See Fig.1 for a typical scenario with Code2. On the other hand, when Code1 is used, the impulse noise cannot be completely removed. However, still the scheme recovers the image fairly
well removing most of the impulses, one reason being its lower code rate of 0.5. Fig.2 shows one typical result for this case after 3 iterations, beyond which the scheme ceases to correct any errors. It can be seen that most of the impulses are removed and the PSNR is 33.2 dB. In essence, we can say that the PDS-based Code2 is a useful candidate with its impulse-noise cancellation ability, higher code rate and computationally efficient decoding.
of the image where decoding failure had occurred. It will be also interesting to investigate different issues (including quantization) for the combination of the proposed CS-based “channel code” and the CS-based source coding of images (that is, when an image is compressed based on CS, either directly by using devices like single-pixel camera or by applying measurement matrix to the conventional images). 6. REFERENCES [1] L. Jacques and P. Vandergheynst, “Compressed Sensing: When sparsity meets Sampling,” Book Chapter in Optical and Digital Image Processing - Fundamentals and Applications, Edited by G. Cristòbal, P. Schelkens and H. Thienpont, Wiley-VCH, April 2011. [2] P. Boufounos, G. Kutyniok and H. Rauhut, “Compressed Sensing for Fusion Frames,” In Proc. SPIE Wavelets XIII, San Diego, Aug. 2009. [3] M.A. Davenport, M.F. Duarte, Y. C. Eldar and G. Kutyniok “Introduction to Compressed Sensing,” Book Chapter in Compressed Sensing: Theory and Applications, Edited by Y. C. Eldar and G. Kutyniok, Cambridge University Press, 2011.
Fig. 1. Typical performance of Code2 for the case of 5% impulse noise; (a)-(d) correspond to iteration number 1-4.
[4] B.S. Adiga, M. Girish Chandra and Shreenivas Sapre, “Guaranteed Error Correction Based on Fourier Compressive Sensing and Projective Geometry,” IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Prague, Czech Republic, pp. 3744 – 3747, May 2011. [5] M. Stojnic, F. Paravaresh and B. Hassibi, “On the Reconstruction of Block-Sparse Signals with an Optimal Number of Measurements,” IEEE Trans. Sig. Proc, Vol.57, No.8, pp. 3075-3085, Aug. 2009. [6] O. Rioul, “A Spectral Algorithm for Removing Salt and Pepper from Images”, IEEE Digital Signal Processing Workshop, Norway, pp. 275-278, Sept. 1996.
Fig. 2. Typical performance of Code1 for 5% impulse noise after 3 iterations; (a)-(c) correspond to iteration number 1-3. 5. CONCLUSIONS The paper presented an effective and computationally elegant scheme based on Compressive Sensing (CS) for impulse noise cancellation in images, comprising of complexfield product codes based on partial Fourier matrix and pursuit decoding algorithms. As a further study, it may be worth considering the removal of residual impulses through other techniques like median filtering by applying the latter only on the portions
[7] A. Gabay, P. Duhamel and O. Rioul, “Real BCH Codes as Joint Source Channel Codes for Satellite Images Coding”, IEEE Global Telecommunications Conference (GLOBECOM), San Francisco, USA, pp. 820-824, Nov. 2000. [8] A. Gabay, P. Duhamel and O. Rioul, “Spectral Interpolation Coder for Impulse Noise Cancellation over a Binary Symmetric Channel”, European Signal Processing Conference (EUSIPCO), Tampere, Finland, Sept. 2000. [9] D Needell and J.A. Tropp, “CoSaMP: Iterative Signal Recovery from Incomplete and Inaccurate Samples”, Appl. Comput. Harmon. Anal., Vol.26, No.3, pp. 301321, 2008.