Texture Features from Mammographic Images and Risk of Breast Cancer

Whether defined by the Wolfe parenchymal pattern or measures of percent breast density, the radiographic appearance of the breast is the strongest risk factor for breast cancer after age, atypia on breast biopsy, and inherited mutations in BRCA1 and BRCA2 (1-4). Wolfe (5) defined four parenchymal patterns: N1, characterized by a fatty breast with only very small amounts of dysplasia (areas of increased density) and no visible duct formation; P1, which consisted of the N1 fatty breast with ducts occupying up to 25% of the breast volume; P2, characterized by a more severe linear or nodular ductal pattern, which occupied >25% of the breast volume; and DY, which had no prominent ducts, but significant densities or dysplasia, usually comprising 50% to 75% of the breast volume. Studies have consistently shown that women with the P2 and DY categories have significantly greater risk of breast cancer than women with N1 or P1 categories (6).

Although the Wolfe pattern may be considered as a texture-based system, it is subjective, which limits its usefulness. To more reliably and accurately quantify the mammographic pattern of the breast, mammographic breast densities have been expressed as a quantitative percentage of total breast area or volume, and computerized techniques were developed to quantitate breast percent density (PD; refs. 7, 8). PD has been well characterized (7), is reproducible (6, 7, 9), and consistently shows a strong association with breast cancer risk (6). The Wolfe pattern predicted breast cancer risk as well as PD in one (10) but not a second (11) study.

Although PD is an established risk factor for breast cancer, it is not sufficiently strong enough to use clinically for risk prediction (12, 13). One strategy to improve risk prediction is to incorporate additional aspects of the mammogram image. The observation that the Wolfe pattern does not correlate perfectly with PD suggests that there is additional information contained within the image observable by radiologists that is not being captured by PD alone. Image texture has long been used in image analysis for segmentation and classification, and multiple approaches to estimate and quantify texture exist. First-order features are those that can be calculated from the histogram of the values in the region of interest. These include mean value, variance, skewness, kurtosis, gray level percentile, and entropy (14). Such features, because they are based on the histogram, ignore information about the spatial relationships of the pixels and consider only their intensities. Second-order features do consider spatial relationships. Perhaps the most commonly used are Markovian features, based on the use of cooccurrence matrices, which measure the probability that a pixel of a certain gray level will be positioned at a particular distance and orientation from a pixel of a certain other gray level (15). Other second-order features include run-length measures, which evaluate texture by measuring “runs” of consecutive pixels with similar gray levels and calculating statistics based on the number and length of such runs (16). Laws features are based on “microtexture masks,” which are specifically designed to detect patterns such as edges or corners at specific orientations (17). Wavelet features measure energy from features at different scales and orientations (18, 19). Finally, Fourier techniques are based on the frequency spectrum of an image or region and can be used to identify preferred periodicities or orientations in a textural pattern or simply to evaluate high-frequency versus low-frequency content (20).

Although several studies have assessed some of these features, we undertook a more systematic and comprehensive assessment of texture features using a large, community-based screening population with the aim of finding features that may capture an element of risk in addition to breast density. We implemented automated routines to calculate a large number of textural features on mammograms from each class described above and trained these on 123 cases and 258 controls (“training” data set). We next examined the strongest features using an additional 123 cases and 264 controls (“validation” data set), and for features that validated, we fit final models in the combined data set. All mammograms were prediagnostic screening mammograms and were available on average 3.7 years before breast cancer diagnosis in the cases.

The study population has been described previously (9). Briefly, subjects were selected from the Mayo Clinic mammography screening practice in Rochester, Minnesota. Patients who did not provide research authorization for medical record studies (3.6%) were not eligible. Cases (n = 373) were women aged 50 y and older who were first diagnosed with breast cancer between 1997 and 2001. Each woman had at least two prior screening mammograms done 2 y before diagnosis and lived within a 120-mile radius of the clinic. Women who had bilateral mastectomies or breast implants before diagnosis were excluded. Controls were also selected from the screening practice, among women with no prior history of breast cancer. Two controls were individually matched to each case on age (within 5 y), final screening exam date (within 4 mo), menopausal status at final exam date, time between initial and final mammogram (within 8 mo), number of prior screening mammograms (within one mammogram), and residence (county).

