2.1.

Light Delivery System

From in vivo measurements with PAM 1, operating at an excitation wavelength of 1064 nm, the absorption contrast could not be attributed to hemoglobin only. Water, lipids, and collagen were also expected to contribute to the contrast.¹² For PAM 2, a dual-wavelength approach was chosen. Excitation with 1064 nm was chosen for its deep penetration into tissue.¹⁸ Water absorption is relatively high at this wavelength,¹⁹ potentially a limiting factor in the penetration depth given the presence of water in the breast and the use of water as a coupling medium between breast and detectors. A second excitation wavelength is needed to distinguish contributions of two chromophores. Excitation with 755 nm was implemented for it has higher hemoglobin absorption and lower water absorption.²⁰

In PAM 2, the woman lies prone on a bed while one of the breasts is pendant through an aperture, hanging inside the imaging tank. The breast is illuminated with a large laser beam from the bottom (the ventral or nipple side of the breast). There is also illumination via nine optical fiber bundles, which are directed at parts of the breast close to the chest wall. A dual-laser system (Qcombo, Quanta System, Milan, Italy) comprising an Nd:YAG and an Alexandrite Q-switched laser is used. The pulse durations are 5 ns at 1064 nm (Nd:YAG) and 60 ns at 755 nm (Alexandrite). The system can be operated in single-wavelength mode or in double-wavelength mode. In single-wavelength mode, one of the wavelengths is emitted with a repetition rate of 10 Hz. In double-wavelength mode, the system alternates between the two wavelengths at a 20-Hz repetition rate.

Each side illumination fiber bundle has a concave lens (focal length of 14.7 mm) to diverge the beam. The fiber bundles are directed slightly upward, under an angle of 5 deg, in order to illuminate the region of the breast close to the chest wall. The beam from the bottom is also diverged by a concave lens (focal length of 25 mm). Inside the laser system, the main laser beams are split into two beams—for bottom and side illumination—by means of a beam splitter (Quanta System, Milan, Italy). By exchanging the beam splitter for another, the energy distribution can be altered. In one setting, 67% of the laser light illuminates the breast through the large bottom beam; the remaining 33% is distributed equally over the nine fiber bundles (67/33 illumination ratio). Another energy division is possible, where 50% of the total energy illuminates from the bottom and the other 50% from the sides through the fiber bundles (50/50 illumination ratio). For practical reasons, only two configurations have been considered so far. The laser pulse energies (measured with StarBright power meter and PE50BF-DIFH-C sensor, Ophir Optronics Solutions Ltd., Jerusalem, Israel) of the bottom illumination beam at a 50/50 configuration are 160 and 130 mJ of the Nd:YAG and Alexandrite laser, respectively. At the 67/33 configuration, the bottom illumination energy was not measured at the same time point but can be expected to be $\sim 220$ and 180 mJ for Nd:YAG and Alexandrite, respectively (with a tolerance of 8%). Laser fluence on the skin is below the safety limit of $20\mathrm{mJ}/{\mathrm{cm}}^2$ at all times (NEN-EN-IEC 60825-1:2014). The smallest safe distances between skin and laser outputs inside the tank were calculated. The design of the imaging tank does not allow the breast to come that close to the laser outputs, ensuring a safe situation. Eye protection goggles are worn whenever the laser is emitting.

2.2.

Ultrasound Detection

The single elements of the ultrasound detector system (Imasonic SAS, Voray sur l’Ognon, France) have a $3.5\times 3.5{\mathrm{mm}}^2$ active surface area. The distance between two adjacent elements is 1.4 mm, resulting in a center to center distance of 4.9 mm. One array consists of 32 piezocomposite elements, placed on an arc with a radius of 120 mm (Fig. S1 in the Supplementary Materials). A total of 12 such arrays were installed in PAM 2.

