# Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination

###### Abstract

In deep tissue photoacoustic imaging, the spatial resolution is inherently limited by acoustic diffraction. Moreover, as the ultrasound attenuation increases with frequency, resolution is often traded-off for penetration depth. Here we report on super-resolution photoacoustic imaging by use of multiple speckle illumination. Specifically, we show that the analysis of second-order fluctuations of the photoacoustic images combined with image deconvolution enables resolving optically absorbing structures beyond the acoustic diffraction limit. A resolution increase of almost a factor 2 is demonstrated experimentally. Our method introduces a new framework that could potentially lead to deep tissue photoacoustic imaging with sub-acoustic resolution.

Light scattering prevents standard optical microscopes to obtain well-resolved images deep inside biological tissues. In the past twenty years, photoacoustic (PA) imaging has been developed to overcome this limitation, by imaging optical absorption deep inside strongly scattering tissue with the resolution of ultrasound [25]. PA imaging relies on the unscattered ultrasonic waves emitted by absorbing structures under pulsed illumination via thermo-elastic stress generation. It therefore provides images at depth in tissue with a spatial resolution limited by acoustic diffraction. Ultimately, the ultrasound resolution for biological soft tissue is limited by the attenuation of ultrasound, which typically increases linearly with frequency. As a result, the depth-to-resolution ratio of PA imaging at depth is around 200 in practice [25, 26]. As an illustration, axial resolution down to 5 µm and lateral resolution down to 10 µm have been reached with high frequency acoustic detectors at depth up to 5 mm [27].

In this letter, we demonstrate that the conventional acoustic-diffraction limit in PA imaging may be overcome by exploiting PA signal fluctuations, building on the super-resolution optical fluctuation imaging (SOFI) technique developed for fluorescence microscopy [28]. SOFI is based on the idea that a higher-order statistical analysis of temporal fluctuations caused by fluorescence blinking provides a way to resolve uncorrelated fluorophores within a same diffraction spot. In this work, we introduce multiple optical speckle illumination as a source of fluctuations for super-resolution PA imaging, inspired by the principle introduced in optics with SOFI [28] or from derived approaches using speckle illumination[29]. In PA imaging, multiple speckle illumination was initially introduced by our group as a mean to palliate limited-view or highpass-filtering artefacts [30]. Here, we demonstrate that a second-order analysis of optical speckle-induced PA fluctuations also provides super-resolved PA images beyond the acoustic diffraction limit.

In this work, we consider PA images reconstructed from a set of PA signals measured with an ultrasound array. A conventional backprojection algorithm is used to reconstruct the images, and it is assumed that the reconstructed PA quantity may be written as a convolution:

(1) |

where is the PSF corresponding to the conventional PA imaging process, the distribution of optical absorption and the optical intensity pattern (see Supplement 1, sec. 1.A, for a detailed justification of Eq. 1). Let us now consider that the region of interest is successively illuminated by many different speckle patterns with mean . The following expression for the mean PA image (estimated from averaging the PA images obtained with many realizations of the speckle illumination)

(2) |

shows that the resolution of the reconstructed image is dictated by the spatial frequency content of . Under the assumption that the optical speckle size is much smaller than that of , the variance image for uncorrelated speckles is given by (see Supplement 1, sec. 1.B)

(3) |

The variance image appears as the convolution of the squared object by the squared PSF, which has a higher frequency content than the PSF itself. As a result, the variance image is expected to have a higher resolution as compared to the mean image.

Our objective was to demonstrate experimentally that the measurement of a two-dimensional (2-D) variance image can provide a super-resolution PA image of the absorption distribution beyond the acoustic-diffraction limit. The experimental setup used is shown in Fig. 1.a. 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. Estimates of the PSF were measured using isolated 50 µm diameter microspheres, while ordered patterns to be imaged were formed using 100 µm diameter microspheres. Fig. 1.c shows the PSF corresponding to our imaging setup. The dimensions of the PSF, defined as the full width at half maximum (FWHM) of its envelope, were µm in the transverse (x) direction and µm in the axial (z) direction. It is worth noting that the anisotropy and bipolarity of the PSF are singular features of limited-view and limited-bandwidth PA imaging, and do not have their counterpart in regular all-optical imaging. The measurement and analysis of the PSF is further described in detail in Supplement 1, sec.3.