Risk factor data, including weight, height, and use of hormone replacement therapy, were obtained from a clinical database or abstracted from the medical record for the dates closest to a particular mammogram. Height and weight were used to construct body mass index (BMI) in kilograms per meters squared. All mammograms during the preceding 10-y period were obtained.

For these analyses, we retrieved the earliest mammogram available on all cases and controls and digitized them at high resolution. Based on the need to digitize at high resolution, mammograms were available on 246 of the cases and 522 of the controls. The analysis was based on film mammograms. All mammograms were digitized on a Lumiscan 85 scanner with 12-bit grayscale depth. The pixel size was approximately 0.100 × 0.100 mm². Most of the images were 18 × 24 cm, yielding images of approximately 1,800 × 2,400 pixels, whereas ∼15% were 24 × 30 cm, yielding images of approximately 2,400 × 3,000 pixels. All four views (left and right mediolateral oblique view and left and right craniocaudal view) were digitized.

Percent breast density (dense area divided by total area × 100) and absolute dense area (cm²) were estimated for each view using a computer-assisted thresholding program (Cumulus) developed specifically to quantify breast density (7). Briefly, two thresholds are set by a single trained technician: one separates the breast from the background and the second separates dense from nondense tissue. We consistently showed high intrareader reliability (r > 0.90) for our technician while reading over 500 duplicate images from this study across varying time frames (9).

Batch files comprised both cases and controls with randomly assigned views and sides within a woman to maximize precision of PD estimates (21). A 5% repeat set of images was included within each batch file for assessment of reliability.

In imaging terms, texture can be described as the spatial arrangement and variation of intensities (gray values) within an image. We assessed the textural properties of the breast tissue using five broad families of features: Markovian cooccurrence matrices, run-length analysis, Laws features, wavelet decomposition, and Fourier analysis. The skewness and kurtosis of the image histograms were also computed. The breast was isolated from the image by using the thresholds and delineations that were defined in the estimation of PD, and texture computations were limited strictly to the breast tissue.

The breast is compressed before image acquisition, which creates two different regions within the breast. The first region consists of the central area of the breast where the thickness is nearly uniform and is referred to as the constant thickness region (CTR). The other consists of tissue near the edge of the breast where the thickness gradually tapers due to the breast geometry. The CTR was approximated by eliminating a margin 160 pixels wide from the perimeter of the breast region. Markovian, run-length, Laws, and wavelet features were calculated for this approximate CTR. A rectangular region within the CTR was also located with a fully automated method developed previously (22). This method finds the largest rectangular box that can be inscribed within the breast region (CTR-box), and this was used to calculate the Fourier texture features, which operate only on rectangular regions. The distribution of total breast area did not vary significantly between cases (138.8 ± 54.4 cm²) and controls (139.5 ± 52.4 cm²; P = 0.68), and this was true for the CTR and the CTR-box.

Markovian Features. Markovian texture measurements quantify image texture by analyzing how frequently pixels of given intensities appear at specific distances from pixels of other intensities. Computation of Markovian features begins by first constructing a cooccurrence matrix. The a,b entry of the cooccurrence matrix specifies the probability that pixels of intensity a and b appear in the image separated by a distance d in either the x or y direction. The distance d is specified before construction of the matrix based on the spatial scale to be explored. To reduce the effects of noise variations in image intensity, the image may be decimated in gray scale by a specified factor g. The decimation is done by dividing the intensity of each pixel by g and retaining the integer portion of the result. Assuming the original image has N gray levels, the resulting image will have N/g gray levels. The gray level reduction is fixed before construction of the matrix based on the expected strength of the texture signal and the degree of noise suppression desired.

Markovian texture features are then computed based on measurements derived from the cooccurrence matrix. Pressman (23) described 21 texture features that could be computed from cooccurrence matrices. All 21 features were used in this study and are referred to as Markov 1-21. Seven additional related features proposed in the literature by various authors or developed here were also calculated and are denoted as Markov 22-28 (24, 25). Two additional features derived from cooccurrence matrices were described in a study of mammographic density patterns in women at high risk of developing breast cancer (26). These two features, termed coarseness and contrast in the latter study, are referred to in this article as Markov 29 and Markov 30, respectively (note that contrast is also the common term for Markov 2, but it is calculated with a different formula).

