3-D Stochastic Numerical Breast Phantoms for Enabling Virtual Imaging Trials of Ultrasound Computed Tomography.

Ultrasound computed tomography (USCT) is an emerging imaging modality for breast imaging that can produce quantitative images that depict the acoustic properties of tissues. Computer-simulation studies, also known as virtual imaging trials, provide researchers with an economical and convenient route to systematically explore imaging system designs and image reconstruction methods. When simulating an imaging technology intended for clinical use, it is essential to employ realistic numerical phantoms that can facilitate the objective, or task-based, assessment of image quality (IQ). Moreover, when computing objective IQ measures, an ensemble of such phantoms should be employed, which displays the variability in anatomy and object properties that are representative of the to-be-imaged patient cohort. Such stochastic phantoms for clinically relevant applications of USCT are currently lacking. In this work, a methodology for producing realistic 3-D numerical breast phantoms for enabling clinically relevant computer-simulation studies of USCT breast imaging is presented. By extending and adapting an existing stochastic 3-D breast phantom for use with USCT, methods for creating ensembles of numerical acoustic breast phantoms are established. These breast phantoms will possess clinically relevant variations in breast size, composition, acoustic properties, tumor locations, and tissue textures. To demonstrate the use of the phantoms in virtual USCT studies, two brief case studies are presented, which addresses the development and assessment of image reconstruction procedures. Examples of breast phantoms produced by use of the proposed methods and a collection of 52 sets of simulated USCT measurement data have been made open source for use in image reconstruction development.

Introduction

ULTRASOUND computed tomography (USCT) is an imaging technique that utilizes tomographic principles to obtain quantitative estimates of acoustic properties, such as speed-of-sound (SOS), density, and acoustic attenuation (AA) [1]–[5]. Because it can produce high-resolution and high contrast images of tissue properties, the development of USCT as a breast imaging modality has received significant attention [5]–[10]. It has several advantages over other breast imaging modalities, such as mammography, including low cost and being radiation- and breast-compression-free [11], [12]. While commercial systems for breast USCT are being actively developed, USCT remains an emerging technology and a topic of active research [13]–[16].

When developing new breast USCT technologies, it is important to assess their clinical utility by use of objective measures of image quality (IQ). Given the large number of system parameters that can impact IQ and variability in the cohort of subjects to-be-imaged, a comprehensive assessment and refinement of modern imaging technologies, such as breast USCT via clinical trials, often is impossible. Furthermore, obvious ethical concerns preclude certain experimental designs that otherwise would be of great benefit toward optimizing imaging systems for diagnostic tasks, such as tumor detection and characterization. As a surrogate for clinical trials, computer-simulation studies of medical imaging technologies, also known as virtual imaging trials (VITs), have been advocated for assessing and optimizing system and algorithm designs [17]–[20]. VITs provide a convenient, safe and cost-effective way to explore system and algorithm designs in the early stages of technology development [21], [22].

For use in computing objective, or task-based, measures of IQ that serve as figures of merit (FOMs) for breast USCT designs, it is critical that VITs employ numerical breast phantoms (NBPs) that accurately convey the anatomical and acoustic properties of the female breast. Moreover, it is known that object variability (i.e., patient-to-patient differences in the breast anatomy and properties) can be viewed as a source of randomness present in image data that limit the performance of human or numerical observers on detection or estimation tasks [23]–[25]. It is therefore important to have the capability of producing ensembles of NBPs that possess prescribed statistical properties associated with a specified to-be-imaged subject cohort; these NBPs can each be virtually imaged and, subsequently, ensemble-averaged objective IQ measures can be computed for use in assessing and refining USCT imaging technologies. However, existing NBPs do not satisfy these requirements and are limited by factors that include oversimplified anatomical structures [2], [26]–[29] or are representative of healthy subjects only [30], [31]. NBPs derived from clinical magnetic resonance images are available [30] but are severely limited in number; as such, they do not accurately depict variability in breast anatomy or acoustic properties that will be present in a prescribed patient cohort. Other tools for generating NBPs [28], [29] rely on digital templates or segmented clinical images with simplified anatomical structures and consider only a limited number of tissue types. In summary, there remains an important need for developing NBPs for use in VITs of breast USCT that comprise realistic structures and acoustic properties, include lesions and/or other pathologies, and are representative of the stochastic variability in breast size, shape, composition, anatomy, and tissue properties observed in a specified cohort of to-be-imaged subjects.

