Wave front analysis for enhanced time-domain beamforming of point-like targets in optoacoustic imaging using a linear array.

Using linear array transducers in combination with state-of-the-art multichannel electronics allows to perform optoacoustic imaging with frame rates only limited by the laser pulse repetition frequency and the acoustic time of flight. However, characteristic image artefacts resulting from the limited view and a lower SNR when compared to systems based on single-element focused transducers represent a burden for the clinical acceptance of the technology. In this paper, we present a new method for the improvement of image quality based on the analysis of the signal amplitudes along summation curves during the delay-and-sum beamforming process (DAS). The algorithm compares amplitude distributions along wave fronts with theoretical patterns from optoacoustic point sources. The method was validated on simulated and experimental phantom as well as in-vivo data. An improvement of the lateral resolution by more than a factor of two when comparing conventional DAS and our approach could be shown (numeric and experimental phantom data). For instance, on experimental data from a wire phantom, a PSF in the range of 0.18-0.22 mm was obtained with our approach against 0.48 mm for standard DAS. Furthermore, the SNR of a subcutaneous vessel 2.5 mm below the skin surface was improved by about 30 dB when compared to standard DAS.

Introduction

Optoacoustic techniques combine the advantages of acoustics and optics by making deep tissue regions beyond the propagation range of ballistic photons accessible to imaging with optical contrast [[1], [2], [3]]. Signal generation is performed with pulsed laser beams that can be applied either in a focused way or through wide field illumination. In the first case, the optical focusing behavior defines the image quality [[4], [5], [6]], while the reception characteristics of the used acoustic detection system define the resolution in the second case [7,8].

When it comes to applications where realtime imaging is needed, optoacoustic systems involving ultrasound arrays are the method of choice [[9], [10], [11], [12], [13], [14]]. Furthermore, wide-field illumination such as used in linear-array optoacoustics is advantageous with respect to the MPE (maximum permissible exposure). Medical ultrasound transducers are often used in such settings due to their availability and their suitability for real-time freehand imaging. In contrast to systems based on single element focusing transducers, arrays allow, furthermore, for an almost homogeneous resolution throughout the whole field of view when images are reconstructed using suitable beamforming algorithms, either in time [[15], [16], [17]] or frequency domain [18,19].

For improvement of the reconstruction quality, different methods have been investigated such as interpolation of complex coefficients in Fourier transform based algorithms [20]. In the time domain, improvement of image quality can be achieved by iterative approaches [21,22], which in particular allow to avoid artefacts that result of simplifications such as homogeneous acoustic properties of the bulk medium. Quantitative optoacoustic imaging, which allows to retrieve parameters such as the optical absorption coefficient [23,24], often relies on iterative approaches [25]. Model based approaches have the advantage of allowing to incorporate knowledge regarding the detection system, such as the influence on the relative positions of the optoacoustic source and the transducer element used for detection on the generated signal´s shape (e.g. its amplitude and frequency spectrum) [[26], [27], [28]]. Other methods for improvement of the resolution and the SNR in optoacoustic reconstruction rely on the modification of standard delay-and-sum beamforming by integration of different coefficients during summation. For instance, coherence weighting [29], which has shown to significantly improve the image quality in confocal optoacoustic imaging using single element systems, can also be applied in linear array optoacoustics [30]. Furthermore, the combination of delay-and-sum and multiply methods [31] can be used for improved reconstruction. Moreover, adaptive beamforming known from ultrasound imaging, which introduces some weighted apodization coefficients, allows enhancement of both resolution and SNR [[32], [33], [34], [35]].

In the present work, a new method for time delay beamforming is introduced, which analyses the signal amplitudes during the reconstruction of individual image points. The improvement of SNR and resolution will be demonstrated on simulation data, phantom data and in-vivo experiments.

Theory

Distribution of optoacoustic signal amplitudes on transducer array elements

