Optimization of volumetric breast density estimation in digital mammograms

All mammograms underwent some steps of preprocessing. First the images were segmented into breast tissue, pectoral muscle (in case of a mediolateral oblique—MLO view) and background (Karssemeijer 1998). Additionally, we performed a thickness correction of the peripheral region. We also determined the Euclidean distance d(r) from each pixel location to the skin-line in the mammographic projection. The distance d(r) is used several times throughout the paper and is needed to define the breast interior.

As described in van Engeland et al, the volume of fibroglandular tissue can be computed from unprocessed (raw) digital mammograms based on a physical model of image acquisition and on the assumption that the breast is composed of two types of tissue: dense tissue and fatty tissue. The dense tissue thickness h_d at each pixel location can be calculated with the following formula:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{h}_{\text{d}}}(r) & = & -\frac{1}{{{\mu}_{\text{d},\text{eff}}}-{{\mu}_{\text{f},\text{eff}}}}\ln \left(\frac{g(r)}{F}\right) \end{array}\end{eqnarray} \tag{ 1 }$$

where, g(r) is the pixel value at position r and F is the pixel value in a fatty tissue reference region where h_d is supposed to be zero. The effective attenuation coefficients, ${{\mu}_{\text{f},\text{eff}}}$ and ${{\mu}_{\text{d},\text{eff}}}$ for fatty and dense tissue respectively, vary with the breast thickness and with the anode / filter combination of the x-ray tube and are described in van Engeland et al.

The total dense tissue volume (VDT) is given by the integral over the projected breast area B:

$$\begin{eqnarray}\begin{array}{*{35}{l}} VDT & = & {{{\int}^{}}_{B}}{{h}_{\text{d}}}(r){{\text{d}}^{2}}r=-\frac{1}{{{\mu}_{\text{d},\text{eff}}}-{{\mu}_{\text{f},\text{eff}}}}{\int}_{B}\ln \left(\frac{g(r)}{F}\right){{\text{d}}^{2}}r \end{array}\end{eqnarray} \tag{ 2 }$$

In this work, we studied three approaches to obtain an estimate for the pixel value that corresponds to the projection of fatty tissue only. Two reference values (F1 and F2) are based on the pixel value distribution in the breast interior, and the third reference value (F3) was calculated by estimating the proportion of dense tissue in the densest location in the breast. Each approach is explained in the following section.

2.3.1. Reference value F1—the maximum pixel value in a small breast interior region.

To obtain the reference value F1, we used the same approach as described in van Engeland et al The pixel value representative for fatty tissue is determined by taking a large quantile (0.99) of the pixel value histogram computed in the breast interior. We used this approach instead of taking the maximum pixel value, because large pixel values may appear in the mammogram due to artifacts or noise. With this approach the breast interior is rather small and depends on the maximum Euclidean distance computed from the breast pixels to the skin-line. To determine this distance, we first calculate the Euclidean distance from all breast pixels (the pixels that are enclosed by the pectoral muscle boundary, the skin-line and the image edge) to the skin-line. Of these distances we take the maximum. In most cases, the point with the maximum distance to the skin is located on the pectoral muscle or breast boundary, but it is noted that depending on the shape of the breast this point may also be located more inward from the boundary. The interior is then defined as the breast pixels that have at least a distance of 0.4 times the maximum Euclidean distance to the skin-line. An example is shown in figure 1(b).

Zoom In Zoom Out Reset image size

2.3.2. Reference value F2—the maximum pixel value in a large breast interior region.

2.3.3. Reference value F3—an estimate from the densest region (minimum pixel value of the mammogram).

When the breast is very dense, the methods above become inaccurate, as it becomes unlikely that the obtained reference value is representative for the projection of only fatty tissue. In this section, we propose a method to obtain a suitable reference value using the densest region in the mammogram. Given the pixel value corresponding to a dense pixel (g_dense) at a certain location r, it is possible to compute the fatty tissue reference value by estimation of the corresponding thickness of dense tissue ${{h}_{\text{dense}}}(r)$. This can be seen if we rewrite formula (1):

$$\begin{eqnarray}\begin{array}{*{35}{l}} F & = & \frac{g(r)}{\exp \left(-\left({{\mu}_{\text{d},\text{eff}}}-{{\mu}_{\text{f},\text{eff}}}\right){{h}_{\text{d}}}(r)\right)} \\ F3 & = & \frac{{{g}_{\text{dense}}}}{\exp \left(-\left({{\mu}_{\text{d},\text{eff}}}-{{\mu}_{\text{f},\text{eff}}}\right){{h}_{\text{dense}}}\right)} \end{array}\end{eqnarray} \tag{ 3 }$$

