OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editors: Andrew Dunn and Anthony Durkin
  • Vol. 8, Iss. 6 — Jun. 27, 2013
« Show journal navigation

Wavelet-based noise-model driven denoising algorithm for differential phase contrast mammography

Carolina Arboleda, Zhentian Wang, and Marco Stampanoni  »View Author Affiliations


Optics Express, Vol. 21, Issue 9, pp. 10572-10589 (2013)
http://dx.doi.org/10.1364/OE.21.010572


View Full Text Article

Acrobat PDF (16231 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Traditional mammography can be positively complemented by phase contrast and scattering x-ray imaging, because they can detect subtle differences in the electron density of a material and measure the local small-angle scattering power generated by the microscopic density fluctuations in the specimen, respectively. The grating-based x-ray interferometry technique can produce absorption, differential phase contrast (DPC) and scattering signals of the sample, in parallel, and works well with conventional X-ray sources; thus, it constitutes a promising method for more reliable breast cancer screening and diagnosis. Recently, our team proved that this novel technology can provide images superior to conventional mammography. This new technology was used to image whole native breast samples directly after mastectomy. The images acquired show high potential, but the noise level associated to the DPC and scattering signals is significant, so it is necessary to remove it in order to improve image quality and visualization. The noise models of the three signals have been investigated and the noise variance can be computed. In this work, a wavelet-based denoising algorithm using these noise models is proposed. It was evaluated with both simulated and experimental mammography data. The outcomes demonstrated that our method offers a good denoising quality, while simultaneously preserving the edges and important structural features. Therefore, it can help improve diagnosis and implement further post-processing techniques such as fusion of the three signals acquired.

© 2013 OSA

1. Introduction

Breast cancer is the most ubiquitous cancer in women. Its detection in early stages increments the survival rate and helps improve the life quality of the patient. Nowadays, absorption-based mammography is the most frequently used detection technique. However, it is not able to detect some important malignancies, due to small differences in x-ray attenuation between them and healthy tissue [1

1. P. Sakellaropoulos, L. Costaridou, and G. Panayiotakis, “A wavelet-based spatially adaptive method for mammographic contrast enhancement,” Phys. Med. Biol. 48(6), 787–803 (2003) [CrossRef] [PubMed] .

].

Phase-contrast and scattering-based x-ray imaging can provide complementary information to the conventional absorption-based method, since they are able to detect subtle electron density differences and measure microscopic density variations in the sample. The acquisition of these types of images can be easily carried out using a grating-based x-ray interferometer, which is able to simultaneously generate absorption, differential phase contrast (DPC) and scattering (DCI) images of the sample and, in addition, performs well with conventional x-ray sources. The first research work with native, non-fixed whole breast samples, including regular and cancerous breast tissue, was presented recently by our group [2

2. M. Stampanoni, Z. Wang, T. Thring, C. David, E. Roessl, M. Trippel, R. Kubik-Huch, G. Singer, M. K. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mmmography,” Invest. Radiol. 46(12), 801–806 (2011) [CrossRef] [PubMed] .

]. In that work, a differential phase contrast mammography (mammoDPC) demonstrator based on a conventional x-ray tube was designed, constructed and used to image whole native breast samples directly after mastectomy. In practice, the recorded DPC and scattering images show a significant amount of noise, compared to the absorption signals, because not as many signal averages as in the latter case are acquired. This noise must be removed, for better visualization, diagnosis and further post-processing, such as the image fusion of the three signals [2

2. M. Stampanoni, Z. Wang, T. Thring, C. David, E. Roessl, M. Trippel, R. Kubik-Huch, G. Singer, M. K. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mmmography,” Invest. Radiol. 46(12), 801–806 (2011) [CrossRef] [PubMed] .

].

In the last two decades, the wavelet transform has gained popularity as a promising denoising technique, due to properties such as sparsity and multi-resolution structure. The discrete wavelet transform generates a scale-space representation of a signal, by successively low-pass and high-pass filtering it [3

3. E. Jerhotova, J. Svihlik, and A. Prochazka, “Biomedical image volumes denoising via the wavelet transform,” in Applied Biomedical Engineering, G.D. Gargiulo and A. McEwan, eds. (Intech, 2011), pp 435–458 (2009).

].

The wavelet-based denoising approach operates as follows: 1) The noisy image is wavelet transformed. 2) Thresholding is applied to the wavelet coefficients in absolute value, assuming that the smallest coefficients correspond to noise. The coefficients above the threshold are left untouched, when using a hard-thresholding method, or shrinked by the value of the threshold, when sticking to a soft-thresholding strategy, which in practice produces more visually agreeable images over the former, as the hard-thresholding method is discontinuous and produces abrupt artifacts. 3) An inverse wavelet transform is applied [4

4. S. Rangarajan, R. Venkataramanan, and R. Shah, “Image denoising using wavelets,” (technical report, 2002).

].

Several methods, based on the orthogonal wavelet transform, have been proposed for threshold selection. Donoho and Johnstone [5

5. D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage,” Biometrika 81(3), 425–455 (2009) [CrossRef] .

] introduced the universal threshold VisuShrink, which is the optimal threshold in the asymptotic sense and minimizes the L2 norm of the difference between the unknown noiseless image and its thresholded version. Although this global threshold might be a good starting point when little is known about the signal, its performance is usually insufficient, because it lacks spatial adaptivity, kills many useful signal coefficients, causing image oversmoothing [6

6. P. Kenterlis and D. Salonikidis, “Evaluation of wavelet domain methods for image denoising,” (technical report).

], and yields visual artifacts [7

7. S. K. Mohideen, S. A. Perumal, and M. M Sathik, “Image denoising using discrete wavelet transform,” J. Comput. Sci. 8(1), 8–11 (2011).

].

