Translator Disclaimer
24 October 2019 Twente Photoacoustic Mammoscope 2: system overview and three-dimensional vascular network images in healthy breasts
Author Affiliations +

We present the Twente Photoacoustic Mammoscope 2, a photoacoustic breast imaging system employing a tomographic configuration. It images one breast pendant inside an imaging tank filled with water while a woman lies prone on a bed. A dual-head laser (755 and 1064 nm) illuminates the breast with one beam directed at the nipple and nine beams directed at the sides. Ultrasound signals are detected using 12 arc-shaped arrays, each curving along the pendant breast. Each array comprises 32 piezocomposite elements each with a center frequency of 1 MHz. The imaging tank and the ultrasound arrays rotate around the breast in steps to obtain additional multiple projections. Three-dimensional images are reconstructed using a filtered backprojection algorithm. The system is described in detail, and measurements on a test object are presented. As part of a preliminary study to assess the system’s in vivo performance, the breasts of two healthy volunteers were imaged. These images show the breast contour, the nipple, and the vascular anatomy within the breast. In the nipple of one case, multiple high-intensity “hot spots” are observed, which we suspect are associated with the lactiferous ducts terminating in the nipple.



Breast cancer is the most frequently diagnosed cancer as well as the leading cause of cancer death among women worldwide.1 Imaging is an important aid in the detection, diagnosis, and treatment monitoring and follow-up of breast cancer.2 Conventional imaging techniques currently used in the clinic are x-ray mammography (MMG), ultrasonography (US), and magnetic resonance imaging (MRI).3 MMG is the only method utilized in organized screening programs in several countries around the globe since it is the only method proven to reduce mortality.4 It does, however, require (often painful) breast compression, makes use of ionizing radiation, and the detection sensitivity in dense breasts is relatively low.5 US is used only in the diagnostic trajectory and is often employed for targeted examinations of palpable lesions and lesions detected in screening. It is also the primary method used for guiding biopsies. It does not use ionizing radiation but has the drawback of being operator dependent.6 MRI is used in the diagnostic trajectory, often as a problem-solving modality. It does not use ionizing radiation, has high resolutions, and possesses high sensitivity. However, it requires the administration of contrast agents, does not have a high specificity, is expensive, and is not widely disseminated.3,7 Each technique has its strengths but also its weaknesses, leaving room for new imaging techniques that can potentially overcome the current limitations and drawbacks.

Photoacoustics (PA) is a technique where the absorption contrast of tissue constituents to light is combined with the resolution of ultrasound imaging.8,9 Upon illuminating the tissue of interest with pulsed light, optical energy absorption at specific locations leads to a local pressure rise due to the effect of thermoelastic expansion. The pressure transient propagates as ultrasound, which can be measured outside the tissue with ultrasound transducers. Hemoglobin as an absorber provides functional contrast to allow visualization of (newly formed) vascularization associated with cancer progression. The production of new blood vessels is called angiogenesis; inducing angiogenesis is thought to be one of the biological capabilities of tumors and a hallmark of cancer.10,11

With the first generation Twente Photoacoustic Mammoscope, the PAM 1 system, we have shown visualization of breast malignancies in vivo using PA.1214 In PAM 1, a woman lay prone on a bed with one breast pendant through an aperture. Below the aperture, the breast was slightly compressed between a glass plate on the cranial side of the breast and a flat detector array on the caudal side of the breast. The planar imaging configuration worked in forward mode, illuminating the breast through the glass plate with 1064 nm and detecting generated ultrasound with a two-dimensional (2-D) ultrasound detector array from the caudal side. In the latest clinical studies, 32 of 33 measured breast malignancies13 were visualized with good contrast. One of the drawbacks of PAM 1 was the limited field of view (90×80  mm2 detector surface area), precluding whole breast coverage in one measurement. Further, the PAM 1 configuration required slight compression of the breast for acoustic coupling between breast and detector and to obtain a uniform breast thickness. This compression is unfavorable with respect to patient comfort. Moreover, the flat detector array in the forward mode configuration made accurate three-dimensional (3-D) visualization of lesions challenging due to the limited views. To overcome these limitations, a tomographic system was designed and built, which we refer to as the second generation Twente Photoacoustic Mammoscope (PAM 2).