For each sample, a set of PA images was reconstructed for 100 uncorrelated speckle patterns, obtained by rotating the diffuser. The mean and variance images were then computed on a per-pixel basis. As described in detail in Supplement 1, sec. 2.B, special care was taken to reduce sources of fluctuations other than the multiple speckle illumination between PA acquisitions. For the image reconstruction, a time-domain backprojection algorithm was used (see Supplement 1, sec. 2.C). Square images of 20 mm side were reconstructed on a grid of square pixels (25 µm side). To demonstrate the resolution enhancement of our technique, we first designed the phantom shown in Fig. 2.a. Pairs of 100 µm diameter beads were positioned along the z and x axis. The distances between beads (center to center) were: 120 µm, 140 µm and 200 µm along the z direction (from top to bottom), and 250 µm, 200 µm and 160 µm along the x direction (from left to right).

Fig. 2.b shows the mean PA image of the sample, obtained from averaging over the 100 speckle realizations, and which corresponds to the PA image obtained under homogenous illumination [30]. Fig. 2.c shows the square root of the variance of the same data set to make a fair comparison with the mean image in term of units and resolution. At this stage, no significant resolution enhancement can be seen from the comparison of the mean and variance images, for both of which the pairs of beads are blurred by the convolution with either (Fig. 2.b) or (Fig. 2.c). To demonstrate that the variance image contains sub-acoustic diffraction information thanks to the higher-frequency content of as compared to , and to remove the peculiar side lobes of the PA PSF observed in Fig. 2.b and Fig. 2.c, image deconvolution was performed based on Eqs. 2 and 3 to retrieve the absorption distribution . The mean and variance images were deconvolved respectively by and . In an ideal noise-free situation, deconvolution of an image should retrieve the absorbing object with no resolution limit. However, in the presence of noise, spatial frequencies are accurately measured up to a certain limit set by the signal-to-noise ratio (SNR). An inversion strategy that allow accounting for the presence of noise in the measurement was therefore carried out to perform the deconvolution. Retrieving from measurements of the variance image modeled as (with accouting for the experimental noise) was carried out by the minimization of the following constrained least-square functional:

(4) |

with the Euclidian norm over the image space, and a regularization parameter. The constrained minimizer provides a regularized solution to the deconvolution problem, which in turn defines an estimation of the absorption distribution . For comparison, the exact same approach was applied to retrieve an estimation of from the measurement of the mean image. In practice, the minimization of the functional above required adjusting the regularization parameter in order to obtain a fair resolution v.s. noise trade-off [31, Sec. 5.6], and choosing in the FISTA numerical method used to perform the minimization (see Supplement 1, sec. 1.C for further details on the minimization approach and the corresponding algorithm).

Deconvolved images are shown in Fig. 3. As mentioned above, the deconvolved variance image is square rooted to retrieve the absorption distribution . The parameters of the deconvolution algorithm were . On the deconvolved mean image (Fig. 3.a), we observe that only absorbers separated by the largest distance in each direction are resolved (z direction: 200 µm, x direction: 250 µm). In comparison, the absorbers appear much sharper on the deconvolved variance image (Fig. 3.b) and almost every pair is resolved on this image. The only exception is the upper one, for which the beads are nearly touching. Cross-sections along both directions are shown respectively on Fig. 3.c and d. These cross-sections show the maximum amplitude projection of the PA image along a band around the white dashed lines plotted in Fig. 3.a.

To further investigate the influence of the PSF anisotropy, a second phantom was designed, in which 100 µm diameter beads are equally spaced from each other with a center-to-center distance of about 150 µm and distributed on a quarter circle. The corresponding results are shown in Fig. 4. The parameters of the deconvolution algorithm were . Not a single bead can be resolved on the deconvolved mean image (Fig. 4.a). In contrast, the 6 beads can easily be distinguished on the deconvolved variance image (Fig. 4.b). The last two beads on the right are still distinguishable but tend to smear, according to the lower resolution in the transverse x direction. Because of the PSF anisotropy (a consequence of the limited view), we observe that the resolution strongly depends on the angle of the beads pair with respect to the probe surface.

