High-throughput volumetric adaptive optical imaging using compressed time-reversal matrix

Deep-tissue optical imaging suffers from the reduction of resolving power due to tissue-induced optical aberrations and multiple scattering noise. Reflection matrix approaches recording the maps of backscattered waves for all the possible orthogonal input channels have provided formidable solutions for removing severe aberrations and recovering the ideal diffraction-limited spatial resolution without relying on fluorescence labeling and guide stars. However, measuring the full input–output response of the tissue specimen is time-consuming, making the real-time image acquisition difficult. Here, we present the use of a time-reversal matrix, instead of the reflection matrix, for fast high-resolution volumetric imaging of a mouse brain. The time-reversal matrix reduces two-way problem to one-way problem, which effectively relieves the requirement for the coverage of input channels. Using a newly developed aberration correction algorithm designed for the time-reversal matrix, we demonstrated the correction of complex aberrations using as small as 2% of the complete basis while maintaining the image reconstruction fidelity comparable to the fully sampled reflection matrix. Due to nearly 100-fold reduction in the matrix recording time, we could achieve real-time aberration-correction imaging for a field of view of 40 × 40 µm2 (176 × 176 pixels) at a frame rate of 80 Hz. Furthermore, we demonstrated high-throughput volumetric adaptive optical imaging of a mouse brain by recording a volume of 128 × 128 × 125 µm3 (568 × 568 × 125 voxels) in 3.58 s, correcting tissue aberrations at each and every 1 µm depth section, and visualizing myelinated axons with a lateral resolution of 0.45 µm and an axial resolution of 2 µm.


Introduction
An arbitrary optical system interacting with light waves can be described by transmission and reflection matrices, as far as the linear light-matter interaction is concerned. The transmission (reflection) matrix of an optical system describes the transmitted (reflected) electric field (E-field) at all the possible detection channels for a set of orthogonal illumination channels. Due to the characterization of the input-output response, the measured matrix can be considered as a replica of a real optical system within the context of the covered illumination/detection channels. Therefore, one can computationally process it as though a real experiment is being conducted. The knowledge of the matrix allows one to find solutions in a variety of applications where lengthy experimental optimizations are required. Examples include focusing light 1,2 , delivering images 3 , and controlling transmission power 4 through scattering media based on the transmission matrix. The reflection matrix, suitable for more realistic in vivo applications for which the detector cannot be placed on the transmission side, has provided exceptional opportunities for deep-tissue imaging 5 . The reflection matrix has also been exploited to focus light on a target embedded deep within strongly scattering media [6][7][8] . The wave correlation of the singlescattered waves in the reflection matrix was tailored to attenuate the effect of multiple light scattering 9 . A wavefront correction algorithm termed closed-loop accumulation of single scattering (CLASS) 10 was developed based on the time-gated reflection matrix for separately identifying the aberrations in the illumination and detection pathways without the need for guide stars and in the presence of strong multiple scattering noise. This offers imaging deep within biological tissues with a subdiffraction-limited resolution 11 . The singular value decomposition (SVD) was applied to a time-gated reflection matrix for retrieving a target image underneath strongly scattering media 12,13 . Recently, it has been demonstrated that the time-gated reflection matrix measured in the space domain made it possible to image a mouse brain through an intact skull inducing extreme aberrations 14 . Indeed, the reflection matrix approaches provide formidable solutions in the context of computational adaptive optics (AO) microscopy 15,16 in that they can deal with extremely severe aberrations with no need for fluorescence labeling and guide stars. In addition, this space-domain reflection matrix study proved that it can serve as a type of wavefront sensorless AO 17,18 that is combined with hardware correction of aberration by wavefront shaping devices such as a spatial light modulator and deformable mirror to realize ideal diffraction-limited multi-photon fluorescence imaging through an intact skull 14 .
Despite these benefits, the matrix-based AO approach has been elusive in real-time bio-medical imaging applications. The recording of the full reflection matrix is a time-consuming process because the E-field map of the reflected wave must be measured for each illumination channel, as opposed to confocal imaging's requirement of point detection. Furthermore, the interferometric detection of the E-field is sensitive to the random phase drift, which can deteriorate the recorded reflection matrix in the dynamic samples. Sparse sampling of the matrix would be a potential solution, but this is accompanied by incomplete sampling of illumination channels. Considering that finding object information embedded within a scattering medium requires identifications of wave distortions in both the illumination and detection pathways, insufficient sampling of the illumination channels can undermine the capability to resolve illumination distortions.
To overcome these issues, we consider a time-reversal matrix RR y instead of the reflection matrix R. Here R y represents the conjugate transpose of R. Unlike the reflection matrix itself, which describes the relationship between the illumination and detection channels, the time-reversal matrix describes the phase-conjugated roundtrip process connecting the detection channels to the same detection channels via the illumination channels. Essentially, this reduces the two-way problem with the reflection matrix to the one-way problem with the time-reversal matrix on condition that the illumination channels are orthogonal. There are two major benefits of dealing with the time-reversal matrix. It can maintain high fidelity in terms of retaining the information on the detection channels even if the illumination channel coverage is much smaller than that of the complete set. Furthermore, it is not even necessary to know the basis of the illumination channels, making it robust to the random phase drift.
Here, we present a high-throughput volumetric AO imaging method termed a compressed time-reversal closed-loop accumulation of single scattering (CTR-CLASS), in which the previously developed CLASS algorithm was extended to a compressed time-reversal matrix constructed by a sparsely sampled reflection matrix for correcting the complex sample-induced aberrations with significantly reduced number of measurements. In this implementation, we took advantage of the time-reversal matrix and made use of dynamically varying unknown speckle patterns as illumination channels. We demonstrated that both the aberration map and object image can be retrieved using the number of speckle patterns as small as 2% of the complete basis while maintaining comparable fidelity to that of the fully sampled matrix. Due to nearly 100-fold reduction of the matrix recording time, the CTR-CLASS has enabled real-time aberration-correction imaging for a field of view (FOV) of 40 × 40 µm 2 (176 × 176 pixels) at a frame rate of 80 Hz. We applied the developed method for the volumetric AO imaging of ex vivo mouse brain and visualized myelinated axons with a lateral resolution of 0.45 µm and axial resolution of 2 µm over a volume of 128 × 128 × 125 µm 3 (568 × 568 × 125 voxels) within a recording time of 3.58 s.

