2.1.

SonixEmbrace ABUS

The SonixEmbrace transducer array consists of 384 elements on a 11.5-cm long circular arc with a curvature of $-12\mathrm{cm}$. It has a center frequency of 10 MHz and 90% bandwidth. The transducer dome has a radius of curvature of 12 cm and a diameter of 18 cm. A photograph of the dome is shown in Fig. 1. This large transducer requires multiplexing (mux) hardware that occupies a significant amount of space directly behind the transducer array under the dome (see Fig. 2), limiting conventional illumination options for breast PAT such as a collimated, expanded beam.^19,29 In further contrast to these systems, the SonixEmbrace uses a thin ($\sim 5\mathrm{mm}$) gel pad for acoustic coupling, as opposed to a flexible membrane or liquid-filled chamber. This means that there is little provision for any divergent optics to expand to cover a significant area of the tissue surface. We explored the possibility of using diode-based illumination, as the smaller footprint of such systems could potentially place the illumination closer to the surface. These systems would provide lower pulse power, so we chose to use a laser-based system in the interest of maximizing fluence and, therefore, the attainable penetration depth. We can also use an optical parametric oscillator (OPO) to perform multi-spectral PAT without changes to the output optics, whereas a diode system would require one illuminator per wavelength.

Fig. 1

The SonixEmbrace ABUS system, consisting of the ultrasound PC and scanning platform. Inset shows a top-down view of the transducer dome, with the blue transducer array visible.

Fig. 2

Schematic diagrams showing (a) the position of the ABUS illuminator and the electromechanical components beneath the ABUS imaging dome and (b) the illuminator, with the cutaway view revealing the optical elements.

The mechanical constraints imposed by the configuration of the ABUS electronics and motor led us to consider a narrow, fan-shaped beam that could fit between the multiplexer boards and the motor axle. We 3D printed an identical transducer dome in poly(lactic acid) (PLA; Afinia, Chanhassen, Minnesota), which has a polycarbonate optical window epoxied parallel to the transducer, as shown in Fig. 2.

2.2.

Laser and Optics

Our illumination system consists of a Continuum Surelite II LASER (Continuum Inc, Santa Clara, California) coupled to an OPO from the same manufacturer to control the output wavelength in the range of 675 to 2500 nm. We use a 80/20 beam splitter (68-376, Edmund Optics, Barrington, New Jersey) to direct 20% of the OPO output to a USB power meter (EnergyMax J-50MB-YAG – Coherent Inc., Santa Clara, California), allowing us to measure per-pulse energy while imaging. The remainder of the beam has a diameter of 9.5 mm, and is subsequently homogenized using a cross cylindrical lens array (Nr.18-00142, SUSS MicroOptics, Hauterive, Switzerland) and coupled by means of a plano-convex spherical lens (LA1608, Thorlabs, Newton) into a 1-mm diameter silica core optical fiber. Our beam homogenization and coupling system are described in further detail in Ref. 30. This system provides 5-ns pulses at 700 nm with a per-pulse energy at the output end of the fiber exceeding 50 mJ.

At the output end of the fiber is an optical assembly contained in a 3D-printed housing to maintain alignment and direct the diffuse output toward the optical window in the ABUS dome. The fiber is connected to the optical assembly by means of a ferrule connector (FC) fiber bulkhead connector (Thorlabs, Newton, New Jersey) and a bearing to decouple the fiber from the dome rotation. The fiber output first passes through an achromatic collimating doublet (65-438, Edmund Optics, Barrington, New Jersey) and is then incident upon an engineered “line diffuser” (EDL-$100\times 0.4$, RPC Photonics, Rochester, New York), which shapes the collimated beam into a fan. This fan has a 90° divergence along one axis and a 0.24° divergence along the other. This shape permits the beam to pass through the optical window, while illuminating its entire 11.5-cm length. A cross-sectional view of this assembly is shown in Fig. 2.

