Passive element enriched photoacoustic computed tomography (PER PACT) for simultaneous imaging of acoustic propagation properties and light absorption

We present a „hybrid‟ imaging approach which can image both light absorption properties and acoustic transmission properties of an object in a two-dimensional slice using a computed tomography (CT) photoacoustic imager. The ultrasound transmission measurement method uses a strong optical absorber of small cross-section placed in the path of the light illuminating the sample. This absorber, which we call a passive element acts as a source of ultrasound. The interaction of ultrasound with the sample can be measured in transmission, using the same ultrasound detector used for photoacoustics. Such measurements are made at various angles around the sample in a CT approach. Images of the ultrasound propagation parameters, attenuation and speed of sound, can be reconstructed by inversion of a measurement model. We validate the method on specially designed phantoms and biological specimens. The obtained images are quantitative in terms of the shape, size, location, and acoustic properties of the examined heterogeneities. ©2010 Optical Society of America OCIS codes: (170.5120) Photoacoustic imaging; (170.0110) Imaging systems; (170.4580) Optical diagnostics for medicine; (170.3010) Image reconstruction techniques. References and links 1. L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nat. Photonics 3(9), 503–509 (2009). 2. S. Manohar, S. E. Vaartjes, J. C. G. van Hespen, J. M. Klaase, F. M. van den Engh, W. Steenbergen, and T. G. van Leeuwen, “Initial results of in vivo non-invasive cancer imaging in the human breast using near-infrared photoacoustics,” Opt. Express 15(19), 12277–12285 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=oe-15-19-12277. 3. S. A. Ermilov, T. Khamapirad, A. Conjusteau, M. H. Leonard, R. Lacewell, K. Mehta, T. Miller, and A. A. Oraevsky, “Laser optoacoustic imaging system for detection of breast cancer,” J. Biomed. Opt. 14(2), 024007 (2009). 4. T. D. Khokhlova, I. M. Pelivanov, V. V. Kozhushko, A. N. Zharinov, V. S. Solomatin, and A. A. Karabutov, “Optoacoustic imaging of absorbing objects in a turbid medium: ultimate sensitivity and application to breast cancer diagnostics,” Appl. Opt. 46(2), 262–272 (2007). 5. J. Jose, S. Manohar, R. G. M. Kolkman, W. Steenbergen, and T. G. van Leeuwen, “Imaging of tumor vasculature using Twente photoacoustic systems,” J Biophoton. 2(12), 701–717 (2009). 6. J. T. Oh, M. L. Li, H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Three-dimensional imaging of skin melanoma in vivo by dual-wavelength photoacoustic microscopy,” J. Biomed. Opt. 11(3), 34032 (2006). #134437 $15.00 USD Received 1 Sep 2010; revised 2 Nov 2010; accepted 8 Nov 2010; published 20 Jan 2011 (C) 2011 OSA 31 January 2011 / Vol. 19, No. 3 / OPTICS EXPRESS 2093 7. R. Nuster, M. Holotta, C. Kremser, H. Grossauer, P. Burgholzer, and G. Paltauf, “Photoacoustic microtomography using optical interferometric detection,” J. Biomed. Opt. 15(2), 021307 (2010). 8. J. Laufer, E. Zhang, G. Raivich, and P. Beard, “Three-dimensional noninvasive imaging of the vasculature in the mouse brain using a high resolution photoacoustic scanner,” Appl. Opt. 48(10), D299–D306 (2009). 9. H. P. Brecht, R. Su, M. Fronheiser, S. A. Ermilov, A. Conjusteau, and A. A. Oraevsky, “Whole-body threedimensional optoacoustic tomography system for small animals,” J. Biomed. Opt. 14(6), 064007 (2009). 10. R. Ma, A. Taruttis, V. Ntziachristos, and D. Razansky, “Multispectral optoacoustic tomography (MSOT) scanner for whole-body small animal imaging,” Opt. Express 17(24), 21414–21426 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-24-21414. 11. M. H. Xu and L. H. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrumn. 77, (2006). 12. X. Jin, and L. V. Wang, “Thermoacoustic tomography with correction for acoustic speed variations,” Phys. Med. Biol. 51(24), 6437–6448 (2006). 13. J. F. Greenleaf, and R. C. Bahn, “Clinical imaging with transmissive ultrasonic computerized tomography,” IEEE Trans. Biomed. Eng. 28(2), 177–185 (1981). 14. N. Duric, P. Littrup, L. Poulo, A. Babkin, R. Pevzner, E. Holsapple, O. Rama, and C. Glide, “Detection of breast cancer with ultrasound tomography: first results with the Computed Ultrasound Risk Evaluation (CURE) prototype,” Med. Phys. 34(2), 773–785 (2007). 15. S. Manohar, R. G. H. Willemink, and T. G. van Leeuwen, “Speed-of-sound imaging in a photoacoustic imager,” SPIE (2007), 64370R. 16. S. Manohar, R. G. H. Willemink, F. van der Heijden, C. H. Slump, and T. G. van Leeuwen, “Concomitant speedof-sound tomography in photoacoustic imaging,” Appl. Phys. Lett. 91 (2007). 17. R. G. H. Willemink, S. Manohar, Y. Purwar, C. H. Slump, F. d. Heijden, and T. G. v. Leeuwen, “Imaging of acoustic attenuation and speed of sound maps using photoacoustic measurements,” in (SPIE, 2008), 692013. 18. P. M. Joseph, and R. A. Schulz, “View sampling requirements in fan beam computed tomography,” Med. Phys. 7(6), 692–702 (1980). 19. F. Marinozzi, D. Piras, and F. Bini, “Spectral analysis of backscattered ultrasound field from hydroxyapatite granules,” in Advances in Medical, Signal and Information Processing, 2008. MEDSIP 2008. 4th IET International Conference on, 2008), 1–4. 20. J. C. Bamber, Acoustical characteristics of biological media (Encyclopedia of acoustics, Wiley, 1997), Vol. 4. 21. R. G. H. Willemink, S. Manohar, T. G. van Leeuwen, and C. H. Slump, “Estimation of integrated ultrasound transmission parameters I,” Speed of Sound under review. 22. R. G. H. Willemink, S. Manohar, C. H. Slump, F. van der Heijden, and T. G. van Leeuwen, “A maximum likelihood method for obtaining integrated attenuation from ultrasound transmission mode measurement,” J. Acoust. Soc. Am. 123(5), 3641 (2008). 23 R. G. H. Willemink, S. Manohar, T. G. van Leeuwen and C. H. Slump, “Estimation of integrated ultrasound transmission parameters II: Acoustic attenuation,” in preparation. 24. R. A. Kruger, P. Liu, Y. R. Fang, and C. R. Appledorn, “Photoacoustic ultrasound (PAUS)--reconstruction tomography,” Med. Phys. 22(10), 1605–1609 (1995). 25. C. R. Crawford, and A. C. Kak, “Multipath artifact corrections in ultrasonic transmission tomography,” Ultrason. Imaging 4(3), 234–266 (1982). 26. R. G. H. Willemink, S. Manohar, J. Jose, C. H. Slump, F. van der Heijden, and T. G. van Leeuwen, Simultaneous imaging of ultrasound attenuation, speed of sound, and optical absorption in a photoacoustic setup,” A. M. Stephen and D. h. Jan, eds. (SPIE, 2009), p. 72650J.


