A computational model to generate simulated three-dimensional breast masses

Purpose: To develop algorithms for creating realistic three-dimensional (3D) simulated breast masses and embedding them within actual clinical mammograms. The proposed techniques yield high-resolution simulated breast masses having randomized shapes, with user-defined mass type, size, location, and shape characteristics. Methods: The authors describe a method of producing 3D digital simulations of breast masses and a technique for embedding these simulated masses within actual digitized mammograms. Simulated 3D breast masses were generated by using a modified stochastic Gaussian random sphere model to generate a central tumor mass, and an iterative fractal branching algorithm to add complex spicule structures. The simulated masses were embedded within actual digitized mammograms. The authors evaluated the realism of the resulting hybrid phantoms by generating corresponding left- and right-breast image pairs, consisting of one breast image containing a real mass, and the opposite breast image of the same patient containing a similar simulated mass. The authors then used computer-aided diagnosis (CAD) methods and expert radiologist readers to determine whether significant di ff erences can be observed between the real and hybrid images. Results: The authors found no statistically significant di ff erence between the CAD features obtained from the real and simulated images of masses with either spiculated or nonspiculated margins. Likewise, the authors found that expert human readers performed very poorly in discriminating their hybrid images from real mammograms. Conclusions: The authors’ proposed method permits the realistic simulation of 3D breast masses having user-defined characteristics, enabling the creation of a large set of hybrid breast images containing a well-characterized mass, embedded within real breast background. The computational nature of the model makes it suitable for detectability studies, evaluation of computer aided diagnosis algorithms, and teaching purposes. C 2015 American Association of Physicists in Medicine. [http: // dx.doi.org / 10.1118 / 1.4905232]


INTRODUCTION
Computational phantom models are widely used in the evaluation of imaging systems and algorithms, where it is often important to have access to large and diverse sets of image data having known characteristics. In this work, we present a new method for generating realistic threedimensional (3D) mass phantoms, expanding on our prior work. 1 We describe an approach for creating simulated masses and a method for embedding the simulated masses digitally within real clinical mammograms to obtain digital hybrid images. 39 We then show that the resulting hybrid images are, for practical purposes, essentially indistinguishable from real mammograms as judged by both expert readers and computer-aided diagnosis (CAD) procedures.
The model introduced in this work employs a Gaussian random sphere (GRS) technique to generate a central mass and an iterative branching algorithm to simulate spicules. The branching process is informed by the principle of minimum work. The random generation of simulated masses is controlled by user-specified parameters. This permits the user to produce an unlimited family mass having particular characteristics.
Several previous approaches to breast tumor simulation found in the literature (e.g., Refs. [2][3][4] are based on simple shapes such as ellipsoids and cylinders. Typically, realism of a simulation method is quantified by the ability of readers to distinguish real tumors from simulated ones, with area under the receiver operating characteristic (ROC) curve (AUC) of 0.5 (pure guessing) being the ultimate goal. In a previous published study by Saunders et al.,2 a proposed mass model yielded AUC values of 0.68 ± 0.07 for the benign masses and 0.65 ± 0.07 for the malignant masses. Thus, the results differed significantly from the target AUC value of 0.5. A similar reader study by Berks et al. 3 was performed for a different model, resulting in a similarly unsatisfactory AUC result of 0.7 ± 0.09. Bliznakova et al. 4 proposed a model based on various geometrical shapes, but readers were easily able to distinguish real masses from the simulated ones, with accuracy exceeding 95%. Rashidnasab et al. proposed several models for mass simulation in mammography using a diffusion limited aggregation algorithm 5 and a random walk algorithm. 6 Their simulated images were generated using the same approach as previously described by us, 1,7 in which healthy mammography images are modified by substituting the effect of a mixture of fatty and fibrous tissue with tumor tissue. The range of structures created by their method was able to be controlled by a set of parameters and produced realistic looking results. However, their work did not provide any distinction in the generation of different mass types (as spiculated or nonspiculated). Further models that have been proposed in the literature [8][9][10] are either simplistic or lack thorough validation. The purpose of this paper is to create breast mass phantoms that are more detailed and realistic than those produced by existing techniques, as well as producing different mass types that can be differentiated by a spiculated or nonspiculated margin.
Our work was originally motivated by a specific need for highly detailed simulations for our research on phase contrast imaging 11,12 where fine details are clearly seen in the image. Note that the purpose of this paper is not to faithfully replicate the biological processes of tumor growth, but simply to produce simulated masses that appear visually similar to what is seen in medical imagery, so as to provide imaging researchers with a tool for evaluating imaging systems and algorithms. The proposed method may also be useful as a training tool.
In Sec. 2, we discuss modeling of the simulated breast mass, computation of its voxelized volume and projection, and its geometrical interpretation. We also describe a method for creation of hybrid digitized mammograms consisting of a simulated mass that is computationally embedded within real breast background. Experiments to measure the realism of the results, using CAD methods and expert readers, are provided in Sec. 3.A, and discussion follows in Sec. 4.

MATERIALS AND METHODS
In the proposed approach, simulated breast masses are modeled according to the steps shown in Fig. 1. Figures  1(a)-1(c) show the process of simulating the central mass, beginning from a shape constructed using GRS, 13 and modifying its surface with various irregularities, first by set of F. 1. Process of generating a simulated mass model and embedding it into an existing mammogram. low-frequency modifications referred to colloquially as "bumps" or "spikes," and then by adding high-frequency modifications to the surface defined as a "fuzzy" surface texture. Next, if a spiculated mass is desired, spiculation structure can be added to the central mass by using an iterative branching algorithm, as shown in Fig. 1(d). Finally, as shown in Figs. 1(e) and 1(f), a projection image of the simulated mass can be embedded within a clinical mammogram to produce a hybrid image (mammogram with simulated mass). The mass is described as a 3D shape model, so it can also be used in simulations of tomographic imaging; however, we did not consider the performance of our method in that setting. The simulated mass is generated in the form of a parametric surface model, so it can be rasterized to form a volumetric image represented on a 3D voxel grid; however, an analytic projection approach has tremendous computational benefits, therefore we describe the analytic solution for planar mammography in the Appendix. Sections 2.A-2.D provide the details of each of the steps of mass simulation shown in Fig. 1.

2.A.1. GRS model [Fig. 1(a)]
In this section, we explain the first step shown in Fig. 1(a) in which the central mass is simulated. We accomplish this by employing a GRS model, which was originally designed to model planets and comets 13 but has also been used in other fields. 14 The GRS is a parametric surface model described by the radial distance of the surface from the origin, r (θ,ϕ), which is given by the following function of spherical coordinates θ and ϕ (see Fig. 2): in which the logarithmic radius s(θ,ϕ) (a random function) is a series of spherical harmonics Y l m (θ,ϕ) truncated to a F. 2. GRS surface with the introduction of two low-frequency modifications: A Gaussian bump (left side) and a spike irregularity (right side), exaggerated here for illustration purposes. maximum order l max , defined as In Eq. (1), α is the mean of r (θ,ϕ), defining the size of the mass, and σ 2 is the variance of r (θ,ϕ), defining the degree of irregularity of the surface of the mass. For s(θ,ϕ) to be normally distributed with zero mean and an angular covariance function appropriate for a closed surface, 11 the expansion coefficients s l m in Eq. (2) are defined by where x G and y G are Gaussian random variables with zero mean and unit variance, i = √ −1 , and δ m0 is the Kronecker delta function. The exponent v is the power-law index of the covariance function. In this work, we fixed v = 4 to correspond to nonfractal shapes without sharp features in their geometry. 14 The coefficients for the negative values of m are defined as follows: where the asterisk denotes complex conjugate and Im(·) denotes the imaginary part. The statistics of the GRS shapes are controlled by l max and σ 2 . For example, weighting the spectrum toward higherdegree harmonics results in Gaussian spheres with finer surface irregularities. Increasing the variance of the logarithmic radius enhances these irregularities radially. Setting the variance to zero, for example, produces a sphere because all surface locations have the same radial distance in this instance.

2.A.2. Introducing low-frequency modifications [spikes and bumps; Fig. 1(c)]
In the next step, the GRS model is given an irregular surface at discrete locations by introducing low-frequency modifications we will refer to colloquially as spikes and bumps [see Figs. 1(b) and 2]. Spikes introduce pointy localized surface changes into the GRS model while bumps introduce localized, lobulated surface changes. These modifications were included to allow the central mass to have a greater degree of surface variation and hence greater realism, especially for nonspiculated masses for which these are the only fine surface structures.
Spikes are introduced as follows. For each modification j of the mass surface, a spherical coordinate pair θ c j ,ϕ c j is randomly chosen, defining the axis of revolution of the modification. The initial GRS model r (θ,ϕ) is then modified in the vicinity of each chosen coordinate pair, according to a specified function as follows. For a given coordinate pair θ c j ,ϕ c j , values of the random spike radius r spike and length l spike are selected and the initial GRS surface is modified using the quadratic ramp function , for d r (θ,ϕ),r θ c j ,ϕ c j ≤ r spike r (θ,ϕ), for d r (θ,ϕ),r θ c j ,ϕ c j > r spike (5) in which r ′ (θ,ϕ) is the new radial position of the mass surface and d (·,·) is the Euclidean distance between two points on the surface. The symbol ± indicates that these variations can be defined to grow outward (+) or inward (−) from the initial GRS surface.
Bumps are introduced in the following manner. For a given coordinate pair (θ ci ,ϕ ci ), values of the random radius r Gauss and length l Gauss are selected. Radial positions r (θ,ϕ) are altered to produce bump features according to in a similar manner as for the spike shapes. Fig. 1

(b)]
A fuzzy surface is created by modifying the surface profile function as follows: where r ′′ is the new radial position of the mass surface, n is a standard normal random variable, and α is a user-defined parameter controlling the variance of the surface variations.

2.B. Modeling spicule structures
The process of mass generation may conclude at this point; but if a spiculated mass is desired, then the algorithm proceeds to the step illustrated in Fig. 1(d), in which spicules are introduced into the central mass structure. This is accomplished by an iterative fractal branching algorithm that recursively creates a set of segments b n , n = 1,..., N based on a set of growing rules. Each of these segments is characterized by a starting location p sn in 3D space, an ending location p en , an initial radius r ini n , and a final radius r fin n , which defines a conical frustum with a hemisphere at the end acting as a "cap," as represented in Fig. 4. Additional variation to the defined geometries is introduced by adding a normal random variable to the distance to the center of each geometry (the revolution axis in the case of the frustum and the center location for the sphere), in a similar way as the high-frequency variations were introduced in the central mass [Eq. (7)]. The growth process is defined by user-selected parameters, including the distribution, bifurcation probability, direction of extension, emerging density, radius, and length of the introduced frustum shapes. Next, we describe the iterative algorithm formulated for the segment set generation.

2.B.1. Iterative branching structure generation
A temporary segment structure at a particular iteration k is defined by the set S k =  s k 1 ,..., s k j ,..., s k L k  , where s k j is the jth segment in the kth iteration of the algorithm, defined by a starting spatial location q k j , a direction of growth defined by angles , a length l k j , and an initial radius r k j . The initial segment set S 0 is generated as described in the flowchart shown in Fig. 3. The user specifies the number of initial segments L 0 and a number M group that indicates the number of initial segments clustered within each neighborhood. The latter is used to allow spicules to emerge in bunches. The starting position coordinates q 0 i are grouped in neighborhoods centered at r ′′ θ ini m ,ϕ ini m , where the coordinates are randomly selected within each neighborhood m. The starting position of segment n is added to neighborhood m as r ′′ θ ini m + log 2 n + γ 1 π/50, ϕ ini m log 2 n + γ 2 π/50 , where γ 1 and γ 2 are uniform random variables. That is, one segment emerges from the center of the neighborhood and the rest are added at increasing angular distances from this neighborhood center, each segment having a random length l 0 i and initial radius r 0 i which are drawn from a Gaussian distribution with user-defined properties. The direction of growth of these segments is defined in spherical coordinates by and where δ 0 and δ 1 are normal random variables with zero mean and a user-defined standard deviation (std) that controls the possible variance expected in the direction of the emerging segments from the radial direction with respect to the center of the mass (in our work, the standard deviation of δ 0 and δ 1 was set to 15(180/π) rad). Figure 4 illustrates an example of the first iteration (iteration 0), showing an initial segment (with parameters having superscript 0) and two child segments (denoted by superscript 1). In show that increasing the number of initial segments L 0 produces a larger number of spiculations, and that increasing the number of maximum emerging segments in a neighborhood M group causes the spicules to emerge from isolated regions in the central mass, while decreasing this parameter results in a more uniform distribution of spicules.
A fractal branching algorithm is applied iteratively to generate the complete set of segments forming a spicule structure. At each iteration k, given the segments S k generated in the previous iteration (parent segments), a new set of segments S k+1 (child segments) is generated following the flowchart shown in Fig. 6. The child segments act as parent segments in the next iteration of the process, which continues until no new segments are generated per a given set of growing rules.
Each parent segment produces zero, one, or two children, with user-defined probabilities ( Fig. 4 shows a two-child generation event). The generated child segments have a starting point that is the ending point of the parent, a radius less than or equal to the final radius of the parent segment, and a direction of growth that deviates randomly from that of the parent. Child segment diameters are scaled such that the flows within a branch follow the physiological principle of minimum work. 9,15 That is, the radii of the parent and child branches can be related by where r k+1 j and r k+1 j+1 are the resulting child segment radii, r k i is the parent segment radius, d r , the dividing ratio, has a value between 0.5 and 1, and τ is a constant known as the diameter exponent 16 (τ = 2.6 is used here per previous work 9 ). Child segments have equal radii when d r = 0.5 and the radius difference increases as d r increases. At the kth iteration, no new segments are added to the structure if a parent segment has a length or initial radius F. 4. Spatial interpretation of an initial parent segment (superscript 0) and child segments (superscript 1) generated in the branching algorithm. smaller than the maximum defined resolution (in this work, a tenth of the pixel size of the imaging system where the breast mass phantom will be used). Otherwise, one of three possible different scenarios in the branching scheme is followed by each parent segment with defined probabilities: 2.B.1.a. Continuing branch. One child segment s k+1 j is generated as a continuing segment from the parent s k i , with slightly reduced length and radius, and a slight change in direction. The radius of the child segment is computed, where decrease r is a chosen factor that models how fast the radius of each branch decreases along its length in the complete structure, the value of which is discussed later. The length of the continuing child segment is computed in a similar way where decrease l is a chosen factor that models how fast the segments' lengths decrease along a branch. The orientation of the continuing segment has a slight random variation from the parents in terms of azimuth and inclination θ k+1 where x 1 and x 2 are independent normal random variables and γ bif is an angle bifurcation factor. The factor γ bif , defined a priori, describes the angular variation in the direction of child segments with respect to their parents, and has a value between 0 and π (smaller values avoid potentially abrupt random changes of the growth direction).
F. 6. Flowchart for the generation of the iterative branching structure.
The factor −|x 2 | ϕ k i /4 appears in Eq. (15) to give the branches a tendency to be stretched along the X-Y plane, which simulates the effect of breast compression. Thus, x 2 is either defined as an independent normal random variable for the compressed-breast case, or set to zero for the uncompressed-breast case. This stretching effect is illustrated in Fig. 7.
After the continuing child is generated, the parent segment s k i is added to the complete structure to form a new element b n = p s n ,p e n ,r ini n ,r fin n of the form 2.B.1.b. Symmetric bifurcation. Two continuing child segments s k+1 j and s k+1 j+1 are generated from the parent, each having similar length and radius, and a similar, but opposed change, in direction. In this case, the dividing ratio d r shown in Eqs. (10) and (11) is computed as a uniform random variable taking values from 0.5 to 0.8. The radii of the two child segments are computed, In the same way, the two child segment lengths are defined, The two different orientations for the child segments are generated with opposed azimuth Gaussian variation in which x i , i = 1,..., 4 are independent, normal random variables. The variables x 5 and x 6 have the same horizontal stretching interpretations as in the previous scenario [parameter x 2 in Eq. (15)], and assume values different from zero when we consider breast compression. Once the two bifurcating children are generated, the parent segment s k i is added to the complete structure forming a new nth element b n = p s n ,p e n ,r ini n ,r fin n . The starting point, ending point, and initial radius are computed in the same way as in the previous scenario [Eqs. (16)-(18)]. The final radius is computed using the dividing ratio for this particular scenario and s k+1 j+1 are generated: the first forms a continuation of the parent branch; the second (smaller) branch diverges from the parent. The radius and length of the child segments are computed in a similar way as in the symmetric bifurcation scenario, but with a dividing ratio d r computed as a uniform random variable assuming values from 0.8 to 1. The continuing segment also has lower angle deviation from the parent segment than the bifurcation. Their orientations are described in a manner similar to the previous scenario, but with azimuth variations given by After the two bifurcating children are generated, the parent segment s k i is added to the complete structure forming a new nth element b n , defined in a manner similar to the previous scenario.
In Figs. 5(c)-5(e), we can observe the effect of adjusting the parameters of the spiculation model, as indicated in the inset table. Comparing Figs. 5(c) and 5(d) shows that larger values for the decrease r and decrease l parameters produce moreextended spicule structures. Comparing Figs. 5(d) and 5(e) shows that defining a higher continuing branch probability reduces the chances of bifurcation, allowing the spicules to grow more extensively.

2.C. Embedding in real tissue images
To allow us to validate the use of our proposed mass simulation model in the context of mammography, we developed a method of embedding the masses within actual clinical digitized mammograms using published values of breast tissue attenuation 17 and a simulated mammography source spectrum. 18,19 Our goal was to computationally modify the mammogram to include the simulated tumor projection by substituting the attenuation effect produced by a mixture of fat and fibrous tissue with that produced by tumor tissue. In these experiments, we used digitized film mammograms from the digital database for screening mammography (DDSM) database. 20,21 We began the embedding process by converting the image values to film optical density values by using the scanner conversion function (as indicated in Refs. 20 and 21). We then computed normalized intensity values from the film optical density values. Assuming a linear relationship between inverse of the recorded optical density and the logarithm of the intensity, we inverted the optical density values and later inverted the logarithm relationship to compute the normalized intensity image (normalized by the integral of the source intensity, that is, values of pixels in air should be 1). According to Beer's law, the normalized intensity image follows: where I d (x, y;ε) is the intensity recorded at the detector at x-ray energy ε, S(ε) is the intensity of the source, ε max is the maximum photon energy produced by the source, L is the thickness of the sample in centimeter, and µ obj (x, y;ε) is the attenuation coefficient of the sample in cm −1 , which is energy-dependent. Since the attenuation coefficients have a nonlinear dependence with energy (Ref. 17), the energy dependence of the integrand I d (x, y;ε) in Eq. (31) is not only on the source intensity but also on the breast thickness and composition, which are not known a priori. Therefore, in order to modify the equation to substitute the attenuation effect produced by healthy tissue with that of the simulated mass, we first need to estimate the composition and thickness of the sample. We approximated the sample composition in the form of a ratio of fat to fibrous tissue (we assume that no lesion is located within the embedding area). Once these factors are estimated, the integrand I d (x, y;ε) in Eq. (31) can be approximated aŝ I(x, y;ε) =Ŝ(ε) · e −(µ fat (ε)ratio fat (x, y)+µ fib (ε)ratio fib (x, y))L(x, y) , (32) whereŜ(ε) is a simulated mammography source spectrum; 18,19 µ fat (ε) and µ fib (ε) are the attenuation coefficients (in cm −1 ) of fat and fibrous tissue, respectively; and ratio fat (x, y), ratio fib (x, y), andL(x, y) are approximations of the fat tissue ratio, fibrous tissue ratio, and tissue thickness observed in the mammogram, respectively, computation of which is described shortly. The simulated tumor projection is then embedded within the original mammogram by substituting the effect caused by the mixture of fat and fibrous tissue with tumor tissue where µ tum (ε) is the attenuation coefficient (in cm −1 ) of tumorous tissue, either benign or carcinoma tissue, and T tum (x, y) is the thickness of the projected simulated tumor at location (x, y) (in cm). The resulting hybrid intensity image I T (x,y) was converted back to optical density values similar to the ones presented in the original digitized mammograms by considering the linear relationship between the logintensity and optical density values and inverting the particular normalization of the scanner used to digitize the original film mammogram. 20,21 The breast thickness considered in Eq. (33) was approximated by considering a secondary image I F (x, y), which approximately describes the intensity recorded by a breast of similar thickness characteristics to the one presented, but composed entirely of fat tissue. That is, in theory We generated the image I F (x,y) by fitting a thin plate spline (TPS) 22 to a selection of candidate locations in I(x,y) where fat tissue was the main breast component. These candidate locations were chosen by selecting the pixels in I(x, y) that are monotonically decreasing from nipple to chest wall. The rationale behind this is that, assuming that breast tissue thickness is not expected to decrease from nipple to chest wall, higher intensity values in this direction should be observed where fat tissue percentage is the most predominant, since fat is the least-absorbing tissue type in breast. We can observe an example of the resulting fat-tissue approximation image in Fig. 8(b), where the logarithm of the intensity I F (x, y) is shown (showing the image in terms of absorption), for easier comparison to the original mammogram shown in Fig. 8(a). Following Eq. (34), we can find the approximated thickness map of the breast sample using a constrained least-squares minimization process 23 Once the approximated sample thicknessL(x, y) is computed, we can obtain the approximate fibrous and fat ratios in Eq. (32) using a similar constrained least-squares process ratio fib (x, y) = arg min ratio fib (x, y) ratio fat (x, y) = 1 − ratio fib (x,y).
A simulated mass can be embedded at any location within any case so long as it can be accommodated within the thickness and spatial extent of the breast at that location. Figure 8(c) shows an example of the visually realistic results obtained by the proposed embedding procedure. Figure 9 shows mammograms for two patients each having one breast containing a real mass and one normal breast. We embedded a simulated mass into the image of the normal breast having similar visual characteristics (general shape, degree of spiculation, size, and location) to that of the actual mass in the opposite breast. In Fig. 9, the normal breast (left), simulated mass in that breast (center), and actual mass in the other breast (right) are shown.

2.D. System performance characteristics and variable values
The modeling software was written in  and executed using a Windows desktop computer (32-bit version Microsoft Windows Vista; 2 GHz dual-core processor with  4 GB RAM). The most challenging task given the limitations of processing speed and memory was the computation of a simulated mass projection and volume. Generating the volume and projection of the simulated tumors directly from a rasterized version proved to be computationally intensive; therefore, we developed a more-efficient analytical solution (described in the Appendix). The computational time and memory requirements for the generation of the simulated T I. Computation times and memory requirements for algorithm steps.  Table I. In our experiments, the simulated masses were defined within a 1000 × 1000 × 1000 voxelized volume in which each voxel had a size of 50 µm in each direction. The tumor simulation described by this model is randomized, following a series of user-defined variables describing the angular undulation of central mass, shape, amount, and size of the included modifications, number and density of emerging spiculation structures, and their growth properties. The parameters were chosen so as to match the appearance of breast tumors found in clinical radiology, and were systematically changed in order to produce a lesion that visually looked to be a mass with distinct features (e.g., lobulated margin versus circumscribed margin). This was accomplished using knowledge from published work, 9,24-28 and guidance from two of the authors (R.A.S. and R.M.N.), both of whom are experts in breast radiology. The process of adjusting the parameters to produce visually realistic simulated masses was conducted independently from creation of the set of hybrid mammo-grams in the evaluation, so as to avoid bias. The resulting parameters of the cases included in the evaluation of our method are summarized in Table II. The radius and length factors included in the table indicate a factor of the mean central mass radial distance value α, selected for the tumor generation in the GRS model [Eq. (1)]. cases, we selected a simulated mass with similar size and spiculation level to the mass that was actually present in the abnormal breast and embedded it in a similar breast location (distribution of parameters shown in Table II). The simulated masses were generated independently from the case and location within the breast in which they were later embedded. We embedded the simulated masses so as to resemble the real breast masses in terms of size, degree of spiculation, and location. (In concept, a simulated mass generated by the proposed approach can be embedded at any location in any case, provided that it fits within the spatial extent and thickness of the breast at that location.) The parameter value distribution was determined, and a set of simulations was generated, in a process separate from that used in the evaluation process, in consultation with two experts in mammography (authors R.A.S. and R.M.N.). We then digitally embedded each simulated mass within the MLO or CC view (chosen randomly for each pair) of the healthy breast using the proposed embedding scheme. By this approach, we created 83 corresponding left-and right-breast image pairs (a total of 166 images), in which the image of one breast depicted an actual mass, while the opposite breast image contained a similar simulated mass, both in the same view (either CC or MLO). Of the 83 cases on which our experiments were based, 31 exhibited benign masses that were classified as nonspiculated in the DDSM database; the remaining 52 cases exhibited malignant tumors that were classified as spiculated. Figure 11 Fig. 12, where 3D representations of the simulated masses are displayed alongside corresponding regions of interest of the simulations embedded in real breast tissue and their corresponding opposite breast containing a real mass.

2.E.1. Validation of hybrid phantom images for CAD analysis
To validate the proposed hybrid phantom images in the context of CAD, we applied a well-established set of quantitative CAD features [29][30][31] that are widely used to charac-terize lesions in digital mammography, the purpose being to observe whether our simulated masses yield the same results as their matched real masses when analyzed using CAD techniques. We only briefly describe the CAD procedures, as the details of these methods can be found in prior publications.
Prior to feature extraction, we segmented each mass in a hybrid image using the region-growing algorithm proposed in Ref. 29. Figure 13 shows an example segmentation of one of our hybrid images, indicating the grown region of the tumor, the tumor margin, an encompassing region, and the surrounding periphery. 30 The surrounding periphery is obtained by a morphological opening applied to the grown region. 30 The encompassing region and surrounding periphery simply extend the region window by 20 pixels along the left, right, top, and bottom edges of the mass.
Five CAD features are next extracted from the detected regions by the methods described in Refs. 29-31. Two of these features-sharpness and full width at half maximum (FWHM) angular deviation (related to spiculation level)-characterize the tumor margin and are important to distinguishing malignant tumors from benign ones; the other three featuresaverage gray level, contrast, and texture-characterize the density of the mass. For more details about the extracted CAD features, we refer the reader to Refs. 29-31.

2.E.2. Reader study design
Our study was based on 83 pairs of corresponding left-and right-breast images, in which each pair consisted of a clinical digitized mammogram of an abnormal breast with a real mass, and a hybrid image (real mammogram with simulated mass) of the opposite breast in the same patient. Thus, the total data set consisted of 166 mammograms. All 166 images were rated independently and sequentially by five expert radiologists. Images containing real or simulated tumors were shown in random order. Each reader assigned a score for each image expressing his or her confidence that the mass shown was real or simulated. Scoring was done on the following seven-point scale: definitely real (0), probably real (1), possibly real (2), unsure (3), possibly simulated (4), probably simulated (5), or definitely simulated (6).
Customized software was developed to conduct the reader study. The software displays the images sequentially in random order to avoid ordering effects, and indicates the position of the tumor to be evaluated in each image. The observers were able to control the pan, zoom, white balance, T III. CAD feature values [mean, (std)] for the mammograms in each category (nonspiculated, spiculated) and p-values from a t-test of the difference in feature values between real and hybrid mammograms. and contrast of each image individually. The readers were entirely free to choose the reading pace; however, after an image was scored, the next image was displayed and the readers were not permitted to revisit previous images. The images were displayed on a 5-megapixel mammographic monitor (Totoku ME551i) with 11-bit grayscale, calibrated by the vendor to the DICOM grayscale standard display function. The ratings were saved for each reader independently, and later processed in a multiple-reader multiple-case (MRMC) analysis. 32-38 Table III summarizes the mean and std values for each of the features extracted from the real and hybrid mammograms, for both nonspiculated and spiculated tumors. Table III also contains the resulting p-values from a t-test comparing feature values from real and hybrid mammograms for each of the five CAD features. Results for nonspiculated and spiculated tumors were computed independently. The result shown in Table II is that none of the comparisons showed a statistically significant difference between results from the real and hybrid images at the level p < 0.05, i.e., p exceeded 0.05 in every comparison, in many cases by a large margin.

3.A.2. Discrimination power of features
We further explored the generated hybrid images in the context of CAD by determining the extent to which the CAD methods behave similarly on the hybrid images and real mammograms. We accomplished this by applying univariate classifiers to discriminate nonspiculated tumors from spiculated ones, following the approach described in Ref. 30, using either real or hybrid images. Classifier performance was then evaluated using ROC curves, with AUC being used to summarize the performance. In our data set, all the real tumors observed in the mammograms that were classified as nonspiculated were found to be benign at histology, and those classified as spiculated were found to be malignant (histology data acquired from the DDSM database 20,21 ). So in this case, the task for discerning tumors according to their spiculation level was equivalent to discerning benign tumors from malignant ones. The resulting measured mean and 95% confidence interval (CI) of the AUCs are summarized in Table IV, together with the results from a similar analysis obtained from previously published work on real breast tumors, 30 demonstrating, in general, good agreement between the real and simulated tumors and the published data. The AUCs for our simulated tumors were all slightly lower than that measured for our real tumors, but within the margin of error.
We compared the performance in terms of AUC obtained from the extracted features using ANOVA analysis, employing the - software, [32][33][34][35][36][37][38] testing the hypothesis that each feature's performance in distinguishing between nonspiculated and spiculated masses is the same for the real masses as for the masses simulated using the proposed model. The mean differences, 95% CI, and p-values obtained by using ANOVA analysis to test this similarity are summarized in Table V. The performance differences were very small and all p-values observed were well above 0.05, showing that no significant differences in performance between the real and simulated masses for the considered extracted features were observed. T IV. AUC (mean ± 95% CI) of the computer-extracted features in distinguishing between benign and malignant tumors for the real tumor set, the simulated tumor set, and previously published data (Ref. 30 Figure 14 shows the distribution of realism scores assigned by each radiologist (on the seven-point scale), displaying their means and 95% confidence intervals. Table VI provides a statistical comparison of the distribution of these ratings for each radiologist by the Mann-Whitney U-test, testing against the null hypothesis that the scores are drawn from equivalent distributions, similar to comparisons in previous studies where breast tumor simulations were proposed. 2,3 The collected data were also used to construct ROC curves for the task of distinguishing real from simulated tumors for each radiologist individually and in a MRMC analysis using the - software available from the University of Chicago. [32][33][34][35][36][37][38] The resulting AUC values in the MRMC scenario analysis, as well as the independent AUC values for each radiologist, are reported in Table VI, showing that AUC = 0.5 was in the 95% CI for each reader and mass type (nonspiculated and spiculated).

DISCUSSION
We have presented a method to generate a collection of random simulated 3D breast masses, with a user-defined spiculation degree. We have also described a method for embedding the simulated masses within real digitized mammograms, yielding hybrid images for which the ground truth about the tumor location, extent, and characteristics is known. We evaluated the realism of the simulated masses by measuring the extent to which CAD features extracted from our proposed hybrid images are similar to those extracted from digitized mammograms containing real masses, both in terms of the features' values as well as their discriminating power (as judged by AUC). We also conducted an expert reader study to evaluate visual realism of the simulated hybrid images, in which we measured the readers' ability to distinguish real and hybrid images, and found that they were not able to achieve significantly better than random guessing.
The values reported in Table III show that CAD features extracted from the hybrid mammograms resulted in a distribution similar to that obtained from the real mammograms. We found no statistically significant difference for any of the features (at level p < 0.05) between the real and hybrid images. Indeed, many of the comparisons yielded very high pvalues, suggesting good correspondence between the features computed from the real and hybrid mammograms. We further evaluated the discrimination power of each feature to discern between nonspiculated and spiculated masses, both for the real and hybrid mammograms, yielding the results reported in Table IV. The features were found to perform similarly on the real and hybrid images, and yielded discrimination power similar to that reported for these specific features in an independent study (Ref. 30). Although differences can be observed in the mean performance in real and hybrid images, further investigation of these differences by ANOVA analysis yielded no evidence of statistically significant differences (Table V).
We also evaluated realism by measuring the ability of five experienced radiologist readers to distinguish the real masses from the simulated ones. The realism scores assigned by the readers (Fig. 14) showed similar distributions for both real and simulated masses, all yielding values around 2-3 on the 0-to-6 realism scale (i.e., the readers rated all masses, whether real or simulated, as "unsure" or "possibly simulated"). This fact, together with comments received from the observers about the difficulty of the discrimination task, suggests that the hybrid phantoms are indeed realistic.
F. 14. Rating distribution for each radiologist in the reader study. Error bars indicate 95% confidence intervals. T VI. U-test comparison of ratings between real and simulated tumors for each observer, testing for the hypothesis that real and simulated tumors have the same score distributions. AUC values are shown for each individual reader and for a MRMC analysis. We quantitatively analyzed the overall performance of the readers given the collected scores for the task of discerning real and simulated masses in a MRMC scenario, which resulted in a mean AUC of 0.544, with a 95% CI of (0.44, 0.65) for nonspiculated tumors; and a mean AUC of 0.588, with a 95% confidence interval of (0.50, 0.68) for spiculated tumors, as indicated in Table VI. The independent AUC values for each observer are also indicated in the same table. The results suggest that the readers were not able to perform significantly better than random guessing in distinguishing between the real and simulated tumors. These AUCs are lower than previously published AUCs (Refs. 2-4) and close to the ideal 0.5. The large overlap of confidence intervals may be also due to the relatively small number of readers and reader variability and statistical significance might have been established if a greater number of simulations or readers were available. There was also some reader variability, which raises the question of whether less-experienced readers may behave differently from more-experienced ones. Reader 2 presented the highest accuracy in distinguishing real from simulated tumors in the spiculated category, while showing low accuracy in the nonspiculated category. Surprisingly, Reader 5, who has substantially more experience than the others, presented the highest accuracy in the nonspiculated category, but moderate accuracy in the spiculated category. One of the previously published methods with highest accuracy was described in Ref. 5, with mean AUC of 0.55 (but without a distinction between masses according to their margin spiculation level). As their computational embedding approach was based on previous work described by our group, 1,7 and which is expanded in this paper, we hypothesize that the proposed embedding process may play an important role in the realistic appearance of the simulations.
The method for generating an array of random simulations presented here depends on a large set of parameters that were user-defined and describe the shape of the simulated masses. The goal of the project was to develop a method that is capable of simulating the wide variation in appearance of breast masses on mammograms. We were not trying to develop a model from first principles, so we took a more pragmatic approach. With our goal in mind, we varied the parameters to change the appearance of the mass to get the desired appearance.
The exact values of the parameters were less important than the appearance of the lesion and how the parameters can be changed to change the appearance of the lesion. The parameter distribution employed in the evaluation of our method was determined by an iterative process guided by authors R.A.S. and R.M.N. Several simulation sets were generated from an array of parameter values, and these were later adjusted after discussion of the appearance of the generated simulations prior to the embedding process, followed by the generation of a new set and repeating the process until the authors were satisfied with the realism of the results. As the parameters reported in Table II correspond to those in our evaluation study, we suggest that they can be used as a guideline to generate successful simulations of breast masses. However, the sensitivity of each of these parameters for the generation of visually realistic results is not presented here and remains a matter to be studied in future work.
It is also important to also note that our methodology of embedding simulated masses projections onto real digitized mammograms does not explicitly take into account film blurring or scatter; 8 however, the manner in which we used the real background appears to have preserved sufficient blur and noise effects so as to obtain very realistic phantoms, as demonstrated from the CAD and human reader studies. But, we acknowledge that further considerations of sharpness and noise effect may improve the results of the embedding process, which will be taken up in future work.
Owing to the availability and ease of use of the DDSM database, we based these first experiments on that source of images. However, we anticipate that our 3D modeling approach can be adapted to other data sets and imaging modalities, which we will study in future research. Further validation studies would be required to prove the effectiveness of our method in these new settings.
The benefit of the proposed 3D mass phantom is that it allows a large number of breast masses to be simulated with known characteristics, extent, and location. We anticipate that the proposed method could have useful application in training of radiologists and in the evaluation and optimization of imaging systems and algorithms.

CONCLUSIONS
In this work, we have presented a realistic, 3D computational breast mass simulation model exhibiting the fine structures and details observed in mammographic images and we have described a method of embedding the simulated masses within real clinical digitized mammograms. The model's versatility allows the creation of a large number of different tumor cases, both benign and malignant, with borders ranging from smooth to highly spiculated. We have demonstrated that computer-aided diagnosis (CA) algorithms yield very similar results when applied to our hybrid mammograms and real ones and we have found that expert readers did not perform significantly better than random guessing in discriminating our hybrid mammograms from real ones. This tool may be helpful in detectability studies for modern breast imaging techniques where a large database of tumors of different sizes and characteristics with knowledge of ground truth of the mass replicated is a requirement. The method may also be useful for training purposes. We plan to make mass model data available via a web site; in the meantime, interested parties can obtain data examples by contacting the corresponding author.

(A14)
In order to find the minor axis of the ellipse, we need to find the corresponding frustum radius defined at that middle point, which we call r rad . This radius is determined by the point in the center axis of the frustum which is the closest to p center in a similar way as described before, p rad = p s +u rad d t = (x rad , y rad ,z rad ), (A16) r rad = ((r fin · u rad ) +r ini (1 −u rad )). (A17) Considering that the base of the frustum is spherical, we can then find the value of the minor axis of the intersecting ellipse, b cyl b cyl =  ( r rad 2 − p center − p rad 2 ) .
Then, the area covered by the frustum in the projection plane can be determined where Int eb (p), Int et (p), and Int ec (p) indicate if the point (p) is inside (1) or outside (0) of the projected area of the base, top, and rest of the conical frustum in the projection plane, respectively, as shown in Fig. 15, with the boundaries indicated in yellow. The top and bottom intersecting points, x xcylt (p), y xcylt (p),z xcylt (p) and x xcylb (p),y xcylb (p),z xcylb (p) , respectively, are then found following the definition of the ellipse formed with the intersection of the frustum and plane: x xcylt (p), y xcylt (p),z xcylt (p) = x p , y p ,z center + a cyl x xcylb (p), y xcylb (p),z xcylb (p) = x p , y p ,z center − a cyl (A23)