1.

Introduction

Optoacoustic computed tomography (OAT), also known as photoacoustic computed tomography, is a non-invasive imaging modality being actively developed for clinical breast imaging and other biomedical applications.^1–8 A unique feature of OAT is the ability to produce an image based on the endogenous optical contrast associated with chromophore concentrations and oxygenation states within tissue, without ionizing radiation and the loss of spatial resolution typically related to purely optical techniques such as diffuse optical tomography.^1,9 This permits the imaging of tissue metabolism and angiogenesis, which have been identified to play a critical role in tumor growth and progression.^7,10 Therefore, optoacoustic imaging is ideally positioned to resolve these two hallmarks of cancer in vivo.^2–8,10 As such, an optimized and validated OAT system can be a powerful tool for the management of breast cancer. By assessing the tumor microvasculature density and blood oxygenation, it can enable the initial evaluation of tumor aggressiveness to inform the treatment plan and prognosis. It also allows for the monitoring of tumor response to treatment over time. However, to realize its full diagnostic potential, OAT should be equipped with the capability of providing quantitative information on true values of the optical absorption coefficient, which is proportional to molecular concentrations.^7,11,12

A large number of different system designs for three-dimensional (3D) breast OAT that deploy varying light delivery and acoustic detection strategies have been proposed.^3–6,13 This is unlike in X-ray mammography, breast magnetic resonance imaging (MRI), and breast ultrasound, in which similar implementations are typically in use per modality. Ideally, candidate designs of new quantitative OAT systems would be evaluated based on clinically relevant objective image quality measures via human subject studies. However, this is not a feasible solution given the large space of possible (technical trade-offs considered) design parameters of OAT imaging systems, the large variety in breast sizes and compositions, and the cost and potential ethical concerns associated with such studies. Instead of clinical trials, virtual imaging trials (VITs), i.e., computer-simulation studies, have been advocated for assessment and optimization of system and algorithm designs in the early stages of technology development. For VITs to be clinically relevant, realistic numerical phantoms must be employed as the to-be-imaged objects.^14–20 Moreover and importantly, because imaging technologies are not typically optimized for use with only a single subject, ensembles of different phantoms that possess clinically relevant variability in anatomy and tissue properties must be virtually imaged. In this way, ensemble-averaged objective measures of image quality can be computed.²¹

Several numerical breast phantoms (NBPs) and numerical lesion phantoms (NLPs) have been proposed.^16–20,22 However, most of the existing NBPs and NLPs were created from a limited number of clinical data, resulting in a lack of variability or oversimplified anatomical structures.^16–19 The virtual imaging clinical trials for regulatory evaluation (VICTRE) project at the U.S. Food and Drug Administration (FDA)¹⁴ released software tools to generate realistic NBPs and NLPs for use in mammography applications. Although the VICTRE NBPs are considerably realistic for VITs of mammography in terms of principal tissue compositions, significant improvements are required for use in VITs of breast OAT. Bao et al.²⁰ reported optoacoustic NBPs based on the VICTRE tools. However, these phantoms do not introduce adaptations of the blood vascular network for use in VITs of OAT and do not account for physiological variability in tissue optoacoustic properties, the values of which are fixed and deterministic for each tissue type.

Although simulation tools for photon transport and acoustic wave propagation^23–25 in general media with spatially varying voxel-wise properties are available, the reported studies that employed NBPs did not leverage this capability. Not only in Ref. 20 but also in the other previous studies,^16–19,22 piecewise-constant optical properties were assigned to each tissue type in the NBPs and NLPs without taking into account spatial variability in the property distribution induced by oxygen transport between tissues. In addition, a correlation between optical properties and both chromophore concentrations and the oxygenation state of the tissue has not been addressed in the existing NBPs and NLPs.^16–20,22 Therefore, further enhancements are needed to establish physiologically realistic NBPs and NLPs for use in VITs of OAT.

In this work, an end-to-end framework for producing ensembles of realistic 3D anatomical, functional, optical, and acoustic NBPs and NLPs for use in VITs of quantitative OAT of breast cancer is established. The generation of the anatomical NBPs is accomplished by extending the VICTRE NBPs with blood vasculature modifications for use in VITs of OAT. Tissue composition within the anatomical NLPs is modeled according to different lesion types and aggressiveness. Depending on the type of cancer of interest, the lesions can contain a necrotic core and/or a tumor angiogenesis region surrounding the viable tumor cell region, which can be modeled because of the proposed modifications to the VICTRE NLPs. Realistic functional, optical, and acoustic properties are stochastically assigned for each breast tissue. The optical absorption coefficient is modeled based on the concentrations of the primary chromophores considered, which are stochastically chosen within their physiological ranges. Our modified version of the VICTRE tool code package,²⁶ which is a fork of the original software, has been released under the same creative commons zero license (CC0) used by the original project. Furthermore, the software framework to generate stochastic distributions of the functional, optical, and acoustic properties based on the modified VICTRE anatomical phantoms has been made publicly available under the GNU general public license version 3 (GPLv3). The software, named stochastic optoacoustic NBP (SOA-NBP),²⁷ is implemented in Python, and it includes capabilities for inserting superficial blood vasculature under the skin layer; creating anatomical NLPs; and assigning functional, optical, and acoustic properties to each breast tissue type. To enable researchers to immediately benefit from this work, 84 datasets were released under the CC0; these consist of anatomical, functional, optical, and acoustic NBPs and the corresponding simulated multi-wavelength optical fluence, initial pressure, and OAT measurement data.^28–30