Reflection matrix description of an imaging system
Let us first start with a mathematical model for an optical imaging system of interest using reflection matrix formalism. We consider the time-gated coherent imaging of a target object through a scattering medium in reflection geometry (Fig. 1a). For convenience, the optical layout is unfolded by flipping the reflection beam path over an object plane, making the layout analogous to transmission geometry. Since the scattering sample serves as a linear system with respect to the E-field in coherent imaging, the reflected wave can be described by a linear superposition of impulse response functions, Here, E o (r o ;r i ) is the time-gated E-field at position r o on the output (detection) plane when a target object is illuminated by a point source located at position r i on the input (illumination) plane. O(r) is the object function that represents complex reflection coefficients of the target object. Both the input and output planes are conjugate to the object plane whose spatial coordinate is r. P i(o) (r i(o) ; r) is the time-gated E-field point-spread-function (PSF) that represents E-field distribution at the input (output) plane generated by a point source located at a position r on the object plane. The PSF describes transmissions of ballistic waves that maintain their propagation directions in propagation through the scattering medium. Therefore, it has the linear shift-invariant property 19 , i.e., P iðoÞ ðr iðoÞ ; rÞ ¼ P iðoÞ ðr iðoÞ À rÞ. E ms (r o ; r i ) represents speckle noise generated by multiple-scattered waves in the scattering medium, whose flight times fall within a finite time-gating window. In scattering matrix formalism, Eq. (1) can be represented by a time-gated reflection matrix R whose element is E o (r o ; r i ) for a column index r i and row index r o . Based on Eq. (1), R can be decomposed as Here, O is a diagonal matrix whose diagonal element is O(r). P i and P o are Toeplitz (diagonal-constant) matrices whose elements are respectively given by the input and output PSFs, P i (r i ; r) and P o (r o ; r). R ms is a matrix composed of E ms (r o ; r i ), and the superscript 'T' denotes matrix transpose operation. The first term on the right-hand side in Eq. (2) is responsible for image reconstruction. Based on the perturbation of the firstorder Born approximation, the first term assumes that the incident wave does not change its propagation direction until it reaches the object plane, but does experiences phase retardation by the scattering For clarity, the reflection pathway is unfolded to the transmission side. O is the object's reflection coefficient matrix, and P iðoÞ ¼ ½P iðoÞ ðr; r iðoÞ Þ is the transmission matrix of the scattering medium for wave propagation from position r to r iðoÞ . The red and blue curves represent incident and reflected waves, respectively. b Reflection matrix R S in the case of speckle illumination. R S ½r o ; j ¼ ½E o ðr o ; jÞ is a set of E-field responses E o ðr o ; jÞ for speckle input channels Sðr i ; jÞ. R S can be written as R S = RS, where S½r i ; j ¼ ½Sðr i ; jÞ is the input illumination matrix. c Geometric interpretation of the time reversal matrix. The time-reversal matrix W½r o ; r ′ o ¼ R S R y S can be interpreted as a reflection matrix describing a roundtrip process for a wave propagating from the output plane r o to the object plane r and reflecting back to the output plane by the target object |O| 2 medium. For identical illumination and detection paths, input and output PSFs are the same due to the reciprocity principle in optics and thus satisfy the relation P i ðr 0 ; rÞ ¼ P o ðr 0 ; rÞ. However, this is not a necessary condition in the present study. As we reported earlier 10 , the CLASS algorithm utilizing R is developed in such a way to find high-fidelity solutions for an unknown object function (O) and two PSFs (P i and P o ) even in the presence of strong multiple scattering noise R ms .