Speckle-induced PA fluctuation imaging is an original method recently proposed by our group to enhance visibility in PA imaging. In this study we demonstrated that fluctuation imaging can also be used to overcome the ultrasound diffraction limit in PA tomography. We showed that this method can separate small absorbers unresolved under standard uniform illumination. This is to the best of our knowledge the first super-resolution PA imaging method which does not require optical focusing, and therefore is applicable beyond the ballistic regime. Super resolution PA microscopy was reported in previous studies, but it relied on light focusing and on nonlinear absorption mechanisms, a situation which cannot be translated to deep tissue PA imaging [32, 33].

The common approach to obtain high resolution at depth in PA imaging was to date to improve the acoustic detection by employing higher frequencies and increasing the detection aperture. Here, we demonstrated that additional efforts can be deployed towards the illumination, in order to induce fluctuations and separate discrete absorbers below the acoustic diffraction limit.

Deconvolution was found essential to reveal the super-resolution properties with the second-order fluctuation images. Deconvolution has already been considered in PA imaging, as part of the reconstruction algorithm [34] and to compensate for smearing induced by the spatial impulse response of the finite-size detectors [35]. However, no super-resolution was yet demonstrated , since there were no physical mechanism extending the high spatial frequency information. Our approach goes beyond past works by considering fluctuations in PA images as a signal that reveal higher spatial frequencies above the noise level. A gain in resolution was obtained in both transverse and axial directions, with a resolution enhancement of at least 1.7 in both directions.

The deconvolution approach was implemented with absorbers of size similar to that of the absorbers used to measure the point spread function. As an immediate drawback, deconvolution was unable to restore the actual size of the absorbers, which lead on the deconvolved image to reconstructed beads smaller than their real size. Nonetheless, this method showed a very good sub-acoustic resolution performance, which was the objective of our work. Although the beads were expected to be almost equally absorbing, a difference in the reconstructed amplitudes was noticed. The proposed method was therefore not shown to be quantitative, which could be attributed to the non-linear deconvolution scheme. Further work should be carried out to investigate the possibility to retrieve quantitative absorption information.

Our proof-of-concept experiments were carried out at low medical ultrasound frequency to simplify the controlled fabrication of absorbing samples (see Supplement 1, sec. 2.A). However, the approach could be scaled down using high frequency detectors, and could also be extended to 3D PA imaging. It must be emphasized that speckle illumination is achievable at depth using a seeded laser [36], whose long coherence length still provides speckle formation after several centimeters of propagation in scattering tissue. PA fluctuations induced by speckle illumination are expected to decrease as the number of speckle grains contained in the absorbers increases [30]. Detection of speckle-induced fluctuations would therefore be challenging deep inside tissue where the speckle grains are on the order of half the optical wavelength, and will require special care on the excitation and detection stability. In the past few years, several techniques combining shaped coherent illumination and PA imaging have been developed. These methods aim at focusing light inside tissue, possibly with sub-acoustic resolution, which could also lead to super-resolved PA imaging [37, 38, 39]. However, these methods require expensive hardware and tedious experimental procedures. We proposed here a very simple imaging technique than does not require any costly equipment and nearly no optical alignment. For biological applications, tissue-induced temporal decorrelation of speckle patterns could even be exploited as a source of fluctuating illumination [40]. Alternatively, this super-resolution method can be extended to fluctuation of the absorption induced by blinking or switchable [41] contrast agents.

#### Acknowledgements.

This work was funded by the European Research Council (grant 278025), the Fondation Pierre-Gilles de Gennes (grant FPGG031) and by the LABEX WIFI (Laboratory of Excellence ANR-10-LABX-24) within the French Program “Investments for the Future” under reference ANR-10- IDEX-0001-02 PSL*. O.K. acknowledges the support of the Marie Curie Intra-European Fellowship. The authors thank Laurent Bourdieu and Jean-François Léger for their helpful assistance on the milling machine.

#### Supplementary materials.

See Supplement 1 for supporting content.

## References