The remainder of this article is organized as follows. Relevant background information including a description of the VICTRE and previous NBPs for OAT is provided in Sec. 2. The proposed framework is described in Sec. 3, and several examples of NBPs generated employing the proposed framework are presented in Sec. 4. A case study that investigates the impact of acoustic breast heterogeneity on image reconstruction quality is provided in Sec. 5, as an example of how the proposed stochastic phantoms can enable important studies. The article concludes with a summary and discussion in Sec. 6.

2.

Background

3.

Methods

To circumvent the limitations described above, adaptations and customizations of the VICTRE NBPs were developed to enable the generation of large ensembles of realistic optoacoustic NBPs that exhibit clinically relevant variability in anatomical structures and functional, optical, and acoustic properties. The flow chart of optoacoustic NBP generation is illustrated in Fig. 1, and details of each step are explained in the following sections.

Fig. 1

Generation of optoacoustic NBPs.

3.1.

Generation of Realistic Anatomical NBPs and Lesion Insertion

In this first step, the construction of ensembles of realistic anatomical NBPs for different BI-RADS breast types [(A) breast is almost entirely fatty, (B) breast has scattered areas of fibroglandular density, (C) breast is heterogeneously dense, and (D) breast is extremely dense] using a modified version of the VICTRE tool is described. Additionally, the generation and insertion of anatomical NLPs into the NBPs are described.

3.1.1.

Definition of breast size and shape

Based on clinical data and constraints by designs of the existing OAT breast imaging systems^3–6 (where a breast should fit within a scanning radius of 85 mm), the distributions of breast size parameters were determined for each BI-RADS breast type.³⁸ The VICTRE tool generates anatomical NBPs according to the parameters specified in the configuration file. Among the parameters, the breast volume extent parameters $a_{1t}$, $a_{1b}$, $a_{2l}$, $a_{2r}$, and $a_3$ and the breast shape parameters ${ϵ}_1$, $B_0$, $B_1$, $H_0$, and $H_1$ were specified in the VICTRE configuration file. These parameters are explained in Ref. 15, and the probability distributions for the parameters that are consistent with the patient’s position during a 3D OAT scan are summarized in Table 2. Here, a truncated Gaussian distribution ($TN$) was chosen among the possible distributions because it is the maximum entropy distribution supported in a bounded interval with a specified mean and standard deviation.³⁹

Table 2

Shape and size parameters.

Parameter	Types A and B	Type C	Type D
$a_{1t}$ (mm)	$\mathrm{TN}(59.70,3.58,50.77,71.5)$		$\mathrm{TN}(50.05,3.58,42.9,57.2)$
$a_{1b}/a_{1t}$	$N(1,0.02)$
$a_{2r}/a_{1t}$	$N(1,0.05)$
$a_{2l}/a_{2r}$	$N(1,0.05)$
$a_3/a_{1t}$	$\mathrm{TN}(0.85,0.14,0.8,1.2)$	$\mathrm{TN}(0.85,0.12,0.7,1.1)$	$\mathrm{TN}(0.85,0.1,0.7,1.1)$
${ϵ}_1$	$N(1,0.1)$
$B_0$	$\mathrm{TN}(0,0.1,-0.18,0.18)$
$B_1$	$\mathrm{TN}(0,0.1,-0.18,0.18)$
$H_0$	$\mathrm{TN}(0,0.15,-0.11,0.11)$
$H_1$	$\mathrm{TN}(0,0.25,-0.3,0.3)$

N(μ,σ): Gaussian distribution with mean μ and standard deviation σ.TN(μ,σ,a,b): truncated Gaussian distribution in interval (a,b).For hemispherical shapes (a1t=a1b=a2r=a2l=a3), the ϵ1 value is set to “1,” and the B0, B1, H0, and H1 values are set to “0.”

The generated anatomical NBPs were discretized in three dimensions using a uniform Cartesian lattice. A voxel size of 0.125 mm was chosen for the discretization, considering computation time and memory usage. This choice allows for avoiding the discretization inverse crime when performing image reconstruction as the typical resolution in OAT breast images is between 0.2 and 0.5 mm. Each voxel value corresponds to an unsigned 8-bit integer tissue label.¹⁴ The anatomical NBPs were rotated to be consistent with a prone position during a 3D OAT scan and cropped to exclude the chest muscle region.

