NOISE FILTERING OF REMOTELY SENSED IMAGES USING ITERATIVE THRESHOLDING OF WAVELET AND CURVELET TRANSFORMS

In recent past, advancements in observing the Earth from space have led to a new class of images with very high spatial resolution. These high resolution images give good quality of detailed information about the properties of the various objects. For different purposes, remote sensing images are used to extract some features, detect the presence various phenomena, or for interpretation. These applications require high signal-to-noise ratio (SNR) to get correct results and better performance. The data that are contaminated with noise can cause a failure to extract valuable information and degrade the interpretation. This article presents some techniques for noise filtering of remotely sensed images based on multiresolution analysis (MRA). Multiresolution techniques provide a coarse-to-fine and scale-invariant decomposition of images for image interpretation. The multiresolution image analysis methods have the ability to analyze the image in an adaptive manner, capturing local information as well as global information. Further, noise being one of the biggest problems in image analysis and interpretation for further processing, is effectively handled by multiresolution methods. The paper aims at the analysis of denoising of image using wavelets and curvelets on high resolution multispectral images acquired by the Quickbird and medium resolution Landsat Thematic Mapper satellite systems. Wavelet transform showed great effect when dealing with one and two-dimensional signal with point singularity features. However, for the two-dimensional image, the main characteristics were characterized by the edges. Wavelets can only capture limited directional information due to its poor orientation selectivity. By decomposing the image into a series of high-pass and low-pass filter bands, the wavelet transform extracts directional details that capture horizontal, vertical, and diagonal details. However, these three linear directions are limiting and might not capture enough directional information in remotely sensed images. Wavelet transform coefficients are not the best and the most sparse to describe the image edge singular features. In order to avoid this shortcoming of wavelet transform and process images of high dimension more effectively, curvelet transforms present in the literature is used. It combines the directional filtering and multiscale ridgelet transform, which makes curvelets useful for denoising of images. Its anisotropic characteristic is advantageous to the edge expression, especially to the curve singularities of signal. To improve the performance of noise filtering an iterative thresholding scheme for wavelets and curvelets is proposed to the problem of restoring an image from noisy image and analyzing the effects of denoising. For each region of the MRA the algorithm is tested to obtain optimum results. The problem is to decide the threshold value to be chosen for noise removal. Each coefficient of the image is considered as a threshold value and the final PSNR of the reconstructed image is calculated. On comparing of the PSNR values obtained, the threshold value giving the highest PSNR is chosen. Two comparative measures are used for evaluation of the performance of the methods for denoising. One of them is the peak signal to noise ratio and the second is the ability of the noise filtering scheme to preserve the sharpness of the edges. By both of these comparative measures, the curvelet has proved to be better than the others. Results are illustrated using Quickbird and Landsat images for fixed and iterative thresholding. Nosiy Image (SNR= 9.54dB) Wavelet with fixed thresholding (PSNR =25.92dB) Wavelet with proposed thresholding (PSNR =27.7dB) Curvelet denoising (PSNR=28.69) ISPRS Technical Commission I Symposium, Sustaining Land Imaging: UAVs to Satellites 17 – 20 November 2014, Denver, Colorado, USA, MTSTC1-103


