Abstract: The rapid development of Extracorporeal Shock Wave Lithotripsy (ESWL) due to its non-invasive advantage calls for efficient image processing methods like the enhancement for visual inspection or a preprocessing for clinic diagnosis. In this study, we proposed a bivariate local Wiener filtering algorithm based on Undecimated Multiwavelet Transform (UMWT), where the dynamically generated windows are used for estimation of the signal variances of noisy wavelet coefficients. The computational complexity of the algorithm is low, since the operation of UMWT is the main computational burden, slight more than single wavelet transform. The performance of the new filter is assessed with simulated data experiments and tested with actual ESWL images. The results show that the proposed technique can help to meet the conflicting requirements of reproducing the low-contrast details with suppressing the overwhelming noise.
INTRODUCTION
Extracorporeal Shock Wave Lithotripsy (ESWL) is a minimally invasive treatment of kidney stones and biliary calculi through the use of shock waves (Srisubat et al., 2009; Khatami et al., 2008). Nevertheless, there are limitations with the success rate in ESWL due to great noise and artifacts present in the images which brings difficulties to clinical analysis. To improve the visual quality and treatment effect, a tradeoff between noise reduction and the preservation of actual image features has to be made to enhance the diagnostically relevant image content (Debakla et al., 2011).
Noise suppression in ESWL images is a challenging problem so that few studies were dedicated to it. A major approach applied in the current lithotripsy equipments is to take the average of the generated images. The median filter is also commonly employed by some manufacturers. New approaches to image enhancement which will maximize its effectiveness and minimize patient trauma are an open issue, object of further research. Although all information is contained in the medical images, obtaining “feature of interest” in a noisy observation is unfortunately difficult sometimes (Selka and Belbachir, 2005; Mueen et al., 2007; Lotfi et al., 2010). On general electric systems, these images are spatially “warped” to correct for gradient imperfections. In ESWL images with high SNR, the noise can be simulated as a Gaussian distribution but when the SNR is low it has a signal dependent mean and deviates slightly from Gaussian. Accordingly, the use of Gaussian mixture model is discussed in this paper (Yin and Zhang, 2011).
In the past few years, researchers have developed powerful wavelet methods for analysis of medical signals (Kirankumar and Devi, 2007; Idi and Kamarudin, 2012; Alnuaimy et al., 2009). Working in wavelet domain is advantageous since signals can be analyzed at multiscales or resolutions. Wavelet-based image denoising include applying different types of filtering in supposedly smooth and supposedly heterogeneous or “edged” image regions (Foucher et al., 2001; Li and Orchard, 2000), spatially adaptive thresholding (Chang et al., 2000) and locally adaptive Wiener filtering (Mihcak et al., 1999). More recently, the spatial dependency of the wavelet coefficients across adjacent scales in the form of multiscale products was developed (Bao and Zhang, 2003; Murtagh and Starck, 2003). Although the idea is simple, the calculation of interscale correlation is often influenced by the tiny excursion of coefficients along the temporal axes. Pizurica et al. (2003) proposed an empirical estimation of the statistical distributions of the coefficients combining these multiscale products, while using a wavelet domain indicator of the local spatial activity.
These-mentioned techniques were based on different thresholding schemes, most of which estimate the “feature of interest” of a noisy wavelet coefficient in its squared neighborhood window with fixed size. An improved algorithm in the undecimated multiwavelet domain is discussed here, providing an adaptive neighborhood windows’ selection and a reformative MMSE (Minimal Mean Square Error) filter-bivariate Wiener filter.
UNDECIMATED MULTIWAVELET TRANSFORM
The Multiresolution Analysis (MRA) method using the wavelets is well adapted to image analysis and very efficient for computation. Unfortunately, the loss of the translation-invariance nature in the wavelet transform leads to visual artifacts when an image is reconstructed after modification of its coefficients. A way to obtain a translation invariant system is to think about the undecimated transform.
An undecimated transform gives more information about the transformed signal compared to a decimated one. This redundant information is especially important when the wavelet coefficients are used for statistically analyzing the original signals. For the undecimated case, the low pass and high pass filter banks are modified by up-sampling at each consecutive level. This up-sampling is done by inserting zero arrays between each of the filter banks’ coefficient matrices. We can realize the Undecimated Multiwavelet Transform (UMWT) in the following way.
For scalar signals, a preprocessing method of mapping the data to the multiple streams is required. This is done by a prefilter which is defined by a sequence of rxr matrices Q and the data is partitioned into a sequence of r-vectors. Thus, for r = 2:
A postfilter P just plays an opposite part, applied to recover the scalar signals after reconstruction. Obviously, P and Q satisfy the condition PQ = I, where I is the identity filter. The commonly used prefilters are identity prefilter, Xia prefilter, Minimal prefilter, etc.
UMWT using the filter banks {H, G} for a prefiltered signal can be executed at scale j and position k by:
| (1) |
where, Cj and Dj are the lowpass and highpass coefficients at scale j, respectively. We compute the filter matrices H(j)(k) and G(j)(k) = G(k) when l/2j is an integer and H(j)(k) = G(j)(k) = 0, otherwise. The reconstruction is:
| (2) |
To ensure perfect reconstruction, the filter banks
| (3) |
THRESHOLDING METHOD
An observed image that is corrupted by additive white Gaussian noise can be represented in the undecimated multiwavelet domain by:
| (4) |
where, s(i, j) and ε(i, j) are the multiwavelet coefficients of the noise-free image and the noise underlying a multivariate normal distribution, respectively. The signal estimation of each noisy wavelet coefficient is considered by an adaptive neighborhood window around the wavelet coefficients to be thresholded. Here, we introduce a selection mechanism associated with adjacent coefficients to determine the windows. For storage of the variable neighborhood windows we take a dynamic array W to save the location information. The primacy condition for providing a window selection is that the current coefficient is sufficiently large, otherwise W includes the coefficient itself. We define the correlation parameter Corr (u, v) of a neighbor coefficient vector Ď (u, v) around an arbitrary coefficient vector D (u, v) to be thresholded as:
| (5) |
where:
Wi, j is regarded as the selected neighborhood window of a coefficient D(i, j). Corr (·) is used to decide whether a coefficient would be added to Wi, j. It can be deduced that the smaller the Corr (·), the stronger dependence between the pair of coefficient vectors there is. The generation of the neighborhood windows is a recursive process. Given parameters λ and δ, obtaining:
| (6) |
where, the initial value of Wi, j is assigned D(i, j). δ is
an empirical threshold related to an observation of coefficients which are similar
enough to the neighborhood window, as we care about the neighbor coefficients
that contain the most information concerned with the thresholded one.
To improve the running efficiency, a flag matrix μ is used which has the equal size to the image. For a MxN image, μ is initialized by 0(MxN). Within each wavelet subband, if the coefficient D (i, j) has been contained in a neighborhood window, μ (I, j) will be set to 1, so that D (i, j) will be out of estimation for the next window selection. This step can be computed in O(MxN) operations.
Once a neighborhood window is established, local filtering by an improved bivariate Wiener filter is performed. All the coefficients D (p, q) inside a window Wi, j are to be thresholded by the following function:
| (7) |
where, (x)+ = max{0, x}, θp, q = DT(p, q)VTq D(p, q), and Vq is the covariance matrix for the error term according to the scale q:
where, #W denotes the size of W. A robust covariance estimation method is used to estimate Vq directly from the observed coefficients (Huber, 1981).
EXPERIMENTAL RESULTS
Simulated data experiments: In the experiments, we present a simulated imaging experiment that compares the performance of two classical multiwavelets, namely GHM and CL, and those existing techniques in terms of signal-to-noise ratio (SNR), as shown in Table 1. The SNR is computed as:
| (8) |
where, σ2s is the signal variance of the noise-free image, and σ2 is the underlying Gaussian noise power.
| Table 1: | SNR comparison |
| Fig. 1(a-d): | Image quality with different filters (a) original image, (b) median filter, (c) spatial wiener filter and (d) proposed filter (CL) |
It is clear that the proposed thresholding technique gives better results than the other methods. Undecimated CL achieves the highest SNR in most cases.
Actual image estimation: To evaluate the medical image quality, the performance of the filtering algorithms is tested with actual ESWL data. Figure 1a-d displays an example of the original noisy and its denoised images using different filters. As expected, our scheme significantly reduces random fluctuations due to noise without loss of image detail and produces superior contrast in the image compared with other filters.
CONCLUSION
This study proposes an ESWL image enhancement scheme using adaptive bivariate Wiener filtering in undecimated multiwavelet domain. Unlike many traditional schemes that directly threshold the wavelet coefficients within fixed-size windows, the proposed scheme adopts adaptive windows of arbitrary shape to determine the signal variance estimation and better differentiate edge structures from noise. Instead of applying univariate thresholding, we experiment with bivariate Wiener thresholding based on UMWT which has proved to be effective in suppressing artifacts and gives better results than the conventional denoising.
It turns out that the idea of restoring image details and important features in adaptive regions of interest may become a powerful post-processing tool in many ESWL applications. Furthermore, it has lower computational demands. Improvements in image quality via the proposed method represent a significant opportunity to advance the state-of-the-art in ESWL techniques.
ACKNOWLEDGMENT
This study was supported by Hunan Provincial Natural Science Foundation of China (Grant No. 10JJ5020), Scientific Research Fund of Hunan Provincial Education Department (Grant No. 11C1215), Scientific Research Fund of Hunan Science and Technology Department (Grant No. 2011GK3205).