Full field detection in photoacoustic tomography

Imaging the full acoustic field around an object by use of an optical phase contrast method is used to accelerate the data acquisition in photoacoustic tomography. Images obtained with a CCD-camera at a certain time show a projection of the instantaneous pressure field in a given direction. In this work a reconstruction method is presented to obtain the two-dimensional initial pressure distribution by back propagating the observed wave pattern in Radon space. Numerical simulations are used to prove the accuracy of the reconstruction algorithm and to demonstrate a method for correcting limited data artifacts. Finally, the overall performance is shown with experimentally obtained data. 2010 Optical Society of America OCIS codes: (170.5120) Photoacoustic imaging; (100.3190) Inverse problems. References and links 1. M. H. Xu and L. V. Wang, "Photoacoustic Imaging in Biomedicine," Review of Scientific Instruments 77, 041101 (2006). 2. 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," Journal of Biomedical Optics 14, (2009). 3. E. Z. Zhang, J. G. Laufer, R. B. Pedley, and P. C. Beard, "In vivo high-resolution 3D photoacoustic imaging of superficial vascular anatomy," Phys. Med. Biol. 54, 1035-1046 (2009). 4. Paul C. Beard, Andrew M. Hurrell, and Tim N. Mills, "Characterization of a polymer film optical fiber hydrophone for use in the range 1 to 20 MHz: A comparison with PVDF needle and membrane hydrophones," IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 47, 256-264 (2000). 5. P. C. Beard and T. N. Mills, "Extrinsic optical-fiber ultrasound sensor using a thin polymer film as a lowfinesse Fabry-Perot interferometer," Appl. Opt. 35, 663-675 (1996). 6. S. Ashkenazi, Y. Hou, S. Huang, T. Buma, and M. O'Donnell, "High Frequency Optoacoustic Transducers for Ultrasonic and Photoacoustic Imaging" (CRC Press, 2009). 7. J. J. Niederhauser, M. Jaeger, and M. Frenz, "Real-Time Three-Dimensional Optoacoustic Imaging Using an Acoustic Lens System," Applied Physics Letters 85, 846-848 (2004). 8. J. J. Niederhauser, D. Frauchiger, H. P. Weber, and M. Frenz, "Real-Time Optoacoustic Imaging Using a Schlieren Transducer," Applied Physics Letters 81, 571-573 (2002). 9. B. Schneider and Shung K. K., "Quantitative Analysis of Pulsed Ultrasonic Beam Patterns Using a Schlieren System," IEEE Trans. Ultrason. Ferroelect. Freq. Control 43, 1181-1186 (1996). 10. N. Kudo, H. Ouchi, K. Yamamoto, and H. Sekimizu, "A simple Schlieren system for visualizing a sound field of pulsed ultrasound," J. Phys.: Conf. Ser. 1, 146-149 (2004). 11. C. I. Zanelli and S. M. Howard, "Schlieren metrology for hight frequency medical ultrasound," Ultrasonics 44, e105-e107 (2006). 12. I. Nunez and J. A. Ferrari, "Bright versus dark schlieren imaging: quatitative analysis of quasi-sinusoidal phase objects," Applied Optics 46, 725 (2007). 13. G. S. Settles, Schlieren and Shadowgraph Techniques (Springer, 2001). 14. Th. Neumann and H. Ermert, "Schlieren visualization of ultrasonic wave fields with high spatial resolution," Ultrasonics 44, e1561-e1566 (2006). 15. E. K. Reichel and B. G. Zagar, "Phase contrast method for measuring ultrasonic fields," IEEE Transactions on Instrumentation and Measurement 55, 1356-1361 (2006). 16. P. Burgholzer, C. Hofer, G. Paltauf, M. Haltmeier, and O. Scherzer, "Thermoacoustic tomography with integrating area and line detectors," IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 52, 1577-1583 (2005). 17. M. Haltmeier, O. Scherzer, P. Burgholzer, R. Nuster, and G. Paltauf, "Thermoacoustic tomography & the circular Radon transform: Exact inversion formula," Mathematical Models and Methods in Applied Sciences 17, 635-655 (2007). 18. G. Paltauf, R. Nuster, M Haltmeier, and P. Burgholzer, "Photoacoustic tomography using a Mach-Zehnder interferometer as acoustic line detector," Appl. Opt. 46, 3352-3358 (2007). 19. S. Helgason, The Radon Transform (Birkhäuser, Boston 1999). 20. F. John, Partial Differential Equautions (Springer Verlag, New York 1982). 21. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, New York 1988). 22. F. Natterer, The Mathematics of Computerized Tomography (SIAM Classics in Applied Mathematics, 2001). 23. X. Pan, E. Y. Sidky, and M. Vannier, "Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction?," Inverse Problems 25, 123009 (2009). 24. M. Haltmeier, O. Scherzer, P. Burgholzer, and G. Paltauf, "Thermoacoustic computed tomography with large planar receivers," Inverse Problems 20, 1663-1673 (2004). 25. G. Paltauf, R. Nuster, M Haltmeier, and P. Burgholzer, "Experimental evaluation of reconstruction algorithms for limited view photoacoustic tomography with line detectors," Inverse Problems 23, S81-S94 (2007). 26. G. Paltauf, R. Nuster, and P. Burgholzer, "Weight factors for limited angle photoacoustic tomography," Phys. Med. Biol. 54, 3303-3314 (2009).