3.1.2.

Realizations of blood vasculature

A specific innovation in our adaptation is the introduction of realistic blood vasculature, as shown in Figs. 2(a) and 2(d). The VICTRE tool iteratively generates sibling and child branches of arterial and venous blood vessel trees within the breast volume, using randomly sampled values based on blood vessel parameters (Table 4) in the configuration file.²⁰ However, the tool assumes four arterial and five venous blood vessel trees with their entry locations predefined by fixed polar angles $\theta$ and a fixed distance from the breast edge $d_{\text{edge}}=20\mathrm{mm}$ [Fig. 2(e)]. This is inconsistent with the anatomy of the breast. Furthermore, the values were hardcoded within the C++ source code of the VICTRE tool, making it inaccessible for users to adjust them through the configuration file. Thus, to more realistically model breast anatomy,⁴¹ the source code of the VICTRE tool was modified to incorporate five sets of blood vessel trees: internal mammary, thoraco-acromial, lateral thoracic, subscapular and thoraco-dorsal, and intercostal arteries and veins [Fig. 2(d)]. Two tunable parameters, vesselEdgeSep1 and vesselEdgeSep2 [$d_{\text{edge}}$ in Table 3], were added to the configuration file to set the distances from the breast edge (1) to the internal mammary, thoraco-acromial, lateral thoracic, and subscapular and thoraco-dorsal arteries and veins and (2) to the intercostal arteries and veins, respectively [Fig. 2(d)]. Table 3 summarizes the entry locations of the blood vessel trees in the left breast, and those in the right breast were set to be bilaterally symmetric.

Fig. 2

Blood vessels in an NBP (type B, left breast) with (a and d) and without (b and e) blood vasculature customization and (c) a clinical OAT image acquired by TomoWave Laboratories employing LOUISA-3D³ at the MD Anderson Cancer Center and postprocessed to extract blood vascular structures.³³ Paraview⁴⁰ was used for volume rendering.

Table 3

Entry locations of blood vessel trees.

Blood vessel	θ (rad), artery	θ (rad), vein	dedge (mm)
Internal mammary	$11\pi /15$, $-8\pi /9$	$41\pi /45$, $-32\pi /45$	2
Thoraco-acromial	$2\pi /5$	$17\pi /30$
Lateral thoracic	$-2\pi /45$	$\pi /5$
Subscapular and thoraco-dorsal	$-\pi /2$	$-13\pi /45$
Intercostal	$16\pi /45$, $-\pi /12$	$\pi /6$, $-3\pi /10$	24

θ: angle in the polar coordinate system.dedge: distance from the edge of the breast along the radial direction.

Among blood vessel tree, branch, and segment parameters in the VICTRE tool, the parameters maxBranch, initRad, minRadFrac, numTry, maxTry, and absMaxTry provided in Table 4 were tuned based on the breast volume percentage occupied by blood vessels estimated from four clinical OAT datasets. The clinical datasets used as references were acquired by TomoWave Laboratories (Houston, Texas, United States) using LOUISA-3D³ at the MD Anderson Cancer Center. The experimental 3D OAT images were reconstructed using a filtered back-projection (FBP) method.⁴² To compute the breast volume percentages occupied by blood vessels from the reconstructed images, the distributions of the total hemoglobin concentration were estimated by employing spectral linear unmixing with optical fluence normalization.³³ Then, a vessel enhancement filter⁴³ was applied to the estimated distributions of the total hemoglobin concentration [Fig. 2(c)], and the ratio of the numbers of the blood vessel voxels and breast voxels was calculated.

Table 4

Blood vessel parameters.

Parameter	Description	Value
maxBranch	Target number of branches	500
initRad	Initial radius of tree (mm)	0.6
minRadFrac	Minimum starting radius as a fraction of parent end radius	0.85
numTry	Number of trial segments to generate	50
maxTry	Maximum number of segments to generate before reducing length	100
absMaxTry	Total number of segment tries	1000

The blood vessels within the depth of 10 mm under the skin layer are dominantly observed with relatively high image contrast in clinical OAT images, i.e., estimated distributions of the initial pressure, but are missing in the VICTRE NBPs [Figs. 2(b) and 2(e)]. To address this, additional blood vasculature [Figs. 2(a) and 2(d)] was stochastically generated using computationally efficient methods (random sampling, Gaussian blurring, Otsu’s thresholding,⁴⁴ and skeletonization), which are commonly used in image processing. The diameter of the blood vessels and their depth from the skin are tunable and were set to 0.75 and 0.375 mm, respectively. Multiple blood vessel segments were formed and alternately divided into arteries and veins depending on the segment locations. Instead of the proposed method, another user-defined method can be employed in our proposed framework if desired. Details of the blood vasculature generation are provided in Appendix A.

