Denoising and enhancing digital mammographic images for visual screening

Denoising and enhancing digital mammographic images for visual screening

Computerized Medical Imaging and Graphics 30 (2006) 243–254 Denoising and enhancing digital mammographic images for visual screening Jacob Scharcansk...

696KB Sizes 0 Downloads 23 Views

Computerized Medical Imaging and Graphics 30 (2006) 243–254

Denoising and enhancing digital mammographic images for visual screening Jacob Scharcanski a,∗ , Cl´audio Rosito Jung b,1 b

a UFRGS, Universidade Federal do Rio Grande do Sul Caixa Postal 15064, 91501-970 Porto Alegre, RS, Brasil UNISINOS, Universidade do Vale do Rio dos Sinos, PIPCA, Graduate School of Applied Computing, Av. UNISINOS, 950. S˜ao Leopoldo, RS 93022-000, Brasil

Abstract Dense regions in digital mammographic images are usually noisy and have low contrast, and their visual screening is difficult. This paper describes a new method for mammographic image noise suppression and enhancement, which can be effective particularly for screening image dense regions. Initially, the image is preprocessed to improve its local contrast and the discrimination of subtle details. Next, image noise suppression and edge enhancement are performed based on the wavelet transform. At each resolution, coefficients associated with noise are modelled by Gaussian random variables; coefficients associated with edges are modelled by Generalized Laplacian random variables, and a shrinkage function is assembled based on posterior probabilities. The shrinkage functions at consecutive scales are combined, and then applied to the wavelets coefficients. Given a resolution of analysis, the image denoising process is adaptive (i.e. does not require further parameter adjustments), and the selection of a gain factor provides the desired detail enhancement. The enhancement function was designed to avoid introducing artifacts in the enhancement process, which is essential in mammographic image analysis. Our preliminary results indicate that our method allows to enhance local contrast, and detect microcalcifications and other suspicious structures in situations where their detection would be difficult otherwise. Compared to other approaches, our method requires less parameter adjustments by the user. © 2006 Elsevier Ltd. All rights reserved. Keywords: Mammography; Contrast equalization; Image denoising; Image enhancement; Multiresolution analysis; Wavelets

1. Introduction Breast cancer currently accounts for more than 30% of cancer incidence, and a significant percentage of cancer mortality in both developing and developed countries. It has been shown that early detection and treatment of breast cancer are the most effective methods of reducing mortality. Despite of advances in resolution and film contrast, screenfilm mammography remains a diagnostic imaging modality where image interpretation is very difficult [1]. Breast radiographs are generally examined for the presence of malignant masses and indirect signs of malignancy such as microcalcifications and skin thickening. A significant effort has been di-



Corresponding author. Tel.: +55 51 3316 7128; fax: +55 51 3319 1576. Tel.: +55 51 3591 1122x1626; fax: +55 51 3590 8162. E-mail addresses: [email protected] (J. Scharcanski), [email protected] (C.R. Jung). 1

0895-6111/$ – see front matter © 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2006.05.002

rected to improve imaging performance, but it is unlikely that improvements will be achieved only by advances in screenfilm radiography. In general, the visualization of mammograms displays a small percentage of the information available [1]. This deficiency of the mammographic technology is caused by the fact that, in general, there are small differences in X-ray attenuation between normal glandular and malignant tissues. Detection of small malignancies is specially difficult in younger women who tend to present denser breast tissue. On the other hand, calcifications have high attenuation properties (because these are denser tissues, similar to bones), but are small in size, and tend to present low local contrast. Therefore, the visibility of small tumors, and any associated microcalcifications, is a problem in the mammography technology based on analog film. Mammographic images are inherently noisy and usually contain low contrast regions. In fact, it is a challenge to improve the visual quality of mammograms by image processing for helping in the early detection of breast cancer. Therefore, two impor-

244

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

tant current problems in mammographic image processing are: (a) improvement of local detail discrimination in low contrast regions, and (b) noise reduction in such images without blurring fine image details (i.e. image edges). Several approaches to improve local contrast, reduce noise and enhance image details (mostly based on the wavelet transform) have been proposed in the literature. Nunes et al. [2] proposed a method for enhancing dense mammograms, that can be useful for the detection of clustered microcalcifications. Their method was tested on scanned images, and for each image, precise equipment calibration and parameter adjustments are required (and also an initial selection of the region of interest). However, nowadays, there is a severe time constraint in medical services provided to the public, where many mammograms must be screened daily. In this work we are interested in less time consuming enhancement methods, that minimize the need for parameter adjustments by the user. Local contrast and image intensity are inter-dependent in mammographic images. In fact, noise tends to increase with pixel intensity in such images, making the discrimination of local details more difficult, specially in dense regions (where the image intensity is higher). Karssemeijer [3] observed this problem and proposed a noise equalization (i.e. isonoise) procedure to obtain images where local contrast is approximately equal at all image intensities. This technique has been improved, and Gurun et al. [4] proposed a method to obtain more reliable estimates of local contrast, and McLoughlin et al. [5] further improved this technique using better noise estimates (specially in the brightest, i.e. most dense, regions). The noise equalization procedure has been proposed as a preprocessing stage to automatic microcalcification detection [5], which can downgrade some image structures for visual screening. However, automatic microcalcification detection algorithms often produce false positives (and false negatives), making direct visual screening necessary to obtain reliable diagnoses. As a consequence, techniques for improving direct visual screening of mammograms are needed in clinical practice. Simoncelli and Adelson [6] proposed a Bayesian estimator to discriminate between image signal and noise (i.e. image denoising), assuming the noise as Gaussian additive and known a priori. Their approach provides a semi-blind noise removal algorithm based on a steerable wavelet pyramid, that can outperform the Wiener filter. Malfait and Roose [7] developed a filtering technique that takes into account two image measurements. The first one measures the local regularity of the image using the H¨older exponent, and the second one takes into account geometric constraints. These two measures are combined in a Bayesian probabilistic framework, and implemented by as a Markov random field model. However, the stochastic sampling procedure, needed in the calculation of the probabilities, is computationally demanding. In the approach of Chang et al. [8], coefficients associated with noise are modelled as Gaussian functions, while coefficients associated to edges are modelled as Generalized Gaussian functions. These probability functions are used to determine a soft threshold, which is applied to the wavelet coefficients. In

a modification of this method [9], context modelling was introduced, and each coefficient is represented by a Generalized Gaussian random variable based on local spatial information. For each coefficient, a threshold is estimated. Mihc¸ak et al. [10] also proposed a spatially-adaptive statistical model for image denoising. Wavelet coefficients are modelled as Gaussian random variables with high local correlation, and a maximum a posteriori probability rule is applied to estimate the original coefficients from noisy observations. Strela et al. [11] described the joint densities of groups of wavelet coefficients as a Gaussian scale mixture (GSM), and developed a maximum likelihood solution for estimating relevant wavelet coefficients (i.e. related to edges) from noisy observations. Portilla and Simoncelli [12] proposed a complex wavelet decomposition and processing for image denoising. The autocorrelations of theoretical noise-free coefficients and their magnitudes are estimated within each subband, and then projected onto spaces having the desired autocorrelations. The image is reconstructed with the modified coefficients, and this process is repeated iteratively, until convergence is achieved. Figueiredo and Nowak [13] proposed a wavelet-based denoising technique without any free parameters, using empirical Bayes estimation based on a Jeffrey’s noninformative prior. The main problem of the approaches mentioned above (except for [13]) is that a noise estimate is needed, which may be difficult to obtain in practical situations, specially for images with inherent noise (e.g. X-ray images, aerial images, etc.). In fact, the reported probabilistic approaches [14,8–13] were not sufficiently tested for these types of images. There are also known approaches designed specifically for mammographic image enhancement, and some of them are discussed next. Yoshida et al. [15] proposed to use the discrete wavelet transform (DWT), multiply every scale by a weight factor, and then reconstruct an image using the inverse DWT. The weights are determined by supervised learning, given a set of training cases. However, the DWT is not translation invariant, meaning that a shift in the image origin leads to results inherently different to the transform applied to the original image. Laine et al. [16,1] used a dyadic wavelet transform [17], and an adaptive enhancement operator on the wavelet coefficient scales. They obtained good contrast improvement for irregular structures such as microcalcifications and masses. The enhancement operation must be defined at each scale separately, and they do not exploit the correlation of wavelet coefficients across scales. Besides, it can be challenging to adjust parameters in order to obtain the desired results. Recently, Heinlein et al. [18] proposed a microcalcification enhancement method based on the continuous wavelet transform (CWT) and a local microcalcification model. Their method requires an empirical selection of appropriate thresholds for image denoising, as well as the specification of an appropriate size range for the structures to be enhanced. This method can not be applied for mammographic image enhancement in general, because the size and shape of the suspicious structures vary significantly in mammograms. In this work, we propose a new adaptive method for mammographic image denoising and enhancement using the wavelet transform, which combines noise equalization, wavelet shrink-

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

age and scale-space constraints. Our approach is flexible enough to allow the user to select the desired image enhancement and scale of analysis, but it does not require the user to adjust any parameters for image denoising. We adopt the approach proposed by Mallat and Zhong [19] to compute a redundant wavelet transform using only two detail images (horizontal and vertical), instead of the classical approach in which three detail images are utilized (horizontal, vertical and diagonal details) [17]. This redundant wavelet transform is known to be virtually shift invariant. The distribution of detail coefficients is modelled by a composition of Gaussian and Generalized Laplacian probability density functions at each scale, a shrinkage function is assembled. Such shrinkage functions are combined in consecutive levels to preserve edges that are persistent over scales and remove residual noise. Finally, an adaptive piecewise linear enhancement function is applied to the denoised wavelet coefficients. The next section describes our contrast enhancement method, the wavelet framework is presented in Section 3, and Section 4 describes the new image denoising and edge enhancement method. Section 5 presents some experimental results, and our conclusions are presented in the final section. 2. The preprocessing stage: contrast enhancement All radiological images contain random fluctuations due to the statistics of X-ray quantum absorption. This noise makes the detection of small and subtle structures more difficult. Usually, the relationship between image intensity and noise variance is nonlinear, and varies significantly from image to image [5]. In this paper, the method proposed by McLoughlin et al. [5] is extended for use on direct screening of digital mammograms. Within a neighborhood  of an image location (x, y), the local contrast is estimated as: c(x, y) = f (x, y) − median (x, y)

(1)

where c(x, y) is the estimated local contrast, f (x, y) is the image gray level at (x, y), and median (x, y) is the median gray level within the neighborhood  of (x, y). Eq. (1) takes the form of a high-pass spatial filter, and local contrast provides a measure of the high frequency image noise. The noise associated with each image gray level I can be measured by the local contrast standard deviation σc (I) ≡ σ{c(I)} (i.e. by the local contrast variability considering all pixels with gray level I). To estimate σc (I) the grayscale is divided into nonoverlapping bins numbered i = 1, . . . , N, where N is the number of bins, and the ith bin has central gray value Ii . For each gray value Ii , the standard deviation σc (Ii ) is calculated for all local contrast values within bin i.1 The gray value variability within bin i is measured by the gray level standard deviation σ(Ii ). Usually, σc (Ii ) and σ(Ii ) increase with Ii , and details in intense regions are severely affected by noise. To improve detail discrimination for visual screening, we correct the contrast in a way that intense details are enhanced w.r.t. noise, and noise is

245

perceived in all gray levels similarly. In other words, details are enhanced with respect to noise in all gray level bins. We approach this problem by defining a contrast enhancement function that is image adaptive, i.e. based on the measured values σ(Ii ) and σc (Ii ). Initially, a transfer function g ≡ σ(Ii )/σc (Ii ) is obtained. At each gray level bin Ii , the transfer function is proportional to the gray level variability σ(Ii ), enhancing the contrast of details. Also, the transfer function provides an equalization of the noise occurrence in the different bins, making the slope at each bin i inversely proportional to the noise intensity σc (Ii ). The gray level variability with the ith bin is captured by σ(Ii ), and increases faster with Ii than σc (Ii ), which accounts essentially for noise. Consequently, when σ(Ii ) increases, 1/σc (Ii ) decreases, and the ratio σ(Ii )/σc (Ii ) tends to increase with Ii . We define the contrast enhancement transfer function as follows: ⎧ ⎨ σ(Ii ) , if σ (I ) > 0; c i fceq (Ii ) = σc (Ii ) (2) ⎩ 0, otherwise. Interpolation of the estimated fceq (Ii ) values provides an estimate of the function fceq (I) for all image intensities I.2 Once the transfer function fceq (Ii ) is obtained, it is normalized in a way that it adds to 1, similar to a discrete probability density function: fceq (I) f¯ ceq (I) = max(I) , (3) j=min(I) fceq (Ii ) where min(I) and max(I) are the minimum and maximum image gray levels, respectively. Finally, we define a contrast equalization function pceq similar to a discrete cumulative density function, which is given by : pceq (I) =

I 

f¯ ceq (j).

(4)

j=min(I)

Therefore, we obtain a contrast equalization function pceq (I) that increases monotonically, and is based on the local contrast and gray level image statistics. The new image gray levels Iceq are then obtained by: Iceq (x, y) = [I(x, y) − min(I)] × pceq (I(x, y)) + min(I), (5) where I(x, y) is the original gray level at the image location (x, y). Our contrast enhancement procedure provides non-linear gray level re-scaling (because pceq (I) is non-linear), which leads to better results compared to a simple gray level normalization.3 The ratio in Eq. (2) increases in average with Ii , implying that our contrast equalization function is non-linear, as Eq. (4) shows. Histogram equalization also is a known contrast enhancement technique, but it does take into account the noise distribution at different gray levels, and can decrease the contrast in intense 2

1

Within each bin, local contrast values are assumed to be normally distributed w.r.t. the median value Ii , and σc (Ii ) describe their dispersion.

In our experiments, we used linear interpolation. The normalized gray levels are given ([I(x, y) − min(I)]2 )/([max(I) − min(I)]) + min(I). 3

by

Inorm (x, y) =

246

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

Fig. 1. Illustration of the data fitting provided by the mixture of Gaussian and Generalized Laplacian model. The wavelet coefficient distributions correspond to vertical details of the contrast enhanced image 148 of the MIAS database: (a) coefficients at scale 21 ; and (b) coefficients at scale 22 (dotted: coefficients distribution; solid: our model).

Fig. 2. Adequacy of the Gaussian noise model for the contrast enhanced image. Normal plots of wavelet coefficients (vertical details) for a simulated Gaussian noise image, that was contrast enhanced with our approach: (a) original Gaussian data; (b) coefficients at scale 21 ; and (b) coefficients at scale 22 .

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

247

image regions (or enhance noise as well as image structures). A comparison between our technique and other approaches is discussed in detail in Section 5. Next, our framework for denoising and enhancing our contrast improved image is presented.

Since we are dealing with digital images f [n, m], we use the discrete version of the wavelet transform [19], and the discrete wavelet coefficients are denoted in this paper by W2i j f [n, m], for i = 1, 2.

3. The wavelet transform

4. Image denoising and enhancement

To compute the redundant wavelet transform with two detail images, a smoothing function φ(x, y) and two wavelets ψi (x, y) are needed. The dilation of these functions are denoted by 1 x y φs (x, y) = 2 φ , , s s s 1 x y , , i = 1, 2, (6) ψsi (x, y) = 2 ψi s s s

Given a digital image f [n, m], we apply the redundant wavelet transform using J dyadic levels. As a result, we obtain 2J detail images W21j f [n, m], W22j f [n, m], for j = 1, . . . , J, and the smoothed image S2J f [n, m]. Next, we describe our method to discriminate coefficients associated with edges from coefficients associated with noise based on a posteriori probabilities. We propagate this information in the scale-space (using consistency along scales), applying an adaptive piecewise linear enhancement function. Finally, we compute the inverse wavelet transform to obtain the processed image.

and the dyadic wavelet transform f (x, y), at a scale s = 2j , has two detail components, given by W2i j f (x, y) = (f ∗ ψ2i j )(x, y),

i = 1, 2,

(7)

and one low-pass component, given by S2j f (x, y) = (f ∗ φ2j )(x, y).

4.1. Wavelet shrinkage (8)

The coefficients W21j f (x, y) and W22j f (x, y) represent the details

in the x and y directions, respectively. Thus, the image gradient at the resolution 2j can be approximated by  W21j f (x, y) W2j f (x, y) = . (9) W22j f (x, y)

Wavelet shrinkage is a known approach for noise reduction, where wavelet coefficients are subject to a non-linearity that reduces (or suppresses) low-amplitude values, and retains highamplitude values [14,20]. For each level 2j , we want to find non-negative nondecreasing shrinkage functions gji (x), 0 ≤ gji (x) ≤ 1, such that wavelet coefficients W21j f and W22j f are updated according to

Fig. 3. Comparative results for the MIAS database mammogram 211, containing a cluster of microcalcifications which is not clearly visible because of the dense tissue, at the bottom right: (a) original image; (b) image enhanced by histogram equalization (local contrast in dense tissues is unsatisfactory); (c) image enhanced using our contrast enhancement approach (local contrast in dense regions is improved, and the cluster of microcalcifications is more visible in the bottom right); and (d) image denoised and enhanced using our method (the cluster of microcalcifications is better defined in the bottom right).

248

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

Fig. 4. Comparative results for the MIAS database mammogram 002, containing two macrocalcifications not clearly visible near the image center: (a) original image; (b) image enhancement by histogram equalization (local contrast is unsatisfactory); (c) image contrast enhancement using our approach (the macrocalcifications are more visible near the image center); and (d) image denoised and enhanced (the macrocalcifications are better defined near the image center).

the following rule: NW i2j f [n, m] = W2i j f [n, m]gji [n, m], i = 1, 2,

(10)

where gji [n, m] = gji (W2i j f [n, m]) are called shrinkage factors. Considering that the same shrinkage procedure (Eq. (10)) will be applied at all scales and subbands of the wavelet decomposition, the indexes j and i will be used only when necessary. To determine a shrinkage function g(x), we analyze the distribution of coefficients Wf [n, m]. In this work, we model the

mammographic image noise by an additive zero-mean Gaussian noise [18]. The detail coefficients Wf [n, m] of an image constituted only by Gaussian noise may be considered Gaussian distributed [21], with standard deviation σnoise .4 Under these circumstances, the function that models the coefficient distribution

4 It shall be noticed that our preprocessing contrast enhancement algorithm preserves the Gaussianity of noise, as illustrated in Fig. 2.

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

is given by p(x|noise) =

1 2 2 √ e−x /2σnoise σnoise 2π

(11)

Therefore, p(x|noise) represents the distribution of coefficients assuming that only noise is present in the image. On the other hand, the distribution of the wavelet coefficients for noise-free images is sharply peaked near the origin (due to homogeneous regions) and has a long tail (due to image edges). In fact, a more suitable function to model such coefficients is a Generalized Laplacian [6,17]. Thus, we assume that the distribution of noise-free wavelet coefficients Wf [n, m] is approximated by a Generalized Laplacian probability density function, given by [22]: x β

e−| α | p(x|edge) = , C(α, β)

for − ∞ < x < ∞, β > 0, α > 0, (12)

where α C(α, β) = 2 Γ β

1 , β

(13)

and Γ (x) is the gamma function. Therefore, p(x|edge) represents the distribution of coefficients assuming that only edges and homogeneous regions are present in the image. The notation

249

p(x|edge) was chosen because we are interested mostly in edges, as will be detailed below. The hypotheses of images containing noise only and edge plus homogeneous regions only are mutually exclusive. This is central to our approach, which models the generic coefficient distribution of images containing noisy edges and homogeneous regions affected by noise. We approximate the above mentioned generic coefficient distribution by a mixture of a Gaussian and a Generalized Laplacian distributions, as discussed next. When we have an image contaminated by zero-mean additive Gaussian noise, the tail of the distribution will not be significantly affected by noise (because typically edge-related coefficients have magnitudes larger than noise-related coefficients). In this work, we assume that coefficients belonging to the tail of the distribution are likely to be edge-related, and coefficients close to the origin of the histogram are likely to be related to noise (and to nearly homogeneous regions). Hence, the origin of the histogram will be Gaussian shaped, and the tail of the distribution will follow a Generalized Laplacian function, and the overall distribution of the coefficients Wf [n, m] (including coefficients related to edges and noise) is given by: p(x) = wp(x|noise) + (1 − w)p(x|edge),

(14)

where w is an unknown parameter of the overall coefficient distribution (0 ≤ w ≤ 1).

Fig. 5. Comparative results for the MIAS database mammogram 148, containing microcalcifications located in dense tissue, which are not clearly visible: (a) original image; (b) image enhancement by histogram equalization (e.g. local contrast is poor in dense tissues, near the microcalcification); (c) image contrast enhancement using our approach (showing details of the microcalcifications, including those located in dense tissue); and (d) image denoised and enhanced (the details of the microcalcifications are visible, including those located in dense tissue).

250

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

The parameters σnoise , α, β and w are estimated by maximizing the likelihood function:  ln L = ln (p(Wf [n, m])), (15) (m,n)inimage

where p(Wf [n, m]) is the function defined in Eq. (14) evaluated at coefficients Wf [n, m]. The second and fourth moments of a Generalized Laplacian corrupted by additive Gaussian white noise, namely m2l and m4l , are given by [6]: α2 Γ (3/β) , Γ (1/β) 6σ 2 α2 Γ (3/β) α4 Γ (5/β) 4 m4l = 3σnoise + . + noise Γ (1/β) Γ (1/β) 2 m2l = σnoise +

(16)

Typically, the variance of noise-related coefficients is smaller than the variance of edge-related coefficients. Also, noise variance decreases as j increases (because of low-pass filtering). j j−1 Therefore, the restrictions σnoise < σedge , σnoise ≤ σnoise also are imposed in the maximization procedure. Besides, the maximization procedure must satisfy the moment equations presented above (see Eq. (16)).

Once the parameters σnoise , α, β and w are estimated, the conditional probability density functions p(x|noise) and p(x|edge) are given, respectively, by Eqs. (11) and (12). The shrinkage function g(x) is given by the posterior probability function p(edge|x), which is computed using Bayes theorem as follows: g(x) =

(1 − w)p(x|edge) . (1 − w)p(x|edge) + wp(x|noise)

(17)

Fig. 1 illustrates how the mixture of Gaussian and Generalized Laplacian model fits the distributions of wavelet coefficients. The wavelet coefficient distributions correspond to vertical details of the contrast enhanced image 148 of the MIAS database, as detailed in Section 5. Also, Fig. 2 shows normal plots of wavelet coefficients (vertical details) of a simulated Gaussian noise image, that was previously contrast enhanced with our approach. These normal plots provide experimental evidence that our Gaussian model for the noise in the WT is still adequate after contrast enhancement. 4.2. Consistency along scales The analysis of the wavelet coefficients at each resolution does not provide sufficient discrimination between edge-related

Fig. 6. Comparative results for the MIAS database mammogram 145, containing a nodule whose boundaries are not clearly visible, located near the image top: (a) original image; (b) image enhancement by histogram equalization (local contrast does not improve much); (c) image contrast enhancement using our approach (the boundaries of the nodule are more visible, near the image top); and (d) image denoised and enhanced (the nodule boundaries are visible, near the image top).

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

and noise-related coefficients, specially in the level 21 , where we find the highest noise contamination. Better discrimination can be achieved by analyzing the information available in different scales. Typically, consistency across scales is explored by computing the Lipschitz exponent of the wavelet coefficients [7], or the direct correlation of coefficients in consecutive scales [23,24]. In this work, we explore consistency across scales simply by combining shrinkage factors in adjacent scales [25], as described below. For each scale 2j , the value gj [n, m] may be interpreted as a confidence measure that coefficient W2j f [n, m] is in fact associated to an edge. To take into account multiple resolutions, we use the information provided by the shrinkage factor gj [n, m], and also by the factors gj+1 [n, m], gj+2 [n, m], . . . , gκ [n, m], where κ − j + 1 is the number of consecutive resolutions that will be taken into consideration for consistency along scales (the values assumed by κ will be described below). We want to combine this information so that, if the value gj [n, m] is close to 1 for

251

several consecutive levels 2j , it is more likely that W2j f [n, m] is associated to an edge. On the other hand, if gj [n, m] decreases as j increases, it is more likely that W2j f [n, m] is actually associated with noise. We update the shrinkage factors according to the following rule:

gj [n, m]γ + · · · + gκ [n, m]γ 1/γ gjscale [n, m] = , (18) κ−j+1 where γ is a positive constant. When γ = 1, the average of the shrinkage factors is obtained. For γ < 1, smaller coefficients carry more weight, and tend to dominate the summation, producing smaller values. This means that when γ is closer to zero, noise reduction is stronger, and when γ is larger, edge preservation is stronger. It was verified experimentally that using γ = 0.8 for different test images (with distinct noise contamination levels), the expected difference in PSNR is smaller than 0.5 dB with respect to the optimal PSNR (determined by exhaustive testing of different γ values). There-

Fig. 7. Comparative results for the MIAS database mammogram 212, containing a large dense mass, located at the image center: (a) original image; (b) image enhancement by histogram equalization (local contrast has not improved much); (c) image contrast enhancement using our approach (details in the dense tissue are more evident); and (d) image denoised and enhanced (details in the dense tissue are even more evident).

252

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

Fig. 8. Comparative contrast enhancement results for the MIAS database mammogram 148: (a) original image; (b) histogram equalization; (c) gray level normalization; and (d) image contrast enhancement using our approach.

fore, we used γ = 0.8 in our experiments with mammographic images. This updating rule is applied from coarser to finer resolutions. The shrinkage factor gJscale [n, m], corresponding to the coarsest resolution 2J , is equal to gJ [n, m]. However, for other resolutions 2j , j = 1, . . . , J − 1, shrinkage factors gjscale [n, m] depend on scales 2j , 2j+1 , . . . , 2κ , where κ = min{J, j + K}, and K + 1 is the maximum number of consecutive scales that will be analyzed. As noticed in [23,24], the support of an isolated edge will increase by a factor of two across scale and the neighboring edges will interfere with each other at coarse scales. Hence, it is sufficient to consider only two adjacent scales (i.e. we use K = 1). Coefficients W2j f [n, m] are then modified according to Eq. (10), using the updated shrinkage factors gjscale [n, m] instead of gj [n, m].5 4.3. Edge enhancement The denoising technique described above can be extended to local contrast enhancement near edges. Linear enhancement is straightforward, and can be obtained by allowing shrinkage factors gjscale [n, m] to be greater than 1, so that local contrast will be enhanced near the edges. Nevertheless, linear enhancement tends emphasize mostly strong edges, and mammograms 5

The subband index i was omitted here to simplify notation.

enhanced by a linear operator containing an obvious macrocalcification (high intensity) will result in gross rescaling with the available dynamic range for image display [1]. Laine et al. [1] proposed a nonlinear enhancement method with the following design guidelines in mind: (a) areas of low contrast are enhanced more than high contrast areas; (b) sharp edges are not blurred; and (c) the enhancement function satisfies the constraints of monotonicity and antisymmetry. These guidelines provide that low contrast features will be more enhanced than the high contrast image features, and that artifacts will not be introduced by the enhancement method. The proposed enhancement function hij , which is used for updating the denoised wavelet coefficients NW i2j f [n, m], is given by the following rule: ⎧ i i ⎪ ⎨ x − (G − 1)T2j , if x < −T2j if |x| < T2ij , hij (x) = Gx, ⎪ ⎩ x + (G − 1)T2ij , if x > T2ij

(19)

where T2ij is the median of coefficient absolute values |NW i2j | at the corresponding scale j and subband i, and G ≥ 1 is a parameter. Laine et al. [1] used a similar nonlinear enhancement method, but in their approach T2ij must be user specified, which can be quite troublesome in practical situations. Our approach is adaptive, and only requires the user to select the desired enhancement level G.

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

Finally, the enhanced wavelet coefficients are computed by: NW i2j enh f [n, m]

=

hij (NW i2j f [n, m]),

(20)

and the inverse wavelet transform is applied to obtain the denoised-enhanced image. 5. Experimental results We have implemented our method in MATLAB, and tested our approach on the mammograms of the MIAS database (http://www.peipa.essex.ac.uk/info/mias.html). These database images are available in reduced resolution, as compared to conventional digital mammograms. Therefore, we used only two dyadic scales in our analysis. Also, we used G = 3 for all images tested (different degrees of enhancement could be used, but this would make it difficult to compare the results). Next, some of our preliminary experimental results are discussed. It shall be noted that our approach can be used in mammographic images with different number of bits per pixel (e.g. contrast resolutions of 8, 12 or 16 bits/pixel). However, computational complexity increases with the number of bits per pixel, as expected. As mentioned in Section 1, the early detection and treatment of breast cancer are the most effective methods of reducing mortality. Therefore, any early sign of abnormality should be detected in mammograms. Fig. 3(a) shows the MIAS database mammogram 211, containing a cluster of microcalcifications in dense tissues (bright image region), at the bottom right. These microcalcifications are not clearly outlined in a first visual evaluation of the original image (Fig. 3(a)). However, as displayed in Fig. 3(c), after image enhancement the cluster of microcalcifications is more evident in the deep portion of the right breast, suggesting a ductal shape; also, Fig. 3(d) shows that after denoising and enhancement, the shape and other details of this microcalcifications cluster are even more outlined. Fig. 4 shows mammogram 002, presenting macrocalcifications in the upper quadrant on the left breast. It is clear in Fig. 4(d) that contrast enhancement and denoising helps outline better the shapes and boundaries of the macrocalcifications. Fig. 5(a) shows tissues of mammogram 148 that apparently are dense. Fig. 5(c) and (d) shows that our method improves the local contrast, making morphological details of the macrocalcifications and the fibre-glandular tissue are more evident (e.g. the annular shape of the microcalcification in the image lower part, left, also is more evident). These details provide relevant diagnostic information, and were not easily seen in the original image. Fig. 6(a) shows the MIAS database mammogram 145, which presents a nodule in dense tissue (near the image top). Our method helps improving the local contrast, and details of the nodule can be seen more clearly in Fig. 6(c); after denoising, the nodule boundaries and surrounding tissues are shown in detail, as illustrated in Fig. 6(d). Also, Fig. 7 shows the MIAS database mammogram 212. In this case, Fig. 7(c) and (d) shows that image denoising and enhancement help seing details despite of the dense tissue, but probably a spot image still would be needed to perceive more details. In some situations (e.g. large X-ray

253

opaque dense masses), the contrast enhancement provided by the method can be limited, and another mammographic image may be required for a reliable screening (which is a common practice). Next, we compare our contrast enhancement approach with other methods often used, based on our experiments. Histogram equalization often provides excessive tissue contrast, affecting the proportionality among local tissue densities, causing an annoying visual effect. However, gray level normalization improves local contrast, without much distortion of the tissue densities, but has a limited effect on high density tissues. Finally, our approach improves the transparency of high density tissues, but keeping enough contrast to characterize adequately the details and the morphology of calcifications. This comparison is clarified, using mammogram 148 as an example, and the results are shown in Fig. 8.6 Our preliminary experimental results indicate that our method can help denoising and improving local contrast in digital mammograms. These image improvements are important for early detection and risk evaluation posed by microcalcifications and suspicious masses found in such images. 6. Conclusions An important research challenge is to improve visual quality of mammograms by image processing in order for helping early detection of breast cancer. This paper describes new methods for mammographic image preprocessing, for noise suppression and edge enhancement based on the wavelet transform. The image preprocessing was designed to enhance the local contrast in dense regions adaptively. The image denoising process also is adaptive, and the selection of a gain factor provides the desired detail enhancement. Our approach was designed to avoid introducing artifacts in the enhancement process, which is very important when analyzing digital mammograms reliably. The preliminary results indicate that our method improves the detection of microcalcifications and other suspicious structures, even in situations where their detection is difficult (e.g. in low contrast image regions, in dense tissues). Compared to other approaches, our method requires less user adjustment parameters. In the continuation of this work, we intend to perform a thorough evaluation of our approach in clinical trials designed by mammography specialists. Despite of Gaussian noise models being widely accepted, we intend to extend our denoising approach to also include non-Gaussian noise. Acknowledgements The authors wish to thank CNPq (Brazilian Research Council) for financial support, and Professor Betel Blum (Hospital de Cl´ınicas de Porto Alegre, Brasil) for advice and useful discussions. 6 The previous examples only compare our approach with histogram equalization, since this technique is very popular

254

J. Scharcanski, C.R. Jung / Computerized Medical Imaging and Graphics 30 (2006) 243–254

References [1] Laine A, Fan J, Yang W. Wavelets for contrast enhancement of digital mammography. IEEE Eng Med Biol 1995;14(5):536–50. [2] Nunes FLS, Schiabel H, Benatti RH, Stamato RC, Escapinati MC, Goes CE. A method to contrast enhancement of digital dense breast images aimed to detect clustered microcalcifications.. In: Proceedings of 9th IEEE International Conference on Image Processing. Thessaloniki, Greece; 2001. p. 305–8. [3] Karssemeijer N. Adaptive noise equalization and recognition of microcalcification clusters in mammograms. Int J Pattern Recog Artif Intell 1993;7:1357–76. [4] Gurun OO, Fatouros P, Khun GM, de Paredes ES. A controlled phantom study of a noise equalization algorithm for detecting microcalcifications in digital mammograms. Medical Physics 2001;28(4):445–54. [5] McLoughlin KJ, Bones PJ, Karssemeijer N. Noise equalization for detection of microcalcification clusters in direct digital mammogram images. IEEE Transact Med Imaging 2004;23(3):313–20. [6] Simoncelli EP, Adelson EH. Noise removal via Bayesian coring. In: Proceedings of 3rd IEEE International Conference on Image Processing. Lausanne, Switzerland; 1996. p. 171–4. [7] Malfait M, Roose D. Wavelet based image denoising using a Markov random field a priori model. IEEE Transact Image Process 1997;6(4):549–65. [8] Chang SG, Vetterli M. Spatial adaptive wavelet thresholding for image denoising. In: Proceedings of the 1997 International Conference on Image Processing (ICIP97); 1997. p. 374–7. [9] Chang SG, Yu B, Vetterli M. Spatially adaptive wavelet thresholding with context modeling for image denoising. IEEE Transact Image Process 2000;9(9):1522–31. [10] Mihc¸ak MK, Kozintsev I, Ramchandram K, Moulin P. Low-complexity image denoising based on statistical modeling of wavelet coefficients. IEEE Signal Process Lett 1999;6(12):300–3. [11] Strela V, Portilla J, Simoncelli EP. Image denoising via a local Gaussian scale mixture model in the wavelet domain. In: Proceedings of the SPIE 45th Annual Meeting. San Diego, CA; 2000. [12] Portilla J, Simoncelli EP. Image denoising via adjustment of wavelet coefficient magnitude correlation. In: Proceedings of the 7th International Conference on Image Processing. Vancouver, BC, Canada; 2000. [13] Figueiredo MAT, Nowak RD. Wavelet-based image estimation: an empirical bayes approach using Jeffrey’s noninformative prior. IEEE Transact Image Process 2001;10(9):1322–31. [14] Simoncelli EP, Adelson E. Noise removal via bayesian wavelet coring. In: Proceedings of the IEEE International Conference on Image Processing. Lausanne, Switzerland; 1996. p. 279–382. [15] Yoshida H, Zhang W, Cai W, Doi K, Nishikawa RM, Giger ML. Optimizing wavelet transform based on supervised learning for detection of microcalcifications in digital mammograms. In: Proceedings of the IEEE International Conference on Image Processing. Lausanne, Switzerland; 1996. p. 152–5.

[16] Laine A, Schuler S, Fan J, Huda W. Mammographic feature enhancement by multiscale analysis. IEEE Transact Med Imaging 1994;13(4):725– 40. [17] Mallat SG. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transact Pattern Anal Machine Intell 1989;11(7):674–93. [18] Heinlein P, Drexl J, Schneider W. Integrated wavelets for enhancement of microcalfications in digital mammography. IEEE Transact Med Imaging 2003;22(3):402–13. [19] Mallat SG, Zhong S. Characterization of Signals from Multiscale Edges. IEEE Transactions on Pattern Analysis and Machine Intelligence 1992;14(7):710–32. [20] Donoho DL. Nonlinear wavelet methods for recovery of signals, densities and spectra from indirect and noisy data. In: Daubechies I, editor. In: Proceedings of the Symposium on Applied Mathematics. Providence, RI; 1993. [21] Donoho DL. Wavelet shrinkage and W.V.D.: A 10 min tour. Progress Wavelet Anal Appl 1993. [22] Srivastava A, Lee AB, Simoncelli EP, Zhu SC. On advances in statistical modeling of natural images. J Math Imaging Vision (18)2003;17–33. [23] Xu Y, Weaver JB, Healy DM, Lu J. Wavelet transform domain filters: a spatially selective noise filtration technique. IEEE Transact Image Process 1994;3(6):747–58. [24] Bao P, Zhang L. Noise reduction for magnetic resonance images via adaptive multiscale products thresholding. IEEE Transact Medical Imag 2003;22(9):1089–99. [25] Jung CR, Scharcanski J. Wavelet transform approach to adaptive image denoising and enhancement. J Electronic Imag 2004;13(2):278–85.

Jacob Scharcanski has a Ph.D. degree in Systems Design Engineering (University of Waterloo, 1993), a M.Sc. degree in Computer Science (1984) and a B.Eng. in Electrical Engineering (1981), both from the Federal University of Rio Grande do Sul (Brazil). His main areas of interest are image processing and analysis, machine learning and medical applications. He has lectured at U. of Toronto, U. of Guelph, U. of East Anglia and U. Manchester (UMIST), as well as in several Brazilian Universities. He authored and co-authored more than 75 publications in Journals and Conferences. He also held research and development positions in the Brazilian and North-American Industry. Currently, he is an Associate Professor at the Institute of Informatics, Federal University of Rio Grande do Sul, Porto Alegre, RS, Brazil. Cl´audio R. Jung received the B.S. and M.S. degrees in Applied Mathematics, and the Ph.D. degree in Computer Sciences, from Universidade Federal do Rio Grande do Sul, Brazil, in 1993, 1995 and 2002, respectively. He is currently an Assistant Professor at Universidade do Vale do Rio dos Sinos, Brazil. His research interests include image denoising, edge detection and enhancement, wavelets, color image processing, computer vision for intelligent vehicles, aerial image analysis, background subtraction techniques and object tracking.