Recently, the Virtual Imaging Clinical Trials for Regulatory Evaluation (VICTRE) project [17], [18] of the Food and Drug Administration (FDA) has validated and released software tools to generate realistic NBPs, as part of an end-to-end simulation framework for virtual mammography imaging studies. The breast size, shape, location, density, and extent of different tissues are tunable parameters, based on which stochastic and physically realistic 3-D numerical phantoms of tissue structures can be generated. The tool also allows to embed a variety of lesions (e.g., circumscribed or spiculated) at physiologically plausible locations.

In this work, a methodology for producing realistic 3-D numerical acoustic breast phantoms for enabling clinically relevant VITs of USCT breast imaging is presented. This will be accomplished by extending the VICTRE NBPs for use in USCT, which will permit virtually imaging of ensembles of NBPs whose physical and statistical properties are representative of clinical cohorts. Modifications to the VICTRE NBPs include the determination of breast shape parameters consistent with a prone imaging position [32], [33], the stochastic assignment of tissue-specific acoustic properties (density, SOS, and AA), as well as the modeling of acoustic heterogeneity within fatty and glandular tissues.

To demonstrate the usefulness of the proposed computational framework, two case studies are presented. Case study 1 assesses the reconstructed SOS IQ using different compensation techniques to account for unknown AA. Case study 2 demonstrates the utility of the proposed framework for generating large-scale ensembles of NBPs for the training of deep learning-based USCT reconstruction methods. To accompany this work, a python library implementing the proposed approach for the generation of 3-D acoustic phantoms has been made publicly available under GPL-2.0 [34]. Furthermore, two datasets have been publicly released under CC-0: the first consists of 52 2-D breast phantom slices and corresponding USCT measurement data [35] and the second contains 4 3-D realizations of NBPs [36].

The remainder of this article is organized as follows. In Section II, the background on USCT breast imaging and the FDA VICTRE project are provided. The stochastic generation of 3-D anatomically and physiologically realistic NBPs for USCT VITs is introduced in Section III. Several examples of NBPs generated with the proposed tool are presented in Section IV. Section V contains the case studies that illustrate possible applications of the proposed phantoms to inform image reconstruction development. Finally, in Section VI, a discussion of the wide range of applications enabled by the proposed framework is provided.

Background

USCT Breast Imaging

In recent decades, a number of research groups have been developing USCT imaging technologies for breast imaging applications [3], [5], [37], [38]. In a typical breast USCT system, the patient lies prone on the imaging table and the breast to be imaged is submerged in water. An array of ultrasound transducers surrounds the breast. Each transducer emits an acoustic pulse one by one until the breast is insonified from all directions. During each shot, all other transducers act as receivers, recording the transmitted, scattered, and reflected wavefield data.

Three types of USCT images are conventionally produced: reflectivity, SOS, and AA [6]. Reflectivity images can be reconstructed by use of integral geometry-based approaches that are similar to the delay-and-sum methods widely employed in conventional B-mode imaging. The majority of the SOS and AA reconstruction methods investigated to date are generally based on two categories: approximated wave equation methods [16], [26], [39] and full-waveform inversion (FWI) methods [1], [2], [15], [27]. Because FWI methods take high-order refraction and diffraction effects into account, they can produce images that possess higher spatial resolution images than those produced by use of linearized or approximate methods [1], [2], [40]. However, FWI is computationally expensive and memory burdensome, especially for 3-D problems, thus hampering the widespread application of FWI to USCT breast imaging. Moreover, FWI suffers from the so-called cycle skipping phenomenon [41], thus requiring an accurate initial estimate of the SOS map to ensure convergence to a useful solution. As a result, there is still an imperative need to systematically investigate and optimize USCT reconstruction methods by means of computer-simulation studies.

Description of VICTRE

The VICTRE project of the FDA has recently released a series of software tools to provide a complete simulated imaging chain for mammography and digital breast tomosynthesis [17]. The VICTRE software includes open-source tools to generate the 3-D random anthropomorphic voxelized phantoms of the human female breast [42]. Using this tool, large ensembles of anthropomorphic NBPs with realistic anatomical structures can be generated by specifying different virtual-patient characteristics that include breast type, shape, granularity, density, and size. By appropriate selection of physical attributes and material coefficients, the VICTRE NBPs can be customized for particular imaging tasks.

The VICTRE software generates NBPs corresponding to the four different levels of breast density defined according to the American College of Radiology’s (ACR) Breast Imaging Reporting and Data System (BI-RADS) [43]: A) breast is almost entirely fat; B) breast has scattered areas of fibroglandular density; C) breast is heterogeneously dense; and D) breast is extremely dense. Each NBP is a 3-D voxelized map consisting of ten tissue types: fat, skin, glandular, nipple, ligament (connective tissue), muscle, terminal duct lobular unit, duct, artery, and vein. Large ensembles of stochastic NBPs with realistic variability in breast volume, shape, fraction of glandular tissue, ligament orientation, and tissue anatomy can be generated by controlling input parameters and selecting the random seed number.