Introduction
Near-infrared photoacoustic (PA) imaging has garnered much attention [1] in imaging of various pathological states related to vascular condition and function. Some important applications of this technique include breast cancer detection [2][3][4][5], skin cancer visualization [6] and small animal imaging [7][8][9][10]. The technique relies on irradiating tissue with nanosecond pulses of visible or near infrared (NIR) light. Optical absorption in tissue causes thermoelastic expansion, which produces broadband pulses (MHz) of acoustic energy. These pulses propagate through the tissue with 2-3 orders lower scattering compared to light and can be detected at the tissue surface [11] at multiple positions. PA Imaging thus combines the typically high optical absorption contrasts of vascular-related pathologies with the high resolutions associated with ultrasound imaging.
The presence of an ultrasound detector in the typical PA imaging configuration, immediately suggests the possibility to perform additional imaging of the acoustic properties of the subject either simultaneous or sequential to the PA studies. Jin et al. [12] used an extra ultrasound transmitter in an additional measurement to perform transmission tomography in a conventional PA imaging system to obtain speed of sound (SOS) maps. This method can be likened to ultrasound computed tomography (USCT) which is being applied to image the female breast [13,14].
Recently, we showed [15,16] the feasibility of obtaining SOS tomograms without the need for an additional ultrasound transmitter and extra measurements. The method is based on creating ultrasound by the interaction of the laser light with a strong absorber fixed in its path in the imaging tank. The propagation times of the ultrasound transients produced by this "passive" element through the object at angles around 360° are measured using a ultrasound detector array. We reconstructed the corresponding spatial distribution of the SOS by crosscorrelating the times-of-flight (TOF) of the measured signal with a reference measurement, acquired with the object retracted from the imaging tank.
In addition to SOS imaging, acoustic attenuation tomograms can also be extracted from the "passive" element signals by analyzing amplitudes with and without the object [17]. Further, since the "passive" element has a small geometric cross-section, a major portion of the incident light illuminates the object allowing conventional PA computed tomography (CT) to be performed as well.
In this paper, we extend our previous work in SOS imaging and perform for the first time hybrid imaging of light absorption, speed-of-sound and acoustic attenuation. We provide a short overview of the imager which we refer to as the passive element-enriched photoacoustic computed tomography (PER-PACT) system. We summarize the methods used for estimating integrated US transmission properties of the subject per projection, and the reconstruction methods used. We finally show the hybrid imaging of specially designed phantoms and biological specimens.

