Modified Curvelet Thresholding Algorithm for Image Denoising

: Problem statement: This study introduced an adaptive thresholding method for removing additive white Gaussian noise from digital images. Approach: Curvelet transform employed in the proposed scheme provides sparse decomposition as compared to the wavelet transform methods which being nongeometrical lack sparsity and fail to show optimal rate of convergence. Results: Different behaviors of curvelet transform maxima of image and noise across different scales allow us to design the threshold operator adaptively. Multiple thresholds depending on the scale and noise variance are calculated to locally suppress the curvelet transform coefficients so that the level of threshold is different at every scale. Conclusion/Recommendations: The proposed algorithm succeeded in providing improved denoising performance to recover the shape of edges and important detailed components. Simulation results proved that the proposed method can obtain a better image estimate than the wavelet based restoration methods.


INTRODUCTION
In the real world signals do not exist without noise, which arises during image acquisition (digitization) and/or transmission (Gonzalez and Woods, 2002). When images are acquired using a camera, light levels and sensor temperature are major factors affecting the amount of noise. During transmission, images are corrupted mainly due to interference in the channel used for transmission. Removing noise from images is an important problem in image processing (Ruggeri and Vidakovic, 1999). This noise removal takes place in the original time space domain or in a transform domain. In transform domain Fourier transform is used in the timefrequency domain and multiresolution transforms like wavelet/curvelet/contourlet transforms are used in the time-scale domain.
Denoising a given noise corrupted signal is a traditional problem in both statistics and in signal processing applications. Linear denoising methods are not so effective when transient non-stationary wideband components are involved since their spectrum is similar to the spectrum of noise (Zhang and Luo, 1999). Nonlinear denoising methods (Smith and Agaian, 2004) rely on the basic idea that the energy of a signal will often be concentrated in a few coefficients in the transform domain while the energy of noise is spread among all coefficients in transform domain. Therefore, the nonlinear methods will tend to keep a few larger coefficients representing the signal while the noise coefficients will tend to reduce to zero. Denoising methods based on multiresolution transforms involves three steps: A linear forward transform, nonlinear thresholding step and a linear inverse transform. Wavelets are successful in representing point discontinuities in one dimension, but less successful in two dimensions. As a new multiscale representation suited for edges and other singularity curves, the curvelet transform has emerged as a powerful tool. The developing theory of curvelets predict that, in recovering images which are smooth away from edges, curvelets obtain smaller asymptotic mean square error of reconstruction than wavelet methods (Candes and Donoho, 2004).
Curvelet transform: An image, when analyzed using a 2-D wavelet transform exhibits large wavelet coefficients, along the edges in the image. At each scale, these edges in the image are seen repeated. This requires many wavelet coefficients to reconstruct the edges in an image and it puts a limit on wavelet denoising. The estimation of these large numbers of coefficients lead to high Mean Square Error (MSE). A comparison of order of error in case of wavelets and that of curvelet reconstruction is presented in the literature (Candes and Donoho, 2004) below.
Non-linear approximation of objects can be considered by treating them as function of two variables, one with discontinuities along edges and second, which are smooth. If an object f is represented in a wavelet basis, then the number of wavelet coefficients of f exceeding the threshold, 1/n, in absolute value increases rapidly as c.n as n→∞. This increase indicates the requirement of many terms to reconstruct a good image. If the best partial reconstruction obtained by selecting the n largest terms is represented by W n f , the best n-term approximation would obey: This result is not optimal and hence wavelets fail to represent objects with edges. If the edge curve is considered to be of length one, then at each scale 2 −j , there are approximately 2 j wavelets interacting with the edge, resulting in coefficient of size 2 −j . This shows that the wavelet coefficient decay only as 1/n. Wavelets do not take advantage of the geometry of the edge and hence about 2 j coefficients are needed to reconstruct the frequency content of an edge up to the subband j 2 ξ ≈ .
Wavelets being nongeometrical, move only in horizontal and vertical directions and not along the curve and thus cannot achieve the optimal rate of convergence.
The curvelet transform has a tight frame which combines multiscale analysis and ideas of geometry and can achieve optimal rate of convergence by simple thresholding (Starck et al., 2002). This multiscale transform has a strong directional character in which elements are anisotropic at fine scales. The support of these elements is according to the parabolic scaling principle length 2 ∼width. Curvelets partition the frequency plane into dyadic coronae. Unlike wavelets, these corona are subpartitioned into angular wedges displaying the parabolic aspect ratio as shown in Fig. 1.
Curvelets at scale 2 −j , are of rapid decay away from a 'ridge' of length 2 −j/2 and width 2 −j and this ridge is the effective support γ µ . If coefficients of only those curvelets, which overlap with the edge and are nearly tangent to it are considered, then for a fixed scale 2 −j , there are at most O(2 −j/2 ) coefficients of such type. The size of each coefficient can be estimated by: Curvelets are L 2 normalized so that 2 L 1 µ γ ≤ and are supported in a box of side length 2 −j/2 and width 2 −j . Therefore they obey: Since, f is a bounded function, the coefficients θ µ , verify the apriori estimate: At each scale 2 −j there are O(2 j/2 ) coefficients which are bounded by C.2 −3j/4 and the nth largest coefficient n θ is bounded by: The above decay gives O(n −2 ) convergence rate for the nonlinear n-term approximation f n defined by keeping the n largest term in the curvelet expansion and it obeys: From the above equation the faster rate of decay of coefficients in a curvelet basis is obvious as compared to the decay of wavelet coefficients in Eq. 1 for an object with arbitrary C 2 singularity.

