Optical coherence refraction tomography.

Optical coherence tomography (OCT) is a cross-sectional, micron-scale imaging modality with widespread clinical application. Typical OCT systems sacrifice lateral resolution to achieve long depths of focus for bulk tissue imaging, and hence tend to have better axial than lateral resolution. Such anisotropic resolution can obscure fine ultrastructural features. Furthermore, conventional OCT suffers from refraction-induced image distortions. Here, we introduce optical coherence refraction tomography (OCRT), which extends the superior axial resolution to the lateral dimension, synthesising undistorted cross-sectional image reconstructions from multiple conventional images acquired with angular diversity. In correcting refraction-induced distortions to register the OCT images, OCRT also achieves spatially resolved refractive index imaging. We demonstrate >3-fold improvement in lateral resolution as well as speckle reduction in imaging tissue ultrastructure, consistent with histology. With further optimisation in optical designs to incorporate angular diversity into clinical instruments, OCRT could be widely applied as an enhancement over conventional OCT.

Introduction

Optical coherence tomography (OCT) is a cross-sectional, micron-scale imaging modality based on coherence gating for depth resolution[1] that has become the clinical standard of care for pathological diagnosis and treatment monitoring in several medical specialties. For long depth of focus OCT systems typical of those in clinical settings, the axial resolution is dominated by the coherence length of the light source and can be sub-micrometre[2]. On the other hand, the lateral resolution of OCT relies on confocal gating of a focused Gaussian beam, characterised axially by its depth of focus. Specifically, there is a tradeoff between a Gaussian beam’s waist diameter (which determines the OCT lateral resolution), and its depth of focus (which limits the OCT imaging depth field of view)[3]. Thus, conventional implementations of OCT commonly accept anisotropic resolutions, in particular with inferior lateral resolutions on the order of 10 μm or greater, in order to obtain long depths of focus on the order of hundreds of microns to millimetres for bulk tissue imaging. Previous studies have addressed this tradeoff between the lateral resolution and depth of focus using techniques such as digital refocusing[4,5], beam-shaping[6-12] (or both of these combined[13]), and dynamic focusing[14-16], but in all cases the lateral resolution was still limited by the diffraction limit. Furthermore, none of these methods explicitly addressed the anisotropic resolution, which can obscure ultrastructural features depending on their orientations with respect to the optical axis.

Another limitation of conventional OCT is that images are distorted due to spatially inhomogeneous refractive index (RI) distributions in tissue. Previous attempts to correct such RI-induced distortions required a priori information on both the sample geometry and RI values[17-19], which severely limits the generalisability of such methods.

To address both shortcomings of conventional OCT, we present a technique called optical coherence refraction tomography (OCRT), which uses multiple OCT cross-sectional images (‘B-scans’) acquired at a diversity of angles to reconstruct isotropic, high-resolution, cross-sectional images with the superior axial coherence gating of conventional OCT extended to the lateral dimension. In the absence of refraction, light rays travel in straight lines and the reconstruction procedure is analogous to that of X-ray computed tomography (CT)[20]. However, due to refraction and optical path delays in the sample, OCT B-scans acquired at multiple angles are distorted differently and hence need to be dewarped and registered. Thus, in correcting for sample-induced refraction, OCRT is simultaneously a novel method for estimating the spatially resolved RI distribution of the sample, which is aligned with the high-resolution reconstruction. Previously published multiangle RI tomography techniques for thick samples[21-23] accounted only for path delay (not changes in ray direction), required access to the other side of the sample, and were not applied to biological samples. Focal shift-based OCT methods for depth-resolved RI measurement typically require experimental repositioning of the focus and physical depth scanning in the sample or reference arm[24-27]. However, OCRT can be implemented using any conventional OCT system, requiring only a method to rotate the sample or an appropriate scanning protocol.

Results

OCRT and CT are Fourier synthesis techniques

