Combining Multifractal Analyses of Digital Mammograms and Infrared Thermograms to Assist in Early Breast Cancer Diagnosis

. We used a 1D wavelet transform modulus maxima (WTMM) method to analyze the temporal fluctuations of breast skin temperature recorded with an infrared (IR) camera from a panel of patients with breast cancer. This study shows that the multifractal complexity of temperature fluctuations observed in healthy breasts, is lost in the region of the malignant tumor in cancerous breasts. Then, we applied the 2D WTMM method to analyze the spatial fluctuations of breast density in the X-ray mammograms of the same patients. Compared to the correlated roughness fluctuations observed in the healthy areas, some clear loss of correlations is detected in malignant tumor foci. These physiological and architectural changes in the environment of malignant tumors detected in both thermograms and mammograms open new perspectives in computer-aided multifractal methods to assist in early breast cancer diagnosis. our computational approach and its potentiality to assist in early breast cancer diagnosis, progression and treatment.


INTRODUCTION
According to World Health Organization, breast cancer is one of the most common causes of cancer deaths among women worldwide. Clinical studies have demonstrated that patients' survivability increases if breast anomalies are detected at early stages [1]. Although X-ray mammography is currently the most effective way of early breast cancer detection, it remains difficult to interpret due to high tissue variability and to two-dimensional projection effects [2,3]. The miss rate in mammography is increased in dense breasts of young women [4,5]. Criticism of the use of screening mammography due to over-diagnosis made clinicians use additional modalities such as ultrasound, magnetic resonance imaging, IR imaging, etc. Since the original observation by R. Lawson [6] that skin temperature over a malignant tumor is higher than its neighborhood, possibly resulting from abnormal increase in metabolic activity and vascular circulation in the tissues beneath [7][8][9][10], IR thermography has been considered as a promising non-invasive screening method of breast cancer detection [11][12][13]. However, the suitability of IR imaging for routine screening has been highly questioned [14], because of insufficient sensitivity for detection of deep lesions and limited knowledge of the relationship between skin surface temperature distribution and thermal diseases. This growing wave of criticism of breast cancer screening techniques has been strengthened by a similar movement criticizing the increasing use of computer-aided detection (CAD) methods that were shown to decrease specificity [15], and to increase recall rates of healthy women [16][17][18][19] resulting in an increase in unnecessary stress in women. Note that most of these CAD methods [20][21][22][23][24] were designed for texture analysis and feature extraction with the prerequisite that background roughness fluctuations of normal breast mammograms are homogeneous and uncorrelated. In this paper, we review the results of recent and ongoing analyses of dynamic infrared thermograms and X-ray mammograms based on wavelet-based multifractal methodologies [25][26][27].
Application of dynamic IR imaging as a diagnostic tool is based on the detection of variations in temperature rhythms generated by the cardiogenic and vasomotor frequencies [28,29]. Human cardiogenic frequency levels lie between 1 and 1.5 Hz and vasomotor frequencies between 0.1 and 0.2 Hz, with higher harmonic frequency rhythms. But, beyond intensity differences in these rhythms between normal and tumor tissues, there is much more to learn from dynamical IR imaging of breast cancer. Using a 1D wavelet-based multiscale analysis of the temperature intensity fluctuations about these physiologic perfusion oscillations, we propose to characterize the multifractal properties of these temperature fluctuations as a novel discriminating method that can be implemented in early screening procedures to identify women with a high risk of breast cancer [25,26].
The main premise in our investigation of mammograms [27] was related to the suggestion that healthy dense and fatty tissues, indeed display persistent long-range correlations and anti-correlations in mammogram spatial fluctuations, respectively [30,31]. Using a 2D wavelet-based multifractal methodology, we recently observed some clear loss of correlations in mammogram area surrounding the malignant tumor and which strongly correlate with regions of anomalous monofractal (as compared to healthy multifractal) temperature fluctuations previously detected in dynamic IR thermograms [27]. This loss of complexity in temperature time-series and mammographic subimages is likely to be the signature of physiological and architectural alterations of the environment of malignant breast tumors. Besides their potential clinical impact, these results open new perspectives in the investigation and understanding of the cellular mechanisms that underlie the loss of tissue homeostasis, cancer initiation and progression.