Preliminary algorithm for image denoising:
The mathematical model of noisy image is as follows: Where: y = The observed image x = The unknown original image n = The contaminating noise Complete curvelet denoising procedure is performed by taking curvelet transform of the image and then applying thresholding to eliminate noisy coefficients. Thus, the inverse curvelet transform of the thresholded coefficients give the denoised image.
The fast discrete curvelet transform (Candes et al., 2006) of the observed image is evaluated as Y using the curvelet transform operator C(.) using following equation: The threshold denoted by λ, for any transform can be expressed in general terms using the operator D(.) as: For wavelet based denoising procedures one such threshold is the universal threshold (Li and He, 2006) given as: Where: σ = The standard deviation of noise N = The size of image Wavelet transform maps white noise in the signal domain to white noise in the transform domain. Thus in the transform domain the signal is concentrated into fewer coefficients but the noise does not concentrate. The principle behind separation of signal and noise is that, when scale 2 −j decreases, wavelet transform maxima of images doesn't increase, but at the same time wavelet transform modulus of white noise increases. Thus different behaviors of wavelet transform maxima of images and noise across different scales allow us to design the operator D(.) adaptively.
An image features a wide variety of characteristics. Hence, instead of using a single value as the global threshold, the operator D(.) can be designed to produce multiple local thresholds λ j adaptively for different scales from fine to coarse. For wavelets an adaptive threshold is (Li and He, 2006) where, j is the decomposition level of wavelet packet transform. This modified multiple local thresholding technique obtained better results than the soft and the hard thresholding methods, which utilizes a single threshold value at every scale. Curvelet transform employs the 1-D wavelet transform as a component step, but, along the radial variable in Radon space. Thus Eq. 11 does not prove to be effective for thresholding the curvelet transform coefficients and requires some modification.

Modified curvelet thresholding for image denoising:
In this study a similar multiple threshold technique for thresholding the curvelet coefficients is proposed. To design the operator D(.), it is proposed to retain all the coefficients at the first scale, since they are the dc values and they provide the average information of the image. For the remaining scales the coefficients, which provide the highest PSNR values seem to be correlated and the curvelet coefficients appear to decay in an exponential manner as shown in Fig. 2. In Fig. 2 Series 1-3 are plots for , j λ corresponding to σ = 20, 30 and 50 respectively. Thus a scale dependent exponential function multiplied by a scale dependent logarithmic function resulted in improvement in PSNR values. Therefore, the multiple local thresholds are proposed as: for j = 2, 3,………J Where: j = The decomposition level of curvelet transform J = The integer corresponding to the last scale After selecting the threshold level, the next step in the process of denoising is applying the threshold operator T(.,.). This can be expressed as: where, Z gives the thresholded coefficients.
This study employs cubic, soft and hard thresholding functions as the threshold operator T(.,.) in Eq. 13, with the difference that the threshold is not single valued λ, but is multivalued λ j as obtained from Eq. 12.
Cubic thresholding function is very flexible and has the capability to adapt to different types of images and threshold operators. The cubic threshold function (Li and He, 2006) given as: The soft (Donoho, 1995) and the hard thresholds (Donoho, 1994) are some simple but powerful shrinkage functions. These thresholds select a single global threshold for all the scales using Eq. 10. Employing , j λ of Eq. 12 with the soft thresholding operator proposes the multiple local soft thresholding as: Similarly, by combining , j λ of Eq. 12 with the hard threshold operator, multiple local hard thresholding can be given by: Finally C −1 (.) takes the fast discrete inverse curvelet transform of the thresholded curvelet coefficients, Z, as: where, x is the reconstructed/denoised image.

RESULTS AND DISCUSSION
The performance of the proposed thresholding methods is evaluated and compared with that of soft, hard and cubic thresholding schemes using wavelets (Li and He, 2006). Gaussian noise was added to the classical Lenna and Saturn images. Multiple local thresholds are obtained using Eq. 12. The curvelet coefficients are processed by thresholding functions in Eq. 14-16. The performance of denoising is evaluated using Peak Signalto-Noise Ratio (PSNR) and Mean Square Error (MSE). PSNR is defined as the ratio of signal power to noise power. It basically obtains the gray value difference between resulting image and original image. MSE is given by the formula:  Table 1 and Table 2 respectively. The curvelet reconstruction using multiple local thresholds enjoys superior performance over the wavelet based reconstructions. The pictorial denoising performance for images 'Lenna' and 'Saturn', using wavelet based soft, hard and cubic thresholds is compared with curvelet based multiple local soft, hard and cubic thresholds in Fig. 3 and 4 respectively. Experiments show that multiple local thresholding based on curvelets outperforms the wavelet based methods on the basis of MSE and PSNR.

CONCLUSION
In this study, a new denoising technique based on adaptive selection of thresholds to suppress noisy curvelet transform coefficients is presented. Due to multiresolutional dictionary, the maxima of the curvelet transform coefficients vary and so the threshold operator is designed to produce as many local threshold values as are the scales. The proposed method efficiently adapts to noise characteristics for different scales and reduces the noise while preserving edges in the image. The thresholding function chosen are the cubic, hard and soft thresholds and the proposed expression is tested against them. From the restored images it can be visually depicted that the edges and texture are well preserved taking the advantage of the fact that curvelets being geometrical very well align themselves to the contours of the edges. Numerical experiments show the good performance of the proposed method in comparison to wavelet based decomposition. Further works involve extending the proposed method to various classes of images which are different from natural images. Another important issue is to test the performance on higher resolution images.