- [1] P. Beard, “Biomedical photoacoustic imaging,” Interface focus 1, 602–631 (2011).
- [2] L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335, 1458–1462 (2012).
- [3] M. Omar, D. Soliman, J. Gateau, and V. Ntziachristos, “Ultrawideband reflection-mode optoacoustic mesoscopy,” Optics letters 39, 3911–3914 (2014).
- [4] T. Dertinger, R. Colyer, G. Iyer, S. Weiss, and J. Enderlein, “Fast, background-free, 3d super-resolution optical fluctuation imaging (sofi),” Proceedings of the National Academy of Sciences 106, 22287–22292 (2009).
- [5] J.-E. Oh, Y.-W. Cho, G. Scarcelli, and Y.-H. Kim, ‘‘Sub-rayleigh imaging via speckle illumination,” Optics letters 38, 682–684 (2013).
- [6] J. Gateau, T. Chaigne, O. Katz, S. Gigan, and E. Bossy, “Improving visibility in photoacoustic imaging using dynamic speckle illumination,” Optics Letters 38, 5188–5191 (2013).
- [7] M. Bertero and P. Boccacci, Introduction to inverse problems in imaging (CRC press, 1998).
- [8] B. Rao, K. Maslov, A. Danielli, R. Chen, K. K. Shung, Q. Zhou, and L. V. Wang, “Real-time four-dimensional optical-resolution photoacoustic microscopy with au nanoparticle-assisted subdiffraction-limit resolution,” Optics letters 36, 1137–1139 (2011).
- [9] L. Wang, C. Zhang, and L. V. Wang, “Grueneisen relaxation photoacoustic microscopy,” Physical review letters 113, 174301 (2014).
- [10] C. Zhang, C. Li, and L. V. Wang, “Fast and robust deconvolution-based image reconstruction for photoacoustic tomography in circular geometry: experimental validation,” Photonics Journal, IEEE 2, 57–66 (2010).
- [11] H. Roitner, M. Haltmeier, R. Nuster, D. P. O’Leary, T. Berer, G. Paltauf, H. Grün, and P. Burgholzer, “Deblurring algorithms accounting for the finite detector size in photoacoustic tomography,” Journal of biomedical optics 19, 056011–056011 (2014).
- [12] J. Bjorkholm and H. Danielmeyer, “Frequency control of a pulsed optical parametric oscillator by radiation injection,” Applied Physics Letters 15, 171–173 (1969).
- [13] T. Chaigne, O. Katz, A. Boccara, M. Fink, E. Bossy, and S. Gigan, “Controlling light in scattering media non-invasively using the photoacoustic transmission matrix,” Nat Photon 8, 58–64 (2014).
- [14] A. M. Caravaca-Aguirre, D. B. Conkey, J. D. Dove, H. Ju, T. W. Murray, and R. Piestun, “High contrast three-dimensional photoacoustic imaging through scattering media by localized optical fluence enhancement,” Optics express 21, 26671–26676 (2013).
- [15] P. Lai, L. Wang, J. W. Tay, and L. V. Wang, ‘‘Photoacoustically guided wavefront shaping for enhanced optical focusing in scattering media,” Nature photonics 9, 126–132 (2015).
- [16] M. Jang, H. Ruan, I. M. Vellekoop, B. Judkewitz, E. Chung, and C. Yang, “Relation between speckle decorrelation and optical phase conjugation (opc)-based turbidity suppression through dynamic scattering media: a study on in vivo mouse skin,” Biomedical optics express 6, 72–85 (2015).
- [17] K. K. Ng, M. Shakiba, E. Huynh, R. A. Weersink, A. Roxin, B. C. Wilson, and G. Zheng, “Stimuli-responsive photoacoustic nanoswitch for in vivo sensing applications,” ACS nano 8, 8363–8373 (2014).
- [18] M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Review of Scientific Instruments 77, 041101 (2006).
- [19] A. Rosenthal, V. Ntziachristos, and D. Razansky, “Acoustic inversion in optoacoustic tomography: A review,” Current medical imaging reviews 9, 318 (2013).
- [20] J. W. Goodman, Speckle phenomena in optics: theory and applications (Roberts and Company Publishers, 2007).
- [21] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM journal on imaging sciences 2, 183–202 (2009).
- [22] D. P. Bertsekas, Nonlinear programming (Athena Scientific, Belmont, ma, USA, 1999), 2nd ed.
- [23] M. Bertero and P. Boccacci, “A simple method for the reduction of boundary effects in the richardson-lucy approach to image deconvolution,” Astronomy & Astrophysics 437, 369–374 (2005).
- [24] M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Physical Review E 71, 016706 (2005).

Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination: supplementary material