3.1.3.

Generation and insertion of lesions

Most common types of malignant lesions, such as invasive ductal carcinoma and invasive lobular carcinoma, have anatomically irregular outlines and develop in milk ducts and lobules.³² Depending on the type of lesion and its aggressiveness, the lesion growth can be accompanied by severe hypoxia in the center region, tumor necrosis, and tumor angiogenesis in the periphery of the viable tumor cells.^45–47 Users can model lesions as being composed solely of a viable tumor cell (VTC) region or containing a necrotic core, VTC region, and peripheral angiogenesis (PA) region [Fig. 3(a)].^47–49 The irregular spiculated boundary of the VTC region was created employing the VICTRE tool. The necrotic core region and PA region were formed via erosion and dilation operations⁵⁰ applied to the surface of the VTC region using spherical structuring elements, i.e., matrices that identify the voxel and its neighboring voxels within a spherical area to be eroded or dilated, respectively. Here, the radii of the spherical structuring elements were 0.75 mm (the thickness of the ring-shaped VTC region)⁴⁹ for the erosion and 5 mm⁴⁸ for the dilation. The generated anatomical NLPs can be optionally inserted into the anatomical NBP. Among the candidate lesion locations predicted via the VICTRE tools, those that do not overlap with other (already inserted) lesions or the skin, nipple, and muscle were randomly chosen as the sites at which to insert the NLP.

Fig. 3

Malignant lesion model: (a) anatomical NLPs without (top) and with a necrotic core and a peripheral angiogenesis region (bottom), and distributions of (b) oxygen saturation $s$ and (c) blood volume fraction $f_b$. The two lesions were inserted at physiologically plausible locations randomly selected among the candidate sites produced by the VICTRE tools. In panels (a)–(c), halves of the lesion volumes are presented to show their cross-sections. In panels (b) and (c), the partial breast volumes clipped at the $y$-coordinate at which both lesions are exhibited are illustrated. The arrows in panel (b) indicate the simulated tumor hypoxia and those in panel (c) indicate the simulated tumor angiogenesis, necrotic tumor core, and relatively high total hemoglobin concentration of the viable tumor cells compared with healthy tissues. These are from a type A breast. Paraview⁴⁰ was used for volume rendering, and color maps were adjusted for better visibility.

4.

Examples of Generated NBPs and Corresponding OAT Images

To provide insight into the visual and quantitative characteristics of NBPs created by the proposed framework and the corresponding OAT images, example NBPs were produced as explained in Sec. 4.1. Subsequently, OAT measurement data were simulated based on a target imaging system described in Sec. 4.2, and 3D estimates of the induced pressure distributions were reconstructed using an FBP method.

4.1.

Examples of Generated NBPs

A total of 84 NBPs of various sizes, 21 for each breast type [(A) breast is almost entirely fatty, (B) breast has scattered areas of fibroglandular density, (C) breast is heterogeneously dense, and (D) breast is extremely dense] were created with a voxel size of 0.125 mm. Of these NBPs, 44 have shapes compatible with the patient’s position during a 3D OAT scan and the others have hemispherical shapes compatible to the use of breast cups to stabilize the breast during imaging.^3,6 Three target wavelengths of 757, 800 (the isosbestic point of deoxy- and oxy-hemoglobin), and 850 nm were selected. Four lesions of different sizes, each smaller than 10 mm in diameter and composed solely of a VTC region, were inserted into each of 80 NBPs (40 natural-shaped and 40 hemispherical-shaped NBPs). Two lesions, one with a necrotic core and PA region and the other without, were inserted into each of the remaining four natural-shaped NBPs [Fig. 4].

Fig. 4

Distributions of functional, optical, and acoustic properties: (a) oxygen saturation $s$, (b) blood volume fraction $f_b$, (c) optical absorption coefficients ${\mu }_a$, and (d) optical scattering coefficients ${\mu }_s$ at a wavelength of 800 nm, (e) sound speed $c$, and (f) density $\rho$ of type A to D breasts (top to bottom). The type A and B breasts (top two rows) have two lesions inserted, one with a necrotic core and PA region and the other without, whereas the type C and D breasts (bottom two rows) have four lesions of different sizes inserted, each smaller than 10 mm in diameter and composed solely of a VTC region. Insets of (a)–(c) show cross-sections, with the lesion location indicated with arrows. Paraview⁴⁰ was used for volume rendering, and color maps were adjusted for better visibility.

