Multi-focus averaging for multiple scattering suppression in optical coherence tomography

Multiple scattering is one of the main factors that limits the penetration depth of optical coherence tomography (OCT) in scattering samples. We propose a method termed multi-focus averaging (MFA) to suppress the multiple-scattering signals and improve the image contrast of OCT in deep regions. The MFA method captures multiple OCT volumes with various focal positions and averages them in complex form after correcting the varying defocus through computational refocusing. Because the multiple-scattering takes different trajectories among the different focal position configurations, this averaging suppresses the multiple-scattering signal. Meanwhile, the single-scattering takes a consistent trajectory regardless of the focal position configuration and is not suppressed. Hence, the MFA method improves the ratio between the single-scattering signal and multiple-scattering signal, resulting in an enhancement in the image contrast. A scattering phantom and a postmortem zebrafish were measured to validate the proposed method. The results showed that the contrast of intensity images of both the phantom and zebrafish were improved using the MFA method, such that they were better than the contrast provided by the standard single focus averaging method. The MFA method provides a cost-effective solution for contrast enhancement through multiple-scattering reduction in tissue imaging using OCT systems.


Introduction
Deep tissue imaging has long been of strong interest in biomedical optics [1].Optical coherence tomography (OCT) provides three-dimensional images of a sample with a-few-millimeter penetration and few-to-few-tens micrometer resolution, and has become widely used for the non-invasive imaging of biological samples.OCT has been successfully adapted to investigations of various tissues, such as human skin [2], ex vivo brain tissue [3,4], and in vitro tumor spheroids [5][6][7].
In general, an OCT image is primarily contributed from the single-back-scattering (SS) signals of the sample, whereas the contribution from the multiple-scattering (MS) signals is trivial.The SS signals are retrieved from the photons that have been back-scattered once at the sample object and hence directly carry the object information.However, when imaging in scattering samples such as tissues, the MS signals become dominant and harm the image contrast especially in deep regions [8].The MS-originated signals at a position in the image are not scattered at the corresponding position in the sample and thus do not convey the correct object information.MS limits the resolvable imaging depth of OCT in tissues and hampers the visualization of deep microstructures [9,10].
In standard OCT with a confocal configuration, most of the MS signals can be rejected by the confocal gating.However, the residual MS signals still disturb the imaging.Techniques for further reducing the MS effect have been explored.One approach is modulating the probing light to decorrelate the MS signals using wavefront manipulation devices, such as a spatial light modulator (SLM) [11,12] or a deformable mirror [13].However, this approach may result in high system complexity and cost.A solution with inexpensive optics that can be easily implemented in the standard OCT scheme is preferable.
In this paper, we propose a method termed "multi-focus averaging" (MFA), to reduce the MS signals and improve the image contrast in deep regions of scattering samples.The MFA method uses a low-cost electrical tunable lens to decorrelate the MS signals over multiple volumetric acquisitions.Here, the depth position of the focus is modulated among the multiple volumes.A computational refocusing is then adopted to cancel out the different defocus in the focus-modulated OCT volumes.A scattering phantom and a postmortem zebrafish were used for validation.

Principle of the MFA method
The MFA method is a combination of sequential volumetric OCT imaging with different defocus, computational refocusing of each volume, and complex averaging of the volumes.In an OCT volume of a scattering sample, the complex en face OCT signal S(x, y; z, z d ) at an imaging depth of z comprises two components, namely S SS (x, y; z, z d ) and S MS (x, y; z, z d ) as S(x, y; z, z d ) = S SS (x, y; z, z d ) + S MS (x, y; z, z d ), where x and y are the lateral positions of the scanning locations, and z d is the amount of defocus at the imaging depth z. S SS and S MS are the signal components originating from SS and MS photons, respectively.According to the formulation of Ralston et al. [14], the S SS component can be expressed as a convolution of a defocus-free OCT SS signal and a depth-and defocus-dependent quadratic phase function ϕ as S(x, y; z, z d ) = S SS (x, y; z, 0) * exp [iϕ(x, y; z, z d )] + S MS (x, y; z, z d ), where S SS (x, y; z, 0) is the defocus-free SS OCT signal.The SS photons experience only a single scattering event, and hence, they take consistent trajectories among the multiple measurements with different defocus, as depicted by the SS trajectories in Fig. 1.Hence, once the defocus is corrected by computational refocusing, the phase of SS components become consistent among the volumes.
Meanwhile, the MS photons can be scattered by different combinations of scatterers if the defocus changes the MS trajectories as shown in Fig. 1.Hence, the trajectories of the MS photons are altered by the defocus.This alteration of the trajectories results in an inconsistent phase of S MS components for the different defocus, even after computational refocusing.These properties of the signal components are formulated as follows.
Computational refocusing can be described as the deconvolution of the quadratic phase function ϕ.After the computational refocusing, the OCT signal becomes where the deconvolution is represented as the convolution with the complex conjugate of ϕ (see the Appendix for details of this representation).No matter the original defocus amount, the computational refocusing provides the consistent defocus-free SS component.In addition, the signal strength of the SS component is enhanced by the refocusing.Meanwhile, the phase of the MS component depends on the original defocus amount because the MS trajectories depend on the defocus amount.Furthermore, the scattering events of the MS photons occurring at a shallower depth is observed at a deeper depth in the image, while the scattering event of an SS photon is observed at the depth of the event.And hence, the MS and SS signals observed at the same image depth may have occurred at the different depths in the sample, and may have different defocus amounts.Therefore, the computational refocusing tailored for the SS component does not enhance the signal strength of the MS component.
In our MFA method, an MS-signal-suppressed OCT volume is obtained by complex averaging the multiple volumes with different defocus amounts after the computational refocusing as where N is the number of volumes being averaged and z d,j is the defocus amount of the j-th volume.As the phase of the MS term (the second term) is unpredictable and practically random, the amplitude of this term is reduced by the averaging.Meanwhile, the phase of SS term (the first term) is consistent among the volumes, and hence, its amplitude is not reduced.In practice, a larger difference in the defocus amounts among the volumes may cause higher mutual randomness (i.e., higher decorrelation) of the phases of the MS components, and a better suppression of MS components can thus be achieved.This hypothesis is experimentally validated in Section 4.2.