Several other research groups are working in the field of photoacoustic tomographic breast imaging. Matsumoto et al.15 showed a (breast) imaging system (PAI-04) where a woman is prone on a bed with one breast inside a breast-holding cup. A hemispherical-shaped detector array (center frequency 4 MHz) is used to measure the ultrasound as the breast is illuminated by a large beam, directed at the ventral side of the breast. The system is scanned in a spiral pattern within a plane to provide full breast images. Lin et al.16 have recently presented an instrument referred to as the single-breath-hold photoacoustic computed tomography breast imaging system, where again a woman lies prone on a bed. The breast to be imaged is slightly compressed upward against the chest wall. A full-ring detector array (2.25 MHz) is scanned elevationally to obtain a 3-D image. Also here, the breast is illuminated by one beam from beneath the breast, directed at the ventral side of the breast. Oraevsky et al.17 co-registered anatomical US images with functional PA images in their LOUISA-3D system, which images a breast inside an imaging bowl while the woman lies in prone position on a bed. An arc-shaped detector array (50 kHz to 6 MHz) is rotated around the breast while multiple optical fibers positioned on an arc-shaped holder, which can rotate around the breast, are used for illumination.

We contribute to the field of photoacoustic breast imaging with our PAM 2 system, which is different from other systems recently presented1517 in terms of illumination and detector specifications. We show that adding more illumination beams and changing the ratio of delivered energy to these beams alters the observed results. Also we present images of absorbing features in the nipple that, to our knowledge, have not been reported earlier in photoacoustic breast imaging. The PAM 2 system contains an imaging tank filled with water, in which a breast is pendant while the woman lies in prone position on a bed. The system has low frequency (1 MHz), broadband detector arrays that are arc-shaped and placed around the breast inside the imaging tank, such that the arcs follow the breast contour. Illumination is done not only with a beam (which we call bottom beam) directed at the ventral side of the breast, as in Refs. 15 and 16, but also with nine smaller beams directed at the region of the breast close to the chest wall. The performance of the system is demonstrated on a test object and on the breasts of two healthy volunteers. Reconstructed images show clearly visible vascularization. Images obtained with two illuminating wavelengths (755 and 1064 nm) are qualitatively compared. We give an indication of the imaging depth. In the breasts of one healthy volunteer, we observe, for the first time, structures that might be linked to lactiferous ducts terminating at the nipple.


System Description


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.12 For PAM 2, a dual-wavelength approach was chosen. Excitation with 1064 nm was chosen for its deep penetration into tissue.18 Water absorption is relatively high at this wavelength,19 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.20

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 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  mJ/cm2 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.


Ultrasound Detection

The single elements of the ultrasound detector system (Imasonic SAS, Voray sur l’Ognon, France) have a 3.5×3.5  mm2 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  μm, to prevent light from being absorbed by the detector surface and thus generating a photoacoustic signal.


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 method21 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×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.



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  μ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.


Signal Processing, Image Reconstruction, and Image Visualization

For the reconstruction, a filtered backprojection based on Algorithm 1 by Haltmeier et al.22 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.23 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.24 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).25 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.


Materials and Methods


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.


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.26 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  μ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.



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.



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.


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.




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.

Maximal value1.0
Minimal value0.84
Mean value0.92
Standard deviation0.03


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×32) individual detector element impulse responses (in gray) in Fig. 4(a). The signals were shifted in time to start at t=0  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-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-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-dB level.




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).



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.20 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.



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.



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.



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.



Discussion and Conclusion

