## Abstract

Acoustic scattering medium is a fundamental challenge for photoacoustic imaging. In this study, we reveal the different coherent properties of the scattering photoacoustic waves and the direct photoacoustic waves in a matrix form. Direct waves show a particular coherence on the antidiagonals of the matrix, whereas scattering waves do not. Based on this property, a correlation matrix filter combining with a time reversal operator is proposed to preserve the direct waves and recover the image behind a scattering layer. Both numerical simulations and photoacoustic imaging experiments demonstrate that the proposed approach effectively increases the image contrast and decreases the background speckles in a scattering medium. This study might improve the quality of photoacoustic imaging in an acoustic scattering environment and extend its applications.

© 2017 Optical Society of America

## Corrections

13 September 2017: A typographical correction was made to the author affiliations.

## 1. Introduction

Photoacoustic imaging is a hybrid biomedical imaging technique, which inherits the rich contrast of optical imaging and the high spatial resolution of ultrasonography in deep tissue [1,2]. During a photoacoustic imaging procedure, a pulse laser is illuminated on the region of interest (ROI). Optical absorbers in the ROI absorb the electromagnetic energy and generate ultrasound due to the photoacoustic effect. Image of the acoustic sources, i.e., the optical absorbers, can be reconstructed from the detected photoacoustic signals. Photoacoustic imaging has shown great potentials in biomedical imaging [3,4], such as osteoarthritis assessment [5], vasculature visualization [6–8], *in vivo* brain monitoring [9], drug delivery monitoring [10], tumor detection [11], guidance of high-intensity focused ultrasound [12,13], and so on [14].

Classical photoacoustic imaging methods are usually based on the principle of coherent beamforming. An ultrasound array detects acoustic signals from the surrounding medium. The measured wave could contain both contributions of the direct wave and the scattering wave. The direct wave propagates from the source to the sensor directly. Taking advantage of the deterministic signal phases, the acoustic source can be located or imaged by achieving coherent beamforming. Particularly, based on the signal coherence, Wang et al. significantly improve the spatial resolution of photoacoustic imaging by using an adaptive coherence weighting technique [15]. However, the scattering wave undergoes one or several scattering events before reaching the sensor. The randomness of scattering breaks the coherence in signals, brings speckles, leads to a markedly decreased penetration depth, and even makes imaging fail [1,16,17]. Therefore, the coherent imaging reconstruction techniques rely on an assumption of homogeneous media where direct wave dominates. Although classical detection methods are capable of detecting optical absorption, they primarily depend on the assumption of homogeneous acoustical properties. Unfortunately, acoustical inhomogeneity widely appears such as in various biological organs. The presence of bones, air cavities, or microbubbles may cause strong acoustic scattering, which may degrade image quality [18]. Acoustic scattering is a nightmare for classical photoacoustic imaging techniques.

This issue has received considerable attention in the past decade. Many efforts have been done to break through the limitation of imaging in acoustic inhomogeneous medium. The statistical reconstruction method is presented to improve image quality [19], which relies on a priori knowledge on the location of acoustic distortions. More approaches that are usually suitable for some specific cases have been demonstrated then, such as time reversal method [20,21], and coherence factor optimization [22], interferometry method [23], and so on [24–26]. For a more generation situation, a photoacoustic wave passing through an acoustically random scattering layer would be seriously distorted. A universal scheme to reconstruct images from the distorted and noise-like signals remains a challenging.

In wave physics especially, the form of the matrix is particularly appropriate to describe the wave transmission and to effectively reveal the inherent deterministic coherence of the acoustic signal in the scattering medium. By treating an inhomogeneous medium as one realization of a random process, some aspects of random matrix theory have been fruitfully applied to optical imaging [27–29], acoustic backscattering imaging [30–32], and telecommunication [33–35] in complex media. Especially, Aubry and Derode [30,36] organize the impulse responses in a matrix form and separate the single scattered echo of the target from the multiple scattering background. Then, the acoustic backscattering images are reconstructed by using a time reversal operator. Their experimental results show that this method is particularly suited to the detection of a target embedded in a strongly scattering medium in comparison with classical imaging techniques based on the Born approximation. These studies open up a new way to imaging the scattering medium.

