1. Introduction

Photoacoustic (or optoacoustic) tomography (PAT) is a non-invasive imaging technique that can provide images with ultrasound resolution related to the distribution of optical absorption within biological tissue. Its increasing significance in medical and biological research arises from the possibility to visualize intrinsic absorbers like hemoglobin. With multispectral PAT also functional information can be extracted such as blood oxygenation. Hence, the field of application is very broad, including preclinical studies of tumor growth, breast imaging, imaging of the microvasculature and small animal imaging [1–3].

In general, PAT relies on the thermoelastic effect, the generation of acoustic waves due to the absorption of pulsed laser light. Local absorption and rapid heating result in generation of broadband ultrasound waves. By recording these waves, the initial absorbed energy distribution can be reconstructed. Thus, photoacoustic imaging is a hybrid technique making use of optical absorption and ultrasonic wave propagation, leading to high contrast and high spatial resolution.

The recording of the broadband, time resolved signals can be done with piezoelectric, optical, or capacitive micromachined ultrasound detectors. Due to the high ultrasound detection sensitivity even for small sensing areas, the uniform amplitude response over a wide frequency range and their optical transparency, optical detection (OD) methods are ideally suited for high resolution, three dimensional PAT [4–10].

Despite the advantages of OD compared to piezoelectric detection, OD will have difficulties entering wide field of applications as long as it is not possible to implement compact optical parallel detection for fast in vivo imaging. However, in the last decade efforts have been made to overcome this drawback by using a charge coupled device (CCD) camera as optical array detector for ultrasound. The CCD array has parallel detection capabilities, but the very limited frame rate precludes recordings of single shot time series. An alternative approach for fast data acquisition with a CCD camera is to take snap shots of the acoustic field and use the spatial information content of a single captured image at a certain time. This information is equivalent to taking time resolved signals recorded at defined detector positions. The proof of principle using a CCD camera in fast three-dimensional (3D) or 2D PAT has been already demonstrated in two different optical arrangements recording either a specific slice or the projection of the acoustic field. While Koestli et. al. [11] and Paltauf et. al. [12] used a camera in combination with a Fresnel reflectance sensor, Lamont et. al. [13] used a camera in combination with a planar polymer film Fabry-Perot interferometer (FPI) sensor to record slices of the acoustic field. The sensitivity of about one percent reflectance change per bar of acoustic pressure of the former robust and simple approach is not sufficient for PAT. The polymer film FPI sensor can achieve in principle an about two orders of magnitude higher sensitivity compared to the reflectance sensor, however the lacking homogeneity of the sensitivity over a large area for a certain interrogation wavelength and film thickness makes the applicability of this method difficult.

Niederhauser et. al. [14,15] used a camera in a Schlieren and a dark field setup to visualize the projections of the acoustic field in real time. We have shown the applicability of an optical phase contrast detection (OPCD) method in PAT [16]. In acoustics, this method has been already used to measure ultrasonic fields [17,18]. However, to our knowledge it had never been used before in PAT. The main advantage compared to Schlieren or dark field techniques is the contrast of the recorded images, which is proportional to the pressure integrated along the light path. Therefore the recorded images are well suited for image reconstruction algorithms.

The information gained with the OPCD method about the optical phase shift is therefore a representation of the pressure field at a given time, projected or integrated along a certain direction. This method is related to the concept of using integrating line detectors in PAT [19–21]. Instead of a single detector that measures time resolved signals from a single line detector (a free or guided light beam), full field detection uses a detector array, interrogating the integrated pressure at a given time from a 2D array of line detectors. To obtain the initial 3D pressure distribution, first of all, the acoustic wave pattern has to be captured from several directions. Each individual pattern is used to calculate a 2D projection of the pressure distribution at time zero. The second step is to apply the inverse Radon transform to the projection data set calculated in the previous step, to obtain the initial 3D pressure distribution. While our earlier work focused on the 2D reconstruction of phantoms from projections of the pressure field [16], we now demonstrate the full 3D reconstruction, provide a thorough characterization of the imaging setup and show first in vivo images of mouse vasculature.

The work is organized as follows: The method section contains information about the working principle of the OPCD setup, the reconstruction method, the data acquisition procedure and the experimental estimation of sensitivity. The results section provides information about the achievable resolution and the applicability of the setup for in-vivo imaging. Finally, the results are discussed and an outlook is given for future developments.