2.3.

Data Acquisition and Pre-Processing

Radio frequency (RF) data from the ABUS transducer were acquired using a SonixDAQ module (BK Medical). The SonixDAQ has a sampling rate at 40 MHz, a 12-bit resolution, and a $-10\text{-}\mathrm{dB}$ sensitivity in the range 2 to 20 MHz. Since the SonixDAQ can only acquire 128 channels at a time, each plane acquisition requires three acquisitions to cover all of the transducer elements. Control of the ABUS motor, and all triggering and synchronization of the optics and DAQ acquisition was accomplished using a custom Arduino-based circuit controlled over USB from a Windows PC.

This circuit allows flashlamp triggers to continually be sent to the laser at its pulse repetition rate of 10 Hz, while illumination pulses (via q-switch triggers) can be requested via hardware or software triggers. We have found that this “steady-state” operation decreases inter-pulse energy variability. For the pulse energies that we employ in this study, this variability is on the order of 10% to 15%. To further minimize the effect of this variability, we use the power meter data to normalize each 128-element acquisition by the relative pulse energy. A diagram summarizing the system configuration is shown in Fig. 3.

Fig. 3

The experimental setup used in this study. Optical elements T, L, H, and S are a translation stage, lens, beam homogenizer, and beam splitter, respectively. Dotted lines represent serial communications, and dashed lines are synchronization triggers.

Several steps are needed to pre-process the RF data prior to photoacoustic reconstruction. First, each A-line is low-pass filtered to remove high-frequency components outside the working range of the transducer. This was performed in software using a fifth-order Butterworth filter with a cutoff at 14.5 MHz. We retain the low-frequency components at this stage as they are essential for the following processing step: the application of singular value decomposition (SVD) denoising to remove cross-channel noise bands.³¹ To perform this denoising, we compute the SVD of a $128\times 2080$ matrix in which each row corresponds to the RF data for one transducer element. In this representation, this band-shaped noise dominates the first few singular value components. We can isolate these components and revert to the original representation, producing a matrix consisting of only the noise. This can then be subtracted from the original data, resulting in the denoised matrix. Since these noise bands are stochastic, this decomposition and denoising must be performed for each 128-element acquisition. Using the GPU implementation of SVD available in the scikit-cuda Python package,³² this denoising takes only 100 ms for each $128\times 2080$ array. Any low-frequency components outside the working range of the transducer will be removed by a ramp filter described in Eq. (3).

We explored several other denoising approaches including template matching, directional filtering³³ (both to the radio-frequency and image data), and temporal averaging (i.e., acquiring multiple frames of RF data per transducer position). Of these, temporal averaging was most successful at reducing the prominence of the noise bands, although it comes at a cost of increasing acquisition time. Our 10-Hz laser system and SonixDAQ create significant overhead for each additional acquisition, especially since the ABUS array requires three 128-element acquisitions per transducer position. For a volumetric scan of 200 angular positions, acquiring a single frame per position resulted in a 20-min scan time, whereas acquiring 10 frames per position increased the scan time to well over 1 h.

The other denoising approaches were either very sensitive to input parameters or unacceptably detrimental to the image quality. We have found that the SVD-based approach works on a wide variety of data with the same input parameters, removes the noise bands with minimal effect on the photoacoustic signals, and is computationally efficient to apply to the data. We have included an example comparing the effectiveness of directional filtering and SVD-based denoising in the Supplementary Materials.

2.4.

Reconstruction

When a laser pulse with spatial energy distribution $\mathrm{\Psi }(\mathbf{r},t)$ at time $t$ and positions $\mathbf{r}$ is incident upon a sample with coefficient of thermal expansion $\beta$, heat capacity $C_v$, acoustic speed $v_s$, and optical absorption coefficient ${\mu }_a(\mathbf{r})$, an acoustic wave is generated with initial pressure distribution