In this study, inspired by the previous studies in the ultrasonic imaging of scattering medium [30,36], we investigate passively detected wide-band photoacoustic waves propagating through complex media by using a matrix form, where direct and scattering wave coexists. We demonstrate that the intrinsic coherence of the direct wave can be revealed in the matrix. Based on this finding, a correlation filter is proposed to reduce the scattering components in the detected photoacoustic signals. As a result, the accuracy and single-to-noise of the imaging are improved. The experiments are used to validate the proposed method.

## 2. Method

As shown in Fig. 1, in order to avoid heavy numerical calculations and without loss of generality, we assume a 2D model in this study. Under the illumination of pulse laser, the optical absorbers in the region of interest (ROI) radiate acoustic waves to the surrounding medium. After the acoustic waves are propagated through a random scattering layer, a passive ultrasound array with *N* elements detects the ultrasound signal. The recorded *N*-channel signals are arranged in a vector **H**(*t*) = [*H*_{1}(*t*), *H*_{2}(*t*), …, *H _{n}*(

*t*), …,

*H*(

_{N}*t*)] with

*n*= 1, 2, …

*N*. After a Fourier transform, the recorded time domain signals can be rewritten in their frequency domain

**P**(

*T*,

*f*) = [

*P*

_{1}(

*T*,

*f*),

*P*

_{2}(

*T*,

*f*), …,

*P*(

_{n}*T*,

*f*), …,

*P*(

_{N}*T*,

*f*)], where

*P*(

_{n}*T*,

*f*) is the short-time Fourier transform of the

*n*-th channel signal in the time window [

*T*– Δ

*t*/2,

*T*+ Δ

*t*/2].

*T*is the time of flight and Δ

*t*is the width of the time window.

**P**(*T*, *f*) contains both the direct wave **P**^{D}(*T*, *f*) from the targets and the scattering wave **P**^{S}(*T*, *f*). That is, *P _{n}*(

*T*,

*f*) =

*P*(

_{n}^{D}*T*,

*f*) +

*P*(

_{n}^{S}*T*,

*f*). In a two-dimensional configuration, under the paraxial approximation, the direct wave

*P*(

_{n}^{D}*T*,

*f*) corresponding to the acoustic source at (

*X*,

*Z*) can be approximated as follows,

*A*

_{0}=

*Z*

^{–1/2},

*k*= 2π

*f*/

*c*is the wave number,

*c*is the speed of sound,

*x*is the

_{n}*x*-coordinate of the

*n*-th element in the array, and

*Z*=

*cT*corresponds to the depth of photoacoustic sources, as shown in Fig. 1. The scattering wave

*P*(

_{n}^{S}*T*,

*f*) is a sum of waves propagating through

*L*scattering paths in the scattering media,

_{n}*A*and

_{l}*s*are the amplitude and phase relating to the

_{l}*l*-th propagating path. For the simple in form, we have extracted the term exp(

*jkZ*) from the term exp(

*jks*) with

_{l}*S*=

_{l}*s*–

_{l}*Z.*Since the nature of random scatterer distribution and scattering paths, all of the amplitude

*A*, the phase

_{l}*S*and

_{l}*s*are random.

_{l}To image a heterogeneous medium, our basic idea is to reduce the influence of the scattering wave. A response matrix **K**(*T*, *f*) can be constructed as **K** = **PP**^{T}. Substituting Eqs. (1-2) into **K** = **PP**^{T}, we can give an insight into the **K** matrix,

*D*= 2

_{m}*Z*+ (

*x*–

_{m}*X*)

^{2}/2

*Z*and

*D*= 2

_{n}*Z*+ (

*x*–

_{n}*X*)

^{2}/2

*Z*.

Observing the Eq. (3), we can find that the **K** matrix contains two components. The first term in the right of the equation, named *K ^{C}_{m,n}*, does not depend on the scatterer distribution and the random scattering path. This implies that whatever the realization of disorder, there is a deterministic phase relation between the elements

*K*in the matrix