The reference value computed in this way is denoted by F3. Note that there is no reference region in the mammogram with this pixel value, but that it can be derived when we can estimate h_dense at the location of g_dense. To find the reference pixel value g_dense, we make use of the large breast interior as used for F2 (see figure 1(c)) and determine the minimum using the 0.01 quantile of the pixel value distribution.

To obtain h_dense we have to estimate the amount of fibroglandular tissue that corresponds to the pixel value g_dense. Since the reference value g_dense corresponds to the densest region in the mammogram, we assume that it represents the projection of dense tissue and the layers of subcutaneous fat. Therefore, this pixel value is also representative for the maximum dense tissue fraction (MDTF) in the compressed breast projection, i.e. the maximum thickness of dense tissue in the path of the x-ray beam divided by the thickness of the compressed breast. Based on this idea h_dense is estimated as $MDTF\times H$ with H the compressed breast thickness. To estimate h_dense, the MDTF in the direction of the x-ray beam is needed. As this information is unavailable, we calculate the MDTF in the image plane and assume that this value is representative for the MDTF in the direction of the x-ray beam. To make this approach work, two key assumptions need to be fulfilled: (1) in a first approximation the MDTF is direction independent with respect to a rotation around the anterior-posterior axis and (2) that the MDTF does not change when the breast is compressed. In section 3.5 we investigate to what extent these assumptions are valid by using MRI data. The MDTF in the image plane is calculated by constructing paths that represent plausible x-ray trajectories through the breast when it would be decompressed, rotated by 90 degrees, and compressed again. We used slightly curved instead of linear paths, to take into account that when imaging a compressed breast the x-ray beam will be approximately perpendicular to the skin and the subcutaneous fat layer.

The curved paths are defined by taking points on a circle centered at a point (P) outside of the breast. A schematic overview is shown in figure 2. To define this point P, we calculate the minimum distance between the nipple (N) and the pectoral muscle (M). In CC views, we used the vertical image edge instead of the pectoral boundary. The point P has a distance to the nipple of two times the distance NM and lies on the straight line with the points N and M. The nipple is assumed to be located on the skin-line at the position with the maximum distance to the pectoral muscle or the vertical image edge for MLO views and CC views, respectively.

Zoom In Zoom Out Reset image size

To estimate MDTF in the plane of the mammogram we first classify all pixels as dense or non-dense with a Random Forest classifier on a rescaled 200 μm image, using features similar to the ones used in Kallenberg et al (2011). Then, for each path, we count the number of dense and non-dense pixels. To minimize variation we group paths to bands with a width of 0.2 mm with a sliding window approach and calculate the dense tissue fraction of these bands. The MDTF of a mammogram is the absolute maximum of the dense tissue fractions. To estimate h_dense, we use the MDTF averaged over all mammograms of an examination (left and right breast, and MLO and CC view), to minimize variation. The MDTF of the MLO and the CC view are assumed to be similar due to the rotation invariance of the MDTF, while the MDTF of the left and right breast are assumed to be similar as the amount of breast density and the breast density distribution of the left and right breast are comparable to each other. Only paths with a minimum distance of NP  +  H/8 are included, to exclude paths that are too close to the nipple. Furthermore, we exclude paths that are too close to the pectoral muscle to prevent that inaccuracies in the pectoral muscle segmentation influence the result. We also use a distance of H/8 here. Using the thickness of the compressed breast H we estimate h_dense by computing $MDTF\times H$.

Previous studies showed that the breast density estimates obtained with F1 provide accurate results for non-dense breasts. However, the reference values F2 and F3 were obtained under the assumption that the breast is dense. Hence, it is likely that the breast density estimations obtained with these reference values do not give accurate results in non-dense breasts. In this study, we developed two hybrid approaches that combine the breast density estimations to obtain suitable results for the complete range of breast densities.

The first combination scheme takes only into account results from F1 and F2. Based on the MDTF a threshold t is set. If MDTF is below that threshold, the fibroglandular tissue volume estimate obtained with F1 is used. Otherwise we use results from reference value F2.

The second combination scheme of the breast density estimates involves three parameters (t, a and b) and results from the three reference values are combined. The first parameter is a threshold (t). If MDTF is below that threshold, the fibroglandular tissue volume estimate obtained with F1 is used. In case of a MDTF above that threshold, a linear combination of the other two fibroglandular tissue volume estimations is used, where the weights depend on the MDTF. The combination scheme can be written as

