1.

Introduction

Optoacoustic tomography (OAT), also known as photoacoustic computed tomography (PACT), is an emerging imaging modality that shows promise in sensitivity for breast cancer detection, especially in dense breasts.^1–4 OAT images exhibit greater detection sensitivity for highly vascularized, i.e., aggressive, breast tumors,⁵ and greater diagnostic specificity for all tumors over other modalities, such as x-ray mammography and ultrasound.^1,5,6 Another advantage of OAT over x-ray mammography is that OAT does not involve ionizing radiation.^1,2 Due to advances in OAT systems design and image reconstruction, a three-dimensional (3D) volumetric scan of the entire breast is now possible.^3,5,7 At a single near-infrared illumination wavelength, natural chromophores in the breast tissue, such as hemoglobin, act as endogenous OAT contrast agents. From external ultrasound measurements of the pressure induced by the laser pulses, the spatial distribution of the chromophores can be estimated; this provides a qualitative, anatomical measure of the blood vasculature.^3,5,7–10 Measurements from multiple illumination wavelengths matching the local maximum, minimum, and isosbestic point of deoxy- and oxyhemoglobin can be reconstructed into quantitative estimates of the blood oxygen saturation.^5,11–16 This technology referred to as quantitative OAT (qOAT), also known as quantitative PACT, when combined with ultrasound, provides both anatomical and functional information of the breast that can facilitate detecting tumor angiogenesis and hypoxia.¹

In 3D OAT breast imaging, it is not feasible to deliver the optical fluence uniformly throughout the whole breast volume, due to optical attenuation in tissue and design constraints of the imaging systems.^{5–7,9,13,17} Whereas distributions of the optical absorption coefficient are determined by tissue types and physiological status, initial pressure distributions decay with depth in the tissue because of light attenuation. Furthermore, the attenuation is wavelength-dependent. Therefore, direct application of spectral linear unmixing methods to reconstructed OAT images, which correspond to the estimated initial pressure distributions, results in inaccurate estimates of blood oxygen saturation.^14,15,18,19

To improve visualization of the reconstructed volumetric images, the optoacoustic imaging community has utilized a commercial tool, AMIRA (Thermo Fisher Scientific),²⁰ and open source, interactive tools such as ImageJ (Wayne Rasband),²¹ ParaView (Kitware),²² and 3D Slicer (Kitware).²³ Also, a 3D PHOVIS (POSTECH, Korea)²⁴ has been recently released that is developed specifically for optoacoustic imaging. However, these tools do not provide physics-informed image processing for contrast enhancement at depth in reconstructed OAT volumes. In addition, these visualization tools are semi-automatic and require substantial manual intervention by the user. Pattyn et al.²⁵ proposed a model-based method to compensate for the optical fluence distribution within a heterogeneous physical phantom that mimics a breast. Monte Carlo (MC) simulation was employed, and known optical properties of the phantom were assumed. However, in practice, the distributions of optical properties within the breast are generally unknown.

A straightforward physics-informed image processing method is proposed to compensate for both the nonuniform incident optical fluence at the breast surface and the wavelength- and depth-dependent optical attenuation within the breast, and the impact of the proposed method on accuracy of the linear unmixing of deoxy- and oxyhemoglobin is investigated. The contributions of this paper are twofold. First, this study establishes an implementation of the linear unmixing method for use with a large object such as a female breast. Second, the proposed method improves sensitivity of the 3D OAT breast imaging by improving contrast at depth.

The paper is organized as follows. Background materials, including existing illumination systems in 3D OAT breast imaging and spectral linear unmixing, are provided in Sec. 2. The proposed method is explained in Sec. 3, and the study description and evaluation metrics are provided in Sec. 4. Results from computer-simulation and experimental studies are presented in Sec. 5, and a discussion is given in Sec. 6. The conclusions of the study are provided in Sec. 7.

2.

Background

In OAT imaging, a short laser pulse is employed to irradiate an object at time $t=0$ and conversion of the absorbed optical energy into the thermal energy results in the generation of an initial pressure distribution $p_0(\mathbf{r},\lambda )$, where $\mathbf{r}=(x,y,z)\in {\mathbb{R}}^3$ and $\lambda$ is a wavelength. The pressure distribution subsequently propagates and is measured by multiple ultrasonic transducers located on a measurement aperture ${\mathrm{\Omega }}_0\subset {\mathbb{R}}^3$ that partially or completely surrounds the object. The propagated pressure wavefields, i.e., the optoacoustic signals, at time $t>0$ are denoted as $p(\mathbf{r},t)$.