In time domain reconstruction, the beamforming process is based on the phase-adjusted summation of signal samples from potential optoacoustic sources. For each pixel of a reconstructed image located at (x,z), the time of flight between the pixel and the individual detectors at position (x is calculated and corresponding signal samples are summed up according towhere b is the pixel value at position (x,z), is the acquired pre-beamformed channel data from N elements and the time is defined aswhere c is the speed of sound. This is performed for each pixel of the image, without the a priori knowledge if this pixel actually corresponds to a source of an optoacoustic wave. The amplitude of the signal depends on the initial pressure distribution in the tissue, which itself results from the optical properties (absorption, scattering and anisotropy coefficients) and the fluency of the applied laser pulse. Furthermore, the propagation through the tissue and the properties of the transducer have an influence on the acquired signal amplitude and temporal shape.

If we consider the case of a homogeneously illuminated spherical optoacoustic source with diameter a leading to an N-shaped signal [36,37], the pre-beamformed channel data for the n-th element of an array can be derived aswhere , with being the position of the source and the position of the m-th sub-element of the n-th element of the transducer array. U is the step function, which equals 1 when the argument is positive and 0 when the argument is negative. is the impulse response of a single array element. For a realistic simulation, the surface of an array element should be discretized in sub-elements with an extent . If the duration of the laser pulse used for generation of the optoacoustic signal is sufficiently short (much shorter than the oscillation period of the used transducer), it can be neglected. Otherwise, Eq. (3) must be convolved with the temporal profile of the laser pulse for a realistic description of the signal of a homogeneously illuminated spherical source.

In the search for a new nonlinear filter coefficient, we consider the maximum amplitudes A of the signals on the different array elements. If we assume that distances R from a source to the different array elements are comparable, the maximum signal amplitudes on all array elements should not vary from element to element and . On the other side, if the source is rather close to the array so that the distance R varies from element to element, while the individual array element is so small that the summation in (3) can be neglected, the amplitudes are modulated with a factor of 1/R so that

However, this only holds true when the transducer elements are small when compared to the wavelength, which mostly is not the case, especially for most medical ultrasound transducers. For a more realistic assessment of the distribution of signal amplitudes from an optoacoustic source along the elements of a transducer array, the angular sensitivity of a transducer element needs to be considered. In fact, in a 2D model, the radial sound field of a rectangular transducer of aperture size L, whose source distribution rect(x) can be described ascan be approximated by a sinc-function. Accordingly, due to the analogy between transmit and receive sound-field, a more general description of the amplitude distribution can be given aswhere

Obviously, the sinc distribution of the angular sound field does only strictly hold true for a monochromatic signal.

Nonlinear filter coefficients based on signal amplitude distribution

In order to better differentiate between pixels which actually correspond to signals sources and pixels which do not, we introduce some new coefficients based on the analysis of the signal amplitude distributions presented in the last section. When a pixel is actually the source of an optoacoustic signal, the amplitude distribution of the signals on the different transducer elements should correspond to . Accordingly, the deviation of the amplitude pattern from the theoretical shape can be used as confidence measure whether the considered point is an optoacoustic source or not. In the simplest case, if we ignore differences in the times of flight and the sound field properties of the individual element, the deviation from the (constant) theoretical pattern corresponds to the standard deviation. However, given that a low deviation must result in a high confidence factor and a large deviation in a low one, we introduce the inverse of the relative standard deviation as confidence parameter withwhere is defined such as described above. Here, is the average of the sample amplitudes summed up during reconstruction of the pixel located at position (x,z). The relative standard deviation was preferred over the absolute standard deviation to differentiate between the cases of an actual optoacoustic source, where high amplitude signals are summed up during beamforming, and the case of a point not corresponding to a source, where only noise values are summed up during beamforming. While the absolute standard deviations can be comparable in these cases, the relative ones take into account the average signal amplitudes, which are much higher for an actual source than for non-source points.

Given that this confidence factor is based on the assumption that sample values are constant, the approach is particularly suitable for distant positions, where differences in amplitude due to varying propagation distances can be neglected. For positions closer to the aperture, the dependency has a greater influence. This can be explained by considering the amplitude distribution for the case of an aperture of length L and two optoacoustic sources at (x,z) = (0,L) and (x,z) = (0,L/4). Under the amplitude distribution assumption, the ratio of the signal amplitudes from the element in the middle of the aperture and from the last element (positioned at (x,z) = (0,L/2)) is for the first source, while it is for the second point being closer to the aperture. We therefore introduced a second coefficient where the amplitude distribution is compared to a best fit function reflecting the mentioned the dependency:with being a best fit to the samples corresponding to the point of interest at according towhere

This confidence measure takes into account the influence of the different acoustic path lengths, however, the sound field properties of the individual transducer element, which become more and more relevant with increasing element size, are still not considered. Accordingly, a third confidence parameter is introduced, in which the signal amplitudes are compared with a best fit function reflecting the sinc pattern corresponding to the sound field of a rectangular transducer element:with being a best fit to the samples corresponding to the point of interest at according towhereand

The confidence parameters , and are calculated for each pixel located at position (x,z) so that a “confidence image” is available in addition to the beamformed image b(x,z). For improvement of the image quality, the beamformed image and the “confidence image” are multiplied pixelwise so that the latter acts as an image filter. For this purpose, Eq. (1) needs to be modified so that the filtered image is reconstructed according towhere i can be 1,2 or 3 depending on the used filter type (standard deviation, 1/r- or sinc-fit, respectively). The influence of the different filter parameters on the image resolution and the contrast is shown in the following chapters for simulated and experimental data.

Simulation

Simulation model

The concept of the new filter coefficients is illustrated in the following with an example from simulated data. Fig. 2 represents a wave front generated after absorption of light by a point source located at (x/z)=(0 mm/ 36.5 mm) in front of a transducer array having 128 elements (pitch = 0.3 mm, element size 0.25 mm × 7 mm in lateral and elevational direction respectively) at a frequency of 7.5 MHz (−6 dB Gaussian bandwidth of 5 MHz). The simulation was performed in Matlab. An illustration of the geometry is given in Fig. 1. For simulation of the elementary signal of a point source, an N-shaped signal with a temporal extent much smaller than one period of the 7.5 MHz oscillation was assumed according to (3). The surface of the transducer element was segmented with a stepsize and the single waves travelling to the individual transducer element segments were summed up in time domain, with a weighting of accounting for the propagation as spherical wave. The final signal was then convolved with the impulse response of the transducer , defined as a Gaussian as described above. The duration of the laser pulse was supposed to be much shorter than the oscillation period of the transducer, so that it could be neglected in the simulation.

Fig. 2

Simulated wave front (blue) and summation curve (yellow) for different pixels (left) and samples corresponding to the summation curves (right). (a) shows the case of a pixel that actually corresponds to the position of an optoacoustic source. (d) shows the amplitudes of the signal values along the yellow curve in (a). (b) shows the case of a pixel that has an axial shift with respect to the optoacoustic source. Accordingly the summation curve (yellow) does not overlap with the wave front (blue). (e) shows the signal values corresponding to the yellow curve in (b). (c) shows the case of a pixel with a lateral offset to the optoacoustic source. The summation curve (yellow) and the actual wave front slightly overlap, which results in the summation of a few high amplitude samples and many low amplitude (noise) samples, as can be seen in (f).

Fig. 1

Geometry of simulated and experimental set-up with linear transducer (128 element, pitch 0,3 mm).

Different levels of random noise could be added to the simulated channel data. For this purpose, a random matrix of the same size than the channel data matrix was generated with the inbuilt Matlab rand-function, scaled to a defined percentage of the maximum of the pre-beamformed channel data matrix and added to the latter. With the simulated data obtained as presented above, the principle of the algorithm is illustrated for the case of the confidence parameter in Fig. 2, where three summation curves corresponding to the actual position of the point source (a), to a position having an axial offset to the point source (b), and to a third point having a lateral offset to the point source (c). As can be seen in this figure, the yellow summation curve corresponds perfectly to the actual wave front of the absorber, when the pixel corresponding to the position of the absorber is reconstructed (Fig. 2(a)). Accordingly, only samples corresponding to actual signals are summed up in this case. Under the assumption that each array element has a comparable sensitivity, and that the amplitude differences due to slightly varying propagation distances between the origin of the wave and the different receivers can be neglected, the summed samples should have comparable values. The distribution of sample values along the theoretical summation curve (shown in the right side in Fig. 2) is therefore characteristic for the position of the reconstructed pixel. If the pixel corresponds to an actual source, a relatively homogeneous distribution can be expected (Fig. 2(d)). If the reconstructed pixel is sufficiently far away from the actual source and the summation curve does not overlap with an actual wave front, only noise signals are summed up during beamforming (Fig. 2(e)). Finally, if the pixel is located near an actual source, the distribution of samples consists of many noise samples and few high amplitude samples, which graphically correspond to the intersection between the summation curve (yellow) and the actual wave front in Fig. 2(c). The corresponding distribution is shown in Fig. 2(f).

In order to enhance pixels corresponding to actual sources and at the same time apply a damping coefficients to other pixels, sample distributions such as in Fig. 2(d) have to be amplified (or at least remain constant) while such as in Fig. 2(e) or (f) need to be suppressed.

To illustrate the ability of the confidence parameter to highlight points corresponding to actual sources of optoacoustic signals, we compared the values for the three cases such as described in Fig. 2.

As shown in Table 1, the absolute standard deviations are comparable in the cases a and c, respectively a pixel corresponding to an actual source and a pixel with a lateral offset to an actual source whose summation curve partially overlaps with an actual wave front. Since the coefficient should damp both the cases b and c while amplifying case a, using the inverse of the relative standard deviation is the method of choice. As can be seen in Table 1, this coefficient amplifies a pixel corresponding to an actual signal source by a factor of 4, while it applies a strong damping to the pixels corresponding to the cases b and c.

Table 1

Wave front homogenity coefficient for different pixels.

	Case a	Case b	Case c
Mean value	712	2.7	2.5
Standard deviation	171	21	110
Relative std dev	0.24	7.79	43.99
1/(rel std dev)	4.2	0.13	0.02

Constraints for sinc-based confidence parameter

The sinc-dependency of the angular sensitivity of a rectangular transducer strictly holds true only for a monochromatic wave. In order to study the influence of a broader signal spectrum, we calculated the angular sensitivity distributions as a function of the bandwidth. For this purpose, an optoacoustic point source leading to an N-shaped signal was assumed at (x/z) = (0 mm/36.5 mm) and the amplitudes of acquired signals on the different elements of an array transducer (7.5 MHz centre frequency, pitch of 0.3 mm, 128 elements) were compared. As can be seen in Fig. 3, the angular sensitivity gets close to a sinc-function, when the bandwidth gets more narrow, which is well in agreement with theory. However, even with bandwidths close to 100%, which is the range of diagnostic probes, the sinc-function is a better fit than a 1/r dependency (see Fig. 2).

Fig. 3

Angular sensitivity distributions of a 128 element transducer with 7.5 MHz centre frequency, a pitch of 0.3 mm transducer, different -6 dB bandwidths for a source positioned at (x,z) = (0/36.5 mm).

Confidence-parameter enhanced DAS

For the assessment of the suitability of the different filter coefficients to enhance the quality of optoacoustic delay-and-sum beamforming, channel data were simulated in Matlab. Again, a point absorber was assumed at (x/z) = (0 mm/36.5 mm) in a non-absorbing background with a sound velocity of 1485 m/s. A linear array transducer with 7.5 MHz centre frequency, a −6 dB bandwidth of 5 MHz, a pitch of 0.3 mm and 128 elements (each 0.25 mm × 7 mm) was assumed. Random noise (maximum noise corresponds to 2% of signal maximum) was added to the simulated channel data prior to beamforming, which resulted in a SNR of 40 dB. The SNR was defined as the ratio between the signal maximum and the average of the signal´s absolute value in a location where there is no target. The sampling rate was set to 80 MSamples/s.

Fig. 4 confirms that the sinc-filter outperforms the other filter functions when it comes to the resolution of the reconstruction point source. In the absence of noise, a (theoretical) FWHM of the point-spread-function of 65 μm is achieved for the sinc-filter, while the STD and 1/r filter lead to a resolution of 120 and 145 μm respectively. In contrast, in this simulation model, the resolution of the standard delay-and-sum beamforming without further filtering is 290 μm.

Fig. 4

1D and 2D PSF after beamforming of the simulated channel data of a single point absorber at (x/z)=(0 mm/36.5 mm) using the different filter coefficients. Dynamic Range (DR) is given for each of the 2D reconstructions.

Fig. 5(b) shows the influence of noise on the performance of the different filters. While the standard beamforming is poorly affected by increasing noise levels, the resolution of the filtered beamforming decreases when more noise is added to the data. The results shown in Fig. 4 are not included in Fig. 5(b), since Fig. 4 was simulated in the absence of noise, which corresponds to an infinite SNR. Fig. 5(b) shows that the performance difference between the different filters vanishes for a SNR lower than 20 dB. This can easily be explained by Fig. 5(a), where the fit of the different functions to the channel data used for reconstruction of a defined pixel is shown. Obviously, with higher noise levels, the quality of the fit decreases. For the same reason, the achievable FWHM is lower in a medium with many signal sources close to each other since data points corresponding to wave fronts of other sources will lead to outliers in the summation data and thereby deteriorate the quality of the fit. Fig. 5(c) and (d) furthermore show the influence of the axial position of the target on the filter´s effect in terms of resolution and SNR. In these simulations, the SNR in the pre-beamformed input data was 26 dB. The results show that all filters perform better when targets are farther away from the aperture than for structures close to the array. This is particularly pronounced in the first 20 mm while the resolution remains approximately constant between 20 and 40 mm.

Fig. 5

Example of sinc and 1/r fit to the channel data summed up for a pixel at (x/z) = (0 mm/36.5 mm) (a), influence of different noise levels on the performance of the filter algorithms (b), influence of the axial position of the target on the resolution and contrast improvement ((c) and (d) respectively).

Experimental results

Materials and methods

In a next step, the performance of the different filter algorithms was evaluated on the basis of experimental data. The chosen setting was analogous to the one from the simulation. A 7.5 MHz linear array transducer (LA-7, VERMON SA) having 128 elements at a pitch of 0.3 mm was used. Optoacoustic signals were acquired with IBMT´s in-house developed multichannel electronics platform DiPhAS (Digital Phased Array System) having 128 independent channels with a digitization rate of up to 80 MSamples/s and a frame rate limited by the pulse repetition frequency of the laser. The system has been tested with respect to compliance with IEC 6060-1-2:2007 (electromagnetic compatibility) and IEC 60601-1:2006 (electric safety). Furthermore, the maximum optical output of the laser at the distal end of the fibre has been tested according to IEC 60825-1:2007 and was measured to 287.5 J/m², which is far below the MPE for skin at the used wavelength of 1064 nm. Accordingly, the compliance with all safety-relevant standards was tested, so that the system could be evaluated in-vivo. When it comes to data acquisition, raw data (pre-beamformed channel data) were transferred from the system to an embedded personal computer and reconstructed in real-time using standard delay-and-sum beamforming in the control and data processing software of the DiPhAS. Real-time reconstruction capabilities were achieved by massive parallel computing on a GPU. After the reconstruction, the data sets were stored in Fraunhofer IBMT´s ORB data format for deeper offline analysis. For the generation of the optoacoustic signals, an Nd:YAG laser (nano S, Litron Lasers, UK) operating at its fundamental wavelength of 1064 nm and providing pulses of 10 ns was used. For synchronization, DiPhAS generated a TTL signal to drive the laser flashlamp while the Q-switch was triggered internally by the laser controller. The phantom consisted of black wires immersed in water with 4% milk for ensuring optical scattering.

Experimental phantom data

Optoacoustic channel data were reconstructed with the different filter options in addition to a standard delay-and-sum beamforming. The channel data from five black wires (diameter of 100 μm) is shown in Fig. 6(a) and a reconstruction with a standard beamforming is given in Fig. 6(c). For comparison of the different algorithms, the 1D and 2D point-spread-functions of a defined point source at (x/z) = (0 mm/ 36.5 mm) are shown as close up in Fig. 6(d).

Fig. 6

(a) optoacoustic channel data after irradiation of 5 black wires in a water phantom. (b) Channel data summed up for (x/z) = (0 mm/ 36.5 mm) and fits corresponding to the filter functions. (c) reconstructed image of the whole field of view without any filter. (d) 1D PSF and close up 2D PSF of the optoacoustic source located at (x/z) = (0 mm/ 36.5 mm) after reconstruction with the different filters. For calculation of the dynamic range DR, the mean of the absolute noise signal in an image area with no optoacoustic source was compared to the maximum.

For the four algorithm types, the FWHM values (mean value +/− standard deviation) are 482 μm, (+/− 50 μm), 236 μm (+/− 44 μm), 222 μm (+/− 54 μm) and 188 μm (+/− 30 μm) for the standard delay-and-sum beamforming, the beamforming with standard deviation filter, with 1/r filter and with sinc filter respectively. The FWHM for the 1/r and sinc filter appear to be identical in Fig. 6(d), which is a comparison of the effect of the different filters applied to the area highlighted by a red square in Fig. 6(c). However, after averaging the FWHM values of the five reconstructed wires in the phantom, the sinc filter performs slightly better than the 1/r one and outperforms all other filter parameters when it comes to the achievable resolution of a reconstructed point source. In addition, given that the diameter of the wire in the phantom was 100 μm, all filters lead to an improvement of the image quality in the sense that the outcome is more similar to the true measurement. When it comes to SNR, it has to be mentioned that the noise levels of the wave front signals are between 10 and 30% of the maximum amplitudes (SNR of 17-26 dB), which corresponds to the noise range in which the sinc filter performs better than the others according to the analysis presented in Fig. 5(b).

In-vivo proband data

In a next step, in-vivo optoacoustic data of a human finger from a healthy volunteer were acquired in order to assess the performance of the filters in noisy datasets with multiple optoacoustic sources instead of single absorbers. The compliance of the system with safety-standards for medical devices was tested by certified laboratories, so that in-vivo data could be acquired from volunteers. For this purpose, the hand of the volunteer was immersed in a water bath and positioned on a hand support for avoiding motion artefacts. Volume data sets were acquired by translation of the probe in elevational direction with a stepsize of 400 μm. Cross-sectional images were reconstructed off-line with an isotropic pixelsize of 25 μm in x and z direction. No apodization was applied during beamforming. Fig. 7(a) shows the channel data that is summed up for reconstruction of a pixel belonging to a subcutaneous vessel (marked with an arrow in Fig. 7(c)–(f)). For comparison, the 1/r and sinc best-fit curves are shown in the graph as well. Despite the noisiness of the data, a strong contrast and resolution enhancement can be achieved such as can be seen in Fig. 7(b)–(f). Fig. 7(b) shows a maximum amplitude projection (MAP) through the marked blood vessel. For a better visualization of the contrast enhancement when using the different filters with respect to the unfiltered beamforming, a dB-scale was chosen in Fig. 7(b). The performances of the standard deviation, 1/r and sinc filter are very comparable and all of them induce a contrast improvement of about 30 dB when compared with the unfiltered beamforming. Furthermore, the FWHM of the marked blood vessel is decreased from 393 μm (no filter) to 197 μm, 180 μm and 177 μm for the standard deviation filter, the 1/r filter and the sinc filter respectively.

Fig. 7

(a) channel data that is summed up for reconstruction of a pixel corresponding to a vessel (marked with the arrow in Fig. 7(c)–(f) as well as the 1/r and sinc best-fit functions. (b) MAP through the vessel marked with the arrow in Fig. 7(c)–(f). Scale is in dB for better visualization of the contrast improvement resulting of the different filter functions. (c–f) optoacoustic B-scan image at the position marked with the white dotted line in Fig. 7(g)–(j) with the different filter functions. (g–j) projection images (maximum projection in z-direction) of 165 reconstructed frames (with the different filters) with an elevational stepsize of 400 μm.

The improvement of the reconstruction quality becomes obvious in the projection images as well, where the background noise vanishes and only signals from the vasculature remain (Fig. 7(g)–(j)).

Discussion

While the optoacoustic contrast between different tissue types is typically much higher than in conventional ultrasound, the absolute signal amplitudes, and thereby the SNR, typically is relatively low. In addition, despite being advantageous with respect to handling, availability and real-time capabilities, optoacoustic systems relying on medical-type linear array probes typically lead to lower resolution than scanning systems based on single focused ultrasound transducers. We therefore investigated different filter algorithms suitable for linear array optoacoustic beamforming with respect to their ability to improve the SNR and the resolution of reconstructed data sets. We report on a new non-linear filter for optoacoustic delay-and-sum beamforming based on confidence measures for each reconstructed pixel which analyse the amplitude distribution of summed signal samples.

If the angle dependent sensitivity of a transducer and the differences in the acoustic path between a pixel position and a transducer element positions can be neglected, signals summed up during the reconstruction of a pixel corresponding to an actual optoacoustic source should have comparable amplitudes. In more sophisticated approaches, the propagation distance or even the sound field of individual receivers can be considered. Our confidence measures reflecting the likelihood of a pixel for being an optoacoustic source are based on the comparison of the sample amplitude distribution with the theoretical pattern, i.e. a sinc function. The approach was first tested on simulated data from point absorbers. Here, we found that the FWHM can be reduced by more than a factor of 4 when taking the sound field into account and by more than 2 in the other cases. Obviously, this only holds true when the actual sample distribution strictly follows the theoretical assumption. Accordingly, noise has a strong influence on the achievable reconstruction quality (Fig. 5(b)).

Furthermore, when using experimental data, overlapping wave fronts or limited knowledge of the actual sound velocity in the medium will have an influence on the sample distribution and lead to a deviation from theoretical assumptions. With an assumed speed of sound not corresponding to the actual one, the wrong samples would be summed up during beamforming and hence analysed for comparison with a theoretical pattern. Effects that are unfavorable for accurate DAS-beamforming such as the presence of fat layers between the optoacoustic source and the transducer will therefore have a detrimental influence on the filter performance as well.

Given that the phantom setting with thin wires in a water bath can be considered quite close to the simulated setting, further tests were performed using in-vivo data, where the multitude of overlapping wave fronts from different subcutaneous vessels, the unknown sound velocity and the lower SNR introduce further challenges. The analysis shows that subcutaneous structures can be reconstructed with a FWHM being half as wide as in the case of non-filtered time-domain beamforming. However, the differences between the filter functions vanish with increasing noise level. For instance, a subcutaneous vessel in a human finger was reconstructed to a width of 390 μm without filters and between 177 and 197 μm for the suggested filter functions. Given the higher computation costs of the algorithms taking into account the acoustic path or the sound field, which require a best-fit of an 1/r or a sinc function to the sample amplitude distribution respectively, the filter based on the analysis of the standard deviation of the sample distribution seems to be the method of choice for most applications.

Furthermore, it has to be mentioned, that the algorithms compares the experimental data to the pattern generated by a point source. Accordingly, it is very suitable for the enhancement of the background contrast of sources of limited spatial extent, but tends to suppress wide structures, whose wave-fronts obviously do not resemble those of a point source. This can be seen in Fig. 7(c)–(f), where both the skin surface and the subcutaneous vessels can be seen in the DAS image (Fig. 7(c)). However, after the non-linear filtering, the contrast of vessels is improved while the signal from the skin surface is suppressed. This has to be taken into account when applying the new algorithm on in-vivo data. While it can be very useful in optoacoustic angiography, where the aim is a high-contrast visualization of the vasculature, it might affect the overall image fidelity when it comes to imaging of larger structures such as organ boundaries. Potential applications of the algorithm would be those where high contrast and high resolution detection and imaging of vasculature is needed. In rheumatoid arthritis, a condition for which optoacoustic imaging has already demonstrated to have great potential [[38], [39], [40]], inflammation-induced angiogenesis in finger joints is an early symptom. Improving the optoacoustic contrast of such vasculature would allow to detect even smaller vessels and thereby potentially allow the identification of earlier disease symptoms.

Conclusion

The presented non-linear filter for optoacoustic DAS reconstruction leads to a contrast improvement of 30 dB and an improvement of the resolution by a factor of two on experimental in-vivo data from subcutaneous vasculature. The analysis shows, that all variants of the filter (STD, 1/r, sinc) perform in a similar fashion when data are relatively noisy, such as in most optoacoustic in-vivo settings, where the optical fluency is limited due to safety considerations. The sinc-variant of the filter, which compares signal distributions to a sinc-function best fit, provides even further improvements with respect to SNR and resolution if the noise level of the input pre-beamformed channel data is low enough. For instance, on experimental data from a wire phantom, the FWHM of the PSF was reduced to about 180 μm using the sinc-fit against 220 to 240 μm for the 1/r- and STD filter respectively (and 480 μm for conventional DAS). However, in most settings, especially in in-vivo ones with typically lower SNR in the pre-beamformed input data, the STD-filter should be the method of choice due to lower computational costs for similar performance. Furthermore, it has to be kept in mind that the filter compares signal patterns to those of optoacoustic point sources and therefore is best suited for amplification of small sources, such as vessels in the context of optoacoustic angiography.

Conflict of interest

The authors declare that there are no conflicts of interest.