OCT setup
In this study, we used a polarization sensitive OCT (PS-OCT) system, which has been described in Refs.[15,16].However, we only used a single-polarization-channel OCT image (i.e., a conventional non-polarization sensitive OCT image) to demonstrate the MFA method.The light source was a sweeping laser source (AXP50124-8, Axsun Technologies, MA) with a central wavelength of 1.31-µm, a bandwidth of 106 nm, and a scanning rate of 50 kHz.The objective (LSM03, Thorlabs, NJ) had an effective focal length of 36 mm and gave a depth of focus (DOF) of 0.36 mm in air.The lateral and axial resolutions were measured to be respectively 18 µm and 14 µm in tissue.The influences from the focus tuning to the resolutions will be discussed in Section 4.3.4.
An electrical tunable lens (ETL, EL-10-30-CI-NIR-LD-MV, Optotune Switzerland AG, Switzerland) was set on the sample arm to modulate the focal position.It is noteworthy that the ETL can be easily aligned and applied to standard OCT systems.

Core methods of the MFA
In our MFA method, multiple complex OCT volumes are acquired with different focal positions.For each volume, the bulk phase error is estimated and corrected using a smart-integration-path [17] method.
The defocus is corrected in each volume by applying computational refocusing, where a series of phase-only deconvolution filters are applied at all depths to correct the defocus.Details of the refocusing method can be found in Section 2.2 of Ref. [18].
When modulating the focal position using an ETL, the deformation of the lens affects the optical path length of the probe beam.This introduces an axial image shift, which sometimes reaches a few tens of microns among the focus-modulated volumes.In our MFA method, the axial shift is estimated and corrected using a sub-pixel intensity cross-correlation method [19].In this process, B-scans from different volumes at the same location are extracted, and four-fold up-sampling along the depth direction is conducted.Intensity-based cross-correlations among the extracted B-scans are performed to estimate shifts with resolution of 1/4 of a pixel along the axial direction.As this image shift is consistent within each volume, the estimated shifts are used for co-registration of all B-scans in each volume.
The inter-volume phase offsets are also corrected.The offsets are computed for each A-line using an intensity-weighted phase difference averaging.In this averaging, one volume is used as a reference, and the product of the reference and the complex conjugate of another volume is computed.The product is then averaged in complex form along the depth direction within a certain depth range, which was 30 pixels with sufficient intensity in the present implementation.The phase of the averaged signal is the phase offset.The phase offsets are then computed for each A-line and corrected.
Finally, all volumes are averaged in complex form to obtain an "MFA volume."The OCT intensity volume is then generated from the MFA volume.

Samples
Two types of sample were measured to validate the performance of the MFA method.One was the scattering phantom illustrated in Fig. 2(a), which comprises glass slips (A, D), a scattering layer (B), and a glass plate buried at the bottom of the scattering layer (C).The glass slip at the top (A) was tilted to prevent a strong reflection from the surface.The scattering layer (B) was a mixture of 0.025 mL of 10%-concentration polystyrene micro-particles (diameter of 10 µm, 72968-10ML-F, Sigma-Aldrich) and 0.6 mL of ultrasound gel (Pro Jelly, Jex, Japan).The glass plate (C) had a thickness of approximately 0.12 mm and was buried at the bottom of the scattering layer to create a space without scattering.The thickness of the scattering layer above the glass plate was approximately 1.5 mm.A photograph of the sample from the top is shown in Fig. 2