In addition, the VICTRE projects include tools to generate 3-D numerical lesion phantoms (NLPs), which can be inserted into the NBPs at clinically plausible locations [44]. Two types of lesions, microcalcification clusters and spiculated masses, can be generated. The size and shape of the lesions can be customized. An example of anatomically realistic NBP and NLP generated using the VICTRE tools is shown in Fig. 1.

Fig. 1.

(a) Volume rendering of fatty breast phantom: partial transparencies are used to highlight anatomical structures and cross-sectional view of this 3-D breast phantom. (b) Volume rendering of spiculated lesion phantom and cross-sectional view of this 3-D lesion phantom.

There exist several challenges that must be addressed in order to extend the VICTRE project to produce NBPs for use in VITs of USCT technologies. These include determination of breast shape parameters consistent with a prone imaging position, the stochastic assignment of tissue-specific acoustic properties (density, SOS, and AA), and the modeling of acoustic heterogeneity within fatty and glandular tissues.

Methods

Several adaptations and customizations of the VICTRE tools were developed that will enable the generation of large ensembles of acoustic NBPs that display clinically relevant variability in both anatomical structures and acoustic properties. The specific procedures for accomplishing this are described in the following.

Generation of Anatomically Realistic Realizations of NBPs and Lesion(s) Insertion

The goal of this step is to generate large ensembles of anatomically realistic NBPs representing four different types of breast (extremely dense, heterogeneously dense, scattered fibroglandular, and fatty). Section III-A1 describes how shape and deformation parameters in the VICTRE NBPs can be set to generate virtual patients with variable breast sizes that are representative of a clinical population and shapes that are consistent with USCT imaging protocols. Section III-A2 describes adaptations to the internal anatomical structures of the NBP to exclude tissues that are not relevant for USCT applications. Finally, Section III-A3 describes how one or more lesions are optionally inserted into the NBPs.

Breast Shape and Deformation Parameters:

Appropriate distributions of breast size parameters were determined for each breast type based on clinical data [45]. In the VICTRE software, the shape of the breast is created by applying a series of transformations to a base superquadratic surface. A detailed description of the breast shape model was presented in [46]. Here, the main parameters affecting size and shape of the breast are discussed. As shown in Fig. 2, the parameters a1, a1 , a2 , and a2 adjust the breast volume in the top, bottom, left, and right hemispheres, respectively. The parameter a3 adjusts the length of the breast. Fig. 3 shows how other parameters affect the final shape of the breast. The parameter z1 is the quadric shape exponent along the polar angle. The ptosis deformation parameters B0 and B1 model the sagging that affects a breast as a subject age. Finally, the turn-pop deformation parameters H0 and H1 change the shape of the top half of the breast laterally. This deformation allows the top part of the virtual breast to point toward the shoulder. The probability distributions assigned to these parameters are summarized in Table I and were set to be consistent with the patient lying prone on the examination table. Here—among all possible distributions with a specified mean, variance, and having support in a bounded interval—a truncated Gaussian distribution is chosen since it represents the maximum entropy distribution that satisfies such constraints.

Fig. 2.

Overview of size parameters: a, a, a, a, and a3.

Fig. 3.

Overview of deformation parameters. Red breast: hemispherical breast phantoms without deformation. Yellow breast: deformed breast phantoms. Left: effect of superquadric exponent deformation (ϵ1) in sagittal plane. Center: effect of ptosis deformation (B0, B1) in sagittal plane. Right: effect of turn top deformation (H0, H1) in coronal plane.

TABLE I

Shape and Size Parameters

Parameter	Types A-C	Type D
a1t (cm)	TN (5.85, 2.3275, 3.85, 7.70)	TN (4.20,1.225,2.80,5.25)
a1b/a1t	N (1, 0.02)	N (1, 0.02)
a2r/a1t	N (1, 0.05)	N (1, 0.05)
a2l/a2r	N (1, 0,05)	N (1, 0,05)
a3/a1t	TN (1.48, 0.18, 1, 1.6)	TN (1.22,0.1,0.75.1.5)
ϵ 1	N (1, 0.1)
B 0	TN (0, 0.1, −0.18, 0.18)
B 1	TN (0, 0.1, −0.18, 0.18)
H 0	TN (0, 0.15, −0.11, 0.11)
H 1	TN (0, 0.25, −0.3, 0.3)

: Gaussian distribution with mean μ and standard deviation σ.