_{m,n}**K**located on the same antidiagonal. Considering the couples of the (

*m*–

*q*)-th and (

*m*+

*q*)-th (

*q*= 1, 2, …,

*N*) channels, which correspond to the elements on the

*m*-th antidiagonal, it has

*x*

_{(}

_{m–q}_{)}+

*x*

_{(}

_{m + q}_{)}= 2

*x*and

_{m}*x*

_{(}

_{m + q}_{)}

*– x*

_{(}

_{m–q}_{)}= 2

*qw*, where

*w*is the array pitch. It can be demonstrated that

*qw*, but not the acoustic source or the scattering paths. The direct waves from the optical absorbers manifest themselves as a particular coherence on the antidiagonals of the matrix. Whereas, the left terms within the braces contain a random phase

*S*, which is related to the random scatterer distribution and the random scattering path. This random phase does not have any coherence in these random terms

*K*.

^{R}_{m,n}Compared to other coherence-based algorithms [15], this coherent property in matrix form allows us to extract the direct waves from the detected signals by using a correlation filter [36,37] and it is effective in a scattering medium. The filtering process can be briefly described as follow: Rotating the matrix **K**(*T*, *f*) 45 degrees with the coordinate transformation (*x _{m}*,

*x*) → (

_{n}*y*,

_{u}*y*): ${y}_{u}=({x}_{m}-{x}_{n})/\sqrt{2}$, ${y}_{v}=({x}_{m}+{x}_{n})/\sqrt{2}$, we can generate a new matrix

_{v}**A**(

*T*,

*f*) =

**A**

^{C}(

*T*,

*f*) +

**A**

^{R}(

*T*,

*f*) with a dimension

*L*, where the matrices

**A**

^{C}and

**A**

^{R}represent the direct contribution and the scattering contribution, respectively. After rotation, the coherence of direct waves appears along the columns of

**A**

^{C}. Then, a filtered matrix

**A**

^{F}can be obtained by projection

**A**

^{F}=

**CC**

^{†}

**A**, where the column vector

**C**with the components

*c*= exp[

_{u}*jky*

_{u}^{2}/2

*Z*] is the characteristic space of direct waves. Its components read:

**K**

^{F}(

*T*,

*f*) with the original coordinates can be obtained by applying the coordinate inversion. In the new

**K**

^{F}matrix, the random components have been reduced without changing the coherence part.

Finally, the image of the optical absorbers can be recovered from the filtered matrix **K**^{F} by using a time reversal operator: the singular value decomposition is applied to the matrix **K**^{F}(*T*, *f*) = **U**^{F}**Λ**^{F}**V**^{F†}. Λ is a diagonal matrix whose non-zero elements *λ _{i}^{F}* are called the singular values of

**K**

^{F}. Each optical absorber is mainly associated to one non-zero singular [38]. Therefore, by numerically back-propagating the singular vector, an image at Z =

*cT*can be achieved by the time reversal operator

**I**

_{i}^{F}(

*Z*,

*f*) =

*λ*

_{i}^{F}|**G***

*V*|, where

_{i}^{F}**G**with the component ${g}_{m}=\mathrm{exp}(jk{r}_{m})/\sqrt{{r}_{m}}$is the propagating operator in a homogeneous media between the array and focal plane and

*r*is the distance between the

_{m}*m*-th element in the array and the point in the focal plane. Since the incoherent scattering components in

**K**

^{F}have been reduced, the quality of the recovered image will be significantly improved.

The method proposed relies on the paraxial approximation and the absence of aberration. Yet, the paraxial approximation is sometimes not valid, in particular at shallow depth. A potential solution of this problem could be to add a time-delay layer in front of the targets in order to increase the equivalent depth of the medium. In addition, ultrasound signals often suffer from aberrations due to the inhomogeneities of the wave speed in the medium. In this situation, the matrix filter is still effective [36]. Moreover, the coherent-weight technique [15] or speed correction techniques [22] could also be involved to reduce the influence of aberrations.

## 3. Results

#### 3.1 Numerical Simulations