As OCRT can be described as a Fourier synthesis technique from multiangle illumination, it is useful to compare it to CT. A single 1D projection in CT can be considered an image resulting from the convolution of the 2D scene with a point-spread function (PSF) that is infinitely wide in projection dimension (x) and narrow in z (Fig. 1a). Thus, in Fourier space, the product of the 2D Fourier transform of the 2D scene and the kernel transfer function is a delta function in the x direction, and of finite width in the z direction (Fourier slice theorem[28]). Hence, by taking projections at multiple directions, the entire 2D Fourier spectrum of the object can be synthesised up to the frequency cutoff defined by the resolution of the 1D detector and the beam size. The inverse 2D Fourier transform then reconstructs the image with isotropic high-resolution. In practice, the backprojection algorithm can be employed, in which each 1D projection is uniformly smeared in the direction of the X-rays, and summed across all angles, which is mathematically equivalent to the Fourier synthesis technique in a continuous framework. Filtered backprojection (FBP) applies a filter to the 1D projections to account for the fact that the centre is overrepresented, the ideal filter being the ramp filter, H(f) = ∣f∣, which arises from the Jacobian of the transformation between Cartesian and polar coordinates[29].

Figure 1.

OCRT, like CT, is a Fourier synthesis technique.

a, CT and OCRT both exhibit anisotropic horizontal and vertical point-spread functions (PSFs) and transfer functions (TFs). By adding contributions from multiple angles, a TF with larger Fourier support is synthesised. In the fourth column, filtering is applied before summation to deemphasise low Fourier components. Both techniques have the same theoretical synthesised PSF, assuming common vertical bandwidths and sufficient angular coverage. Because OCRT also has horizontal Fourier support, coarser angular spacing suffices, unlike for CT, as apparent in the fourth column. b, Gaps in the TF decrease with angular spacing. Red circles denote the minimum radial Fourier support; black circles denote the maximum support. ‘Min support’, defined as the ratio of the red to black radius, approaches 1 as angular spacing decreases. The assumed OCT axial-lateral PSF ratio is 0.30. c, The worse the lateral relative to the axial resolution, the finer the required angular spacing for a particular ‘min support’. The dark red point corresponds to our experimental settings (spacing = 6°, axial-lateral ratio ~0.3, minimum Fourier support ~0.986). d, Lateral resolution increases with the max illumination angle θmax, approaching the axial resolution as θmax approaches full coverage (±90°).

OCRT and CT are conceptually similar (Fig. 1a), except that in OCRT the PSF in the projection dimension (x), defined by the lateral resolution, is not infinite but still wider than the PSF in the propagation dimension (z). Note that the projection dimension here is not necessarily the propagation dimension, but rather in general it is the dimension of the lowest resolution (in this case, the lateral dimension). Thus, if we assume that the lateral resolution is approximately constant with depth and that the wavefront curvature away from the focus of the Gaussian beam is negligible, the image is the convolution between the scene and an anisotropic PSF. In Fourier space, the transfer function corresponding to the PSF is anisotropic, with high frequency information diminished in one dimension. Analogously to CT, by taking B-scans at multiple angles, the missing high frequency information can be filled in. However, since OCT already has some lateral frequency support, OCRT requires data from fewer angles than CT. The angular spacing (assumed equal) depends upon the axial-to-lateral resolution ratio and the amount of tolerance on the gaps in the Fourier space at intermediate angles (Fig. 1b-c). As full angular coverage at ±90° is approached, isotropic resolution in the image plane, xz, is achieved (Fig. 1d).

Refraction correction

While RIs of biological tissue in the X-ray regime are close to 1, in the optical regime a spatially varying RI distribution distorts OCT images by changing the path lengths and the directions of the rays, preventing rigid-body registration among the multiangle B-scans. If the RI distribution were known, the rays could be propagated accordingly to dewarp the image prior to FBP application to generate a high-resolution reconstruction. To infer the unknown RI distribution, in OCRT we solve the inverse problem using the ray equation as the forward model, which in 2D (x and z) is where n(x, Z) is the RI distribution parameterised by and s is the position along the 1D ray trajectory[30]. The parameterisation we chose was a sum of a regularly spaced grid of Gaussian kernels such that n(x, z) is differentiable everywhere and minimises the effects of the ‘staircase’ artefacts stemming from discretisation onto a Cartesian grid. In particular, we chose the index distribution given by the Nadaraya-Watson kernel parameterisation[31]. This parameterisation also avoids having to use finite differences to compute spatial gradients, as analytic expressions are available.