Later on, Chang, Yu and Vetterli brought in BayesShrink [8

8. S. G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. Image Process. 9(9), 1532–1546 (2000) [CrossRef] .

]. This threshold is derived in a Bayesian framework, is adaptive to each wavelet decomposition subband and depends on data-driven estimates of the parameters. More recently, Chen, Bui and Krzyzak [9

9. G. Y. Chen, T. D. Bui, and A. Krzyzak, “Image denoising using neighbouring wavelet coefficients,” in Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing (IEEE, 2004), pp. 917–920.

] introduced NeighShrink, which thresholds the wavelet coefficients by considering the magnitude of the sum of the squares of the wavelet coefficients within a pre-defined neighboring window, under the assumption that a large wavelet coefficient will have large coefficients as its neighbors. BayesShrink and NeighShrink certainly outperform VisuShrink, mainly because they are subband-adaptive.

Although not as robust as the wavelet-based denoising methods, there have been some useful spatial-domain based denoising strategies. The first one is the Wiener filter, a method that relies on the assumption that the noisy image f can be expressed as the sum of the noiseless image x and the noise n, which are both stationary Gaussian processes. It aims to find a linear estimate of x that minimizes the mean square error. When x is also a white Gaussian process, the filter has a very simple scalar form that only depends on the noise variance and local estimators of the signal variances and means, which can be easily computed [10

10. Y. Jin, E. Angelini, and A. Laine, “Wavelets in medical image processing: denoising, segmentation and registration,” in Handbook of Biomedical Image Analysis, Volume 1: Segmentation models, Part A, J. S. Suri, D.L. Wilson, and S. Laxminarayan, eds. (Kluwer Academic/Plenum Publishers, 2005), pp. 305–358 [CrossRef] .

]. Chambolle et al [11

11. J.F. Aujol and G. Gilboa, “Constrained and SNR-based solutions for TV-Hilbert space image denoising,” J. Math. Imaging Vision 26,217–237 (2006) [CrossRef] .

] proposed a denoising algorithm based on the total variation minimization. They defined the total variation as a weighted sum of a term favoring smoothness, based on the local derivatives, and a term accounting for the fidelity to the input image. Through a fixed point algorithm, they were able to find the optimal weighting factor and compute the noiseless image simultaneously. Even though their approach results certainly interesting from the theoretical point of view, it has shown to yield oversmoothed outcomes.

The denoising methods described above can be classified as “blind” methods, because they assume that the noise distributes in a Gaussian manner along the whole image and use certain estimators to calculate the variance of this noise. However, in mammoDPC, the noise models of the absorption, DPC and scattering signals have been studied and it has been found that the noise variance, instead of being the same for the whole image, varies in a pixelwise manner and can be calculated [12

12. V. Revol, C. Kottler, R. Kaufmann, U. Straumann, and C. Urban, “Noise analysis of grating-based x-ray differential phase contrast imaging,” Rev. Sci. Instrum. 81,073709 (2010) [CrossRef] [PubMed] .

]. Therefore, a denoising method considering these pixelwise noise maps is expected to be more reliable than one assuming uniform Gaussian noise.

One potential choice to solve this denoising problem is the spatial-domain method proposed by Kisilev, Shaked and Hwan Li [13

13. P. Kisilev, D. Shaked, and S. H. Lim, “Noise and signal activity maps for better imaging algorithms,” in Proceedings of IEEE International Conference on Image Processing (IEEE, 2007), pp.117–120.

], which claims that the visual perception of noise depends on both the noise and the underlying image content or signal activity. Based on this argument, known as the masking effect, they introduced a pixelwise noise threshold which is proportional to the ratio of the noise standard deviation to the signal activity. The signal activity was calculated by using local directional derivatives. Another option to solve this particular problem is the Wiener filter, violating the assumption of uniform Gaussian noise and using a different noise variance value for each pixel; a more robust strategy is a wavelet-domain based method incorporating the available noise model, because this transform has shown remarkable advantages. For the design of this method, two issues have to be addressed: 1. The noise variance must be translated from the spatial domain to the wavelet domain. 2. A pixelwise threshold in the latter domain has to be defined.

The results obtained show that the method proposed manages to denoise the images, while preserving the edges and improving the image quality simultaneously. Therefore, this algorithm can certainly be useful for additional post-processing, visualization of structures and diagnosis.

2. Materials and methods

2.1. Signal acquisition

A differential phase-contrast mammography prototype [2

2. M. Stampanoni, Z. Wang, T. Thring, C. David, E. Roessl, M. Trippel, R. Kubik-Huch, G. Singer, M. K. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mmmography,” Invest. Radiol. 46(12), 801–806 (2011) [CrossRef] [PubMed] .

] was used for signal acquisition in the present work. This prototype consists of a Talbot-Lau interferometer, that is composed of a Seifert ID 3000 x-ray generator, an unfiltered tungsten line focus tube (operated at 40 kVp with mean energy of 28 keV and a current of 25 mA), a 3-grating interferometer (with periods p0 = 14μm, p1 = 3.5μm and p2 = 2μm) and a Hamamatsu C9732DK flat panel CMOS detector with a 12cm × 12cm field of view and a 50μm × 50μm pixel size. The interferometer was tuned at the fifth Talbot distance, yielding an angular sensitivity of 0.15 μrad. The distance between the source and the G1 grating was 140 cm and the geometry could be considered parallel; thus, flat gratings were employed. The resulting field of view, limited by the size of the G1 and G2 gratings, was 5cm × 5cm, so to get an image of the whole breast various acquisitions were stitched together. The exposure time was set to 9 seconds. To separate absorption, phase, and scattering signals, a phase-stepping approach, where G2 is translated perpendicular to the grating lamina by a fractional distance of the grating period, was followed. For this study, minimum 8 phase steps were applied, resulting in a total exposure time of at least 72 seconds for the 5cm × 5cm field of view. This fact explains the very low noise level in the absorption images.