With PAM 2, we have noninvasively acquired 3-D images of vasculature in the human breast. Photoacoustic signals are detected in a tomographic scan by 384 single elements placed around the breast with a center frequency of 1 MHz and a 6-dB fractional bandwidth of 100%. The imaging capabilities were explored by measuring on a test object comprising of highly absorbing beads on a thread and by measuring on two healthy volunteers. A resolution of 1  mm along all axes was reached.

PAM 2 was developed for noninvasively imaging vascularity in the female breast and thereby visualizing breast cancer. It is well-known from the literature that tumor growth induces new vessels to grow in the process of angiogenesis.11 These blood vessel networks are disorganized; tumor vasculature contains vessels that are immature, tortuous, and hyperpermeable. The vessels are characterized by bulges and blind ends.28,29 So far, we have shown images of healthy volunteers only. Since we are able to visualize these vessel networks in great detail, the expectation is to be able to identify regions with higher vessel densities and further find abnormal structures and thereby visualize breast cancer.

Additionally, photoacoustic imaging can potentially extract information about blood oxygen saturation using the knowledge that the absorption spectra of hemoglobin and oxyhemoglobin are markedly different at certain wavelengths. It is clinically relevant to assess oxygen saturation of blood vessels in or near a breast tumor since this is expected to be lower in tumor tissue than in healthy tissue.3032 Matsumoto et al.15 adopted a method to extract the so-called S-factor from measurements at two wavelengths (756 and 797 nm), in which a measure of oxygen saturation could be evaluated. Oraevsky et al.17 have shown binary functional images with oxygenated and deoxygenated blood obtained by illumination with 757 and 797 nm.

We have shown that images obtained with our PAM 2 system at two illumination wavelengths show a good anatomical overlap of visualized structures. The fact that our laser system alternates pulses of the two wavelengths leads to this correspondence between the images obtained at both wavelengths. The isosbestic point of oxy-and deoxyhemoglobin lies around 800 nm. Since our system employs illumination at 755 and 1064 nm, this opens possibilities for analysis of oxygen saturation. For oxygen saturation assessment, it is important to know (at least an estimate of) the light fluence distribution throughout the breast. Image intensities have not been related to absolute pressures nor have been corrected for fluence variations yet. We are planning to implement a method where a fluence distribution estimation for both wavelengths based on the diffusion approximation can be used to evaluate contributions of oxy- and deoxyhemoglobin to the absorption contrast throughout the reconstructed volume. Apart from the possibility for extracting functional information, the two wavelengths are complementary since they result in slightly different images of the same breast volume, for example, the different appearances of the skin surface in images obtained with the two wavelengths.

The breast is not only illuminated with a large bottom beam, but also with multiple side beams directed at the chest wall. Slightly more structures were visualized when delivering a higher amount of energy to the side illumination (at the cost of bottom illumination). It is, however, not easy to directly compare images with both illumination configurations [Figs. 7(a)7(b) versus Figs. 7(c)7(d)] since the breast was repositioned in between the measurements and thus its position was not the same. Moreover, the amount of movement of the volunteer in both measurements was difficult to control and complicates comparison. It would be of interest to study the added value of the side illumination and the optimal distribution of light energy and compare it to a situation with only a bottom illumination beam.