Passive element enriched photoacoustic computed tomographic (PER-PACT) imager
A schematic of the PER-PACT set-up is depicted in Fig. 1. The light source is a Q-switched Nd:YAG laser (Brilliant B, Quantel) with a frequency doubler delivering 6 ns pulses. For the experiments described here, 532 nm was used as the exciting wavelength. The object is mounted on a rotary stage and immersed in an imaging tank with water. The curvilinear array (Imasonic, Besançon, France) consists of 32 piezocomposite elements arranged in a curved aperture covering an angle of 85° with radius of curvature 40 mm. The elements have a central frequency of 6.25 MHz with a receiving bandwidth greater than 80%. Individual elements have a dimension of 10 mm to 0.25 mm with an inter-element spacing of 1.85 mm. Each element in the curved surface is shaped to produce an elevation plane focus of 1 mm at a distance of 48 mm from the detector surface. The signals from each detector element are amplified by 60 dB with a 32-channel pulse-receiver system (Lecoeur, Paris) with a sampling rate of 80 MHz. Prior to each set of measurements, a calibration measurement was performed using a horsetail hair phantom, to ascertain the imaging geometry such as the centre of rotation and position of detector elements.
The passive element was a horsetail hair with a diameter of 250 µm positioned at a distance of 90 mm from the array (Fig. 1). The arrangement was such that the object lies completely within the fanbeam traced by the line of sight from passive element to the edge detector elements for all projections. The ultrasound detector measures conventional photoacoustic signals from the object in addition to the passive element signals. The two sets of signals have different times-of-flight (TOF); the interference of these two signals at the detector can be avoided by choosing a proper time window during the analysis. A reference measurement is made by retracting the object, allowing the ultrasound to propagate in water alone. The signals from reference and object measurements are analyzed yielding time-offlight (TOF) differences and acoustic attenuation changes compared with the reference. Obtaining multiple projections around the object permits reconstruction of the cross-sectional map of the acoustic attenuation and speed-of-sound.
The spatial resolution in fan beam geometry is dependent on the spacing between the "rays" in each "view". Here, "rays" refers to the single measurement of the relevant density function integral along a narrow acoustic beam from the passive element. While a "view" can be defined as a collection of such rays which together form the projection of the density as seen from a particular angular projection. The minimum number of detector elements N min or the angular "views" required to obtain a faithful image function without aliasing can be calculated using [18]: where m  is the required highest spatial frequency, R is the distance from the center point to the outermost pixel in the area of interest and  is the opening angle of the fan beam from passive element to the detector array.
In our experimental geometry ( m  = 2.5 lines mm 1 , R = 28 mm, and  = 40°), 1336 angular "rays" on the circle of diameter 80 mm are required. Considering 32 rays in one angular "view", leads to a minimum of 42 angular "views" required to obtain a good spatial resolution. In the measurements we typically acquire 45 views around the object using 15 mJcm 2 of energy at 532 nm with 100 averages.

