Photoacoustic clutter reduction by inversion of a linear scatter model using plane wave ultrasound measurements

Hans-Martin Schwab; Martin F. Beckmann; Georg Schmitz

doi:10.1364/BOE.7.001468

1. Introduction

In medical photoacoustic (PA) imaging, the observed tissue is illuminated by short light pulses, and the evolving sound waves, generated by photoacoustic energy conversion, are detected with conventional ultrasound transducers. The resulting images contain information about the optical absorption of the observed regions. By incorporating different laser wavelengths, PA imaging can provide spectral information. For that reason, major attention has been paid to PA imaging in several disciplines, such as micro-vascular imaging and cancer diagnosis [1–4].

Compared to pure optical methods, PA imaging features enhanced imaging depth. However, compared to ultrasound (US) imaging, the penetration depth of PA imaging is rather poor. Besides strongly decreasing light fluence with depth, one main limitation with high impact on penetration depth in clinical PA imaging is clutter. Indirect clutter, also known as skin-echo, is a distortion of the received PA signal that originates in strong light absorption at the tissue surface, generating high amplitude waves that propagate into the tissue and are reflected by acoustical heterogeneities. Due to the high fluence at the irradiation area on the skin surface and the strong absorption of melanin in human skin, the amplitudes of the reflected clutter waves can overcome the amplitudes of the waves carrying the desired signal information. In general, direct PA waves cannot easily be distinguished from clutter waves, which is a severe detriment on the validity of images that suffer from clutter artifacts [5,6].

Most current approaches for clutter reduction rely on decorrelation of the clutter signal by variations of the imaging geometry. Jaeger et al. follow an approach of varying the axial position of the transducer in subsequent images and tracking the deformation with ultrasound imaging [7, 8]. Decorrelation by variations of the illumination geometry in subsequent images has been proposed by Nicholas et al. [9]. The fact that reflected waves are shifted twice as much as direct PA waves during a tissue displacement induced by radiation force is exploited in [10]. These approaches, based on decorrelation, have been shown to reduce the clutter artifacts effectively, but require additional mechanical setups or well-trained clinicians. Other approaches in literature assume the clutter waves to be of lower spatial coherence than direct PA waves and have been shown to improve the contrast in PA images by using a short-lag weighted reconstruction [11,12].

The contrast in images that are affected by clutter strongly depends on the distance between light irradiation and the transducer. The best illumination of the imaging plane and, consequently, the highest signal amplitudes are achieved by an irradiation underneath the transducer, for example by employing an acoustical standoff or a setup in water. However, the desired in-plane irradiation is usually problematic, due to even stronger clutter amplitudes. These cause the irradiation position with maximum contrast to be shifted away from the transducer, if clutter is not corrected for [5,6].

The approach presented here focuses on setups with in-plane irradiation. The idea is to estimate the wave field of the reflected PA wave on the basis of reflected US waves that are generated using a linear array transducer. Therefore, the primary clutter source needs to be localized in the reconstructed image, which requires in-plane irradiation. To reproduce the clutter wave, plane wave ultrasound acquisitions are incorporated as they can easily be superposed in post-processing to imitate arbitrary wave fronts. Employing plane waves also allows for simple relations and fast algorithms in the frequency domain.

The idea to exploit ultrasound measurements to retrieve information about the scattering behavior for clutter reduction has also been assessed by Singh et al. in [13]. This approach uses focused ultrasound that is emitted onto high intensities in the PA image to identify clutter artifacts.

In the following, a relation in the frequency domain is derived that links the plane wave US measurements and a clutter-free PA measurement to the disturbed PA measurement. This relation is inverted in order to attain the clutter-free PA measurement. Followed by a standard reconstruction, a PA image with distinctly reduced clutter artifacts can be computed. The proposed method has been implemented for both simulated and experimental measurement data and the capabilities of the method are presented below.

2. Theory

