We present a method for enhancing the spatial resolution of optical-absorption-based photoacoustic imaging through or within highly scattering media. The optical speckle pattern that emerges as light propagates through diffuse media provides structured illumination to an object placed behind a scattering wall. The photoacoustic signal produced by such illumination is detected using a focused ultrasound transducer. We demonstrate through both simulation and experiment that, by acquiring multiple photoacoustic images, each produced by a different random and unknown speckle pattern, an image of an absorbing object can be reconstructed with a spatial resolution far exceeding that of the ultrasound transducer. The imaging approach and reconstruction method exploiting joint sparsity will be valuable for biomedical applications and other fields where high-frequency structural information is scrambled by diffusive processes.
© 2016 Optical Society of America
Optical imaging of objects inside or behind turbid media is a challenging endeavor, as multiple scattering scrambles the object information and hides features from view. While techniques such as optical coherence tomography  and multiphoton microscopy  have been developed to extract ballistic light from scattered light, thus enabling diffraction limited imaging, these approaches are suitable only for imaging depths of the order of 5 times the scattering mean free path, beyond which the ballistic photon population is too weak with respect to the scattered light to provide useful information. Recently, there have been significant developments in wave front shaping techniques that use scattered light to image objects in turbid media [3–8]. Here, light is manipulated using a spatial light modulator or nonlinear optical crystal such that scattering is compensated. If the transmission matrix of the sample is unknown, as is often the case in dynamic scattering media, focusing and imaging using wave front shaping requires a feedback signal to guide the optical field to a given location . Noninvasive feedback approaches have been developed that combine optics with acoustics, exploiting either acousto-optic interactions [10–12], where an ultrasound field is used to frequency modulate light in the interaction region, or the photoacoustic effect, where the intensity of light at a given location is inferred by the photoacoustic signal amplitude [13–17]. Wave front shaping in living tissue is also complicated by the fact that it must be done very quickly, as the speckle decorrelation time is of the order of milliseconds. While techniques for rapid wave front correction continue to evolve, imaging approaches that eliminate the need for implicit or explicit determination of the transmission matrix and wave front correction may prove advantageous both in terms of reducing system complexity and potentially increasing speed.
In this paper, we propose a feedback-free imaging method that uses the random optical speckle patterns that naturally emerge as light propagates through strongly scattering media as a structured illumination source for photoacoustic imaging. Our approach, termed blind structured illumination photoacoustic microscopy (BSIPAM), was inspired by recent work in fluorescence microscopy where super-resolution imaging was demonstrated using multiple unknown speckle illumination patterns [18,19]. We extend this concept to the multiple scattering domain using photoacoustics, with the speckle pattern serving to generate ultrasound, which is detected by a focused transducer. In structured illumination microscopy, the resolution enhancement is limited to about a factor of 2 because the maximum spatial frequency of the illumination pattern is constrained by the optical transfer function of the microscope . In BSIPAM, however, the resolution enhancement can potentially be higher. The lateral resolution in conventional acoustic-resolution photoacoustic microscopy at depth in scattering media is determined by the point spread function (PSF) of the transducer. If we consider structured illumination of an object using a speckle pattern, the maximum spatial frequency in the illumination pattern is dictated by the speckle size and can be well beyond that typically accessible by the transducer. However, by structuring the illumination, the high spatial frequencies in the object are downshifted into the transducer passband with the resolution of BSIPAM ultimately limited by the speckle size rather than the PSF of the transducer. To account for the unknown speckle patterns, we develop a new reconstruction algorithm, named block-FISTA, that uses joint sparsity constraints and exploits the fact that all speckle patterns illuminate the same absorber distribution.
We note that while speckle patterns have been used to enhance photoacoustic imaging in the past, the approaches are markedly different than that proposed here. Multiple speckle illumination was initially used to enhance the visibility in photoacoustic imaging and address limited view issues by introducing an additional smaller granular structure caused by speckles . In addition, taking multiple photoacoustic images using different speckle patterns and imaging the variance of the photoacoustic response rather than the mean shows a resolution enhancement of about 1.4 or 1.6 with subsequent deconvolution . The technique is conceptually similar to super-resolution optical fluctuation imaging . Our approach is more closely related to structured illumination microscopy, where the physical origin of the resolution enhancement is frequency mixing between the illumination and object spatial frequencies .
Next we formally describe how to recover the absorber distribution at a spatial resolution close to the speckle size.
Suppose that different speckle patterns are used. After proper discretization, we assume that the speckle patterns and the absorber distribution are represented by discrete vectors , where the components denote values at equidistant points in the imaging domain. The measured photoacoustic data are then given by1). The product is the sound source corresponding to the th speckle pattern. Because the intensity distribution and the speckle patterns are unknown, the system in Eq. (1) consists of equations for unknown scalar parameters and is therefore underdetermined. On the other hand, the sound sources are (theoretically) uniquely determined by the deconvolution equations for . However, due to the ill-conditioned nature of deconvolution with a smooth kernel, these uncoupled equations are sensitive to errors and further provide only low-resolution reconstructions when solved independently and without proper regularization. In this work we propose to use a block-sparsity penalty to obtain high-resolution reconstructions.
It is now well accepted that exploitation of sparsity allows super-resolution signal recovery. As opposed to classical linear models, sparse recovery algorithms assume the estimated signal to be a superposition of a few (but unknown) elements from a large collection of base signals containing high-frequency as well as low-frequency information. If the performed measurements allow differentiating sufficiently well between the base signals, appropriate sparse recovery algorithms can recover sparse combinations of base elements [24,25]. A naive approach incorporating sparsity would be to look, separately, for the sparsest solution of any individual equation . However, such an approach still misses the main feature of the reconstruction problem, namely, that all products come from the same density distribution . Such a joint-sparsity (or block-sparsity) situation can be implemented by using the joint sparsity term2) consists of the -norm of the -norms of the pixel values over all blocks . The inner -norm couples the different sound sources, whereas the outer -norm accounts for the required sparsity. By exploiting the joint-support assumption, the -term has higher super-resolution capacity than the pure -term . In Supplement 1 we develop an efficient reconstruction algorithm, the block-FISTA, that realizes such a joint-sparsity approach for solving Eq. (1).
3. NUMERICAL RESULTS
We now provide an example of how BSIPAM can be implemented to improve the resolution of photoacoustic imaging. We start with a simple one-dimensional sample [constant along the vertical axis in Fig. 1(a)] where the absorber distribution ρ is a set of eight lines of 10 μm thickness. The distance between lines varies between 40 and 150 μm. The absorbers are illuminated with the speckle pattern shown in Fig. 1(b), where the speckle size is 25 μm. The product of the speckle intensity and absorber distribution gives the spatial distribution of the photoacoustic sources shown in Fig. 1(c). Next, we assume that the absorber lies in the focal plane of an ultrasonic transducer with a Gaussian PSF (full width at half-maximum ). When scanning the transducer across the imaging plane, the photoacoustic response is found by convolving the source distribution with the transducer PSF, and the result is shown in Fig. 1(d), where only lines with the largest separations are clearly resolved.
The experiment is now repeated times by applying a new random speckle pattern each time and using the procedure above to calculate the photoacoustic response to a line scan across the absorber distribution [a single horizontal line from Fig. 1(d)]. Six such line scans, each from a different speckle pattern, are shown in Fig. 2(a). It is apparent that the structured illumination patterns, via the random speckle illumination, influence the line scan response. It is precisely this variation that is exploited for super-resolution imaging. Figure 2(b) shows a comparison of several techniques to reconstruct the absorber distribution. Gaussian noise with a standard deviation of 1% of the photoacoustic signal amplitude was added to the data prior to reconstruction. The mean response corresponds to the photoacoustic response obtained through averaging all of the line scans. The variance reconstruction is obtained through plotting the square root of the signal variance at each spatial position. Given a Gaussian PSF, we expect this to result in a increase in the spatial resolution over the mean response and the improvement is evident in the simulation . The Richardson–Lucy deconvolution (RLD) plot shows the deconvolution of the mean response given the associated PSF. Finally, the BSIPAM response shows the result obtained through our block-FISTA algorithm. Here, even the smallest feature spacing (40 μm) is resolved giving a clear resolution advantage over the other approaches. Additional one-dimensional simulation results and discussion regarding the maximum achievable resolution and required number of speckle patterns are included in Supplement 1.
Our reconstruction strategy is directly applicable to two-dimensional image reconstruction, as well. To investigate the capabilities of our block-FISTA algorithm for two-dimensional image reconstruction, we considered a star as an absorber consisting of several lines, as shown in Figure 3. We use a sample of size where the closest distance of the lines is 17 μm. The PSF has been taken as two-dimensional Gaussian kernel with a FWHM of 35 μm. In this case, Gaussian noise with a standard deviation of 10% of the photoacoustic signal amplitude was added to the data prior to reconstruction. In the mean image the individual lines clearly cannot be separated. If we define the resolution as the distance where two lines are clearly separated, then the resolution in the mean image is close to 42 μm. Also, after applying a (regularized) deconvolution (as in Ref. ), the lines cannot be separated, although the resolution increases by a factor 1.4. To further increase the resolution, we apply the proposed approach and simulate random speckle patterns with a speckle size of 9 μm. The result of our algorithm clearly separates the lines having a minimal distance of 17 μm. In this case, our approach increases resolution by at least a factor of 2.4.
4. EXPERIMENTAL VALIDATION
We now turn to a discussion of the experimental validation of the proposed approach. A schematic of the experimental setup is depicted in Fig. 4. The sample is placed in a water-filled test tank with transparent walls. Light from a pulsed, frequency-doubled Nd:YAG laser () is sent through a ground glass diffuser and focused onto the sample using a lens with a focal length . The illumination beam diameter is set by an adjustable aperture placed adjacent to the lens. We estimate the speckle size at the focal plane of the lens as , where is the wavelength of the illumination source and is the aperture diameter . Ultrasonic signals are recorded by a focused ultrasonic transducer (Panametrics V3512, 90 MHz, element size 6 mm, focal length 13 mm) connected to an ultrasonic pulser/receiver (Panametrics 5073PR, bandwidth 1 kHz–75 MHz) and sampled with a 12 bit digital oscilloscope (Lecroy HDO6104). A silicon photodetector is used to sample the illumination pulses in order to correct the photoacoustic signal amplitudes for shot-to-shot variations in the laser intensity. The same photodetector is also used to trigger the oscilloscope. The ultrasonic transducer is mounted on a motorized linear scanning stage, which moves the transducer in the y-direction, perpendicular to the sample. A manual linear stage in the z-direction is used for adjusting the distance between sample and transducer. A second motorized linear scanning stage is used to rotate the diffuser via a v-belt.
Photoacoustic signals were acquired with a spatial step size of 10 μm over a total scan length of 1.0 mm. The measured time traces were processed using a Fourier transform, and the magnitude spectra were summed between 15 and 85 MHz to get a net photoacoustic response at each spatial location. Using this processing approach, the system was first characterized by scanning over a single alpaca fiber with a diameter in the 18–25 μm range. To avoid speckle artifacts, the measurement was repeated 17 times at different diffuser positions and the results were averaged to simulate plane wave optical illumination. The mean response was then deconvolved with a step function (25 μm width) approximating the hair to yield the line spread function (LSF) shown in Fig. 5(a). The LSF determined in this manner has a FWHM of approximately 110 μm.
A sample consisting of five alpaca fibers (see Fig. 4 inset) was fabricated, and the relative location of the fibers with respect to each other at the ultrasound scan location (based on the optical image) is shown as the object in Figs. 5(a)–5(c). We note that the absolute position of the fiber array is only an estimate. The aperture diameter was set to 3.0 mm, yielding a speckle size of approximately 10 μm, and 101 scans across the fiber array were acquired. The diffuser was rotated after each scan to provide a unique speckle pattern. The mean photoacoustic response from all scans, simulating uniform illumination, is shown in Fig. 5(a), along with five individual scans acquired with different speckle patterns. The weighting between individual features changes in response to the structured illumination associated with each scan. The mean, RLD, and BSIPAM photoacoustic reconstructions are shown in Fig. 5(b). Here, BSIPAM is the only technique able to clearly resolve all five fibers, and the relative positions of the fibers in the reconstructed response agree well with those measured using optical microscopy. Given the measured LSF, one can just begin to resolve two lines when the line spacing reaches approximately 80 μm. The minimum fiber spacing in our experiment, resolved using BSIPAM, is estimated to be less than 40 μm, giving a resolution enhancement of at least 2.
Figure 5(c) shows a simulation of the experiment. Here, the fiber array we use consists of 18 μm fibers, with the fiber spacing matching that from the optical image. The transducer LSF is calculated for a broadband transducer with a 40 MHz central frequency and a FWHM of 110 μm . Gaussian noise with a standard deviation of 1% of the photoacoustic signal amplitude was added to the data prior to reconstruction (the experimentally measured data had a noise standard deviation of the signal amplitude), and the speckle size was set to 10 μm. One-hundred line scans with different speckle patterns were simulated, and the individual line scans again show fluctuations indicative of the random spackle illumination. The BSIPAM approach allows for resolution of all five lines in agreement with what was found experimentally. We note that extension of the experimental results to two-dimensional imaging presents no additional challenges, and can be facilitated by a transducer array for data collection. Finally, we emphasize that BSIPAM utilizes the fact that the photoacoustic signals generated by different random speckle patterns are indeed different. If the speckle size is too small with respect to the PSF/LSF, then we expect that the differences will be small, as the received signal is integrated over the detection area. On the other hand, large speckle size may lead to large differences in the photoacoustic signals, but will also eliminate the high spatial frequency components and ultimately limit resolution. Additional results showing the speckle size effects on the variation in photoacoustic responses across speckle patterns are included in Supplement 1.
In conclusion, if structural information is scrambled by diffusive or dissipative phenomena, like scattering, absorption, or diffusion, we can use joint sparsity (called also block sparsity) to gain better spatial resolution from blind structured illumination measurements. We have demonstrated super-resolution experimentally in one-dimensional photoacoustic imaging, and numerically also in two-dimensional imaging, using multiple speckle illuminations. The differences in the photoacoustic signals generated with random speckle patterns are utilized in BSIPAM using the block-FISTA reconstruction algorithm. This algorithm reconstructs the absorbing structure from measured PA signals with a resolution close to the speckle size and has the potential for reconstructing structural information from signals measured with other modalities using blind structured illumination, like photothermal imaging or lateral structures in optical coherence tomography.
Austrian Science Fund (FWF) (P25584-N20); European Regional Development Fund (ERDF) (IWB2020: MiCi); Federal Government of Upper Austria; Christian Doppler Research Association.
See Supplement 1 for supporting content.
1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, and C. A. Puliafito, “Optical coherence tomography,” Science 254, 1178–1181 (1991). [CrossRef]
2. F. Helmchen and W. Denk, “Deep tissue two photon microscopy,” Nat. Methods 2, 932–940 (2005). [CrossRef]
3. R. Horstmeyer, H. Ruan, and C. Yang, “Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue,” Nat. Photonics 9, 563–571 (2015). [CrossRef]
4. A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nat. Photonics 6, 283–292 (2012). [CrossRef]
5. Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang, “Optical phase conjugation for turbidity suppression in biological samples,” Nat. Photonics 2, 110–115 (2008). [CrossRef]
6. J. Tang, R. N. Germain, and M. Cui, “Superpenetration optical microscopy by iterative multiphoton adaptive compensation technique,” Proc. Natl. Acad. Sci. USA 109, 8434–8439 (2012). [CrossRef]
7. K. Si, R. Fiolka, and M. Cui, “Fluorescence imaging beyond the ballistic regime by ultrasound pulse guided digital phase conjugation,” Nat. Photonics 6, 657–661 (2012). [CrossRef]
8. I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. 32, 2309–2311 (2007). [CrossRef]
9. S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. 104, 100601 (2010). [CrossRef]
10. X. Xu, H. Liu, and L. V. Wang, “Time-reversed ultrasonically encoded optical focusing into scattering media,” Nat. Photonics 5, 154–157 (2011). [CrossRef]
11. K. Si, R. Fiolka, and M. Cui, “Breaking the spatial resolution barrier via iterative sound-light interaction in deep tissue microscopy,” Sci. Rep. 2, 748 (2012). [CrossRef]
12. B. Judkewitz, Y. M. Wang, R. Horstmeyer, A. Mathy, and C. Yang, “Speckle-scale focusing in the diffusive regime with time reversal of variance-encoded light (TROVE),” Nat. Photonics 7, 300–305 (2013). [CrossRef]
13. F. Kong, R. H. Silverman, L. Liu, P. V. Chitnis, K. K. Lee, and Y. C. Chen, “Photoacoustic-guided convergence of light through optically diffuse media,” Opt. Lett. 36, 2053–2055 (2011). [CrossRef]
14. T. Chaigne, O. Katz, A. C. Boccara, M. Fink, E. Bossy, and S. Gigan, “Controlling light in scattering media non-invasively using the photoacoustic transmission matrix,” Nat. Photonics 8, 58–64 (2014). [CrossRef]
15. 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,” Opt. Express 21, 26671–26676 (2013). [CrossRef]
16. P. Lai, L. Wang, J. W. Tay, and L. V. Wang, “Photoacoustically guided wavefront shaping for enhanced optical focusing in scattering media,” Nat. Photonics 9, 126–132 (2015). [CrossRef]
17. D. B. Conkey, A. M. Caravaca-Aguirre, J. D. Dove, H. Ju, T. W. Murray, and R. Piestun, “Super-resolution photoacoustic imaging through a scattering wall,” Nat. Commun. 6, 7902 (2015). [CrossRef]
18. E. Mudry, K. Belkebir, J. Girard, J. Savatier, E. Le Moal, C. Nicoletti, M. Allain, and A. Sentenac, “Structured illumination microscopy using unknown speckle patterns,” Nat. Photonics 6, 312–315 (2012). [CrossRef]
19. J. Min, J. D. Jang, D. M. Keum, S. W. Ryu, C. H. Choi, K. H. Jeong, and J. C. Ye, “Fluorescent microscopy beyond diffraction limits using speckle illumination and joint support recovery,” Sci. Rep. 3, 2075 (2013). [CrossRef]
20. M. G. L. Gustafsson, “Surpassing the lateral resolution by a factor of two using structured illumination microscopy,” J. Microsc. 198, 82–87 (2000). [CrossRef]
21. J. Gateau, T. Chaigne, O. Katz, S. Gigan, and E. Bossy, “Improving visibility in photoacoustic imaging using dynamic speckle illumination,” Opt. Lett. 38, 5188–5191 (2013). [CrossRef]
22. T. Chaigne, J. Gateau, M. Allain, O. Katz, S. Gigan, A. Sentenac, and E. Bossy, “Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination,” Optica 3, 54–57 (2016). [CrossRef]
23. T. Dertinger, R. Colyer, G. Iyer, S. Weiss, and J. Enderlein, “Fast, background-free, 3D super-resolution optical fluctuation imaging (SOFI),” Proc. Natl. Acad. Sci. USA 106, 22287–22292 (2009). [CrossRef]
24. S. Mallat, A Wavelet Tour of Signal Processing: The Sparse Way (Academic, 2008).
25. J. Yang, J. Wright, T. Huang, and Y. Ma, “Image super-resolution via sparse representation,” IEEE Trans. Image Process. 19, 2861–2873 (2010). [CrossRef]
26. M. Haltmeier, “Block-sparse analysis regularization of ill-posed problems via ℓ1,2-minimization,” in 18th International Conference on Methods and Models in Automation and Robotics (MMAR) (IEEE, 2013) pp. 520–523.
27. J. E. Ward, D. P. Kelly, and J. T. Sheridan, “Three-dimensional speckle size in generalized optical systems with limiting apertures,” J. Opt. Soc. Am. A 26, 1855–1864 (2009). [CrossRef]
28. S. Canumalla, “Resolution of broadband transducers in acoustic microscopy of encapsulated IC’s: transducer selection,” IEEE Trans. Compon. Packag. Technol. 22, 582–592 (1999). [CrossRef]