2.2. Noise calculation

The noise model of the absorption, DPC and scattering signals was investigated by Revol et al [12

12. V. Revol, C. Kottler, R. Kaufmann, U. Straumann, and C. Urban, “Noise analysis of grating-based x-ray differential phase contrast imaging,” Rev. Sci. Instrum. 81,073709 (2010) [CrossRef] [PubMed] .

]. The detector quantum noise and the phase-stepping jitter noise are the two main noise sources in differential phase contrast imaging. In our case, the phase stepping jitter noise is insignificant because the mechanical uncertainty of the piezo motor is in the order of nm[16

16. PIHera Piezo Lineartisch manual PI. P-620.1 P-629.1.

], which yields a noise variance that is around 1% of that caused by quantum noise. Therefore, we take the quantum noise as the main noise source.

The detector quantum noise variance was found to be linearly proportional to the correspondent mean intensity:
σI,det2=f1I
(1)

The slope f1 has to do with the signal and noise transfer of the incoming x-rays to the output in arbitrary digital units. Due to the beam hardening effect, f1 is different for the reference (flat field) scan ( f1r) and the sample measurement ( f1s). The variances σI translate into errors in the transmission(T), DPC and dark-field (V) images σT,det, σDPC,det and σV,det, respectively, that can be computed employing the error propagation formula and Eqs. (2)(4):
σT,det2T2=f1rNpsa0r(1+f1sTf1r)
(2)
σDPC,det2=f1r2π2vr2Npsa0r(1+f1sTf1rV2)
(3)
σV,det2V2=f1rvr2Npsa0r[vr2(1+f1sTf1r)+2(1+f1sTf1rV2)]
(4)
where:

T

Mean intensity of the transmission image.

V

Mean intensity of the dark field image.

a0r

Mean intensity of the reference measurement.

Nps

Number of phase steps acquired over one period p2.

f1r

f1 for the reference measurement.

f1s

f1 for the sample measurement.

vr

Visibility of the phase stepping curve for the reference measurement.

The dark field signal is computed as the ratio of the visibility of the phase stepping curve in the reference measurement to its visibility in the sample measurement.

To calculate f1r and f1s, 50 images were acquired for a fixed tube voltage and 13 different anode currents, ranging from 1 mA to 25 mA in steps of 2 mA, in the presence and absence of a 1mm Al filter. The mean and variance of the recorded intensity histogram for both the filtered and unfiltered cases and each anode current, were approximated in a pixel-wise manner. Afterwards, the average mean intensity and average variance were computed over a 100 × 100 pixel region. Finally, a linear regression was performed. The values found for f1r and f1s in our system were 0.33 and 0.27, respectively.

2.3. Wavelet based Noise model driven denoising algorithm (WND)

This method translates the noise variance from the spatial to the wavelet domain and uses the latter to define a pixelwise soft threshold.

2.3.1. Noise variance translation from the spatial to the wavelet domain

Let us assume that σs(t) corresponds to the known noise standard deviation of pixel t in the spatial domain calculated using Eqs. (3) and (4). Therefore, the noise signal n(t), having that the noise can be assumed pixelwise independent in mammoDPC, can be expressed as:
n(t)=σs(t)ε(t),
(5)
where ε(t) is white Gaussian noise. Take A(ω) as the DTFT of σs(t), so that:
σs(t)=12πππA(ω)ejωtdω,
(6)

Thus, combining Eqs. (5) and (6) and from [19

19. B.M. Priestley, “Evolutionary spectra and non-stationary processes,” J. R. Stat. Soc. Series B 27(2), 204–237 (1965).

],
n(t)=12πππA(ω)ejωtdF(ω),
(7)
with F(ω) the DTFT of ε(t). The variance σs2(t) can be computed as [19

19. B.M. Priestley, “Evolutionary spectra and non-stationary processes,” J. R. Stat. Soc. Series B 27(2), 204–237 (1965).

]:
σs2(t)=E[n(t)n(t)]
(8)
=1(2π)2ππππE[dF(ω1)dF(ω2)]A(ω1)A(ω2)e(jω1t)e(jω2t).
(9)

Let A*(ω) denote the complex conjugate of A(ω). Due to the fact that n(t) is real-valued, A(ω) = A*(−ω)),
σs2(t)=1(2π)2ππππE[dF(ω1)dF(ω2)]A*(ω1)A(ω2)e(jω1t)e(jω2t).
(10)
and since ε(t) is a white noise process, F(ω) has orthogonal increments. Therefore:
E[dF(ω1)dF(ω2)]=2πδ(ω1ω2)dω1dω2,
(11)
Taking Eq. (11) into Eq. (10), we have:
σs2(t)=12πππA*(ω)A(ω)e(jωt)dω.
(12)

Now, let H(s,o)(ω) be the frequency response of the wavelet filters at scale s and orientation o, and H*(s,o)(ω) be its complex conjugate. Then, using a derivation similar to that for Eq. (12), the wavelet domain noise variance is approximately given by [14

14. B. Goossens, A. Pizurica, and W. Philips, “Em-based estimation of spatially variant correlated image noise,” in Proceedings of IEEE International Conference on Image Processing (IEEE, 2008), pp. 1744–1747.

]:
[σwavelet(s,o)(t)]2=12πππH(s,o)(ω)H*(s,o)(ω)A*(ω)A(ω)e(jωt)dω.
(13)