Because in general there is no closed solution to the ray equation, we employed numerical integration (fourth-order Runge-Kutta[32]). To specify optical path delay, we used a constant step size for the numerical solver, which was scaled by the inverse of the RI value at the current position of the ray. To initiate ray propagation, we specified the initial conditions, which were the initial ray positions and directions. Assuming telecentric scanning and uniform sampling of A-scans within each B-scan, the initial directions were all parallel, and the initial positions were equally spaced for each angle. In this way, we propagated the A-scans from each B-scan from each angle as individual rays, which caused the overall B-scan image to dewarp. In theory, if the optimal RI distribution is found, the images would be perfectly registered. We assumed that the group and phase indices are approximately the same in biological tissue, such that a common index distribution governs both path trajectory and delay. Thus, in this paper we use the terms group (refractive) index and RI interchangeably.

Optimisation: joint registration, isotropic, high-resolution reconstruction, and RI map estimation

To provide feedback on the accuracy of n(x, z) to aid its optimisation, we required a differentiable metric that quantifies the degree of joint registration among all the B-scans. Here we describe an intensity-based metric, by which we computed the mean squared error (MSE) between the raw B-scans and a forward prediction of the B-scans based on the estimated high-resolution reconstruction (Fig. 2). The forward model started with the current estimate of the high-resolution reconstruction based on FBP along the ray trajectories governed by the current estimate of the RI distribution. For each B-scan orientation, that current estimate was then rotated, warped according to the same trajectories, blurred by the OCT PSF, depth-wise attenuated, and locally intensity-rescaled according to the orientation of the structure to the illumination. The result constituted the forward prediction of the B-scan, the MSE between which and the raw B-scan data was to be minimised with respect to the forward model parameters (e.g., the RI distribution, attenuation parameters, etc.). We also included regularisation on the index distribution to impose spatial smoothness, stabilising the solution and optimisation. After registration, we optimised the filter for FBP. More details about the forward model are provided in the supplementary Sec. 1.

Figure 2.

Overview of the iterative OCRT reconstruction algorithm.

Multiangle B-scans served as the input, which were backprojected along trajectories derived from ray propagation (ray prop) according to the current estimate of the RI distribution to form an estimate of the high-resolution reconstruction. The same trajectories sampled the current high-resolution estimate to generate forward predictions of the B-scans. The MSE (mean square error) between the B-scan data and B-scan estimates was iteratively minimised with respect to the RI distribution. For a more detailed description, see supplementary information Sec. 1.

We used a modified version of gradient descent called Adam[33] to minimise this regularised MSE, which jointly registered the B-scans and estimated both the RI distribution and the high-resolution, isotropic image. We computed the gradient of the intensity-based registration metric through the numerical differential equation solver with respect to the forward model parameters using TensorFlow[34], a software library that employs automatic differentiation, a widely used technique in the deep learning community requiring only the specification of the differentiable forward model. Based on this forward model, the gradient is computed through recursive application of the chain rule.

In summary, the OCRT optimisation procedure simultaneously registers all B-scans, generates an undistorted, isotropic, high-resolution reconstruction, and a co-aligned estimate of the RI distribution of the sample. An overview of this intensity-based implementation of OCRT is given in Fig. 2. A more detailed description of this implementation can be found in the supplementary Sec. 1.

Experimental results

To validate the isotropic resolution of OCRT, we imaged 560-nm polystyrene beads embedded in a 2% (w/v) agarose gel inside a glass microcapillary tube (Fig. 3). We used a custom-designed inverted rotation stage (supplementary Fig. S4) to rotate the tube sample, which was immersed in water to reduce the RI contrast with the glass tube wall. Because the index contrast was still large, our optimisation procedure started with gross estimation of the index distribution by first assuming a circular tube geometry and optimising with respect to only the index of the medium and the glass, with the tube diameter, wall thickness, and centre calibrated in a separate step (see supplementary Secs. 1 and 2). The subsequent optimisation step refined this gross estimate by allowing the entire 2D index distribution to vary arbitrarily, as well as small translations of the OCT B-scans. Excision of tissue samples and their insertion into cylindrical tubes were used in this first report to aid in sample mounting and rotation, however these steps are not in general required for OCRT (see Discussion).