$$\begin{eqnarray}VDT=\left\{\begin{array}{*{35}{l}} VDT1 & \quad \text{if}\,\,MDTF<t \\ (1-w)\times VDT2+w\times VDT3 & \quad \text{if}\,\,MDTF\geqslant t \end{array}\right. \end{eqnarray} \tag{ 4 }$$

with $w=a\times MDTF+b$ and VDT1, VDT2 and VDT3 the dense tissue volumes obtained with the reference values F1, F2 and F3 respectively.

In the evaluation, five fold cross validation was used to obtain optimal values for the parameters. The data set was divided into five groups, or folds, of equal size. To avoid bias all exams of a woman were in the same fold. Then, four folds were used to optimize the parameters, which were then applied on the fifth fold. This process was repeated such that the obtained parameters were once applied to each fold. The parameters were obtained by minimizing the sum of the squared residuals in the comparison of percent density between mammography and MRI averaged per breast.

In this work, we used two approaches to estimate the breast volume. The first one is based on the semi-circle model (Snoeren and Karssemeijer 2004). In that approach, which was applied in van Engeland et al, the interior is assumed to consist of parallel planes, while the peripheral zone is approximated by semi-circles. The thickness as a function of the compressed breast thickness (H) is given by:

$$\begin{eqnarray}h(r)=\left\{\begin{array}{*{35}{l}} 2{{\left[{{(H/2)}^{2}}-{{\left((H/2)-\text{d}(r)\right)}^{2}}\right]}^{1/2}} & \quad \text{if}\,\,d(r)<H/2 \\ H & \quad \text{if}\,\,d(r)\geqslant H/2 \end{array}\right. \end{eqnarray} \tag{ 5 }$$

where d(r) is the Euclidean distance between the breast pixel location and the skin-line.

The second breast volume estimate algorithm used in this study is based on the work of de Groot et al (2014). That method uses a model where the peripheral zone of the breast in the lateral direction, i.e. the region where tissue is not in contact with the compression paddle, is assumed to have half the width of the peripheral zone in the posterior-anterior direction, due to forces related to the attachment of the breast to the chest wall.

To segment the peripheral zone, de Groot performs a sequence of scaling with a factor of two along the rows, erosion, dilation and rescaling of the segmented mammogram. The pixels in the lateral direction that are closer to the skin-line than H/4 are considered to belong to the peripheral zone, while in the posterior-anterior direction pixels closer than H/2 are included. The different steps are shown in figure 3.

Zoom In Zoom Out Reset image size

Since the description proposed by de Groot et al is for CC mammograms only, we adapted this approach for MLO views. We introduced a coordinate system in which the y-axis is defined by the pectoral muscle boundary, and the x-axis is perpendicular to the y-axis and goes through the nipple position. For each point within the breast segmentation, we determined the absolute angle ($|\alpha |$) between the x-axis and the line going through that point and the origin. Points that were closer to the skin-line than

$$\begin{eqnarray}&&\frac{H}{2}\left(1-\frac{|\alpha |}{\pi}\right)\end{eqnarray} \tag{ 6 }$$

were considered to belong to the region that is not in contact with the compression plates. Figure 4 shows a MLO mammogram with the coordinate system. Like in de Groot et al, we used H, obtained from the DICOM header, as thickness for the fully compressed breast region, and a thickness of $H\pi /4$ (average thickness of the peripheral zone) for the peripheral zone. The breast volume is the integral over the breast thickness over the breast area.

Zoom In Zoom Out Reset image size

For evaluation of the methods we used breast MRI scans. First, the MRI scans were segmented into breast tissue and background. Each breast voxel was then labeled as fibroglandular tissue or fatty tissue. For most MRI's this was performed automatically with a previously developed method (Gubern-Mérida et al 2015). Some images were segmented manually. A manual segmentation was necessary, when the automatic segmentation was not accurate enough judged by a radiologist. Manual segmentations were performed by a trained researcher and were reviewed by a radiologist with expertise in breast MRI.

Manual segmentation was done as follows: in the axial view, in every 5–10 slices a contour was drawn around the breast outline using spline interpolation. Using these contours, the breast outline in the remaining slices was computed using a spline surface function to get the breast mask. These masks were checked and in case of an inaccurate segmentation, additional contours were drawn until the masks were accurate enough. Within the breast mask, contours for the segmentation into fibroglandular and fatty tissue were drawn in the same slices as the contours to segment the breast from background. These contours were extended to the whole breast with the method mentioned above. Within this fibroglandular tissue mask a threshold was applied to distinguish fibroglandular tissue from fatty tissue. Voxels outside the fibroglandular tissue mask but inside the breast mask were considered to be fatty. Segmentations were performed on bias field corrected images. The N4 bias field correction algorithm (Tustison et al 2010) was used.

We made use of 202 pairs of DM and MR examinations of 119 women recorded between December 2000 and December 2011 in the Radboudumc. MRI exams and mammograms were performed within two months of each other. All digital mammograms were acquired on GE Senographe systems using standard clinical settings with a non flexible compression paddle. For all examinations, a complete exam consisting of MLO and CC images of the left and right breast was available. We used the unprocessed (raw) data.

The MR examinations were performed on either a 1.5 or a 3 Tesla Siemens scanner (Magnetom Vision, Magnetom Avanto or Magnetom Trio), with a dedicated breast coil (CP Breast Array, Siemens, Erlangen, Germany). The segmentation of breast and fibroglandular tissue was performed on pre-contrast T1-weighted MR volumes without fat saturation. For 159 MRI examinations the automated segmentation was approved by a radiologist, 43 MRI's were segmented manually.

The fibroglandular tissue volume estimations obtained with the three different reference values and the two combined approaches were compared to estimations based on MRI data. Results were averaged over MLO and CC views to obtain an estimate per breast. For an estimation per exam results were averaged over the left and right breast. Pearson's correlation coefficients and 95% confidence intervals (CI) were calculated using the statistical software package R. Because of the log-normal distribution of the data, correlation coefficients were computed after converting the measurements using the natural logarithm (Wang et al 2013, Gubern-Mérida et al 2014).

Mammographic breast volumes obtained with the two different geometrical approaches, the semi-circle model and the adaption of de Groot et al, were compared to breast volumes based on MRI data. Results were averaged over both mammographic views and over the left and right breast to obtain a single score for each exam. The MRI based breast volumes of the left and right breast were averaged as well. Pearson's correlation coefficients and 95% CI were calculated. Based on the performance of the two algorithms, we decided to proceed with one algorithm for further computations.

The percent density estimations obtained with the three different reference values and the two combined estimates were compared to the estimations based on MRI data. Percent density is defined as fibroglandular tissue volume divided by breast volume. We made use of the breast volume estimations that are based on the work of de Groot et al In the original work of van Engeland the semi-circle model was used to obtain the breast volume. For comparison, percent density is given based on the fibroglandular tissue estimations with F1 and the breast volumes estimations of the semi-circle model as well. The estimations were again averaged over MLO and CC views to obtain an estimate per breast. To get an estimation for each exam, estimations were averaged over the left and right breast. Pearson's correlation coefficients and 95% CI were calculated.

To obtain a fatty tissue reference value with the dense tissue reference region approach we assume that the maximum dense tissue fraction (MDTF) in the direction of the mammographic compression (the direction of the x-ray beam) is the same as the MDTF in the image plane. This assumption is based on the idea that the MDTF does not depend on the direction of the measurement and is not affected by deformation of the breast due to compression. To investigate the validity of this idea we used the previously described data set of 202 MRI exams of 119 women. Additional, we included exams of 20 women who underwent MR guided biopsies and had breast density categorized as either BI-RADS 3 or 4. For each of these 20 women two MRI exams were available: an MRI recorded with the breast under compression in the MR guided biopsy procedure and a regular diagnostic MRI. The compression of the biopsied breast was always in the mediolateral direction. All MR images were manually segmented into fibroglandular tissue, fatty tissue and background.

To verify the rotation invariance of the MDTF measurements, we took samples in regular (uncompressed) MRI exams along lines in multiple directions in the coronal planes. Samples in each direction were taken on a regular grid in a square region of $6\times 6\quad \text{c}{{\text{m}}^{2}}$, were the center of the square region projected to the center of the breast. By doing so, we prevented taking samples close to the pectoral boundary or near the breast edge. Samples had a cross section of $2\times 2\,\text{m}{{\text{m}}^{2}}$. In each of these samples the fraction of dense tissue was computed from the segmented MRI. The MDTF was computed as a function of the projection angle by taking the maximum dense tissue fraction in the square region while rotating the breast. We took samples using 36 different angles, ranging between 0 and 175 degrees in steps of five degrees, simulating 36 different projection directions. The mean MDTF and the coefficient of variation when changing the projection angle were computed and displayed with a scatter plot. For this experiment we used both the set of 20 MRI's of uncompressed dense breasts and the 202 segmented MR images. In figure 5(a) a coronal slice of a MRI volume is given. The sample direction is craniocaudal.

Zoom In Zoom Out Reset image size

In a second experiment we looked at the effect of compression using the exams of the patients who underwent MR-guided breast biopsies. We determined the maximum dense tissue fraction in the MRI's of the compressed breasts in two directions, one in the compression direction and one perpendicular to it, both along lines in parallel with the chest wall. These correspond to the directions in which we compute the maximum dense tissue fractions in the mammograms, with the difference that the MRI's are compressed in mediolateral directions while mammograms are compressed in the mediolateral oblique or craniocaudal direction. A visual example is shown in figure 5(b).

Figure 6 shows the results of the five methods for dense tissue volume estimation using MRI as a reference. The Pearson's correlation coefficients for the comparisons with MRI were 0.79, 0.80, 0.73, 0.83 and 0.86 when using F1, F2, F3, the combination of F1 and F2, and the combined estimates of F1–F3, respectively. In table 1 the correlation coefficients are given for the estimations averaged per breast and per exam.

Zoom In Zoom Out Reset image size

The parameters t, a and b that are needed to combine the results were obtained through cross validation. The optimal values for the parameters were determined in each fold and varied only a little. For both approaches, when combining F1 and F2, and when combining all three estimates, the optimal threshold was set to 0.35 (in four of the five folds). Therefore, the fibroglandular tissue volume as estimated with F1 is used if the MDTF  <  0.35. The parameters a and b are needed when combing all three fibroglandular tissue volume estimates. In the majority of folds, the optimal values for a and b were 1.4 and  −0.8 respectively. The MDTF is a first indication for the density. So in non-dense breasts the estimates as obtained with F1 are used while for dense breasts a combination of F2 and F3 works best.

Figure 7 shows results of the two breast volume estimation methods in comparison to MRI. Breast volume based on mammography has a linear relationship with the breast volume calculated from the MR images with both approaches. A Pearson's correlation coefficient of 0.96 (0.95–0.97 95%CI) was obtained for the semi-circle model and for the model in which the breast region that is in contact with the compression plates has a direction dependency. We decided to use the second approach, which is based on the idea of de Groot et al, for further calculations. This method appears to have less bias as the data points are closer to the identity line.

Zoom In Zoom Out Reset image size

Figure 8 shows the comparison between volumetric percent density estimates from DM and MRI. The Pearson's correlation coefficients between the volumetric percent density estimations per study were 0.86, 0.89, 0.74, 0.87 and 0.90 when using F1, F2, F3, the combination of F1 and F2 and the combination of F1  −  F3, respectively. When using the fibroglandular tissue estimations with F1 and the semi-circle model for the breast volume a correlation of 0.88 was found. Pearson's correlation coefficients with 95% CI are shown in table 2 for comparisons per breast and per exam.

Zoom In Zoom Out Reset image size

We determined MDTF values for different projection angles in the MR images of 20 dense breasts and the 202 segmented MR volumes. With a rotation interval of 5 degrees, this yielded to MDTFs of 36 different angles. The mean MDTF ranged between 0.11 and 0.95 and the coefficient of variation of the measurements ranged between 0.01 and 0.54. In figure 9(a) a scatter plot of the mean MDTF and the coefficient of variation is shown. We see that the coefficient of variation decreases with an increase of the mean MDTF. The low coefficient of variation of dense breasts (high mean MDTF), indicates that the MDTF is almost direction invariant in extremely dense breasts. On the other hand, the coefficient of variation is larger of small mean $MDT{{F}^{\prime}}s$. As the MDTF is a first indication for the breast density, we can see that the idea of direction invariance is violated in non-dense breasts. Hence, it is no surprise that the MDTF used to estimate F3 becomes unreliable, resulting in unrealistic low dense volume measurements (see also figure 6(c)).

Zoom In Zoom Out Reset image size

The comparison of the maximum dense tissue fractions measured in coronal planes of the compressed breast MRI's is shown in figure 9(b). Results obtained in the mediolateral direction, the direction of the compression, are on the vertical axis, while results of the craniocaudal direction are plotted horizontally. It can be observed that the maximum dense tissue fraction (MDTF) in the two directions are quite similar, which supports the validity of our approach.