Compressive sensing of the reflection matrixS
The matrix R in Eq. (2) can be directly measured by scanning the position r i of the focused illumination and wide-field detection of the backscattered wave field E o ðr o ; r i Þ across r o . To obtain the full time-gated reflection matrix for a given FOV, it is necessary to scan a focused beam over the FOV with a lateral sampling interval of the diffraction-limited resolution, Δx = λ/(2NA), where λ is the wavelength of light source, and NA is the objective numerical aperture. For a 2-dimensional (2D) FOV of size L × L, the required number of sampling points for a complete sampling is given by N = L/Δx which is the total number of orthogonal spatial modes for the given FOV, NA, and λ. In fact, the full reflection matrix R 2 C N N can be measured by sending any complete set fE i ðr i ; jÞg N j¼1 of N illumination fields, instead of point-by-point scanning with a focused beam. One can measure the respective output field E o ðr o ; jÞ for each j th illumination and construct a reflection matrix R m 2 C N N whose columns are assigned by the measured output E-fields fE o ðr o ; jÞg N j¼1 . In this case, the column and row indices of R m are j and r o , respectively. Then, the measured R m is expressed as R m ¼ RE i , where E i 2 C N N is an illumination matrix constructed by fE i ðr i ; jÞg N j¼1 in the same way as R m . The reflection matrix of the sample can be obtained by multiplying the measured matrix R m by the i . This requires knowledge of the illumination fields.
Complete sampling of the reflection matrix for a large FOV can be time-consuming and resource-intensive, and thus limits its practical applicability. Main objective of this study is to reduce the data acquisition time. Here, we propose the use of a set of M unknown random speckle illumination patterns, fSðr i ; jÞg M j¼1 for compressive sensing of R. In particular, we consider the case in which M is significantly smaller than N (M ≪ N). The time-gated E-field image E o ðr o ; jÞ is recorded for each j th speckle illumination, and the sparsely-sampled reflection matrix R S 2 C N M is then constructed using fE o ðr o ; jÞg M j¼1 as a matrix element (Fig. 1b). Therefore, the matrix R S with column index j and row index r o can be expressed as where S 2 C N M is the sensing matrix constructed by fSðr i ; jÞg M j¼1 . To realize aberration correction and image reconstruction without a prior knowledge of the illumination pattern, we consider a CTR matrix, where W ms is the noise matrix associated with the multiple scattering R ms , and the superscript '*' denotes the complex conjugate. Note that the covariance matrix SS † ∈ C N×N is almost an identity matrix I for sufficiently large M. However, finite overlaps among random speckles can cause non-zero off-diagonal elements, i.e., SS † = I + σ, where σ denotes the additive complex random noise whose matrix elements σ½r i ; r 0 i are given by correlations of two series of random speckle fields illuminating different positions, r i and r 0 i . Therefore, W can be expressed as where O I denotes |O| 2 , a diagonal matrix with its diagonal elements given by the reflectance of the object, |O(r)| 2 . The first term in Eq. (4) uses the relation P T i P Ã i ¼ P i P y i ¼ I, which is valid when P i ðr i ; rÞ is a shift-invariant PSF induced by a phase-only pupil aberration. The last term W σ ¼ P o OP T i σP Ã i O y P y o denotes the sparse sampling-induced noise associated with σ. Statistically, mean amplitude of σ½r i ; r 0 i is given by hj where the bracket notation denotes an average over the variables in the subscript. Sðr i ; jÞ is normalized such that h P M j¼1 jSðr i ; jÞj 2 i r i ¼ 1. Therefore, the magnitude of the matrix elements of W σ becomes smaller as M increases, which makes W σ negligible for sufficiently large M.
Physical interpretation of the time-reversal matrix W is given in Fig. 1c. In this discussion, we excluded the noise terms in Eq. (4) to focus our attention more on a successful time-reversal process. The noise degrades the fidelity of the aberration correction and the signal-to-noise ratio (SNR) of the reconstructed image, which will be discussed in detail in "Analysis of image SNR" section. By the successive timereversal operation R y S , a spherical wave (red curves) emitted from a point source at a position r 0 o on the output plane propagates in the backward direction through the object (P Ã i O y P y o ) followed by a fictious scattering layer whose transmission matrix is S † . Afterwards, the R S is applied such that the reflected wave (blue curves) returns to the scattering layer (S) and the object (P o OP T i ) in the forward direction to arrive at the output plane. Here, the important point is that the operation indicated by the shaded gray box (P T i SS y P Ã i ) serves as a phase-conjugation mirror when SS y % I, i.e., the illumination speckles are sufficiently orthogonal. In other words, a point source emanating from an object plane comes back to the same position via its travel through P T i SS y P Ã i . This eliminates the need to consider the input aberration P i and illumination patterns S. As a result, the matrix W ¼ R S R y S can be interpreted as a time-gated reflection matrix describing an imaging system that images a reflective object with the reflectance O I through a scattering medium with an input transmission matrix of P Ã o and an output transmission matrix of P o . The whole process becomes an N-by-N square matrix with its column and row indices both corresponding to r o .
There are a few major benefits of considering W instead of R. First, the sensing matrix S describing the illumination patterns is removed in W, thereby eliminating the need to know the illumination speckle patterns. This makes it possible to send an arbitrary choice of illuminations, such as dynamic speckle patterns generated by a rotating diffuser, and thus it is no longer necessary to scan pre-defined positions of point illumination using scanning mirrors. Second, W is greatly simplified such that only P o and O I remain to be identified. Imperfection in illumination and detection optics often causes discrepancy between P i and P o in the reflection geometry. P i and P o are intrinsically different in the transmission geometry. However, in W, it is not necessary to consider P i . Another critical benefit is the possibility of sparse sampling. The condition SS † ≈ I satisfies even when M is extremely small. In contrast, if there is significant downsampling in the focused illumination, both the ability to identify the aberration and the imaging fidelity are significantly degraded.