2. Method

2.1 Experimental setup

The OPCD setup for capturing pictures of acoustic wave patterns (the phase object) is shown in Fig. 1. Laser pulses with 10 ns pulse duration coming from the frequency doubled output (532 nm) of a 10 Hz Nd:YAG laser system were used for photoacoustic excitation. All PAT experiments shown in this work were performed in backward mode by illuminating the sample from below with a radiant exposure below the American National Standards Institute (ANSI) limit of 20 mJ/cm². The sample was arranged in the center of the sample holder (petri dish) filled with water or ultrasound gel. The sample holder was mounted on a rotary table and only slightly dipped in water for acoustic coupling and positioned in the rear focal plane (object plane) of the Fourier transforming lens (L1, f₁ = 200 mm). The pulsed probe laser beam coming from a diode pumped solid state laser system (wavelength 527 nm, pulse duration 8 ns) was collimated and expanded to a beam diameter of about 35 mm, matched to the size of the field of view (FOV) of the detection system. This ensured that the detection aperture was sufficiently large to capture the outgoing acoustic waves coming from the optically absorbing regions inside the sample and that the total area of the CCD element was used.

 

Fig. 1 Phase contrast detection setup to capture projection images of acoustic fields. PP: partially absorbing phase plate, OW, optical window, NBPF: narrow bandpass filter, L: lens, M: mirror

Download Full Size | PDF

The working principle of the detection system is based on the interaction of the acoustic and optical fields (elasto-optic interaction). The latter attains a phase variation proportional to the pressure integrated along the probe beam path. This phase variation is converted into a measurable intensity modulation by arranging an optical signal processing element (PP: phase plate) in the Fourier plane of L1. Like in optical phase contrast microscopy an additional relative phase shift of either plus or minus π/2 between the non-diffracted and diffracted portions of the probe laser beam leads to constructive or destructive interference between both parts. This results in a so called positive or negative phase contrast image. The lens L2 (f₂ = 90 mm) is arranged in a way that its rear focal plane and the phase plate coincide, producing a reversed image of the phase object on the CCD element of the camera. A narrow bandpass filter (NBPF) is placed directly in front of the camera to avoid exposure with stray light from the excitation laser pulse. Snapshots of the pressure distribution are taken with a CCD camera (pco.2000s, 14 bit dynamic range) at certain time delay with respect to the excitation laser pulse. The time delay can be adjusted by using a delay generator.

The pressure induced optical phase variation is given by the integral of the pressure field over the probe beam path in y direction [17],

In principle, the resolution is adjustable with MO. The natural limit is the optical diffraction limit determined by the wavelength and the numerical aperture of imaging system. Note that the acoustic field is integrated along the propagation direction of the probe laser beam. Therefore, also the depth of field (~8 mm) of the optical imaging system is an important parameter.

2.2 Reconstruction

A recorded camera image reveals a map of the acoustic field around the object, integrated or projected in direction of the detection beam for a certain orientation $p_{t=T_0}^{\alpha }\left(x,z\right)$ of the sample. A very fast and easy to implement reconstruction algorithm is the Fourier algorithm. This algorithm is used to obtain projection images by back propagating the waves to time zero$p_{t=0}^{\alpha }\left(x,z\right)$, when they were produced by pulsed laser excitation. The 2D Fourier reconstruction algorithm is given by

2.3 Data acquisition method

The measurement system is sensitive to any phase object along the optical path of the probe beam. Besides the desired phase object, the acoustic field, also undesired phase objects like impurities on and in optical components, temperature fluctuations and floating particles in water appear on recorded images. Hence, it is almost impossible to record an image with homogenous background. To solve this problem a “background image” (BG, image without wave pattern) and a “foreground image” (FG, image with wave pattern) have to be recorded for each orientation α of the sample (see Fig. 2(a) and Fig. 2(b)). The subsequent subtraction of those images yields wave pattern images with almost constant background. To improve the contrast to noise ratio (CNR), this procedure can be repeated several times for averaging. In the final preprocessed image (Fig. 2(c)) the background pixel values fluctuate around zero due to the noise level of the CCD camera. On closer examination of Fig. 2(c) it can be seen that the time varying part of the background, like floating particles (bright spots in Fig. 2(a) and Fig. 2(b)), is not completely removed due to the movement of the particles during the time delay of 0.1s (determined by the frame rate of the camera) between the recordings. The stripe with negative pixel values along x direction at position z = 7.5 mm also results from a time varying background that could not be removed by the subtraction method. The reconstructed projection image is shown in Fig. 2(d). Spurious bright spots and stripes in the preprocessed image cause arc and line shaped artifacts in the reconstructed projection image due to the divergence after reconstruction. This underlines the necessity of preventing artifacts by using optically clean water as acoustic propagation medium and by avoiding turbulences in water so that the background subtraction method can work efficiently.

 

