Photoacoustic imaging beyond the acoustic diffraction-limit with dynamic speckle illumination and sparse joint support recovery

In deep tissue photoacoustic imaging the spatial resolution is inherently limited by the acoustic wavelength. Recently, it was demonstrated that it is possible to surpass the acoustic diffraction limit by analyzing fluctuations in a set of photoacoustic images obtained under unknown speckle illumination patterns. Here, we purpose an approach to boost reconstruction fidelity and resolution, while reducing the number of acquired images by utilizing a compressed sensing computational reconstruction framework. The approach takes into account prior knowledge of the system response and sparsity of the target structure. We provide proof of principle experiments of the approach and demonstrate that improved performance is obtained when both speckle fluctuations and object priors are used. We numerically study the expected performance as a function of the measurements signal to noise ratio and sample spatial-sparsity. The presented reconstruction framework can be applied to analyze existing photoacoustic experimental datasets containing dynamic fluctuations.


Introduction
Optical microscopy is an invaluable tool in biomedical investigation and clinical diagnostics. However, its use is limited to depths of not more than a fraction of a millimeter inside tissue due to light scattering. At depths beyond a few hundred microns light scattering in tissue prevents the ability to focus light to its diffraction limit. While non-optical imaging techniques, employing non-ionizing radiation such as ultrasound, allow deeper investigations, they typically possess inferior resolution and generally do not permit microscopic studies of cellular structures at depths of more than a millimeter 1 . One of the leading approaches for deep tissue optical imaging is photoacoustic imaging/tomography 2, 3 . Photoacoustic imaging relies on the generation of ultrasonic waves by absorption of light in a target structure under pulsed optical illumination. Ultrasonic waves are produced via thermo-elastic stress generation, and propagate to an externally placed ultrasonic detector-array without being scattered. Photoacoustic imaging thus provides images of optical contrast with a spatial resolution limited by acoustic diffraction. Ultimately, the ultrasound resolution in soft tissue is limited by the attenuation of high frequency ultrasonic waves. As a result, the depth-to-resolution ratio of deep-tissue photoacoustic imaging is ~100-200 in practice 2, 3 . For example, at a depth of 5 mm one can expect a resolution of around 20 m at best, more than an order of magnitude above the optical diffraction limit.
Recently, it has been demonstrated that the conventional acoustic diffraction-limit can be overcome by exploiting temporal fluctuations in photoacoustic signals originating from illuminating the sample with dynamically varying optical speckle patterns 4 . This work was inspired by the notion of super-resolution optical fluctuation imaging (SOFI) 5 , developed for fluorescence microscopy. In SOFI, a resolution beyond the diffraction limit is obtained via highorder statistical analysis of temporal fluctuations of fluorescence, recorded in a sequence of images. To apply SOFI in photoacoustics a set of random, unknown, optical speckle illumination patterns was used as a source of fluctuations for super-resolution photoacoustic imaging 4 ( Fig.1a-b). Using this approach, an effective resolution enhancement of 1.6 beyond the acoustic diffraction-limit was obtained by analyzing temporal fluctuations second moment (variance) using a set of 100 photoacoustic images. Obtaining higher resolution by analyzing higher statistical moments with such a limited number of images results in strong artifacts due to statistical noise, caused by the insufficient number of analyzed frames (Fig.1b, right inset).
Here, we show that by adapting an advanced computational reconstruction algorithm based on a compressed-sensing framework it is possible to obtain an enhancement in resolution and reconstruction fidelity in photoacoustic imaging beyond that possible with the basic statistical fluctuations analysis of SOFI 4 , while using the same experimentally obtained dataset (Fig.1c). Specifically, we recognize that photoacoustic imaging under dynamic unknown speckle illumination 4, 6 is an instance of blind structured illumination microscopy (blind-SIM) 7,8 . Since the photoacoustic signal generation and detection is a linear process, reconstructing the target object from the measured set of photoacoustic images is formulated as a linear inverse problem (Fig.1c), studied at depth in many other instances of imaging in optics and other domains [8][9][10][11][12] . In principle, a reconstruction approach for solving such inverse problems should exploit all available information. This includes, in addition to the acquired images and detection system response, any prior information on the statistics or structure of the unknown illumination patterns, the non-negativity of the illumination intensity and the object absorption, and any inherent structural correlations or sparsity.
In this work, we employ a reconstruction approach based on compressed sensing (CS) 13-15 . CS has been demonstrated to enable super-resolved optical imaging of microscopic structures 11, 12, 16 , and imaging using sub-Nyquist sampling 14, 17 , i.e. imaging using a number of measurements that is lower than the number of image pixels. Our use of CS for super-resolved photoacoustics combines the high-resolution information contained in the temporal fluctuations, with the super-resolution recovery capability of CS to retrieve the maximum amount of information using a minimum number of acquired photoacoustic frames. Unlike the use of dynamic speckle illumination for enhanced resolution and optical sectioning in optical microscopy 7,8,18,19,[20][21][22] , in photoacoustics the optical speckle grain size is not limited by the same diffraction-limit as the imaging point spread function (PSF), which is acoustic. In practice, the speckle grain size can be orders of magnitude smaller than the PSF dimensions 6 . This suggests that the resolution increase is not limited by the usual factor of two as in structuredillumination microscopy 7 , even without nonlinearities in the imaging process 23 . Given the differences in dimensions, the speckle grains can be thought of as playing a similar role to blinking sources of signal ('molecules') with dimensions much smaller than that of the imaging system PSF. A situation analogous to the one considered in SOFI 5 , PALM 24 and STORM 25 super-resolution microscopy techniques utilizing blinking fluorescent molecules. In the presented approach a compressed-sensing reconstruction algorithm takes advantage of available prior information on the object structure (here, sparsity of the absorbing structure, given by the diagonal matrix D), acoustic system response (given by the convolution matrix H), and optical speckle properties (the matrix U whose columns are the illumination speckle patterns) to perform the recovery. The result is a reconstruction with improved fidelity and resolution. The reconstructions in (b-c) are numerical results obtained using 2000 speckle patterns with measurement noise of 5% of signal peak.