INTRODUCTION
In recent past, advancements in observing the Earth from space have led to a new class of images with very high spatial resolution.These high resolution images contain detailed information about the properties of various objects.For different purposes, remote sensing images are used to extract some features, detect the presence of various phenomena, and for interpretation.These applications require high signal-tonoise ratio (SNR) to get correct results and better performance.The data that are contaminated with noise can degrade the interpretation and hamper extraction of useful information.
A large number of noise filtering algorithms are present in the literature.There are a number of benchmark non-linear filters like median filters, Wiener filters, Min/Max filters; adaptive filters like Lee filters (Lee, 1981), Frost filters (Frost et al., 1980) and sigma filters, along with linear filters like low-pass and high-pass filters.While many techniques are present for denoising, wavelet denoising is the technique used due to its added advantage of localization in frequency as well as spatial domain and multi-scale, multi-resolution analysis (Chen and Bui, 2003;Chen and Kegl, 2007;Dabov et al. 2007).Amongst the numerous wavelets available to use, the most commonly used are the family of Daubechies wavelets (Daubechies, 1992), the bi-orthogonal wavelets.
Wavelet transform showed great effect when dealing with one and two-dimensional signals with point singularity features.However, for the two-dimensional image, the main characteristics were characterized by the edges.Wavelets can only capture limited directional information due to its poor orientation selectivity (Cand`es and Donoho, 2000).By decomposing the image into a series of high-pass and low-pass filter bands, the wavelet transform extracts directional details that capture horizontal, vertical, and diagonal details.However, these three linear directions are limiting and might not capture enough directional information in remotely sensed images.Wavelet transform coefficients are not the best and the most sparse to describe the image edge singular features.In order to avoid this shortcoming of wavelet transform and process images of high dimension more effectively, Cand`es and Donoho (2000) introduced curvelet transform, combining the directional filtering and multi-scale ridgelet transform (Cand`es and Donoho, 1999), which makes curvelets useful for noise filtering of images.Its anisotropic characteristic is advantageous to the edge expression, especially to the curve singularities of two dimensional signals.
In this paper, first some basic ideas of wavelet, ridgelet and curvelet transforms are discussed.Next, a model noise filtering problem is considered for remotely sensed imagery with Gaussian noise, and existing methods of thresholding followed by proposed method of iterative thresholds for denoising are described.Finally the results with various performance measures are discussed with interpretations and direction of future work is presented.

Wavelet Transforms
The application of multi-resolution analysis (MRA) for image analysis and interpretation has become very popular in recent past.A multi-resolution technique can provide a coarse-to-fine and scale-invariant decomposition for image interpretation.Therefore it is effective to analyze an image starting from the global view of the coarse resolution and then gradually increasing the resolution to a more local view.
According to algorithm in (Mallat, 1989), the input image is convolved with low-pass and high-pass filters associated with a mother wavelet and down-sampled.Four images (each one with half the size of the original image) are obtained, corresponding to high frequencies in the horizontal direction and low frequencies in the vertical direction (HL), low frequencies in the horizontal direction and high frequencies in the vertical direction (LH), high frequencies in both directions (HH) and low frequencies in both directions (LL).The last image (LL) is a low-pass version of the original image, and is called as the approximation image.This procedure is repeated for the approximation image at each resolution 2j (for dyadic analysis and synthesis).If the wavelet transform is applied up to the scale 2j, the original image can be reconstructed using the above direction images as shown in Fig. 1.Also, since a lowpass filter is involved, noise suppression is implicit to this approach.

Ridgelet Transforms
Wavelets are suitable for dealing with objects with point singularities.Wavelets can only capture limited directional information due to its poor orientation selectivity.By decomposing the image into a series of high-pass and low-pass filter bands, the wavelet transform extracts directional details that capture horizontal, vertical, and diagonal details.However, these three linear directions are limiting and might not capture enough directional information in remotely sensed images.
Wavelets do not isolate the smoothness along the edges, and thus more suitable for reconstruction of sharp point singularities than lines or edges.These shortcomings are addressed by ridgelet transform introduced in (Cand`es and Donoho, 1999), where a line singularity is mapped into a point singularity using Radon transform.Then, the wavelet transform is used to handle the point singularity.The result is an efficient representation for 2-D functions with piecewise smooth regions separated by a line singularity.Ridgelet analysis may be considered as wavelet analysis in the Radon domain because Radon transform translates singularities along lines into point singularities, for which the wavelet transform is known to provide a sparse representation (Fadili and Starck, 2007) The algorithm of the Discrete Ridgelet Transform (DRT) is depicted in Fig. 2. The DRT of an image of size n × n is an image of size 2n × 2n, introducing a redundancy factor equal to 4 (Welland, 2003).

Curvelet Transforms
The ridgelet transform is optimal at representing straight-line singularities.This transform with arbitrary directional selectivity provides a key to the analysis of higher dimensional singularities.But, the ridgelet transform is only applicable to objects with global straight-line singularities, which are rarely observed in remotely sensed images (Woodcock and Strahler 1987).In order to analyze local line or curve singularities, a natural idea is to consider a partition for the image, and then to apply the ridgelet transform to the obtained sub-images.This block ridgelet based transform, which is named curvelet transform, was first proposed in (Cand`es and Donoho, 2000).Curvelet basis functions can be viewed as a local grouping of wavelet basis functions into linear structures so that they can capture the smooth discontinuity curve more efficiently as illustrated in Fig. 3. Curvelets partition the frequency plane into dyadic coronae and (unlike wavelets) sub-partition those into angular wedges which again display the parabolic aspect ratio.Hence, the curvelet transform refines the scale-space viewpoint by adding an extra factor, orientation, and operates by measuring information about an object at specified scales and locations but only along specified orientations.The curvelet transform has gone through two major revisions (Donoho and Duncan, 2000).The first version (Cand´es and Donoho, 2000) used a complex series of steps involving the ridgelet analysis of the radon transform of an image.Their performance was very slow; therefore, an improved version was developed known as Fast Discrete Curvelet Transform (FDCT) (Cand`es et al., 2006).In this new method, the use of the ridgelet transform as a preprocessing step of curvelet was discarded, thus reducing the amount of redundancy in the transform and increasing the speed considerably.
According to Cand`es et al. (2006), two implementations of FDCT are proposed: (i) Unequally Spaced Fast Fourier transforms (USFFT) (ii) Wrapping Function Both implementations of FDCT differ mainly by the choice of spatial grid that is used to translate curvelets at each scale and angle.Both digital transformations give a table of digital curvelet coefficients indexed by a scale parameter, an orientation parameter, and a spatial location parameter.For numerical computations we need to discretize the continuous curvelet transform, since we usually work with discrete data sets in applications.As shown by Cand`es and Donoho (2003), a discrete version of the continuous curvelet transform can be derived by a suitable sampling at the range of scales, orientations and locations.In this paper, an implementation of this algorithm with wrapping function is used.

Fast Discrete Curvelet Transform via Wrapping
The new implementation of curvelet transform based on Wrapping of Fourier samples takes a 2D image as an input in the form of a Cartesian array f [m, n], where 0 ≤ m<M, 0 ≤ n<N where M and N are the dimensions of the array.The following are the steps of applying wrapping based FDCT algorithm (Cand`es et al. 2006) as shown in Fig. 4. 1) Apply the 2D FFT to an image f to obtain Fourier samples F[m, n] 2) For each scale j and angle l, form the product The support of Uj is contained the rectangle Rj 3) Wrap this product around the origin and calculate ~~ where the range for m, n and θ is 0 ≤ m< 2 j , 0 ≤ n< 2 j/2 , and −π/4 ≤ θ<π/4.4) Apply IFFT to the processed curvelet coefficients Fig. 4: Construction of Fast Discrete Curvelet Transform (Cand`es et al., 2006) Discrete curvelet transform in the spectral domain utilizes the advantages of FFT.During FFT, both image and curvelet at a given scale and orientation are transformed into the Fourier domain.The convolution of the curvelet with the image in the spatial domain then becomes their product in the Fourier domain.After this step, a set of curvelet coefficients are obtained by applying IFFT to the spectral product.This set contains curvelet coefficients in ascending order of the scales and orientations.

NOISE FILTERING
Noise represents the unwanted information which deteriorates the quality of the image.A remotely sensed optical image is mainly corrupted by noise which is modelled as additive Gaussian noise.Multiscale transform based noise filtering is achieved by hard-thresholding of the transform coefficients.

Hard thresholding
The denoising method with hard thresholding is as follows (Starck, Candès, and Donoho 2002);  Estimate the noise standard deviation σ in the input image. Calculate the wavelet and curvelet transforms of the input image.Get a set of bands w j , each band contains N j coefficients and corresponds to a given resolution level. Calculate the noise standard deviation σ j for each band j of the transformed coefficients. For each band j: Calculate the maximum of the band and multiply each coefficient. Reconstruct the image from the modified coefficients.
As described in (Starck, Cand`es and Donoho 2002), threshold at 3σ jl for all levels and 4σ jl at finest scale is selected.σ jl is the noise level of a coefficient at scale j and angle l (equal to the noise level times the L 2 norm of the corresponding transform).

Proposed method of iterative threshold
The fixed thresholding (3σ or 4σ) method does not give the maximum Signal to Noise Ratio (SNR) while donoising, to obtain the optimum results in terms of SNR and Mean Squared Error (MSE), an iterative thresholding method is proposed.
The problem is to decide the threshold value to be chosen for noise removal at each scale.Each coefficient of the transformed image is considered as a threshold value and the final SNR of the reconstructed image is calculated by performing hard thresholding on noisy image in the transform domain followed by respective reconstruction.
On comparing of the SNR values obtained, the threshold value giving the highest SNR is chosen.Therefore at a given resolution (scale) of the transformed image, a unique set of threshold values is created from which the optimum threshold is selected.
The selection of coefficient value as threshold is done at each resolution in different high frequency sub-bands so that for SNR calculations, coefficients in the particular detailed subband only are considered instead of all the coefficients in all the sub-bands at a given resolution.Therefore for a given sub-band of the wavelet or curvelet transformed image, a unique set of threshold values is created from which the optimum threshold is selected.Since lesser are the number of elements in the unique set, the lesser number of computations are required to successfully denoise the image.Thus, the unique set is formed of only the existing coefficient values present in the specific band of the wavelet or curvelet transformed image.Then hard thresholding is done as Inverse wavelets, curvelets transforms are then taken.This algorithm is repeated for different wavelet families of Daubechies and Bi-orthogonal filters.For selected remotely sensed data, bi-orthogonal filters of wavelets are used in this paper.This algorithm is repeated for high and medium resolution images with different levels of noise.
In this approach computational load increases as the algorithm runs more number of times (in contrast to fixed thresholding) to obtain the thresholds, however it gives the optimum SNR.

RESULTS AND ANALYSIS
To evaluate the performance of these algorithms for noise filtering, sample high resolution Quickbird image (1024x1024) of Powai area of Mumbai city and Landsat Thematic Mapper (TM) image (512x512) of Mumbai city are used.To create the effect of noise of varying amounts, zero mean white Gaussian noise with various standard deviations (σ = 10, 15, 20, 25 and 30) was added to the input images, and input to the denoising filters are developed and then denoised by fixed threshold and optimal threshold by iterative method using wavelet and curvelet transforms.In addition to visual observations, objective measures like Signal to Noise Ratio (SNR), Peak SNR (PSNR) and Mean Squared Error (MSE) are used to evaluate the performance of the filters.
Further, a few edges are interactively selected in the images, and the edge contrast is measured to evaluate the performance of these filters for edge preserving capacity while noise removal.It is expected that as denoising is performed, edges may get blurred, but due to removing the high frequency content related to noise, it is expected that mean difference across the edge should remain the same before and after denoising.This analysis is done for both the images with two different levels of noise.It is found that curvelets outperform the wavelet transform in terms of SNR, PSNR, MSE and edge preservation capability (Figures 5 and 6).
It is observed that for Quickbird image, with σ = 10, curvelet with iterative threshold (CvI) gives SNR of 17dB (3.86 dB higher than wavelet with iterative threshold (WvI), and 1.71 dB higher than curvelet with fixed threshold (CvF)) and MSE of 88 (14 less than wavelet) (Table 1).For Landsat TM image also, curvelet performs better than wavelet with improvement in SNR, PSNR and MSE (Table 2).
It is expected that as denoising is performed, edges get blurred, so the denoising mechanism should be such that while removing the high frequency content related to noise, the actual edge information should be retained, so it is expected that mean difference across the edge should remain the same before and after denoising.Also, to have intra-class variance minimum and interclass variance maximum, mean to standard deviation ratio should be improved after denoising, or in worst case it should remain same as before denoising.This analysis is done for both the images with two different levels of noise.
From the original image, three edges are chosen to measure the edge preservation while noise filtering.For this purpose one very strong edge (edge 1 in Quickbird image (Fig. 5g)), one line edge (edge 2 in Fig. 5g) and a weak edge (edge 3 in Fig. 5g) are considered as testing points for edge preservation.Mean and standard deviation across edges are considered as objective measures of performance for edge preserving capacity while denoising.It is observed from Table 3, for edge 1 in Quickbird image, in curvelet based denoising, standard deviation across the edge is reduced and mean difference remained almost same as expected, which in turn has improved the mean to standard deviation ratio in both the sides of the edge.Similar improvement on both the sides of edges 2 and 3 is observed.
For edge 1, wavelet also gives better result (as compared to original) but curvelet outperforms the wavelets by giving more mean to standard deviation ratio.
When noise level is increased to σ = 30, from Table 4, for edges 1 and 2, it is observed that curvelet gives higher mean to standard deviation ratio than original whereas wavelet degrades the edges.For side 1 of edge 3, again curvelet performs better than wavelet.For side 2, though as compare to original (with 35.61) mean to standard deviation ratio reduces slightly, curvelet with 29.96 is better than wavelet which gives the ratio of 29.5.
For Landsat TM image with σ = 15, it is observed from Table 5, curvelet gives better results as compare to wavelets for edges 1 & 3, and for side 1 of edge 2. But when the noise is increased to σ = 25, it is observed that curvelet gives better results on only side of all the edges.Though on the other side of edges, mean to standard deviation ratio is less as compared to original but more than that of obtained by wavelets, therefore curvelet with proposed thresholding performs better in this case as well.

CONCLUSIONS
From experimental results, it is found that the curvelet transform outperforms the wavelet transform in terms of SNR, PSNR, MSE and the curvelet denoised images appear visually sharper than the wavelet denoised images.The curvelet transform with proposed iterative thersholding provides high SNR and PSNR values and can filter out Gaussian white noise from remotely sensed images more efficiently than the fixed thresholding method.
For high resolution image of Quickbird, curvelet gives promising results in terms of edge preservation even when noise level is increased.For medium resolution Landsat TM image, curvelets preserve more edge information at low noise levels but when noise level is increased edge quality degrades.The proposed iterative thresholding (in both wavelets and curvelets) significantly improves the SNR and PSNR in comparison with fixed thresholding method but marginally improves edge preservation capability.Therefore, curvelets with iterative thresholding brings two fold advantages both in terms of higher SNR and better edge preservation From the experimental results it is observed that the denoising by the curvelet with proposed scheme is more effective than wavelets, and it is more applicable to high resolution images as low and medium resolution images are less prone to sensor noise.As future work, these methods can be combined with adaptive scales of MRA and variance stabilization of thresholds to further improve the performance of noise filtering.

Table 1 :
Performance of Noise Filtering for Quickbird image