Truncated Gaussian distribution in interval (a, b).

Relabeling of Tissue Types Invisible to USCT:

The generated NBPs are high-resolution volumes with a voxel size as small as 50 μm. Each voxel is assigned a label corresponding to one of the ten tissue types (fat, skin, glandular, nipple, ligament, muscle, terminal duct lobular unit, duct, artery, and vein). Of these tissues, only four are typically visible in USCT imaging: fatty, glandular, skin, and ligament. Voxels corresponding to tissue types for which there is not enough clinical evidence that they can be well resolved in USCT imaging are relabeled as fatty or glandular based on the type of the neighboring voxels. An ad hoc inpainting algorithm was designed to ensure consistent anatomical structures when relabeling voxels. The first step in the algorithm marks all voxels to be relabeled. Marked voxels are assigned to regions based on connectivity (two voxels are connected if they share a face) and process each connected region independently. For each region, the algorithm selects voxels near the boundary of the region (i.e., all voxels that share at least one face with unmarked voxels), reassigns their labels to the most occurring label among those of neighboring (unmarked) voxels, and unmarks them. This step is repeated until all voxels in all regions have been relabeled. An example of the result of replacement of USCT-invisible tissues is shown in Fig. 4.

Fig. 4.

Illustration of relabeling of tissues type invisible to USCT. (a) Anatomical phantom (tissue labels) generated by VICTRE. (b) Phantom after tissue relabeling. The different colors represent distinct tissue types.

Lesion Insertion:

To generate NBPs that contain tumors, synthetic lesions can be inserted in the healthy NBPs as follows. First, an ensemble of numerical tumor phantoms (NTPs) with various sizes and irregular (spiculated) shapes can be generated by the use of the VICTRE tool. One or more NTPs can then be inserted in each NBP at locations among those suggested by the VICTRE phantom tools as candidate tumor locations. Additional location constraints are included to ensure that tumors do not overlap each other or skin layer and are not inserted too close to the chest or nipple.

Assignment of Acoustic Properties

By use of the anatomical breast maps generated in Section III-A, 3-D acoustic NBPs can be established via stochastic assignment of acoustic properties. The acoustic properties considered are the SOS c (m/s), density ρ (kg/m3), and AA coefficient α0 (Np/m/MHz) with power-law exponent y. The 3-D acoustic property maps are constructed as follows. First, acoustic properties values are stochastically assigned to each phantom voxel based on the tissue type label as described in Section III-B1. Next, to model variations in the acoustic properties across voxels of the same tissue type, SOS and density maps are perturbed by additive colored noise with a prescribed correlation structure as described in Section III-B2. Finally, the choice of the power-law exponent y is presented in Section III-B3.

Stochastic Assignment of Acoustic Properties to Each Tissue Type:

Acoustic properties (SOS, AA, and density) are assigned to each voxel of the anatomical NBPs generated in Section III-A as follows. For each tissue type, values of SOS, AA, and density are sampled from a predefined probability distribution and assigned to all voxels of that tissue type. Table II shows the probability distributions of the acoustic parameters assigned to each tissue type. These were chosen based on a comprehensive literature survey to represent anatomically realistic values. The SOS values of healthy breast tissues were based on the clinical studies reported in [10] and [49]. The distributions of density and AA in healthy breast tissues were set according to [48], a database providing comprehensive estimates of material properties of several human tissues, as well as statistical information about the spread of those properties. This information was based on a meta-analysis of over 150 references. The variance of AA values for each tissue type was set to 10% of the respective mean values. Finally, tumor acoustic properties were also chosen from clinical literature of breast pathology [4], [50], [51].

TABLE II

Acoustic Property Values of Different Tissue Types

Medium	SOS [m/s]	AA [N p/m H zy]	Density [k g/m3]
Water	1500 @ 26°C [47]	0.025328436023 [48]	994 [48]
Skin	TN (1555.0, 10.0, 1530, 1580) [10]	N (21.158, 2.16) [48]	TN (1109, 14, 1100, 1125) [48]
Fat	TN (1440.2, 20.9, 1412, 1485) [10], [49]	N (4.3578, 0.436) [48]	TN (911, 53, 812, 961) [48]
Glandular	TN (1540.0, 15.0, 1517, 1567) [10], [49]	N (8.635, 0.86) [48]	TN (1041, 45.3, 990, 1092) [48]
Ligament	TN (1457, 18.5, 1422, 1565) [10], [49]	N (14.506, 1.45) [48]	TN (1142, 45, 1110, 1174) [48]
Tumor	TN (1548, 10.3, 1531, 1565) [4]	N (31, 2.3)	TN (945, 20, 911, 999) [51]