Principle
The principle of our approach is presented in Fig 1.a,c. We consider the conventional photoacoustic imaging experimental setup given in Fig.1a, employing an ultrasound transducerarray for detection, and a pulsed laser source with sufficiently large coherence length to produce speckles. A rotating diffuser in the illumination path produces random unknown optical speckle illumination patterns on the absorbing object that we wish to image. For each of the m=1..L unknown speckle intensity patterns, u m (x,z), a single photoacoustic image, y m (x,z), is measured, where x and z are spatial coordinates (Fig.1a). For simplicity of the analysis we consider a linear transducer array which images a small two-dimensional absorbing structure. The structure spatial absorption pattern is given by o(x,z). We assume the object lies at the center of the transducer array FOV, effectively having a shift-invariant acoustic PSF 4 . Under these considered assumptions, each acquired photoacoustic image, y m (x,z), is a convolution of the acoustic detection PSF, given by h(x,z) with the object structure, o(x,z), multiplied by the unknown speckle illumination pattern intensity u m (x,z), which generates the photoacoustic signals: In principle, the same analysis can be performed in three dimensions and for shift-varying PSF. The imaging challenge is thus to find o(x,z) given the known/measured system PSF h(x,z), and the set of photoacoustic images y m (x,z), without knowing the speckle patterns u x (x,z). Importantly, the photoacoustic images are of considerably lower resolution (spatial frequency bandwidth) than both the object and speckle patterns due to the convolution with the acoustic PSF h(x,z). While in conventional structured illumination and ghost-imaging 17 the speckle illumination patterns are known and the reconstruction is straightforward, here, as in blind-SIM the speckle patterns are unknown. However, many of the speckle patterns properties, such as their non-negativity, intensity statistics and correlations, are universal and can be used in the reconstruction algorithm 7 .
While using SOFI to analyze the N th -order statistical cumulant of y m yields, in principle, a N resolution increase without deconvolution 5 , it is accompanied by strong artifacts when an insufficient number of frames is available 5 (Fig.1b, rightmost inset). However, since SOFI's simple statistical analysis does not take into account all available information besides the temporal fluctuations, its performance can be surpassed through a model-based approach that considers prior knowledge of the object and the system. This is exactly the design goal of CS: to recover the maximum amount of information from a minimal number of measurements. In a nutshell, a CS algorithm can solve a set of underdetermined linear equations, such as those given by equation 1, by exploiting the inherent sparsity of natural objects in an appropriate transform basis. Remarkably, such sparsity is a general property of most natural images, and is at the core of modern lossy image compression algorithms, such as JPEG 14 . Here we exploit the sparsity of the object O.
To establish a CS framework for our problem, we formulate the problem that is given by equation (1) in a continuous coordinate space by its representation in an adequately sampled discrete space: The intensity pattern exciting photoacoustic signals is then given by a vector u m where its entries represent the intensity of the m th illumination pattern at all relevant spatial positions (x,z). The acoustic spatial emission pattern can then be written as: vm=Dum, where D=diag (O 1 ,…,O N ) is a non-negative diagonal matrix representing the object pattern on an N pixels grid. In the case where the object pattern is sparse in real space, D has a sparse diagonal. Each photoacoustic image ym is a result of a convolution of the photoacoustic emission pattern with the acoustic detection system PSF, given by the vector h. We denote by H its convolution matrix. With these notations the m th measured photoacoustic image can be written as (see Fig  1c.): Stacking the entire series of measurements ym as columns in a matrix Y, and um as columns in a matrix U, the entire image set acquisition process can be written as (see Fig 1c.): where Y and H are the measured photoacoustic image-set and system response, correspondingly, and are assumed to be known. The matrices D and U are the unknown object absorption pattern (on the diagonal of D) and the unknown speckle illumination patterns (as columns of U).
Let X=DU. Since D is diagonal and sparse, the support of X and D will be the same, where the support is equal to the nonzero rows of a matrix. Thus, our recovery problem can be stated as: argmin || ||, . . 0, is row sparse.
Once we recover X=DU, we can obtain the support of X and thus of D. Since for fully developed speckles the ensemble average of the speckle intensity is the same for all spatial positions 7, 26 , D can be found by averaging the columns of X (average over the illumination patterns). See detailed explanation of (4) in Methods section.
The challenge of recovering a common support from the acquisition of multiple correlated sparse signals is referred to in the CS literature as the multiple measurement vector (MMV) problem 15, 27 . Several methods have been proposed for estimating the support 15 . Here we use a Bayesian approach to perform the recovery 8, 28, 29 , referred to as multiple sparse Bayesian learning (M-SBL). We chose M-SBL since it led to superior reconstruction fidelity compared to other MMV methods we have tested. The M-SBL algorithm implements a maximum aposteriori estimate (MAP) to find the optimal X while defining some constraints on X so as to encourage solutions which match the prior knowledge. In our setting, we employed a spatial sparsity prior (see Methods).