Figure 3.

Experimental validation of the isotropic resolution and RI estimates of OCRT.

a, 1,200-averaged OCT B-scan (also averaged across 18 frames in the out-of-plane dimension), histogram-matched to b, the OCRT reconstruction. The two sets of close-up views (red and blue boxes; dashed, averaged B-scan; solid, OCRT reconstruction) show clear lateral resolution enhancement. The arrows in the zoomed-in views indicate beads resolved in the OCRT reconstruction, but not in the averaged B-scan. c, Bar plots of the 2D Gaussian width (FWHM) fits (median, with 90% data coverage intervals after filtering out poor fits; numerical labels are medians) demonstrate a >3-fold enhancement in lateral resolution, and that the spatial resolution of OCRT is isotropic. For OCRT, the x and z axes were defined arbitrarily. d, The left column shows the OCRT RI map for the polystyrene bead sample embedded in a 2% agarose gel reconstructed in b. The right column shows the RI map for a separate polystyrene bead phantom, which was instead embedded in PDMS, which has a higher RI. The estimated RI values for the embedding media fall within the uncertainty range (dotted lines represent one standard deviation) estimated using OCT (see methods section). Scale bars, 100 μm.

To quantify the resolution of OCRT, we localised the beads of both the reconstruction (Fig. 3b) and the averaged volume (Fig. 3a), and fit a 2D Gaussian function with axes oriented with the Cartesian coordinate system, where z and x are the axial and in-plane lateral dimensions, respectively (for the reconstruction, the x and z axes are arbitrary). The resulting distributions of the PSF width fits in x and z (Fig. 3c) demonstrate a factor of >3 enhancement in the lateral resolution, and confirm that the resolution afforded by OCRT is isotropic, given by the axial resolution of the original OCT images.

To validate our RI distribution estimates, we created another bead phantom substituting the agarose gel for polydimethylsiloxane (PDMS). Fig. 3d shows that the OCRT index estimate was in good agreement for both agarose and PDMS, falling within one standard deviation of the bulk RI estimated independently from OCT pathlength measurements (see Methods).

To demonstrate the performance of our method, we imaged and applied OCRT to a variety of ex vivo biological samples, including several mouse organs (vas deferens, femoral artery, bladder, and trachea), a human donor cornea, and a marsh crane fly (Tipula oleracea) leg. All samples were inserted into microcapillary tubes for convenience during imaging, and haematoxylin and eosin (H&E)-stained histological sections of neighbouring tissue samples were obtained for comparison. Masson’s trichrome stain was additionally obtained for the artery sample. Overall, we observed that lateral resolution was significantly improved in the OCRT reconstruction across all samples imaged, and that the features in the RI maps matched features in the OCRT reconstructions. In the first vas deferens sample (Fig. 4a-h), the structures of the adventitia, smooth muscle, lamina propria, and the transitional epithelium are clearly resolved. For example, it is apparent in the OCRT reconstruction but not in the B-scan (Fig. 4a,b) that the smooth muscle layer consists of more circumferential structure, in contrast to the adventitia, which has a different texture; this is consistent with the histology (Fig. 4d). Furthermore, the lamina propria is very prominent in the OCRT reconstruction as a darker, more strongly scattering layer, unlike in the histology where the delineation is more subtle, and the clover-leaf-shaped transitional epithelium is much more recognisable in the OCRT reconstruction than in the B-scan. In the second vas deferens sample (Fig. 4i-p), in addition to the features in the first vas deferens sample, we observe clear evidence of resolution enhancement in the artefactual ribbon-like detachments of the adventitia, consistent with the histology (Fig. 4l).

Figure 4.

OCRT of mouse vas deferens and femoral artery.