The 1D Wavelet Transform (Time-Frequency Analysis)
The WT is a space-scale analysis which consists in expanding a signal in terms of wavelets that are constructed from a single function, the "analyzing wavelet" , by means of translations and dilations. The WT of a function ( ) t is defined as [37]: where 0 t and ( 0) a are the time and scale (inverse of frequency) parameters, respectively. By choosing a wavelet whose (n + 1) first moments are zero ( ( )d 0,0 ), m t t t m n one makes the WT microscope blind to order-n polynomial behavior, a prerequisite for multifractal fluctuations analysis [54][55][56][57].

The 2D Wavelet Transform (Space-Scale Analysis)
With an appropriate choice of the analyzing wavelet, one can reformulate Canny's multiscale edge detection [58] in terms of a 2D wavelet transform [37]. The general idea is to start smoothing the discrete image data by convolving it with a filter and then to compute the gradient of the smoothed image. Let us consider two wavelets that are, respectively, the partial derivatives with respect to x and y of a 2D-smoothing function ( , ) x y : x y x and 2 ( , ) ( , ) . x y x y y The WT of I (x, y) with respect to 1 and 2 has two components and therefore can be expressed in a vectorial form [37]: from which we can extract the modulus and argument of WT:

The Wavelet Transform Modulus Maxima Method
The wavelet transform modulus maxima (WTMM) method consists in computing the WT skeleton defined, at each fixed scale a, by the local maxima L(a) of the WT modulus M (·, a). In 1D, these WTMM are disposed on curves connected across scales called maxima lines which defines the WT skeleton [54-56, 59, 60]. In 2D, these WTMM lie, at a given scale, on connected chains called maxima chains. When considering the points along these maxima chains where M (·, a) is locally maximum, we define the so-called WTMMM L(a). The WTMMM are then linked through scales to form the WT skeleton [31,[61][62][63]. Along these maxima lines the WTMM (resp. WTMMM) behave as a h(t) (resp. a h(x) ), where h(t) (resp. h(x)) is the Hölder exponent characterizing the singularity of (resp. I) at time t (resp. position x). The multifractal formalism amounts to characterize the relative contributions of each Hölder exponent value via the estimate of the D(h) singularity spectrum defined as the fractal dimension of the set of points t (resp. x) where h(t) (resp. h(x)) = h. This spectrum can be obtained by investigating the scaling behavior of partition functions defined in terms of WT coefficients [31,[54][55][56][57]: where q . Then from the scaling function (q), D(h) is obtained by a Legendre transform: In practice, to avoid instabilities in the estimation of the singularity spectrum D(h) through the Legendre transform, we instead compute the following expectation values [31,[54][55][56][57]: where ˆ( , , ) ( , ) ( , ) W q l a M a Z q a are the equivalent of the Bolzmann weight in the analogy that links the multifractal formalism to thermodynamics [56]. Then from the slopes of h(q, a) and D(q, a) vs ln a, we get h(q) and D(q) and therefore the D(h) and, consequently, singularity spectrum as a curve parametrized by q.

Monofractal versus Multifractal Functions
The (q) data are fitted by the quadratic approximation The corresponding singularity spectrum has a characteristic single humped shape is the fractal dimension of the support of singularities of (resp. I), c 1 is the value of h that maximizes D(h) and 2 , c the so-called intermittency coefficient [64], characterizes the width of the D(h) spectrum. Homogeneous (mono) fractal functions, i.e. functions with singularities of unique Hölder exponent H are characterized by a linear (q) curve 1 ( ) . h q c H The D(h) singularity then reduces to a single point 2) as the signature of no change in WTMM statistics across scales [31, 39, 54-57, 61-63, 65, 66]. A nonlinear (q) is the signature of nonhomogeneous functions exhibiting multifractal properties, meaning that the Hölder exponent h(t) [54][55][56] (resp. h(x) [31,[61][62][63]) is a fluctuating quantity that depends on t (resp. x). Then D(h) has a single humped shape 2 ( 0) c that traduces some change in WTMM statistics across scales [31,55,63,65].

Experimental Protocols
In this work [25][26][27], we analyzed the thermograms of 33 patients (from 37 to 83 years old) of Perm Regional Oncological Dispensary with breast cancer (invasive ductal and/or lobular cancer) and 14 healthy volunteers (from 23 to 79 years old) [67]. Tumors (1.2 to 6.5 cm in size) were found at 1 to 12 cm depths. Both breasts of each patient and volunteer were imaged with an IR camera (FLIR SC5000) at 50 Hz frame rate during the 10 min immobile patient phase. Each image set contained 30 000 image frames. The mammographic procedure involved two X-ray images of both compressed breasts of each patient, referred to as the medio-lateral oblique (MLO) and the craniocaudal (CC) views. Spatial resolution of images was 200 pixels per 1 cm. For mammogram analysis, our database includes mammograms collected from 30 of the mentioned 33 patients with breast cancer.

WTMM Multifractal Analysis of Breast Infrared Thermograms
We grouped single-pixel temperature time-series into 8×8 pixel 2 squares covering the entire breast [25,26]. The reported results are averaged over 64 temperature time-series in these 8×8 subareas. We analyzed 1-pixel temperature time-series taken from 8×8 pixels 2 -squares covering patient entire breasts. As expected, these timeseries generally fluctuate at a higher temperature when recorded in the tumor region of a cancerous breast (Fig. 1a) than in a symmetrically localized square on the contralateral unaffected breast (Fig. 1b), as well as on a healthy breast (Fig. 1c). As commonly done for noise signals [38,55] and previously experienced when analyzing rainfall timeseries [68], we applied the WTMM method to the cumulative of these temperature time-series using the secondorder compactly supported version (2) (3) of the Mexican hat [26,69] (hence the singularities with a possible negative Hölder exponent 1 0, h become singularities with cum 0 1 1 h h in the cumulative). We found that the partition functions Z(q, a) (Eq. (6)) display nice scaling properties for 1 to 5, q over a range of time-scales limited to [0.43 s, 2.30 s] for linear regression fit estimates in a logarithmic representation. The (q) so-obtained are well approximated by quadratic spectra (Fig. 1d). For the cancerous breast of patient 20, (q) is nearly linear as quantified by a very small value of the intermittency coefficient  (Fig. 1e). In contrast, the (q) spectra obtained for the contralateral unaffected breast of patient 20 and the right healthy breast of volunteer 14 (Fig. 1d)  respectively. This hallmark of multifractal scaling is confirmed by the computation of the corresponding D(h) spectra in Fig. 1e that both display a single-humped shape [26].
The results of our comparative wavelet-based multifractal analysis of (cumulative) temperature fluctuations over the 33 cancerous breasts, the 32 contralateral unaffected breasts (no right breast for patient 6) and the 28 volunteer healthy breasts are shown in Fig. 2 [25,26]. In Fig. 2a, the corresponding histograms of c 1 values extend over a rather wide range 1 0. 6 1.8, c but turn out to be quite similar with the mean values 1 1.066 0.002 c (cancerous), 1.104 0.002 (contralateral unaffected) and 1.103 0.002 (healthy). In contrast, the intermittency parameter c 2 has a definite discriminatory power. The histogram for cancerous breasts 2 ( 0.045 0.001) c is definitely shifted towards smaller values relative to the ones for contralateral unaffected 2 ( 0.056 0.001) c and normal breasts of healthy volunteers 2 ( 0.058 0.001) c (Fig. 2b). The small values of c 2 dominate in cancerous breasts compared to normal breasts, confirming that the cancerous breasts are enriched in squares where temperature fluctuations display significantly reduced multifractal properties. This justifies that we considered 2 0.03 c as the threshold below (resp. above) which a square was qualified as monofractal (resp. multifractal) (Figs. 3d-3f).
According to the proposed classification criteria, 8×8 pixel 2 squares covering breast thermograms were color coded in the following way: squares with 2 0.03 c were shown in black, squares with 2 0.03 c in white and squares without scaling in gray. A high percentage of squares in the cancerous breast displays monofractal temperature fluctuations with small intermittency coefficient values (49.7% in Fig. 3a), while only few of those squares are found in the contralateral unaffected breast (7.7% in Fig. 3b) and in the healthy volunteer breast (11% in Fig. 3c). Consistently, in the thermograms of the contralateral unaffected and volunteer's normal breasts, a majority of squares (89.4% and 65%, respectively) manifests multifractal scaling with an intermittency coefficient 2 0.03 c (Figs. 3b and 3c). Note that 43.1% of the squares in the cancerous breast also display multifractal temperature fluctuations, as observed for normal breasts (Fig. 3a). These squares indeed cover regions of the breast that are far from the tumor area (left upper quadrant) mostly covered by monofractal square. Let us point out that for each breast, a small percentage (<20%) of (gray) squares was removed from our analysis because of lack of scaling.
A comparative statistical analysis of monofractal and multifractal square distributions in breast thermograms confirmed that the percentage of monofractal (black) squares covering cancerous breasts (26.8 ± 3.5%) is about twice higher than the ones covering contralateral unaffected breasts (13.1 ± 2.3%) and breasts of healthy volunteers (11.3 ± 2.2%). This excess is compensated by a smaller percentage of multifractal (white) squares in cancerous breasts (56.9 ± 4.4%) than in contralateral unaffected breasts (68.5 ± 3.8%) and normal breasts (69.4 ± 4.3%) in healthy volunteers (Table 1).
We have compared the percentages of monofractal squares in cancerous and contralateral breasts [26] and found that 25 patients have higher quantity of monofractal squares in the cancerous breast. 18 cancerous breasts have a percentage of monofractal squares greater than the mean value (26%). Among the other 15 cancerous breasts, 7 have a smaller percentage of monofractal squares but these squares are well localized in the tumor region. The remaining 8 cancerous breasts were detected as false negatives due to small percentage of monofractal squares and their distant location with regard to the tumor region. Four of them correspond to rather deep tumors (6 cm, 7 cm, 8 cm and 12 cm depth) in the breasts with massive adipose tissue that can explain the absence of significant changes in skin temperature dynamics. Importantly, in five of the contralateral unaffected breasts a high (>26%) percentage of monofractal squares was detected which could be an early sign of pathology and require thus further check-ups. As a control, we reproduced this comparative analysis on both normal breasts of the 14 healthy volunteers. Only 4 volunteers manifested high (>26%) percentage of monofractal squares, while in most of them this value did not exceed 10% [26]. Breast-wide multifractal analysis of skin temperature temporal fluctuations. As estimated from the (q) spectra computed with the WTMM method (Fig. 1d), the 8×8 pixels 2 squares covering the entire breast are coded (a-c) and represented as dots in the 1 2 ( , ) c c plane (d-f  To evaluate the efficiency of the test based on the multifractal analysis of breast infrared thermograms, we computed sensitivity (Se) and specificity (Sp). The results reported in Table 2 for our rather small set of patients, namely 76% sensitivity and 86% specificity, are quite encouraging for future applications to larger cohorts of females with breast cancer.

WTMM Multifractal Analysis of X-ray Mammograms
We divided each mammogram view into 360×360 pixel 2 squares spanning 18×18 mm 2 [27]. From the central 256×256 pixels 2 part of the WT skeleton, we computed the partition functions Z(q, a) (Eq. (6)) that were found to display scaling properties over a range of scales 2 1 log 3, a corresponding to [0.7 mm, 2.8 mm], for linear regression fit estimates in a logarithmic representation. Within the range 1 3, q statistical convergence was achieved and the so-obtained (q)-spectrum was found remarkably linear (Fig. 4e) for a great majority of squares in both the cancerous and contralateral unaffected breasts, as the signature of monofractal roughness fluctuations [31,61]. We systematically run our 2D WTMM algorithm for all mammographic views in our database. We confirmed the result of our preliminary study [30] that normal tissue regions display monofractal scaling properties 1 2 ( , 0) c H c as characterized by the Hurst exponent values H < 0.5 in fatty areas which look like anti-persistent self-similar random surfaces, while H > 0.5 in dense areas which exhibit persistent long-range correlations in mammographic roughness fluctuations [31,61]. Importantly, our multiscale method provides an objective way to identify areas with some clear loss of correlations in roughness fluctuations ( 0.5).

H
We extended our waveletbased multifractal analysis of CC and MLO mammographic images to the entire sets of 256×256 pixels 2 squares that cover each breasts of the 30 patients. As illustrated in Figs. 4a-4d for patient 20, we color-coded each of these squares according to the monofractal diagnosis obtained with the 2D WTMM method. Besides confirming the existence of areas of correlated spatial roughness fluctuations, our methodology reveals that peculiar areas of uncorrelated roughness fluctuations are mainly found in the neighborhood of malignant tumors and strongly correlate with the regions of anomalous monofractal (as compared to healthy multifractal) temperature fluctuations previously detected in dynamic IR thermograms (Fig. 3a).
We computed the percentages of monofractal uncorrelated H = 0.5 squares in the mammograms of the two breasts of 30 patients with breast cancer. Good agreement was observed between the percentages of uncorrelated H = 0.5 squares in the CC and the corresponding MLO mammograms of both breasts of each patient (Fig. 5a). Obviously, the percentage of uncorrelated H = 0.5 squares in the cancerous mammograms is much larger than in the opposite breast mammograms. When comparing the two breasts of patients with breast cancer (Fig. 5b), we definitely evidenced some important asymmetry with a significant excess of H = 0.5 squares in the mammograms of the cancerous breast compared to the contralateral unaffected breast. This is a very promising indication that the loss of correlations in the mammographic spatial density fluctuations can be used as an indicator of malignancy. We correlated our findings in IR thermograms and X-ray mammograms of the 30 patients (Fig. 5c).   It is clear that the majority of patients have consistent large numbers of "risky" uncorrelated H = 0.5 squares in mammograms and of "risky" monofractal squares in thermograms of their cancerous breast as compared to their contralateral unaffected breast. We further investigated the spatial distribution of uncorrelated H = 0.5 squares over the entire breast defining a H = 0.5 cluster as an aggregation of H = 0.5 squares sharing a common edge [27]. Isolated H = 0.5 squares were considered as cluster-singlet. As reported in Table 3, the mean size of a cluster in a cancerous breast is larger than in the contralateral unaffected breast. This clustering effect is even more evident for the size of the largest H = 0.5 square that is twice bigger for the cancerous breasts than for the opposite ones. These findings suggest that uncorrelated H = 0.5 squares conglomerate in rather big clusters in cancerous breasts likely defining "risky" area that can be used as a new indicator of malignancy.

CONCLUSIONS
By combining a 1D WTMM multifractal analysis of IR temperature time-series and a 2D WTMM multifractal analysis of X-ray mammograms, we have revealed some change in both temporal and spatial scaling properties in the region of the malignant tumor in cancerous breasts. This is a clear indication of some physiological alterations and architectural disorganization of the "niche" surrounding breast tumor that may be the signature of a transition from an antitumorigenic to a tumorigenic microenvironment underlying the process of oncogenic transformation, tissue invasion and metastasis invasion during cancer progression [70][71][72][73]. Despite the exploratory nature of this preliminary study involving a rather small cohort of females who all went through surgery to remove the histologically confirmed malignant tumor, the obtained results are quite encouraging in the perspective of generalizing our approach to the analysis of longitudinal studies. This will allow us to quantify the statistical performance (sensitivity, specificity) of our computational approach and its potentiality to assist in early breast cancer diagnosis, progression and treatment.

ACKNOWLEDGMENTS
This work was supported by INSERM, ITMO Cancer for its financial support under contract PC201201-084862 "Physiques, mathématiques ou sciences de l'ingénieur appliqués au Cancer", the Russian Foundation for Basic Research (grant No. 16-41-590235) and the Maine Cancer Foundation.
The study reported in this article was conducted according to accepted ethical guidelines involving research in humans and/or animals and was approved by an appropriate institution or national research organization.
The study is compliant with the ethical standards as currently outlined in the Declaration of Helsinki. All individual participants discussed in this study, or for whom any identifying information or image has been presented, have freely given their informed written consent for such information and/or image to be included in the published article.