Experimental results
To experimentally demonstrate the advantage of the proposed reconstruction strategy over conventional photoacoustic reconstruction and statistical based fluctuations analysis, we have used the above described M-SBL algorithm to analyze a set of experimental photoacoustic images of test samples made of absorbing beads of diameters 50m-100m, measured under dynamic speckle illumination 4 . The experiments have been performed using the experimental system sketched in Fig.1a, and described in detail in reference 4 (see Methods).
The experimental results are presented in Fig.2. The leftmost column of Fig.2 presents the microscopic optical image of the absorbing beads for three different samples as imaged directly without the presence of any scattering, which can be considered as the 'ground truth' of o(x,z). The other four columns present images reconstructed from photoacoustic acquired data: (1) a conventional photoacoustic image, given by the mean of the photoacoustic image set ∑ ; (2) a 2 nd order SOFI fluctuation analysis followed by a Richardson-Lucy deconvolution using the squared PSF for deconvolution. These may be considered as the best results of the photoacoustic SOFI approach, as achieved in reference 4; (3) an M-SBL reconstruction using as input only the mean photoacoustic image, i.e. providing the increase in resolution relying only on sparsity without additional information from speckle fluctuations; and (4) an M-SBL reconstruction using the entire image dataset, providing the main result of this work. One may consider the results of SOFI as exploiting only the temporal fluctuations of the signal, the results for M-SBL on the standard photoacoustic image as optimal recovery using only the sparsity prior, and the rightmost column as exploiting simultaneously the sparsity priors, the common support and the fluctuations information for all of the speckle realizations. As expected, exploiting more information yields superior reconstruction fidelity, recovering most accurately the number and positions of the absorbing beads, and reducing imaging artifacts. Since the PSF used to form the deconvolution matrix, H, was measured using beads of size similar to the imaged beads (see Methods), some of the reconstructed beads appear smaller than their real size.