Fig. 2 Representative recorded background image (a) and foreground image (b). Four times averaged wave pattern image after background subtraction (c). Reconstructed projection image (d) from the preprocessed wave pattern image (c).

Download Full Size | PDF

For both, the BG and the FG image, the sample is irradiated with the excitation laser pulse. The difference is the chosen delay time between excitation and probe laser pulse. The delay time, set by the delay generator, is switched for each exposure between the values t = 0 µs and t = T₀ > 0µs. The traveling time T₀ should be smaller than 10 µs, which is the time acoustic waves need to cross half of the FOV with its size of about 30 mm. The traveling time and the resulting traveling distance c·T₀ are directly linked to the aperture γ of the visible wave pattern from a point source. If γ approaches 180°, a half circular pattern is visible in the image, and the horizontal resolution is comparable to the vertical resolution. However, smaller angles will lead to a deterioration of lateral resolution (limited view problem). Since such large detection angles are only possible for superficial structures a compromise has to be found concerning optimum lateral resolution and imaging depth.

2.4 Sensitivity and dynamic range

The sensitivity of the OPCD system was estimated using a well-defined photoacoustic source. This source was generated by coupling the excitation laser pulses into a 1 mm glass fiber. The flat fiber exit was submerged in a cuvette filled with an absorbing dye solution (OrangeG/water mixture). Knowing the excitation pulse energy, the absorption coefficient of the dye solution and the shape of the acoustic source the pressure field could be simulated for a certain propagation distance. For a distance of 15 mm, matched to half the size of the FOV, the simulated pressure fields were integrated along y-direction, yielding maximum pressure-length-product-values $p\cdot L_{\max }$ for six different laser pulse energies used in the experiment. By plotting the corresponding CNR values from the recorded preprocessed wave pattern images against the $p\cdot L_{\max }$values and applying a linear regression to the data points, the noise equivalent pressure length product (NEPLP) of 5.1 kPa mm could be obtained from the inverse of the slope k of the regression line. The CNR values of the preprocessed images were estimated using the relation:

 

Fig. 3 Thirty-two times averaged preprocessed wave pattern image originating from a defined acoustic source (a). Profile taken from the wave pattern image along z direction at x = 0 position (b). Linear dependency of the contrast to noise ratio (CNR) on the pressure length product (c). The reciprocal value of the slope corresponds to the noise equivalent pressure length product (NEPLP).

Download Full Size | PDF

The theoretically achievable sensitivity of the imaging system can be derived via the CNR of the CCD camera. It is given by the ratio of generated signal charge carriers (signal electrons) to the number of unwanted charge carriers (noise electrons).

The upper limit of the dynamic range can be estimated from the requirement of a linear relation between the phase shift caused by the integrated pressure and the intensity modulation. Assuming that the linearity holds up to a phase shift of π/10, the related pressure value is 200 kPa (from Eq. (1)), integrated over a length of 1 mm. With the NEPLP from above, this gives a theoretical dynamic range of 56 dB.

3. Results

3.1 Phantom imaging experiment

To determine the spatial resolution and to test whether the proposed method is suitable for fast imaging an experiment was conducted using a phantom sample containing black human hair loops with diameters ranging from 55 to 70 µm and black polystyrene microspheres with 110 µm diameter embedded in a mixture of 2% agar and 5% intralipid in water. The sample was arranged in the setup as described in section 2.1. The chosen radiant exposure was 12 mJ/cm², about two times below the ANSI limit for 10 ns laser pulses with a wavelength of 532 nm. A delay time of T₀ = 7 µs was chosen, corresponding to a wave propagation distance of 10.5 mm, to maximize the detection aperture.