Image reconstruction
The concept of the time-reversal matrix was initially introduced for selective focusing in acoustics [20][21][22] and then has been intensively studied in microwaves 23,24 and optics 7,12,25 . In these previous studies, either iterative operation or the singular value decomposition (SVD) of the time-reversal matrix W = R † R was used for selective focusing on a few highly reflecting targets embedded in a scattering medium. Each input singular vector with a nonzero singular value of W corresponds to a specific wavefront of the incident light focusing on one of the targets whose reflectivity is proportional to the eigenvalue. This approach has also been applied to deep optical imaging of an extended target in a scattering medium, where the image is reconstructed by using the dominant singular values and corresponding singular vectors 12,13 .
In contrast, we consider the CTR matrix W ¼ R S R y S and introduce a matrix decomposition of W into a product of three matrices, W ¼ P o O I P y o in order to find the unknown object function (O I ) embedded within a scattering medium inducing optical aberrations (P o ). In conventional compressive sensing 26,27 , the sampling process is modeled as y = Sa + n, where y 2 R M 1 is a vector of sparsely sampled data, a 2 R N 1 is the original signal to be recovered, S 2 R M N is a known sensing (or measurement) matrix with M ≪ N, and n denotes a measurement noise. Here, M is the number of measurements and N is the number of signals of interest. The problem of recovering the signal a is underdetermined or ill-posed because there are more unknowns than equations. To reliably solve the problem, we need a prerequisite that the degree of signal sparsity S (the number of nonzero elements in a) is smaller than the number of measurements, M: M ! OðS logðN=SÞÞ.
Conventional compressive sensing concerns the sparsity of the target a. In our study, the reflection matrix R 2 C N N is sparsely sampled by the sensing matrix S 2 C N M : R S = RS. Thus, the degree of sparsity of R is of concern. Essentially, our model R ¼ P o O I P T i allows us to treat R as a highly sparse matrix. Since O I is diagonal and P oðiÞ is a Toeplitz matrix, R contains only 3N unknowns. Since the measured R S contains the number of elements N × M, it forms a system of N × M equations with 3N unknowns. Theoretically, there needs M ! Oð3 logðN=3ÞÞ to accurately estimate the solutions. However, we used the unknown speckle illuminations, i.e. unknown sensing matrix S, and thus converted the problem to a decomposition of time reversal matrix, where W ms is due to multiple scattering noise, and W σ is the additional sparsity-induced noise. Since W σ scales with 1= ffiffiffiffiffi M p , the reduction of sampling M gives rise to the increase of noise. This makes it necessary to set M larger than the estimated value to properly decompose W.
To identify the aberration and reconstruct an aberration-corrected image, the basis of W is changed to the spatial-frequency domain (k-space) by taking the Fourier transform for both the column and row bases. The CTR matrix in k-space is expressed as where a tilde (~) above the variable denotes the Fourier transformation of the variable into spatial-frequency domain. The matrixP o ½k 0 o ; k o represents the planewave based transmission matrix between the object and output planes. When the output PSF is space-invariant,P o becomes a diagonal matrix whose elements are given by its complex pupil functionP o ðk o Þ ¼ e iϕ o ðk o Þ , where ϕ o ðk o Þ is the output pupil phase map. We consider here a phase-only pupil function that has amplitude of unity. The matrixÕ I is the target spectrum matrix in which each column consists of a shifted special-frequency spectrum of the target: 10 , acting on the CTR matrixW , identifies an aberration correction matrixP c that maximizes the total intensity of the object imagẽ The final output correction phase function is found by accumulating all the preceding correction phases,