a, 1,200-averaged B-scan of mouse vas deferens, histogram-matched to b, the OCRT reconstruction. c, RI map that registers the multiangle B-scans. d, Histology of vas deferens from a separate animal. e-h, Zoomed-in regions of interest in a and b. i-p shows the same information for another vas deferens sample from a different animal. In both samples, structures in the adventitia (Ad) are better resolved in the OCRT reconstruction, particularly where the adventitia is artefactually detached (d). The adventitia has a different texture than the smooth muscle layer (SM) in b, consistent with histology. The transitional epithelium (TE), lumen, and lamina propria (LP) are more apparent in the OCRT reconstruction. q-x panels show the same information for a collapsed mouse femoral artery. Panel t shows two adjacent histological slices of femoral artery from a separate animal—the top slice was stained with H&E, the bottom slice with Masson’s trichrome stain to accentuate the internal and external elastic lamina (IEL and EEL) lining the SM. Unlike in the B-scan, in the OCRT reconstruction, the IEL and EEL are clearly resolved, the former exhibiting undulatory patterns, consistent with histology. Scale bars, 100 μm.

In the mouse femoral artery sample (Fig. 4q-x), the smooth muscle layer and the external and internal elastic lamina are much more apparent in the OCRT reconstruction, while essentially invisible in the B-scan. Note in the OCRT reconstruction the undulatory pattern of the internal elastic lamina (Fig. 4v,x), consistent with the Masson’s trichrome-stained histology, which accentuates the elastic lamina as with a blue colour (Fig. 4t, bottom slice). The internal elastic lamina is lined by nuclei in the OCRT reconstruction, appearing as dark dots (Fig. 4x), which is consistent with the H&E-stained histological slide (Fig. 4t, top slice).

In the mouse bladder sample (Fig. 5a-h), the layered structures of the smooth muscle are much clearer in the OCRT reconstruction than the B-scan, where they are less apparent due to lower lateral resolution and presence of speckle noise. The connective tissue in the lamina propria is also better resolved in the OCRT reconstruction. In particular, comparison of Fig. 5e and 5g to 5f and 5h, respectively, reveals connective tissue features of the lamina propria present in the OCRT reconstruction, but not in the B-scan. In the mouse trachea sample (Fig. 5i-p), the basement membrane and ciliated epithelium are well resolved across the entire sample as dark lines, while only partially visible in the B-scan where they are near normal to the optical axis, suggesting that the anisotropic resolution of the B-scan is the culprit. The round cartilaginous structures are also clearer in the OCRT reconstruction.

Figure 5.

OCRT of mouse bladder and trachea.

a,1,200-averaged B-scan (inverted) of a mouse bladder sample, histogram-matched to b, the corresponding high-resolution OCRT reconstruction, which shows the layers of the smooth muscle (SM) more clearly than in the B-scan. Furthermore, connective tissue structures in the lamina propria (LP) are much better resolved. c, RI map that registers the multiangle B-scans, which exhibits two distinct layers, corresponding to the smooth muscle (SM) layer and the transitional epithelium (TE). d, H&E-stained histological section from the same animal. e-h, close-up views of regions of interest of the OCT B-scan and OCRT reconstruction denoted in a and b. i-p shows the same images except for mouse trachea (B-scan also inverted). The histology in l was obtained from a separate animal. The basement membrane (BM) and ciliated epithelium (CE) are well resolved throughout the entire sample, while partially for the averaged B-scan. The cartilaginous (C) structures are also clearer in the OCRT reconstruction. The adventitia (Ad) is delineated by another dark line that is not readily apparent in the B-scan. Scale bars, 100 μm.

The human donor cornea sample is shown in Fig. 6a-h. The lamellar structures of the corneal stroma are parallel to the front and back surfaces of the cornea. Thus, the features of the stroma are largely axial and hence the increase in resolution for OCRT over the B-scan is not immediately appreciable, but rather the visibility of the lamellae increases due to speckle reduction (see also Fig. S5). More interestingly, the RI maps generated by OCRT range from 1.38-1.40 (Fig. S6), which agrees with previous bulk group index measurements at similar near-infrared wavelengths, considering the influence of dehydration on the index[35-37]. We obtained similar results on a different human cornea sample (Supplementary sec. 3), giving additional validation to the RI-mapping capability of OCRT.

Figure 6.