Since pendant breast diameters (at the chest wall) can be as high as 13.7 cm (mean for D/DD cup sizes),33 reaching an imaging depth of around 6 to 7 cm would be required to image the entire breast. Imaging depth was estimated from reconstructed images obtained from healthy volunteer 2 (Fig. 9) to be 22 mm, which is in the same range as recent results from Toi et al.34 The estimation is based on the continuity of a vessel, so it does not necessarily reflect the deepest retrieved signal and is, therefore, meant as an indication. We believe that our estimated imaging depth is smaller than what our system would be capable of. Currently, the breast hangs freely inside the imaging tank and this may lead to (part of) the breast not being imaged optimally, when it is outside the region that is being illuminated and/or from where photoacoustic signals are detected. Moreover, the SOS utilized in the reconstruction algorithm is not physically correct. Currently, one SOS is assumed, based on the temperature-dependent SOS in the water surrounding the breast inside the imaging tank. We know, however, that the SOS in breast tissue is different from that of water and that within the breast, different tissue components all have their own SOS. The signals coming from relatively deep in the breast are most affected by a physically incorrect SOS; hence, the image quality deeper in the breast will be decreased most by assuming only one SOS, ultimately limiting our imaging depth. Finally, breathing and motion artifacts may be of influence given the measurement duration. In PAM 2, the breast is pendant in water. From investigating individual time-domain signals, it was observed that the breast drifts during the measurement duration of 4 min. Perhaps volunteers are tense at the beginning of a measurement and slowly relax toward the end of it. Small periodic movements due to breathing were found as well. Abrupt motion of the volunteer is also possible but is always registered by the PAM 2 operator and taken into account when analyzing reconstructed images. In order to enhance depth performance, individual time-domain signals might be further studied and employed to make corrections for motion in the future. Another option would be to immobilize the breast during a measurement. Oraevsky et al.17 used a hemispherical, optically clear and acoustically thin plastic cup-stabilizer. Lin et al.16 used an agar pillow to slightly compress the breast against the chest wall from below and measure within a single-breath hold to minimize breathing motion distortion. Matsumoto et al.15 used a breast holding cup to immobilize the breast. When choosing any form of breast immobilization, a challenge would be to find a shape that supports every breast shape enough to reduce motion artifacts while at the same time having all the breast tissue anatomically available to be imaged. Especially imaging the areas of the breast close to the chest wall and axilla might be compromised when a suboptimal immobilization strategy is chosen. In designing such an immobilization strategy, it is important that a material is used that distorts the optical illumination and acoustic wave propagation as little as possible. Immobilizing the breast could also contribute to a larger portion of the breast being visualized by forcing the breast to always be in the center of the imaging tank. The position of the breast is currently not reproducible [compare breast positions in Figs. 7(a)7(b) and 7(c)7(d)] since the breast hangs freely inside the imaging tank. Furthermore, with an immobilization strategy, it will be easier to employ a simple two SOS model in the reconstruction algorithm once the water-tissue boundary is known. This would improve our reconstruction algorithm and the visibility of deep features. Also having recognized the skin surface will simplify calculation of the depth of features with respect to the skin. In the long run, it will be beneficial to obtain a 3-D SOS map of the breast tissue simultaneously with the photoacoustic measurement, such that the map can be used to further improve the image reconstruction of that particular measurement. Efforts are made in this field, for example, by Matthews and Anastasio,35 who proposed a method for a joint reconstruction of photoacoustic initial pressure and SOS distributions from photoacoustic and few-view ultrasound tomography measurements. PAM 2 depth performance will also be enhanced by employing an iterative reconstruction algorithm in addition to the currently used filtered backprojection algorithm. The relatively weak signals from deeper inside the breast will be retrieved from the background signal and enhanced. Both the improvement of the reconstruction algorithm as well as the immobilization of the breast are topics of current investigation.

In volunteer 2, we observed intriguing features inside the nipples of both breasts. Figures 7(c)7(d) show the features in the right breast, and Fig. 10 shows the same features with more detail in the left breast. We hypothesize that these features belong to the lactiferous ducts terminating from the nipple through orifices. The lactiferous ducts in the nipple are surrounded by smooth muscle cells,36,37 which are vascularized. In general, ducts in the human breast are surrounded by a high number of microvessels.38 We can observe a total of either seven or eight of these high-intensity features. From previous work,39 the average breast was estimated to contain five to nine ductal orifices. The features are within an area with a diameter of 12 mm. The nipple (excluding the areola) of volunteer 2 was measured to be 13 mm. Given that these features are thus inside the nipple and the number matches the estimated number of ductal orifices, it is plausible that they indeed belong to the lactiferous duct orifices. This is worth exploring further in the near future. Moreover, it could be valuable to learn more about the anatomy of the nipple in a noninvasive manner, with PA. Knowledge about (ductal and vascular) nipple anatomy is of importance for ductal approaches to diagnosis, planning nipple-sparing mastectomy40 and for understanding normal breast development.41 Finally, in the literature, there is still some discussion about and a large spread in the number of orifices in the nipple.39,40,42,43 More measurements, if possible on the same volunteer but also on other volunteers, will be carried out to further investigate these findings.