Upon completion of this step, piecewise constant acoustic maps are constructed that present variability both in their values, which are randomly sampled, and spatial distribution, which is dictated by the NBP stochastic anatomical structure. Fig. 5(a) shows an example of a slice through a piecewise constant 3-D SOS phantom generated by the described procedure.

Fig. 5.

Illustration of texture generation on SOS phantoms. (a) Piecewise constant SOS phantom. (b) SOS phantom after texture generation.

Modeling Spatial Heterogeneity Within Fatty and Glandular Tissues:

Acoustic scattering in breast tissues arises not only from jumps in acoustic impedance across tissue types but also from spatial heterogeneity within each tissue [52]. The latter is a predominant effect in fatty and glandular tissues. To account for the spatial heterogeneity within these tissues, random textures are introduced into the SOS and density maps. SOS and density textures in glandular tissue are modeled as a spatially correlated Gaussian random field with zero mean and Gaussian covariance function. SOS and density textures in fatty tissue are modeled as truncated (plus or minus 0.9 standard deviations) spatially correlated Gaussian random field with zero mean and Gaussian covariance function account for the lower acoustic scattering observed in fatty tissues [53]. The pointwise standard deviations σ and correlation lengths ℓ are shown in Table III and are based on reflectivity tomography studies [53]. SOS and density textures are sampled independently one from the other. Each voxel in the generated textures maps is added to the corresponding voxel in the piecewise constant property maps described in the preceding paragraph; this results in NBPs that display random heterogeneity with the glandular and fatty tissues. Fig. 5(b) shows an example of a slice through a 3-D SOS phantom that contains tissue texture generated by the described procedure. Note that acoustic heterogeneity is stronger in glandular tissue (gray regions) than in fatty tissue (black regions).

TABLE III

Pointwise Standard Deviations and Correlation Lengths Uses to Model Texture in Fatty and Glandular Tissues

Property	SOS		Density
	σ (m/s)	ℓ (mm)	σ (kg/m3)	ℓ (mm)
Fatty tissue*	28.8	0.21	18.22	0.21
Glandular tissue	30.4	0.21	20.82	0.21

Random texture in fatty tissues is truncated within the ±0.9σ range.

Power-Law Attenuation Model:

To model frequency dependence in AA, a fractional power-law model [54] is assumed. Specifically, frequency-dependent AA α (Np/m) is defined as where α0 (Np/m/MHz−1) is the AA coefficient, y is the fractional power-law exponent, and f (MHz) is the acoustic wave frequency. In general, the exponent y varies for different tissue types and estimates for several breast tissues can be found in the IT’IS database [48]. However, several widely employed time-domain wave propagation solvers [55], [56] assume a spatially homogeneous exponent y.

To address this, a homogenization technique based on the solution of a nonlinear least squares problem is proposed. The proposed technique considers wave propagation in one spatial dimension for which an analytical model of AA can be constructed. Specifically, for a monochromatic wave with frequency f propagating through a heterogeneous medium with thickness L, the log amplitude ratio between the transmitted A and incident A wave is where is the tissue-dependent fractional power-law exponent. Since attenuation in water is negligible and the volume of skin, tumor, and ligament tissues is small compared to the whole breast, a medium consisting of only fatty and glandular tissues is considered. Under this simplifying assumption, is determined as a function of the fatty tissue volume fraction only. The spatially homogeneous fractional power-law exponent y is then defined as where is the average value of α0(x) and the frequencies f (k = 1, . . . , K) are uniformly distributed over the range of frequencies typically employed in USCT imaging. Table IV reports the estimated power-law exponent y as a function of the fatty tissue volume fraction vfat when K = 22 frequencies evenly spaced between 0.2 and 2.3 MHz are used to evaluate (3).

TABLE IV

Homogeneous Power-Law Exponent y as a Function of the Fatty Tissue Volume Fraction V

Breast type	A	B	C	D
v fat	~95%	~85%	~66%	~40%
y	1.1151	1.1642	1.2563	1.3635

Examples of Generated NBPs

Fig. 6 shows four 3-D visualization examples, one for each breast type, of 3-D acoustic NBPs produced by the proposed framework. The generation of the anatomical structures for these NBPs using the VICTRE tools took about 80 ~ 240 mins on a single node of the Golub cluster at the University of Illinois at Urbana–Champaign campus cluster (Two ten-core Intel E5–2670v2 CPUs and 64 GB of memory per node). Tissue relabeling, tumor insertion, and assignment of spatially varying acoustic properties took between 20 and 50 mins on the same machine. Time variations depended on the volume of phantoms. Paraview [57] was used for volume rendering to highlight internal structures. Note the variability in size, shape, internal structures, and values of acoustic properties among the four NBPs. Fig. 7 shows the examples of 2-D cross-sectional slices extracted from the phantoms, one for each breast type. Yellow rectangles indicate the location of the inset zoom region where a lesion was inserted. The diameters of the inserted lesions were sampled from a uniform distribution between 1.5 and 5 mm, to mimic small lesions in early breast cancer.