Experimental setup of CTR-CLASS microscopy
The schematic of the experimental setup is shown in Fig. 2a for recording a CTR matrix W. The basic configuration is a low-coherence wide-field interferometric microscope, but a sample is illuminated by random speckle fields while the reference wave is a planar wave. A  custom-built wavelength-tunable Ti:Sapphire pulsed laser (center wavelength of 800-900 nm, bandwidth of 30 nm) was used as a low-coherence light source. An optical diffuser mounted on a motorized rotation stage was inserted at a conjugate image plane in the illumination path to produce uncorrelated random speckle fields for the sample wave. Backscattered sample waves from the target were captured by an objective lens (Nikon, ×60, NA 1.0) and delivered to a high-speed CMOS camera (Photron, FASTCAM mini UX100) placed at a conjugate image plane. A reference plane wave was introduced at the camera to generate the off-axis low-coherence interferogram, from which we obtained the time-gated E-field of the backscattered sample wave (see Supplementary Information Note 1 for the detailed setup). Figure 2b shows three representative E-field amplitude and phase images of the sample under different speckled illuminations. To obtain a complete N-by-N reflection matrix for a FOV having N orthogonal modes, N E-field images must be acquired for an orthogonal set of N illumination fields, where each image has a total of N orthogonal pixels. However, we recorded only M (<N) E-field images using random speckle illuminations to reduce the acquisition time and reconstruct a sparsely sampled N-by-M reflection matrix R S as shown in Fig. 2c. For efficient sparse sampling of the reflection matrix, it is important to minimize the correlation between the speckle fields as much as possible. For a given camera exposure time and frame rate, the angular velocity of the rotating diffuser was carefully selected to minimize the spatial correlation between two consecutive speckle patterns ( Supplementary Information Note 2). To justify the validity of SS † ≈ I, we experimentally measured the S matrix by placing a mirror at the sample plane. Figure 2d shows SS † matrix obtained using M = 700 speckled illuminations for a FOV of 40 × 40 µm 2 with a diffractionlimited resolution of 450 nm, resulting in the total number of orthogonal modes, N = 88 × 88 = 7744. The SS † was nearly diagonal matrix with the ratio between offdiagonal and diagonal elements of ≲0.1.

Proof-of-concept experiment
To demonstrate the high-throughput data acquisition and aberration correction capabilities of the CTR-CLASS microscopy, we imaged a homemade Siemens star target covered by a 600-µm-thick plastic layer introducing strong optical aberrations. The laser operated at a center wavelength of 900 nm, and had a coherence length of 12 µm. For a FOV of 40 × 40 µm 2 (N = 88 × 88 pixels), M speckled E-field images of the target were imaged by the high-speed camera operating at a frame rate of 12,500 Hz with an exposure time of 20 µs.
We define the compression ratio that indicates the degree of sparse sampling as CR = M/N. Figure 3a shows the CTR matrix W constructed by speckle patterns for CR = 1. In the absence of aberration (P o = I), the matrix W % P o O I P y o is almost diagonal because it is reduced to O I . In the presence of aberration and scattering, the signal in the diagonal spreads out to the off-diagonal elements. The matrixW in spatial-frequency domain and the identified correction phase ϕ c ðk o Þ by the CTR-CLASS algorithm are shown in Fig. 3b, c, respectively. The corrected CTR matrix W c ¼ P y c W P c is shown in Fig. 3d. The intensity images before and after the aberration correction were reconstructed fromW andW c , respectively. The uncorrected image shown in Fig. 3e is blurry and hardly recognizable, while object structures are sharply resolved in the corrected image in Fig. 3f. Both the imaging resolution and the signal-to-background ratio (SBR) were significantly improved compared to those in the uncorrected image. We quantified the degree of aberration correction by measuring the normalized intensity profiles of PSFs before and after the aberration correction ( Fig.  3g-i). The Strehl ratio, defined by the peak intensity of the PSF, is a measure of AO performance. The Strehl ratio α o before aberration correction is given by α o ¼ jP o ð0Þj 2 . Likewise, the Strehl ratio α c after correction is given by α c ¼ jP res ð0Þj 2 , where P res (r) is the residual PSF after the correction (P res ¼ P y c P o ). The enhancement in the Strehl ratio, (α c /α o ) measured from the line profiles of the PSFs was at least 20. The measured full-width at half-maximum (FWHM) of the aberration-corrected PSF was about 450 nm, which is the diffraction-limit spatial resolution of the system.