OCRT of human cornea and crane fly leg reveals additional RI information.

a, 1,200-averaged B-scan of human cornea, histogram-matched to b, the high-resolution OCRT reconstruction. d, Histology from the same sample. e-h show zoomed-in regions of interest in a and b. The lamellar structures of the stroma (St) are more clearly defined in the OCRT reconstruction and consistent with histology. The epithelium (Ep) and Descemet’s membrane (DM) are also apparent, although partially damaged from processing (Ep damage is consistent with histology). The RI map in c, further analysed in Fig. S6, agrees with literature values[35-37] (RI ~1.38-1.40). i-l shows information analogous to that of a-d for crane fly femur (Fem) and tibia (Tib). The outer wall of the femur (cuticle; Cu) is apparent in the OCRT reconstruction, consistent with histology (j). Small hair-like protrusions of the cuticle (setae; Se) are apparent in the OCRT reconstruction (l). The circumferential orientation of the longitudinal muscle fibres (LMF) are clearer in the OCRT reconstruction than B-scan. The RI map (i) shows an index close to 1 inside the tibia, while close to that of water inside the femur, which is not immediately apparent in the B-scan (k). Scale bars, 100 μm.

Finally, the OCRT reconstructions of the tibia and femur cross sections of a crane fly (Fig. 6i-l) exhibit superior resolution of the outer walls (cuticle) and the ultrafine features that surround them. The cuticle, apparent as a thin orange structure in the H&E, manifests at the bottom of femur as two thin dark lines in the OCRT reconstruction. Note also that although the larger hair-like structures (setae) present in the B-scan and OCRT reconstruction are not visible in the H&E, likely lost during tissue processing, the much finer setae lining the wall in the H&E are still visible in the OCRT reconstruction. The circumferential orientation of the internal longitudinal muscle fibre is also much more apparent in the OCRT reconstruction than in the B-scan. Furthermore, the RI map indicates that the lumen of the tibia was filled with air (n ≈ 1) while the lumen of the femur was filled with water. Although this is disparity is likely due to the sample preparation procedure, it is a feature not at first apparent from the B-scan. However, upon closer inspection, we note that the region indicated by the green arrow in Fig. 6k is displaced upward relative to the regions in the neighbouring A-scans that do not include the tibia lumen, which indicates that the tibia lumen has a lower RI (i.e., air). Furthermore, we note that the tissue-to-lumen interface is more intense in the tibia than in the femur, indicating a larger index contrast. Thus, the RI mapping capabilities of OCRT can reveal additional information about the sample and simultaneously dewarp the image.

Discussion

We introduce optical coherence refraction tomography (OCRT), a method that uses OCT B-scans from multiple angles to simultaneously reconstruct an undistorted, isotropic, high-resolution image and estimate the RI distribution. In particular, the superior axial resolution of OCT is extended to the inferior lateral resolution to form an isotropic PSF. We demonstrate resolution enhancement and RI estimates of OCRT not only in sub-resolution scattering bead phantoms, but more importantly in several biological samples.

OCRT differs from optical projection tomography[38], which although similar to CT uses projections along rather than perpendicular to the optical axis, hence uses weakly focused beams and suffers from the same depth-of-focus tradeoff. Our method also differs from rotational imaging OCT[39], which uses multiple angles to overcome only the limited depth of penetration of OCT.

OCRT is advantageous over other techniques that address the depth-of-focus tradeoff[4-16] because it completely replaces the lateral resolution dependence on diffraction with the axial resolution dependence on coherence. For example, recent work in applying OCT to the extreme ultraviolet regime has yielded axial resolutions on the order of 10s of nm, but with lateral resolutions three orders of magnitude lower[40]. We also note that unlike these other techniques OCRT by design renders the resolution isotropic. Furthermore, since OCRT does not require high-NA objectives, much longer working distances associated with low NAs are possible, which make in vivo imaging more practical. Moreover, OCRT can tolerate lateral aberrations, in particular those induced by the sample. While the resolution of OCRT can be degraded by chromatic dispersion induced by the sample, in principle they can be corrected computationally. Another advantage of OCRT over other inverse optimisation approaches such as interferometric synthetic aperture microscopy (ISAM) is that OCRT does not require phase stability[4,13], which for near infrared wavelengths typical of OCT requires nanometre-scale stability of both the scanner and the sample over the acquisition time. On the other hand, OCRT, which employs intensity images, only requires stability on the order of the axial resolution, and even this constraint may in principle be relaxed through computational corrections during the registration step. Another advantage is improved image contrast from the speckle reduction due to angular compounding of independent speckle patterns[41-43]. Finally, OCRT does not sacrifice signal-to-noise ratio (SNR), unlike Bessel beam-based OCT[6,7,11-13] because a substantial fraction of the power in Bessel beams is contained in the side lobes, and unlike ISAM[4,13] which suffers SNR fall off with distance from the nominal focus.