By solving the associated acoustic inverse problem,²⁶ an estimate of the absorbed energy density distribution within the object can be obtained. Functional quantities such as the vascular blood oxygenation can be reconstructed via qOAT from multiwavelength measurements.^5,11–15,27

In most implementations of 3D OAT breast imaging, the patient lies prone on a bed with their breast located inside a water-filled bowl just below the surface plane of the bed. The breast is illuminated at a near-infrared wavelength with a short laser pulse.^{3,5,7,9,17,28} The induced pressure waves $p(\mathbf{r},t)$ propagate out of the breast and are measured with the ultrawide-band ultrasonic transducers.¹⁶

2.1.

Illumination in 3D OAT Breast Imaging

Several different 3D OAT breast imaging systems have been proposed and established, but their data-acquisition principles are similar.^{3,5,7,9,17,28} Existing systems for 3D OAT breast imaging are equipped with two subsystems: an illumination system and an optoacoustic signal detection system. The focus here is on the illumination system. A common feature in the existing illumination systems is that the light is delivered in a radially symmetric pattern to the breast surface from either single^3,7,9,17,28 or multiple⁵ light sources. The laser fluences involved are well below the maximum permissible exposure for skin defined by the American National Standards Institute (ANSI).²⁹

Table 1 presents salient design features of the light-delivery systems employed by the five OAT breast imaging systems shown in Fig. 1. These systems deliver the light to the entire breast surface. However, less optical energy per unit area (optical fluence) is delivered to regions near the chest wall compared with the center of the breast, and this imbalanced distribution of the incident optical fluence causes lower voxel brightness near the chest wall in the reconstructed OAT images.

Table 1

Light-delivery systems employed by the five OAT breast imaging systems in Fig. 1.

Ref.	λ (nm)	Salient design features
7,9	755, 795	• A conical laser beam is emitted from below and then reflected through a planoconvex lens and a holographic diffuser [Fig. 1(a)].
		• The breast is contained by a 0.5-mm-thick polyethylene terephthalate glycol-modified (PETG) cup which is optically and acoustically transparent.
3	1064	• A donut-shaped laser beam is emitted from below and then reflected through an axicon lens and an engineered diffuser [Fig. 1(b)].
		• An agar pillow is used to slightly compress the breast.
28	755, 1064	• A laser beam is split into a bottom beam and side beams, and then those are diverged via concave lenses [Fig. 1(c)].
		• The side beams are emitted slightly upward from nine optical fiber bundles that rotate around the breast in discrete steps.
17	532	• A ring-shaped light beam is formed via a cone-shaped reflector, stationary conical reflector, and mobile conical reflector [Fig. 1(d)].
		• The mobile conical reflector vertically moves together with the ring-shaped detector array during the scan.
5	755 to 1064	• Laser beams are emitted from five fiber-optic segments that are constrained to the surface of a light paddle that rotates around the breast in discrete steps [Fig. 1(e)].
		• The PETG breast cup is 0.1 mm thick.

Fig. 1

Illumination in 3D OAT imaging systems for the breast: (a) Kruger et al.⁹ and Toi et al.,⁷ (b) Lin et al.,³ (c) Schoustra et al.,²⁸ (d) Alshahrani et al.,¹⁷ and (e) Oraevsky et al.⁵

3.

Normalization of Optical Fluence Distribution

The proposed method seeks to estimate and compensate for the nonuniform optical fluence distribution. This is referred as normalization of the optical fluence distribution hereafter. The method was designed based on a common feature of the existing illumination systems; specifically, that radially symmetric light delivery to the breast surface is employed. In this section, the details of the method are described based on the reference imaging system shown in Fig. 2, where a patient lies prone on the examination bed and the patient’s breast is located inside the breast cup. A spherical coordinate system is assumed, with the origin corresponding to the center of the breast cup [see Figs. 2(a) and 3(d)]. Here, $\theta$ is defined to be a polar angle from the positive $z$ axis with $0\le \theta \le \pi$.

Fig. 2

3D OAT scan using LOUISA-3D: (a) breast scan schematic and (b) photograph of phantom scan. (b) Nonuniform illumination is observed on the surface of the tissue-mimicking physical phantom.

Fig. 3