Before reconstruction, the recorded wave pattern images were preprocessed as described in 2.3. As the number of averages was set to four, we recorded eight images (four FG and four BG images) for one projection image, requiring less than one second with the 10 Hz pulse repetition and camera frame rate. To obtain the 3D image, two hundred averaged projection images were recorded over a half rotation. Consequently, the recording of the 3D data took just about 3 minutes. From this experimental result (Fig. 4) we estimated the vertical/axial (z-direction) and horizontal/lateral resolution (xy-planes) to be better than 73 µm and 80 µm, respectively.

 

Fig. 4 Photoacoustic maximum amplitude projection images of the 3D-image of the phantom sample along three different directions (a), (b) and (c). Horizontal and vertical profiles taken from (a) and (b) are shown in (d), where pixel values along the line HP are displayed and in (e) for values along line VP. The full width half maximum (FWHM) values of the peaks, labeled by the black arrows, are representative for the horizontal and vertical resolution.

Download Full Size | PDF

3.2 In-vivo imaging experiment of mouse vasculature

The experiment started with the preparation of a 14 week old female CD1 mouse. The mouse was anesthetized by a single intraperitoneal injection of ketamine (70 mg/kg; Pharmacia & Upjohn, Germany) and xylazin (12 mg/kg; Bayer, Leverkusen, Germany) to avoid unwanted extensive movement of the mouse during the experiment over a time period of about 30 minutes to avoid motion artifacts. The musculus Biceps femoris was defined as the region of interest (ROI) and all hair was removed from the hindlimb using depilation cream prior to scanning. Then the mouse was carefully placed in an upright position in the petri-dish filled with water and gently fixed to the sample holder with a soft elastic tape. The sample holder was mounted to the rotary table and slightly tipped in water for acoustic coupling (see Fig. 1).

The chosen radiant exposure was 15 mJ/cm², which is still below the ANSI limit. The used excitation wavelength in the visible range (532 nm) limited the imaging depth to a few millimeters. The wave pattern images were recorded at delay time T₀ = 9 µs and due to eight times averaging the recording time of each projection image shown in Fig. 5 was 1.6 seconds. The projection images can be obtained in principle in real-time (without averaging) and reveal already the complex vessel structure. However, a single projection image does not directly show the absorbed energy density. For instance, a bright region in a photoacoustic projection image (few of them are indicated by arrows in Fig. 5) can either mean an enhanced local light absorption or that an absorbing structure is oriented along the projection direction. In applications where the projection images are sufficient, like the monitoring of biopsies, it is a fast alternative to full 3D imaging.

 

Fig. 5 Photoacoustic projection images from four different sample orientations α of the left hind leg of a mouse. The peculiarity of photoacoustic projection images is that blood vessels (extended absorbing structures) oriented in projection direction appear as bright spots (white arrows).

Download Full Size | PDF

The 3D absorbed energy distribution in the hind leg of the mouse was obtained by recording one hundred projection images over a half rotation. Consequently, the recording time of the 3D data took again 3 minutes. In Fig. 6 maximum amplitude projection images are depicted from the top and the side view showing the vasculature network in the muscle of the hind leg of the mouse. The appearance of the blood vessels visible in the zoom image (Fig. 6(c)) clearly indicates a lateral resolution below 100 µm.

 

Fig. 6 Photoacoustic maximum amplitude projection (MAP) images from side and top view of the left hind leg of the imaged mouse in (a) and (b). A zoom of small blood vessels in the dashed box in (b) is shown in (c) indicating the image resolution <100 µm. Scale bar: 300 µm.

Download Full Size | PDF

Approval to this study by the Animal Protocol Review Board of the City of Vienna was obtained prior to conducting experiments. They were carried out in accordance with the guideline for “Care and Use of Laboratory Animals” of the National Institute of Health.

3. Discussion

The presented photoacoustic tomography technique takes advantage of the high detector density and number offered by a CCD array. We could demonstrate that it can generate 3D photoacoustic images with high resolution in a short time. Its working principle is similar to a high density line detector array, but contrary to an arrangement of single ultrasound detectors it captures snapshots of the acoustic wave pattern at given delay time after the excitation laser pulse, rather than time dependent pressure signals at the predefined detector positions. Consequently, the information on the acoustic field depends on the field of view of the camera and on the chosen delay time between excitation and probe laser pulses.