The underlying theory for this clutter reduction approach is based on two dimensional wave propagation. This implies that the pressure is assumed constant along the x-dimension in a three dimensional space. The accuracy of this simplification in photoacoustic imaging has been discussed in [14].

The quantities used in the following can be divided into functions of r = (y, z) ∈ Ω in object space Ω ∈ ℝ² and functions of (y, t) ∈ δΩ × [0, T] in measurement space δΩ × [0, T] ∈ ℝ². For the latter, measurements will be restricted to pressure signals on the line z = 0, which is assumed to be the location of a linear array sensor. This implies quantities in object space to equal 0 for z ≤ 0. The measurement geometry is depicted in Fig. 1.

Fig. 1 Geometrical arrangement, object region Ω at z > 0 with acoustical heterogeneities γ_κ and γ_ρ, as well as initial pressure distribution p₀, infinite line sensor at z = 0, plane US waves from ϑ ∈ [0, π].

Download Full Size | PDF

Quantities in object space are the initial pressure distribution p₀(r), which is aimed for in the PA reconstruction, as well as the acoustic medium properties compressibility κ(r) and mass density ρ(r). These are hereafter expressed by their variations γ_κ(r) = κ(r)/κ₀ − 1 and γ_ρ(r) = 1 − ρ₀/ρ(r), respectively, with κ₀ and ρ₀ being their respective steady parts. The average speed-of-sound in Ω is $c_{0} = 1 / \sqrt{ρ_{0} κ_{0}}$ .

In the measurement space, the photoacoustic measurement including reflected waves is described by p_pa(y, t). The reflections are caused by the heterogeneous medium properties γ_κ and γ_ρ. p_pa(y, t) is treated as a superposition of the pure photoacoustic measurement p_h(y, t) that would have been received in a homogeneous medium, and the measurement of the scattered wave p_sc(y, t):

p_{pa} (y, t) = p_{h} (y, t) + p_{sc} (y, t)

The last quantity in the measurement space is the plane wave ultrasound measurement, which is described by p_us(y, t, ϑ) with ϑ being the angle of incidence relative to the y–axis. It will be shown that p_sc can be expressed as a function of p_h and p_us. Since p_pa and p_us can actually be measured, p_h is the only unknown quantity in (1) and can thus be solved for.

Even though the object functions are considered unknown in the following derivation, it will be shown that the quantities in measurement space contain all necessary information. This is due to the following direct relations between the object functions and the measurement functions in the spatial and temporal frequency domain. The arguments of the functions in frequency domain correspond to the arguments in time and space domain by y ↔ k_y, z ↔ k_z, t ↔ ω. For better readability, the temporal frequency ω is hereafter expressed by its corresponding wave number, using the dispersion relation: $k_{t} = \frac{ω}{c_{0}}$ .

It has been shown in several publications that, in a homogeneous medium, the photoacoustic source p₀ relates to the photoacoustic measurement p_h as [15–18]:

p_{h} (k_{y}, k_{t}) = \frac{k_{t}}{2 c_{0} κ_{z}} p_{0} (k_{y}, κ_{z})

where κ_z is a wave number defined as

κ_{z} : = sgn (k_{t}) \sqrt{k_{t}^{2} - k_{y}^{2}}

with sgn(k_t) being the signum-function.

Analogous to that, in a plane wave setup in reflection mode with ϑ ∈ [0, π] being the angle of incidence in relation to the y–axis, the acoustic scatter properties γ_κ and γ_ρ relate to the ultrasound measurement p_us as [19–22]:

p_{us} (k_{y}, k_{t}, ϑ) = \frac{j k_{t}^{2}}{2 κ_{z}} [γ_{κ} (k_{t} e_{φ} + k_{t} e_{ϑ}) + e_{φ} \cdot e_{ϑ} γ_{ρ} (k_{t} e_{φ} + k_{t} e_{ϑ})]