A waterproof stainless steel case encloses an array and its integrated preamp (PA Imaging R&D B.V., Enschede, The Netherlands). The 32-channel low-noise ultrasound amplifiers with a 38-dB voltage gain employ a high-input impedance and filtering to prevent aliasing. The outer surface of the elements is covered with aluminum tape with a thickness of $30\mu \mathrm{m}$, to prevent light from being absorbed by the detector surface and thus generating a photoacoustic signal.

2.3.

Imaging Configuration

The fiber bundles and detector arrays are fixed inside the imaging tank. Detector arrays were placed in a circular arrangement inside the imaging tank, with an angle of 30 deg between two adjacent arrays. Fiber bundles and detector arrays are installed in an interleaved manner, but the flexible design of the tank allows other configurations also (e.g., transmission mode with all detectors on one side). The current configuration is depicted in Fig. 1. During a measurement, the tank with the fiber bundles and detectors rotates in steps around the breast to acquire more views (projections) around the breast for a tomographic scan. This increases the number of sampling locations and thereby improves the reconstructed image quality. At each position, an averaged signal is obtained from a number of laser pulses. This number is a parameter that can be chosen per measurement. After acquiring the average signal, the tank rotates a few degrees before stopping and obtaining another measurement. The choice for such a stop-and-go scanning method over a continuous scanning method²¹ was made to limit the time needed for reconstruction (linearly increasing with the number of projections). The centers of the arcs of all detector arrays coincide, such that all single elements lie on the surface of a virtual hemisphere. The elements of six arrays are shifted 2.45 mm on the virtual hemisphere with respect to the elements of the other six arrays, to enhance the spatial sampling of the breast volume. The positions of all 384 detector elements were primarily deduced from the designed geometry of each arc-shaped array and the mounting positions inside the imaging tank (30 deg between two adjacent arrays). Subsequently, small corrections to the elements’ positions were made based on photoacoustic calibration measurements with a small absorbing dot made from blank epoxy (G14250, Thorlabs, Inc., Newton, New Jersey) mixed with carbon black, suspended in the imaging tank. In every measurement, the exact angular rotation of the imaging tank between consecutive scanning steps is measured at each projection (MagLine system, SIKO GmbH, Buchenbach, Germany) since the actual motor rotation was found to deviate from the programmed rotation. In Fig. S2 in the Supplementary Materials, the positions (where tank rotation was measured with the MagLine system) of the centers of all $12\times 32$ elements are plotted (Lambert azimuthal equal-area projection), for a tomographic scan obtained from 21 projection angles, covering a total angle of 60 deg.

Fig. 1

(a) Top view drawing of the imaging tank and (b) photograph of the imaging tank.

Prior to a measurement, the subject climbs onto the bed via a step stool to lie prone on a bed with one of her breasts through the aperture with a diameter of 21 cm and pendant in the imaging tank. The imaging tank is filled with water to achieve acoustic coupling between the breast and the ultrasound detectors. The bed originally belonged to a prone breast stereotactic biopsy system (Hologic, Inc., Marlborough, Massachusetts). The imaging tank is double-walled with an inner diameter of 35.5 cm and an outer diameter of 40.5 cm. The gap between the two walls carries water heated externally in a heating circulator (F12-ED, Julabo GmbH, Seelbach, Germany). By heat transfer through the inner wall, the coupling water in the imaging tank is warmed up to maintain a temperature of 25°C to 30°C during the measurement. Figure 2 shows the PAM 2 system.

Fig. 2

The PAM 2 system with (a) the bed on which a woman lies prone, (b) a foot rest, and (c) the breast aperture with (d) a schematic of the shape and position of the imaging tank below the bed, (e) the power supply and cooling unit, (f) laser head (behind the panel), (g) DAQ unit, and (h) step stool.

2.4.

Data Acquisition

Photoacoustic data are acquired through a 12 times 32 channels analog to digital converter as part of a data acquisition (DAQ) system (PA Imaging R&D B.V., Enschede, The Netherlands). Sampling is done with a frequency of 25 MHz at a depth of 14 bits and for $167\mu \mathrm{s}$ following each laser pulse, to detect signals from over a distance of 25 cm [assuming a speed of sound (SOS) of 1500 m/s]. Average signals as well as raw signals following each laser pulse are acquired, transferred to and saved on a PC. Raw signals are not used for reconstructing images but are saved to enable more detailed analysis of the data, for example with respect to movement of the breast.