(b).
A 30-dpf wild-type postmortem female zebrafish [Fig.2(c)] was imaged to demonstrate the measurement of a biological sample.The zebrafish was anesthetized by tricaine and sacrificed by low-temperature treatment (placed in ice for 2 minutes) and then placed in a petri dish and immersed in saline solution for measurement.A piece of black tape was placed on the bottom of the petri dish to prevent strong specular reflection.A wide-field-of-view en face OCT-intensity projection of the zebrafish is shown in Fig. 2(d), where the orange box indicates the measured area.The postmortem zebrafish measurement was performed following the animal experiment guidelines of the University of Tsukuba.The measurement protocol was approved by the Institutional Animal Care and Use Committee (IACUC) of University of Tsukuba.

Measurement protocol
In both the phantom and zebrafish measurements, the samples were placed on a linear translation stage and measured.Seven volumes with different focal positions were acquired for averaging.The interval between focal positions is referred to as the focus shifting step ∆z and was set at 0.12 mm (approximately 1/3 of the DOF).The overall shifting distance of the focus is denoted as D = (N − 1)∆z and was 0.72 mm in the study.These parameters were determined empirically, and details of the parameter selection strategy are given in Section 4.2.For comparison, another set of raw volumes were acquired and averaged without focus modulation.This averaged volume without focus modulation is referred as the "single focus averaging" volume, in short "SFA".
In the measurements, the focal positions were set in the deep regions of the samples to enhance the light collection efficiency.The lateral scanning range was 1.5 mm × 1.5 mm sampled with 512 × 512 A-lines, which gave an isotropic lateral pixel separation of 2.93 µm (around 1/6 of the lateral optical resolution).The acquisition time for each volume was approximately 6.5 s.

Quantitative image contrast analysis
We defined the signal-to-background ratio (SBR) to quantify the image contrast in the scattering phantom measurement.An en face image at a depth slightly beneath the top surface of the buried glass plate was extracted.Several scatterers with similar intensities were manually selected to compute the SBR.The signal intensity of each scatterer was defined as the averaged image intensity within a 3-pixel × 3-pixel window centered at the scatterer.The signal intensity used to compute the SBR was then defined as the mean of the signal intensities of the selected scatterers.The background intensity was defined as the mean intensity of pixels in the buried glass plate region (i.e., the scattering-free region).
It is noteworthy that the imaging depth is difficult to be defined and quantified for practical samples.Hence, here we used SBR to quantify the performance of MFA instead of directly quantifying the imaging depth.The focal positions were estimated by finding the depth with zero defocus from the linear regression of defocus amounts estimations.The defocus estimations were performed at each depth by finding the phase-only defocus-correcting deconvolution filter that minimizes the information entropy of the en face image.An example of such linear regression can be found in Section 5.6.1 in Ref. [20].Among the B-scans, the MFA image provides a lower background noise in the scattering layer than the other images, especially in the deep regions (denoted by the yellow dashed boxes).

Scattering phantom
Figure 3(g) presents the intensity depth profiles obtained using the three methods.The depth intensity profile of the MFA (solid orange line) shows a lower intensity than the other profiles as depth increasing (the region with yellow background).Since the MS signal is expected to be more pronounced in deeper regions, this profile may indicate the MS suppression by MFA method.Meanwhile, the single acquisition (dotted black line) and SFA (dashed blue line) profiles have intensities similar to each other, which indicates that the SFA method does not appreciably reduce MS.At the depths of the buried glass plate region (blue background region) in Fig. 3(g), we see that the MFA and SFA methods show appreciably lower signal intensity than that of the single acquisition method.It may be because of the suppression of measurement noise by the averaging.At a depth close to the top surface of the buried glass plate (indicated by an orange arrowhead), the SFA curve has a lifted "tail" with intensity higher than that in the air.As the glass plate region is scattering-free, this tail might be due not to scattering in the glass but to the MS in the superior scattering part of the phantom.A similar result was reported by Yadlowsky et al. [8].This lifted tail is reduced when using the MFA method, which further supports the MS-suppression ability of MFA.Three more depth profiles using the same number of A-lines but at different locations are shown in Supplementary Fig. S1, and all of them gave consistent results.These findings suggest that the SFA method reduces the system noise, and the MFA method reduces both system noise and MS signals.
Figure 3(d)-(f) show the en face images at the same depth indicated by the orange arrowheads in Figs.3(a) and (g).The SFA image has slightly better contrast of the scatterers than the single acquisition image, but this contrast is not as good as that of MFA image (magnifications in the yellow boxes).The MFA method provides better image contrast than the SFA method because the former benefits from the MS reduction.In addition, we also noticed some "low-intensity scatterers" appeared in the single acquisition and SFA images [Figs.3(d) and (e), respectively] (arrow in the magnification insets), and became dimmer by applying MFA [Fig.3(f)].We suspect these "low-intensity scatterers" might be the speckle caused by the MS signal, because the SFA method only suppresses the system noise but not the MS signal, whereas the MFA method suppresses both.
The SBR was computed from Fig. 3(d)-(f) for quantitative comparison.We selected ten scatterers that have the highest intensities [red circles in Fig. 3(d)] from all the scatterers in the en face images to compute the signal intensity.Since these scatterers appeared clearly, we expected the signals of these scatterers were dominated by SS.The background intensity was computed from the pixels in the glass plate region [blue line in Fig. 3(d)].The computed SBRs are 30.12,32.13, and 34.64 dB for single acquisition, SFA, and MFA images, respectively.Adopting the current protocol, the MFA method provides an SBR improvement of 4.5 dB relative to the single acquisition, and 2.5 dB relative to the SFA method.  .This low-scattering-intensity structure is expected to be the notochord according to its anatomic features because a similar transparent structure of the notochord has been observed at 5 dpf [21].At depth Z2 [Fig.4(g)-(i)], the contrast between the myosepta (black arrows) and the surrounding muscle tissues is again highest for the MFA method among the three methods.Some hyper-scattering spots in the belly region (yellow dashed circles) are sharp and recognizable in the MFA image, whereas they are blurred in the SFA image.These hyper-scattering spots are almost unrecognizable in the single acquisition image.