The Markovian features were computed for gray-scale reduction values (g) of 1, 2, 4, 8, 16, 32, and 64 and for distance values (d) of 1, 2, 4, 8, and 16. These reduction values were selected to cover a wide range of spatial scales and feature strengths relative to noise. Thus, a total of 7*5*30 = 1,050 Markovian texture features were generated for each image.

Run-Length Features. As the name implies, run-length texture features examine runs of similar gray values in an image. Runs may be labeled according to their length, gray value, and direction (either horizontal or vertical). Long runs of the same gray value correspond to coarser textures, whereas shorter runs correspond to finer textures.

Run information from an image is recorded in the run distribution matrix. The a,b entry of the matrix specifies the number of runs of length b with gray value a that occur in the image. Similar to the Markovian texture features, the number of gray levels may be decimated to reduce noise variations in image intensity. Both the gray level decimation and run direction are fixed before construction of the matrix.

Texture content is quantified by computing eight features derived from the run-length distribution matrix. The first five features, described by Galloway (16), examine the length and uniformity of the runs and are invariant with respect to intensity. These features are referred to as RL1 through RL5. Two additional features were computed that weight the texture content according to pixel intensity and are referred to as RL6 and RL7 (27). A final feature, referred to as RL8, simply counts the number of runs within the breast region.

Run-length textures were computed for gray level reduction values of 1, 2, 4, 8, 16, 32, and 64 to explore a range of feature strengths relative to the noise level. Runs were tallied and analyzed separately for the horizontal and vertical directions. Thus, a total of 7*2*8 = 112 run-length features were generated for each image.

Laws Features. Laws features (17) are constructed from a set of five one-dimensional filters, each designed to respond to a different type of structure in the image. These one-dimensional filters are denoted E5 (edges), S5 (spots), R5 (ripples), W5 (waves), and L5 (low pass, or average gray value). By applying a one-dimensional filter in the horizontal direction followed by a (possibly different) one-dimensional filter in the vertical direction, results from 25 different two-dimensional filters can be calculated. These filters were applied to the image at full resolution and to the image downsized by factors of 2, 4, 8, and 16 using cubic spline interpolation. The downsized images allow the Laws filters to detect increasingly broader and coarser textures. The second and fourth moments within the CTR for each of the 25 filters were then calculated. Thus, a total of 25*5*2 = 250 Laws features were calculated for each image.

Wavelet Features. The discrete wavelet transform iteratively decomposes an image into four components based on frequency content and orientation (18, 19). Each iteration splits the image both horizontally and vertically into low-frequency (low pass) and high-frequency (high pass) components. Thus, four components are generated: a high-pass/high-pass component consisting of mostly diagonal structure, a high-pass/low-pass component consisting mostly of vertical structures, a low-pass/high-pass component consisting mostly of horizontal structure, and a low-pass/low-pass component that represents a blurred version of the original image. Subsequent iterations then repeat the decomposition on the low-pass/low-pass component from the previous iteration. These subsequent iterations highlight broader diagonal, vertical, and horizontal textures. The frequency content of a given iteration, termed an octave, is one half that of the previous iteration.

Five iterations of the symlet 12 wavelet transform (18) were done, yielding 15 diagonal, horizontal, and vertical components across five octaves. Each component was then isolated in the image by setting the other components to zero and inverting the transform, resulting in a filtered version of the original. The second and fourth moments of the intensities within the CTR region were then computed. Thus, a total of 5*2*3 = 30 wavelet features were generated for each image.