However, as a technique that relies on nonlinear, nonconvex optimisation, OCRT may be sensitive to initialisation. For completeness, we also note three limitations stemming from the fact that this present study is a proof of concept of OCRT. First, we only varied the angles in the 2D plane, such that there was no out-of-plane (y) resolution enhancement. However, in theory this method could be extended to 3D by taking volumes at orientations described by two angles in spherical coordinates[44]. Second, all samples in this study were subject to destructive preparation using glass capillary tubes to assist in full 360° data collection, which led to some sample artefacts (i.e., adventitial separation in the vas deferens and disruption of the corneal epithelium and endothelium). Both of these limitations could be solved by replacing sample rotation with an angularly scanning probe, as pioneered for metrology applications[44]. Alternatively, we have previously described OCT system designs in which dual scanners placed in planes conjugate and anti-conjugate to the sample under study allow independent control of the beam’s entry position and angle to the sample[45]. These approaches would make OCRT extensible to in vivo medical imaging of exposed biological surfaces such as the cornea, ocular anterior chamber, skin, or luminal surfaces. Such imaging conditions would place constraints on the range of angles available, and thus future work will include extending OCRT to limited angle tomography[46], which we discuss below. The third limitation is the relatively long acquisition times stemming from B-scan averaging (each of the 60 angular B-scans was 20 times averaged). However, we show in supplementary Sec. 5 that 20 times averaging is conservative, specifically that even with no averaging OCRT still obtains substantial improvement over conventional OCT. Nevertheless, there are numerous avenues to speed up acquisition to make OCRT compatible with in vivo imaging. Although our OCT system operated at a 20-kHz A-scan rate, multi-MHz-rate systems have been reported[47] that would acquire the same number of raw B-scans (1,200) used even in our highly averaged OCRT reconstructions in <0.5 seconds. We also note that OCRT includes a registration step such that even at slower acquisition rates OCRT can in principle correct for motion. Moreover, OCRT can take advantage of advanced denoising methods, which have been demonstrated to mitigate the need for averaging[48]. Another approach to reduce the number of averages and therefore acquisition time is to use higher incident powers to increase the SNR, subject to appropriate laser exposure safety limitations. Finally, the number of angles from which images are acquired could be reduced, since according to Fig. 1c, we oversampled the angles required. Thus, future work will explore limited angle tomography and other ray subsampling strategies, in analogy with compressive sensing in CT[49]. This is analogous to the ‘missing cone’ problem[46,50-54] in diffraction tomography, by which strictly along the axial dimension there is zero Fourier support. Nevertheless, various strategies have been employed to fill in this missing cone using regularisation[46,52], which could readily be adapted to limited angle OCRT.

In summary, OCRT is a general framework that leverages multiangle OCT images to computationally synthesise both an enhanced-resolution reconstruction and an RI map, which is not available in conventional OCT. We demonstrated substantial improvements of OCRT over OCT in various biological samples, and found that there were structures readily apparent in OCRT reconstructions but missing in conventional OCT images. Finally, there are readily achievable multiangle scanning strategies and extensive prior work on limited angle tomography in other domains such as CT and diffraction tomography that would realistically propel OCRT towards in vivo medical imaging applications.

Methods

OCT imaging