(a)–(c) 3D OAT breast image ($\widehat{\alpha }$) of a healthy volunteer at a wavelength of 755 nm, scanned by TomoWave Laboratories using LOUISA-3D⁵ at MD Anderson Cancer Center, and (d) breast shapes in a given breast cup. Maximum voxel brightness projection (MVBP) along (a) $x$ axis and cross-sections on (b) $y\text{-}z$ plane at $x=0$ mm and (c) $x\text{-}y$ plane at $z=-30\mathrm{mm}$. The slice is indicated with a white dotted line in (a). The image ($\widehat{\alpha }$) was reconstructed using filtered backprojection (FBP)⁴³ method. The brightness range of the images was adjusted for better visibility. In panels (b) and (c), a white solid circle indicates the location of the brightest voxel in the cross-section, and a cyan dashed line represents the approximated breast boundary.

In the proposed method, the nonuniform distribution of the optical fluence within the breast is estimated from the voxel values in the reconstructed 3D OAT image $\widehat{\mathit{\alpha }}$ that is an estimate of $p_0(\mathbf{r},\lambda )$ discretized employing a uniform Cartesian lattice. The following reasonable assumptions are made:

The normalization of the optical fluence distribution is conducted in the following order: (1) compensation for nonuniform incident optical fluence, (2) estimation of breast surface and depth of each voxel relative to the breast surface, and (3) compensation for the effective optical attenuation. The location of the breast surface must be known for optical attenuation compensation, but breast surface detection is highly challenging because the top skin layer (epidermis) can appear dimmer than the noise surrounding the breast due to the nonuniform incident optical fluence [see Figs. 3(a)–3(c)]. Therefore, the incident optical fluence needs to be normalized in advance of the breast surface detection.

In the proposed method, a hemispherical breast stabilizer cup is assumed, as it is employed in several 3D OAT breast imaging systems.^5,7,9 The cup is selected for each breast size, so it maintains the breast shape radially symmetric. Therefore, the only possible breast shape is a partial spheroid, as shown in Fig. 3(d). Whereas the nipple and areola absorb more light than the other breast tissues due to their high concentration of pigment, their impact on the optical fluence distribution within the breast is insignificant because of their relatively small volume.⁴⁵ Thus, in the proposed method, the region $160\mathrm{deg}<\theta \le 180\mathrm{deg}$, which was determined by an average diameter ratio of the breast and areola,⁴⁵ is excluded from consideration.

The flowchart of the proposed method is provided in Fig. 4, and the details of each step are given in the following sections.

Fig. 4

Flowchart of normalization of optical fluence distribution in 3D OAT breast images.

3.1.

Compensation for Nonuniform Incident Optical Fluence

Under assumption A1, the distribution of the incident optical fluence is radially symmetric, which means it can be interpreted as a function of polar angle [Fig. 3(d)]. If strong optical absorbers that have similar ${\mu }_a$ values, such as blood vessels, are densely located near the object surface ($d=0$), the spatial distribution of $p_0$ on the surface is proportional to the distribution of the incident optical fluence ($p_0\propto \varphi$). As mentioned earlier, the blood vessels in subdermal regions appear as the brightest voxels at any polar angle in the reconstructed image $\widehat{\mathbf{\alpha }}$ [assumption A3; see Figs. 3(b) and 3(c)]. Thus, the distribution of the incident optical fluence can be approximated by the voxel brightness of the blood vessels near the breast skin layer (subdermal), i.e., maximum voxel brightness projection (MVBP) at discretized polar angles $[\mathbf{\theta }{]}_n\in ({\theta }_i-\frac{\mathrm{\Delta }\theta }{2},{\theta }_i+\frac{\mathrm{\Delta }\theta }{2})$ with an increment $\mathrm{\Delta }\theta$ of 1 deg. The polar angle $[\mathbf{\theta }{]}_n$ is calculated as ${\cos }^{-1}(z_n/r_n)$, where $r_n=\sqrt{x_n^2+y_n^2+z_n^2}$ is the distance of the n’th voxel from the origin, and $x_n$, $y_n$, and $z_n$ are the $x$, $y$, and $z$-coordinates of the n’th voxel in the uniform Cartesian grid, respectively.

An $L$-degree polynomial curve $q_L(\theta )$ is fitted to the extracted maximum voxel brightness according to the discretized polar angles within the range of interest (90 deg to 160 deg), where the region containing the nipple and areola is excluded, as shown in Fig. 5(a). The estimated distribution of the incident optical fluence ${\widehat{\mathbf{\varphi }}}_{\mathbf{0}}$ is described as

Eq. (2)