Fourier Features. The two-dimensional Fourier transform decomposes an image into the sum of sinusoidal waves of increasing frequency. A Fourier transform was done on the CTR-box region for each image. The transform was then broken up into 40 equally spaced annular regions (discarding the central circle) and also a “corner” region consisting of the region outside of the largest inscribed circle. The energy (second moment) within each region was calculated, resulting in 41 Fourier features for each image. The innermost annular region was labeled region 1 and corresponds to the lowest frequency content of the image. Higher numbered regions represent annuli closer to the outside of the transform and correspond to higher frequency content in the image. Large energy values for lower-numbered regions indicate the presence of broader, coarser textures, whereas large energy values for higher-numbered regions correspond to finer textures. These quantities were calculated both as raw numbers and as relative numbers normalized to the total energy in all the annuli adding to unity (22).

Frequency content was also summarized by estimating the power law spectrum of the Fourier transform as previously suggested (28). The power law spectrum relates the power (energy, or second moment) of the Fourier transform to frequency by finding a coefficient β such that Power ∝1/(Radial Freq)^β. The coefficient β was estimated using ordinary least-squares regression. Comparisons of previous work in mammographic spectral analysis show that the power law model holds for mammograms acquired from various sensors and different data representations (22).

Histogram Measurements. In addition to the textural features specified above, the sample skewness and kurtosis of the image histograms were also computed (20). The computations were limited to the CTR.

Statistical Analyses. The study population was randomly divided into training (123 cases and 258 controls) and validation sets (123 cases and 264 controls) while retaining the matched nature of the data. Although controls had been matched to the cases, unconditional logistic regression was used to test for associations between the mammographic features and breast cancer case-control status. This approach provided a more natural structure under which to do our entire suite of analyses, including the calculation of concordance statistics and repeated bootstrap sampling and analysis. Breaking the matches, and doing analysis using unconditional logistic regression, provided results that were comparable with what would have been observed in the matched analyses because the distributions of the matching variables remained consistent between cases and controls even after the explicit matching was dissolved. We verified the comparability of results from conditional and unconditional logistic regression approaches on selected analyses.

A series of logistic regression–based analyses were carried out to identify features that were associated with breast cancer case-control status. Initial analyses were done within each class of features (i.e., Markovian, run length, Laws, wavelets, and Fourier), and later analyses were done across these classes. Due to the strong correlation among some features within each class, a bootstrap-based variable selection procedure was used to identify the subset of features most strongly associated with breast cancer. Five hundred bootstrap samples of the same size as the original collection of cases and controls in the training set were selected with replacement from the training set. For each bootstrap sample, a stepwise logistic model selection procedure was carried out. This identified the subset of features most strongly associated with breast cancer in the bootstrap sample. Features in this model were recorded and the process was repeated for each bootstrap sample. The number of times each feature appeared in the final stepwise-selected logistic regression model was tabulated.

The features that were most frequently retained in the models specific to each bootstrap sample were considered to be candidate variables for the specific class of image features being considered. These features were then entered into a multiple logistic regression model, which was subsequently simplified using a backward elimination procedure: the least significant features were removed, one at a time, until all variables in the final model were significantly (P < 0.05) associated with breast cancer case-control status. The features identified through these model selection procedures were then assessed in the validation data set. This was done by applying the regression coefficient, estimated from the training data set, from the identified feature to the measurements of that feature from the images in the validation data set. The resulting score was assessed for significance in a simple logistic regression model, and the validation concordance statistic was computed.

For the features that validated within each class, models were then refit on the entire data set. These models were used to obtain estimates of odds ratios (OR), confidence intervals (CI), and concordance statistics. Because age, BMI, and PD are strongly associated with risk of breast cancer, ORs and c-statistics were also estimated after adjusting for these factors. To measure the predictive precision of each feature, the concordance or c-statistic, also known as the area under the receiver operator characteristic (ROC) curve (AUC), for a logistic regression model was calculated. This statistic measures how often the model correctly identifies the case in a random case-control pair as having a higher risk. Values of this statistic range from 0.0 (perfect incorrect prediction) to 0.5 (chance) to 1.0 (perfect prediction).

In an attempt to combine all features across the different classes into one model, a final bootstrap procedure similar to the one used within each class of features described above was done where the final validated features from each class were considered to be the only candidate features. As in the individual classes of features, the top selected features were removed by a backward elimination procedure until only those with P < 0.05 remained. Two-sided P values were obtained for all tests. Statistical Analysis System software (SAS) was used for all analysis.