Expected performance as a function of experimental parameters
The resolution enhancement of the proposed nonlinear recovery approach depends on many experimental parameters. While it depends most critically on the size of the acoustic PSF (given laterally by the acoustic diffraction limit and axially by the transducer impulse response), it is also highly sensitive to the experimental signal to noise ratio (SNR), the absorbing sample structure/sparsity, and the speckle grain dimensions compared to the PSF dimensions. Qualitatively, the best performance would be expected for the narrowest PSF, highest measurement SNR and the sparsest object. In addition, a smaller speckle grain size will result in a lower fluctuations to mean signals ratio, and thus a lower SNR for resolving the fluctuations 6 .
To quantitatively analyze the expected reconstruction fidelity as a function of the above parameters we have numerically investigated a large set of imaging scenarios involving different PSF size, object sparsity, and SNR. The results of this study are presented in Fig 3.  Figure 3a-b display the correlation between the reconstructed images and the object for each of the considered scenarios. The vertical axis represents the sample sparsity/complexity, taken here as the number of absorbing pixels contained in the area enclosed by the PSF, where the pixel size is taken as the optical diffraction limit (a speckle grain dimension). The horizontal axis provides a measure of the PSF size, taken as the ratio between the width of the PSF and the optical diffraction limit. All simulations were performed with a pixel grid of 70 by 70 pixels, with a pixel size equal to the speckle grain dimensions, i.e. no structures with dimensions below the optical diffraction-limit are considered. The PSF used in the simulations was generated by simulating the acoustic response of a 50um bead being uniformly illuminated by a 1ns laser pulse and recorded by a linear transducer array with 256 elements, upper frequency limit of 8MHz, inter-element pitch of 0.125mm, and element width of 0.125mm.
As expected from the intuitive qualitative description given above, both the PSF dimensions and the object sparsity play a crucial role in obtaining a high fidelity reconstruction. One can immediately appreciate that all of the results in the top row in the plots of Fig.3a-b display a near perfect reconstruction. This result is not very surprising when considering that the top row presents very sparse samples that contain nearly a single absorber inside a resolution cell given by the PSF dimensions. This is close to the scenario considered in localization microscopy techniques such as PALM/STORM. In these cases, a centroid analysis can provide a good estimation of the absorber location.
To illustrate what recovery errors appear when the reconstructed images fail to perfectly recover the object structure, we provide in Fig.3l-o several examples for reconstructed images side by side with the original objects (h-k), for several cases presented in Fig.3b. It can be seen that when the concentration of absorbers is too high (or sparsity too low) the algorithm fails to identify their exact positions but instead delivers some continuous line-like structure connecting them, which can be understood as blurring of the original image. This gradual 'break-down' of the recovery algorithm is encouraging for practical imaging purposes, as even the reconstructions with low calculated correlation to the object carry relevant information on the object. This is expected to be advantageous when continuous structures such as blood vessels are considered. While different measures for the sample sparsity, SNR, and PSF size relative to the object structure can be chosen, we expect the graphs of Fig.3a-b to serve as a reference to the expected performance given a specific experimental imaging scenario.

Discussion
We have proposed an advanced reconstruction algorithm for photoacoustic imaging which efficiently exploits dynamic temporal fluctuations, joint sparse support constraints, and known system response for improved resolution and reconstruction fidelity. Our approach provides superior performance compared to the recently proposed SOFI-based photoacoustic speckle fluctuation analysis 4 .
We formulated the photoacoustic imaging problem in the case of dynamic speckle illumination as an instance of blind structured illumination. As such, other algorithms that were developed to solve this problem could be employed and their performance compared to the specific algorithm used here. From an estimation theory point-of-view, it is clear that improved performance will be obtained for a reconstruction algorithm that takes into account all available information. Here, we have used CS to take into account some of this information by exploiting the object sparsity and the multiple random measurements provided by the random speckle illumination. Improved algorithms could be developed by incorporating also the non-negativity of the speckle intensity and object structure, and the known universal Rayleigh statistics of fully developed speckles 26 , which we have used only implicitly here to reconstruct the object from the reconstructed matrix X (see Methods).
We have demonstrated our approach using two-dimensional objects and sparsity constraints in real space, however one may consider applying our proposed compressed-sensing recovery approach to reconstruct three-dimensional (3D) objects utilizing sparsity in any other 3D-sparse transform basis representations, which better matches the object structure, e.g. wavelet or minimum total variance (min-TV). In this work we have made use of the measured system PSF. However, when the system PSF is not measured it may be possible to develop an advanced reconstruction algorithm that estimates the PSF and the object simultaneously, as is done in blind deconvolution 30 .
As noted earlier, an important practical challenge for applying the approach for deep tissue photoacoustic imaging arises from the large difference between the speckle grain dimensions and the acoustic PSF dimensions. At large imaging depths the speckle grain dimensions, which are given by the optical diffraction limit, would be orders of magnitude smaller than the imaging acoustic PSF 6 . In this scenario, the measured value at each pixel in the raw photoacoustic frames is the result of a sum over a large number of fluctuating uncorrelated speckle grains (the PSF convolution kernel being much larger than the speckle grain). This results in an overall small fluctuation to mean value in each image pixel between the different frames 6 . Resolving the small fluctuations over a large background may be challenging under low SNR conditions. Choosing a long optical wavelength and a high ultrasound frequency would be advantageous for this task.