$${[{\widehat{\mathbf{\varphi }}}_{\mathbf{0}}]}_n=q_L({[\mathbf{\theta }]}_n),n=1,\dots N,$$

where $N$ is the total number of voxels. The $L$ is set depending on the illumination pattern. It was found that $L=1$ and $L=2$ were sufficient for accurately fitting the data in the experimental studies (Sec. 4.2) and computer-simulation studies (Sec. 4.1), respectively.

Fig. 5

Compensation for nonuniform incident optical fluence: (a) estimated incident optical fluence as a function of polar angle and (b) MVBP of 3D OAT breast image after the compensation (${\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}$) along $x$ axis. The results were obtained from Fig. 3(a). In panel (a), a black solid line indicates the maximum voxel brightness according to polar angles, and a red solid line represents a first-degree polynomial curve $q_1({[\mathbf{\theta }]}_i)$ fitted to the maximum voxel brightness. The polar angles to the right of the blue dashed line correspond with nipple and areola in (a).

Elementwise division of ${\widehat{\mathbf{\varphi }}}_{\mathbf{0}}$ is applied to the reconstructed image $\widehat{\mathbf{\alpha }}$, to compensate for the nonuniform incident optical fluence as shown in Fig. 5(b):

Eq. (3)

$${[{\widehat{\mathbf{\alpha }}}_{\mathit{N}\mathbf{0}}]}_n=\frac{{[\widehat{\mathbf{\alpha }}]}_n}{{[{\widehat{\mathbf{\varphi }}}_{\mathbf{0}}]}_n},n=1,\dots ,N.$$

A comparison of the images before and after the compensation is shown in Fig. 6

Fig. 6

Comparison of images before and after compensation for non-uniform incident optical fluence: MVBP of a 5-mm-thick $y$ slice at $y=0$ mm (a) before ($\widehat{\mathbf{\alpha }}$) and (b) after (${\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}$) the compensation along $y$ axis and (c) MVBP of their difference (${\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}-\widehat{\mathbf{\alpha }}$) along $y$ axis. The results were obtained from Fig. 3(a). The voxel brightness near the chest wall ($\theta =90\mathrm{deg}$) in (a) is lower than in the region near the areola ($\theta \ge 160\mathrm{deg}$) and, accordingly, the compensation procedure leads to a higher gain near the chest wall as shown in (b) and (c).

3.2.

Estimation of Breast Surface and Depth

Under assumption A3, the blood vessels located in subdermal regions can be employed to infer the breast surface in the reconstructed 3D OAT image $\widehat{\mathbf{\alpha }}$. The average range of skin thickness of healthy female human breasts is between 0.5 and 2.4 mm.⁴⁶ This is a relatively thin layer that attenuates light negligibly in comparison to attenuation in the bulk. To extract the voxels that belong to blood vessels in the close proximity of the breast surface, first a 3D median filter is applied to reduce the noise. Subsequently, the contrast of the resulting image is increased by elementwise square operation:

Eq. (4)

$${[{\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}^{\prime }]}_n={([\mathrm{med}\{{\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}\}{]}_n)}^2,n=1,\dots ,N,$$

where $\mathrm{med}\{·\}$ is a 3D median filter function with a 3-by-3-by-3 filter.

The voxels corresponding to the blood vessels near the breast surface are extracted using Otsu thresholding⁴⁷ applied to ${\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}^{\prime }$. The set of the extracted voxels is defined as $V=\{n\le N:{[{\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}^{\prime }]}_n\ge T\}$ where $T$ is an Otsu’s threshold. For each polar angle, the voxels in $V$ that are the farthest from the $z$ axis are selected to estimate the breast boundary. The estimated radius of the breast in the cross-section of slice ($j$’th $z$-slice), is given as

Eq. (5)

$$[{\mathbf{\rho }}_{\mathit{z}}{]}_j=\max \{[\mathbf{\rho }{]}_n\}n\in S_j\cap V,$$

where $[\mathbf{\rho }{]}_n$ is the distance of the $n$’th voxel from the $z$-axis calculated as $\sqrt{x_n^2+y_n^2}$, and $S_j$ denotes the set of all voxels, on the $j$’th $z$-slice, i.e., the $x\text{-}y$ plane.