2.5.

Signal Processing, Image Reconstruction, and Image Visualization

For the reconstruction, a filtered backprojection based on Algorithm 1 by Haltmeier et al.²² was implemented in Python 3.4. The time interpolation step mentioned in Ref. 22 has not been included in our algorithm, given the relatively small sampling interval of 40 ns compared to the period of the acoustic waves involved. Deconvolution of the signals was performed, as a preprocessing step, in order to reconstruct the pressure waves received by the transducers as accurately as possible. Deconvolution is implemented as a regularized division of the measured signal by the average frequency response of the transducers. Regularization refers to an attenuation of the effect of division for frequencies where the average frequency response is very low. The relative sensitivities of all transducer elements are used to scale the measurement data before being passed to the reconstruction algorithm. The SOS is assumed to be homogeneous inside the measurement volume (breast and surrounding water) and calculated from the measured temperature of the water surrounding the breast.²³ To keep reconstruction time within feasible limits, the core of the algorithm was implemented as a C++ extension module with support of OpenMP, a programming interface for parallel computation.²⁴ A single-image reconstruction, i.e., for one breast or phantom and one laser wavelength (300 MB of averaged data), takes about 5 to 10 min on a dual-socket server platform with two octa-core Intel Xeon processors with 32 GB RAM. Without the C++ extension, calculation time would be over 5 h per reconstruction.

Reconstructed 3-D volumes are by default visualized as maximum intensity projections (MIPs). Another way of visualizing a reconstructed volume is by means of local maximum intensity projections (LMIPs).²⁵ The first pixel above a set threshold is plotted, where the pixels are color-coded for distance from the observer (white is superficial and orange is further away). The depth corresponding to the orange color is different per image, depending on the size of the reconstructed volume along each axis.

3.1.

Detector Relative Sensitivity

All detector elements were insonified with a 1-MHz V303 ultrasound transducer (Olympus NDT, Inc., Waltham, Massachusetts), driven at 1 MHz with a fixed pressure in water at room temperature. The V303 transducer was driven by an AFG3102 function generator (Tektronix, Inc., Beaverton, Oregon). The 32 detecting elements of each array were individually insonified by rotating the array around the V303 transducer, keeping the distance between detector and V303 fixed (120 mm). Alignment was optimized by translating the detector array to maximize the amplitude of the detected signal. The average was taken of 100 individual signals. From all measured signals, the amplitudes were extracted and normalized to the maximal value measured in any of the detectors to obtain relative sensitivities.

3.2.

Detector Frequency Response

The frequency responses of the detecting elements were measured by acquiring signals induced by a custom-built laser-induced ultrasound (LIUS) setup.²⁶ This technique of LIUS is reviewed in Ref. 27. The frequency content of the LIUS signal was broad enough (frequencies up to 10 MHz were measured when the source was characterized) to assume the signal as an ideal acoustic calibration source. Hence, the measured PAM 2 detector response is considered to only represent the detector response and we did not need to correct for the source frequency content. The LIUS transmitter comprised polydimethylsiloxane mixed with carbon black and applied in a thin layer ($29\mu \mathrm{m}$) onto a 5-mm substrate of poly(methyl methacrylate). The LIUS transmitter was illuminated with a 10-mm diameter beam from a frequency-doubled Quanta-Ray Pro 250-10-Hz laser (Spectra-Physics, Santa Clara, California) operating at a wavelength of 532 nm. The 32 elements of each array were individually characterized by rotating the array around the LIUS transmitter, the distance between detector and LIUS transmitter being fixed at 120 mm. Detector arrays, LIUS transmitter, and illumination fiber were immersed in water at room temperature during these measurements. Each signal was averaged 100 times. The acquired LIUS signal was Fourier transformed to obtain the frequency response of the detector.

3.3.

Resolution