Table 1 presents descriptive information on the 246 cases and 522 controls included in this report. The time interval between the baseline mammogram and the date of cancer diagnosis in cases or date of last mammogram in controls was 3.7 ± 1.0 years on average (range, 2.0-5.9) and was >3 years for over 73% of the participants. The matching algorithm was quite effective, as evidenced by the similarity of cases and controls with regard to the design variables. The data set was randomly split into training (n = 381, 123 cases and 258 controls) and validation (n = 387, 123 cases and 264 controls) data sets. There were no differences between the two data sets with respect to age or BMI at time of mammography, PD, or time from mammography to diagnosis (cases) or selection (controls; data not shown).

Table 2 shows the main results for PD and the strongest texture features that validated within each class (details of these features are given in the Appendix), and Table 3 shows descriptive characteristics and correlations with age and BMI for each feature. In the combined data set, the unadjusted ORs were in the range of 1.3 to 1.4 per SD change, which is consistent with previous studies. Adjustment for age and BMI either did not change or actually strengthened the ORs and AUCs. AUCs for the age- and BMI-adjusted models of PD and the individual texture features were similar and in the range of 0.60 to 0.62. There were many other features within each feature class that validated and had AUC statistics nearly as large as those listed in Table 2. These features were highly correlated with the top validated feature in each class and were therefore not selected.

As shown in Table 4, the validated texture features were correlated with PD (Spearman correlation coefficients, 0.39-0.76) and with each other (0.46-0.91). When each validated texture feature was individually included in a model with PD, all ORs remained statistically significant, whereas the ORs for PD attenuated from 6% to 14% and the ORs for texture features attenuated from 11% to 18% (Table 5). In these models, the AUC increased only slightly (to 0.62 to 0.63) over models with just the texture feature or PD alone. Finally, in an attempt to combine across the feature families, the bootstrap procedure was repeated using all validated features. The results showed that the Fourier (normalized, region 2) and wavelet (order 5 diagonal, 2nd moment) features were the two most frequently selected. When these two features were fit together in the same model, only the Fourier feature remained statistically significant (OR, 1.31; 95% CI, 1.07-1.60).

The top validated features were all at or near the lowest spatial resolution (coarsest scale) used in these calculations, and a notable general trend for all features of a given class was that the differences between cases and controls (larger intensity variations in the cases) progressively declined toward finer scales. For example, the single best Laws feature measured differences in vertical edge strength at the lowest resolution considered; in this case, the differences were measured at a width of ∼40 pixels. As shown in Fig. 1, for all Laws features, the average AUC progressively decreased from 0.60 to 0.50 with progressively finer scales. Notice also that, as mentioned above, many Laws features have AUC results similar to the single best Laws feature. Similar behavior is seen for the wavelet and Fourier features. The wavelet feature showed differences in diagonally oriented tissue structures at the coarsest scale considered, corresponding to 32 to 64 pixels. The Fourier feature is measuring energy at a very low spatial frequency annulus, which varies with the size of the CTR-box but corresponds roughly to 27 to 40 pixels, and the Fourier power law exponent shows that there is relatively more energy at lower frequencies in the cases compared with controls. Figures 2 and 3(for the wavelet and Fourier features, respectively) show that as one moves to lower wavelet orders or higher frequency bins (i.e., toward finer scales in both cases), the ability to discriminate between cases and controls degrades gradually. The Markovian feature, which measures a variant of image contrast, examined differences 8 pixels apart, which was the coarsest scale considered for this family. The run-length features were calculated only at one scale, so no general statement can be made about scale, but the relationship found indicates a similar trend: the cases have a less uniform image, with shorter runs of adjacent pixels with similar gray levels.

Following a comprehensive evaluation of >1,000 individual measurements from five different classes of texture features, we identified texture features derived from mammographic images that predict breast cancer risk at the same magnitude as PD. When individual features were included in the same model as PD, both the feature and PD remained significant predictors, although the strength of the association (OR) weakened slightly for both variables, likely due to their moderate correlation. Including variables from across feature classes did not improve prediction ability. More importantly, none of these features added significantly to risk prediction beyond PD alone.