In this way, the connection between the spatial noise variance, σs2(t), and the noise variance in the wavelet domain, [σwavelet(s,o)(t)]2 is set up by Eqs. (6), (12) and (13). For simplicity, from now on, the index (s, o) will be omitted.

2.3.2. Thresholding

Each wavelet coefficient is assumed to follow a Generalized Gaussian Distribution (GGD) corrupted with noise; this GGD has a standard deviation σx(t). Therefore, the following soft threshold is a good approximation to the optimal threshold [15

15. S. G. Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Trans. Image Process. 9(9), 1522–1531 (2000) [CrossRef] .

]:
T(σx)=σwavelet2(t)σx(t)
(14)

For the calculation of σx(t), context modeling, a widespread procedure employed to adapt the coder to varying image characteristics in the field of image compression, is used [15

15. S. G. Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Trans. Image Process. 9(9), 1522–1531 (2000) [CrossRef] .

]. The context Z[t] of each wavelet coefficient Y(t) is defined as a weighted average of the absolute value of its p neighbors. If these p neighbors are placed in a vector ut, then the context becomes:
Z[t]=wTut
(15)

The weighting factor w is found through least squares:
wLS=(UTU)1UT|Y|,
(16)
where U is a N × p matrix with each row being utT for all t and Y is the N × 1 vector containing all the wavelet coefficients Y (t) belonging to a specific wavelet sub-band (N is the number of coefficients per sub-band) In this case, nine neighbors were considered, eight of them taken from the same subband and the other one from the parent subband [15

15. S. G. Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Trans. Image Process. 9(9), 1522–1531 (2000) [CrossRef] .

].

Afterwards, the variance of the random variable is estimated from other coefficients whose context variable are close in value to Z[t], as follows:
σ^x2[t]=max(12L+1[k]BtY[k]2σwavelet2(t),0),
(17)
where Bt comprises the coefficients correspondent to the L closest points above Z[t], the coefficients correspondent to the L closest points below and Y(t). The term σwavelet2(t) shall be subtracted, because {Y(t)} are noisy coefficients and the variance is independent of the noise.

The proposed method can be used either with the orthogonal or the undecimated wavelet transform. However, the former has been found to yield images with evident ring artifacts. Therefore, in this work, the latter was used. Since the wavelet coefficients produced by this kind of transform are correlated, it is necessary to decorrelate them beforehand. For the sth level of decomposition, the coefficients can be separated into 22s sets of uncorrelated coefficients. These sets are given by Y [2st + k1]t, k1 = 0, 1,...,2s− 1 [15

15. S. G. Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Trans. Image Process. 9(9), 1522–1531 (2000) [CrossRef] .

].

2.4. Data sets

Fig. 1 a. Noiseless DPC image, b. Noiseless DCI image
Fig. 2 a. Noisy DPC image, b. Denoised using VisuShrink, c. Denoised by BayesShrink, d. Denoised by NeighShrink, e. Denoised using the median filter, f. Denoised by Signal-activity-maps, g. Denoised using the Wiener filter and h. Denoised by WND.
Fig. 3 a. Noisy DCI image, b. Denoised using VisuShrink, c. Denoised by BayesShrink, d. Denoised by NeighShrink, e. Denoised using the median filter, f. Denoised by Signal-activity-maps, g. Denoised using the Wiener filter and h. Denoised by WND.

2.5. Performance metrics

To evaluate the performance of the implemented denoising methods in the simulations, the Mean Squared Error(MSE) and the Complex Wavelet Structural Similarity (CWSSIM) index were used.

The MSE between two images x and y with the same number of pixels is defined as:
MSE(x,y)=1Mi=1M(xiyi)2,
(20)
where M is the number of pixels of the images.

The CWSSIM index is an extension of the Structural Similarity index (SSIM) to the complex wavelet domain [17

17. M. P. Sampat, Z. Wang, S. Gupta, A. C. Bovik, and M. K. Markey, “Complex wavelet structural similarity: a new image similarity index,” IEEE Trans. Image Process. 18(11), 2385–2401 (2009) [CrossRef] [PubMed] .

]. The principle behind the SSIM is that the Human Visual System (HVS) is mainly adapted to extract the structure of the objects from the images. Consequently, a measure of the structural similarity constitutes a good approximation of perceptual image quality. SSIM tries to ignore non structural distorsions, such as local intensity patterns and is defined in the spatial domain as:
SSIM(x,y)=(2μxμy+C1)(2σxy+C2)(μx2+μy2+C1)(σx2+σy2+C2),
(21)
where C1 and C2 are two small positive constants; μx and μy are the mean of images x and y, respectively; σx and σy are the standard deviation of images x and y, correspondingly; and σxy is simply:
σxy=1Mi=1M(xiμx)(yiμy)
(22)

The maximum value SSIM can achieve is 1, and it is only achieved when x and y are identical.