The spatial resolution of the system was estimated by measuring the full-width at half-maximum (FWHM) of the point-spread function (PSF) of a subresolution strongly absorbing sphere. A batch of 26 drops made of two component blank epoxy (G14250, Thorlabs Inc., Newton, New Jersey) mixed with carbon black was produced. These drops were suspended on two perpendicular nonabsorbing nylon threads to create a test object (Fig. 3). The smallest drop (as identified by eye) had a diameter of 0.46 mm, the largest drop 1.65 mm and the remaining drops ranged in diameter size from 0.96 to 1.56 mm. The test object was placed in the center of the imaging tank and illuminated with 1064 nm with 21 projection angles (total covered angle of 60 deg) and 100 averages per laser pulse. The 67/33 illumination ratio was used. Images were reconstructed with a grid size of 0.3 mm and are shown as MIPs. Gaussians were fitted to the PSFs of line profiles through the center of the smallest (subresolution) sphere with a diameter of 0.46 mm. The lateral and elevational resolutions were estimated as the FWHM of the PSF’s Gaussian fits.

Fig. 3

Photograph of the test object with 26 different sized epoxy drops on nonabsorbing nylon.

3.4.

Healthy Volunteer Measurement

All human volunteer measurements were carried out as part of a study approved by an institutional review board (METC Twente, Medisch Spectrum Twente, Enschede) and registered in the Netherlands National Trial Register (TC 6508). To all subjects participating in our study, the procedure was explained and consent was obtained.

To test the system and the measurement protocol, a measurement on a 29-year-old woman with a light skin tone was performed (healthy volunteer 1). The 67/33 illumination ratio was employed. A tomographic scan with 21 projection angles (total covered angle of 60 deg) and 100 laser pulses per average was done. Water temperature inside the imaging tank was 30°C. The scan took 4 min per breast. Images were reconstructed with a 0.4-mm grid size and are plotted as MIPs. Differences between images obtained with 755 nm versus 1064 nm were compared.

3.5.

Optimizing the Illumination Scheme

In the PAM 2 system, illumination of the breast is spread over a total of 10 beams. The main beam illuminates the nipple side of the breast with a diverging beam from the bottom of the tank. The remaining energy is spread over nine fiber bundles that are installed at the top of the imaging tank, serving to illuminate breast tissue close to the chest wall. As mentioned, a beam splitter inside the laser splits the main beam and sets the energy division. In the original configuration, 67% of the energy was delivered to the nipple side of the breast. In order to find out whether a more homogeneous intensity distribution could be obtained, another configuration with more energy going to the side illumination was considered. Both configurations were compared by measuring on the breasts of a 25-year-old healthy volunteer with a light skin tone (healthy volunteer 2).

After the first measurement set with a 67/33 configuration, the beam splitter inside the laser system was replaced with a 50/50 beam splitter and a second measurement set was performed. The subject left the bed between these measurements and thus the breast was repositioned in between measurements. Leaving the bed was necessary since the system needed to be partially disassembled to reach the laser and replace the beam splitter, a procedure that took more than an hour. In all measurements, a tomographic scan over 21 projection angles (total covered angle of 60 deg) was done where signals were averaged over 100 laser pulses. The water temperature was 30°C. The acquisition time for one breast was 4 min. Reconstructed images (grid size = 0.4 mm) are plotted as MIPs. Another reconstruction was performed with a 0.25-mm grid size and plotted as LMIPs. Images were compared, where the focus was on the visibility of blood vessels near the chest wall and the homogeneity of signal intensity throughout the breast. Measurements on this volunteer were also used to estimate the imaging depth relative to the skin surface.

4.1.

Detector Relative Sensitivity

The sensitivity of all 384 single elements was measured. All values were normalized and the interelement variability in terms of sensitivity is depicted in Table 1.

Table 1

Sensitivity interelement variability of the detectors.

4.2.

Detector Frequency Response