Phantoms used
The base material used for the phantom was 3% agar gel in water. Acoustic inhomogeneities were prepared according to Marinozzi et al [19] using 4% agar gel in a 80% dilution of milk in water. To acoustically characterize the samples, prior to imaging we used the insertion technique [20] and obtained SOS and acoustic attenuation values as 1498 ms 1 and 0.15 dBcm 1 MHz 1 for the phantom, with 1505 ms 1 and 0.4 dBcm 1 MHz 1 for the inhomogeneities at room temperature (22° C). The acoustic values are consolidated in Table  1, and detailed descriptions of the phantoms are provided below: where s c is the SOS of the object, w c is the SOS of the water (1490 ms 1 at room temperature), d is the thickness (25 mm transverse dimension of mouse), and t  is the sampling interval of the data acquisition system (12.5 ns sampling frequency 80 MHz). The above equation leads to a minimum detectable SOS contrast ( c  ) of 1 ms 1 with the system, which is the integrated value between the passive element and detector element. To investigate the feasibility of mapping such an SOS contrast lumped into small targets, we developed a phantom embedded with acoustic inhomogeneities of various dimensions. Figure 2(a) depicts the layout of the phantom with its acoustic inhomogeneities; Fig. 3(a) shows the photograph of the phantom in preparation: Rod shaped agar-milk inserts have been placed on a pre-formed agar gel cylinder in a 25 mm diameter mould. The inserts are then submerged in aqueous agar solution which gels and forms the phantom. The cross-section of the inhomogeneities ranges from 2.6 mm to 4.6 mm. The smallest object of 2.6 mm (SOS 1505 ms 1 ) induces a TOF delay of 8 ns with respect to the phantom measurement, and is considered to be a sub-resolution target. The largest object of 4.6 mm cross-section induces a normalized TOF delay of 15 ns, greater than the sampling period Δt used. (b) Elevation plane resolution: To quantify the slice thickness, a cylindrical phantom comprising two regions of different SOS was made, as shown in Fig. 2(b). Imaging slices are made of the object on either side and across the interface.
(c) Hybrid imaging: To investigate the feasibility of hybrid imaging of light absorption, simultaneous with acoustic transmission parameters, a biological specimen composed of porcine muscle and adipose tissue was prepared. Figure 5

Estimation of integrated ultrasound transmission parameters
The first step to reconstruct SOS and acoustic attenuation tomograms is to estimate integrated TOF delays and attenuation along paths connecting the passive element with a detector element.
The sampled version of the signal in a measurement window, is then: Since the measurement function is a non-linear function of the time parameter, we can iteratively linearize it with a Gauss-Newton optimization using (1) , where x (1) is the initial value obtained from the matched filter approach. A full description of the timedelay estimator is provided in Willemink et al [21].

(ii) Frequency-domain estimation
In tissue due to frequency dependent attenuation, dispersion is present, which causes various phase components of the broadband pulse to travel at different velocities. This could lead in cases to distortion in the signal shapes making a time-domain analysis inaccurate.
In such cases, the phase spectral method may be used which is performed in the frequency domain. Here the FFTs (Fast Fourier Transform) of the sample and water signal are obtained; the phases are subtracted from each other. From this the phase velocities may be calculated as: where   is the phase difference between sample and water signals at a specific frequency  , l is the ray-path defined at the rotation angle  and n is the detector position.
Of the methods outlined above, time-domain maximum likelihood (ML) method has the advantage that all frequencies in the signal participate in the estimation which results in a low variance in the estimated value. The disadvantage is that the bias is higher for certain samples due to dispersive effects which are not modeled. The phase-spectral method while addressing the problem of dispersion does have the disadvantage that the variance in estimation is higher. Only a single frequency component contributes to the estimate and further no ML implementation has been developed. For the phantoms used, we decided to retain the timedomain method. While this is not optimum, it provides a good first approximation to the integrated time delay estimate on the phantoms we used.