The original motivation to extend the SSIM to the complex wavelet domain was its high sensitivity to rotation, translation and scaling of the images, in the spatial domain [17

17. M. P. Sampat, Z. Wang, S. Gupta, A. C. Bovik, and M. K. Markey, “Complex wavelet structural similarity: a new image similarity index,” IEEE Trans. Image Process. 18(11), 2385–2401 (2009) [CrossRef] [PubMed] .

]. Additionally, the fact that the HVS is organized in a multirresolution manner and that phase contains more structural information than magnitude in natural images, inspired this wavelet-based index. The CWSSIM is defined as:
CWSSIM(cx,cy)=2|i=1Ncx,icy,i*|+Ki=1N|cx,i|2+i=1N|cy,i|2+K,
(23)
with cx,i and cy,i being the wavelet coefficients of x and y, respectively, at the decomposition sub-band i (N is the number of coefficients per sub-band), and c* denoting the complex conjugate. K is a small positive constant, whose task is to improve the robustness of the method in regions where the Signal-to-Noise Ratio (SNR) is poor [17

17. M. P. Sampat, Z. Wang, S. Gupta, A. C. Bovik, and M. K. Markey, “Complex wavelet structural similarity: a new image similarity index,” IEEE Trans. Image Process. 18(11), 2385–2401 (2009) [CrossRef] [PubMed] .

].

CW SSIM(cx,cy) can be re-written as [17

17. M. P. Sampat, Z. Wang, S. Gupta, A. C. Bovik, and M. K. Markey, “Complex wavelet structural similarity: a new image similarity index,” IEEE Trans. Image Process. 18(11), 2385–2401 (2009) [CrossRef] [PubMed] .

]:
CWSSIM(cx,cy)=2i=1N|cx,i||cy,i|+Ki=1N|cx,i|2+i=1N|cy,i|2+K2|i=1Ncx,icy,i*|+K2i=1N|cx,icy,i*|+K.
(24)

The first component of this product is determined by the magnitude of the wavelet coefficients and, therefore, is equivalent to SSIM. The second component depends on the consistency of the phase changes between cx and cy and achieves its maximum value 1 when their difference in phase is a constant for all i. It is assumed that the structural information of local image features is principally located in the relative phase changes and a constant phase shift does not cause any structural distorsion [17

17. M. P. Sampat, Z. Wang, S. Gupta, A. C. Bovik, and M. K. Markey, “Complex wavelet structural similarity: a new image similarity index,” IEEE Trans. Image Process. 18(11), 2385–2401 (2009) [CrossRef] [PubMed] .

]. Thus, the more similar x and y are, the higher the value of CWSSIM. The maximum value achievable is 1.

When having real data, it is not possible to calculate any of the metrics presented above, because the ”ground truth” is unknown. Therefore, metrics such as the SNR can be employed to assess the denoising algorithms. For the experiment data, the SNR was calculated in the background of the image, as a reference. The SNR was computed as:
SNR=μσ,
(25)

μ represents the mean and σ the standard deviation. Due to the characteristics of the breast anatomy, it is hard to find constant intensity regions inside it and, therefore, not possible to calculate the object SNR. In this case, the visual appearance of the images and some intensity profiles were employed as the main indicator.

3. Results and discussion

3.1. Simulations

The denoised simulated images corresponding to the highest noise level (worst case scenario) are shown in Figs. 2 and 3. The metrics for the whole simulated data set are reported in Table 1, Table 2, Table 3 and Table 4.

Table 1. MSE values for each denoising method and all the simulated noise levels for the DPC images. The different noise leveles were generated by changing the value of a0r.

table-icon
View This Table
| View All Tables

Table 2. CWSSIM for each denoising method and all the simulated noise levels for the DPC images

table-icon
View This Table
| View All Tables

Table 3. MSE values for each denoising method and all the simulated noise levels for the DCI images

table-icon
View This Table
| View All Tables

Table 4. CWSSIM values for each denoising method and all the simulated noise levels for the DCI images

table-icon
View This Table
| View All Tables

It is worth mentioning that the Signal Activity Maps (SAM) [Figs. 2(f) and 3(f)] and WND [Figs. 2(h) and 3(h)] methods are the best at preserving the edges. However, the former slightly overenhances them, while the latter preserves better their natural appearance. Conversely, the median filter attenuates them strongly [Figs. 2(e) and 3(e)], while VisuShrink completely loses those corresponding to the two oval lobes located inside the phantom and distorts the remaining ones [Figs. 2(b) and 3(b)]. Aditionally, this global threshold yields an oversmoothed image with ring artifacts. This issue constitutes a big advantage of the present denoising algorithm, because for samples such as the breast it is crucial that edges and local features are maintained after denoising.

According to the MSE metric, NeighShrink appears to be the most appropriate technique for both signals [Figs. 2(d) and 3(d)]. Nevertheless, it can be seen that the output yielded by this method is a little oversmoothed and edges are not that well preserved. In addition, it is evident that theVisuShrink outcome [Figs. 2(b) and 3(b)] is much worse than the WND result [Figs. 2(h) and 3(h)], although the MSE value associated to the latter is 1.24 times higher than that associated to the former. Even though, for more than 50 years, the MSE has been the dominant quantitative performance metric in the signal processing field, is simple to calculate and has a clear physical meaning, it has proven to fail to predict human perception of signal quality and fidelity, because it does not take into account any dependency between pixels, such as ordering, textures and patterns. Additionally, it completely ignores the signs of the error and the effects of these signs are significant from a perceptual point of view [18

18. Z. Wang and A. C. Bovik, “Mean-square error : Love it or leave it?,” IEEE Signal Process Mag. 26(1), 98–117 (2009) [CrossRef] .

]. Thus, the optimal denoising method can not be selected solely based on the MSE and more reliable metrics should be used.