In Fig. 7(a), blue circle markers correspond to an estimate ${\mathbf{\rho }}_{\mathit{z}}$ for each $z$-slice. An elliptical curve is fit to ${\mathbf{\rho }}_{\mathit{z}}$ to obtain a smooth representation of the breast boundary according to assumption A4. The ellipse equation is $\frac{(z-z_C{)}^2}{a^2}+\frac{{\rho }^2}{b^2}=1$. Here, $a$ and $b$ are lengths along semimajor and semiminor axes, respectively, $z_C$ is a $z$-coordinate of the ellipse center. These parameters are determined by the elliptical curve fitting. The surface formed by rotation of the estimated elliptical curve [a red solid line in Fig. 7(a)] around the $z$ axis is then used to approximate the breast surface as shown in Fig. 7(b). From this, a breast mask $M_{\text{Breast}}$ is built by assigning the value “1” to voxels inside the surface and “0” outside.

Fig. 7

Breast surface estimation: (a) estimated radii on $x\mathit{\text{-y}}$ planes ${\widehat{\mathbf{\rho }}}_{\mathbf{z}}$; (b) estimated breast surface; and (c) estimated breast boundary on $y\text{-}z$ plane at $x=0$ overlaid on the MVBP of ${\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}$ along the $x$ axis.

Finally, the depth $\mathbf{d}$ at each voxel that is used in the Beer–Lambert law is approximated as the minimum distance from the breast surface, i.e., the estimated spheroid surface, using Newton’s method.⁴⁸

3.3.

Compensation for the Effective Optical Attenuation

With consideration of assumptions A2 and A5, optical attenuation can be approximated as a function of depth from the decay of the voxel brightness inside the blood vessels in the reconstructed 3D OAT image. Specifically, such an approximation uses the maximum voxel brightness value in a certain depth range of $(d_m-\frac{{\mathrm{\Delta }}_d}{2},d_m+\frac{{\mathrm{\Delta }}_d}{2})$, where ${\mathrm{\Delta }}_d$ is an increment of 1 voxel:

Eq. (6)

$${[{\widehat{\mathbf{\alpha }}}_{\mathbf{BV}}]}_m={\max }_{n\in C_m}\{{[\widehat{\mathbf{\alpha }}]}_n\}.$$

Here, $C_m$ denotes a set of voxels in the $m$’th depth range. Figure 8 shows ${\widehat{\mathbf{\alpha }}}_{\mathbf{BV}}$ according to depth.

Fig. 8

Estimation of optical attenuation at a wavelength of 755 nm. A black solid line indicates maximum voxel brightness ${\widehat{\mathbf{\alpha }}}_{\mathrm{BV}}$ at a certain depth range of $(d_m-\frac{{\mathrm{\Delta }}_d}{2},d_m+\frac{{\mathrm{\Delta }}_d}{2})$, and a blue curve is the estimated optical attenuation ${\widehat{\mathbf{\varphi }}}_{\mathbf{a}}$. The ${\mu }_{\mathrm{eff}}$ estimated from the 3D OAT breast image was $1.0041{\mathrm{cm}}^{-1}$.

An exponential curve based on the Beer–Lambert law (assumption A5) is fit to ${\widehat{\mathbf{\alpha }}}_{\mathbf{BV}}$ as shown in Fig. 8. The estimated optical attenuation is expressed as

Eq. (7)

$$[{\widehat{\mathbf{f}}}_{\mathit{a}}{]}_n=c^{*}\exp (-{\mu }_{\mathrm{eff}}^{*}{[\mathbf{d}]}_n),n=1,\dots ,N,$$

where $c^{*}$ and ${\mu }_{\mathrm{eff}}^{*}$ are the prefactor and the effective optical attenuation coefficient estimated from the curve fitting, respectively.

Elementwise division of ${\widehat{\mathbf{\varphi }}}_{\mathbf{a}}$ is applied to the image after normalization of nonuniform incident optical fluence ${\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}$ as follows:

Eq. (8)