Assuming the medium properties $\beta$, $C_v$, and $v_s$ are spatially homogeneous and the laser pulse is temporally short, $\mathrm{\Psi }(\mathbf{r},t)=\mathrm{\Psi }(\mathbf{r})\delta (t)$, we have

Finding $p_0$ in Eq. (3) is typically accomplished using back-projection algorithms such as delay-and-sum and is a well-studied problem.³⁵ More sophisticated reconstruction algorithms take into account the fact that the measured pressure data are not exactly equal to the actual pressure incident on the detector due to factors such as the electromechanical and spatial impulse responses (SIRs) of the detector.^35,36 We can generalize the measured pressure data as

Since the SonixEmbrace transducer was designed for ultrasound imaging, it includes an acoustic lens for beam focusing, and as such it would be a poor assumption to fully neglect the SIR, as is common in many PAT reconstruction schemes. It is known, however, that accounting for the SIR is computationally difficult.^35,36 In an attempt to strike a balance, we separate the SIR into a spatial portion that describes the directional sensitivity and a temporal portion that describes the averaging effect due to the finite size of the transducer element

Under this assumption, $S_0({\mathbf{r}}_{\mathbf{D}},\mathbf{r})$ has no time dependence and as such does not depend on $|\mathbf{r}-{\mathbf{r}}_{\mathbf{D}}|$ when substituted into Eq. (3). This assumption also implies that the transducer has a flat frequency response. In the coordinate system of one transducer, we define ${\theta }_L$ and ${\theta }_E$ as the angles between the element normal (the axial direction) and the sample point in the lateral (in-plane) and elevational (out-of-plane) directions, respectively, such that

Since we are already computing the total angle between the element normal and the sample point to calculate the solid angle element $d{\mathrm{\Omega }}_0$ in Eq. (3), $S_0({\theta }_L,{\theta }_E)$ does not appreciably change the computational effort required. Since ${\theta }_L$ and ${\theta }_E$ must be computed for every pair of ${\mathbf{r}}_{\mathbf{D}}$ and $\mathbf{r}$ ($\sim {10}^5·{10}^7={10}^{12}$), this approach lends itself well to parallelization using GPUs.

2.4.1.

Non-uniform illumination

If the illumination source is spatially uniform within the sample, i.e., $\mathrm{\Psi }(\mathbf{r})=\mathrm{\Psi }$, Eq. (2) implies that ${\mu }_a(\mathbf{r})$ is directly proportional to the initial pressure $p_0(\mathbf{r})$. If the illumination is non-uniform, we must deconvolve these two spatial distributions. Since strong light attenuation in tissue limits imaging depth, it is common in PAT to apply an exponential weighting to the reconstructed $p_0(\mathbf{r})$ to enhance signals further from the surface. Our illumination is not collimated (uniform) at the tissue surface, but it can still be accounted for in a similar manner.

With sufficient geometric constraints, we model $\mathrm{\Psi }(\mathbf{r})$ as a solution of the RTE for homogeneous medium with temporally short pulses.³⁷ For an illumination source given by ${\mathbf{r}}_{\mathbf{i}}$, the optical fluence at a point $\mathbf{r}$ will be a function of the distance $\mathrm{r}=|\mathbf{r}-{\mathbf{r}}_{\mathbf{i}}|$ from the source:

Eq. (9)

$$\mathrm{\Psi }(\mathbf{r}-{\mathbf{r}}_{\mathbf{i}})=\frac{1}{4\pi Dr}\exp (-{\mu }_{\mathrm{eff}}r),$$

where the effective attenuation coefficient ${\mu }_{\mathrm{eff}}$ and the photon diffusion coefficient $D$ are defined as

Eq. (10)

$${\mu }_{\mathrm{eff}}=\sqrt{\frac{{\mu }_{a0}}{D}}D=\frac{1}{3({\mu }_{a0}+{\mu }_{s^{\prime }})}{\mu }_{s^{\prime }}={\mu }_s(1-g),$$