Fig. 6.

3-D rendering of acoustic phantoms from four breast types. From up to bottom: (a) almost entirely fatty, (b) scattered areas of fibroglandular density, (c) heterogeneously dense, and (d) extremely dense. From left to right: the SOS (m/s), AA (Np/m/MHz), and density (kg/m3) volumes.

Fig. 7.

Realization of a cross-sectional slice from four types breast. From left to right: the SOS (mm/μs) image, AA (Np/m/MHz) image, and density (kg/m3) image. Tumor region is zoomed in. From up to bottom: (a) almost entirely fatty, (b) scattered areas of fibroglandular density, (c) heterogeneously dense, and (d) extremely dense.

Fig. 8 compares the distributions of the breast diameter and depth in a virtual population of 1000 NBPs to that observed in a sample of 219 women with age ranging between 35 and 82 years and median age of 54 years [45]. The proportion for each breast type in the virtual population was set to 10% for breast types A and D and 40% for breast types B and C [43]. The figure shows good qualitative agreement in the diameter and depth distributions between the virtual population and the clinical sample. It is worth noting that the distributions of the virtual population are skewed toward slightly larger breast sizes compared to those of the clinical sample. This is intentional and aims to address a limitation of the sample in [45], which is biased toward denser—and therefore smaller—breast types (23% type A, 40% type B, 28% type C, and 9% type D).

Fig. 8.

Breast size distribution comparison. (a) Breast diameter distributions. (b) Depth distributions. Blue: distributions of the generated NBPs. Orange: distributions estimated from clinical data [45].

Case Studies

Two case studies were conducted to demonstrate the usefulness of the proposed framework for generating acoustic NBPs. Case study 1 (Section V-C) assesses reconstructed SOS IQ when heuristic procedures for compensating for unknown AA are employed. Case study 2 (Section V-D) demonstrates the utility of the proposed framework for the training and assessment of deep learning-based USCT reconstruction methods. In both studies, 2-D cross-sectional slices extracted from the 3-D NBPs are (virtually) imaged using the stylized 2-D imaging system described in Section V-A.

Virtual Imaging System

A stylized 2-D virtual imaging system was modeled to generate USCT measurement data. It comprised 1024 idealized, point-like, transducers that were evenly arranged in a circular array with a radius of 110 mm. The excitation pulse employed in this study was assumed to be spatially localized at the emitter location. The central frequency and duration of the pulse were set to 1 MHz and 10 μs, respectively. The pulse profile s(t) was defined as the sum of three sinusoidal functions tapered by a Gaussian kernel as where σ = 1.6 μs is the standard deviation of the Gaussian kernel, t = 3.2 μs is a constant time shift, and f = 1 MHz is the central frequency. The maximum frequency of s(t) is 2.3 MHz.

Cross-sectional slices were extracted from the 3-D NBPs and centered within the field of view of the imaging system. Bilinear interpolation was employed to downsample the maps of acoustic properties to a computational grid comprised of 0.1 mm isotropic pixels. To emulate the imaging process, the propagation of the pressure waves through the object was modeled by solving the lossy acoustic wave equation with power-law frequency-dependent AA [58] by use of a time-explicit pseudo-spectral k-space method [59]–[61]. Further details regarding the wave solver and its implementation are presented in Appendix C. The simulated measurement data were corrupted with Gaussian independent and identically distributed (i.i.d.) noise that had zero mean and a standard deviation of 0.02% of the maximum pressure amplitude at the emitting transducer.

Computation of USCT measurement data for a single slice took about 110 GPU hours using a single NVIDIA GK110 Kepler GPU on the Blue Water cluster at the National Center for Supercomputing Applications.

SOS Images Reconstructed Under Favorable Conditions