Our results consistently show that texture features at low spatial frequencies (i.e., coarser mammographic textures) provided the strongest predictors of future breast cancer risk. The general trend is that the cases seem to have stronger intensity variations (i.e., more energy) across coarse scales than the controls. Many features from different classes have similar ability to discriminate case and control status, and this discrimination degrades as features are calculated at finer scales. Indeed, the results suggest that a wide variety of methods of measuring intensity differences over scales of several mm or so may yield predictive power comparable with PD. We hypothesize that this may primarily be a reflection of the simple presence of breast density itself. Breast density is a measure of what percentage of the breast has intensity above some threshold that is meant to separate dense fibroglandular tissue from darker fatty tissue. Because the breast density for most women is well under 50%, an increase in density usually means that there are more bright areas in an image that is still primarily dark. In this scenario, texture features would show increased energy at low frequencies as density increases because pixel differences would increase at broad spatial scales. This would explain why so many features have predictive power roughly equivalent to PD, but the predictive power does not significantly increase when they are combined with PD. As noted by Pepe et al. (29), for the AUC to be increased substantially by the addition of new risk factors, the magnitude of association for at least one of the added variables needs to be quite large. Unfortunately, the ORs corresponding to the features identified here were not of a magnitude that could lead to a marked improvement in the AUC.

It is also possible that the textural differences between mammogram images from cases and controls may be correlated with the Wolfe classifications. Both the high-risk P2 and DY categories consist of broad, coarse patterns that are likely to exhibit strong energy at lower frequencies. Furthermore, our texture measurements, similar to the Wolfe patterns, offer little discrimination capability above and beyond PD. However, we did not have the Wolfe classification on our study population, and additional study would be required to establish a definitive link between the texture measurements and the Wolfe classifications.

Several previous studies have assessed texture features as predictors of breast cancer risk. Torres-Mejia (10) found a 39% decreased risk of breast cancer per 1 SD increase in lacunarity after adjusting for SD and other confounders. A decrease in lacunarity at a given density corresponds to a more regular or homogeneous arrangement of objects versus gaps in the image (i.e., a more homogeneous mammographic parenchymal pattern) and is associated with a higher risk of breast cancer. Because this analysis was done in conjunction with a fractal dimension analysis and calculated a single lacunarity coefficient across all scales, it is difficult to compare this result to ours in terms of strength of features at different scales.

Byng et al. (30) and Yaffe et al. (31) reported a negative correlation between skewness and cancer risk, but not as strong as PD itself, and probably due to the high correlation between skewness and PD. Nagao et al. (32) studied an automated measure based on the location of the peak of the histogram and also saw a result similar to but not as strong as PD.

Other studies have also used mammographic texture analyses to identify differences among women at different levels of risk for developing breast cancer. Notably, Huo et al. (26) used a linear combination of four mammographic texture and histogram features to distinguish between mammograms of a small sample of BRCA1/BRCA2 mutation carriers (n = 30), of whom the majority (n = 17) also had breast cancer, and age-matched noncarriers (n = 60). The latter study achieved a very high AUC (0.91) for distinguishing mutation carriers from low-risk noncarriers. Both their Markovian features (coarseness and contrast) and their histogram features (skewness and kurtosis) were analyzed in our data set but did not show a significant ability to predict breast cancer risk in our average-risk population. Huo et al. (26) concluded, however, that BRCA1/BRCA2 carriers tend to have mammograms that are low in contrast and are coarse in texture, which is consistent with our findings for discriminating case and control status.

The texture measurements presented here may be useful in a multivariate risk prediction model. Although the predictive capabilities of these features are moderate at best, they are similar to the widely used Gail model, including the more recent models that incorporate mammographic density (12, 13). However, density estimation is both subjective and laborious, which may limit its use in clinical practice. The texture measurements presented in this work, though, are both objective and automated, which is an important advantage. However, numerical values of texture features tend to vary with differences in acquisition, so standardizing across institutions and scanners may be difficult. It might be possible to standardize for variables such as compression force, thickness, kVp, and mAs with techniques such as the representation proposed by Highnam and Brady (33). Accounting for these influences may not only improve predictive capability but may also generalize the measurements to films from other institutions. Unfortunately, acquisition variables were not available for the vast majority of the films in this data set. It will also be important to confirm the significance of these features in the full-field digital mammography environment, where these features could automatically be assessed in the digital data at the time of examination.

Our study design had several key strengths. We used a relatively large case-control study that allowed us to directly assess the association of PD and texture features from screening mammograms with subsequent breast cancer risk. The case-control study included careful control of age, mammogram characteristics (timing and frequency), menopausal status, and residence. We also had data on weight, height, and other breast cancer risk factors at the time of mammogram. We used a quantitative estimate of PD that has consistently shown association with risk (6) and we had excellent intrareader reliability. Our approach to evaluating texture features was comprehensive and systematic, with >1,000 texture features from five different classes analyzed. Our statistical analysis used a robust bootstrapping algorithm that was able to handle the large number of correlated independent variables. Both a split sample validation as well as cross-validation in the training set were used to minimize overly optimistic predictions.

There are some limitations that affect interpretation of our results. Our study sample had limited ethnic diversity, and whether these results will hold in other groups is not known. However, PD seems to be a strong risk factor in several populations, including Asian Americans, African-Americans, and Native Hawaiians (34-36). Although this was not a population-based study, by restricting to the 120-mile radius and requiring serial mammograms on all women, the study population was close to a community-based rather than a referral or high-risk population and the internal validity is high. Finally, all mammograms were acquired from a single institution, and therefore, future studies will need to assess the robustness of these results across different institutions and imaging platforms.

In summary, our results show that breast cancer cases tend to have stronger textural properties at coarse scales in their screening mammograms than controls. The ability to predict individual breast cancer risk from these features was modest, although it is similar to both PD and the widely used Gail model. Further analyses that account for the image acquisition variables may improve predictive ability. The objective and automated nature of these features may make them useful within a multivariate risk prediction model. Finally, this work only shows the relationships with the breast cancer risk but does not attempt to relate the findings to the underlying biology or disease etiology. Such knowledge would be useful in developing better risk prediction models and perhaps identifying novel interventions to reduce risk of breast cancer.

The details of the specific features cited in Tables 2 to 4 are described below. Note, however, that we do not consider these features to be unusually predictive; as discussed above, many other features from each of the classes at coarse scales are almost equally predictive.

The Laws feature (E5L5 M2 scale 16) is the second moment (energy) in an image obtained by decimating the original image by a factor of 16 in each direction and convolving the result with the E5 “edge” filter (1 -2 0 2 1) in the x direction and the L5 “low-pass” filter (1 4 6 4 1) in the y direction. This can be interpreted as a measure of vertical edge strength at a coarse scale.

The Markovian feature (M23 S8 level 1) is calculated from the original intensities (no gray level decimation). We term the quantity we calculate the “gray level difference variance,” and it is simply the variance of the absolute values of the differences between pixels and their neighbors; in this case, neighbors 8 pixels away in x and y. This is related to the classic Markovian quantities M2 and M11 but differs from the “difference variance” feature M10 (23). It is also related to the analysis by Chandrasekaran et al. (25): if p(i) is the histogram whose ith component is the probability that two neighboring pixels have an absolute intensity difference of i, then the mean of this histogram is

The run-length feature (X direction #3 level 64) calculates the run distribution matrix p(i,j) for runs of identical pixels in the x direction after the image intensities are divided by 64 and calculates the gray level nonuniformity GLN from this matrix according to Galloway (16):

The wavelet feature (O5 D M2) applies five iterations of the symlet 12 wavelet transform (18), keeps the values for the diagonal band of this fifth octave, sets all other wavelet coefficients to zero, inverts the wavelet transform, and calculates the 2nd moment of the image intensities within the CTR.

The Fourier feature (norm 2) is the normalized energy in the second annular bin of the Fourier transform, as described in Materials and Methods. The Fourier power law feature is fully described there as well.

No potential conflicts of interest were disclosed.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

We thank Fang Fang Wu who worked tirelessly to digitize and estimate density on thousands of films for this study and Sondra Buehler for secretarial support.