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

Hojun Lee,Seokchan Yoon,Pascal Loohuis,Jin Hee Hong,Sungsam Kang and Wonshik Choi✉

1Center for Molecular Spectroscopy and Dynamics,Institute for Basic Science,Seoul 02841,Korea

2Department of Physics,Korea University,Seoul 02841,Korea

Abstract 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 light1,2,delivering images3,and controlling transmission power4through 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 imaging5.The reflection matrix has also been exploited to focus light on a target embedded deep within strongly scattering media6-8.The wave correlation of the singlescattered waves in the reflection matrix was tailored to attenuate the effect of multiple light scattering9.A wavefront correction algorithm termed closed-loop accumulation of single scattering (CLASS)10was 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 resolution11.The singular value decomposition (SVD) was applied to a time-gated reflection matrix for retrieving a target image underneath strongly scattering media12,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 aberrations14.Indeed,the reflection matrix approaches provide formidable solutions in the context of computational adaptive optics (AO) microscopy15,16in 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 AO17,18that 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 skull14.

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 matrixRR†instead of the reflection matrixR.HereR†represents the conjugate transpose ofR.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µm2(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µm3(568×568×125 voxels) within a recording time of 3.58 s.

Results

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,

Fig.1 Schematics of imaging geometry and time-reversal process.

Here,Eo(ro;ri) is the time-gated E-field at position roon the output (detection) plane when a target object is illuminated by a point source located at position rion 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.Pi(o)(ri(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 property19,i.e.,Pi(o)(ri(o);r)=Pi(o)(ri(o)-r).Ems(ro;ri) 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 matrixRwhose element isEo(ro;ri) for a column index riand row index ro.Based on Eq.(1),Rcan be decomposed as

Here,Ois a diagonal matrix whose diagonal element isO(r).PiandPoare Toeplitz (diagonal-constant)matrices whose elements are respectively given by the input and output PSFs,Pi(ri;r) andPo(ro;r).Rmsis a matrix composed ofEms(ro;ri),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 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 relationPi(r′;r)=Po(r′;r).However,this is not a necessary condition in the present study.As we reported earlier10,the CLASS algorithm utilizingRis developed in such a way to find high-fidelity solutions for an unknown object function (O) and two PSFs (PiandPo) even in the presence of strong multiple scattering noiseRms.

Compressive sensing of the reflection matrixS

The matrixRin Eq.(2) can be directly measured by scanning the position riof the focused illumination and wide-field detection of the backscattered wave fieldEo(ro;ri) across ro.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 sizeL×L,the required number of sampling points for a complete sampling is given byN=L/Δxwhich is the total number of orthogonal spatial modes for the given FOV,NA,andλ.In fact,the full reflection matrixR∈CN×Ncan be measured by sending any complete setofNillumination fields,instead of point-by-point scanning with a focused beam.One can measure the respective output fieldEo(ro;j) for eachjthillumination and construct a reflection matrixRm∈CN×Nwhose columns are assigned by the measured output E-fieldsIn this case,the column and row indices ofRmarejand ro,respectively.Then,the measuredRmis expressed asRm=REi,whereEi∈CN×Nis an illumination matrix constructed byin the same way asRm.The reflection matrix of the sample can be obtained by multiplying the measured matrixRmby the inverse ofEi,i.e.,R=RmE-1i.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 ofMunknown random speckle illumination patterns,for compressive sensing ofR.In particular,we consider the case in whichMis significantly smaller thanN(M≪N).The time-gated E-field imageEo(ro;j) is recorded for eachjthspeckle illumination,and the sparsely-sampled reflection matrixRS∈CN×Mis then constructed usingas a matrix element (Fig.1b).Therefore,the matrixRSwith column indexjand row index rocan be expressed as

whereS∈CN×Mis the sensing matrix constructed byTo realize aberration correction and image reconstruction without a prior knowledge of the illumination pattern,we consider a CTR matrix,W=RsR†s.By inserting Eq.(3) inW,the matrix is expressed asW=(PoOPTiS)(S†P*iO†P†o)+Wms,whereWmsis the noise matrix associated with the multiple scatteringRms,and the superscript‘*’denotes the complex conjugate.Note that the covariance matrixSS†∈CN×Nis almost an identity matrixIfor sufficiently largeM.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σ[ri,r′i] are given by correlations of two series of random speckle fields illuminating different positions,riand r′i.Therefore,Wcan be expressed as

whereOIdenotes |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 relationPTiP*i=PiP†i=I,which is valid whenPi(ri;r)is a shift-invariant PSF induced by a phase-only pupil aberration.The last termWσ=PoOPTiσP*iO†P†odenotes the sparse sampling-induced noise associated with σ.Statistically,mean amplitude ofσ[ri,r′i] is given bywhere the bracket notation denotes an average over the variables in the subscript.S(ri;j) is normalized such thatTherefore,the magnitude of the matrix elements ofWσbecomes smaller asMincreases,which makesWσnegligible for sufficiently largeM.

Physical interpretation of the time-reversal matrixWis 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 operationR†S,a spherical wave(red curves)emitted from a point source at a position r′oon the output plane propagates in the backward direction through the object(P*iO†P†o) followed by a fictious scattering layer whose transmission matrix isS†.Afterwards,theRSis applied such that the reflected wave(blue curves)returns to the scattering layer(S)and the object(PoOPTi)in the forward direction to arrive at the output plane.Here,the important point is that the operation indicated by the shaded gray box (PTiSS†P*i)serves as a phase-conjugation mirror whenSS†≈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 throughPTiSS†P*i.This eliminates the need to consider the input aberrationPiand illumination patterns S.As a result,the matrixW=RSR†Scan be interpreted as a time-gated reflection matrix describing an imaging system that images a reflective object with the reflectanceOIthrough a scattering medium with an input transmission matrix ofP*oand an output transmission matrix ofPo.The whole process becomes anN-by-Nsquare matrix with its column and row indices both corresponding to ro.

There are a few major benefits of consideringWinstead ofR.First,the sensing matrixSdescribing the illumination patterns is removed inW,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,Wis greatly simplified such that onlyPoandOIremain to be identified.Imperfection in illumination and detection optics often causes discrepancy betweenPiandPoin the reflection geometry.PiandPoare intrinsically different in the transmission geometry.However,inW,it is not necessary to considerPi.Another critical benefit is the possibility of sparse sampling.The conditionSS†≈Isatisfies even whenMis 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 acoustics20-22and then has been intensively studied in microwaves23,24and optics7,12,25.In these previous studies,either iterative operation or the singular value decomposition (SVD) of the time-reversal matrixW=R†Rwas used for selective focusing on a few highly reflecting targets embedded in a scattering medium.Each input singular vector with a nonzero singular value ofWcorresponds 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 vectors12,13.

In contrast,we consider the CTR matrixW=RSR†Sand introduce a matrix decomposition ofWinto a product of three matrices,W=PoOIP†oin order to find the unknown object function (OI) embedded within a scattering medium inducing optical aberrations (Po).In conventional compressive sensing26,27,the sampling process is modeled as y=Sa+n,where y ∈RM×1is a vector of sparsely sampled data,a ∈RN×1is the original signal to be recovered,S∈RM×Nis a known sensing (or measurement) matrix withM≪N,and n denotes a measurement noise.Here,Mis the number of measurements andNis 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 sparsityS(the number of nonzero elements in a)is smaller than the number of measurements,M:M≥O(Slog(N/S)).

Conventional compressive sensing concerns the sparsity of the target a.In our study,the reflection matrixR∈CN×Nis sparsely sampled by the sensing matrixS∈CN×M:RS=RS.Thus,the degree of sparsity ofRis of concern.Essentially,our modelR=PoOIPTiallows us to treatRas a highly sparse matrix.SinceOIis diagonal andPo(i)is a Toeplitz matrix,Rcontains only 3Nunknowns.Since the measuredRScontains the number of elementsN×M,it forms a system ofN×Mequations with 3Nunknowns.Theoretically,there needsM≥O(3 log(N/3)) to accurately estimate the solutions.However,we used the unknown speckle illuminations,i.e.unknown sensing matrixS,and thus converted the problem to a decomposition of time reversal matrix,W=PoOIP†o+Wms+Wσ,whereWmsis due to multiple scattering noise,andWσis the additional sparsity-induced noise.SinceWσscales withthe reduction of samplingMgives rise to the increase of noise.This makes it necessary to setMlarger than the estimated value to properly decomposeW.

To identify the aberration and reconstruct an aberration-corrected image,the basis ofWis 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 matrix[k′o;ko] represents the planewave based transmission matrix between the object and output planes.When the output PSF is space-invariant,becomes a diagonal matrix whose elements are given by its complex pupil function(ko)=eiϕo(ko),where ϕo(ko) is the output pupil phase map.We consider here a phase-only pupil function that has amplitude of unity.The matrixis the target spectrum matrix in which each column consists of a shifted special-frequency spectrum of the target:

The CTR-CLASS algorithm10,acting on the CTR matrix,identifies an aberration correction matrixthat maximizes the total intensity of the object image(k)reconstructed from a corrected CTR matrix=Here,(k) is reconstructed by the sum ofalong diagonals:The algorithm iteratively finds the solution(ko)in the following way.At thenthiteration,thenthcorrection pupil functiontarget spectrum(Δk),and timereversal matrixare calculated as

where Δk=k′o-ko.Thenthcorrection phase angles(ko) is found by taking inner products of the angular spectrum images[ko+Δk;ko] and the corrected image(Δk).Note that the target spectrum(Δk)is reconstructed by synthesizing all the spatial frequency spectra covered by the input and output angles (|k/(2π)|,|/(2π)|≤2NA/λ),resulting in a special-frequency band of |Δk/(2π)|≤2NA/λ.The iteration starts with the initial conditions of(ko)=1 and=and continues until the root-mean-square(RMS)of the phase,σ2ϕ=〈|ϕ(n)c(ko)|2〉kobecomes less than a predefined value.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 matrixW.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 completeN-by-Nreflection matrix for a FOV havingNorthogonal modes,NE-field images must be acquired for an orthogonal set ofNillumination fields,where each image has a total ofNorthogonal pixels.However,we recorded onlyM(＜N) E-field images using random speckle illuminations to reduce the acquisition time and reconstruct a sparsely sampledN-by-Mreflection matrixRSas shown in Fig.2c.

Fig.2 Experimental setup for recording a compressed time-reversal matrix.

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 ofSS†≈I,we experimentally measured the S matrix by placing a mirror at the sample plane.Figure 2d showsSS†matrix obtained usingM=700 speckled illuminations for a FOV of 40×40µm2with 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µm2(N=88×88 pixels),Mspeckled 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 asCR=M/N.Figure 3a shows the CTR matrixWconstructed by speckle patterns forCR=1.In the absence of aberration (Po=I),the matrixW≈PoOIP†ois almost diagonal because it is reduced toOI.In the presence of aberration and scattering,the signal in the diagonal spreads out to the off-diagonal elements.The matrixin spatial-frequency domain and the identified correction phase ϕc(ko) by the CTR-CLASS algorithm are shown in Fig.3b,c,respectively.The corrected CTR matrixWc=P†cWPcis shown in Fig.3d.The intensity images before and after the aberration correction were reconstructed fromand,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αobefore aberration correction is given byαo=|Po(0)|2.Likewise,the Strehl ratioαcafter correction is given byαc=|Pres(0)|2,wherePres(r) is the residual PSF after the correction (Pres=P†cPo).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.

Fig.3 Wavefront aberration correction using time-reversal matrix.

Analysis of image SNR

We evaluated the performance of image recovery depending onCRto determine the minimum achievableCR.Reconstructed images and aberration phase maps for variousCRare shown in Fig.4a.The first row shows the reconstructed intensity images normalized byMforCR=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,althoughCRwas significantly reduced.Diffraction-limited high-resolution images could be successfully restored forCR≥0.02 (M≥155).Considering the camera frame rate of 12,500 Hz,it took only 12.4 ms to record the CTR matrixWforCR=0.02,setting the highest achievable aberration-correction image frame rate to 80 Hz for a FOV of 40×40µm2(88×88 pixels).If theCRwas 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 ofCR,so far as the image reconstruction was successful.This means that even the severe decrease inCRdoes 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 ofCR.This indicates that aberrations were properly corrected even with the smallM.The pupil phase retardation ϕ(cn)(ko) is estimated from the inner product of the angular spectrum image for the illumination angle of ko,[ko+Δk;ko] and the angular synthesized image(Δk)in Eq.(6),where the summation is taken over totalNspectral frequencies,Δk.Thus,the standard error of the estimated pupil phase scales withEssentially,the performance of aberration correction is determined byN,the number of sampled spectral frequencies or equivalently the image size.Therefore,as long asNis sufficiently large,the pupil phase retardation can be identified with high fidelity even for a small number of measurements,M.

Fig.4 Image reconstruction fidelity depending on the compression ratio(CR).

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 matrixWc=P†cWPc=P†cPoOIP†oPc+P†cWmsPc+P†cWσPc.Note that the incident speckle field is normalized as=1,which means the average intensity of speckle illumination fields at each pixel riis 1/M.For intuitive understanding,let us consider sending speckle illumination fields whose average intensity per pixel is 1.We define the signalscas the intensity of the target reconstructed after the correction,which corresponds to the main diagonal of the first termP†cPoOIP†oPc.Therefore,scis given bysc=MαcOI(r).There are mainly two kinds of noise sources:(1) measurement noiseWmsand (2) sparse sampling-induced noiseWσ.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,σmsremains almost unchanged after aberration correction and is proportional to the square root ofM:σms=Heremdenotes 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 byWσ=PoOσO†P†o.The standard deviation of the sparse sampling-induced noise,σsalso scales proportionally with the square root ofM:σs=where 〈OI(r)〉ris the average reflectance of the target.In fact,theσsis analogous to the measurement noise,except that it is proportional to the target’s average reflectance.Then,the SNRccan be estimated to be SNRc=

In the strong multiple scattering regime (OI(r)≪m),the SNRcof the reconstructed image from the CTR matrix is approximated as SNRc≈sc/σms=where SNR1=αoOI(r)/mis the SNR of a single wide-field image obtained without correction.In the weak scattering regime of SNR1≫1,the SNRcis mainly determined by the sparse sampling-induced noise and given by SNRc≈sc/σs=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 (OI(r)=1) and a fill factor value of～50% (〈OI(r)〉r=1/2),giving the ratioOI(r)/〈OI(r)〉r≈2.In this specific case of a binary target with 50% fill factor,the SNRcis expressed as SNR=In terms of compression ratio,it is expressed as SNR=We plotted the SNR as a function ofCRin Fig.4c (red dots)and observed that the SNR fits well withThe ratio between the diagonal and off-diagonal elements ofWcis the signal to background noise ratio,which in turn is equal to the SNR of the reconstructed image.We found that ratios between the diagonal to offdiagonal elements of bothSS†andWcalso fit well withas shown in Fig.4d,confirming the validity of our analysis.We found that this ratio fits well withconfirming the validity of our SNR analysis.Image reconstruction was successful up to aCRvalue of 0.02 (M≥155),but it was shown to fail with a further decrease ofCR.As explained above,the reduction inCRintroduces random noise.Below a certain thresholdCRvalue,the noise level becomes too high for the algorithm to retrieve the pupil phase map.It is difficult to generalize the minimum achievableCRvalue because it depends on various factors such as the SNR1of the uncorrected image,the distribution of the target spectrum,and the complexity of the aberration.Therefore,appropriate choice ofCRvalue 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 thez-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µm2(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µm2(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 depthz,E-field images taken within a range ofz±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 depthz.The reconstructed 3D volume image over 128×128×125µm3(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 images28.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 atz=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 atz=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 atz=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 atz=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 alongz-axis was～2µm.

Fig.5 Aberration-free volumetric high-speed imaging of a mouse brain.

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µm2at a frame rate of 80 Hz.Furthermore,we realized the volumetric imaging of a mouse brain over a volume of 128×128×125µm3in 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 approach29,30using 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 CTRCLASS with methods that can reconstruct super-resolution images without prior knowledge of the illumination patterns,such as blind structured illumination microscopy31and random illumination microscopy32,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 limit34.

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,fcam.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 byT=M/fcam.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µm2in 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µm2in 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).

Acknowledgements

This work was supported by the Institute for Basic Science (IBS-R023-D1).

Author details

1Center for Molecular Spectroscopy and Dynamics,Institute for Basic Science,Seoul 02841,Korea.2Department of Physics,Korea University,Seoul 02841,Korea.3Department of Applied Mathematics,University of Twente,Drienerlolaan 5,7522 NB Enschede,Netherlands.4Achmea Holding BV,Handelsweg 2,3707 NH Zeist,Netherlands

Author contributions

S.Y.conceived the idea.W.C.supervised the project.H.L.and S.Y.performed the experiments.P.L.,H.L.,and S.Y.wrote the Matlab code and analyzed the experimental data.J.H.H.prepared biological samples.W.C.,S.Y.,H.L.,and S.K.wrote the manuscript.All authors discussed the results and commented on the paper.

Data availability

All relevant data are available from the authors upon request.

Conflict of interest

The authors declare no competing interests.

Supplementary informationThe online version contains supplementary material available at https://doi.org/10.1038/s41377-021-00705-4.