A safer metric would be the CWSSIM, because it focuses on the distortions that are meant to cause important perceptually alterations [17

17. M. P. Sampat, Z. Wang, S. Gupta, A. C. Bovik, and M. K. Markey, “Complex wavelet structural similarity: a new image similarity index,” IEEE Trans. Image Process. 18(11), 2385–2401 (2009) [CrossRef] [PubMed] .

]. Correspondingly, for the DPC signals, our proposed method yields the highest CWSSIM [Fig. 2(h)]. However, in the DCI case it indicates that VisuShrink is the best denoising algorithm for the four highest noise levels, although [Fig. 3(b)] shows that this method causes oversmoothing, edge loss and ring artifacts. This result could be explained by the fact that the SNR of the DCI images is significantly lower than that of the DPC images, and therefore the CWSSIM works better for the latter kind of signals (see Table 5). In order to make the CWSSIM safer for DCI signals, an SNR threshold could be established, so we could make sure that above certain value the calculated CWSSIM would be trustable.

Table 5. SNR for the DPC and DCI mammograms calculated in a background ROI

table-icon
View This Table
| View All Tables

For the DCI signal, WND has the largest CWSSIM for the lowest noise level, and the second largest for the rest; this result is more coherent, because as [Fig. 3(h)] proves, this technique causes no structural distortions at all and has a notable edge preservation ability.

3.2. Mammography data

A ROI surrounding a tumor was selected, to illustrate better the denoising effects (Figs. 4 and 5), and an intensity profile of one of the rows was selected to make more clear these effects (Figs. 6 and 7). The mammogram’s background SNR values computed are reported in Table 5.

Fig. 4 a. Noisy DPC image, b. Denoised using VisuShrink, c. Denoised by BayesShrink, d. Denoised by NeighShrink, e. Denoised using the median filter, f. Denoised by Signal-activity-maps, g. Denoised using the Wiener filter and h. Denoised by WND. The red line marks the position of the row profile plotted in Fig. 5.
Fig. 5 a. Noisy DCI image, b. Denoised using VisuShrink, c. Denoised by BayesShrink, d. Denoised by NeighShrink, e. Denoised using the median filter, f. Denoised by Signal-activity-maps, g. Denoised using the Wiener filter and h. Denoised by WND. The red line marks the position of the row profile plotted in Fig. 6.
Fig. 6 Mammography DPC image row intensity profile showing the denoising effects of the methods implemented. The row selected traverses the tumor and the big spike corresponds to a salient edge. Blue and red represent noisy and denoised respectively. a. Denoised using VisuShrink, b. Denoised by BayesShrink, c. Denoised by NeighShrink, d. Denoised using the median filter, e. Denoised by Signal-activity-maps, f. Denoised using the Wiener filter and g. Denoised by WND.
Fig. 7 Mammography DCI image row intensity profile showing the denoising effects of the methods implemented. The row selected traverses the tumor and the big spike corresponds to a salient edge. Blue and red represent noisy and denoised respectively. a. Denoised using VisuShrink, b. Denoised by BayesShrink, c. Denoised by NeighShrink, d. Denoised using the median filter, e. Denoised by Signal-activity-maps, f. Denoised using the Wiener filter and g. Denoised by WND.

As can be grasped from the tables, VisuShrink yields the highest SNR value for both signals. This is understandable, since these values were computed in constant intensity regions with no edges, so the oversmoothing caused by this technique boosts the SNR. Additionally, it is worth to point out that the WND method yielded the highest SNR among the three noise model-based methods in the DPC mammograms. Nevertheless, a disadvantage of the SNR is that it does not assess the preservation of edges and small features. Since these characteristics are important in mammography, this metric is not enough, although it can be a good starting evaluation point, e.g. to confirm that the noise level was indeed reduced. Consequently, to be able to perform a more reliable evaluation, we plotted the intensity profile of one of the rows traversing the tumor (Figs. 6 and 7). These plots confirm that VisuShrink kills many important signal components and attenuates the edges. They also evidence that the median filter causes the largest reduction to the edge strength. Conversely, the WND method manages to denoise the signal while simultaneously preserving the borders. The latter effect is easily visualized, but cannot be quantified precisely.

Regarding all the results obtained, it can be stated that the blind denoising techniques are able to provide an acceptable solution as long as they take into account statistical characteristics of the image and do not rely exclusively on the noise standard deviation. For instance, VisuShrink is the least optimal threshold of this kind, because it ignores these characteristics and strongly depends on the estimated variance value. On the contrary, BayesShrink and NeighShrink offer better solutions in general, because they consider this information.

Another important conclusion is that it is difficult to provide a reliable assessment of a denoising method, relying solely on the existent automatic performance metrics. As it was proved, there is still a long way to go over to develop a robust and exact metric.

A relevant and unique characteristic, and also advantage, of this denoising task is that the variance of the noise associated to each pixel can be calculated. Even so, it is fundamental to keep in mind that in image acquisition systems, there is always random noise that is not possible to calculate. In fact, due to this additional noise, the denoising process yielded much better outcomes in the simulated scenario. Therefore, albeit the noise a-priori knowledge driven denoising methods are more trustable, they are still not perfect. A good way to reduce the uncertainty of these methods is to improve the acquisition hardware, e.g. gratings and detectors.

4. Conclusion

Several other blind as well as model-based existing denoising methods were also implemented. All these denoising techniques were tested with simulated and real mammography DPC and DCI data. Two performance metrics, the Mean-Squared error (MSE) and the Complex Wavelet Structural Similarity Index (CWSSIM), were used for evaluation in the simulated scenario. Since these metrics demand knowing the ground truth, they can not be employed when dealing with real data. Therefore, the SNR and was used to assess the real data denoising results.

