Evaluation of the BreastSimulator software platform for breast tomography

The aim of this work was the evaluation of the software BreastSimulator, a breast x-ray imaging simulation software, as a tool for the creation of 3D uncompressed breast digital models and for the simulation and the optimization of computed tomography (CT) scanners dedicated to the breast. Eight 3D digital breast phantoms were created with glandular fractions in the range 10%–35%. The models are characterised by different sizes and modelled realistic anatomical features. X-ray CT projections were simulated for a dedicated cone-beam CT scanner and reconstructed with the FDK algorithm. X-ray projection images were simulated for 5 mono-energetic (27, 32, 35, 43 and 51 keV) and 3 poly-energetic x-ray spectra typically employed in current CT scanners dedicated to the breast (49, 60, or 80 kVp). Clinical CT images acquired from two different clinical breast CT scanners were used for comparison purposes. The quantitative evaluation included calculation of the power-law exponent, β, from simulated and real breast tomograms, based on the power spectrum fitted with a function of the spatial frequency, f, of the form S(f)  =  α/f  β. The breast models were validated by comparison against clinical breast CT and published data. We found that the calculated β coefficients were close to that of clinical CT data from a dedicated breast CT scanner and reported data in the literature. In evaluating the software package BreastSimulator to generate breast models suitable for use with breast CT imaging, we found that the breast phantoms produced with the software tool can reproduce the anatomical structure of real breasts, as evaluated by calculating the β exponent from the power spectral analysis of simulated images. As such, this research tool might contribute considerably to the further development, testing and optimisation of breast CT imaging techniques.

f, of the form S( f ) = α/f β . The breast models were validated by comparison against clinical breast CT and published data. We found that the calculated β coefficients were close to that of clinical CT data from a dedicated breast CT scanner and reported data in the literature. In evaluating the software package BreastSimulator to generate breast models suitable for use with breast CT imaging, we found that the breast phantoms produced with the software tool can reproduce the anatomical structure of real breasts, as evaluated by calculating the β exponent from the power spectral analysis of simulated images. As such, this research tool might contribute considerably to the further development, testing and optimisation of breast CT imaging techniques.
Keywords: breast CT, breast model, software simulation, anatomical structure (Some figures may appear in colour only in the online journal)