Experimental setup
The setup used to perform the experiments is drawn schematically in Fig.1. The beam of a nanosecond pulsed laser (Continuum Surelite II-10, 532 nm wavelength, 5 ns pulse duration, 10 Hz repetition rate) was focused on a ground glass diffuser (Thorlabs, 220 grit, no significant ballistic transmission). The scattered light illuminated a 2-D absorbing sample embedded in an agarose gel block. This phantom was located 5 cm away from the diffuser, leading to a measured speckle grain size of 30 µm. The absorbing sample was placed in the imaging plane of a linear ultrasound array (Vermon, 4 MHz center frequency, >60% bandwidth, 128 elements, 0.33 mm pitch), connected to an ultrasound scanner (Aixplorer, Supersonic Imagine, 128-channel simultaneous acquisition at 60 MS/s). A collection of black polyethylene microspheres (Cospheric, 50 µm and 100 µm in diameter) was used to fabricate phantoms with isotropic emitters. The PSF of the system was measured for each sample by concentrating light on one single 50 µm-diameter bead. The diffuser was removed from the light path during this step, to ensure a homogenous illumination of the bead.
For each sample, a set of photoacoustic signals for 100 uncorrelated speckle patterns was obtained by rotating the diffuser. Special care was taken to reduce sources of fluctuations other than the multiple speckle illumination between photoacoustic acquisitions. The raw recorded (RF) acoustic signals were processed by a low-pass filter with a sharp cutoff that eliminated all frequencies above 2.4MHz in sample 1, and 2.8MHz and 5.3MHz for samples 2 and 3 correspondingly. Different cutoff frequencies were used to create more challenging recovery scenarios for the different algorithms. The dimensions of the resulting PSF, defined as the full width at half maximum (FWHM) after the low-pass filtering were 613/1643um, 537/1231um and 393/562um in the transverse/axial directions for samples 1, 2 and 3 correspondingly.
For each sample, 100 photoacoustic images were reconstructed from the raw acoustic signals, for each of the 100 speckle patterns, using a time-domain backprojection algorithm on a grid of 814 by 814 pixels with a pixel pitch of 25um. The time domain backprojection reconstruction is based on summing the photoacoustic signals taken at appropriate retarded times 31 . The reconstructed photoacoustic images ym were downsampled to half of their original size by bilinear interpolation, to reduce the required computational resources and run time.

Recovery algorithm
The MMV recovery algorithm we used in this work is M-SBL 28, 29 . As mentioned above, the M-SBL algorithm implements a maximum a posteriori estimator. Through the application of the Bayes law it searches for the value of which maximizes the joint probability of , , see details below. The M-SBL algorithm we have adapted to use in this work assumes fluctuations having Gaussian statistical distribution with zero mean 29 . Dynamic speckle illumination having a speckle grain size that is considerably smaller than the reconstruction grid provides indeed a Gaussian statistical distribution for the temporal fluctuations (as a direct consequence of the central limit theorem and the large number of summed speckle grains in each pixel), but with a mean that is not zero. Thus, when running M-SBL for the MMV case, the pixel-wise calculated temporal mean of the fluctuation images was subtracted from them pixel-wise. This ensures that the prior zero mean Gaussian distribution of X used in M-SBL matches the provided measured data. Improving the algorithm by employing a prior containing the exponential intensity statistics of fully developed speckles may increase the reconstruction fidelity. To provide a fair comparison with the 2 nd order SOFI reconstruction, the M-SBL reconstructions presented in Figures 1-2 show the variance over each row of the recovered matrix X. Providing the temporal variance at each reconstructed image pixel. In Figure 3 the displayed M-SBL results are the standard deviation of each reconstructed image pixel, since the standard deviation provides a measure that is linearly related with the mean absorption in each spatial position (pixel). In the case where the reconstruction grid pixel size is smaller than the speckle grain dimensions, the standard deviation of each pixel provides a quantitative estimate of the mean absorption since for the exponential statistics of fully developed speckle 26 the temporal standard deviation is equal to the mean. Since the speckles fluctuations are uncorrelated, in the case where the grid pixel size is larger than the speckle grain dimensions, the standard deviation provides the mean absorption times the square root of the number of speckles contained in the absorbing area inside the reconstructed pixel.
The reconstruction M-SBL algorithm we have used in this work is an adaptation of the algorithm of Zhang et al. 29 . Briefly, under certain prior assumptions of Gaussian distributions of signal and noise, the posterior density of the j-th column of X becomes 28 :