This document provides supplementary information to “Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination”. The first section contains theoretical informations: the conditions are given to model the photoacoustic images reconstruction as a convolution, the variance image is shown to result from the convolution of the squared PSF, and details on the deconvolution algorithm used are given. The second section provides more detailed information on the experimental methods, including samples preparation, measurements and processing of photoacoustic signals, and the backprojection algorithm used for the photoacoustic reconstruction. The last section analyze some properties of the photoacoustic point spread function.

## Appendix A Theory

### a.1 Photoacoustic imaging as a convolution process

In its simplest form, the generation and propagation of photoacoustic (PA) pressure waves , generated by a distribution of optical absorption illuminated by a pulsed light with intensity , may be described in a homogenous acoustic medium (with speed of sound and Grüneisen coefficient ) by the following equation

(S1) |

Eq.S1 shows that for a given temporal pulse shape , the pressure field is linearly related to the distribution of optical absorption times the optical intensity pattern. Various algorithms exist that aim at reconstructing from a set of PA signals measured at various points located on some boundary around the sample to be imaged [42, 43]. In the context of our work, we consider a conventional back-projection algorithm based on summing values of the PA signals taken at appropriate retarded times [42]. As such, the reconstructed images remain linearly related to . Moreover, if this linear reconstruction process may also be considered as translation-invariant, i.e. a spatial translation of the object simply translates the reconstructed object, then the reconstruction process may be written as a convolution as assumed by Eq. 1 of our work. Because the PA signals are measured on some boundary around the region to be imaged, the translation invariance cannot be strictly verified. However, for a field of view small enough such that every point that it contains is reconstructed with the same point spread function (PSF), it becomes possible to reasonably assume that the reconstruction process is translation-invariant. In our case, this assumption was validated by measuring PSF at four different locations in the field of view. The four PSF appeared identical (see sec. C.2 below), and the results shown in the manuscript were independant of the particular PSF that was used to perfom the deconvolutions, thus validating our assumption that the reconstruction process may be written as a convolution.

### a.2 The variance image as a convolution with

From (valid only for objects lying in a small enough field of view, see sec. A.1 above), the variance image reads

where

is the autocovariance of the intensity fluctuation in the speckle patterns. Under the common assumption that the speckle is wide-sense stationary [44, p.67], the autocovariance function is a function of only one variable . This function is sharply peaked around its center, and its characteristic length is by definition the speckle grain size [44, p.72]. If the speckle grains are small enough compared to the PSF, may be considered in the double integral above proportional to a delta function , yielding the following convolution expression:

This expression is strictly equivalent to that found for the second-order analysis in SOFI [28], apart from the origin of the fluctuations.

### a.3 Deconvolution algorithm: FISTA

As introduced in the main manuscript, deconvolution of the variance image was performed by minimizing the following constrained least-square functional:

(S2) |

From a practical viewpoint, an iterative minimization algorithm is required for the numerical evaluation of . Since is a strictly convex functional if , its global minimizer is asymptotically reached by any locally convergent iteration, whatever the initial-guess of the algorithm. For this constrained optimization task, a natural candidate is the standard projected-gradient (PG) since its computational burden is very low. The PG is however rather slow to converge and the FISTA iteration [45] that achieves a faster convergence is presently considered as a popular alternative and was used in our work. In comparison to the standard projected-gradient algorithm [46], the FISTA method [45] only requires the additional storage of an auxiliary component . Let be an initial-guess and , the FISTA update for is:

(S3) |

with , the projection operator over , and the gradient of the penalized least-square functional above that reads

(S4) |

where stands for the adjoint of the operator . It is worth noting that all the convolution operations that appear above can be computed efficiently in the Fourier domain via the fast Fourier transform (FFT) and appropriate boundary conditions [47]. Finally, the global convergence of (S3) is granted provided that the constant step-size is adjusted within the interval with

where stands for the Fourier transform of .

## Appendix B Experimental methods

### b.1 Samples preparation

A collection of black polyethylene microspheres (Cospheric, 50 µm and 100 µm in diameter) was used to fabricate phantoms with isotropic emitters. Estimates of the PSF were measured using 50 µm diameter microspheres, while ordered patterns to be imaged were formed using 100 µm diameter microspheres. To precisely control the relative position of the beads, melted gel was first poured on a mold drilled with micrometer precision (Mini Mill, Minitech). The beads were then manually placed in the molded holes of the solidified gel.

### b.2 Measurements and signal processing

To avoid other sources of fluctuation apart from the multiple speckle illumination, appropriate care was taken to eliminate noise and triggering-induced temporal jitters. The intensity of the laser pulses was monitored with a photodiode to compensate for the laser pulse-to-pulse energy fluctuations. In addition, for each speckle illumination (i.e. each diffuser position), the PA signals were averaged over 25 laser pulses to improve the signal-to-noise ratio. The signals were then filtered between 1 MHz and 8 MHz ( order butterworth filter) for noise removal outside of the array bandwidth. A set of PA images was reconstructed for 100 uncorrelated speckle patterns. The mean and variance images of this data set were then computed on a per-pixel basis. The same procedure was repeated with one single speckle realization (static diffuser). The resulting variance image was then subtracted from the previous variance image (with rotating diffuser). This procedure was found to correct for systematic variance noise.

### b.3 Backprojection algorithm

In the context of our work, we considered a conventional back-projection algorithm based on summing values of the PA signals taken at appropriate retarded times. More specifically, the backprojection algorithm described in [48, eq. (20)] was implemented while keeping only the first time-derivative of the processed signal, assuming point-like detectors located in the centers of the elements of the ultrasound array.

## Appendix C Photoacoustic point spread function

### c.1 Measurement

The PSF of the imaging setup was measured by concentrating light on a single 50 µm diameter bead located in the vicinity of the structured 100 µm bead pattern (see Fig. 2.a). This ensured that the 50 µm microsphere was the only PA source. The diffuser was removed from the light path during this step.

### c.2 Translation-invariance of the point spread function across the field-of-view

Following the measurement approach described above, PSFs were reconstructed for 4 different locations in the field of view to confirm that it can be assumed as constant, as required to model the reconstruction by a convolution process. The corresponding results are summarized in Fig. S1. With our current ultrasonic detection and the back-projection algorithm, we showed that this approach was at least valid for a 3 mm x 3 mm field-of-view. On the deconvolved images, a slight difference in the PSF shapes would results in additional artefacts, for instance some rebounds around the reconstructed objects that are located far from the place where the PSF is measured.

### c.3 Analysis

In Fig. S2, we display the PSF (Fig. S2.a) of the PA system and its envelope (Fig. S2.b, computed using Hilbert transform), and the squared PSF (Fig. S2.c). The conventional resolution of the imaging system (when using uniform illumination) is given by the full width at half maximum (FWHM) of the enveloppe of the PSF. The FWHM measured on the data shown on Fig. S2.b was µm in the transverse x direction and µm in the axial z direction. On Fig. S2.c, we observe that the squared PSF is sharper than the PSF itself, which is expected to lead to a higher measurable frequency content in the fluctuation image. The spatial frequency content of the PSF and the squared PSF of the PA imaging setup are shown in Fig. S3. We observe that for a given noise level, higher spatial frequencies are measurable with the squared PSF, hence on the variance image.

## References