The mean impulse response of all LIUS signals, used for measuring the frequency response, is plotted (in red) together with the 384 ($12\times 32$) individual detector element impulse responses (in gray) in Fig. 4(a). The signals were shifted in time to start at $t=0\mathrm{s}$. The 384 corresponding Fourier transformed frequency responses (again in gray tones) can be seen in Fig. 4(b). The mean frequency response is plotted in red. The detector is most sensitive at a frequency of 1 MHz. The average $-6\text{-}\mathrm{dB}$ fractional frequency bandwidth of all detectors is 100%, ranging from 500 kHz up to 1.5 MHz. There is little variation in response among the individual detector elements (especially in the $-6\text{-}\mathrm{dB}$ frequency bandwidth). Therefore, the average frequency response was used in the deconvolution step of the reconstruction algorithm.

Fig. 4

Response to an acoustic impulse. (a) Measured photoacoustic signal traces generated by LIUS. The red line represents the mean impulse response and the gray lines are the individual impulse responses of the 384 elements. Normalized amplitude has arbitrary units. (b) Corresponding frequency spectrum: the red line represents the mean response and the gray lines are the corresponding individual frequency responses of the 384 elements. The dashed line represents the $-6\text{-}\mathrm{dB}$ level.

4.3.

Resolution

Reconstructed images plotted as MIPs of the test object with epoxy beads can be seen in Figs. 5(a)–5(c). Figure 5(d) [line drawn in Fig. 5(b)] estimates the elevational resolution to be 0.96 mm. In Fig. 5(e) [line drawn in Fig. 5(a)], the lateral resolution was measured to be 1.06 mm.

Fig. 5

MIPs of the resolution test object. (a) Top-view $XY$ projection: the line profile was drawn along the blue line, (b) side-view $XZ$ projection: the line profile was drawn along the red line, (c) side-view $YZ$ projection, (d) Gaussian fit to the line profile along $z$ in (b), and (e) Gaussian fit to the line profile along $y$ in (a).

4.4.

Healthy Volunteer Images

Figure 6 shows MIPs (coronal and sagittal views) of reconstructed data of the left breast of healthy volunteer 1, illuminated with both 755 nm [Fig. 6(a)–6(b)] as well as 1064 nm [Figs. 6(c)–6(d)]. With both wavelengths, the nipple is clearly visible (yellow arrowheads) and blood vessels around it. There is good anatomical correspondence between the images, due to the laser system alternating pulses of both wavelengths within one measurement. Some differences can, however, be observed. When illuminating with 755 nm, the portion of the breast illuminated by the bottom beam seems a bit larger than at 1064 nm [compare the structures indicated with blue arrowheads in Figs. 6(b) and 6(d)]. It is known that the bottom beam divergence of both wavelengths is not exactly the same, explaining this observation. In the 755-nm images, the skin surface around the nipple-areolar complex is better visible than at 1064 nm. This can be explained by the higher absorption coefficient of cutaneous melanosomes at 755 nm compared to 1064 nm.²⁰ Finally, illumination with 1064 nm shows structures close to the chest wall with more detail, indicated with green arrowheads in Figs. 6(b) and 6(d).

Fig. 6

Reconstructed MIP images of the left breast of healthy volunteer 1. (a) Coronal image obtained at 755 nm, (b) corresponding sagittal image, (c) coronal image obtained at 1064 nm, and (d) corresponding sagittal image. Yellow arrowheads indicate the nipple. Blue arrowheads indicate a part of the breast better visible when illuminating with 755 nm compared to 1064 nm. Green arrowheads indicate structures close to the chest wall, better visible when illuminating with 1064 nm compared to 755 nm.

4.5.

Optimization of the Illumination Scheme

In Fig. 7, MIPs of reconstructed data of the right breast of healthy volunteer 2 using 1064-nm light with two different illumination schemes (illumination ratio bottom/side 67/33 and 50/50) are shown. From both measurements, coronal [Figs. 7(a) and 7(c)] and sagittal images [Fig. 7(b) and 7(d)] are shown. In all images, a network of blood vessels (both arteries and veins) as well as the nipple [indicated with a yellow arrowhead in Figs. 7(b) and 7(d)] are visible. In general, images obtained with the 50/50 configuration show more structures. In the 67/33 configuration images, reconstructed signal intensity is more localized.