In the setup the sample was separated from the coupling liquid by a thin plastic film (base of the petri-dish). This ensured that no impurities of the sample penetrated the optical imaging volume. Other than in a previous study where the full acoustic field was detected after leaving the object [16], in this setup only the acoustic waves entering the half space adjacent to the sample surface are captured. This gives rise to some limited data artifacts and reduced lateral resolution with increasing imaging depth and normal distance to the rotation axis but more closely corresponds to a real imaging situation where only one side of an object is accessible. This effect can be clearly seen in the MAP image shown in Fig. 6(b). Blood vessels in the central region are clearly resolved but those located towards the periphery of the imaging area become increasingly blurred.

3D imaging with the integrating detection approach needs to record projection images from different orientations of the sample. Hence, the detection aperture varies for off-axis absorbing structures due to varying normal distances to the rotation axis while rotating the sample. This circumstance improves the off-axis 3D image resolution compared to the resolution in a single projection image since all projection images from the different orientations of the sample contribute to the 3D image. This is an advantage of the integrating line detection concept compared to approaches where a point detector scans over a plane [4,7–9], where the detection aperture of a certain reconstructed image voxel is fixed. The estimated horizontal and vertical 3D-image resolution from the phantom experiment (see Fig. 4) in the central region of the imaging volume is better than 80 µm and 73 µm, respectively, in accordance with the theoretical limit of 66 µm. Currently, only the PAT system based on the thin film Fabry-Perot sensor achieves a comparable resolution <100µm [4].

The imaging time of the proposed OPCD tomograph depends on the operating mode. The system can operate in the projection or 3D imaging mode. The projection mode without averaging allows recording of projection images in real-time. The currently achievable frame rate for projection imaging is 10Hz, limited by the CCD-camera and the repetition rate of the pulse laser system. This practically real-time mode could be used for applications where a full 3D image is not needed, such as guided needle biopsy. The 3D imaging mode requires the rotation of the sample while recoding the projection images. For a required resolution Λ_min within a section of the image with side length SL, the angular step size should not exceed $\Delta {\alpha }^{\max }=90°\cdot {\Lambda }_{\min }/SL$, derived from the Fourier-Slice-Theorem and the Nyquist-criterium [23]. Assuming a resolution of 100 µm and a cubic imaging volume with side length of 10 mm the maximum angular step size is 0.9°. Hence, the data acquisition for 3D imaging without averaging takes less than 1 minute with the 10 Hz frame rate considering that 200 projection images must be recorded over an angular range of 180°.

The experimentally and theoretically estimated sensitivity, the NEPLP value, is 5.1 kPa mm and 0.31 kPa mm, respectively. The large discrepancy of these values is a clear hint that still many things can be optimized. One is the improvement of the background subtraction method by measuring the probe pulse energy and using it to normalize the FG and BG image before image subtraction, the improvement of the optical phase plate by adding an anti-reflection coating avoiding ghost images and the reduction of unwanted phase objects (floating particles) that cause a limitation of the dynamic range due to the saturation of the CCD-camera.

The experimentally estimated NEPLP value of 5.1 kPa mm for acoustic wave monitoring is about 10 times larger (less sensitive) compared to the free beam optical line detection system (0.5 kPa mm) based on two beam interference (Mach-Zehnder interferometer) [24]. Compared to a 1 mm diameter piezoelectric PVDF (28µm) sensor terminated with a high quality preamplifier (0.66 kPa) [25] the NEPLP is about 8 times worse. Despite the somewhat poorer performance of the sensitivity of the OPCD method in comparison with other detection methods, high quality images (see Fig. 4 and Fig. 6) can be obtained without exceeding the permissible exposure levels. The reason is the large number of detection positions (pixels) contributing to the reconstruction compared to other detection approaches. For instance, the wave pattern of a point like source recorded 7 µs after the photoacoutic excitation appears as a semicircle in the recorded image. Dividing the arc length of the semicircle (~33 mm) by the pixel size in the object plane (~33 µm) the number of pixels (detection positions) is 1000 and all these pixels with corresponding pixel value contribute to one specific voxel in the reconstructed projection image improving the CNR in the final image.

4. Conclusion

The presented technique provides real time 2D and fast 3D photoacoustic imaging using a purely optical, parallel detection method working in backward imaging mode. It combines the advantages of optical detection with the possibility to acquire the required high amount of data within a relatively short time using parallel detection with a CCD-camera. Next steps will be multispectral in-vivo imaging to answer scientific questions from biological and medical research.