Numerical simulations are utilized to validate the proposed method. As shown in Fig. 1, a passive ultrasonic array with 101 elements and an array pitch of 0.5 mm is carried out to pick up the photoacoustic signals from the surrounding environment. The acoustic scattering layer consists of 30 randomly distributed optical absorbers with a diameter of *d _{1}* = 0.8 mm are located at a distance of

*a*= 40 mm from the ultrasonic array. The thickness of the scattering layer is

*L*= 20 mm. The concentration of the scatterers is 3 rods/cm

^{2}. The frequency-averaged scattering mean free path

*l*is 27.9 ± 1.33 mm for this medium between 1.26 and 2.68 MHz [39]. Whereas, the speed of sound and the density of these optical absorbers are 5200 m/s and 7870 kg/m

_{s}^{3}, respectively. Two targets (speed of sound 1500 m/s and the density 1000 kg/m

^{3}) with a diameter of

*d*= 0.6 mm in ROI are placed behind the scattering layer. These acoustic sources (in the scattering layer and in the ROI) emit ultrasonic pulses with a central frequency of 2.0 MHz and a - 3 dB bandwidth of 1.26 ~2.68 MHz to the surrounding medium. The surrounding medium has the speed of sound 1500 m/s and the density 1000 kg/m

_{2}^{3}, which is closed to water or soft tissues. The strong impedance mismatch will generate complex acoustic scattering. The photoacoustic signals are detected by the ultrasonic array and recorded with a sampling frequency of 20 MHz. The ROI in the dashed box is the furthest away from the array. This region corresponds to the worst signals polluted by random acoustic scattering. Therefore, in the following discussion, we will mainly focus on imaging the two targets in the ROI.