Fig. 7

Reconstructed MIP images of the right breast of healthy volunteer 2, illuminated with 1064 nm. Zoomed versions of regions within the colored boxes can be found in Fig. 8. (a) Coronal image with illumination ratio bottom/side at 67/33, (b) corresponding sagittal image, (c) coronal image with illumination ratio bottom/side at 50/50, and (d) corresponding sagittal image.

Four subsections of the coronal images were chosen and are shown enlarged in Fig. 8. They were not placed at the same coordinates in both images but were placed in a way that they show the same structures. All boxes intend to show portions of the breast close to the chest wall. The yellow box [Figs. 8(a) and 8(e)] shows a small loop of blood vessels, which is resolved in the 50/50 configuration but less so in the 67/33 configuration. When we look at the region bounded by the blue box, we again see more structures with the 50/50 configuration [Fig. 8(g)] compared to with the 67/33 configuration [Fig. 8(c)]. This could be partially explained by the fact that the $x$-position of the breast in the 50/50 configuration is more centrally located within the imaging volume. The $y$-position is, however, less favorable compared to in the 67/33 configuration. In the green box, the differences between the two configurations are small and both zoomed images are comparable. In the magenta box, the 67/33 configuration shows more structures, probably due to the more favorable coordinates both along $x$ as well as $y$.

Fig. 8

Zoomed subsections of reconstructed coronal MIP images of the right breast of healthy volunteer 2, illuminated with 1064 nm. Box colors indicate regions in Figs. 7. (a)–(d) Sections imaged with illumination ratio bottom/side at 67/33 and (e)–(h) sections imaged with illumination ratio bottom/side at 50/50.

4.6.

Imaging Depth

In the reconstruction of the left breast of healthy volunteer 2, we assessed the imaging depth by tracing a vessel toward the center of the breast. Figure 9 shows a reconstruction done with a smaller grid size, namely 0.25 mm. Images were plotted as LMIPs. Markers (green crosses) were placed along a vessel in the transverse view [Fig. 9(b)] up to the depth where the vessel was visible. The $x$-coordinates of those pixels were retrieved and subsequently the same vessel was marked in the sagittal view [Fig. 9(a)]. From this, the imaging depth in these images was estimated to be 22 mm with respect to the skin surface.

Fig. 9

Color-coded LMIPs of the left breast of healthy volunteer 2, illuminated with 755 nm. (a) Sagittal image, where max depth on the colorbar corresponds to 110 mm and (b) transverse image, where max depth on the colorbar corresponds to 100 mm. Green crosses track a blood vessel down to a depth of 22 mm from the skin surface.

4.7.

Photoacoustic Features Inside the Nipple

The measurements on healthy volunteer 2, carried out to assess the illumination configuration, led to some interesting observations on the photoacoustic appearance of the nipple. In the images in Fig. 7 [particularly in Figs. 7(c)–7(d)], a number of small high-intensity features in the nipple show up. Figure 10 shows images of the left breast illuminated of the same volunteer (2) at 755 nm with a 50/50 illumination ratio. Here these features are even more clear. The reconstruction was done with a smaller grid size, namely 0.25 mm, to further investigate the phenomenon. Images were plotted as LMIPs. Figure 10(a) shows the coronal LMIP. The area of interest is shown enlarged in Fig. 10(b). Sagittal and transverse LMIP images are plotted in Figs. 10(c) and 10(d), respectively.

Fig. 10

Color-coded LMIPs of the left breast of healthy volunteer 2, illuminated with 755 nm. (a) Coronal image, where max depth on the colorbar corresponds to 55 mm, (b) the zoomed subsection shows the features that are present in a circular pattern inside the nipple (three of them are pointed out), (c) sagittal image, where max depth on the colorbar corresponds to 110 mm, and (d) transverse image, where max depth on the colorbar corresponds to 100 mm.