The generation of anatomical structures of the breast using the VICTRE tool with the modifications in terms of the blood vessel trees described in Sec. 3.1.2 took approximately 100 to 310 mins using a 16-core Intel Xeon Gold 6130 CPU and 256 GB of memory. The computation time varied depending on the volume of the phantoms. Note that breast sizes, shapes, and structures varied among NBPs [see Fig. 4]. The computation time to generate and insert the blood vasculature under the skin layer and the lesions was 1 to 6 mins (depending on the phantom volume), and the assignment of functional, optical, and acoustic properties took approximately 40, 15, and 2 mins, respectively, using the same machine.

Figure 5(a) presents box plots of the breast volume percentages occupied by blood vessels in which 40 proposed NBPs and 40 unmodified VICTRE NBPs are compared. Figure 5(b) shows box plots of the spatially averaged effective optical attenuation ${\mu }_{\mathrm{eff}}$ estimated from 40 proposed NBPs at the three wavelengths via the Beer–Lambert law-based estimation method in Ref. 33.

Fig. 5

Box plots of (a) breast volume percentages occupied by blood vessels in which 40 proposed NBPs (left) and 40 unmodified VICTRE NBPs (right) are compared and (b) spatially averaged effective optical attenuation coefficients ${\mu }_{\mathrm{eff}}$ estimated from 40 proposed NBPs at wavelengths of 757, 800, and 850 nm (left to right). The used NBPs have hemispherical shapes with radii stochastically sampled from the probability distribution ($a_{1t}$) in Table 2. In panel (a), the black dotted horizontal line indicates the estimated mean of four clinical OAT images acquired by TomoWave Laboratories (Houston, Texas, United States) using LOUISA-3D³ at the MD Anderson Cancer Center. In panel (b), the reference values were calculated from the measurements in Refs. 51 and 8586.–87.

As shown in Fig. 5(a), the reference value (0.439%) of the breast volume percentage occupied by blood vessels, i.e., an estimated mean of the clinical OAT images (a dotted horizontal line), is within the interquartile range of the proposed NBPs (0.336% to 0.48%, median of 0.438%), but it is not within the interquartile range of the VICTRE NBPs (0.249% to 0.39%, median of 0.31%). In Fig. 5(b), the ${\mu }_{\mathrm{eff}}$ estimates of female breasts reported in Ref. 51 and 8586.–87 all fall within the interquartile range of the proposed NBPs for light excitation wavelengths of 800 and 850 nm. For a light wavelength of 757 nm, the two ${\mu }_{\mathrm{eff}}$ estimates from Ref. 51 are slightly out of the interquartile range of the proposed NBPs (the 84th and 85th percentiles). Once more data on the functional and optical properties of breasts, specifically the $f_w$ and ${\mu }_a$ values of breast fat and glandular tissues for different breast types, become available and are incorporated into NBPs, the ${\mu }_{\mathrm{eff}}$ values are expected to better match the reference values found in the existing literature.

4.2.

Example OAT Images from Simulated Data

OAT images that depict the induced pressure distributions were reconstructed from simulated measurement data corresponding to the considered NBPs. The light delivery subsystem of the virtual OAT imaging system consisted of 20 arc-shaped illuminators uniformly distributed along the azimuthal angle. Each illuminator arc had a radius of 145 mm and was modeled using 25 cone beam sources as illustrated in Fig. 6. The acoustic detection subsystem consisted of 51,472 transducer elements uniformly distributed on a hemispherical aperture with a scanning radius of 85 mm. Each virtual transducer element recorded 3,720 time samples of pressure data at a sampling frequency of 20 MHz. Additional details on the virtual OAT imaging system are presented in Appendix C.

Fig. 6

Virtual light delivery system: (a) five linear fiber-optic segments (red lines) are attached on the surface of an arc-shaped illuminator (radius of 145 mm, central angle of 80 deg), and (b) five cone beams (half-angle of 12.5 deg) are emitted from the locations uniformly distributed on each linear fiber-optic segment. Panels (a) and (b) show light delivery at an illumination view. The total number of illuminator arcs is 20, so 500 cone beams illuminate the breast.

OAT data were virtually acquired via multiphysics simulation of the photoacoustic effect and subsequent wave propagation. First, the simulation of photon transport in biological tissues was conducted using the GPU-accelerated MCX version 1.9.0^23,24 to compute the induced initial pressure distribution at the three optical wavelengths of interest (757, 800, and 850 nm). Figures 7(a) and 7(b) illustrate the simulated distributions of optical fluence $\varphi$ and true initial pressure $p_0$, and Figs. 7(c) and 7(d) show the $p_0$ ratios of the breast blood vasculature and lesions between wavelength pairs of 800 and 757 nm and of 850 and 800 nm, respectively.

Fig. 7

Simulated distributions of optical fluence $\varphi$ and initial pressure $p_0$: (a) $\varphi$ and (b) $p_0$ of the breast in the water-filled imaging bowl at a wavelength of 800 nm, and the $p_0$ ratio of the breast blood vasculature and lesions between wavelengths (c) of 800 and 757 nm and (d) of 850 and 800 nm. These are from a type C breast. Paraview⁴⁰ was used for volume rendering, and color maps were adjusted for better visibility.