Reconstruction of SOS images under favorable conditions (namely, AA map known exactly) is considered here. However, this study does not represent an inverse crime as it includes three sources of model mismatch: 1) a constant density map was employed in the reconstruction; 2) measurement data were corrupted with additive Gaussian noise; and 3) reconstruction was performed on a coarser grid. By use of the procedures described above, 2-D slices from 52 NBPs (13 for each of the 4 breast types) were extracted and virtually imaged to produce USCT measurement data. From these data, SOS images were reconstructed on a grid with pixel size of 0.2 mm by use of a previously published waveform inversion with source encoding (WISE) method [2]. The reconstruction was initialized by use of a blurred version (Gaussian blur with 8-mm correlation length) of the true SOS. Because true values of AA were considered here, the reconstructed SOS estimates are expected to generally be of higher quality than would be obtained if attenuation properties had to be concurrently estimated with the SOS or if incorrect fixed values of AA were employed. In this sense, it will be useful to compare these reconstructed SOS estimates against the images reconstructed in the two case studies mentioned next.

Fig. 9 presents the examples, one for each breast type, of the ground truth and reconstructed SOS images assuming that the AA distribution and density are known. Table V reports the average mean square error (MSE)[1] and structural similarity index measure (SSIM) [62] for each breast type.

Fig. 9.

Reference images: ground truth (top row) and reconstructed (bottom row) SOS maps. From left to right: breast type A-D. The unit is mm/μs.

TABLE V

Reference SOS Reconstructions: MSE and SSIM

Breast type	MSE (std)	SSIM (std)
A	1.786e-04 (3.923e-5)	0.9835 (0.0056)
B	2.571e-04 (1.087e-4)	0.9788 (0.0069)
C	2.459e-04 (1.797e-4)	0.9732 (0.0102)
D	2.258e-04 (1.301e-4)	0.9835 (0.0066)
all types	2.269e-04 (1.210e-04)	0.9799 (0.0083)

Case Study 1: Heuristic Compensation of AA in SOS Reconstruction

In this study, two heuristic approaches to compensating for AA when reconstruction SOS estimates were compared: a two-region attenuation model (TRAM) and a data domain attenuation compensation (DDAC). The TRAM assumes that the breast boundary is known (reflectivity imaging could be possibly used to estimate it) and assigns one constant AA value to the water bath (α0 = 0) and another to the breast region. The attenuation coefficient of the breast region was set to 5.20 [Np/m/MHz] , which corresponds to a weighted average (80%–20% split) of the mean values of AA in fatty and glandular tissues, as reported in Table II. The heuristic DDAC procedure seeks to compensate for AA by modifying the amplitudes of the recorded pressure data, rather than explicitly modeling attenuation in the wave propagation forward model. Specifically, for each pair of emitting/receiving transducers, the maximum amplitude of the recorded signal was rescaled to match that of the corresponding measurement when only the water bath was present [63]. The generation of synthetic data and reconstruction method used in this case study is the same as described in Section V-B.

Fig. 10 shows the examples of reconstructed images of four breast types using the proposed techniques (TRAM and DDAC) to compensate for unknown AA properties. The corresponding ground truth images and reference reconstructions are shown in Fig. 9. Table VI shows the quantitative evaluations of all reconstructed images on each breast type from A to D. Fig. 11 shows the variation with respect to MSE and SSIM in all image samples and breast types A–D. In all the cases tested but 9, TRAM led to measurable improvement (i.e., smaller MSE) compared to DDAC. Remarkably, of the nine cases in which DDAC led to a smaller MSE, six were for breast type C. Furthermore, TRAM achieved an MSE smaller than that of the reference reconstruction (Section V-B) for 25 out of the 52 tested cases, thus suggesting that TRAM is an effective approach to compensate for unknown AA in FWI reconstruction of SOS maps.

Fig. 10.

Case study 1: reconstructed SOS images corresponding to the same phantoms shown in Fig. 9 using TRAM (top row) and DDAC (bottom row). From left to right: breast type A-D. The unit is mm/μs.

TABLE VI

Case Study 1: MSE and SSIM of SOS Images Using the TRAM and DDAC Approaches to Compensate for Unknown AA

		A	B	C	D	All types
MSE (std)	TRAM	1.825e-4 (3.863e-5)	2.853e-4 (1.199e-4)	2.686e-4 (1.449e-4)	2.300e-4 (2.006e-4)	2.416e-4 (1.414e-4)
	DDAC	3.132e-04 (7.309e-5)	3.521e-4 (1.210e-4)	2.668e-4 (1.666e-4)	2.720e-4 (2.059e-4)	3.006e-4 (1.513e-4)
SSIM (std)	TRAM	0.9819 (0.0080)	0.9766 (0.0083)	0.9700 (0.0126)	0.9843 (0.0066)	0.9785 (0.0103)
	DDAC	0.9777 (0.0058)	0.9732 (0.0073)	0.9709 (0.0094)	0.9818 (0.0071)	0.9761 (0.0085)

Fig. 11.