With PAM 2, we have demonstrated 3-D images of vascular structures in human breasts. This vascular-morphological information combined with the potential of extracting functional information on oxygen saturation in the future could be of clinical use not only in breast cancer management but also in visualizing different arterial and venous disorders and abnormalities of the breast.44 Although most of them uncommon, accurate diagnosis often demands multimodality imaging and pathologic correlation. Some benign tumors are vascularized (such as hemangiomas), which might complicate differentiation from malignant lesions. It is important to investigate different vascular disorders and abnormalities of the breast and whether PA could have an added value in this field.

PAM 2 system’s parameters are unique in terms of the illumination scheme and detector characteristics, and it is worth exploring how it performs in measuring on breast tumors in the near future and how results compare with the other systems around the world. So far, promising results have been obtained in the form of 3-D images of vascular patterns in the breasts of healthy volunteers. The knowledge gained from exploring the system and measuring on healthy women will be used for imaging on breast cancer patients in a noninvasive and painless manner.


Authors Schoustra, Piras, Steenbergen, and Manohar have no relevant financial interests in this paper and no potential conflicts of interest to disclose. PA Imaging R&D B.V. is a 100% subsidiary of PA Imaging Holding B.V. PA Imaging Holding B.V. is a privately held company. Authors Huijink, op ‘t Root, Alink, and Muller Kobold (employees of PA Imaging R&D B.V.) do not have financial interests in the company.


Authors acknowledge Stichting Achmea Gezondheidszorg for funding in Project No. Z620, the Health∼Holland LSI-TKH PPP project SACAMIR (No. LSHM17007), and the PAMMOTH project funded by the EU’s Horizon 2020 RIA H2020 ICT 2016–2017 under Grant Agreement No 732411. Authors would like to thank Johan van Hespen, Rutger Pompe van Meerdervoort, and Daan Sprünken for their contributions.



L. A. Torre et al., “Global cancer statistics, 2012,” CA Cancer J. Clin., 65 (2), 87 –108 (2015). CAMCAM 0007-9235 Google Scholar


P. A. T. Baltzer et al., “New diagnostic tools for breast cancer,” Mag. Eur. Med. Oncol., 10 (3), 175 –180 (2017). Google Scholar


L. W. Bassett, S. Lee-Felker, “Breast imaging screening and diagnosis,” The Breast, 337361.e2 5th ed.Elsevier(2018). Google Scholar


D. A. Berry et al., “Effect of screening and adjuvant therapy on mortality from breast cancer,” N. Engl. J. Med., 353 (17), 1784 –1792 (2005). NEJMAG 0028-4793 Google Scholar


D. S. Al Mousa et al., “What effect does mammographic breast density have on lesion detection in digital mammography?,” Clin. Radiol., 69 (4), 333 –341 (2014). Google Scholar


R. J. Hooley, L. M. Scoutt and L. E. Philpotts, “Breast ultrasonography: state of the art,” Radiology, 268 (3), 642 –659 (2013). RADLAX 0033-8419 Google Scholar


B. E. Adrada, R. Candelaria and G. M. Rauch, “MRI for the staging and evaluation of response to therapy in breast cancer,” Top. Magn. Reson. Imaging, 26 (5), 211 –218 (2017). TMRIEY 0899-3459 Google Scholar


S. Manohar and D. Razansky, “Photoacoustics: a historical review,” Adv. Opt. Photonics, 8 (4), 586 –617 (2016). AOPAC7 1943-8206 Google Scholar