$$[{\widehat{\mathbf{\alpha }}}_{\mathbf{N}}{]}_n=\{\begin{array}{ll}\frac{{[{\widehat{\mathbf{\alpha }}}_{\mathbf{N}\mathbf{0}}]}_n}{{[{\widehat{\mathbf{\varphi }}}_{\mathbf{a}}]}_n},& \text{if\hspace{0.17em}}n\in M_{\text{Breast}}\\ 0,& \text{otherwise}.\end{array}$$

4.

Study Description and Evaluation Metrics

5.

Results

5.1.

Computer-Simulation Studies Results

5.1.1.

Physical measures of image quality

Figure 9 shows results obtained by applying CLAHE and the proposed optical fluence normalization method to the reconstructed estimate of ${\mathbf{p}}_{\mathbf{0}}$. In the results from CLAHE [Fig. 9(c)] and the proposed method [Fig. 9(d)], more structures at depths deeper than 5 mm (green to red color) were revealed compared to the reconstructed initial pressure distribution in Fig. 9(b). The vasculature in the image produced using the proposed optical fluence normalization in Fig. 9(d) is visually more similar to that in the ground truth [Fig. 9(a)] than that produced by CLAHE in Fig. 9(c).

Fig. 9

Comparison between distributions of (a) the optical absorption coefficient ${\mu }_{\mathbf{a}}$, (b) initial pressure estimate reconstructed from the noisy measurements, simulated at a wavelength of 800 nm, using FBP with no normalization, and (c) images processed via CLAHE and (d) optical fluence normalization method. The images are presented as MVBP along $y$ axis and color-encoded by depth. A depth range of 0 to 30 mm was visualized. A Jet color map in MATLAB was used to illustrate the breast tissues at different depths.

Figure 10 shows PSNR [Eq. (9)] and SSIM [Eq. (10)] comparisons between CLAHE [Fig. 9(c)] and the proposed method [Fig. 9(d)]. As shown in Fig. 10, the results of the proposed method showed higher PSNR and SSIM than those of no normalization and CLAHE for all three wavelengths.

Fig. 10

Comparison on PSNR and SSIM between no normalization (black), CLAHE (cyan), and the proposed method (red).

5.1.2.

Task-based measures of image quality

Figure 11 shows the detected blood vasculature and the estimated blood oxygenation using two- and three-wavelength linear unmixing methods with no optical fluence normalization [Fig. 11(a)], CLAHE [Fig. 11(b)], and the proposed method [Fig. 11(c)].

Fig. 11

Estimates of vascular blood oxygenation obtained using (a) no optical fluence normalization, (b) CLAHE, and (c) the proposed method. The used wavelength pairs are 757 and 800 nm (first column), 757 and 850 nm (second column), 800 and 850 nm (third column), and 757, 800, and 850 nm (fourth column). The vascular blood oxygenation of the numerical phanton (ground truth) is shown on the top right. Paraview²² was used for volume rendering.

With respect to the blood vasculature detection, the majority of the voxels corresponding to the blood vessels were not detected without normalization of the optical fluence [Fig. 11(a)]. Although many of the blood vessel voxels near the breast surface were detected via CLAHE, the voxels in the region deeper than 15 mm [orange to red color in Fig. 9(a)] were not detected [Fig. 11(b)]. The proposed method significantly improved the blood vasculature detectability [Fig. 11(c)]. With respect to estimation of the vascular blood oxygenation, the proposed method [Fig. 11(c)] enhanced the estimation accuracy in regions deeper than 18 mm [orange to red color in Fig. 9(a)] for all choices of wavelengths [Fig. 11(c)]. The results of vascular oxygenation estimation from CLAHE [Fig. 11(b)] were relatively inaccurate regardless of the wavelength pairs and voxel location, compared to the proposed method [Fig. 11(c)].

Table 2 provides detectability (DET), detectability-classification accuracy, and classification accuracy with respect to arteries (TAR and DTAR), veins (TVR and DTVR), and both (ACC and DACC).

Table 2

Artery/vein detectability and classification accuracy (%).

λs (nm)	Normalization	DET	TAR	TVR	ACC	DTAR	DTVR	DACC
757, 800	None	12.00	89.81	75.97	81.58	3.62	17.69	11.99
	CLAHE	33.53	72.15	79.16	76.32	16.71	38.83	29.86
	Proposed method	78.89	77.88	97.21	89.38	66.85	76.69	72.71
757, 850	None	11.39	74.78	88.71	83.06	2.78	17.08	11.32
	CLAHE	31.04	62.54	90.82	79.36	16.21	37.12	28.65
	Proposed method	76.19	71.62	96.52	86.43	58.98	74.38	68.14
800, 850	None	11.87	30.41	85.61	63.24	2.88	15.71	10.51
	CLAHE	34.75	25.81	86.07	61.65	16.39	30.37	24.70
	Proposed method	79.51	48.66	77.47	65.79	37.35	65.33	53.99
757, 800, 850	None	11.58	78.35	88.06	84.12	3.07	17.28	11.52
	CLAHE	32.13	64.91	90.03	79.85	16.96	38.05	29.51
	Proposed method	77.11	73.12	96.80	87.21	61.01	75.24	69.47

Note: For each metric and each choice of wavelengths, the entry corresponding to the best performing method is in bold.

As presented in Table 2, the DET of the proposed method was, on average, 6.66 and 2.37 times greater than no optical fluence normalization and CLAHE, respectively. The proposed method showed slightly better ACC compared with the others in Table 2. The proposed method increased the DACC by 5.81 and 2.34 times on average compared to no optical fluence normalization and CLAHE, respectively. The distribution of artery/vein voxels was not uniform with respect to depth. There were more voxels that correspond to the blood vessels (veins in particular) near the surface. Thus, further analysis of the classification accuracy according to depth will be presented henceforward [Fig. 12(b)].

Fig. 12

(a) Artery/vein detectability and (b) classification accuracy of no optical fluence normalization (black color), CLAHE (cyan color), and the proposed method (red color), according to 10 mm-depth ranges. The used wavelength pairs are 757 and 800 nm (dashed lines); 757 and 850 nm (dash-dotted lines); 800 and 850 nm (dotted lines); and 757, 800, and 850 nm (solid lines).

Figure 12 shows (a) vasculature detectability and (b) artery/vein classification accuracy of no optical normalization, CLAHE, and the proposed method as a function of depth. In this analysis, depth was quantized using 5 bins with a width of 10 mm. At all depth ranges, the proposed method [red color in Fig. 12(a)] outperformed the other two [black and cyan colors in Fig. 12(a)] in blood vessel detectability. The artery/vein voxels located deeper than 20 and 30 mm were not detected when no optical fluence normalization [black color in Fig. 12(a)] and CLAHE [cyan color in Fig. 12(a)] were applied, respectively. As shown in Fig. 12(b), in the depth ranges of 10 to 20 mm and 20 to 30 mm, the ACC of the proposed method (red color) was higher than the others (black and cyan colors), and the ACC largely dropped in the results of all three methods at a depth deeper than 30 mm. It is speculated that this is because the strength of the attenuated optoacoustic signals at depths deeper than 30 mm is similar to or lower than that of the noise. The ACC of CLAHE [cyan color in Fig. 12(b)] was either lower or slightly higher, up to 2.53%, than that of no optical fluence normalization [black color in Fig. 12(b)].

5.1.3.

Results from Experimental Studies

Figure 13 shows reconstructed 3D OAT images with no optical fluence normalization [Figs. 13(a) and 13(b)], CLAHE [Figs. 13(c) and 13(d)], and the proposed method [Figs. 13(e) and 13(f)]. A depth range of 0 to 30 mm was visualized. In Fig. 13, the region deeper than 15 mm (orange to red color) is nearly invisible to the human eye in the CLAHE results [Figs. 13(c) and 13(d)] while it is clearly visible in the results of the proposed method [Figs. 13(e) and 13(f)]. Additional visualization of the comparison in Figs. 13(a), 13(c), and 13(e) is provided in Video 1, showing a $z$-slice ($x\text{-}y$ plane) of the breast image at each descretized $z$ location with an increment of 0.25 mm (from $-46$ to $-2.25\mathrm{mm}$). In Video 1, the visibility of the blood vessels seated deeper than 15 mm (orange to red color) is consistent with the results in Fig. 13. The effective optical attenuation coefficient estimated from the left breast [Figs. 13(a), 13(c), and 13(e)] was $1{\mathrm{cm}}^{-1}$ and that from the right breast [Figs. 13(b), 13(d), and 13(f)] was $0.98{\mathrm{cm}}^{-1}$.

Fig. 13

Comparison between reconstructed images with [(a), (b)] no optical fluence normalization, [(c), (d)] CLAHE, and [(e), (f)] the proposed method. The used wavelength was 755 nm. Images in the left column [(a), (c), and (e)] are from the left breast and those in the right column [(b), (d), and (f)] are from the right breast. The images (a) to (f) are presented in MVBP of the entire breast volume along $y$ axis. The still images in panel (g) of Video 1 (Video 1, MP4, 944 kB [URL: https://doi.org/10.1117/1.JBO.27.3.036001.1]) illustrate a $z$ slice ($x\text{-}y$ plane) at $z=-10\mathrm{mm}$. The images were color-encoded by depth using the Jet color map in MATLAB.

6.

Discussion

In spite of the method’s simplicity, the numerical results demonstrated that the proposed method significantly improved vasculature detectability by compensating for optical attenuation and increased estimation accuracy of the vascular blood oxygenation by mitigating the spectral coloring effect (Figs. 11 and 12, and Table 2). Voxel brightness in the reconstructed estimate of $p_0(\mathbf{r},{\lambda }_i)$ decreased with depth due to optical attenuation. This resulted in severely underestimating total hemoglobin concentration at depths deeper than 10 mm when applying spectral linear unmixing directly to $p_0(\mathbf{r},{\lambda }_i)$, rather than ${\mu }_a(\mathbf{r},{\lambda }_i)$. On the other hand, in the estimation of oxygen saturation ($C_{{\mathrm{HbO}}_2}(\mathbf{r})/C_{t\mathrm{Hb}}(\mathbf{r})$), the effect of the optical attenuation could not be canceled out because of its dependence on wavelengths, i.e., the spectral coloring effect. Thus, the classification accuracy constantly decreased with depth without optical fluence normalization in Fig. 12(b). The proposed method ameliorated such reduction [Fig. 12(b)].

Furthermore, the value of the effective optical attenuation coefficient (Fig. 8), which was estimated from the in vivo 3D OAT breast images using the proposed method, correlates well with experimental measurements ($\approx 1{\mathrm{cm}}^{-1}$) that were reported in previous studies.^44,68–70 The proposed method is completely measurement-data-driven, therefore, a prior knowledge of the optical properties of the breast tissues, anatomy of the vascular network, and precise characterization of the illumination pattern and incident fluence is not required.

The proposed method was specifically implemented for the 3D OAT breast imaging system presented in Fig. 2. However, the general framework for the normalization of the optical fluence distribution is not limited to breast imaging and to this specific system. For example, the curve fitting for incident optical fluence estimation can be opportunely modified to account for different optical illumination patterns.

Although these studies demonstrated qualitative and quantitative enhancement achieved via use of the proposed method, there remain limitations. First and foremost, the proposed method assumes a constant effective optical attenuation coefficient when estimating the fluence map within the breast. Errors in the estimation of the fluence map due to neglecting spatial variations of effective optical attenuation coefficient may introduce bias in the optical energy absorption estimates.

Besides, to obtain 3D quantitative images of the vascular blood oxygenation from in vivo data, further investigations should address acoustic heterogeneity of breast tissue and noise removal (thermal acoustic noise from the medium and thermal noise from ultrasound transducer and electromagnetic interference). Because the proposed method compensates for the depth-dependent optical attenuation by amplifying the image brightness at each voxel in the reconstructed 3D OAT images as a function of depth, existing noise and artifacts are also amplified depending on the depth. Application of the proposed method to images reconstructed using advanced regularization techniques can reduce such noise and artifacts, thus extending imaging depth.

Future directions include investigation of the proposed method to imaging of breasts with benign and malignant lesions and other 3D OAT imaging applications, such as transcranial imaging and small animal imaging (whole or partial body). It is expected that the performance of the proposed method largely depends on the distributions of the optical properties within the target. For example, in whole mouse imaging, hemoglobin-concentrated organs, such as liver, kidneys, and colon, locally occupy a certain extent as bulk. This causes a locally varying imbalance in the optical fluence distribution. In such case, the assumptions of the proposed method are invalid, thus, further investigation is required, including the use of more sophisticated numerical models to estimate the fluence distribution, such as MC photon transport simulation or simplified spherical harmonics approximation of radiative transfer equations.¹²

7.

Conclusion

In this work, a straightforward physics-based method to normalize optical fluence distributions in 3D OAT breast images was proposed. The method is based on generally accepted assumptions on breast anatomy and optical properties as well as common features of light delivery in existing 3D OAT breast imagers. In the proposed method, both distributions of incident optical fluence and optical attenuation within the breast tissues are estimated solely from the voxel brightness in the reconstructed images, thus, a prior knowledge of the breast and specific geometry of the light-delivery system is not required.

Numerical studies demonstrated that the proposed method—in conjunction with spectral linear unmixing—significantly enhanced blood vasculature detectability and improved estimation accuracy of vascular blood oxygenation down to a depth of 30 mm, when compared with no optical fluence normalization and a general-purpose image contrast enhancement technique called CLAHE. In addition, the proposed method outperformed CLAHE, in terms of PSNR and SSIM. It was also demonstrated that the proposed method can be applied to in vivo data. In particular, the effective optical attenuation coefficients estimated from the reconstructed 3D OAT breast images via the proposed method were found to be consistent with those experimentally measured in in vivo studies. With further investigations on acoustic heterogeneity, noise removal, and vascular segmentation, the use of the proposed method can potentially achieve 3D in vivo functional OAT images of the whole breast.