Analysis of image SNR
We evaluated the performance of image recovery depending on CR to determine the minimum achievable CR. Reconstructed images and aberration phase maps for various CR are shown in Fig. 4a. The first row shows the reconstructed intensity images normalized by M for CR = 0.5, 0.1, 0.02, and 0.017, and the second row shows the corresponding aberration phase maps. The identified aberration maps were almost identical all the way into the high spatial frequency range, although CR was significantly reduced. Diffraction-limited high-resolution images could be successfully restored for CR ≥ 0.02 (M ≥ 155). Considering the camera frame rate of 12,500 Hz, it took only 12.4 ms to record the CTR matrix W for CR = 0.02, setting the highest achievable aberration-correction image frame rate to 80 Hz for a FOV of 40 × 40 µm 2 (88 × 88 pixels). If the CR was further reduced to less than 0.02, the image reconstruction failed to find the correct aberration map and object image. The line profiles along the white dotted lines on the reconstructed images in Fig. 4a are compared in Fig. 4b to quantify the image quality. Interestingly, neither the image contrast nor the spatial resolution of the reconstructed image was diminished by the reduction of CR, so far as the image reconstruction was successful. This means that even the severe decrease in CR does not hamper the performance of aberration correction. We compared the residual root-mean-square (RMS) wavefront errors of the identified wavefront maps relative to the aberration map obtained from the full reflection matrix (blue dots in Fig. 4c). The residual RMS wavefront errors remained nearly constant regardless of CR. This indicates that aberrations were properly corrected even with the small M. The pupil phase retardation ϕ ðnÞ c ðk o Þ is estimated from the inner product of the angular spectrum image for the illumination angle of k o , W ½k o þ Δk; k o and the angular synthesized imagẽ O I ðΔkÞ in Eq. (6), where the summation is taken over total N spectral frequencies, Δk. Thus, the standard error of the estimated pupil phase scales with 1= ffiffiffiffi N p . Essentially, the performance of aberration correction is determined by N, the number of sampled spectral frequencies or equivalently the image size. Therefore, as long as N is sufficiently large, the pupil phase retardation can be identified with high fidelity even for a small number of measurements, M.
An important figure of merit in imaging is the signal-tonoise ratio (SNR), which is defined by the ratio of the mean intensity of the target image and the standard deviation of random background noise. To estimate the image SNR, we consider the aberration-corrected CTR matrix W c ¼ P y c W P c ¼ P y c P o O I P y o P c þ P y c W ms P c þ P y c W σ P c . Note that the incident speckle field is normalized as h P M j¼1 jSðr i ; jÞj 2 i r i ¼ 1, which means the average intensity of speckle illumination fields at each pixel r i is 1/M. For intuitive understanding, let us consider sending speckle illumination fields whose average intensity per pixel is 1.
We define the signal s c as the intensity of the target reconstructed after the correction, which corresponds to the main diagonal of the first term P y c P o O I P y o P c . Therefore, s c is given by s c ¼ Mα c O I ðrÞ. There are mainly two kinds of noise sources: (1) measurement noise W ms and (2) sparse sampling-induced noise W σ . First, the measurement noise includes multiple scattering noise, photon shot noise, and dark count noise of a camera sensor. In the case of imaging through a scattering medium, the measurement noise is dominated by multiple scattering noise, which is much larger than photon shot noise and dark count noise. Due to the random nature of noise, the standard deviation of the measurement noise, σ ms remains almost unchanged after aberration correction and is proportional to the square root Here m denotes the measurement noise in a single wide-field image, which is mainly given by the average intensity of multiple scattering noise. Second, as we discussed earlier, the sparse sampling-induced noise is inherently caused by the overlaps among speckle illumination patterns and is given by The standard deviation of the sparse sampling-induced noise, σ s also scales proportionally with the square root of M: σ s ¼ ffiffiffiffiffi M p hO I ðrÞi r , where hO I ðrÞi r is the average reflectance of the target. In fact, the σ s is analogous to the measurement noise, except that it is proportional to the target's average reflectance. Then, the SNR c can be estimated to be SNR c ¼ s c = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi σ 2 ms þ σ 2 s p . In the strong multiple scattering regime (O I ðrÞ ( m), the SNR c of the reconstructed image from the CTR matrix is approximated as of a single wide-field image obtained without correction.
In the weak scattering regime of SNR 1 ≫ 1, the SNR c is mainly determined by the sparse sampling-induced noise and given by SNR c % s c =σ s ¼ ffiffiffiffiffi M p α c (O I ðrÞ=hO I ðrÞi r ). In this proof-of-concept experiment, the plastic aberrating layer caused strong aberrations, but not much the multiple scattering noise. The Siemens star target used in this demonstration consists of binary patterns with nearly 100% reflectivity (O I ðrÞ ¼ 1) and a fill factor value of 50% (hO I ðrÞi r ¼ 1=2), giving the ratio O I ðrÞ=hO I ðrÞi r % 2. In this specific case of a binary target with 50% fill factor, the SNR c is expressed as SNR ¼ 2 ffiffiffiffiffi M p α c . In terms of compression ratio, it is expressed as SNR . We plotted the SNR as a function of CR in Fig. 4c (red dots) and observed that the SNR fits well with ffiffiffiffiffiffi ffi CR p . The ratio between the diagonal and off-diagonal elements of W c is the signal to background noise ratio, which in turn is equal to the SNR of the reconstructed Uncorrected CR = 0.02 image. We found that ratios between the diagonal to offdiagonal elements of both SS † and W c also fit well with ffiffiffiffiffi M p as shown in Fig. 4d, confirming the validity of our analysis. We found that this ratio fits well with ffiffiffiffiffi M p , confirming the validity of our SNR analysis.
Image reconstruction was successful up to a CR value of 0.02 (M ≥ 155), but it was shown to fail with a further decrease of CR. As explained above, the reduction in CR introduces random noise. Below a certain threshold CR value, the noise level becomes too high for the algorithm to retrieve the pupil phase map. It is difficult to generalize the minimum achievable CR value because it depends on various factors such as the SNR 1 of the uncorrected image, the distribution of the target spectrum, and the complexity of the aberration. Therefore, appropriate choice of CR value is necessary depending on the sample for an optimal image acquisition speed and quality.