All OCT images were acquired using a commercial spectral domain OCT system (Bioptigen Envisu R4110 XHR SDOIS) at a 20-kHz A-scan rate, with an incident power of 1 mW, an 820-nm centre wavelength, and a nominal axial resolution of 1.2 μm in air and lateral resolution of 8.5 μm. Our measured resolution using sub-resolution beads is given in Fig. 3a. The B-scans consisted of 500 A-scans over a field of view of 1.5 mm; each A-scan had 2,048 pixels with a maximum imaging depth of 2.22 mm in air.

For each sample, we acquired 20-averaged OCT B-scans or volumes every 6° across 360°. For comparison, we acquired, registered using a previously developed subpixel registration algorithm[55], and averaged 1,200 B-scans or volumes at one angle (Fig. 3a). For the polystyrene bead phantoms only, we performed out-of-plane averaging of B-scans across 18 frames spaced 3 μm apart to account for slight tube misalignment, due to which the beads close to the wall of the tube away from the axis of rotation, because of their small size, were not present in every B-scan. We use these y-averaged B-scans as the raw data for OCRT reconstruction. For all biological samples this was not an issue and so we did not average B-scans across the y direction.

Preparation of samples

All samples were inserted into capillary microtubes (Drummond Scientific Microcaps) with an inner diameter of 797.6 μm and immersed in phosphate buffered saline (PBS). 560-nm polystyrene beads (Thermo Scientific) were embedded in 2% agarose gel (Sigma Aldrich) or PDMS (Dow Corning).

The animal organs used in this study were from mice (C57BL/6 wild-type) euthanised and discarded from unrelated experiments and thus were not subject to Institutional Animal Care and Use Committee review. The organs were fixed in 10% neutral buffered formalin (NBF) and kept at 4°C for 24 hours. After fixation, all organs were micro-dissected using an upright dissection scope (Zeiss). The dissected tissues were inserted into the glass microtubes using dissection forceps. In some instances, the tissues were placed into the glass tube following aspiration with a 23G truncated needle on the opposite side of the tube. This procedure allowed the specimen to gently slide into the inner wall of the glass microtubule. To avoid dehydration after fixation, all specimens were transferred with the glass tube in PBS. The remaining organ tissue was used for paraffin block preparation followed by H&E staining or trichrome staining for histological examination. Histological images were acquired at 10× or 20× magnification using an Olympus microscope and digitally white balanced.

The donor human cornea samples were obtained from Miracles in Sight and were declared exempt from review by the Duke University Health System Institutional Review Board. The samples were stored in Optisol at 4°C until it was transferred to and fixed in 10% NBF for 24 hours. Cornea samples underwent the same procedure for tube embedding and histology as for the mouse tissue.

Estimating the ground truth RI of uniform embedding media using OCT

We randomly sampled n = 8 tubes and estimated their inner diameters using OCT by measuring the difference between the luminal reflections, and obtained 797.3 ± 0.9 μm (standard error of the mean), in very good agreement with the manufacturer’s specification of 797.6 μm. We then estimated the inner diameter pathlength of the bead samples (with either agarose or PDMS embedding media) using the same method for all angles and volume (y) slices, across which the standard deviation values in Fig. 3d were calculated.

Sample stage

To acquire images at multiple angles, we used a custom-built, inverted rotation stage (Thorlabs) designed to mount a single tube vertically. The tube has two orientational degrees of freedom and two translational degrees of freedom perpendicular to the axis of rotation, which we use to manually align the tube. The tube is immersed in water or PBS in a cuvette, which remains stationary during sample rotation. The entire above setup is mounted on xyz-translation stage (Newport). See supplementary Fig. S4 for a schematic layout. Sample rotation and data acquisition were automated using a custom Python script. For all experiments, we acquired B-scans at 60 orientations, spaced evenly across 360°.

Numerical optimisation

Optimisation was conducted on TensorFlow 1.8 using Python 2.7 on an Intel Core i7-3930K with 48 GB of RAM or in Google Cloud Platform. Because each iteration required 10s of GB of memory, we used the CPU version of TensorFlow. For all samples we ran gradient descent for 200-500 iterations, on the order of a 2-3 minutes per iteration. Filter optimisation was run for an additional 100 iterations. We used the same regularisation hyperparameters for all biological samples.