- [25] P. Beard, “Biomedical photoacoustic imaging,” Interface focus 1, 602–631 (2011).
- [26] L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335, 1458–1462 (2012).
- [27] M. Omar, D. Soliman, J. Gateau, and V. Ntziachristos, “Ultrawideband reflection-mode optoacoustic mesoscopy,” Optics letters 39, 3911–3914 (2014).
- [28] T. Dertinger, R. Colyer, G. Iyer, S. Weiss, and J. Enderlein, “Fast, background-free, 3d super-resolution optical fluctuation imaging (sofi),” Proceedings of the National Academy of Sciences 106, 22287–22292 (2009).
- [29] J.-E. Oh, Y.-W. Cho, G. Scarcelli, and Y.-H. Kim, ‘‘Sub-rayleigh imaging via speckle illumination,” Optics letters 38, 682–684 (2013).
- [30] J. Gateau, T. Chaigne, O. Katz, S. Gigan, and E. Bossy, “Improving visibility in photoacoustic imaging using dynamic speckle illumination,” Optics Letters 38, 5188–5191 (2013).
- [31] M. Bertero and P. Boccacci, Introduction to inverse problems in imaging (CRC press, 1998).
- [32] B. Rao, K. Maslov, A. Danielli, R. Chen, K. K. Shung, Q. Zhou, and L. V. Wang, “Real-time four-dimensional optical-resolution photoacoustic microscopy with au nanoparticle-assisted subdiffraction-limit resolution,” Optics letters 36, 1137–1139 (2011).
- [33] L. Wang, C. Zhang, and L. V. Wang, “Grueneisen relaxation photoacoustic microscopy,” Physical review letters 113, 174301 (2014).
- [34] C. Zhang, C. Li, and L. V. Wang, “Fast and robust deconvolution-based image reconstruction for photoacoustic tomography in circular geometry: experimental validation,” Photonics Journal, IEEE 2, 57–66 (2010).
- [35] H. Roitner, M. Haltmeier, R. Nuster, D. P. O’Leary, T. Berer, G. Paltauf, H. Grün, and P. Burgholzer, “Deblurring algorithms accounting for the finite detector size in photoacoustic tomography,” Journal of biomedical optics 19, 056011–056011 (2014).
- [36] J. Bjorkholm and H. Danielmeyer, “Frequency control of a pulsed optical parametric oscillator by radiation injection,” Applied Physics Letters 15, 171–173 (1969).
- [37] T. Chaigne, O. Katz, A. Boccara, M. Fink, E. Bossy, and S. Gigan, “Controlling light in scattering media non-invasively using the photoacoustic transmission matrix,” Nat Photon 8, 58–64 (2014).
- [38] A. M. Caravaca-Aguirre, D. B. Conkey, J. D. Dove, H. Ju, T. W. Murray, and R. Piestun, “High contrast three-dimensional photoacoustic imaging through scattering media by localized optical fluence enhancement,” Optics express 21, 26671–26676 (2013).
- [39] P. Lai, L. Wang, J. W. Tay, and L. V. Wang, ‘‘Photoacoustically guided wavefront shaping for enhanced optical focusing in scattering media,” Nature photonics 9, 126–132 (2015).
- [40] M. Jang, H. Ruan, I. M. Vellekoop, B. Judkewitz, E. Chung, and C. Yang, “Relation between speckle decorrelation and optical phase conjugation (opc)-based turbidity suppression through dynamic scattering media: a study on in vivo mouse skin,” Biomedical optics express 6, 72–85 (2015).
- [41] K. K. Ng, M. Shakiba, E. Huynh, R. A. Weersink, A. Roxin, B. C. Wilson, and G. Zheng, “Stimuli-responsive photoacoustic nanoswitch for in vivo sensing applications,” ACS nano 8, 8363–8373 (2014).
- [42] M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Review of Scientific Instruments 77, 041101 (2006).
- [43] A. Rosenthal, V. Ntziachristos, and D. Razansky, “Acoustic inversion in optoacoustic tomography: A review,” Current medical imaging reviews 9, 318 (2013).
- [44] J. W. Goodman, Speckle phenomena in optics: theory and applications (Roberts and Company Publishers, 2007).
- [45] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM journal on imaging sciences 2, 183–202 (2009).
- [46] D. P. Bertsekas, Nonlinear programming (Athena Scientific, Belmont, ma, USA, 1999), 2nd ed.
- [47] M. Bertero and P. Boccacci, “A simple method for the reduction of boundary effects in the richardson-lucy approach to image deconvolution,” Astronomy & Astrophysics 437, 369–374 (2005).
- [48] M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Physical Review E 71, 016706 (2005).