Volumetric aberration-correction imaging of a mouse brain
With the CTR-CLASS microscope, we demonstrated the high-throughput aberration-corrected volumetric imaging of myelinated axons in an ex vivo mouse brain. In this demonstration, the pulsed laser operated at a center wavelength of 848 nm with a coherence length of 8 µm. Typically, a set of E-field images constituting a reflection matrix are recorded for each fixed depth to maintain the input and output planes throughout the measurements. This depthwise matrix recording slows down the volumetric imaging. Here, we employed continuous depth scanning to speed up the volume coverage. Since the output planes are continuously varying, we added numerical propagation steps to synchronize the output planes prior to the application of CTR-CLASS algorithm. The synchronization of the input planes is not necessary as there is no need to know the used input channels. To experimentally implement this concept, a whole mouse brain was mounted on a motorized stage and continuously scanned along the z-axis at a constant speed of 35 µm/s while dynamically varying speckle patterns illuminated the specimen (Fig. 5a). The E-field images of the backscattered waves were recorded by the high-speed camera at a frame rate of 5000 Hz. Therefore, there was 7 nm depth difference between the neighboring frames. The camera exposure time (50 µs) and the angular velocity of the rotating diffuser (210 degree/s) were carefully selected to ensure that the spatial correlation of the speckle patterns between the consecutive frames was less than 0.1 (See Supplementary Information Note 2). For the total image acquisition time of 3.58 s, we obtained a series of 17,891 E-field images with a frame size of 128 × 128 µm 2 (300 × 300 pixels), spanning a depth range of 125 µm from the surface of the brain. The angular spectra of the obtained E-field images are filtered by applying a binary pupil mask with NA = 0.94. Then, we finally obtained E-field images with a size of 128 × 128 µm 2 (284 × 284 pixels) and diffraction-limited resolution of 0.45 nm.
To reconstruct a 3D volume image, we first prepared depth-corrected E-field images at 125 depths spaced 1 µm apart. To retrieve an aberration-corrected image at each given depth z, E-field images taken within a range of z ± 4 µm numbering 1147 were numerically propagated to the depth z by adding appropriate quadratic spectral phase factors in their angular spectra. To deal with positiondependent aberrations, the depth-corrected E-field images were divided into 11 × 11 subregions, and the CTR matrix for each subregion was separately constructed using these images. Finally, we retrieved aberration-corrected 2D images by applying the CLASS algorithm to the CTR matrices at individual subregions in all depth z. The reconstructed 3D volume image over 128 × 128 × 125 µm 3 (568 × 568 × 125 voxels) is shown in Fig. 5b (See Supplementary Video 1). Note that the depth-dependent defocus due to refractive index mismatch between immersion water and the tissue causes the separation of the objective focus and coherence volume resulting in blurred images 28 . We could find and compensate the depthdependent defocus by numerically propagating the E-field images such that the total intensity of reconstructed images without aberration correction was maximized. The depth-dependent defocus was about 3 µm at z = 100 µm, which was less than the coherence length (8 µm) of light source. Representative section images at various depths are shown in Fig. 5c. The left-hand column shows maximum intensity projection (MIP) of the reconstructed images without aberration correction, whereas the middle column shows corresponding MIP images after aberration correction. The right-hand column shows identified output pupil phase maps for 11 × 11 subregions. There were no significant aberrations at z = 15 µm except for spherical aberration. Up to a depth of 50 µm, the aberration-corrected images were almost identical to those without correction. As the imaging depth was increased, the aberration maps became more complex, and myelinated fibers in the uncorrected images became blurred due to the inhomogeneity within the tissue. Specifically, there was a blood vessel with a diameter of~30 µm located close to the surface of the brain. The white dashed curves in the uncorrected image at z = 15 µm in Fig. 5c indicate wall boundaries of the blood vessel. We observed that the blood vessel induced pronounced aberrations such that the aberration maps in areas under the vessel were more complex than those in other areas. In addition, correlation between aberration maps of neighboring subregions decreased rapidly, suggesting that the isoplanatic patch size was reduced with depth. At the depth of 115 µm, myelinated axons were almost invisible without aberration correction. Intensity line profiles along the white dotted lines in the images at z = 115 µm are shown in Fig. 5d. Comparing the line profiles, we observed that the CTR-CLASS can recover a nearly diffraction-limit resolution of~0.45 µm (the minimum thickness of myelinated fiber in FWHM) and high-contrast images of myelinated fibers (up to~7-fold increase in signal-to-background ratio). The axial resolution measured from cross-sections of myelin fibers along z-axis was~2 µm.

