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

: Photoacoustic imaging aims to visualize light absorption properties of biological tissue by receiving a sound wave that is generated inside the observed object as a result of the photoacoustic effect. In clinical applications, the strong light absorption in human skin is a major problem. When high amplitude photoacoustic waves that originate from skin absorption propagate into the tissue, they are reﬂected back by acoustical scatterers and the reﬂections contribute to the received signal. The artifacts associated with these reﬂected waves are referred to as clutter or skin echo and limit the applicability of photoacoustic imaging for medical applications severely. This study seeks to exploit the acoustic tissue information gained by plane wave ultrasound measurements with a linear array in order to correct for reﬂections in the photoacoustic image. By deriving a theory for clutter waves in k-space and a matching inversion approach, photoacoustic measurements compensated for clutter are shown to be recovered.


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][2][3][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 inplane 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.

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 Ω ∈ R 2 and functions of (y,t) ∈ δ Ω × [0, T ] in measurement space δ Ω × [0, T ] ∈ R 2 . 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.
Quantities in object space are the initial pressure distribution p 0 (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)/κ 0 − 1 and γ ρ (r) = 1−ρ 0 /ρ(r), respectively, with κ 0 and ρ 0 being their respective steady parts. The average speedof-sound in Ω is c 0 = 1/ √ ρ 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. Figure 1. Geometrical arrangement, object region Ω at z > 0 with acoustical heterogeneities γ κ and γ ρ , as well as initial pressure distribution p 0 , infinite line sensor at z = 0, plane US waves from ϑ ∈ [0, π].
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 = ω c 0 . It has been shown in several publications that, in a homogeneous medium, the photoacoustic source p 0 relates to the photoacoustic measurement p h as [15][16][17][18]: where κ z is a wave number defined as κ z := sgn(k t ) k 2 t − k 2 y with sgn(k t ) being the signumfunction.
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][20][21][22]: 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 t e ϕ + k t e ϑ ) 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.

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]: with k = k 2 y + k 2 z 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: 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]: with g(r) being the two dimensional Green function of free space and * being the twodimensional 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: 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: 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 −1 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 0 . 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: The angle of the plane wave is determined by the frequency vector of p sc as cos −1 ( 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.
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.

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 : 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. 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: for each kt 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 p * h appear. Dividing the equation into real part and imaginary part yields a directly invertible equation: 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 0 (y, z), such as an algorithm based on (2). The processing chain of the complete inversion process is displayed in Fig. 3.

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/ ρ(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. 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 speedof-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.
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.

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

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.