The results yielded by WND were satisfactory. Its outcomes demonstrated a reduction of the noise while showing preservation of the edges in both the simulated and real scenarios. The CWSSIM and SNR metrics also gave a good grade to this method. However, as was largely discussed, these metrics are not equivalent to the human visual system and, therefore, it would be ideal to carry out a reader study with experts, who can provide a more reliable assessment of the proposed method. Another important issue is the noise calculation. The noise model for the three signals has been studied and the noise variance associated to each pixel can be estimated. However, this model only considers the detector and phase-jitter noise and, in this work, the latter was not included because, according to the manual of the motor used to move the grating G2, it was insignificant [16

16. PIHera Piezo Lineartisch manual PI. P-620.1 P-629.1.

]. Even so, it would be worth to get to measure this uncertainty and try to include other possible sources of noise, disregarded by this model, in the noise calculation, so it can be more reliable and the denoising can yield better outcomes.

As to the beam hardening effect, it depends on both the sample composition and thickness, and the x-ray energy or tube voltage. However, since the grating-based x-ray interferometry system is designed for a fixed tube voltage, in our case only the sample properties have an effect on beam hardening. To image different objects, in principle, f1s should be measured and calibrated with different materials and sample thicknesses in advance. Nonetheless, our major interest here is mammography, and for imaging breast samples, the beam hardening effect is not pronounced, as presented in [20

20. J. Kaufhold, J. A. Thomas, W. Eberhard, C. E. Galbo, and D. E. Trotter, “A calibration approach to glandular tissue composition estimation in digital mammography,” Med. Phys. 29,1867–1880 (2002) [CrossRef] [PubMed] .

] and as we can grasp from the similarity between f1r and f1s. Breast samples are usually assumed to have a simple composition (fat and glandular tissue, which have close attenuation properties), therefore f1s is estimated with only one material in our case. For situations with severe beam hardening and complex sample compositions with high atomic number materials, we suggest to do the calibration for f1s as described above, in order to get a good estimation of the noise. Since the noise models we are using take into account the difference in f1 between the reference and sample measurements, we believe that our denoising approach would still be useful in these situations.

In conclusion, a working noise model based denoising algorithm for Differential phase contrast mammography was designed, implemented and evaluated, yielding good results and showing potential in this emerging field.

Acknowledgments

The authors thank the Paul Scherrer Institut (PSI) for providing them with the facilities to develop this project and the TOMCAT group, for their valuable inputs and comments. Part of this work has been supported by the ERC Grant ERC-2012-StG 310005-PhaseX.

References and links

1.

P. Sakellaropoulos, L. Costaridou, and G. Panayiotakis, “A wavelet-based spatially adaptive method for mammographic contrast enhancement,” Phys. Med. Biol. 48(6), 787–803 (2003) [CrossRef] [PubMed] .

2.

M. Stampanoni, Z. Wang, T. Thring, C. David, E. Roessl, M. Trippel, R. Kubik-Huch, G. Singer, M. K. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mmmography,” Invest. Radiol. 46(12), 801–806 (2011) [CrossRef] [PubMed] .

3.

E. Jerhotova, J. Svihlik, and A. Prochazka, “Biomedical image volumes denoising via the wavelet transform,” in Applied Biomedical Engineering, G.D. Gargiulo and A. McEwan, eds. (Intech, 2011), pp 435–458 (2009).

4.

S. Rangarajan, R. Venkataramanan, and R. Shah, “Image denoising using wavelets,” (technical report, 2002).

5.

D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage,” Biometrika 81(3), 425–455 (2009) [CrossRef] .

6.

P. Kenterlis and D. Salonikidis, “Evaluation of wavelet domain methods for image denoising,” (technical report).

7.

S. K. Mohideen, S. A. Perumal, and M. M Sathik, “Image denoising using discrete wavelet transform,” J. Comput. Sci. 8(1), 8–11 (2011).

8.

S. G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. Image Process. 9(9), 1532–1546 (2000) [CrossRef] .

9.

G. Y. Chen, T. D. Bui, and A. Krzyzak, “Image denoising using neighbouring wavelet coefficients,” in Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing (IEEE, 2004), pp. 917–920.

10.

Y. Jin, E. Angelini, and A. Laine, “Wavelets in medical image processing: denoising, segmentation and registration,” in Handbook of Biomedical Image Analysis, Volume 1: Segmentation models, Part A, J. S. Suri, D.L. Wilson, and S. Laxminarayan, eds. (Kluwer Academic/Plenum Publishers, 2005), pp. 305–358 [CrossRef] .

11.

J.F. Aujol and G. Gilboa, “Constrained and SNR-based solutions for TV-Hilbert space image denoising,” J. Math. Imaging Vision 26,217–237 (2006) [CrossRef] .

12.

V. Revol, C. Kottler, R. Kaufmann, U. Straumann, and C. Urban, “Noise analysis of grating-based x-ray differential phase contrast imaging,” Rev. Sci. Instrum. 81,073709 (2010) [CrossRef] [PubMed] .

13.

P. Kisilev, D. Shaked, and S. H. Lim, “Noise and signal activity maps for better imaging algorithms,” in Proceedings of IEEE International Conference on Image Processing (IEEE, 2007), pp.117–120.

14.

B. Goossens, A. Pizurica, and W. Philips, “Em-based estimation of spatially variant correlated image noise,” in Proceedings of IEEE International Conference on Image Processing (IEEE, 2008), pp. 1744–1747.

15.

S. G. Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Trans. Image Process. 9(9), 1522–1531 (2000) [CrossRef] .

16.

PIHera Piezo Lineartisch manual PI. P-620.1 P-629.1.