Discussion
The reflection matrix containing full optical input-output response of a scattering medium has offered robust image reconstruction in comparison with conventional adaptive optics approaches relying on partial information. In particular, it enables the correction of extremely complex aberrations in stringent conditions where there are strong multiple scattering noise and no guide stars available. As a trade-off, the matrix recording is too time-consuming to perform real-time imaging.
Throughout our study, we demonstrated that the use of a time-reversal matrix, instead of the reflection matrix, can be a solution for the high-throughput volumetric imaging equipped with all the benefits of the reflection matrix approaches. We proved that the time-reversal matrix approach can maintain the fidelity of aberration correction and image reconstruction using as small as 2% of the full basis sampling. Due to nearly 100-fold reduction of the matrix recording time, we could achieve aberrationcorrection imaging for a 2D FOV of 40 × 40 µm 2 at a frame rate of 80 Hz. Furthermore, we realized the volumetric imaging of a mouse brain over a volume of 128 × 128 × 125 µm 3 in 3.58 s with a lateral resolution of 0.45 µm and an axial resolution of 2 µm throughout all the voxels including the areas underneath a blood vessel. The proposed method presents a noteworthy conceptual advance. It is a new discovery that the timereversal matrix can be highly compressed in terms of illumination channel coverage. We found that it is not even necessary to know what the illumination channels were. These conceptual findings naturally led to the advances in practicality. In addition to the reduction of illumination channel coverage, there is no need to perform time-consuming pre-calibration to gain prior knowledge on illumination field. It is no longer necessary to concern the phase stability among the E-field images. This enabled us to use dynamically varying random speckle patterns for illumination, instead of laser beam scanning by carefully aligned scanning mirrors, which greatly simplifies the experimental setup. We also presented novel volumetric image processing algorithm that replaces previous depth-wise angular scanning with continuous depth scanning in conjunction with dynamic speckle illuminations. We introduced the depthcorrection step where all E-field images taken at different depths within the coherence length of the light source were numerically propagated to the target depth. This increases the number of images to be used for constructing a time-reversal matrix at each target depth, which effectively increases the volumetric imaging speed.
All these benefits of using the compressed time-reversal matrix come with a price to pay. A finite overlap between random illumination channels introduces additive noise in addition to multiple scattering noise. Therefore, achievable imaging depth is reduced relative to the full sampling by the amount of sparse sampling-induced noise. Using orthogonal illumination channels such as the Hadamard patterns instead of unknown speckles can minimize the sparse sampling-induced noise at the expense of hardware simplicity. In case when a priori knowledge of the scene is known, the number of required measurements could be drastically reduced by introducing a learned sensing approach 29,30 using optimized illumination channels. Another drawback is that the achievable imaging resolution with the CTR-CLASS algorithm is diffraction limited. This is because, without knowledge of the illumination channels, the spatial cut-off frequency is solely determined by that of detection channels. The above shortcoming can be overcome by introducing a new image reconstruction algorithm combining the CTR-CLASS with methods that can reconstruct super-resolution images without prior knowledge of the illumination patterns, such as blind structured illumination microscopy 31 and random illumination microscopy 32,33 . In this study, ballistic waves scattered once by an object are used for image reconstruction, and multiple-scattered waves inside a scattering medium are considered as background noise. However, multiple-scattered waves do also carry spatial information of the object. CTR-CLASS algorithm can potentially be extended to make the deterministic use of multiple-scattered waves in image reconstruction for further reducing measurement time or lowering the achievable spatial resolution well below the diffraction limit 34 .
High-throughput volumetric imaging equipped with aberration correction capability for every depth section allows detailed mapping of microstructures deep within tissues. This will lead to accurate quantification of structural and molecular information in various biological systems. Therefore, the presented method will find its use for a wide range of studies in life science and medicine including the myelin-associated physiology in neuroscience, retinal pathology in ophthalmology and endoscopic disease diagnosis in internal organs. Due to the high-speed measurement of tissue aberration, it can also serve as wavefront sensing AO to provide aberration information for the hardware aberration correction. This will help to improve the imaging depth of fluorescence and nonlinear imaging modalities such as multi-photon microscopy, superresolution microscopy, and coherent Raman microscopy.

Materials and methods
Acquisition time for the CTR-matrix As long as the laser power is enough, the major factor determining the acquisition of the CTR-matrix is the frame rate of the camera, f cam . In this study, a laser power of about 40 mW was illuminated on mouse brain sample for volumetric imaging, which was sufficient to obtain high SNR camera images with an exposure time of 50 µs. The acquisition time set by the frame rate of the camera is given by T = M/f cam . The size of the camera pixel corresponds to 0.128 µm in sample space. In the proof-ofconcept experiment with a phantom aberrating layer, a series of camera images (1280 × 400 pixels) were acquired with a frame rate of 12,500 Hz and cropped to 312 × 312 pixels, which corresponds to a FOV of 40 × 40 µm 2 in sample space. For volumetric imaging of the mouse brain, a series of camera images (1280 × 1000 pixels) were acquired with a frame rate of 5000 Hz, which corresponds to a FOV of 164 × 128 µm 2 in sample space.

Animal preparation
Adult (over 8 weeks) C57BL/6 mice were deeply anesthetized with an intraperitoneal injection of ketamine/xylazine (100/10 mg/kg) and decapitated. After the scalp and skull were removed, the brain was fixed with 4% paraformaldehyde at 4°C overnight and washed with phosphate-buffered saline (PBS) three times. For imaging, the fixed brain was stuck to a plastic dish and immersed in PBS. All animal experiments were approved by the Korea University Institutional Animal Care & Use Committee (KUIACUC-2019-0024).