Introduction
Photoacoustic tomography (PAT) is a hybrid biomedical imaging technique combining optical contrast with ultrasound resolution.Acoustic waves are excited via the thermoelastic effect.By illuminating a sample with nanosecond laser pulses, each light absorbing inclusion acts as an individual acoustic source.The distribution of the initial pressure in these sources is determined by the optical properties of the sample.After propagating from their initial source position towards the surface the acoustic waves are detected at defined detector positions around the sample.Subsequently the recorded time resolved signals are used to reconstruct the initial three dimensional pressure distribution with a resolution that is determined by the acoustic part of the imaging procedure [1][2][3].
Piezoelectric detectors are commonly used to record the acoustic waves due to their high sensitivity, their flexibility to manufacture arbitrary shapes, and their applicability to parallel detection.However, in the last decade efforts have been spent to develop optical ultrasound detectors, which are preferable in some applications because of their uniform amplitude response over a wide frequency range and their optical transparency for the excitation laser pulse.Furthermore, in applications where the active sensing area should be as small as possible ("point like") optical sensors based on interferometry are more sensitive than piezoelectric sensors [4][5][6].
Despite their advantages, PAT devices using optical detection will have difficulties entering commercial applications as long as it is not possible to realize compact optical parallel detection.The main problem is the need for some kind of array detector.Charge coupled device (CCD) cameras have parallel detection capabilities but at very limited camera frame rate.For real-time detection in PAT the desired frame rate would have to be in the range of megahertz, which is not possible at reasonable costs.An alternative approach to speed up the data acquisition is to use the spatial information content of a single captured image at a certain time instead of using time resolved signals recorded at defined detector positions.Under certain conditions a single snapshot of the acoustic wave pattern outside an object contains all information to reconstruct a two-dimensional (2D) projection of the initial pressure distribution inside the object.The proof of principle for use of a CCD camera in real-time PAT has already been demonstrated [7,8].
There are different ways for visualizing the acoustic wave pattern and for the reconstruction of the initial pressure distribution.Commonly the Schlieren or dark field method is used to visualize acoustic waves or similar optical phase objects.The fields of application vary from the observation of flow dynamics in liquids or gases to the estimation of the focusing properties of ultrasound transducers [9][10][11][12][13][14].However, as these methods do not deliver a contrast proportional to the phase variations they are not ideally suited for PAT.As an alternative, if the pressure induced phase shift of the probe laser beam is much smaller than π, interferometry delivers a linear relation between the detected light intensity and the pressure values.This relationship is based on a linear dependence of the change of optical refractive index on acoustic pressure.The phase contrast technique, developed by Fritz Zernike, uses the interference between diffracted and non-diffracted light in an optical imaging setup and also delivers a contrast in the captured images which is proportional to the pressure amplitude [13,15].All the imaging methods described above use illumination of a clear coupling liquid surrounding a sample with a collimated light beam.The gained information about the optical phase shift in the coupling liquid 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 [16][17][18].Instead of a single detector that measures time resolved signals from a single line detector (a free or guided light beam), full field detection (FFD) uses a detector array, interrogating the integrated pressure at a given time from a 2D array of line detectors.To obtain the initial three dimensional (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 (RT) to the projection dataset calculated in the previous step, to obtain the initial 3D pressure distribution.
The projection images of the initial pressure distribution can be obtained either by convolving or deconvolving the captured acoustic wave pattern with the Green function of the 2D wave equation.These methods, however, yield a limited accuracy of the resulting reconstructed images.The convolution is equivalent to simultaneous forward and backward propagation in time of the 2D acoustic wave pattern, leading to distributions which are not completely spatially separated.The deconvolution on the other hand is very sensitive to noise.Furthermore both methods are not ideally suited to obtain accurate reconstruction results from limited data, when parts of the pressure distribution are missing.
In this work a reconstruction method is presented to obtain the 2D initial pressure distribution from projection images of the pressure field.Numerical examples are shown to point out the accuracy of the proposed reconstruction method.Next, a strategy is presented for improving the image quality in the limited data case, and finally the overall performance is shown with experimentally obtained data.

The formal problem in full field detection
In PAT, the propagation of an acoustic wave, originating at the absorbing inclusions in an object after illumination by a short laser pulse, is governed by the wave equation .0 0 Here c is the speed of sound and ) ( p x 0 denotes the initial pressure distribution.In the FFD, linear projections of the pressure T) p( , x , for some fixed time T , are recorded.The measurement time T has to be chosen large enough, such that essential parts of the acoustic wave have left the domain { } We denote by σ E the plane that is orthogonal to the direction ( ) , q q = q denoting the coordinates in the plane σ E and s denoting the coordinate orthogonal to σ E .Finally, we denote by

(
) As illustrated in the left image of Fig. 1, the X-ray transform is the linear projection of ) (x h onto the imaging plane σ E , obtained by integrating that function along lines orthogonal to σ E and passing through the point ) 0 , (q x .
In the FFD, pictures of the instantaneous projections of the pressure field at a fixed positive time 0 > T are captured, given by The imaging problem in FFD consists in reconstructing the initial pressure ) ( p x 0 from the observed data ) , ( T P q σ .In the following section, we decompose this imaging problem into several 2D problems, and present an explicit solution method for the resulting 2D problems.

Reconstruction process
The X-ray transform and the Laplace operator satisfy the commutation relation [19] ( denoting the Laplace operator in two and three spatial dimensions, respectively.Because the Laplace operator is invariant with respect to rotations around the ) 1 , 0 , 0 ( axis, Eq. ( 1) also holds for the new coordinates ) , , ( This implies, that the linear projections The FFD provides data ) , ( T P q σ for a fixed time T .Below we will demonstrate how to recover the projection of initial pressure ) ( 0 , q σ P in (3) from such data.Having found the projection images ) )( (X : ) ( 0 0 , q q σ σ p P = the 3D initial pressure distribution ) ( p x 0 can be reconstructed by inverting the X-ray transform, which is well known from classical X-ray computed tomography. In the following we assume that the direction σ is fixed, and denote by from the origin (see the right image in Fig. 1).
The Radon transform and the Laplace operator satisfy a commutation relation, similar to the one in Eq. ( 2), see [19].Therefore, we can argue as before, and conclude that ) , , )( ( Equation (5)

Numerical examples
All numerical calculations were performed on a 200 200× grid, using the MATLAB programming environment with the embedded imaging toolbox.First, to show the performance of the developed reconstruction method the acoustic wave pattern at time T=8 µs was computed (see Fig. 2b) from a structure consisting of five squares as shown in Figure 2a.For this purpose we used the solution formula of the 2D wave equation [20], and assumed the speed of sound in water with 1.5 mm/µs.Figure 2b shows an image which would be captured under ideal conditions in the experiment.Applying the Radon transform to the acoustic wave pattern results in an image consisting of two separated structures (bands), as mathematically outlined in Eq. (5).Each column of this Radon transform denotes plane waves counterpropagating in direction ) sin , cos ( φ φ − , separated by twice the propagation distance cT (black arrow in Fig. 2c).Equation ( 6) performs the sum of the two bands after left and the right translation by an amount of cT and yields three bands from which only the one in the middle is significant, representing the Radon transform of the initial 2D source.This is shown in Fig. 2d where, according to Eq. ( 7), the two bands localized around cT 2 ± have already been removed.Finally, by applying the inverse Radon transform, the initial source is obtained (Fig. 2e).Fig. 2f shows a comparison of horizontal and vertical profiles through the center of the initial source and the reconstruction.Both, the amplitude and the shape of the profiles taken from the reconstructed image and from the original source fit perfectly together, proving the accuracy of the reconstruction method.
In order to test the stability of the outlined reconstruction method against noise we added 20% noise (RMS relative to the image maximum) to the computed acoustic wave pattern.The noisy wave pattern and the resulting reconstruction are shown in Fig. 3. Comparing these two images it can be clearly seen that the reconstruction algorithm increases the signal to noise ratio, and consequently enhances the contrast.This is mainly due to the fact that the Radon transform involves an integration that tends to smear out the statistical fluctuations of the noise.
Equation ( 4) and Fig. 2(c) show that there is redundancy in the Radon transform of the measured data, because each of the two counter propagating solutions (or bands in the Radon transform) contains the same information.Only one band would be sufficient to recover the initial source, offering the possibility to improve the reconstruction when only limited data are available.For instance, in the experiment this is the case when the sample holder blocks a certain region of the imaging field of view.In Fig. 4a the Radon transform is shown when the acoustic wave pattern is available only outside the region indicated by the dashed line in Fig 2b .Mainly the lower band is distorted by the missing data.Consequently, removing the lower band by setting the pixel values of the lower half space, corresponding to 0 > d to zero and weighting the image with a factor of two leads to a significant improvement of the reconstructed image.Using limited data without any corrections the reconstruction of the squares appears distorted, as shown Fig. 4c.Applying the correction greatly improves the reconstruction (see Fig. 4d).This is shown in more detail in Fig. 4b, by comparing the horizontal and the vertical profiles crossing the square localized in the center.
Also mentionable in this context is that the Radon transform of a fully accessible ideal 2D acoustic wave pattern contains only nonnegative pixel values (see Fig. 2c) representing the two counter-propagating plane waves.This is not the case in the realistic situation, where data are missing in the captured image (see Fig 4a).Consequently, negative pixel values appear after applying the Radon transform and cause artifacts in the reconstructed images.These artifacts are significantly reduced by applying the correction method described above.

Experimental details
The experimental setup for capturing a picture of the acoustic wave pattern (the phase object) is shown in Fig. 5.To excite the initial acoustic source, the sample was illuminated from three directions.Laser pulses with 10 ns duration and with variable wavelength were generated by an optical parametric oscillator (OPO) pumped by a frequency tripled Nd:YAG laser system.
The OPO pulses illuminated the sample from two sides.To achieve a reasonably homogenous light distribution in the sample, an additional illumination beam came from the front side.Due to limitations of the OPO pulse energy from our system this beam was taken from the frequency doubled output (532 nm) of the Nd:YAG pump laser.The sample itself was positioned in the water tank, in the center of the probe beam and in the rear focal plane (object plane) of the Fourier transforming lens (L2).Two lenses (Obj., L1) and a pinhole were used as a telescopic beam expander.The probe beam coming from a continuous, diode pumped solid state laser with a wavelength of 532 nm (DPSS laser) was collimated and expanded to a beam diameter which was at least three times the size of the sample.This ensured that the outgoing acoustic waves, although the central part of the image was obscured by the sample, could be fully captured in the clear coupling medium (water) outside the sample.Since it is not possible to visualize the pressure distribution inside a light scattering sample, the factor of three is necessary for letting the waves generated in all parts of the sample propagate into the coupling liquid.
The interaction of the ultrasonic field with the collimated probe beam causes diffraction of the latter and consequently higher spatial frequency components in the Fourier plane of L2 (backside focal plane).A phase plate arranged at the position of this plane performs an optical filtering process.For phase contrast imaging it should introduce a phase shift of a quarter of the optical wavelength between the diffracted light and the central spot, containing the nondiffracted light.These components are then able to interfere with each other in the image plane, giving rise to an intensity modulation proportional to the accumulated phase shift [13,15].The lens L3 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.Snapshots of the pressure distribution were taken with an image intensified, time gated CCD camera (PCO DiCAM, 8 Bit dynamic range).With the control unit the exposure time and the time delay with respect to the laser pulse could be adjusted.The inset in Fig. 6b shows the sample which was used in the experiment.It was a thread of black plastisol with a knot and was attached to a plastic frame.The total radiant exposure of the sample was 100 mJ/cm² from all sides as mentioned before.The exposure time of the camera set to 40 ns.Hence, the propagation distance of the sound wave in water (speed of sound: 1.5 mm/µs) within this time limited the spatial resolution to 60 µm.Ten times averaged snapshots of the acoustic wave pattern were taken 3 µs after the excitation laser pulse.A second image was taken at a delay time at which the acoustic waves had already left the field of view and was subtracted from the first image.This removed most of the disturbing background in the image caused by the impurities of the non ideal optical components along the optical path and by the phase change caused by heat flow around the sample.What remained was an image showing the acoustic wave pattern (Fig. 6a).Finally, the linear projection of the initial pressure distribution was reconstructed, by applying the proposed Radon reconstruction method (Fig. 6b).Apart from a small region which was not illuminated and is therefore not visible, the reconstructed image shows quite well the shape of the phantom.

Discussion
A phase contrast optical imaging method providing full field data of a photoacoustically generated pressure distribution and a corresponding image reconstruction algorithm have been proposed.The method is based on capturing projection images of the acoustic wave that has propagated from the imaging object into a surrounding, clear coupling liquid.The reconstruction yields a 2D image of a projection of the initial pressure distribution in a selected direction.Since only two laser pulses are necessary to acquire the data for the projection image, this method has real time imaging capability.Provided the imaging data are complete, the proposed reconstruction algorithm is exact.It uses a reduction to a 1D problem, in a similar way as the previously proposed photoacoustic tomography method using large, planar detectors [24].In the current study we restricted ourselves to the 2D imaging problem, however, the extension to a 3D imaging procedure only requires additional rotation of the sample during data acquisition and an inversion of the standard Radon transform in the reconstruction.Considering the short time to obtain imaging data for one projection, it should take about one minute or less to gain information for a 3D image, only limited by the laser repetition rate and the frame rate of the CCD-camera.
The reconstruction step involving the Radon transform (Fig. 2c) makes the redundancy in the imaging information visible, which can be used to correct for artifacts that arise due to incomplete data.This is an important property of the imaging method, since in practice it is difficult to obtain an image of the complete acoustic field around a sample.Interestingly, it is the omission of data (zero/one weighting) which clearly improves the result of the reconstruction.A similar positive effect of omission or weighting of redundant data has recently been demonstrated for photoacoustic tomography with point or line detectors [25,26].
The phase contrast technique for acquiring the pressure field images combines several advantages.First, it gives rise to an imaging contrast proportional to the relevant quantity, the projected pressure variations.Furthermore, although it is an interferometric technique, it is inherently stable because it uses properties (separation of diffracted and undiffracted portions) of the same illumination beam and therefore does not require any stabilization as other interferometric methods.A disadvantage seems to be the limited sensitivity.In the current study we used a quite high radiant exposure and also a high contrast imaging object.However, several improvements of the setup can be made to increase the sensitivity.One obvious improvement is to use a camera with higher dynamic range.The camera used in this study can resolve intensity modulations in the range of 0.5 % due to its resolution of 8 bit.In further studies, a camera with 12 bit resolution will be used.Cooling the CCD chip and using a pulsed light source instead of the gated image intensifier will improve the signal to noise ratio and enable a higher dynamic range.Finally, the phase contrast method can be optimized by use of a partly absorbing phase plate [13,15].
In conclusion, a concept for photoacoustic tomography using optical array detection of acoustic waves has been demonstrated.With the envisaged improvements of the sensitivity, fast 3D imaging will be possible.

Fig. 1 .E
Fig.1.Graphical representation of the reconstruction procedure.FFD measures linear projections of T) p( , x onto planes σ E (left).The Radon transform of the projection image reduces the 2D wave propagation to sets of two counter-propagating plane waves (right).
a ball with radius r around the origin in3  R .In such a situation, the Radon transform ) principle for the 1D wave equation (which immediately follows from D'Alemberts formula (5)) implies that ) separated, and Eq.(6) implies the following explicit inversion formula, ( ) the inverse Radon transform in each of these planes gives the desired 3D reconstruction of )

Fig. 2 .
Fig. 2. Numerical examples proving the performance of the reconstruction method.(a) Mathematical phantom representing the projection of the initial pressure distribution.(b) Computed acoustic wave pattern at time T=8 µs.(c)-(d) Radon transform of the acoustic wave pattern and the reconstruction procedure in Radon space.(e) Reconstruction of the initial 2d source.(f) Horizontal and vertical profiles through the center of images (a) and (e) .

Fig. 4 .
Fig. 4. Numerical examples for limited data reconstruction.(a) The effect of missing data on the Radon transformed image.(c) and (d) Reconstruction from limited data with and without proposed correction method.(b) Comparison of horizontal and vertical profiles.Red line: without correction, blue line: with correction, black dotted line: original source distribution.

Fig. 5 .
Fig. 5. Experimental setup to capture a snap shot of the acoustic wave pattern at time T. M: mirror, BS: beam splitter, PH: pin hole.

Fig. 6 .
Fig. 6.(a) Captured acoustic wave pattern.(b) Reconstructed projection of the initial pressure distribution.The inset shows the plastisol phantom.