with ${\mu }_{a0}$ and ${\mu }_s$ being the bulk absorptive and scattering attenuation coefficients, respectively, and $g$ being the tissue anisotropy. Equation (9) holds in the diffusion limit, when ${\mu }_{a0}\ll {\mu }_s$. For human breast tissue at 700 nm, ${\mu }_s=10.0{\mathrm{cm}}^{-1}$ and ${\mu }_{a0}=0.2{\mathrm{cm}}^{-1}$.³⁸ Assuming constant optical properties in the tissue constitutes a first order approximation to the true fluence, which would require either prior knowledge or joint reconstruction of the scattering and absorption coefficients.

The geometry of the illuminator and the resultant fluence at the sample surface is shown in Fig. 4. For a given position in the sample, $\mathbf{r}$, the fluence can be computed by integrating Eq. (9) over values of ${\mathbf{r}}_{\mathbf{i}}$ on $B$, the section of the dome surface bounded by the angles $\theta$ and $\varphi$ corresponding to the two diverging axes of the fan-shaped beam:

Eq. (11)

$$\mathrm{\Psi }(\mathbf{r})={\iint }_B\frac{1}{4\pi Dr}\exp (-{\mu }_{\mathrm{eff}}r)\mathrm{d}\theta \mathrm{d}\varphi .$$

Fig. 4

The coordinate system for the illumination system, indicating the divergence angles of the beam and the intersection with the sample surface.

2.5.

Reconstruction Implementation

Figure 5 shows a summary of our reconstruction scheme. We first filter and back-project the measured signals to recover the initial pressure, taking into account the spatial sensitivity of the transducer array. We then compute the illumination map and use it to normalize the initial pressure, recovering the optical absorption map for the current illuminator position. It is important to reiterate that, since the illuminator is fixed relative to the transducer and not the sample, $\mathrm{\Psi }(\mathbf{r})$ must be recalculated for each acquisition angle. By applying this method for each angle, and combining the results, the optical absorption coefficient in the entire 3D tissue volume can be reconstructed. This accumulation over acquisition angles means that the final estimate of ${\mu }_a$ cannot be linearly separated into $\mathrm{\Psi }$ and $p_0$ as would be the case with systems employing stationary illumination.

Fig. 5

A schematic outline of our reconstruction scheme. Green boxes indicate operations computed on GPUs.

For all results in the present study, we collected RF data at 200 planes with a constant angular spacing of 1.8°. We structured our GPU reconstruction into a number of blocks equal to the number of image voxels, $i=1,\cdots ,N$, with each 384-element acquisition plane divided into 384 threads. This was performed sequentially for each of the planes, with the results being accumulated in GPU memory.

We computed the illumination map given by Eq. (11) at each of the $N$ image voxels by dividing the surface $B$ into small patches and summing their contributions to the fluence within the sample. These contributions are independent and can thus be computed in parallel on a GPU. For the present study, we divided $B$ into 200 patches, each of which was computed in a separate GPU thread, resulting in a $200\times$ reduction in computation time.

Reconstruction was implemented in Python, with the following computationally-intensive portions implemented in CUDA: back-projection, filtering, SVD denoising, and calculation of the illumination map $\mathrm{\Psi }(\mathbf{r})$. Reconstruction was performed on a 64-bit Windows 10 PC with 16 GB of memory, an Intel Core i7-7700 processor (Intel Corporation, Santa Clara, California), and a GeForce GTX 1060 3GB GPU (Nvidia Corporation, Santa Clara, California).