Next, wave propagation in acoustically heterogeneous media was simulated using the GPU-accelerated k-Wave toolbox version 1.3,²⁵ and pressure traces measured by the virtual transducer elements were recorded. Finally, measurement noise was modeled as independent and identically distributed Gaussian noise with zero mean and a standard deviation equal to 1% of the maximum optoacoustic signal strength of the entire ensemble for all three wavelengths, as was empirically determined based on the in vivo breast OAT data.^3,33 Using the simulated noisy acoustic measurements, the images, i.e., $p_0$ estimates, were reconstructed employing the FBP method with a voxel size of 0.25 mm. Figure 8 presents a visual comparison of 3D OAT images reconstructed from clinical measurement data and simulated measurement data produced using the proposed NBP and the unmodified VICTRE NBP. The OAT image generated using the proposed NBPs [Fig. 8(b)] plausibly resembles the 3D clinical OAT image [Fig. 8(a)], compared with the unmodified VICTRE NBP [Fig. 8(c)].

Fig. 8

Visual comparison of 3D OAT breast images reconstructed from (a) clinical measurement data and simulated measurement data produced using (b) the proposed NBP and (c) the unmodified VICTRE NBP. In panel (a), the clinical data were acquired by TomoWave Laboratories using LOUISA-3D.³ In the unmodified VICTRE NBP, tissue properties were assigned to each tissue type in a piecewise constant manner. All three images were reconstructed using the FBP method.⁴² The images were visualized using Paraview’s volume rendering technique, which accumulates intensities based on the selected color and opacity maps.⁴⁰

5.

Case Study: Acoustic Heterogeneity in Image Reconstruction

To demonstrate the usefulness of the proposed framework in virtual OAT studies, a case study was conducted to explore the impact of acoustic heterogeneity in OAT images of the breast. One phantom for each breast type (A, B, C, or D), four in total, was used for this case study. Acoustic measurements were simulated via virtual OAT data acquisition based on the target imaging system described in Sec. 4.2.

To explore the impact of acoustic heterogeneity in OAT images of the breast, the images, i.e., $p_0$ estimates, were reconstructed from the noisy acoustic measurements simulated in Sec. 4.2 with a voxel size of 0.25 mm. Three different approaches were used, each with the following assumptions: homogeneous acoustic properties of water in Table 7 for the entire computation domain (Approach 1); two sound speed values for water ($c_{\text{water}}$) and breast tissue ($c_{\text{breast}}$) each, and homogeneous density and acoustic attenuation properties of water for the whole domain (Approach 2); and the true distributions of acoustic properties (Approach 3). For image reconstruction, the LSQR⁸⁸ algorithm was implemented⁸⁹ to compute a least-squares estimate of $p_0$. As the stopping rule, the iteration was terminated when the squared error between the true $p_0$ and estimated $p_0$ started to increase to avoid overfitting the noisy data. The goal of this case study was not to illustrate a practical image reconstruction algorithm but to investigate the effect of acoustic heterogeneity on image quality. Thus, a stylized stopping criterion, which requires knowledge of the true object, was purposely chosen to isolate the effect of acoustic heterogeneity from other sources of error or bias in OAT images, such as the design of an appropriate regularization functional. Figure 9 shows the true $p_0$ (first column) and the $p_0$ estimates (second to fourth columns) of Approaches 1 to 3, respectively. To tune $c_{\text{breast}}$ in Approach 2, sound speed values between 1.45 and $1.54\mathrm{mm}/\mu \mathrm{s}$ were swept with a step size of $0.005\mathrm{mm}/\mu \mathrm{s}$, and the selected values for each breast type are presented in Table 8.

Fig. 9

True distributions of initial pressure $p_0$ at a wavelength of 800 nm (first column) and $p_0$ estimates reconstructed assuming homogeneous acoustic properties of water (second column), two sound speeds for water and breast tissue regions each and homogeneous density and acoustic attenuation properties of water (third column), and true distributions of acoustic properties (fourth column). The whole breast regions are presented from a side view in panels (a) and (c), and the regions at depths between 10 and 20 mm from the breast surface are illustrated from a front view in panels (b) and (d). The images in panels (a) and (b) are from a type A breast, and those in panels (c) and (d) are from a type D breast. Paraview⁴⁰ was used for volume rendering, and color maps were adjusted for better visibility.

Table 8

Average relative squared error (RSE) and structural similarity index (SSIM).