For better readability, in (3), two unit vectors were introduced: e_ϑ := (cos(ϑ), sin(ϑ)) and e_φ := (k_y/k_t, κ_z/k_t). e_ϑ points into the direction of the incoming plane wave, as depicted in Fig. 1, and e_φ points into the direction of a wave-vector with magnitude k_t and y-component k_y. Thus, γ_κ(k_te_φ + k_te_ϑ) might also be written as γ_κ(k_y + k_t cos(ϑ), κ_z + k_t sin(ϑ)). These equations are essential for the following mathematical description of the scattered photoacoustic wave measurement.

2.1. PA wave reflection in the frequency domain

The first step in this derivation is to introduce the wave that originates in a photoacoustic excitation as an incoming wave for the scattering process. It is assumed that the wave is only scattered once (Born approximation). In the spatial frequency domain, the temporally evolving photoacoustic wave for t ≥ 0 can be described by [16,17]:

p_{in} (k_{y}, k_{z}, t) = p_{0} (k_{y}, k_{z}) cos (k c_{0} t)

with

k = \sqrt{k_{y}^{2} + k_{z}^{2}}

being the magnitude of the spatial frequency vector. The wave can be decomposed into plane mono-frequent waves in the spatial domain and in the temporal frequency domain by Fourier transforms and a transform into polar coordinates. Plane waves that propagate into negative z–direction, which implies angles of ϑ ∈] − π, 0[, are considered not to contribute to the reflected clutter wave and are therefore neglected. The incoming wave is then described by:

p_{in} (r, k_{t}) = - \frac{k_{t}}{4 c_{0} π} \int_{0}^{π} p_{0} (- k_{t} e_{ϑ}) exp (- i k_{t} e_{ϑ} \cdot r) d ϑ

Here, ϑ accounts for the propagation direction of the individual plane waves in relation to the y–axis. This incoming sum of plane waves is plugged into the following general solution to the inhomogeneous acoustic wave equation [22]:

\begin{array}{l} p (r, k_{t}) = & - k_{t}^{2} (γ_{κ} (r) p_{in} (r, k_{t})) * g (r) \\ + (γ_{ρ} (r) \nabla p_{in} (r, k_{t})) * \nabla g (r) \end{array}

with g(r) being the two dimensional Green function of free space and * being the two-dimensional convolution operator in Ω. Assuming Born approximation, the scattered wave in the measurement space can be derived in the frequency domain in line with Schiffner et al. [22] as:

\begin{array}{l} p_{sc} (k_{y}, z = 0, k_{t}) = & \frac{j k_{t}^{3}}{8 c_{0} π κ_{z}} \int_{0}^{π} p_{0} (- k_{t} e_{ϑ}) γ_{κ} (k_{t} (e_{φ} + e_{ϑ})) d ϑ \\ + \frac{j k_{t}^{3}}{8 c_{0} π κ_{z}} \int_{0}^{π} p_{0} (- k_{t} e_{ϑ}) e_{φ} \cdot e_{ϑ} γ_{ρ} (k_{t} (e_{φ} + e_{ϑ})) d ϑ \end{array}

which defines the forward problem that relates the quantities in object space to the space of measurement. By substituting the integration over ϑ with an integration over k_y,ϑ := k_t cos(ϑ) with dk_y,ϑ = −κ_z,ϑdϑ and κ_z,ϑ := k_t sin(ϑ), the quantities in object space can be expressed by their representations in the measurement space according to (2) and (3). This results in a representation relying only on quantities in measurement space:

p_{sc} (k_{y}, k_{t}) = - \frac{1}{2 π} \int_{- k_{t}}^{k_{t}} p_{h} (- k_{y, ϑ}, - k_{t}) p_{us} (k_{y}, k_{t}, ϑ (k_{y, ϑ})) d k_{y, ϑ}

The measurement of the scattered wave for one temporal frequency can consequently be described as an inner product of the non-scattered measurement with plane wave ultrasound measurements resulting from emissions at certain angles at the same temporal frequency.

The angles of the plane waves in the ultrasound measurements depend on the integration variable k_y,ϑ by ϑ = cos⁻¹ (k_y,ϑ/k_t). Hypothetically, this requires the acquisition of one plane wave angle for each pair of k_y,ϑ and k_t. In a practical implementation, a limited number of angles is acquired and the missing values are approximated by interpolation.

There is a more intelligible representation of the relation in (8). As e_φ and e_ϑ are always linked as sum or as scalar product in (3), which are both commutative operations, the two unity vectors are interchangeable with no impact on the result. Thus, the integration in (7) can be performed as an integration over φ instead of ϑ, if e_ϑ is replaced by e_φ in the argument of p₀. This yields a relation in measurement space, where for each frequency vector of p_sc, only one plane wave measurement p_us is read out and its product with p_h is multiplied and integrated along the lateral frequency component:

p_{sc} (k_{y}, k_{t}) = - \frac{1}{2 π} \int_{- k_{t}}^{k_{t}} p_{h} (- k_{y, φ}, - k_{t}) p_{us} (k_{y, φ}, k_{t}, {cos}^{- 1} (\frac{k_{y}}{k_{t}})) d k_{y, φ}

The angle of the plane wave is determined by the frequency vector of p_sc as

{cos}^{- 1} (\frac{k_{y}}{k_{t}})

and missing angles can, again, be attained by interpolation. A descriptive illustration of this representation of the relation is shown in Fig. 2.

Fig. 2 K-space relation between functions in measurement space for scattered photoa-coustic measurements. To get the value of p_sc at (k_t,0, k_y,0), the measurement of p_us under the plane wave angle $ϑ = {cos}^{- 1} (\frac{k_{y, 0}}{k_{t, 0}})$ and the measurement p_h are read out on the respective lines, multiplied and summed up; White areas refer to imaginary lateral wave numbers and contain no data.

Download Full Size | PDF

The relation in (8), or (9) respectively, does not sufficiently describe a scattered photoacoustic wave in general. Neither is transmission scattering considered nor does the theory account for the fact that (4) is only valid for t ≥ 0 and returns a converging wave at times t < 0. This does not resemble a photoacoustic wave propagation and falsely creates additional scattering of the converging wave. However, both of these simplifications have low impact on the application for clutter reduction, since the skin, which is the major source for the incoming wave in (5), is located between the sensor and the scatterers.

2.2. Inversion

Equation (9) describes the forward problem of a scattered PA wave in measurement space. The aim of this paper, however, is to derive p_h as a function of known quantities. Therefore, (9) needs to be inverted.

In an actual combined ultrasound and photoacoustic acquisition, the measurable quantities p_pa and p_us are sampled in both space and time. The resulting two dimensional data sets need to be Fourier transformed in both directions, which can be done efficiently by using fast Fourier methods. In this discrete representation, (8) can be expressed as one matrix-vector multiplication for each temporal frequency k_t:

{\underline{p}}_{sc} = - {\underset{̳}{P}}_{us} \cdot {\underline{p}}_{h}^{*}

Here, p̱_sc and p̱_h are column vectors in the two dimensional measurement matrix in Fourier domain, which separates temporal frequencies into columns and spatial frequencies into rows. P̳_us is a matrix that consists of a collection of respective column vectors in the ultrasound measurement matrices for different angles.

{\underline{p}}_{h}^{*}

determines the complex conjugate of p̱_h. Since, under the assumption of real values in time- and space-domain, p_h(−k_y,ϑ, −k_t) equals p_h(k_y,ϑ, k_t)^*, the conjugate accounts for the negative arguments of p_h in (8). It should be mentioned that, due to the interchangeability of k_ϑ,y and k_y as integration variable, as presented in (8) and (9), P̳_us is actually a symmetric matrix.

Substituting (8) in (1) and expressing the result in matrix-vector-notation yields:

{\underline{p}}_{pa} = \underset{̳}{I} \cdot {\underline{p}}_{h} - α \cdot {\underset{̳}{P}}_{us} \cdot {\underline{p}}_{h}^{*}

where I̳ is the unity matrix and α is a scalar scaling factor that accounts for the different signal amplitudes of the photoacoustic measurement and the plane wave measurements.

To be more precise, (2) and (3) do not contain the temporal system responses of the two imaging systems, such as the acousto-electric system response of the transducer and the transmit pulse shapes of both the laser pulse and the plane wave. If these had been considered, α would be an inverse filter that amplifies damped frequency ranges of the PA and US measurements. Here, α is simplified to be a scalar factor that needs to be found empirically and can be used for fine tuning of the result. The optimal α is found by minimizing the signal strength at the location of a known artifact. It does not depend on the content of the image but rather accounts for the unknown transmission gain of the US measurement that scales (10) linearly. Hence, the determination of α only needs to be done once for each imaging system with constant ultrasound transmission gain.

In equation (11), both p̱_h and ${\underline{p}}_{h}^{*}$ appear. Dividing the equation into real part and imaginary part yields a directly invertible equation:

(\begin{matrix} ℜ {\underline{p}}_{pa} \\ ℑ {\underline{p}}_{pa} \end{matrix}) = (\begin{matrix} ℜ [\underset{̳}{I} - α {\underset{̳}{P}}_{us}] & - ℑ α {\underset{̳}{P}}_{us} \\ ℑ α {\underset{̳}{P}}_{us} & ℜ [\underset{̳}{I} + α {\underset{̳}{P}}_{us}] \end{matrix}) \cdot (\begin{matrix} ℜ {\underline{p}}_{h} \\ ℑ {\underline{p}}_{h} \end{matrix})

Here, ℜ is the operator that returns the real part and ℑ is the operator that returns the imaginary part. In order to account for noise, inaccurate setups, inaccurate interpolation and model deviations, a least-square approximation is implemented to retrieve an approximate solution for p̱_h. This solution requires only the measurable quantities p̱_pa and P̳_us. This least-square approximation needs to be applied once for each contributing temporal frequency in order to return a complete data set for p_h(k_y, k_t). Subsequently, any standard photoacoustic reconstruction algorithm can be applied to p_h(k_y, k_t) to yield p₀(y, z), such as an algorithm based on (2). The processing chain of the complete inversion process is displayed in Fig. 3.

Fig. 3 Block diagram of clutter compensation algorithm.

Download Full Size | PDF

3. Results

3.1. Simulation

The clutter reduction method was tested on simulated data gained by wave field computations using the k-wave toolbox for Matlab [23]. While the clutter reduction algorithm is derived only for linear scattering, the simulation yields non-linear solutions including multiple scattering. A speed-of-sound (SOS) phantom was developed in combination with a photoacoustic source phantom containing a strongly radiating source at the top and weakly radiating sources inside the tissue, as depicted in Fig. 4. In general, the SOS distribution relates to the aforementioned acoustic medium properties by $c (r) = 1 / \sqrt{ρ (r) κ (r)}$ . A transducer with 128 elements, an element pitch of 0.245 mm and a center frequency at 10% of the sample frequency with 60% relative bandwidth was simulated for detection. The same transducer geometry was used for the generation of 100 plane waves at equidistant angles between 45° and 135°. For reference, an additional simulation with a homogeneous SOS distribution was acquired that does not contain any clutter due to the lack of reflection.

Fig. 4 Phantom for simulation, weak sources are evaluated in green ROI.

Download Full Size | PDF

During post-processing, missing data was interpolated linearly between the data sets of neighboring angles. Since shallow angles cannot accurately be generated with a limited size linear array, the values of measurements at angles smaller than 45° and larger than 135°, which had not been acquired, were set to zero. In addition, the ultrasound data was weighted with a Gaussian function over the angle to suppress the contribution of US data from shallow angles without cropping the data too strongly.

Figure 5 shows the results of the clutter reduction approach by comparing the reconstruction of the recovered measurement data to the reconstruction of unprocessed measurement data and to the reconstruction of measurement data gained from a simulation with a homogeneous speed-of-sound distribution, which represents the best achievable image. The results are displayed in a region of interest (ROI) as defined in Fig. 4 that crops the area with strong skin absorption. While, in the uncompensated image, clutter cannot be distinguished from PA sources, the amplitudes of the sources overcome the clutter signal distinctly in the image with applied clutter reduction. In comparison to the reconstruction results from the homogeneous simulation, the background signal is still higher even after clutter reduction, but the desired structures can easily be identified. The reconstruction artifacts that can be seen in the reference image are completely obscured by clutter in the other two images.

Fig. 5 Comparison of PA reconstructions of simulated data, axes in mm.

Download Full Size | PDF

A PSNR was determined by the ratio of the squared maximal intensity of the strongest absorber and a mean power in a region with no absorbers. The image with homogeneous medium exhibits a PSNR of 31.4 dB, which is only limited by reconstruction artifacts. Due to clutter, the image with heterogeneous medium properties has a PSNR of only 15.4 dB. In the image with applied clutter reduction, the PSNR could be recovered up to 20.6 dB, which implies an artifact reduction of 5.3 dB.

Even though the quality in the compensated image has been improved, the clutter artifacts could not be suppressed completely. We see two main reasons for this: Firstly, as the wave field simulation accounts for multiple scattering, the fact that the derived theory relies on a linear scatter model can limit the success. Secondly, the interpolation between the measurement data of neighboring plane wave angles, which is necessary to reduce the amount of data sets to acquire, is inaccurate. Such an interpolation resembles two weighted plane waves rather than one plane wave from a different angle. Using a higher number of plane wave angles or deriving an improved interpolation strategy might help to improve the suppression.

The entire processing time of the reconstruction compensated for clutter ranges between two to three seconds in a non-parallelized Matlab implementation. Since each temporal frequency can be treated independently, the algorithm can be sped up by parallelized implementation.

3.2. Experiment

In addition to an assessment based on simulated measurement data, an experiment was conducted that employs a simple phantom to exhibit the general functionality of the clutter reduction method. The experimental setup contained a blackened plastic foil to simulate the skin surface and a blackened metal rod of 1 mm diameter underneath the foil that was orientated perpendicular to the image plane, serving as both an absorber and a scatterer. The phantom was placed in water. For the data acquisition, a hand held diode laser PA/US transducer prototype was used, which is being developed in the EU research project Fullphase. The prototype is similar to the the one described in [24] and has a 808 nm diode stack integrated on the ultrasound transducer that can supply light pulses of up to 1 mJ.

For the plane wave US acquisition, one measurement per integer angle was acquired in a range from 45° to 135°, yielding a data set of 91 measurements. The acquired PA measurements were averaged over 50 frames. Both the US and PA measurements were band-pass filtered to suppress noise and to match up the respective frequency ranges of both modalities. Measurements relating to shallow plane wave angles were suppressed by applying a Tukey window onto the rows of P̳_us in (11).

Figure 6 a displays the US reconstruction incorporating all acquired angles. The foil appears as a bright bar structure at 6.2 mm depth. Underneath, at a depth of 9.8 mm, the metal rod can be located as a bright point. Figure b shows the PA reconstruction without compensation. The foil and the metal rod can both be seen in the PA image. Due to reflection of the wave originating in the foil, a clutter artifact appears in the PA image at 13.4 mm depth. As expected, the distance between artifact and foil is exactly twice the distance between metal rod and foil. In Fig. c, the PA reconstruction including clutter reduction is displayed. The clutter exhibits a distinctly reduced amplitude, even though it does not vanish completely. In contrast to that, the amplitude of the absorbing rod remains on an equal level and so does the background noise. In order to provide a more detailed impression of the clutter reduction, a region of interest is magnified and displayed in Fig. 7. Here, the background noise is more visible. In the compensated image, the background signal is maintained at an equal level, while the clutter amplitude is decreased. The relation between maximum intensity of the absorber to maximum intensity of the artifact was reduced by 5.0 dB.

Fig. 6 Comparison of Reconstructions of experimental data.

Download Full Size | PDF

Fig. 7 Comparison of clutter artifacts in uncompensated and compensated image.

Download Full Size | PDF

4. Discussion and conclusions

We have predicted a reduction of clutter artifacts in medical photoacoustic images by incorporating plane wave ultrasound measurements and processing the acquired data in combination with the distorted PA data to reproduce a clutter-reduced PA image. The methods can be applied to all devices or setups with in-plane light irradiation.

The general applicability of the proposed method was shown in a simple experimental setup that provides a high degree of ground truth. The contrast between clutter artifacts and PA sources has clearly been improved.

In addition, measurement data gained by a non-linear wave field simulation were acquired to further assess the clutter reduction method. In this assessment, more clinical relevant tissue structures could be investigated than in the experimental setup, while ground truth could still be assured. The results of this simulation based study promise a practicality of the proposed method for medical applications. PA sources can be clearly distinguished from clutter artifacts.

The influence of out-of-plane clutter on the performance of the method needs to be investigated in further studies, as neither the experimental study nor the simulation study involved scatterers that are located outside the imaging plane. Assuming a strong elevational focus of the ultrasound transducer, the influence of out of plane clutter, however, is considered low. In vivo studies will be performed in the future to validate the algorithm in a more clinically relevant context.

The scaling factor α, which accounts for different signal amplitudes in PA and US measurement, might be replaced by an inverse filter that considers the transfer functions of the PA and US imaging systems. Currently, α is found empirically. Even if the system transfer functions are not considered, the perfect choice for α might be determined by a system calibration or an adaptive algorithm in the future.

Currently, the number of plane wave angles that need to be emitted to ensure stable results is very high, which is a limiting factor in terms of real-time processing. With a smaller number of angles, both the time for data acquisition and the time for data processing might decrease significantly. There is a high potential for reducing the number of required angles by employing a more advanced interpolation method. A better interpolation might also affect the necessary amount of zero padding, which has a strong influence on computation time.

The fact that the inversion process can be calculated independently for each temporal frequency encourages the real-time capability of the described method, since it implies a high degree of parallelizability. Assuming real-time processing, this automated clutter reduction method that has no demand on additional hardware or training, might be a beneficial contribution to the introduction of photoacoustics as a bedside diagnosis tool in medical imaging.

Acknowledgments

The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 318067.

References and links

1. L. V. Wang and H.-i. Wu, Biomedical Optics (John Wiley & Sons, 2007).

2. L. V. Wang, Photoacoustic imaging and spectroscopy (Taylor & Francis Group, LLC, 2009). [CrossRef]

3. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77, 041101 (2006). [CrossRef]

4. P. Beard, “Biomedical photoacoustic imaging: a review,” Interface Focus pp. 602–631 (2011). [CrossRef]

5. M. Frenz and M. Jaeger, “Optimization of tissue irradiation in optoacoustic imaging using a linear transducer: theory and experiments,” Proc. SPIE 6856, 68561Y (2008). [CrossRef]

6. G. Held, S. Preisser, S. Peeters, M. Jaeger, and M. Frenz, “Effect of irradiation distance on image contrast in epi-optoacoustic imaging of human volunteers,” Biomed. Opt. Express 5, 3765–3780 (2014). [CrossRef] [PubMed]

7. M. Jaeger, L. Siegenthaler, M. Kitz, and M. Frenz, “Reduction of background in optoacoutic image sequences obtained under tissue deformation,” J. Biomed. Opt. 14054011 (2009). [CrossRef]

8. M. Jaeger, D. Harris-Birtill, A. Gertsch, E. O’Flynn, and J. Bamer, “Deformation-compensated averaging for clutter reduction in epiphotoacoutic imaging in vivo,” J. Biomed. Opt. 17066007 (2012). [CrossRef]

9. T. N. Erpelding, H. Ke, and L. Wang, “In-place clutter reduction for photoacoustic imaging,” (2014).

10. M. Jaeger, J. C. Bamber, and M. Frenz, “Clutter elimination for deep clinical optoacoustic imaging using localised vibration tagging (LOVIT),” Photoacoustics 1, 19–29 (2013). [CrossRef] [PubMed]

11. E. J. Alles, M. Jaeger, and J. C. Bamber, “Photoacoustic clutter reduction using short-lag spatial coherence weighted imaging,” in Proceedings of IEEE International Ultrasonics Sympossium (IEEE, 2014), pp. 41–44.

12. B. Pourebrahimi, S. Yoon, D. Dopsa, and M. C. Kolios, “Improving the quality of photoacoustic images using the short-lag spatial coherence imaging technique,” Proc. SPIE 8581, 85813Y (2013). [CrossRef]

13. M. Kuniyil Ajith Singh and W. Steenbergen, “Photoacoustic-guided focused ultrasound (PAFUSion) for identifying reflection artifacts in photoacoustic imaging,” Photoacoustics 3, 123–131 (2015). [CrossRef]

14. K. P. Köstli and P. C. Beard, “Two-dimensional photoacoustic imaging by use of fourier-transform image reconstruction and a detector with an anisotropic response,” Appl. Opt. 42, 1899–1908 (2003). [CrossRef] [PubMed]

15. N. Baddour, “A multi-dimensional transfer function approach to photo-acoustic signal analysis,” J. Franklin Inst. 345, 792–818 (2008). [CrossRef]

16. K. P. Köstli, M. Frenz, H. Bebie, and H. P. Weber, “Temporal backward projection of optoacoustic pressure transients using Fourier transform methods,” Phys. Med. Biol. 46, 1863–1872 (2001). [CrossRef] [PubMed]

17. B. Cox and P. Beard, “Fast calculation of pulsed photoacoustic fields in fluids using k-space methods,” J. Acoust. Soc. Am. 117, 3616 (2005). [CrossRef] [PubMed]

18. Y. Xu, D. Feng, and L. Wang, “Exact frequency-domain reconstruction for thermoacoustic tomography. I. Planar geometry,” IEEE Trans. Med. Imaging 21, 823–828 (2002). [CrossRef] [PubMed]

19. A. Devaney, “A filtered backpropagation algorithm for diffraction tomography,” Ultrason. Imaging 4, 336–350 (1982). [CrossRef] [PubMed]

20. A. Devaney, Mathematical Foundations of Imaging, Tomography and Wavefield Inversion (Cambridge University Press2012). [CrossRef]

21. J.-y. Lu, “Experimental study of high frame rate imaging with limited diffraction beams,” IEEE Trans. Ultrason. Ferr. 45, 84–97 (1998). [CrossRef]

22. M. F. Schiffner and G. Schmitz, eds., “Plane wave pulse echo ultrasound diffraction tomography with a fixed linear transducer array,” in Acoustic Imaging (Springer, 2012). [CrossRef]

23. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave-fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef]

24. K. Daoudi, P. van den Berg, O. Rabot, A. Kohl, S. Tisserand, P. Brands, and W. Steenbergen, “Handheld probe combining laser diode and ultrasound transducer array for ultrasound/photacoustic dual modality imaging,” Opt. Express 22, 26365–26374 (2014). [CrossRef] [PubMed]

Photoacoustic clutter reduction by inversion of a linear scatter model using plane wave ultrasound measurements

Abstract

Corrections

1. Introduction

2. Theory

2.1. PA wave reflection in the frequency domain

2.2. Inversion

3. Results

3.1. Simulation

3.2. Experiment

4. Discussion and conclusions

Acknowledgments

References and links

Cited By

Figures (7)

Equations (12)

Biomedical Optics Express