Postmortem zebrafish
At a deeper depth Z3, a thin hyper-scattering tissue is noted in the MFA image [green arrowheads in Fig. 4(l)] but difficult to recognize in the single acquisition and the SFA images [Figs.4(j) and (k)].This tissue is also observed at the upper depth Z1 [green arrowheads in Fig. 4(d)-(f)].A similar structure has been visualized in adult zebrafish by OCT [22], and it is expected to be the outer layer of the swim bladder.Unidentified structures elongated along the dorsal side of the zebrafish [white ellipses in Fig. 4(j)-(l)] are also visualized with the highest contrast in the MFA image.We also computed the SBR from Fig. 4(j)-(l) for image contrast quantification, where the swim bladder layer structure was selected as the signal and the air region was the background.The signal region was selected by a window with a size of 15 × 15 pixels as indicated by the boxes in Fig. 4(j)-(l).The computed SBRs are 11.7, 14.4, and 16.1 dB for single acquisition, SFA, and MFA images, respectively.
These findings suggest that the MFA method is applicable to scattering tissue imaging, reduces the MS signal, and enables better visualization of structures.

Validation of the induced defocus and its correction by the MFA method
To validate the defocus correction over multiple acquisitions, we compared the lateral resolutions of a set of seven OCT volumes with different focal positions.These volumes are the same phantom volumes used for the MFA method demonstration in Section 3.1 (Fig. 3).Figures 5(a  The lateral resolution at each depth was evaluated by the speckle size of the en face intensity image, which was estimated using a linear-intensity-based auto-correlation.The speckle size was defined as the full width at half-maximum (FWHM) of the auto-correlation function.This estimation was performed along both fast-scanning (denoted as "X") and slow-scanning (denoted as "Y") directions.Examples of the auto-correlation functions are shown in Figs.5(c) and (d), whose FWHMs are used to determine the speckle size along both directions.A depth range of 1.1 mm (150 pixels) beneath the surface of the scattering layer was used for this validation.
Without refocusing, the seven volumes with different defocus have speckle sizes (X and Y FWHMs) varying along the depth [Fig.5(e), first row].After applying refocusing to each volume, the FWHMs become almost constant at each depth [Fig.5(e), second row].Note that the horizontal lines in the plot are not regression or theoretical curves but plotted dots appearing as lines because the FWHMs are highly constant.Most of the FWHMs converge on a value of 14.6 µm, whereas the system optical lateral resolution (defined as the focus radius where the intensity falls to 1/e 2 of the maximum) is 18 µm.This suggests that computational refocusing corrects the defocus over a depth range exceeding 1 mm, regardless of how much defocus is applied in the measurement.

Experimental optimization of MFA parameters
To find the optimal protocol, we experimentally explored the relationship between the number of averaged volumes N, the overall focus shifting distance D, the focus shifting step ∆z, and the SBR.In this experiment, we performed five sets of sequential volume measurements with different ∆z as summarized in Table 1.For each set of measurements, we extracted subsets to examine several configurations of N. D is determined from ∆z and N as D = (N − 1)∆z.The sample used for this measurement was the phantom shown in Fig. 3.As there are two independent parameters, we plot the SBRs in a three-dimensional space of the SBR, D, and N as shown in Fig. 6.The color gradient of the plotted points represents the SBR.The SBR takes a minimum value when N and D are close to zero and increases with both N and D. The black plot is the projection along the N-axis that shows the relationship between the SBR and D. We see that the SBR saturates at D approximately twice the DOF and starts to drop at D exceeding approximately 4 times the DOF.The green plot shows another projection along the D-axis that gives the relationship between the SBR and N. It is seen that the datasets with different ∆z saturate at different N.For example, the datasets of configuration #1 (∆z = 1 × DOF) and configuration #4 (∆z = 1/4 × DOF) saturate at approximately N = 4 and N = 8, respectively.These results suggest that both the parameters N and D contribute to the SBR improvement in the MFA method.
The optimal parameters were selected based on the above results and used in the experiments described in the previous sections.We first selected D as twice the DOF, where the projection curve starts saturating (black plots in Fig. 6).We then performed measurements of the zebrafish sample with ∆z of 1, 1/2, 1/3, and 1/4.As D was fixed to be twice the DOF, N was 3, 5, 7, and 9, respectively.To quantify the image contrast, we computed the SBRs from the region of a scattering spot and a region with low-scattering intensity (indicated by the red boxes in the magnified insets of Fig. 7, with a size of 15 × 15 pixels).Results show that the parameters correspond to ∆z = 1/3, D = 2, N = 7 (third column) provide the highest SBR (17.2 dB) in comparison with the other protocols.We therefore selected this parameter set as our optimal protocol and used it to perform the measurements presented in the results section.This configuration was thus selected as the optimal parameter set in the present study.

Interpretation of the SBR drop
In the protocol optimization, we note that the SBR drops when D becomes too large.A possible reason of this SBR drop is the uncorrected confocality; i.e., the signal amplitude drops when the sample plane is far from the focal plane.When D exceeds a few DOFs, the focus shifts far from the depth of interest.As a result, the amplitudes of both SS and MS signals decrease even with the computational refocusing, whereas the amplitude of the measurement noise is unchanged regardless of the focal position.The solution to overcome this limit remains an open issue and requires further investigation.Here, the astigmatism was introduced to the illumination beam, and its astigmatic angle was modulated to decorrelate the MS among the measurements [13].In addition, a deformable membrane mirror has been combined with a full-field SS-OCT, and fast three-dimensional volumetric cross-talk-free imaging has been demonstrated [23,24].Borycki et al. used a spatial phase modulator for digital aberration correction and MS reduction [25].
An advantage of the ETL over these wavefront manipulation devices is the low cost of the ETL.A simple comparison of the aforementioned devices is shown in Table 2. Another advantage of the ETL is that the ETL is easily integrated with standard OCT systems because it is a refractive (not reflective) device.The optical design schematics of the sample arm without and with the mounting of the ETL in our system are shown in Fig. 8(a) and (b), respectively.Figure 8(c) is a photograph of our sample arm, where the ETL (shown by red dashed box) was inserted between the collimator and the mount of the galvanometric scanner.In our implementation, the replacement of the ETL and the realignment of the system usually takes only a few minutes.This simple implementation of the ETL makes our probe arm switchable between the standard OCT mode and MFA-compatible mode.Wavelength-dependent transmission of the ETL should be considered.According to the device specification sheet, the ETL provides approximately 85% transmission (−0.7 dB) for current probe wavelength of 1,310 nm.In practice, the single-pass power loss was measured to be −0.5 dB, leading to an overall −1 dB power loss for the double-passing scheme.In future applications, this loss can be minimized by using another ETL optimized for our probe wavelength.

Temperature dependency
The ETL's accuracy for focal position modulation could be affected by the change of temperature.The guaranteed tuning range, according to the official technical documents of the device, may be reduced by 0.07 diopter per degree.To maintain the tuning consistency of the ETL, we roughly controlled the room temperature to be around 23 to 24 degree Celsius by a standard air conditioner.Hence, the temperature-related effect of the ETL to our system was negligible.Additional calibration might be required if we perform the measurement in some extreme room temperatures.

Influence on axial and lateral resolutions
To investigate the effect of focus modulation on the axial resolution, we measured the axial resolution by analyzing the surface signal of a gold-coated reflection mirror set at the focus.We found the focus variation by ETL does not have significant influence on the axial resolution across the full tuning range of ETL.
To investigate the effect of focus modulation on the lateral resolution, a numerical investigation on the focus spot size was conducted.The focus spot size (1/e 2 -width) was evaluated by numerical simulation using Optic Studio (version 19.4 SP1 Professional, Ansys Zemax, USA).The focus spot size was simulated with several focal positions varying in the range of approximately 6 mm.In addition, the spot size without ETL implementation was also computed.The results are plotted in Fig. 8(d), which shows a reduction in spot size as the focal position is shifting away from the objective.
However, in our current implementation, the focal position was shifted over a short total distance, i.e., 0.72 mm (as indicated by the blue background area).It results in only a slight change of the spot size that is less than 0.5 µm.Hence, we consider that the effect of focus modulation on the lateral resolution was negligible in the current configuration.Careful examination and potential correction of focal spot variations might be necessary for some future cases.

Future extension for PS-OCT
PS-OCT is a functional extension of OCT that provides additional polarization contrast for the quantitative evaluation of the optical properties of tissues [26][27][28].It is also known that MS affects the polarization measurement of PS-OCT [27].Adie et al. reported that the presence of MS randomized the polarization states [29], leading to so-called depolarization.Several researchers suspected that MS generates artifacts in accumulative or local phase retardation measurements in the imaging of biological samples [30,31].
However, the MS issue in PS-OCT has rarely been addressed.One example was demonstrated by Gao et al. in FF-OCT, who used a Muller-matrix-based MS subtraction process to remove the MS-induced local phase retardation [32].
The MFA method is applicable also to PS-OCT.Our preliminary study showed that in a scattering phantom, the measured degree-of-polarization uniformity (DOPU) decreases in deep regions, which can be mitigated by adopting the MFA method [33].We are currently working on adapting the MFA method to PS-OCT imaging.This adaption may improve the accuracy of quantitative polarization measurements in biological tissues.

Acquisition and computation time of current MFA
For future application of MFA to in vivo measurements, we need to consider acquisition and computation times.
The acquisition time of MFA can significantly affect the feasibility of in-vivo application.In the current implementation, the total acquisition time for seven volumes being averaged is approximately 45.5 s.It is hard to suppress the sample motion and achieve phase-stable OCT measurement for such a long measurement time.One possible solution is MFA based on one-dimensional computational refocusing as we will discuss in the next Section 4.5.2.
The data processing time is less significant than the acquisition time.However, shorter computation time is preferable for several applications.In the current implementation, the total computation time for the refocusing, axial shift correction, and phase offset correction for seven volumes are approximately 385 s, 175 s, and 270 s, respectively, with a CPU based Python program on a laptop PC (CPU: Intel Core i7-8750H, memory: 32 GB).Since both the computational refocusing and the shift correction heavily use numerical Fourier transform, the computation time can be significantly shortened by introducing GPU based processing.

Feasibility of in vivo application
The ETL enables fast and accurate focus modulation with a mechanical response time of 15 ms [34].However, the current MFA method captures multiple volumes for averaging, which is time consuming and a challenging task for in vivo imaging.
Nevertheless, the long acquisition time can be overcome by several means.One possible solution is to use ultra-high-speed acquisition OCT [35][36][37] and a motion tracking module to reduce the effect of sample motion.Another possible solution is to use B-scan-based focus modulation instead of modulating the focus for each volume.In other words, we can capture multiple B-scans at the same location with different focal positions, applying a one-dimensional version of computational refocusing [38], and average the focus-corrected B-scans to generate an MS-reduced B-scan.Finally, the three-dimensional volume can be obtained by repeating this process for a volume.In this case, the practical acquisition time, in which the measurement should be stable, can be significantly shorter than that of the original MFA method.Namely, the practical acquisition time of the new method is the time to acquire a set of B-scans, while that of the original method is a time to acquire a set of volumes.A proof-of-concept experiment showed that this method achieves a level of MS suppression similar to that achieved by the current MFA method [39].

Conclusion
We have developed a new method termed "MFA", to suppress the MS signals and improve the image contrast in OCT.A scattering phantom was measured to validate the MS reduction of the MFA method.The contrast improvement of the MFA image was compared with single acquisition and single focus averaging images quantitatively, showing the MFA method achieved a 4.5-dB enhancement of the signal-to-background ratio (SBR) to the single acquisition method.A postmortem zebrafish was measured to demonstrate the ability of the MFA method in biological imaging.Images of the zebrafish showed the MFA method better visualized anatomic structures in deep regions.We expect this proposed method will help better visualize the anatomic features.MFA method has great potential of being adopted in in vivo imaging and other imaging modalities to reduce the multiple-scattering effect.

Fig. 1 .
Fig. 1.Schematics of trajectories of SS and MS photons with different focus depths; i.e., different defocus amounts.No matter the focal position, the SS photon is scattered only once and hence has the same path length.This results in the consistent phase of the SS signal after computational refocusing.Meanwhile, the trajectories of MS photons are scrambled by changes in the focal position.Hence, the phase of the MS signal is randomized even after the computational refocusing.

Fig. 2 .
Fig. 2. Schematic (a) and photograph (b) of a scattering phantom.The phantom comprises glass slips A and D, a scattering layer B, which is a mixture of polystyrene micro-particles and ultrasound gel, and a glass plate C embedded in the scattering layer.The glass plate C provides a scattering-free area.A postmortem zebrafish at 30 days post fertilization (dpf) was used as a biological sample.The sample is shown in a color photograph (c) and a wide-field-of-view OCT intensity projection (d).The orange box around the belly region denotes the measurement area in the validation study of the MFA method.

Figure 3 (
Figure 3(a)-(c) shows the OCT intensity B-scans of the single acquisition, SFA, and MFA methods, respectively.The approximate focal positions are denoted by the orange arrowheads.Multiple orange arrowheads in Fig.3(c) denote different focal positions of each image used to generate the MFA image.The focal positions were estimated by finding the depth with zero defocus from the linear regression of defocus amounts estimations.The defocus estimations were performed at each depth by finding the phase-only defocus-correcting deconvolution filter that minimizes the information entropy of the en face image.An example of such linear regression can be found in Section 5.6.1 in Ref.[20].Among the B-scans, the MFA image provides a lower background noise in the scattering layer than the other images, especially in the deep regions (denoted by the yellow dashed boxes).Figure3(g) presents the intensity depth profiles obtained using the three methods.The depth intensity profile of the MFA (solid orange line) shows a lower intensity than the other profiles as depth increasing (the region with yellow background).Since the MS signal is expected to be more pronounced in deeper regions, this profile may indicate the MS suppression by MFA method.Meanwhile, the single acquisition (dotted black line) and SFA (dashed blue line) profiles have intensities similar to each other, which indicates that the SFA method does not appreciably reduce MS.At the depths of the buried glass plate region (blue background region) in Fig.3(g), we see that the MFA and SFA methods show appreciably lower signal intensity than that of the single acquisition method.It may be because of the suppression of measurement noise by the averaging.At a depth close to the top surface of the buried glass plate (indicated by an orange arrowhead), the SFA curve has a lifted "tail" with intensity higher than that in the air.As the glass plate region is scattering-free, this tail might be due not to scattering in the glass but to the MS in the superior scattering part of the phantom.A similar result was reported by Yadlowsky et al.[8].This lifted tail is reduced when using the MFA method, which further supports the MS-suppression ability of MFA.Three more depth profiles using the same number of A-lines but at different locations are shown in Supplementary Fig.S1, and all of them gave consistent results.These findings suggest that the SFA method reduces the system noise, and the MFA method reduces both system noise and MS signals.Figure3(d)-(f) show the en face images at the same depth indicated by the orange arrowheads in Figs.3(a) and (g).The SFA image has slightly better contrast of the scatterers than the single acquisition image, but this contrast is not as good as that of MFA image (magnifications in the yellow boxes).The MFA method provides better image contrast than the SFA method because the former benefits from the MS reduction.In addition, we also noticed some "low-intensity scatterers" appeared in the single acquisition and SFA images [Figs.3(d)and (e), respectively] (arrow in the magnification insets), and became dimmer by applying MFA [Fig.3(f)].We suspect these "low-intensity scatterers" might be the speckle caused by the MS signal, because the SFA method only suppresses the system noise but not the MS signal, whereas the MFA method suppresses both.The SBR was computed from Fig.3(d)-(f) for quantitative comparison.We selected ten scatterers that have the highest intensities [red circles in Fig.3(d)] from all the scatterers in the en face images to compute the signal intensity.Since these scatterers appeared clearly, we expected the signals of these scatterers were dominated by SS.The background intensity was computed from the pixels in the glass plate region [blue line in Fig.3(d)].The computed SBRs are 30.12,32.13, and 34.64 dB for single acquisition, SFA, and MFA images, respectively.Adopting the current protocol, the MFA method provides an SBR improvement of 4.5 dB relative to the single acquisition, and 2.5 dB relative to the SFA method.

Fig. 3 .
Fig. 3. (a)-(c) and (d)-(f) show the B-scans and en face images of the phantom, respectively.In (a)-(c), the orange arrowheads denote the approximate focal positions, and the yellow boxes denote a deep region of the scattering layer.In (d)-(f), the magnification insets are the en face images of a small region in the scattering layer containing several scatterers.In (d), the red circles denote the scatterers selected to compute the signal intensity.The blue line indicates a region with 2 × 200 A-lines in the glass plate region used to compute the background intensity for the SBR analysis, and their locations in the B-scans are indicated by the blue bracket in (a).White arrows in (d)-(f) indicate a low-intensity scatterer signal, which is dimmed in (f).It could be the speckle caused by the MS signal.(g) shows the intensity depth profiles that are averaged by the A-lines at the en face locations denoted by the blue line in (d).The 0-mm depth refers to the top surface of the cover glass on the top, and the green arrowhead denotes the depth in air where the intensity is 0 dB.The orange arrowheads in (a) and (g) denote the depth where the en face images are taken.

Figure 4 (
Figure 4(a)-(c) shows the B-scans of the sample.The MFA images better visualize structural boundaries [orange magnification inset in Fig. 4(c)], which are smeared by strong noise in Figs.4(a) and (b).This low-scattering-intensity structure is expected to be the notochord

Fig. 4 .
Fig. 4. Cross sectional and en face images of the postmortem zebrafish.In (a)-(c), the orange arrowheads denote the approximate focal positions, and the magnification insets are supposed to show notochord structure.The black arrows in (a)-(c) and (g)-(i) may indicate the myosepta.In (d)-(f), white dashed areas denote the belly region, and blue dashed ellipses mark fine structures observed in the muscle region.The green arrowheads in (d)-(f) and (j)-(l) may indicate the outer layer of the swim bladder.In (g)-(i), yellow dashed circles denote several hyper-scattering spots observed in the belly region.In (j)-(l), white ellipses denote unidentified structures are better contrasted in the MFA image, and small red boxes indicate two selected regions for computing SBR.At all depths, the MFA method has the images with best contrast among the three methods.

Figure 4 (
Figure 4(d)-(f) shows the en face images at depth Z1.The MFA image provides better contrast for some fine structures in the muscle region (within blue dashed ellipses) and some hyper-scattering spots in the belly region (top-left dashed areas).At depth Z2 [Fig.4(g)-(i)], the contrast between the myosepta (black arrows) and the surrounding muscle tissues is again highest for the MFA method among the three methods.Some hyper-scattering spots in the belly region (yellow dashed circles) are sharp and recognizable in the MFA image, whereas they are blurred in the SFA image.These hyper-scattering spots are almost unrecognizable in the single acquisition image.At a deeper depth Z3, a thin hyper-scattering tissue is noted in the MFA image [green arrowheads in Fig.4(l)] but difficult to recognize in the single acquisition and the SFA images [Figs.4(j)and (k)].This tissue is also observed at the upper depth Z1 [green arrowheads in Fig.4(d)-(f)].A similar structure has been visualized in adult zebrafish by OCT[22], and it is expected to be the outer layer of the swim bladder.Unidentified structures elongated along the dorsal side of the zebrafish [white ellipses in Fig.4(j)-(l)] are also visualized with the highest contrast in the MFA image.We also computed the SBR from Fig.4(j)-(l) for image contrast ) and (b) show representative B-scans of the phantom without and with refocusing, respectively.Note that Fig. 5(b) is identical to Fig. 3(a).

Fig. 5 .
Fig. 5. (a) and (b) are B-scans at the same location of a phantom volume without and with refocusing, respectively.(c) and (d) show representative spatial auto-correlation functions of the linear-scaled images at a depth [dashed vertical lines in (e)] along both directions without and with applying refocusing, respectively, from which the FWHMs are measured.(e) shows the FWHMs along the depth, which are considered to be proportional to the speckle size.The first and second rows present the results without and with computational refocusing, respectively.Each column presents the results for a volume measured with different focal positions.The blue and red plots show the FWHMs for the fast-scanning (X) and slow-scanning (Y) directions.The B-scans (a) and (b) were taken from the volumes corresponding to the plots highlighted by green boxes.The used volumes are identical to those used for the MFA images in Fig. 3.

Fig. 6 .
Fig. 6.SBR plot against the number of averaged volumes N and the overall focus shifting distance D. The type of plotted dot indicates the focus shifting step ∆z.The plot was produced to select optimal parameters of N, D, and ∆z.The measured sample is the phantom shown in Fig. 3, and the SBR was computed by the process described in Section 3.1.The black and green plots are projections showing the relationships between the SBR and D and between the SBR and N, respectively.

Fig. 7 .
Fig. 7. En face images of the zebrafish measured with different MFA parameters at the same depth and the same lateral position.Red boxes indicate two selected regions for computing SBR.Among the four sets of parameter configurations, configuration (c) gives best contrast.This configuration was thus selected as the optimal parameter set in the present study.

4. 3 .
Implementation of ETL 4.3.1.Advantages of ETL Other wavefront manipulation devices have recently been used to reduce MS for improving the OCT image contrast.One such device is a deformable mirror, which was used by Liu et al. in aberration-diverse OCT.

Fig. 8 .
Fig. 8. Optical design schematics of the sample arm without (a) and with (b) the ETL.(c) Color photograph of the sample arm with the ETL.GS: galvanometric scanner, OBJ:objective, ETL: electrical tunable lens, and CM: collimator.Since the ETL is a refractive optical element and easy to align, switching between the two configurations without and with the ETL takes only a few minutes.(d) shows the simulation which reveals that the focus spot size becomes smaller as the focus is shifted away from the objective.The focus position is denoted by the distance from the bottom surface of the objective to the focal plane, the blue background denotes the focus shifting range D applied in the optimal protocol, and the star mark denotes the simulated spot size without ETL.