The relative squared error $\mathrm{RSE}={\Vert x_{\text{true}}-x_{\mathrm{est}}\Vert }_2^2/{\Vert x_{\text{true}}\Vert }_2^2$ and structural similarity index measure (SSIM)⁹⁰ between the true $p_0$ ($x_{\text{true}}$) and the reconstructed image ($x_{\mathrm{est}}$) were calculated to quantify the accuracy of the reconstructed images. Table 8 provides the average RSE and SSIM calculated from the data of breast types A to D for all three wavelengths. The RSE and SSIM improved for all breast types when acoustic heterogeneity was accounted for in the image reconstruction. As shown in Fig. 9, the mitigation of clutter artifacts using the two-sound speed model (Approach 2) with respect to the homogeneous model (Approach 1) was more significant in the type A breast (almost entirely fatty breast) than the type D breast (extremely dense breast). Here, clutter refers to background artifacts caused by uncompensated scattering or refraction effects. For the fatty type breast, the sound speed mismatch between fatty tissue and water is the main cause of clutter; therefore, the two-sound speed model (Approach 2) can effectively mitigate such artifacts. However, for denser type breasts, glandular tissue has a sound speed similar to water, and acoustic impedance heterogeneity across different breast tissues is the primary cause of clutter; thus the two-sound speed model is less effective.

6.

Conclusion

In this work, a framework to generate realistic 3D NBPs that can be used in large-scale VITs of OAT breast imaging was established. For the first time, anatomically, physiologically, and optoacoustically realistic 3D NBPs and NLPs in varying sizes and shapes could be produced using the proposed framework. This was achieved by extending VICTRE tools to OAT imaging with substantial modifications. Because the proposed framework was implemented in a modular form, it facilitates further customization of the NBPs. Users can replace the relatively simple lesion model applied in this study with more biologically realistic models of tumor growth and metabolism. For example, heterogeneous distributions of multiple necrotic regions and blood vessels proliferating toward the lesion could be included. Future studies may also include explicit geometric modeling of vascular growth associated with tumors, functional and optical modeling of a specific type of breast cancer, analysis of lesion detectability as a function of their size and depth, and mechanical deformations of the breast to simulate slight breast compression in the existing OAT system. In summary, the produced NBPs and NLPs enhance the authenticity of virtual OAT studies, and the proposed framework can be widely employed for the investigation and development of advanced image reconstruction methods and machine learning-based methods, as well as the objective evaluation and optimization of the OAT breast imaging systems. Code packages of the proposed tools and 84 datasets that were made publicly available will enable researchers to immediately benefit from this work.

7.

Appendix A: Generation of Blood Vasculature under the Skin Layer

Blood vasculature under the skin layer [Figs. 2(a) and 2(d)] was stochastically generated visually similar to the blood vessels observed in the clinical OAT images^3–6,33 via relatively computationally inexpensive image processing methods (random sampling, Gaussian blurring, Otsu’s thresholding,⁴⁴ and skeletonization) that are commonly used.

In Algorithm 1, the inputs $N_x$, $N_y$, and $N_z$ are the sizes of the anatomical NBP volume along the $x$-, $y$-, and $z$-axes; inputs ${\sigma }_1$ and ${\sigma }_2$ of Gaussian filters control the distance between blood vessel branches in the vessel tree; and input ${\sigma }_3$ determines the thickness of blood vessels under the skin layer. Based on the breast volume percentage occupied by blood vessels estimated from the clinical OAT images, the values of ${\sigma }_1$, ${\sigma }_2$, and ${\sigma }_3$ were set to 3.875, 6.125, and 0.219 mm (corresponding to the vessel diameter of 0.75 mm), respectively. The erosion in step (7) of Algorithm 1 shrinks bright regions in the breast mask without skin, i.e., voxels under the breast skin. Thus, the degree of erosion determines the distance between the inserted blood vasculature and the skin. For example, three erosion iterations at a voxel size of 0.125 mm yield the blood vessels at a depth of 0.375 mm under the skin. The proposed method does not produce a specific number and volume of arteries and veins under the skin layer. Instead, it stochastically creates visually plausible vascular structures employing computationally efficient methods, so the number and volume of the blood vessel segments produced in step (11) of Algorithm 1 are not always even for arteries and veins. The blood vessel segments were assigned to the artery mask and vein mask alternately, depending on the segment locations, in step (12) of Algorithm 1.

Algorithm 1

Blood vasculature under the skin layer

The produced blood vasculature was embedded into the anatomical NBPs by assigning the artery label to the artery mask region and the vein label to the vein mask region in the anatomical NBPs [Figs. 2(a) and 2(d)].

8.

Appendix B: Modeling Distributions of Oxygen Saturation and Blood Volume Fraction

To model a spatially varying distribution of oxygen saturation $s(\mathbf{r})$ at each spatial coordinate $\mathbf{r}$ within the phantom, the following diffusion-reaction PDE⁸¹ was solved:

Here, $\mu >0$ is a numerical parameter controlling the smoothness of the oxygen saturation $s(\mathbf{r})$; $\mathrm{\Delta }$ is the Laplace operator; ${\mathcal{T}}_s=\{$ skin, artery, vein, lesion $\}$ is a set of tissues of interest; ${\chi }_i(\mathbf{r})$ and ${\overline{s}}_i$ are the indicator function and the target value of oxygen saturation, respectively, (randomly sampled according to the tissue-specific probability distributions in Table 5) assigned to the tissue $i$ ($i\in {\mathcal{T}}_s$).

Tumor angiogenesis was modeled in the blood volume fraction distribution $f_b(\mathbf{r})$ at a spatial coordinate $\mathbf{r}$ within the phantom by solving the following diffusion-reaction PDE:

The PDEs above were solved using a linear finite element on an unstructured tetrahedral mesh that is fitted to tissue interfaces. Unstructured meshes were generated using the pygalmesh software,^91,92 and the PDE was solved using the parallel finite element library FEniCS.⁸²

9.

Appendix C: Details of Virtual OAT Imaging System and Data Acquisition

The light delivery mechanism and the measurement geometry of the target virtual OAT imaging system were chosen to emulate one of the existing OAT systems for breast imaging, LOUISA-3D^3,33 developed by TomoWave Laboratories (Houston, Texas, United States). The detailed illumination geometry of the virtual light delivery system using 20 arc-shaped illuminators [see Fig. 6] is described in Table 9.

Table 9

Light delivery geometry of arc-shaped illuminators.

For the virtual OAT data acquisition, idealized point-like transducers, uniformly distributed on the arc-shaped optoacoustic probe specified in Table 10, were assumed; the probe collects a total of 3,720 time samples of pressure data at a sampling frequency of 20 MHz. For the simulation of OAT data acquisition, the k-Wave GPU code was used as the acoustic wave propagation simulator. With the assumption of idealized point-like transducers, the possible option to define transducer locations in the simulation supported by the k-Wave GPU code is a binary matrix, where the value at the corresponding transducer locations is “1” and at the others is “0.” Therefore, 51,472 transducer locations were defined as a target measurement geometry, excluding the overlaps due to the discretization of Cartesian coordinates into the binary matrix with the given voxel size.

Table 10

Measurement geometry of the arc-shaped optoacoustic probe.

Acoustic OAT measurements were virtually acquired through the simulation of photon transport in breast tissues, calculation of initial pressure distribution $p_0$, and simulation of acoustic wave propagation. In the simulation of the photon transport using the GPU-accelerated MCX, the light source and direction parameters were set according to the target system design [Fig. 6 and Table 9], and the simulation domain size was set to $340\times 340\times 170\text{voxels}$ with a spatial step size of 0.5 mm. Among the optical NBPs employed as inputs of the simulation, the distributions ${\mu }_a$ and ${\mu }_s$ for wavelengths of 757, 800, and 850 nm were downsampled using linear interpolation (from a voxel size of 0.125 to 0.5 mm), and constant values $g$ and $n$ were averaged within breast tissues. With these inputs, the distributions of optical fluence $\varphi$ were simulated using ${10}^8$ photons per beam for a time duration of 5 ns.

The true $p_0=\mathrm{\Gamma }{\mu }_a\varphi$ at a voxel size of 0.25 mm [Fig. 7(b)] was calculated via elementwise multiplication of ${\mu }_a$ and the simulated $\varphi$, assuming a constant Grüneisen parameter $\mathrm{\Gamma }=1$, as is commonly done for soft tissues.^16,33 The ${\mu }_a$ and the simulated $\varphi$ were downsampled (from a voxel size of 0.125 to 0.25 mm) and upsampled (from a voxel size of 0.5 to 0.25 mm) via linear interpolation, respectively. The MCX simulation at a coarser grid accompanied by the upsampling not only is equivalent to smoothing the $p_0$ in the acoustic wave propagation simulation to reduce the amplitudes of the high spatial frequency components but also significantly reduces the computation time of the MCX simulation.

In the simulation of acoustic wave propagation using the k-Wave toolbox and its GPU code,²⁵ the measurement geometry was configured based on the target system design [Table 10] as explained above. The voxel size was set to the smallest possible, i.e., 0.25 mm, given the measurement geometry and the use of an NVIDIA Tesla V100 GPU with 32 GB memory. The simulation domain size was set to $700\times 700\times 350\text{voxels}$. The acoustic NBPs $c$, $\rho$, and ${\alpha }_0$ employed to simulate acoustic measurements were downsampled using linear interpolation (from a voxel size of 0.125 to 0.25 mm). An anisotropic absorbing boundary layer, known as a perfectly matched layer, was set to be outside the simulation domain to prevent undesired acoustic reflection at the boundaries. A pressure wavefield generated by the simulated $p_0$ was virtually propagated based on the acoustic NBPs and measured by the virtual transducer elements in the target imaging system.