17.

M. P. Sampat, Z. Wang, S. Gupta, A. C. Bovik, and M. K. Markey, “Complex wavelet structural similarity: a new image similarity index,” IEEE Trans. Image Process. 18(11), 2385–2401 (2009) [CrossRef] [PubMed] .

18.

Z. Wang and A. C. Bovik, “Mean-square error : Love it or leave it?,” IEEE Signal Process Mag. 26(1), 98–117 (2009) [CrossRef] .

19.

B.M. Priestley, “Evolutionary spectra and non-stationary processes,” J. R. Stat. Soc. Series B 27(2), 204–237 (1965).

20.

J. Kaufhold, J. A. Thomas, W. Eberhard, C. E. Galbo, and D. E. Trotter, “A calibration approach to glandular tissue composition estimation in digital mammography,” Med. Phys. 29,1867–1880 (2002) [CrossRef] [PubMed] .

OCIS Codes
(100.2000) Image processing : Digital image processing
(100.7410) Image processing : Wavelets
(110.3000) Imaging systems : Image quality assessment
(170.3830) Medical optics and biotechnology : Mammography

ToC Category:
Image Processing

History
Original Manuscript: January 29, 2013
Revised Manuscript: March 31, 2013
Manuscript Accepted: April 8, 2013
Published: April 23, 2013

Virtual Issues
Vol. 8, Iss. 6 Virtual Journal for Biomedical Optics

Citation
Carolina Arboleda, Zhentian Wang, and Marco Stampanoni, "Wavelet-based noise-model driven denoising algorithm for differential phase contrast mammography," Opt. Express 21, 10572-10589 (2013)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-21-9-10572


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. P. Sakellaropoulos, L. Costaridou, and G. Panayiotakis, “A wavelet-based spatially adaptive method for mammographic contrast enhancement,” Phys. Med. Biol.48(6), 787–803 (2003). [CrossRef] [PubMed]
  2. M. Stampanoni, Z. Wang, T. Thring, C. David, E. Roessl, M. Trippel, R. Kubik-Huch, G. Singer, M. K. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mmmography,” Invest. Radiol.46(12), 801–806 (2011). [CrossRef] [PubMed]
  3. E. Jerhotova, J. Svihlik, and A. Prochazka, “Biomedical image volumes denoising via the wavelet transform,” in Applied Biomedical Engineering, G.D. Gargiulo and A. McEwan, eds. (Intech, 2011), pp 435–458 (2009).
  4. S. Rangarajan, R. Venkataramanan, and R. Shah, “Image denoising using wavelets,” (technical report, 2002).
  5. D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage,” Biometrika81(3), 425–455 (2009). [CrossRef]
  6. P. Kenterlis and D. Salonikidis, “Evaluation of wavelet domain methods for image denoising,” (technical report).
  7. S. K. Mohideen, S. A. Perumal, and M. M Sathik, “Image denoising using discrete wavelet transform,” J. Comput. Sci.8(1), 8–11 (2011).
  8. S. G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. Image Process.9(9), 1532–1546 (2000). [CrossRef]
  9. G. Y. Chen, T. D. Bui, and A. Krzyzak, “Image denoising using neighbouring wavelet coefficients,” in Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing (IEEE, 2004), pp. 917–920.
  10. Y. Jin, E. Angelini, and A. Laine, “Wavelets in medical image processing: denoising, segmentation and registration,” in Handbook of Biomedical Image Analysis, Volume 1: Segmentation models, Part A, J. S. Suri, D.L. Wilson, and S. Laxminarayan, eds. (Kluwer Academic/Plenum Publishers, 2005), pp. 305–358. [CrossRef]
  11. J.F. Aujol and G. Gilboa, “Constrained and SNR-based solutions for TV-Hilbert space image denoising,” J. Math. Imaging Vision26,217–237 (2006). [CrossRef]
  12. V. Revol, C. Kottler, R. Kaufmann, U. Straumann, and C. Urban, “Noise analysis of grating-based x-ray differential phase contrast imaging,” Rev. Sci. Instrum.81,073709 (2010). [CrossRef] [PubMed]
  13. P. Kisilev, D. Shaked, and S. H. Lim, “Noise and signal activity maps for better imaging algorithms,” in Proceedings of IEEE International Conference on Image Processing (IEEE, 2007), pp.117–120.
  14. B. Goossens, A. Pizurica, and W. Philips, “Em-based estimation of spatially variant correlated image noise,” in Proceedings of IEEE International Conference on Image Processing (IEEE, 2008), pp. 1744–1747.
  15. S. G. Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Trans. Image Process.9(9), 1522–1531 (2000). [CrossRef]
  16. PIHera Piezo Lineartisch manual PI. P-620.1 P-629.1.
  17. M. P. Sampat, Z. Wang, S. Gupta, A. C. Bovik, and M. K. Markey, “Complex wavelet structural similarity: a new image similarity index,” IEEE Trans. Image Process.18(11), 2385–2401 (2009). [CrossRef] [PubMed]
  18. Z. Wang and A. C. Bovik, “Mean-square error : Love it or leave it?,” IEEE Signal Process Mag.26(1), 98–117 (2009). [CrossRef]
  19. B.M. Priestley, “Evolutionary spectra and non-stationary processes,” J. R. Stat. Soc. Series B27(2), 204–237 (1965).
  20. J. Kaufhold, J. A. Thomas, W. Eberhard, C. E. Galbo, and D. E. Trotter, “A calibration approach to glandular tissue composition estimation in digital mammography,” Med. Phys.29,1867–1880 (2002). [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited