Image Anal Stereol 2001;20:169-174 Original Research Paper
WAVELET-BASED MULTIFRACTAL FORMALISM TO ASSIST IN DIAGNOSIS IN DIGITIZED MAMMOGRAMS
Pierre Kestener1, Jean Marc Lina2, Philippe Saint-Jean2 and Alain Arneodo1
1Centre de Recherche Paul Pascal, Av. A. Schweitzer, 33600 Pessac cedex, France, 2Centre de Recherche Mathematiques, Universite de Montreal, CP 6128, Suc. centre ville, H3C 3J7 Montreal, Canada e-mail: kestener@crpp.u-bordeaux.fir
(Accepted August 24, 2001)
ABSTRACT
We apply the 2D wavelet transform (WTMM) method to perform a multifractal analysis of digitized mammograms. We show that normal regions display monofractal scaling properties as characterized by the so-called Hurst exponent H = 0.3±0.1in fatty areas which look like antipersistent self-similar random surfaces, while H = 0.65 ±0.1 in dense areas which exibit long-range correlations and possibly multifractal scaling properties. We further demonstrate that the 2D WTMM method provides a very efficient way to detect tumors as well as microcalcifications (MC) which correspond to much stronger singularities than those involved in the background tissue roughness fluctuations. These preliminary results indicate that the texture discriminatory power of the 2D WTMM method may lead to significant improvement in computer-assisted diagnosis in digitized mammograms.
Keywords: breast tissue, fractional Brownian motions, Hurst exponent, image analysis, mammogram, microcalcifications, multifractal formalism, rough surface, scale invariance, wavelet transform.
INTRODUCTION
In the past 20 years, several national mass screening mammography programs, Health Insurance Plan of Greatest New-York (1982) and the Swedish 2-county Program of Mammography Screening for Breast Cancer (1992) have shown that early diagnosis can significantly decrease breast cancer mortality of about 23 to 31 % in women of age 49 to 69. Thus mammo-graphy (X-ray examination) has become the most reliable imaging technique for the early diagnosis of breast cancer which is still the leading cause of cancer-related death in women. Indeed mammography plays a vital role in diagnosis of the disease as well as in pretherapeutic management and in control during and after treatment. But the radiological interpretation of mammograms is a rather difficult task since the mammographic appearance of normal tissue is highly variable. This explains that 10 — 30% of cancers which could have been detected are missed while a high percentage of patients called back at screening turn out not to have cancer. Recently, much research has been devoted to developing reliable computer aided diagnosis (CAD) methods (see Doi et al., 1993, for a general review). Many of these methods are based on multiresolution analysis, global and local thresholding, difference image techniques, stastistical approaches, neural networks, fuzzy logic, and the wavelet transform (WT) and related techniques (Heine et al., 1997; Netsch et al., 1999; Qian et al., 2000). Currently most of these methods are often combined
to detect and classify clusters of microcalcifications (MC) which is an important mammographic sign of early (in situ) breast cancer despite the fact that several benign diseases show MC as well. In the middle nineties, fractal methods have also been applied to the analysis of radiographic images with some success in improving the performances of previous CAD schemes (Priebe et al., 1994; Lefebvre et al., 1995; Thiele et al., 1996). But most of these methods have been intrinsically elaborated on the prerequesite that the background roughness fluctuations of normal breast texture are statistically homogenous (i.e. monofractal) and uncorrelated. Regions that contain statistical aberrations that deviate from this monofractal picture are considered as abnormal regions where tumors or MC are likely to be found. Our goal here is to propose an alternative wavelet-based method to perform multifractal analysis of digitized mammograms. The so-called 2D WTMM method was originally designed to describe statistically the roughness fluctuations of fractal surfaces (Arneodo et al., 2000). In the present work, we briefly summarize how the 2D WTMM method provides an efficient framework to study synthetic and natural fractal images (Arneodo et al., 2000; Decoster et al., 2000; Roux et al., 2000). Then we report the results of a preliminary application of this method to mammogram analysis, which are very encouraging for potential use to classifying tissue types as well as detecting tumors and clustered MC.
169
Kestener P et al: Multifractalformalism in mammography
THE 2D WAVELET TRANSFORM MODULUS MAXIMA (WTMM) METHOD
Most of the fractal methods used to analyze digitized mammograms simply rely upon the estimate of the fractal dimension DF which is related to the so-called Hurst exponent H that characterizes statistically the global roughness of the mammogram landscape (Arneodo et al, 2000). The multifractal formalism accounts for possible fluctuations of the local regularity of a rough surface as defined by the Holder (local roughness) exponent h(r) of the function f(r) whose graph defines the rough surface under study (f(r + l) -f(r) ~ |l|hr, |l| -)• 0). The 2D WTMM method (Arneodo etal, 2000) provides a way to estimate the so-called D(h) singularity spectrum defined as the Hausdorff dimension of the set of points r where the local roughness exponent h(r) is h. This method consists in performing a wavelet-based multiscale Canny edge detection. Let us define two wavelets: i/^(x,_y) = d6(x,y)/dx and \j/2(x,y) = d6(x,y)/dy , where 0(x,y) is a 2D smoothing function well localized around x = y = 0. For any function f(x,y) G L2(R), the WT defined with respect to i/A1 and \j/2 can be expressed as a vector:
TVLfKb a)
v{Te[f](b,a)},
(1)
where Te\f](b,a) = a~2 JJd2r 0(a"1(r - b))f(r). If 0 is just a gaussian 0(r) = exp(—r2/2), then Eq. 1 defines the 2D WT as the gradient vector of f(r) once smoothed by dilated versions 0 (r/a) of this filter. At a given scale a > 0, the WTMM are defined by the positions b where the WT modulus ^v[f\(b:a) is locally maximum in the direction «e^,[f](b,a) of the gradient vector Ty,[f]. When analyzing rough surfaces, these WTMM lie on connected chains (Figs. 1b and 1c). Then the WTMM maxima (WTMMM) are identified as the local maxima of M along the WTMM chains. These WTMMM are disposed along connected curves across scales. The WT skeleton (Fig. 1d) defined by these maxima lines contains a priori all the information about the hierarchical organization of the singularities of the function f(r). In particular, one can prove that, provided the first n₯ moments of \j/ are zero, then .Mw\f\ ~ ahr 0 along a maxima line pointing to the point r0 in the limit a —>• 0+, where h(r0) (< n₯) is the local Holder exponent off. The 2D WTMM method consists in defining the following partition functions directly from the WTMMM that belong to the WT skeleton:
a)    =       £     I       sup      .J(w\f\{x,a')
&e&(a)  \(x,a)E&,a<a                          /
(2)
where q G R. Note that the scaling exponent x{q) has well-known meaning for some specific values of q: (i) — t(0) is the fractal dimension of the set of singularities off; (ii) t(1) is related to the fractal dimension of the rough surface § (z = f(x, y)): dfS) = max(2,1 — t(1)); (iii) t(2) is related to the scaling exponent ί of the spectral density S(k) = |f ˆ(k)|2 ~ |k|-cwith/3=4 + T(2).
From the deep analogy that links the multifractal formalism to thermodynamics, the D(h) singularity spectrum can be determined from the Legendre transform of T(q): D(h) = minq[qh — T(q)). Homogenous (monofractal) self-affine functions involve singularities of unique Holder exponent h = dxjdq, i.e. the %{q) spectrum is a linear function of q. On the contrary, a nonlinear %{q) curve is the signature of nonhomogenous functions that display multifractal properties, in the sense that the Holder exponent h(r) is a fluctuating quantity that depends upon the spatial position r.
From a practical point of view, one often prefers to avoid Legendre transforming. Instead one computes the Boltzmann weights Wy\f\(q,J£,a) =
sup(x,a)e^,a<a-^VLf](x'a/) l^{qa).which yield
to the following expectation values:
h(q,a) D(q,a)
1    ln
sup      e^V,[f](x,a)
(x,a')e.Ža<a
xWv[f](q,3?,a), X    Wv[f](q,S?,a)
(3)
&ES?(a)
xln(Wv[f\(q,J?7a)) 7    (4)
from which one extracts h{q) = lima_^0h{q7a)jlna and D(q) = lima^0+D(q,a)/lna, and therefore the D(h) singularity spectrum by eliminating q.
TEST APPLICATIONS OF THE 2D WTMM METHOD TO ISOTROPIC FRACTIONAL BROWNIAN SURFACES
In the eighties, the fractional Brownian motion (fBm) has become a very popular model in signal and image processing (Peitgen et al, 1987). 2D fBm BH(r) indexed by H g]0, 1[ are Gaussian stochastic processes with stationary increments that can be used to generate random self-affine surfaces with known
170
Image Anal Stereol 2001;20:169-174
statistical properties: z(q) = qH — 2. By Legendre transforming T(q), one gets the following D(h) singularity spectrum: D{h) = 2 if h = H, and D(h) = —°° for any other h. Indeed one can prove that almost
_20   .....................               h....................1        0.
12345                12345               0                  0.5
log2(a)                                             log2(a)                                                  h
Fig. 2. Determination of the z(q) and Dih) spectra of 2D fBm with the 2D WTMM method. a) log23?(q,a) vs log2 a; the solid lines correspond to the theoretical spectrum z(q) = qH — 2 with H = 1/3. b) h(q,a) vs log2a; the solid lines correspond to the theoretical slope H = 1 /3. c) z(q) vs qfor H = 1 /3 (•), 1/2 (m) and 2/3 (k); the solid lines correspond to linear regression fit estimates ofH. d) D(h) vs h as obtained from Eqs. 3 and 4. Results correspond to annealed averaging over 32 (1024x1024) fBm images. a is expressed in oW units.
In Figs. 1 and 2 are reported the results of a test application of the 2D WTMM method to 32 (1024 x 1024) realizations of a 2D fBm withH = 1/3, when using a first-order analyzing wavelet (Arneodo et al, 2000). In Fig. 1 is illustrated the computation
all realizations of the fBm process are continuous, everywhere non-differentiable, isotropically scale-invariant and monofractal as characterized by a unique Holder exponent h(r) = H1 Vr.
of the maxima chains and the WTMMM for an individual image (Fig. 1a) at two different scales (Figs. 1b and 1c). As seen in Fig. 1b, at a given scale, the maxima chains correspond to well defined edge curves of the original image once smoothed by the filter 0(r/a). In Fig. 1d is illustrated the corresponding WT skeleton defined by the maxima lines obtained by connecting the WTMMM across scales. In Fig. 2 are summarized the results of the computation of the x(q) and D(h) spectra using the Eqs 2 to 4. As shown in Fig. 2a for B1,3, the annealed average partition function 3f(q,a) displays a remarkable scaling behavior over more than 3 octaves when plotted versus a in a logarithmic representation, and this for a rather wide range of values of q G [—4,6]. When proceeding to a linear regression fit of the data over the first two octaves, one gets the x{q) spectra shown in Fig. 2c for the three values of the fBm index H = 1/3 (anticorrelated), 1/2 (uncorrelated) and 2/3 (positively correlated increments). Whatever H, the data systematically fall on a straight line, the signature of homogenous (monofractal) scaling properties. These results are corroborated by the computation of h(q1a) in Fig. 2b according to Eq. 3. In the limit of small a values, the slope of h(q,a) vs log2(a) is found not to depend on q: h(q) = H = 1/3,Mq. In Fig. 2d are reported the corresponding estimates of D(h) using Eqs. 3 and 4. Independently of the value of q G [—4,6], B1 ,3, B1j2 and B2,3 display the same behavior, namely D(h = H) = 2.00 ± 0.02. The results obtained for both the r(q) and D(h) spectra
Fig. 1. 2D wavelet transform analysis of BH^1j3{r).  9 is an isotropic 2D-Gaussian function. a) 32 grey-scale coding of
the central 512x512 portion of the original image. In b) a=2aW and c) a=201oW, are shown the maxima chains; the local maxima of',/M-y along these chains are indicated by dots (•); from each dot originates an arrow whose length is proportional to ^$y and its direction (with respect to the x-axis) is given by s>/y. Inb), the smoothed image 0b a*B1,3 (Eq. 1) is shown as a grey-scale coded background from white (min) to black (max). In d) is shown the corresponding wavelet transform skeleton defined by the set of maxima lines obtained after linking the WTMMM across scales. oW=13 (pixels) is the characteristic size ofw at the smallest resolved scale.
171
Kestener P et al: Multifractalformalism in mammography
are thus in remarkable agreement with the theoretical predictions. We refer the reader to Decoster et al. (2000), for similar applications of the 2D WTMM methodology to synthetic random multifractal rough surfaces.
APPLICATION OF THE 2D WTMM METHOD TO DIGITIZED MAMMOGRAMS
As we want to study scaling properties of digitized mammograms, we chose to use full-breast images from the Digital Database for Screening Mammography (DDSM) project (Heath et al, 1998) which provides online more than 2600 studies (http:// marathon.csee.usf.edu/Mammography/Database.html) sorted into three categories: normal, cancer, benign. Mammograms were digitized using a 12 bits scanner with both a good spatial resolution of 43.5mm. Full-breast images enable us to select about ten overlapping 1024 x 1024 pixels squares; indeed, in order to master edge effects, only cores of the images were used for the computation of the WT skeleton and partition functions.
Mammographic tissue classification: dense and fatty tissues
Many statistical studies devoted to mammography analysis actually use fractal techniques or models. Our aim here is to analyze normal mammary parenchyma with our multifractal 2D WTMM method. We have selected a set of 10 images in the DDSM database according to ACR breast density rating with some index ranking from 1 to 4, as assigned by an experienced mammographer: 5 fatty (rated 1 on ACR density scale) and 5 dense (rated 4) breasts. The main steps of the 2D WTMM computations are illustrated in Fig. 3 on two (1024 x 1024 pixels) images selected respectively in a dense-glandular and in a fatty breasts. Figs. 3a and 3d show the original (1024 x 1024) images cut out of these mammograms. The corresponding smoothed images and WT maxima chains computed at the scale a = 39 pixels are shown in Figs. 3b and 3e, respectively. Figs. 3d and 3f represent, at a smaller scale, the location of the WTMMM (•) from which originates an arrow which represents the WT vector Ty[f](b,a). In Fig. 4 are reported the results of the computation of the partition functions ^(q,a),h(q,a) and D(q,a) obtained when averaging over 12 overlapping (1024 x 1024) images cut out of the original dense and fatty mammograms. As shown in Figs. 4a and 4b, both dense and fatty tissues display rather good scaling properties over two and a half octaves. The scaling
actually deteriorates progressively when considering large scales, due to finite size effects. When proceeding to a linear regression fit of log2(3?(q1a)) vs log2(a) over the range of scales extending from amin = 1.6s W to amax = 4s W, one obtains the tq) spectra reported in Fig. 4c. From a simple visual inspection, one realizes that dense and fatty breast tissues display quite different scaling properties. The latter presents a tq) spectrum which is remarkably linear in the range q G [-3,3] with a slope H = 0.25 ± 0.05, while the former presents a larger slope H = 0.65 ± 0.05 with some possible non-linear departure which might be the signature of multifractality. This monofractal vs multifractal discrimination between fatty and dense breast tissues is also evidenced by the computation of the corresponding D(h) singularity spectra in Fig. 4d. However the multifractal diagnosis for dense tissues requires further numerical analysis to ensure statistical convergence of the tq) exponents for large values of \q\. Nevertheless, what seems to be robust, considering the whole set of processed images, is the fact that fatty tissues display monofractal scaling behavior with a Hurst exponent H taking value in the range [0.20,0.45] as the signature of anti-persistent roughness fluctuations while dense tissue display (possibly multifractal) scaling with H G [0.55,0.75] as the signature of persistent long-range correlations.
Fig. 3. 2D wavelet transform analysis of 2 mammograms: a) dense breast tissue and d) fatty breast tissue. q is an isotropic 2D-Gaussian function. In b) and e) is shown the WT modulus at scale a = 3s W with the same grey level coding as in Fig. 1b; the maxima chains are shown for comparison. In c) and f) only the maxima chains and the local maxima of \My along these chains are represented (•) at scale a=2.5s W.
172
Image Anal Stereol 2001;20:169-174
(a>	6 ^ 5 4	(b)     .••: .v • •V •• • • • •. • • • . • • • •••¦ • • • • • • • . • • * * •  «noOOno ¦ • mr,m°o*~a m  1°  • no00o° ¦ _o.   8o°  o0o0 _ • 8°°*0o°o0° 8%S° o0co°°°°. o8 o°° °» 0°°   Oo°oo0 . °°° ° o°°° o o°	0 h -5	¦(c)	f^\
' °S». °°r> •.. °8«y __V...,					
			4 3 2 O 0	-3	0    3    6 1
				(d)         ; ---°—^^g —	
2    3    4    5 log-(a)		2    3    4    5 logp(a)		)	0.5 1 h
Fig. 4. Determination of the tq) and D(h) spectra of dense (•) and fatty (o) breasts with the 2D WTMM method. a) log2^{q,a) vs log2a. b) h(q,a) vs log2a. c) tq) vs q. d) D(h) vs h obtained from Eqs. 3 and 4. Same analyzing wavelet as in Fig. 3. These results correspond to annealed averaging over 12 (1024x1024) squares cut out of full-breast images. a is expressed in s W units.
singularities are seen by our mathematical microscope as Dirac singularities; thus the corresponding maxima lines pointing to the MC are likely to display scaling properties with a local Holder exponent h = — 1 (¦y[f\ ~ a~1) down to scales of the order of the MC size where one should observe a cross-over to the value h = 0 ( y\f\ ~ cst) as the signature of the discontinuity induced by the MC boundary. The behavior of the WT modulus along several maxima lines pointing to background points and to MC is illustrated in Fig. 5b. One can thus perform a classification of these lines according to the behavior of ^y[f\ along these lines, and then separate MC (h ~ — 1) from dense background tissue (h ~ 0.65 ± 0.05). Figs. 5c and 5d show the maxima chains that are found to correspond to MC at two different scales. We see that these maxima chains can be used not only to detect MC at the smallest resolved scale (Fig. 5c), but also to perform MC clustering when investigating largest scales (Fig. 5d). Work in this direction is in current progress.
Detecting microcalcifications through WT skeleton segmentation
The presence of clustered MC is one of the most important and sometimes the only sign of cancer in a mammogram. As a potential computer aided-diagnosis tool, let us show how our WT methodology can identify MC which are small calcium deposits in tissue, appearing as clusters of bright spots. Fig. 5 shows how one can actually detect MC by inspecting the WT maxima chains. Indeed, at the smallest scale resolved by our WT microscope (s W = 13 pixels), MC which can be considered as strong singularities, are contour-shaped by some maxima chains. Since the average size of MC is about 200 mm (5 pixels), these
CONCLUSION
We have presented a new space-scale methodology for studying, within the same algorithmic framework, background tissue properties and abnormal singularities associated with breast cancer. For its ability to revealing and distinguishing persistent and anti-persitent long-range correlations, the 2D WTMM method looks very promising in classifying tissues by quantifying breast density in an very accurate way. Furthermore, we plan to improve detection and segmentation of MC by mixing and combining the 2D WTMM method with neural networks techniques, with the ultimate goal of discriminating benignancy from malignancy.
(a)
(c;
c,  V
(d)

V
^71
Fig. 5. Detection and characterization of microcalcifications. a) Original 726x 726 image of dense breast tissue containing MC. b) Scaling behavior of the WT modulus ,/M-y along some maxima lines pointing towards dense tissue background (a) and microcalcifications (a). The dashed (resp. solid) straight line corresponds to the slope h=0.65 (resp. —1) characteristic of background tissue roughness fluctuations (resp. MC). (c) and (d) show the maxima chains obtained after eliminating background tissue maxima chains at scales a = s W c) and 2.5s W d), using the WT skeleton space-scale information.
173
Kestener P et al: Multifractal formalism in mammography
ACKNOWLEDGMENTS
We are very grateful to Dilhuydy MH (Institut Bergonie΄, Bordeaux, France) and to Lalonde L (CHUM, Montre΄al, Canada) for stimulating discussions and for the permission to use some of their mammographic data.
REFERENCES
Arneodo A, Decoster N, Roux SG (2000). A wavelet-based method for multifractal image analysis. I. Methodology and test applications on isotropic and anisotropic random rough surfaces. Eur Phys J B 15:567-600.
Decoster N, Roux SG, Arneodo A (2000). A wavelet-based method for multifractal image analysis. II. Applications to synthetic multifractal rough surfaces. Eur Phys J B 15:739-64.
Doi K, Giger ML, Nishikawa RM, Hoffmann KR, MacMahon H, Schmidt RA, Chua KG (1993). Digital Radiography. Acta Radiol 34:426-39.
Heath M, Bowyer KW, Kopans D (1998). In Karssemeijer N, ed. Digital Mammography. Dordrecht: Kluwer Academic Publishers, 457-60.
Heine JJ, Deans SR, Cullers DK, Stauduhar R, Clarke LP (1997). Multiresolution statistical analysis of high resolution digital mammograms. IEEE Trans Med Imaging 16:503-15.
Lefebvre F, Benali H, Gilles R, Kahn E, Di Paola R (1995). A fractal approach to the segmentation of microcalcifications in digital mammograms. Med Phys 22:381-90.
Netsch T, Peitgen HO (1999). Scale-space signatures for the detection of clustered microcalcifications in digital mammograms. IEEE Trans Med Imaging 18:774-86.
Peitgen HO, Saupe D, eds. (1987). The Science of Fractal Images. New York: Springer Verlag.
Priebe CE, Solka JL, Lorey RA, Rogers GW, Poston WL, Kallergi M, Qian W, Clarke LP, Clark RA (1994). The application of fractal analysis to mammographic tissue classification. Cancer Lett 77:183-9.
Qian W, Li L, Sun X, Clark RA (2000). Wavelet-based image processing for digital mammography. In: Aldroubi A, Laine AF, Unser MA, eds. Wavelets Applications in Signal and Image Processing VIII. SPIE San Diego, CA, 596-604.
Roux SG, Arneodo A, Decoster N (2000). A wavelet-based method for multifractal image analysis. III. Applications to high resolution satellites images of cloud structure. Eur Phys B 15:765-86.
Thiele DL, Kimme-Smith C, Johnson TD, McCombs M, Bassett LW (1996). Using tissue texture surrounding calcification clusters to predict benign vs malignant outcomes. Med Phys 23:549-55.
174