Figure 2 gives the results of the simulation with the acoustic thickness (*L*/*l _{s}*) of 0.7168. The time window length (

*∆t*) is set to 5 μs to promise that the direct waves within ROI are in the same time window [40]. Figure 2(a) presents the photoacoustic field detected by the passive array. As shown, the wavefronts are fragmentized and appear many noise-like patterns because of the scattering. Figure 2(b) gives the

**K**(

*T*,

*f*) matrix at time

*T*= 32 μs [the window 1 in Fig. 2(a)] and the frequency

*f*= 2.0 MHz. This region is closest to the array. Although there is still some interference of the scattered waves, the direct waves play a leading role. As the above theory has predicted, the K matrix displays a deterministic coherence along its antidiagonals, as indicated by the arrow marks in Fig. 2(b). Figures 2(c) and 2(d) respectively present the

**K**(

*T*,

*f*) and

**K**(

^{F}*T*,

*f*) matrix corresponding to the ROI in this study, where time

*T*= 50 μs [the window 2 in Fig. 2(a)] and the frequency

*f*= 2.0 MHz. The array is far from this region and the scattered waves dominate. As a result, the matrix

**K**shows random pattern [Fig. 2(c)]. However, after the matrix filtering, the matrix

**K**

^{F}displays a deterministic coherence along its antidiagonals again, which implies the reduced scattering components. We can quantify the reduction of the scattering components by calculating the average background intensities in the ROI. The percentage of the reduction between Fig. 2(c) and Fig. 2(d) is 48% and it can reach 61% while the acoustic thickness

*L*/

*l*

_{s}increases to 1.1947.

We use the first singular value to reconstruct the image because it contains 96.53% of the total energy for the target. The better photoacoustic images could be reconstructed from the filtered matrix, as shown in Fig. 3. And to show comparison, Figs. 3(a) and 3(b) are the images obtained by the classic delay-and-sum beamforming method and the proposed method, respectively. A deconvolution algorithm has been applied to reconstruct the image by the beamforming method and the image is plotted with a normalized scale. For the classic beamforming method, the image has low contrast and shows strong speckles in the background because of the strong scattering. The proposed method provides a much better image than the traditional one. The background is very clear since the scattered waves are reduced. The image has a good contrast.

Figures 4(a) and 4(b) quantify the comparison of the signal-to-noise ratio (SNR) and the resolution (axial and transverse) of these two methods. The axial and transverse resolution is estimated by extracting the diameter of the targets from the full width at half-maximum (FWHM) of the image. More numerical simulations are utilized to compare the proposed method with the classical beamforming method. We keep the thickness of the scattering layer *L* = 20 and the number of array elements *N* = 101 constant and change the number of the scatterers from 10 to 50 in the scattering layer. The corresponding concentration of the scatterers varies from 1 rods/cm^{2} to 5 rods/cm^{2}. And the frequency-averaged scattering mean free paths are 83.75 ± 4.00 mm, 41.88 ± 2.00 mm, 27.9 ± 1.33 mm, 20.9 ± 1.00 mm, 16.75 ± 0.80 mm with a corresponding range of frequency from1.26MHz to 2.68 MHz, respectively.

In the image, the signal region is estimated considering the half maximum contour of the images of two optical absorbers under the condition of no scattering. The mean image intensity in the two regions (the length of the axial direction is 2.2 mm, the length of the transverse direction is 5 mm, the centers of the two regions are at the actual central position of the optical absorbers) is considered as the signal intensity, and the noise intensity is defined by the mean value in the background region of the ROI. The image SNR is obtained as a ration of the signal intensity and the noise intensity. As shown in the Fig. 4(a), the horizontal axis is the acoustic thickness of the scattering layer. The scattering becomes stronger as the acoustic thickness increases. When the scattering is weak, both methods provide high SNR. With the increase of the acoustic thickness, the image SNR of the traditional beamforming methods is decreased due to the interference of scattering waves while the SNR of the proposed method is always higher than that of the traditional beamforming, which means the proposed method could provide high image quality in scattering medium.

Figure 4(b) compares the resolution of the proposed method and the beamforming method. Both the axial and transverse resolution becomes worse with the increase of the acoustic thickness. The axial resolution of the proposed method is slightly worse than the beamforming method. For photoacoustic imaging, the axial resolution depends on the bandwidth of the detected signals. The larger axial resolution of the proposed method might be related to the loss of a major part of the initial frequency bandwidth. The target is detected over a frequency bandwidth narrower than that of the emitted signal. Besides, the information overlap in the axial direction caused by a large time window may also lead to a bad axial resolution. Whereas, the beam-forming method shows a bad transverse resolution, because the limited-view scanning will lead to the strong distortion along the transverse direction [41].

The other parameters in Figs. 4(c) and 4(d) are same as the parameters described in the Fig. 1. We only change the number of the array elements to 17, 33, 49, 101, respectively. (The position of the first and the last transducers remain the same, changing the pitch between the transducer.) The SNR of both methods can be improved by simply adding the number of the transducers. When the number of transducers is large enough, the increasing trend of the SNR is no longer as obvious as before. [Fig. 4(c)]. As shown in the Fig. 4(d), the number of the array elements has limited influence on the resolution of photoacoustic imaging for both of the mentioned ways. Given a certain number of the transducers, the proposed method can obtain a better transverse resolution than the classical beamforming method, while the axial resolution of the two methods maintains almost the same. The above simulation results validate the proposed theory and the imaging method.

In practical situations, acoustic scattering may also cause wave mode conversion, which turns the longitudinal wave to shear wave. Since the direct waves does not undergo acoustic scattering and the mode conversion, their coherence on the antidiagonals of the matrix are always valid. While for the scattering waves, the mode conversion during scattering will not bring additional coherence to them. Therefore, the mode conversion will not affect the theory.

#### 3.2 Experiment

Finally, experiments have also been done to examine the performance of the proposed method in the real situation. Our photoacoustic experimental setup is given in Fig. 5(a). Two pencil leads with a diameter of 0.5 mm in the ROI are irradiated by pulse laser with a wavelength of 532 nm and a pulse width of about 8 ns. After absorbing the laser energy, the pencil leads emit ultrasound to the surround media. The pencil leads are served as acoustic sources in the experiments.

A scattering layer consists of 100 randomly distributed steel rods with a diameter of 0.8 mm is located between the ROI and the transducer. The thickness of the scattering layer is 10 mm. The concentration of the scatterers is 16.67 rods/cm^{2}. The scattering mean free path is 3.13 ± 0.14 mm for this medium. The time window length (*∆t*) is set to 5 μs to promise that the direct waves within ROI are in the same time window. The ultrasonic transducer (WS1-10, ULTRAN) with a central frequency of 10 MHz detects the photoacoustic signal propagating through the scattering layer. The transducer is shifted along the *x*-direction with 0.5 mm step to simulate a 49-element array. The detected signals are amplified (SA-230F5, NF) and recorded at a sampling frequency of 60 MHz (PCI-5105, NI). Figure 5(b)-5(d) illustrate the experimental results. Figure 5(b) presents the measured signals, where the random pattern implies the domination of the scattered waves. The axial resolution of the classical beamforming method in the Fig. 5(c) is 1.37 mm. And the strong acoustic scattering causes the two targets unable to be clearly distinguished in the transverse direction by a traditional beamforming method [Fig. 5(c)]. Then, we use the first two singular values which contain 80.15% of the total energy, to reconstruct the by the proposed method. As shown in Fig. 5(d), the speckles in the background are significantly reduced. Two pencil leads can be clearly identified in the image. The axial and transverse resolution in the Fig. 5(d) is 1.51 mm and 2.13 mm, respectively. The experiments validate the performance of the proposed method.

## 4. Conclusion

In summary, we reveal the different coherent properties of the scattering wave and the direct wave by using a matrix form. The direct waves show a particular coherence on the antidiagonals of the matrix, while the scattering waves do not. Based on this finding, we propose a correlation filter combining with a time reversal operator to achieve photoacoustic image behind a scattering layer. Numerical simulations show that the proposed approach effectively increases the image contrast and decreases the background speckles in a scattering medium. In comparison with the traditional delay-and-sum method, the method can improve the image SNR by about 6 dB. Finally, the experiments demonstrate the proposed method recovers the photoacoustic image successfully in a complex scattering environment, while the transitional delay-and-sum method failed. This study might be valuable in improving the quality of photoacoustic imaging and broadening its applications in scattering tissue.

## Acknowledgement

This work was supported by the National Basic Research Program of China under Grant No. 2016YFC0102300, the NSF of China under Grant Nos. 11422439.

## References and links

**1. **L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science **335**(6075), 1458–1462 (2012). [CrossRef] [PubMed]

**2. **L. V. Wang and J. Yao, “A practical guide to photoacoustic tomography in the life sciences,” Nat. Methods **13**(8), 627–638 (2016). [CrossRef] [PubMed]

**3. **L. Xi and H. Jiang, “High resolution three-dimensional photoacoustic imaging of human finger joints in vivo,” Appl. Phys. Lett. **15**, 18076–18081 (2015).

**4. **M. S. Singh and H. Jiang, “Estimating both direction and magnitude of flow velocity using photoacoustic microscopy,” Appl. Phys. Lett. **104**(25), 1145–1151 (2014). [CrossRef]

**5. **J. Xiao, L. Yao, Y. Sun, E. S. Sobel, J. He, and H. Jiang, “Quantitative two-dimensional photoacoustic tomography of osteoarthritis in the finger joints,” Opt. Express **18**(14), 14359–14365 (2010). [CrossRef] [PubMed]

**6. **A. Dima and V. Ntziachristos, “Non-invasive carotid imaging using optoacoustic tomography,” Opt. Express **20**(22), 25044–25057 (2012). [CrossRef] [PubMed]

**7. **Z. Yang, J. Chen, J. Yao, R. Lin, J. Meng, C. Liu, J. Yang, X. Li, L. Wang, and L. Song, “Multi-parametric quantitative microvascular imaging with optical-resolution photoacoustic microscopy in vivo,” Opt. Express **22**(2), 1500–1511 (2014). [CrossRef] [PubMed]

**8. **B. Ning, M. J. Kennedy, A. J. Dixon, N. Sun, R. Cao, B. T. Soetikno, R. Chen, Q. Zhou, K. Kirk Shung, J. A. Hossack, and S. Hu, “Simultaneous photoacoustic microscopy of microvascular anatomy, oxygen saturation, and blood flow,” Opt. Lett. **40**(6), 910–913 (2015). [CrossRef] [PubMed]

**9. **L. Yao, L. Xi, and H. Jiang, “Photoacoustic computed microscopy,” Sci. Rep. **4**(1), 4960 (2014). [CrossRef] [PubMed]

**10. **J. R. Rajian, M. L. Fabiilli, J. B. Fowlkes, P. L. Carson, and X. Wang, “Drug delivery monitoring by photoacoustic tomography with an ICG encapsulated double emulsion,” Opt. Express **19**(15), 14335–14347 (2011). [CrossRef] [PubMed]

**11. **L. Xi, S. R. Grobmyer, L. Wu, R. Chen, G. Zhou, L. G. Gutwein, J. Sun, W. Liao, Q. Zhou, H. Xie, and H. Jiang, “Evaluation of breast tumor margins in vivo with intraoperative photoacoustic imaging,” Opt. Express **20**(8), 8726–8731 (2012). [CrossRef] [PubMed]

**12. **X. Xiong, Y. Sun, A. Sattiraju, Y. Jung, A. Mintz, S. Hayasaka, and K. C. Li, “Remote spatiotemporally controlled and biologically selective permeabilization of blood-brain barrier,” J. Control. Release **217**, 113–120 (2015). [CrossRef] [PubMed]

**13. **Y. Sun and B. O’Neill, “Imaging high-intensity focused ultrasound-induced tissue denaturation by multispectral photoacoustic method: an ex vivo study,” Appl. Opt. **52**(8), 1764–1770 (2013). [CrossRef] [PubMed]

**14. **J. Xia, G. Li, L. Wang, M. Nasiriavanaki, K. Maslov, J. A. Engelbach, J. R. Garbow, and L. V. Wang, “Wide-field two-dimensional multifocal optical-resolution photoacoustic-computed microscopy,” Opt. Lett. **38**(24), 5236–5239 (2013). [CrossRef] [PubMed]

**15. **D. Wang, Y. Wang, Y. Zhou, J. F. Lovell, and J. Xia, “Coherent-weighted three-dimensional image reconstruction in linear-array-based photoacoustic tomography,” Biomed. Opt. Express **7**(5), 1957–1965 (2016). [CrossRef] [PubMed]

**16. **J.-L. Robert and M. Fink, “Green’s function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging,” J. Acoust. Soc. Am. **123**(2), 866–877 (2008). [CrossRef] [PubMed]

**17. **X. Jin, C. Li, and L. V. Wang, “Effects of acoustic heterogeneities on transcranial brain imaging with microwave-induced thermoacoustic tomography,” Med. Phys. **35**(7Part1), 3205–3214 (2008). [CrossRef]

**18. **M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. **77**(4), 041101 (2006). [CrossRef]

**19. **X. L. Dean-Ben, V. Ntziachristos, and D. Razansky, “Statistical optoacoustic image reconstruction using a-priori knowledge on the location of acoustic distortions,” Appl. Phys. Lett. **98**(17), 412 (2011). [CrossRef]

**20. **D. Wu, C. Tao, and X. Liu, “Photoacoustic tomography in scattering biological tissue by using virtual time reversal mirror,” J. Appl. Phys. **109**(8), 084702 (2011). [CrossRef]

**21. **D. Wu, C. Tao, and X. Liu, “Photoacoustic tomography extracted from speckle noise in acoustically inhomogeneous tissue,” Opt. Express **21**(15), 18061–18067 (2013). [CrossRef] [PubMed]

**22. **C. Yoon, J. Kang, S. Han, Y. Yoo, T. K. Song, and J. H. Chang, “Enhancement of photoacoustic image quality by sound speed correction: ex vivo evaluation,” Opt. Express **20**(3), 3082–3090 (2012). [CrossRef] [PubMed]

**23. **J. Yin, C. Tao, P. Cai, and X. Liu, “Photoacoustic tomography based on the Green’s function retrieval with ultrasound interferometry for sample partially behind an acoustically scattering layer,” Appl. Phys. Lett. **106**(23), 503 (2015). [CrossRef]

**24. **Y. Liu, Y. Shen, C. Ma, J. Shi, and L. V. Wang, “Lock-in camera based heterodyne holography for ultrasound-modulated optical tomography inside dynamic scattering media,” Appl. Phys. Lett. **108**(23), 231106 (2016). [CrossRef] [PubMed]

**25. **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**(1), 58–64 (2013). [CrossRef]

**26. **S. Kang, S. Jeong, W. Choi, H. Ko, T. D. Yang, J. H. Joo, J. S. Lee, Y. S. Lim, Q. H. Park, and W. Choi, “Imaging deep within a scattering medium using collective accumulation of single-scattered waves,” Nat. Photonics **9**, 253–258 (2015).

**27. **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**(10), 100601 (2010). [CrossRef] [PubMed]

**28. **A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, “Smart optical coherence tomography for ultra-deep imaging through highly scattering media,” Sci. Adv. **2**(11), e1600370 (2016). [CrossRef] [PubMed]

**29. **A. Badon, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, “Retrieving time-dependent Green’s functions in optics with low-coherence interferometry,” Phys. Rev. Lett. **114**(2), 023901 (2015). [CrossRef] [PubMed]

**30. **A. Aubry and A. Derode, “Random matrix theory applied to acoustic backscattering and imaging in complex media,” Phys. Rev. Lett. **102**(8), 084301 (2009). [CrossRef] [PubMed]

**31. **S. Shahjahan, A. Aubry, F. Rupin, B. Chassignole, and A. Derode, “A random matrix approach to detect defects in a strongly scattering polycrystal: how the memory effect can help overcome multiple scattering,” Appl. Phys. Lett. **104**(23), 046607 (2014). [CrossRef]

**32. **A. Aubry, L. A. Cobus, S. E. Skipetrov, B. A. van Tiggelen, A. Derode, and J. H. Page, “Recurrent scattering and memory effect at the Anderson localization transition,” Phys. Rev. Lett. **112**(4), 043903 (2014). [CrossRef] [PubMed]

**33. **A. Derode, A. Tourin, J. de Rosny, M. Tanter, S. Yon, and M. Fink, “Taking advantage of multiple scattering to communicate with time-reversal antennas,” Phys. Rev. Lett. **90**(1), 014301 (2003). [CrossRef] [PubMed]

**34. **A. M. Tulino, S. Verd, *Random matrix theory and wireless communications* (Now Publishers Inc., 2004).

**35. **R. Sprik, A. Tourin, J. de Rosny, and M. Fink, “Eigenvalue distributions of correlated multichannel transfer matrices in strongly scattering systems,” Phys. Rev. B **78**(1), 012202 (2008). [CrossRef]

**36. **A. Aubry and A. Derode, “Detection and imaging in a random medium: A matrix method to overcome multiple scattering and aberration,” J. Appl. Phys. **106**(4), 044903 (2009). [CrossRef]

**37. **W. V. Meyer, D. S. Cannell, A. E. Smart, T. W. Taylor, and P. Tin, “Multiple-scattering suppression by cross correlation,” Appl. Opt. **36**(30), 7551–7558 (1997). [CrossRef] [PubMed]

**38. **C. Prada and J. L. Thomas, “Experimental subwavelength localization of scatterers by decomposition of the time reversal operator interpreted as a covariance matrix,” J. Acoust. Soc. Am. **114**(1), 235–243 (2003). [CrossRef] [PubMed]

**39. **A. Derode, V. Mamou, and A. Tourin, “Influence of correlations between scatterers on the attenuation of the coherent wave in a random medium,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **74**(3), 036606 (2006). [CrossRef] [PubMed]

**40. **A. Aubry and A. Derode, “Singular value distribution of the propagation matrix in random scattering media,” Waves Random Complex Media **20**(3), 333–363 (2010). [CrossRef]

**41. **D. Wu, C. Tao, X. J. Liu, and X. D. Wang, “Influence of limited-view scanning on depth imaging of photoacoustic tomography,” Chin. Phys. B **21**(1), 014301 (2012). [CrossRef]