Case study 1: boxplots of MSE and SSIM value with respect to TRAM and DDAC. From left to right: breast types A-D and all breast types together.

Case Study 2: Deep Learning Reconstruction Method

There remains an important need to lessen the computational burden of FWI. A supervised learning-based method is proposed to reduce the number of FWI iterations and drastically lessen the computational burden, as well as enhancing the reconstruction by using ground truth images as training data. Fig. 12 shows the proposed learning framework, in which a deep neural network (a five-level U-Net [64]) is trained to minimize the MSE of the reconstructed SOS images. The input to the network is an intermediate SOS estimate obtained by early stopping of the WISE method after 35 iterations. The rationale of this method is that early stopped reconstructed images capture structural information of the SOS map but lack quantitative accuracy. The network was trained for 220 epochs on a dataset consisting of 622 2-D slices, each extracted from a different NBP (312 type B NBPs and 310 type C NBPs). The Adam optimizer was used with a batch size equal to 32 and the initial learning rate 0.001. The learning rate was reduced by a factor of 0.9 after each epoch. Several models corresponding to different architecture hyper-parameters (e.g., number of layers at each level) were trained. The selected model was that achieving the highest mean MSE on the validation set consisting of 100 slices (equally split between types B and C).

Fig. 12.

Case study 2: supervised deep learning framework for SOS reconstruction.

The testing set consisted of the 52 2-D phantoms described in Section V-B, thus allowing us to evaluate the accuracy of the network for both in-distribution (types B and C) and out-of-distribution data (types A and D).

Fig. 13 reports the examples of learned reconstructed images of four breast types. The top row corresponds to the early stopped WISE reconstruction after 35 iterations, and the bottom row shows the output of the neural network. The corresponding ground truth images and reference reconstructions are shown in Fig. 9. The proposed learning approach improved the visual quality of the images, leading to sharper tissue interfaces. Table VII and Fig. 14 show the quantitative evaluations on the test dataset. The reported MSE and SSIM values are stratified by breast types: breast types A and D (out of distribution) have a larger median MSE and smaller median SSIM than types B and C (in distribution). While the reported MSE and SSIM are comparable (or sometimes even better) than those reported for the reference reconstructions in Table V, the learned reconstruction method may mistakenly introduce some fine structures (hallucinations[2]) that are not existing in the ground truth image [65]. An example is shown in Fig. 15. This case study demonstrates that, while deep learning methods can be used to enhance perceived and quantitative IQ, their results must be interpreted with particular care due to the possibility of introducing hallucinated features in the image.

Fig. 13.

Case study 2: reconstructed SOS images of the phantoms shown in Fig. 9 using a machine learning-based method. Top row: input to the neural network. Bottom row: corresponding estimated image. From left to right: breast type A–D. The units are mm/μs.

TABLE VII

Case Study 2: MSE and SSIM

Breast type	MSE (std)	SSIM (std)
A	2.165e-04 (8.413e-5)	0.9675 (0.0081)
B	1.973e-04 (7.9898e-5)	0.9707 (0.0079)
C	2.160e-04 (1.149e-4)	0.9788 (0.0069)
D	2.887e-04 (2.009e-4)	0.9651 (0.0073)

Fig. 14.

Case study 2: boxplot of MSE and SSIM value of learned reconstructed results for breast type A–D. The subscript ∗ denotes out of distribution breast types.

Fig. 15.

Case study 2: false structures in the reconstructed image. (a) Ground truth image. (b) Neural network-based reconstructed image that contains small hallucinated features.

Conclusion

In this work, procedures were established by which 3-D NBPs can be computed for use in large-scale VITs of 2-D or 3-D breast USCT. This will, for the first time, permit 3-D realistic NBPs to be computed that possess varying shapes, acoustic properties, tissue texture, and tumors. This was accomplished by adapting VICTRE tools to USCT imaging. While some modeling choices and simplifications had to be made, the modular and flexible implementation of the phantom generation procedures allows for additional customizations of the NBPs. For example, future studies may include additional tissue type, use of different lesion models, analyze the detectability of microcalcifications, or develop advanced biomechanical models to capture deformations of the submerged breast due to buoyancy. In summary, the generated NBPs improve the authenticity of USCT virtual imaging studies and can be employed widely for the investigation of advanced image reconstruction methods, objective evaluation of the USCT breast imaging systems, and the development of machine learning-based methods.

TABLE VIII

Imaging System Discretization Parameters

Computational grid	2560 by 2560 pixels(0.1 mm pixel size, ~6.5 ppw)
Time step size	1/50 μs, CFL number =0.3
Simulation time	170 μs, 8500 time steps