All reconstructions in this paper have a resolution of 0.5 mm and assume a constant acoustic velocity of $1500\mathrm{m}{\mathrm{s}}^{-1}$. To calculate the illumination profile using Eqs. (10) and (11), we used the scattering and absorptive attenuation coefficients of milk at 700 nm³⁹—${\mu }_s=30{\mathrm{cm}}^{-1}$ and ${\mu }_{a0}=0.015{\mathrm{cm}}^{-1}$ for the calibration phantom and of human breast tissue at 700 nm³⁸—${\mu }_s=10.0{\mathrm{cm}}^{-1}$ and ${\mu }_{a0}=0.2{\mathrm{cm}}^{-1}$ for the chicken breast phantom.

2.6.

Phantoms

To provide a robust and reproducible way to test our entire imaging system and reconstruction scheme, we designed and built a modular photoacoustic phantom.

As an absorber, we used 0.25-mm diameter monofilament fishing line, coated with black spray paint. This provided a strong photoacoustic signal, and the small diameter served as a test of the system resolution. To suspend the fishing line in the imaging volume, we designed and 3D printed a scaffold offering multiple mounting points between which the line was extended. The scaffold consisted of a base that precisely aligns the phantom with the imaging volume and several cylindrical columns with regularly-spaced holes through which the fishing line was threaded. This scaffold is shown in Fig. 6. As shown in Fig. 6, the fishing line was threaded through multiple points, providing a complex inclusion geometry that is nonetheless simple to describe and reproduce. For imaging the calibration phantom, the ABUS dome was filled with a 50/50 mixture of water and milk to provide optical scattering. To account for this mixture, we scaled the scattering and absorptive attenuation coefficients in our reconstruction by a factor of two, giving ${\mu }_s=15{\mathrm{cm}}^{-1}$ and ${\mu }_{a0}=0.0075{\mathrm{cm}}^{-1}$.

Fig. 6

Images of the calibration phantom used in this study. The white template is 3D printed plastic, and black painted fishing line serves as a photoacoustic source. The circular portion of the template has an outer diameter of 19 cm.

An additional advantage to using a 3D printed template is that we can easily use the design file to generate the inclusion geometry for use in simulations and validation of our reconstructed data.

To test our system in a more highly-scattering medium, we imaged a piece of 0.5-mm diameter pencil graphite embedded 2-cm deep in a piece of chicken breast tissue.

3.1.

SVD Denoising

Figure 7 shows an example of the effect of the denoising method described in Sec. 2.3 on representative RF data. For our data, discarding the first 10 singular values produced the best balance between removing noise and maintaining signal, as characterized by measuring the SNR for an A-line containing a known inclusion signal.

Fig. 7

The effect of SVD denoising on the RF data. Right panels show the signal along the red line for each of the three images. Noise signal corresponds to the first 10 singular values. (a) The raw data without processing. (b) The isolated noise. (c) The subtraction of (a) and (b).

3.2.

Transducer Directivity

We measured the spatial sensitivity of the ABUS transducer, $S_0({\theta }_L,{\theta }_E)$, by attaching a 2 mm long, 0.7 mm diameter graphite rod on a piece of fishing line to a 3-axis translation stage such that it was suspended in the water-filled ABUS dome. We coupled our laser system into the optical fiber as described in Sec. 2.2 and mounted the output end on the translation stage, directed toward the graphite inclusion. This ensured relatively constant illumination as the inclusion was translated. Photoacoustic data was acquired at a total of 200 points, over 1 cm in the elevational direction and 1.6 cm in the axial direction. The lateral position was fixed near the midpoint of the transducer since there are sufficient elements to constitute a range of lateral positions.

We observed a Gaussian dependence of the signal intensity on the angle between the transducer normal and the inclusion position, both in the lateral and elevational directions. The measured standard deviations of these distributions, ${\sigma }_L$ and ${\sigma }_E$, corresponding to the lateral and elevational angles ${\theta }_L$ and ${\theta }_E$ was 18.5° and 6.2°, respectively. These data and the associated fits are included in Fig. 8.

Fig. 8

The sensitivity of the SonixEmbrace transducer elements as a function of the elevational and lateral angle between the element normal and the imaging point. Red lines indicate Gaussian fit.