P. Beard, “Biomedical photoacoustic imaging,” Interface Focus, 1 (4), 602 –631 (2011). Google Scholar


D. Hanahan and J. Folkman, “Patterns and emerging mechanisms of the angiogenic switch during tumorigenesis,” Cell, 86 (3), 353 –364 (1996). CELLB5 0092-8674 Google Scholar


D. Hanahan and R. A. Weinberg, “Hallmarks of cancer: the next generation,” Cell, 144 (5), 646 –674 (2011). CELLB5 0092-8674 Google Scholar


M. Heijblom et al., “Photoacoustic image patterns of breast carcinoma and comparisons with magnetic resonance imaging and vascular stained histopathology,” Sci. Rep., 5 11778 (2015). SRCEC3 2045-2322 Google Scholar


M. Heijblom et al., “The state of the art in breast imaging using the Twente Photoacoustic Mammoscope: results from 31 measurements on malignancies,” Eur. Radiol., 26 (11), 3874 –3887 (2016). Google Scholar


S. Manohar et al., “Initial results of in vivo non-invasive cancer imaging in the human breast using near-infrared photoacoustics,” Opt. Express, 15 (19), 12277 –12285 (2007). OPEXFF 1094-4087 Google Scholar


Y. Matsumoto et al., “Visualising peripheral arterioles and venules through high-resolution and large-area photoacoustic imaging,” Sci. Rep., 8 (1), 14930 (2018). SRCEC3 2045-2322 Google Scholar


L. Lin et al., “Single-breath-hold photoacoustic computed tomography of the breast,” Nat. Commun., 9 (1), 2352 (2018). NCAOBW 2041-1723 Google Scholar


A. A. Oraevsky et al., “Full-view 3D imaging system for functional and anatomical screening of the breast,” Proc. SPIE, 10494 104942Y (2018). PSISDG 0277-786X Google Scholar


B. J. Tromberg et al., “Non-invasive in vivo characterization of breast tumors using photon migration spectroscopy,” Neoplasia, 2 (1–2), 26 –40 (2000). Google Scholar


G. M. Hale and M. R. Querry, “Optical constants of water in the 200-nm to 200-μm wavelength region,” Appl. Opt., 12 (3), 555 –563 (1973). APOPAI 0003-6935 Google Scholar


S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol., 58 (11), R37 (2013). PHMBA7 0031-9155 Google Scholar


A. Sharma, S. K. Kalva and M. Pramanik, “A comparative study of continuous versus stop-and-go scanning in circular scanning photoacoustic tomography,” IEEE J. Sel. Top. Quantum Electron., 25 (1), 1 –9 (2019). IJSQEN 1077-260X Google Scholar


M. Haltmeier, T. Schuster and O. Scherzer, “Filtered backprojection for thermoacoustic computed tomography in spherical geometry,” Math. Methods Appl. Sci., 28 (16), 1919 –1937 (2005). Google Scholar


W. Marczak, “Water as a standard in the measurements of speed of sound in liquids,” J. Acoust. Soc. Am., 102 (5), 2776 –2779 (1997). JASMAN 0001-4966 Google Scholar


, available from: (2018) October ). 2018). Google Scholar


Y. Sato et al., “Local maximum intensity projection (LMIP): a new rendering method for vascular visualization,” J. Comput. Assisted Tomogr., 22 (6), 912 –917 (1998). JCATD5 0363-8715 Google Scholar


M. Olsman, “A new variant of plane wave imaging using laser-induced ultrasound,” University of Twente, (2016). Google Scholar


S.-L. Chen, “Review of laser-generated ultrasound transmitters and their applications to all-optical ultrasound transducers and imaging,” Appl. Sci., 7 (1), 25 (2017). Google Scholar


D. W. Siemann, “The unique characteristics of tumor vasculature and preclinical evidence for its selective disruption by tumor-vascular disrupting agents,” Cancer Treat. Rev., 37 (1), 63 –74 (2011). CTREDJ 0305-7372 Google Scholar