Introduction
Mammography, and in particular digital mammography (DM), is a fundamental imaging technique in breast cancer screening and diagnosis. DM returns a two-dimensional (2D) digital representation of a compressed three-dimensional (3D) object. Therefore, tissues belonging to different planes are all projected onto the same x-ray image plane, making it difficult to detect possible abnormalities. This condition may be even worse in dense breasts, characterized by a high fraction of glandular tissue. In recent years, digital breast tomosynthesis (DBT) has been introduced as a form of 3D imaging for screening and clinical diagnosis of the compressed breast in order to overcome this limitation (Sechopoulos 2013). On the other hand, cone-beam computed tomography (CT) scanners dedicated to the uncompressed breast (breast CT, BCT) are available both experimentally and commercially, characterized by the use of quasi-mono-energetic (McKinley et al 2005) or poly-energetic x-ray beams (Lindfors et al 2008, O'Connell et al 2010, Mettivier et al 2011, Kalender et al 2012, Sarno et al 2015. Parallel-beam synchrotron radiation mono-energetic BCT is also under investigation (Longo et al 2016. But before BCT can become a clinical procedure a number of issues should be optimized, such as the source and the detector design (Kalender et al 2012), the acquisition strategy (Lindfors et al 2008, Mettivier et al 2012, McKinley et al 2012, and the reconstruction methods. To perform such investigations, there is a strong need of large databases of clinical images. Alternatively, images may be simulated from computational 3D digital breast models. They can be classified as digital phantoms based on patient data or mathematical data. In BCT, simple mathematical breast phantoms, usually in the form of cylinder, half-ellipsoid or slabs of homogeneous material with a given glandular to adipose breast ratio, are widely used in simulations particularly for dosimetry and optim ization of acquisition geometry (Boone et al 2004, Lanconelli et al 2013. However, when it is necessary to investigate parameters such as the detectability of lesions, the performance of image processing algorithms or the reconstruction algorithms, the use of a homogeneous background is a limitation, since the anatomical structure is not reproduced.
Mathematical breast phantoms for BCT may be produced also with the BreastSimulator software tool (Bliznakova et al 2003), which is a software application dedicated for research in x-ray breast imaging (Bliznakova et al 2015). This research tool allows the creation of realistic 3D uncompressed breast models. The simulation of the breast compression adopted during DM and DBT, is also included (Zyganitidis et al 2007). With this software, it is possible to simulate mammographic, tomosynthesis and fully tomographic breast imaging geometries. The BreastSimulator tool was previously validated and evaluated as a reliable tool for the simulation of DM systems (Bliznakova et al 2010, Mettivier et al 2016b. The purpose of this study was to validate this software tool as an appropriate x-ray simulator for dedicated BCT imaging. This investigation is based on the quantification of the anatomical noise, evaluated by calculating the β exponent deduced from the power spectral analysis of the simulated CT images. The closeness of the power spectrum coefficient β (calculated from simulated CT images) to that calculated from clinical CT images of the uncompressed breast was the criterion for validating BreastSimulator for 3D breast imaging studies.

BreastSimulator
The main components of the BreastSimulator (Bliznakova et al 2012) are the module for the creation of breast composition models and the module for the simulation of x-ray imaging.
2.1.1. Breast model generation module. This module was used to generate breast models of a given size and composition. The simulated breast features include the breast shape, the duct system, the Cooper's ligaments, the pectoralis muscle, the 3D mammographic texture, the skin, the lymphatic and blood systems and breast abnormalities (masses and calcifications). Additionally, a 3D texture matrix was created to simulate breast structures not explicitly model led, such as nerves and blood vessels, as well as to increase the realism of the simulated ligaments. The user can increase the complexity of the breast model by including any such features and by increasing their number (e.g. simulating a large number of Cooper ligaments, or of lactiferous ducts) or size (e.g. by simulating short-sized or long-sized ducts). The aim was to simulate the complexity of the actual breast anatomy, as found for women with different glandular fraction, breast size and shape, and anatomical texture. An example of a simulated breast phantom with and without Cooper ligaments is shown in figures 1(a) and (b), respectively.
The duct system is simulated as a network of cylinders (figure 1(a)), marked as a fibrous tissue and probabilistically arranged in the breast in a tree-like arrangement. The duct model includes the major ducts and the lactiferous ducts. Cooper ligaments are simulated as thin ellipsoidal shells, originating at randomly sampled positions in the breast model ( figure 1(b)). Their linear attenuation coefficient (at the given x-ray photon energy) is equal to that of the ducts, while the compartments enclosed by them are assigned the attenuation coefficient of the breast adipose tissue. A mixture of adipose, fibrous and connective tissues as well as other non-glandular tissue types not explicitly modelled, simulates the mammographic texture. The algorithm for generating this texture is based on the use of a random walk, following the concept of the 'fractional Brownian motion model'. The pectoralis muscle (figure 1(a)) is approximated as a cone-shaped object. Breast abnormalities are modelled with round, ovoid, elongated or irregular shapes.
2.1.2. X-ray imaging module. The x-ray imaging module contains information for the acquisition geometry and allows for setting acquisition parameters like source-to-detector (SDD) and source-to-isocenter (SID) distances, number of projection images, gantry angles, beam energy and detector type. X-ray projection images were obtained by simulating the transport of mono-energetic x-ray photons in the breast model. Image formation was based on the Lambert-Beer's law for x-ray attenuation. Figure 2 shows the BCT acquisition geometry used in the present simulation study. X-ray tracing was based on Siddon's algorithm to calculate the exact radiological path through voxels. Poisson quantum noise was added to the noise-free projection images, using a Gaussian random number generator, with a variance equal to the number of photons incident on each detector pixel. To simulate poly-energetic beams, the images obtained at each spectral energy were subject to weighted sum based on the corresponding photon fluence of the x-ray spectrum adopted for the simulation.

Setup description
For the imaging setup, we simulated a cone-beam irradiation geometry, with SID and SDD set to 458 mm and 866 mm, respectively. The detector was modelled with a size of 700 × 700 pixels, while the pixel size was set to 0.333 mm, in the order of the effective pixel size adopted in projection images with clinical BCT scanners (Boone et al 2005). Photon scatter in the breast and detector response were not simulated. In particular, scatter in uncompressed breasts, though significant, introduces a low-frequency trend on the spatial variation of the background signal, at  frequencies below the minimum spatial frequency (0.05 mm −1 ) considered here in the analysis of the power spectrum (PS) for evaluation of the power-law exponent, β .
We note that considering scatter-free views of the simulated breast may introduce a decrease in the noise power at low frequencies in CT slices, so producing a decrease of the power-law magnitude, α. However, for the purpose of this work, only the power-law exponent, β, will be analysed. On the other hand, both α and β will be considered as free parameters in the linear fit of the log(PS) curve log[S( f )] = (log α) − β · f in the analysis of the radial PS from simulated or measured CT slices.
For each projection and photon energy, 10 7 photon histories were run. Simulated images were generated for all incident mono-energetic beams with energy in the range 15-80 keV, with 1 keV increments. In the case of poly-energetic beams, for each energy we run a code to produce three projection images from the initial x-ray spectra and for each gantry angle. These three projection images have pixel values representing distances in mm (d adipose , d gland , d skin ) as the travelled distance of the x-ray through (a) the adipose tissue, (b) the gland tissue and (c) the skin, respectively. A total of 70 200 projections were obtained for the 80 kVp polyenergetic spectrum, while for the 49 kVp and 60 kVp spectra the number were 36 720 and 48 600, respectively. Then, the corresponding image per gantry angle was obtained as follows: where w E , μ(E) adipose , μ(E) gland, and μ(E) skin are the weighting coefficient calculated, based on the incident x-ray spectra and the attenuation coefficients for the corresponding tissue and energy E, respectively.
CT scans with 360 angular views were simulated, with mono-energetic (27, 32, 35, 43 and 51 keV) as well as poly-energetic x-ray spectra (49, 60, or 80 kVp). These last spectra correspond to the spectra used for patient scans performed with the breast CT scanners at University of California Davis Medical Center (UCDMC) (Lindfors et al 2008, Gazi andBoone 2014), and the commercial scanner of Koning Ltd (O'Connell et al) also used at Radboud University Medical Center (RUMC). The spectra were calculated using the spectral model by Boone et al (1997). The beam qualities were: (49 kVp, 1.39 mm Al) (Sechopoulos et al 2010); (60 kVp, 4.26 mm Al) (80 kVp, 5.74 mm Al) (Boone et al 2005). CT scans for mono-energetic beams with photon energies of 32 keV, 43 keV and 51 keV were simulated, since these energies correspond to the mean energy of the poly-energetic spectra at 49, 60 and 80 kVp, respectively. On the other hand, mono-energetic beams at 27 keV and 35 keV have a photon energy close to those used in the experimental studies ongoing at the Elettra synchrotron radiation facility on phase-contrast BCT (Longo et al 2016, Delogu et al 2017. Image pixel dimensions were 0.333 mm × 0.333 mm. The simulation of 360 x-ray projections took 24 h for the mono-energetic beams on a dedicated personal computer (see below).

Breast phantoms
For the purpose of this work, we used BreastSimulator to create eight different 3D digital breast phantoms with realistic anatomical features (breast models from #1 to #8, with different sizes and compositions). This research tool runs on an Intel Core 2 Quad Processor Q8200 2.33 GHz, with 8 GB RAM and 64 bit Linux operating system. The attenuation coefficients of the different breast components were derived from the XCOM Program (NIST Database). Table 1. Main features of the eight breast models. The breast size is expressed as the product of three half axes.  Table 1 shows the characteristics of each breast phantom (background matrix size, voxel dimensions, glandular fraction and breast size) and the corresponding simulated tissue features. The size (in pixels) of the 3D breast matrix ranges from (560) 3 for breast #2, to (1000) 3 for breast models #4 to #7. The voxel size of these matrices ranges from 0.10 mm to 0.25 mm (in each dimension) and the volume from 1.4 × 10 5 mm 3 (#7) to 8.5 × 10 5 mm 3 (#3). As a result of simulating a certain breast anatomy, each breast model turns out to have a different glandular fraction by volume which varies from 10% to 35% (i.e. the glandular fraction is an output datum of the software). The breast models were generated free of lesions in order to avoid any bias in the results. Based on the breast model complexity, the CPU time for the realization of a digital breast model ranged from 15 to 30 min.

Clinical data
In order to make a comparison with measured clinical data, the same procedure for the evaluation of β was applied on clinical BCT images acquired by the team at UCDMC and at RUMC. The data acquired with the UCDMC scanner involved a set of 180 breasts. The reconstructed slices have different dimensions with a voxel size of 0.20 × 0.20 × 0.35 mm 2 . Figure 3 shows clinical CT views and 3D renderings of one of the samples acquired with the UCDMC scanner. For comparison, figure 3 shows a simulated breast model (#5) with the same glandular fraction. The data acquired with the RUMC scanner include CT scans of eight different real breasts. The number of projections is 300 over 360° and the tube voltage was 49 kVp. The reconstructed slices have different dimensions with a voxel size of 0.27 × 0.27 × 0.27 mm 3 .

Assessment of the β parameter
It has been previously reported (Burgess et al 2001) that the anatomical structure in mammograms may represent a major impediment to lesion detectability and that it shows an isotropic power law spectrum of the form S( f ) = α/( f β ). Here, α describes the magnitude of fluctuation in the signal power, f represents the radial frequency and β (the log spectrum slope) has a critical role in determining the size at which a lesion reaches the threshold for detection (Chen et al 2012). Its value (for natural scenes and for mammography) ranges between 1.5 and 3.5 (Burgess and Judy 2007) and its mean value has been reported to be approximately 3 for DM (Heine andVelthuizen 2002, Burgess et al 2001) and approximately 2 for BCT (Metheany et al 2008, Chen et al 2012. For the analysis of the β parameter, the projections of the simulated uncompressed breast were reconstructed with commercial software (COBRA, Exxim Computing Corporation, Pleasanton, CA, USA) implementing the Feldkamp-Davis-Kress algorithm, providing axial, coronal and sagittal views. In order to evaluate the impact of slice thickness, the dimension of the reconstructed voxel varied from 0.250 × 0.250 × 0.250 mm 3 (512 × 512 × 512 voxels) to 0.250 × 0.250 × 2.00 mm 3 (512 × 512 × 64 voxels). The power-law exponent β was derived from fitting the PS functional form S( f ) = α/( f β ) evaluated on ROIs selected randomly in reconstructed slices (Chen et al 2013). A linear fit was applied to the radial log PS. Assuming uniform tissue characteristics in the various regions of the breast volume, one thousand ROIs have been selected randomly inside a single coronal slice or in the whole breast. A rejection method was used to ensure that all ROIs were located within the breast anatomy on the image. Following the method of Metheany et al (2008), for each ROI, the mean pixel value was subtracted and then a Hanning window was applied. The ROI was chosen empirically while ensuring that it was large enough to allow for an accurate estimate of β but not too large that it could emphasize non-uniformities in the image. Then, the 2D PS was computed by means of the Fast Fourier Transform for each ROI and the mean 2D PS was determined by averaging the PS from the 1000 ROIs. In order to obtain a 1D PS, a radial profile was evaluated. Finally, the β coefficient was calculated as the negative slope of the fitting line returned by computing a linear fit of log(1D PS) versus log( f ). The optimal frequency window (from f 1 to f 2 ) was selected for each PS( f ), as the one which maximized the R 2 fit statistic, and the values of α and β were then recorded for each breast data set. We tested this method by evaluating the β coeffi-     Figure 5 shows sample CT slices and the corresponding 3D renderings of all eight breast models created with BreastSimulator for this work, for a simulated 80 kVp spectrum. In this figure, the different anatomical complexities explained in the various models are reflected in the different appearance of the simulated CT slices. In particular, the inclusion in model breasts of the Cooper's ligaments and of the various sized duct systems are clearly distinguishable features of the various breasts. The glandular tissue is well distinguished from the flat tissue as well as the skin tissue (when present).

Analysis of the power spectrum
After CT reconstruction of breast projections created with BreastSimulator, we evaluated the β coefficient from the resulting CT coronal slices (figures 6(a) and (b)); an example of PS plots and corresponding power-law fits is shown in figure 6(c), together with PS data evaluated from clinical BCT scans at RUMC and UCDMC. The three illustrative cases shown in figure 6(c) were randomly selected between the simulated phantoms; in the case of the real breast data from the UCDMC set we selected a case with a β parameter equal to the measured mean value. We observe that the BreastSimulator data show a different intercept, log(α), of the fitting curve: indeed, the simulated data show a lower power, at all frequencies, than patient data.
Measured β values for all the simulated and real breasts are reported in table 2. Table 2 reports the β coefficients calculated in the central coronal slice of each breast model described in table 1, reconstructed with four thicknesses of the tomographic slice and for the different beam qualities. For instance, the results for a slice thickness of 0.250 mm and 0.500 mm are shown in figure 7, for both mono-energetic (27,32,35,43 and 51 keV) (figures 7(a) and (b)) and poly-energetic (49, 50 and 80 kVp) spectra (figures 7(c) and (d)).
We observe that the simulated projection images from computer breast models show a different value of the β coefficient, but for a given breast model, there is limited variation in the β values when the energy of the x-ray beam is varied. When considering the range of values found in UCDMC data (β = 2.11 ± 0.55), the data provided by the BreastSimulator partially overlap this range, with breast models from #2 to #6 falling inside this range. With regards to the change in β deriving from a different choice of the slice thickness in the CT reconstruction, the comparison of figures 7(a)-(d) illustrates that there is a slight increase in β when the slice thickness increases. These considerations are also confirmed by the data in figure 8 where the β coefficients are reported for all digital phantoms as a function of the beam energy for monoenergetic (figure 8(a)) and poly-energetic (figure 8(b)) spectra.
Figures 9(a)-(d) show the β coefficients for each phantom for the 80 kVp poly-energetic spectrum as a function of the phantom glandularity. The continuous and dashed lines in the figure show the fitted value and range of ±1 std. dev. reported in Chen et al (2013), respectively. Figure 10 shows the β coefficients as a function of the slice thickness, calculated in the central coronal slice for breast models #1 to #8 for the 80 kVp poly-energetic setup. The area shaded in light grey represents the average value (dashed red line, β = 1.86) ±1 std. dev. of the parameter β measured on clinical CT data reported by UCDMC on 44 breast images (Chen et al 2012). It is observed that breast models #1 to #4 are totally outside this area and that models #5 to #8 are totally or partly within it. We also note that Chen et al (2013) at UCDMC evaluated the β coefficients from coronal CT slices on a cohort of 185 patients and they found a mean value β = 1.96 ± 0.46.   ure 11(a), the values of β show some fluctuation with the slice position from chest wall to the nipple; however, independent parabolic fits to the data point relative to the six different beam qualities reveal a common trend, with β decreasing markedly toward the nipple and slightly decreasing also toward the chest wall. The simulated data for 80 kVp and for the corresponding 51 keV mono-energetic spectrum are replotted in figure 11(b), where it is seen that the two datasets of β values share the same mean value and standard deviation, for the slices along the longitudinal position from the chest wall to the nipple. A direct comparison between the β values within the simulated breast model #5 and a real patient scan (randomly selected), acquired with the RUMC scanner is shown in figure 12. The x-ray spectra utilized for the simulation is the same used in the commercial scanner. Good agreement within the experimental uncertainty is obtained.

Discussion
The power spectrum analysis of simulated and clinical images suggests that the BreastSimulator tool generates breast models with a tomographic characteristic close to the real breast tomograms. Eight 3D uncompressed breast models characterized with different content, size and voxel resolution were created with the BreastSimulator software package. While the agreement between real and simulated β distributions is not taken as a paradigm for 'correctness' of the anatomical structure of the 3D digital breast models, it is a positive aspect of BreastSimulator that it is capable of reproducing the range of variability of this parameter found in patients' BCT scans. In this respect, we point out that the detailed distribution of the glandular tissue in the breast volume has no influence on the evaluation of the β parameter, since in the evaluation procedure (as reported in section 2.5 and shown in figures 6(a) and (b)) a large number of ROIs were selected randomly in the slice for computation of the average pixel value. Just for sake of example, the average β value plays the same role as the average glandular dose in breast dosimetry, where two different 3D distributions of the glandular dose might share the same MGD. It is the fraction of the glandular tissue that has an influence on the β value as reported in the literature (e.g. Metheany et al (2008) and Chen et al (2013)) and shown in figure 9. This is reflected in figure 3, where the spatial distributions of the glandular tissue in real and simulated breasts looks quite different, though the average β's are comparable.
The creation of fatty, glandular and dense breast models is influenced by the simulated breast features: ducts, Cooper ligaments, adipose, fibrous and other connective tissues. For each phantom, we simulated the tomographic acquisition with eight different setups  (5 mono-energetic and 3 poly-energetic setups) and obtained 64 CT reconstructions (figure 5) with four slice thicknesses (0.250, 0.500, 1.000 and 2.000 mm) for a total of 256 CT reconstructions. The β coefficient for the slices in each reconstruction was calculated and reported in table 2. The mean value and the standard deviation of calculated β values of simulated images (β ~ 2 for this 3D imaging modality) are in general in the range of the calculated β values obtained for the clinical CT images from dedicated breast CT scanners. Published BCT data from a relatively large cohort of patients indicate an average β exponent between 1.86 ± 0.38 and 1.96 ± 0.46 (Metheany et al 2008, Chen et al 2013. The increase of the thickness of the reconstructed slice results in an increase of the β values for each breast model and setup (table 1 and figure 8) as expected due to the increased complexity of the anatomical structure related to the greater thickness. Those results are in agreement with the data reported in Chen et al (2012) and in Chen et al (2013) where an increase of the β value from 1.96 to 3.15 was due to the increase of the slice thickness (from 0.23 mm to 44 mm).
The detailed comparison of β values calculated from simulated and real images (in figure 7) shows that, with the exception of breast models #1, 7 and 8, all other created breast phantoms have β values close to the β value calculated on the real images used in this study and are in the range of the β values reported in the literature. In particular, we can note that the β value depends on the breast matrix size (from about 600 3 to 1000 3 ) and the voxel side (from 0.25 to 0.1 mm) as well as on the number of Cooper ligaments. Lower voxel resolution results in averaging of adipose and glandular tissues in one voxel. Lower voxel resolution combined with higher glandularity leads to slices with highly irregular structures which result in higher β values, as in the case of breast #1.
Increasing the voxel side has the opposite effect on β. If the parameters breast size and voxel size are well chosen and set, the β values decrease with the decrease of the breast glandularity. These results are in agreement with the results reported also in Metheany et al (2008) and Chen et al (2013): in particular the data reported in figures 2 and 8 of these works, respectively. It was also observed that the β value depends on the number of simulated structures and anatomical details in the breast model. In particular, breast models #5 and #6 differ only in the number of Cooper's ligaments and their initial radius. The increase in the number of these ligaments leads to the increased β value. On the contrary, the β value is independent of the photon energy value (figure 10(a)) and tube kilovoltage ( figure 10(b)) in the mono-energetic and poly-energetic setup, respectively. This is expected, since the β value quantifies the anatomical structure in images.
Another aspect is highlighted by figure 11. In the literature, there is no indication about the appropriate choice of the slice used for evaluation of the β parameter; in our study we used the central coronal slice, following Boone (Chen et al 2013). Results for β values calculated in all the slices are reported in figure 11. A good agreement can be observed between these β values and those obtained for a real breast (figure 12) as well as data reported by Engstrom Figure 11. Power-law exponent β calculated on coronal slices (0.250 mm thick) of simulated CT scans for the breast model #5, as a function of the distance from the chest wall, for (a) all six beam qualities, and (b) for a mono-energetic (51 keV) and a poly-energetic (80 kVp) spectrum. The six quadratic fit lines (in red) in (a) show a decreasing trend of β at the sides of the chest wall and nipple. The order of the six fitting curves is indicated on the right side of the plot. In (b), the horizontal dashed line and continuous lines (in red) represent, respectively, the mean value and the (mean ± standard deviation) values of β calculated for the two datasets. The dashed vertical line (in blue) in the two plots marks the longitudinal position of the central slice. In this work, we used eight computer models of the breast with glandularity ranging between 10% and 35%. In the glandularity calculations, we included the skin. Four of the breasts had a glandular fraction of 10%, two had a glandular fraction of 24%, and two other breasts have 35% and 30% glandularity, respectively. The average glandularity was 23%, comparable to the value (19.3%) obtained by Yaffe et al (2009) for the average glandularity obtained from 191 patient breast-CT images. The current efforts of our group aim towards generating a very large number of computer models of the breast with different glandularity, sizes and dimensions, to be stored in a dedicated database and to be used for further research.
Regarding the computational time of about 24 h, it includes (a) generation of breast models, (b) generation of 360 projection images and (c) volume reconstruction. The most timeconsuming part of the algorithm for generation of breast models is the part related to the generation of the Cooper ligaments. The locations of the Cooper ligaments are sampled randomly within the breast volume; then, a procedure verifies that the generated ligament does not have any intersection with other ligaments. At present, this process has not been optimised computationally, but we are working for reducing the time needed for their generation. Simulation for CT imaging was performed on two 6-core processor workstations, with 24 GB RAM. One projection image was simulated in about 4 min. Image formation simulation may also be optimised if cloud technology is used. On the other hand, the time for reconstructing tomograms was comparably negligible.
We are also working on the use of the BreastSimulator to simulate noise in projection images. Noise sources to be considered are the photon noise, the noise in the detector as well as the scattered x-rays. In these studies, a Monte Carlo code will be exploited to simulate x-ray interactions in the breast models and in the detectors. It is expected that simulated tomographic images with improved noise description will have characteristics closer to those from clinical scans. This future work will further establish the practical importance of the BreastSimulator in carrying out feasibility studies in the field of breast x-ray imaging and in particular in advancing research and optimisation studies for breast CT.

Conclusions
This paper described the evaluation of a software package called BreastSimulator for research purposes in breast CT. Breast models of different size and content were simulated and the anatomical structure properties were evaluated by calculating the β exponent from the power spectral analysis of the simulated images. The good agreement between simulated and measured β in clinical scans indicated the potential of BreastSimulator in devising a digital phantom for describing the complex anatomy of the female breast. It is expected that the increase in the complexity of the present models (e.g. increase of the number of tissue structures) for breast CT with BreastSimulator will produce an even better description of the corre sponding complexity of the anatomical structure of the breast; this also depends on the computing power available for simulation. This would be particularly important for the description of dense breasts.