(b) Acoustic attenuation estimation
In the estimation of the ultrasound propagation parameters of the object, we also perform a reference measurement in water without the object. By this we can eliminate instrumental unknowns such as the transfer function of the transducer and its spatial dependency due to diffraction effects introduced by its finite aperture size. In terms of transfer functions, the dependence of the object signal and the reference signal can be written as: The reconstructed images of acoustic propagation clearly show the gross features of the hidden inhomogeneities. In the SOS image the 5 mm milk-agar inserts show a SOS of 1505 ms 1 , which is in good agreement with independently characterized values ( Table 1). The subresolution targets show a SOS of 1501 ms 1 with a spread over a slightly larger area than the original dimension (3.2 mm compared with the actual 2.6 mm). This is likely due to the resolution limit of the system, the detection of sub-sampling period TOF delays is strongly influenced by SNR and leads to inaccurate SOS estimation.
In the frequency dependent attenuation image (Fig. 3(c)), the inhomogeneities are faithfully reconstructed with the expected attenuation values (0.4 dBcm 1 MHz 1 ) and dimensions. However, there are artifacts at the boundary of the phantom with the presence of concentric rings. This is most likely due to ray refraction at the boundaries leading to multipath propagation [25]. This occurs when more paths than the line-of-sight path exist between passive element and detector due to refraction. The TOF of multipath signals will have small propagation time differences: the signals interfere so that the assumed model of measuring a single signal is not valid. Therefore the results exhibit artifacts of high attenuation and amplification at the boundaries. Fig. 4. Average values of speed-of-sound tomograms obtained by imaging the phantom (see Fig. 2(b)) at various heights.
We obtained multiple SOS tomograms spaced 1 mm apart in the z direction of the phantom described in Fig. 2(b). Figure 4 is the average SOS value of each image plotted against the corresponding slice number. The SOS values obtained are 1506 ms 1 in the agar-milk layer and 1497 ms 1 in the agar-only layer (See Table 1). As the interface between the two layers is spanned, the relative influence of the layers is seen in intermediate values of SOS. The transition between the two values extends over 5 z-axis positions (5 mm). Thus the slice thickness or elevation plane resolution of the system can be defined as close to 5 mm, and this value is expected from the imaging geometry and the dimension of the individual piezocomposite elements.  The photoacoustic image shows only the outer surface of the muscle tissue; this is due to the high absorption and thus limited light penetration depth at 532 nm. Adipose tissue has low absorption in the green and exhibits a low photoacoustic response. The SOS image resolves the muscle from fatty tissue clearly with values obtained for speed-of-sound (1572 ms 1 and 1525 ms 1 respectively) close to previously measured values. The mixed layer is indistinguishable from the neighboring muscle tissue.

Hybrid imaging
The value of acoustic attenuation is higher in adipose tissue compared to muscle, making the acoustic attenuation image appear complementary to the SOS image. The values are close to the values obtained in the characterization measurements. (see Table 1.) Significantly, the mixed muscle-fat layer is clearly visible in the image due to the higher content of fat in the tissue. This shows the sensitivity of the system and also the relevance of performing hybrid imaging to distinguish the different tissue structures; the mixed layer would otherwise not be uncovered from only photoacoustic or speed-of-sound images.

Summary
In summary, we have demonstrated a "hybrid" imaging approach which can simultaneously image both optical absorption properties and acoustic transmission properties of an object in a two-dimensional slice by adding a "passive" element to a computed tomography photoacoustic imager. The reconstructions are correctly showing the structure of the inhomogeneities. Also the reconstructed speeds of sound and acoustic attenuation values are in good agreement with the actual values obtained from the characterization measurements. The images obtained show the feasibility of our technique to resolve the time of flight difference as low as 12.5 ns with a slice thickness of 5 mm. The obtained photoacoustic image however shows only the outer surface of the muscle tissues, which is due to the wavelength used in the experiment. It is evident that the light penetration depth is limited at 532 nm due to the high optical absorption and scattering in the tissue.
The PER-PACT system provides an intrinsically hybrid tool to image ultrasound transmission parameters without any additional measurement to the conventional photoacoustic imaging protocol. Such an approach permits not only functional information from conventional photoacoustics to be extracted but also anatomic and morphological information from ultrasound parameter depiction. Since the speed of sound in cancerous tissue (1559 ms 1 ) is different from that in surrounding fat (1470 ms 1 ) and glandular tissue (1515 ms 1 ) [14], we believe that this hybrid imaging modality has provided enhanced potential in tumor detection than photoacoustic imaging alone.
The samples used in the experiment possess small SOS differences, for which the ray model used to recover speed of sound and acoustic attenuation is satisfactory. However, for higher speed of sound differences, in future we will use an approximated solution to the wave equation [13] with a spatially varying speed of sound function. This solution incorporates refraction effects of the propagating wave due to the speed of sound inhomogeneities. Further, the obtained speed of sound (SOS) values can be used as a feedback to correct photoacoustic reconstruction in the presence of SOS inhomogeneities thereby improving resolution and contrast [26]. Future work using the PER-PACT approach is in performing hybrid imaging of small animals using near-infrared wavelengths.