A. R. Pries et al., “Structural adaptation and heterogeneity of normal and tumor microvascular networks,” PLoS Comput. Biol., 5 (5), e1000394 (2009). Google Scholar


M. Höckel and P. Vaupel, “Tumor hypoxia: Definitions and current clinical, biologic, and molecular aspects,” J. Natl. Cancer Inst., 93 (4), 266 –276 (2001). JNCIEQ Google Scholar


P. Vaupel, F. Kallinowski and P. Okunieff, “Blood flow, oxygen and nutrient supply, and metabolic microenvironment of human tumors: a review,” Cancer Res., 49 (23), 6449 –6465 (1989). Google Scholar


V. Ntziachristos et al., “MRI-guided diffuse optical spectroscopy of malignant and benign breast lesions,” Neoplasia, 4 (4), 347 –354 (2002). Google Scholar


S. Y. Huang et al., “The characterization of breast anatomical metrics using dedicated breast CT,” Med. Phys., 38 (4), 2180 –2191 (2011). MPHYA6 0094-2405 Google Scholar


M. Toi et al., “Visualization of tumor-related blood vessels in human breast by photoacoustic imaging system with a hemispherical detector array,” Sci. Rep., 7 41970 (2017). SRCEC3 2045-2322 Google Scholar


T. P. Matthews and M. A. Anastasio, “Joint reconstruction of the initial pressure and speed of sound distributions from combined photoacoustic and ultrasound tomography measurements,” Inverse Probl., 33 (12), 124002 (2017). INPEEY 0266-5611 Google Scholar


S. A. Hoda et al., Rosen’s Breast Pathology, Wolters Kluwer Health, Philadelphia (2014). Google Scholar


M. Tezer et al., “Smooth muscle morphology in the nipple-areola complex,” Nipple-Areolar Complex Reconstr., 28 171 –175 (2011). Google Scholar


A. G. Naccarato et al., “Definition of the microvascular pattern of the normal human adult mammary gland,” J. Anat., 203 (6), 599 –603 (2003). JOANAY 0021-8782 Google Scholar


S. M. Love and S. H. Barsky, “Anatomy of the nipple and breast ducts revisited,” Cancer, 101 (9), 1947 –1957 (2004). CANCAR 0008-543X Google Scholar


J. E. Rusby et al., “Breast duct anatomy in the human nipple: three-dimensional patterns and clinical implications,” Breast Cancer Res. Treat., 106 (2), 171 –179 (2007). BCTRD6 Google Scholar


J. J. Going and T. J. Mohun, “Human breast duct anatomy, the ‘sick lobe’ hypothesis and intraductal approaches to breast cancer,” Breast Cancer Res. Treat., 97 (3), 285 –291 (2006). BCTRD6 Google Scholar


D. T. Geddes, “Inside the lactating breast: the latest anatomy research,” J. Midwifery Women’s Health, 52 (6), 556 –563 (2007). Google Scholar


F. Hassiotou and D. Geddes, “Anatomy of the human mammary gland: current status of knowledge,” Clin. Anat., 26 (1), 29 –48 (2013). CLANE8 1098-2353 Google Scholar


R. A. Jesinger et al., “Vascular abnormalities of the breast: arterial and venous disorders, vascular masses, and mimic lesions with radiologic-pathologic correlation,” Radiographics, 31 (7), E117 –E136 (2011). Google Scholar

Biographies of the authors are not available.

© The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Sjoukje M. Schoustra, Daniele Piras, Roeland Huijink, Tim J. P. M. op't Root, Laurens Alink, Wouter Muller F. Kobold, Wiendelt Steenbergen, and Srirang Manohar "Twente Photoacoustic Mammoscope 2: system overview and three-dimensional vascular network images in healthy breasts," Journal of Biomedical Optics 24(12), 121909 (24 October 2019).
Received: 4 May 2019; Accepted: 22 August 2019; Published: 24 October 2019

Back to Top