3.3.

Combined Imaging

Figure 9 shows representative in-plane images of the calibration phantom, demonstrating registered B-mode and photoacoustic imaging. For clarity of display, the photoacoustic data were thresholded such that only pixels with amplitudes in the top 10% of the image histogram were visible. The full dynamic range of the non-thresholded photoacoustic image was 10 dB. Figure 9(c) shows the model inclusion data derived from the design file for the 3D printed wire phantom, which further confirms that the designed structure of the phantom matches what we see when imaging. It is worth noting that the diameter of the inclusions in the model data in Fig. 9(c) are exaggerated $10\times$ to be more easily visible. There are two circular features in the model data, corresponding to wires oriented normal to the imaging plane and therefore generating photoacoustic signals too weak to be above the display threshold. We used two measures to quantify the agreement between any two of these registered modalities. The first was the root-mean-squared error (RMSE), defined as

Fig. 9

B-mode and photoacoustic data for a representative plane of the calibration phantom, indicating registration with modeled inclusion positions derived from the wire phantom design file.

Table 1

Registration errors between imaging modalities.

3.4.

Illuminator Spatial Profile

Figure 10 shows the calculated illumination structure within the sample, as per Eq. (11). We also measured the fluence at the sample side of the optical window to estimate the overall losses in our fiber coupling and diffusing system. The measured data are shown in Fig. 10(b), alongside the modeled fluence for the same location. The broader shape of the measured data is due to the 1-cm aperture of the power meter used. Integrating the fluence over the entire window indicates a total energy per pulse of around 35 mJ, down from 100 mJ at the OPO output. If this same laser pulse were collimated to the size of the dome (surface area $250{\mathrm{cm}}^{-2}$), the average fluence at the tissue surface would be $0.14\mathrm{mJ}{\mathrm{cm}}^{-2}$.

Fig. 10

Optical fluence for one acquisition plane. (a) The two-dimensional structure of the illumination, relative to the transducer face, which is shown in black. (b) The measured fluence at the sample side of the optical window in black, compared with its predicted value, as a function of the position along the arc of the window.

3.5.

Reconstruction Comparison

To test our reconstruction scheme, we reconstructed our measured data both with and without accounting for transducer directivity and non-uniform illumination. We then measured SNR of a known inclusion for each method. We measured SNR here by extracting the photoacoustic amplitude along a line perpendicular to the long axis of the inclusion. We found the maximum of this signal and defined the “peak” as the region where the amplitude is greater than 50% of this maximum. We then defined the SNR as the mean amplitude of the signal in this region divided by the standard deviation of the amplitude over the entire image.⁴² The width of this peak region is the full width at half maximum (FWHM) and provides a measure of the resolving power of the system. Figure 11 shows an example of this measurement from the chicken breast phantom. A comparison of the measured values is given in Table 2. Figures 11(a)–(d) were scaled to the same dynamic range for display to make qualitative differences in the images clearer.

Fig. 11

Comparison of reconstruction methods applied to chicken breast phantom. Insets show the PA amplitude along the dashed white line, normalized to the same dynamic range across all four images, with vertical dashed red lines indicating the FWHM of the peak. Corresponding measurements are summarized in Table 2.

Table 2

Comparison of reconstruction methods applied to chicken breast phantom (see Fig. 11).

To further quantify the performance of our reconstruction scheme, we compared the image error in this reconstruction as a function of imaging depth in our calibration phantom. Using the known model for the inclusions, $M_i$, we define image error as

Fig. 12

Image error as a function of imaging depth, compared for reconstruction schemes with and without transducer directivity and/or non-uniform illumination compensation.

Fig. 13

Maximum intensity projections of a tomographic reconstruction of the calibration phantom. (a) The data reconstructed with directivity and illumination compensation; (b) the model data; and (c) an overlay of the two. White scale bar indicates 1 cm.