Multifocal interferometric synthetic aperture microscopy

Yang Xu; Xiong Kai Benjamin Chng; Steven G. Adie; Stephen A. Boppart; P. Scott Carney

doi:10.1364/OE.22.016606

1. Introduction

Optical coherence tomography (OCT), is an optical imaging method that reveals the 3D structure of semitransparent samples with micron-resolution [1–6]. OCT has proved useful in many fields including ophthalmology [7–12], dermatology [13–15], and cardiology [16], among many others. We present a method of 3-D optical imaging that combines multifocal OCT and interferometeric synthetic aperture microscopy (ISAM).

In OCT, two of the more important figures of merit are DOF and transverse resolution. However, using Gaussian beam optics, it is impossible to improve both at the same time. For a given OCT system, the transverse resolution at each depth is equivalent to the width of the scanning beam at that depth. The transverse resolution is best on the focal plane, where the scanning beam has minimum beam width. For given free space wavelength λ₀, free space numerical aperture NA₀ and media refractive index n, the beam waist radius is w₀ = λ₀/(nπNA₀). The transverse resolution degrades as the position of reconstruction is further away from focal plane, where the scanning beam diverges in width. For a Gaussian beam propagating in the z direction, with focal plane at z₀, the beam radius can be expressed as

w (z - z_{0}) = \sqrt{w_{0}^{2} + {[\frac{λ_{0} (z - z_{0})}{π n w_{0}}]}^{2}} .

Because of the beam-spreading, the reconstructed image becomes defocused at positions far away from the focal plane. Therefore, the DOF of the OCT system is usually taken to be within a Rayleigh range (Z_R) of the focal plane on both sides. The Rayleigh range is defined as the distance from the focal plane to the plane where

w (z) = \sqrt{2} w_{0}

. It can be expressed as

Z_{R} = n π w_{0}^{2} / λ_{0},

which is proportional to the square of the transverse resolution.

Thus, the DOF is a function of resolution as shown in Fig. 1, where the wavelength is taken to be 1 μm and the sample refractive index is 1.4. For example, if the desired transverse resolution is 0.76 μm on the above mentioned condition, which can be achieved by using a 0.3 NA scanning beam, the achievable DOF is 5.1 μm, as can be seen at the lower sample point in Fig. 1. On the other hand, if the DOF is increased to 182 μm, which can be implemented using a 0.05 NA scanning beam, the best possible transverse resolution degrades to 4.55 μm, as can be seen at the upper sample point.

Fig. 1 Trade-off between the (effective) DOF and transverse resolution in OCT, multifocal OCT and ISAM, at wavelength of 1 μm and with refractive index of 1.4.

Download Full Size | PDF

Although OCT is limited by the above-mentioned trade-off, it is possible to mitigate the trade-off by improving either the physical setup or the reconstruction algorithm. Recent solutions include Bessel beam OCT, interferometric synthetic aperture microscopy (ISAM), and multifocal OCT, each have their own advantages and challenges.

Bessel beam OCT [17,18] uses Bessel beam for illuminating the sample to achieve all-optical DOF extension. The Bessel beam is generated using an axicon lens. The scattered signal is collected using either the axicon lens or a conventional objective lens. This method trades for larger DOF at the cost of higher side lobes.

Multifocal OCT, also called multibeam OCT or multichannel OCT [19, 20], is a version of OCT where multiple apertures are used simultaneously in scanning. Each aperture is focused at a different depth during the scan. One set of OCT data is collected from each aperture. An image stitching method is used to combine each set of data, where the part with the best transverse resolution from each set of data are stitched into the final reconstruction. Compared with each set of OCT data, the multifocal OCT image shows an improved resolution over the whole depth, or from another point of view, enables a larger DOF. Consider the case of a multifocal OCT system with 3 identical apertures. The distance between two adjacent focal planes is taken to be 2 Rayleigh ranges. The blue curve in Fig. 1 plots the trade-off curve of multifocal OCT. As three stacked apertures are used, for any value of transverse resolution, multifocal OCT provides 3 times the DOF compared with OCT; or on the other hand, for the same desired DOF, multifocal OCT can achieve improved transverse resolution compared to OCT. Even more apertures can be used at the cost of additional setup complexity. Multifocal OCT also faces challenges such as intensity discontinuities along stitching edges in the reconstructed image [19, 20]. Although artificial blending can be used to mitigate this effect [21], it is not physically derived and thus does not always provide accurate image reconstruction.

The ISAM experiment is similar to that of OCT. ISAM takes the OCT data and corrects for the spreading of the beam using the solution to the inverse problem of OCT [22–27]. The algorithm can be efficiently implemented using a frequency domain resampling technique [25, 28], which can be expressed by

〈 r | η 〉 = w (z - z_{0}) \int d^{2} k_{‖} 〈 r_{‖} | k_{‖} 〉 \int d β 〈 z | - β 〉 H^{- 1} 〈 k_{‖}, \frac{1}{2} \sqrt{β^{2} + k_{‖}^{2}} | S 〉,

where r is the space coordinate; k is the wavevector; k_‖ and r_‖ are the vector projection onto the plane normal to the z axis, i.e. the beam axis; and β is the longitudinal spatial frequency. The expression 〈k_‖, k₀|S〉 is the 2-D Fourier transform of 〈r_‖, k₀|S〉. 〈r_‖, k₀|S〉 represents the OCT data |S〉 collected at transverse position r_‖ and wavenumber k₀. In the equation, k₀ is replaced with

\frac{1}{2} \sqrt{β^{2} + k_{‖}^{2}}

, which represents a set of calculated grid points to be resampled onto each β. The inner product 〈r_‖|k_‖〉 = (2π)⁻¹exp(−ik_‖ · r_‖), and 〈z| − β〉 = (2π)^−1/2exp(izβ). 〈r|η〉 is the sample susceptibility at position r. The expression w(z − z₀) is the beam radius at depth z with the focal plane at z = z₀, as is defined in Eq. (1). The factor H depends on k_‖, k₀, and position of the focal plane z_n of the n^th beam focus. This equation is similar to Eq. (19.32) in [28] except that the expression |z − z₀| in [28] are replaced by w(z − z₀), which generalizes the asymptotic results at both near focus and far from focus, and provides a smooth transition in between.

Because the DOF of OCT is infinitely extended by applying the algorithm, ISAM is able to provide spatially invariant resolution, as opposed to the decaying transverse resolution of OCT along depth axis. ISAM has been validated and is proved to provide accurate reconstruction with focal-plane resolution [24, 29].

However, ISAM faces the challenge of signal-to-noise ratio (SNR) decreasing with 1/z as does for standard OCT, which limits its effective DOF [24]. Assuming the SNR at the focal plane is SNR₀ (linear), and the lowest tolerable SNR is SNR_L, the effective DOF can be expressed as

Effective DOF = \frac{2 π n ω_{0}^{2}}{λ_{0}} \sqrt{{(\frac{{SNR}_{0}}{{SNR}_{L}})}^{2} - 1} \approx \frac{2 π n ω_{0}^{2}}{λ_{0}} \frac{{SNR}_{0}}{{SNR}_{L}} .

The red curve in Fig. 1 plots the trade-off curve of ISAM when SNR₀/SNR_L = 10 as an example. The curve would shift upward or downward linearly with the ratio of SNR₀/SNR_L for SNR₀/SNR_L ≫ 1. Compared with multifocal OCT, ISAM further shifts the curve to the upper-right quadrant of the plot, mitigating the trade-off.

2. Multifocal interferometric synthetic aperture microscopy

Seeing the trade-off OCT is facing and the advantages and challenges of ISAM and multifocal OCT, we have developed a new imaging technique which combines the advantages of both multifocal OCT and ISAM. By applying a coherent combination of the reconstruction of each channel with appropriate regularization, higher lateral resolution and signal-to-noise ratio can be achieved without artifacts from the synthesis of the data sets. The experiment is similar to conventional OCT, except multiple channels, each with a different focus, are used.

2.1. Forward model

The sample is again described by a susceptibility |η〉, and the data set by |S〉, but |S〉 is now in a larger Hilbert space. The position of the n^th focal spot is denoted r_n. Transverse wavevector states are denoted |k_‖〉 so that $〈 r_{‖} | k_{‖} 〉 = \frac{1}{2 π} e^{- i k_{‖} \cdot r_{‖}}$ . We start from Eq. (19.26) in [28], which describes the multifocal OCT signal component acquired with illumination focal spot at r_m and detection focal spot at r_n with wavenumber k₀, that is

\begin{array}{r} 〈 r_{m}, r_{n}, k_{0} | S 〉 = \int d^{2} k_{‖}^{″} \int d^{2} k_{‖}^{'} \int d^{2} r_{‖} d z e^{i k_{‖}^{″} \cdot r_{n}} {\tilde{g}}_{n} (k^{″}, k_{0}) 〈 k_{‖}^{″} | G_{z_{n}} | r 〉 \\ \times 〈 r_{‖} | K_{z_{m} - z} | k_{‖}^{'} 〉 e^{- i k_{‖}^{'} \cdot r_{m}} {\tilde{g}}_{m} (k^{'}, k_{0}) 〈 k_{0} | U_{0} 〉 〈 r | η 〉, \end{array}

where again 〈r|η〉 is the sample susceptibility at position r; 〈k₀|U₀〉 is the field amplitude as a function of wavenumber;

{\tilde{g}}_{m, n} (k^{″}, k_{0}) = exp (- π^{2} k_{‖}^{2} / 2 k_{0}^{2} {NA}_{m, n}^{2})

is the Fourier transform of the beam waist of the illumination and the detection aperture; K and G are the field propagator and Green’s function operator respectively.

We then assume the transverse-scanning position of the illumination aperture is infinitely close to that of the detection aperture and have the same NA. Therefore, we can set r_‖m = r_‖n = r_‖, g̃_m = g̃_n = g̃, so that the problem is converted to the single static case with transverse-scanning position r_‖, illumination focal plane at depth of z_m and detection focal plane at z_n. The datum at transverse-scanning position r_‖, wavenumber k₀, illumination focal plane z_m and detection focal plane z_n is denoted 〈r_‖, k₀, z_m, z_n|S〉. By performing a change of basis and expressing the received signal as the superposition of the contribution from each illumination aperture, the multi-foci version of Eq. (19.29) in [28] can be derived. The signal received at aperture n can be expressed as

\sum_{m = 1}^{M} 〈 k_{‖}, k_{0}, z_{m}, z_{n} | S 〉 = 2 π 〈 k_{0} | U_{0} 〉 \int d z \tilde{h} (k_{‖}, z, k_{0}, z_{m}, z_{n}) 〈 k_{‖}, z | η 〉,

where 〈k_‖, z|η〉 is the Fourier transform of 〈r|η〉 with respect to the transverse coordinates. Here

\tilde{h} (k_{‖}, z, k_{0}, z_{m}, z_{n}) = \int d^{2} k_{‖}^{″} \frac{\tilde{g} (k^{″}, k_{0}) \tilde{g} (k_{‖}^{″} - k_{‖}, k_{0})}{\sqrt{k_{0}^{2} - k_{‖}^{″ 2}}} \sum_{m = 1}^{M} e^{i \sqrt{k_{0}^{2} - k_{‖}^{″ 2}} (z_{n} - z) + i \sqrt{k_{0}^{2} - {(k_{‖}^{″} - k_{‖})}^{2}} (z_{m} - z)}

, which is the two dimensional Fourier transform of the point spread function of the system, which consists of M stacked components, as can be visualized in Fig. 2.

Fig. 2 Illustration of the different components of the detected signal and the stacked PSF in multifocal three-channel OCT.

Download Full Size | PDF

In order to eliminate the stacked point spread function which is computationally too expensive to deconvolve, we assume that the path-length difference between channels is much larger than the source coherence length so that cross-talk components 〈k_‖, k₀, z_m, z_n|S〉 where m ≠ n may be neglected. This can be achieved experimentally by inserting different lengths of fiber before each aperture [30]. With this assumption, the signal detected from the n^th aperture only consists of the back-scattering caused by the illumination from itself, resulting a simplified point spread function, as is described by

〈 k_{‖}, k_{0}, z_{n}, z_{n} | S 〉 = 2 π 〈 k_{0} | U_{0} 〉 \int d z {\tilde{h}}^{'} (k_{‖}, z, k_{0}, z_{n}, z_{n}) 〈 k_{‖}, z | η 〉,

where

{\tilde{h}}^{'} (k_{‖}, z, k_{0}, z_{n}, z_{n}) = \int d^{2} k_{‖}^{″} \frac{\tilde{g} (k_{‖}^{″}, k_{0}) \tilde{g} (k_{‖}^{″} - k_{‖}, k_{0})}{\sqrt{k_{0}^{2} - k_{‖}^{″ 2}}} e^{i [\sqrt{k_{0}^{2} - k_{‖}^{″ 2}} + \sqrt{k_{0}^{2} - {(k_{‖}^{″} - k_{‖})}^{2}}] (z_{n} - z)} .

Since only the signals with the same illumination and detection apertures will be considered, 〈k_‖, k, z_n, z_n|S〉 will be denoted as 〈k_‖, k, z_n|S〉 from now on for brevity. The method of stationary phase can be applied when |z_n − z| is large compared to Rayleigh range, where the phase is stationary around $k_{‖}^{″ (stat .)} = k_{‖} / 2$ [28], resulting in an simplified asymptotics result. The near-focus case can also evaluated by expanding the integral around the k′_‖ = k_‖/2 [28]. Generalizing both asymptotics results in the far-from-focus regime and the near-focus regime, the forward problem can be expressed as,

〈 k_{‖}, k_{0}, z_{n} | S 〉 = H (k_{‖}, k_{0}) \int d z 〈 - 2 \sqrt{k_{0}^{2} - {(k_{‖} / 2)}^{2}} | z - z_{n} 〉 w_{r} {(z - z_{n})}^{- 1} 〈 k_{‖}, z | η 〉,

where

H (k_{‖}, k_{0}) = {\begin{array}{r} (2 i k_{0}^{3} {NA}^{3}) / [π (k_{0}^{2} - k_{‖}^{2})] 〈 k_{0} | U_{0} 〉 {\tilde{g}}^{2} (k_{‖} / 2, k_{0}) & near focus \\ 4 π^{3} 〈 k_{0} | U_{0} 〉 {\tilde{g}}^{2} (k_{‖} / 2, k_{0}) & far from focus \end{array}

and w_r(z − z_n) is defined as the relative beam radius compared to w₀,

w_{r} (z - z_{n}) = \frac{w (z - z_{n})}{w_{0}} = \sqrt{1 + \frac{λ_{0}^{2} {(z - z_{n})}^{2}}{n^{2} π^{2} w_{0}^{4}}} .

2.2. Regularized inverse

At this stage we have formed a complete forward model of multifocal OCT which relates the susceptibility distribution |η〉 of the sample to the signal |S〉 detected from the aperture. The next step is to derive its inverse based on the forward problem, which enables us to reconstruct the susceptibility distribution from the detected signal. We start from Eq. (12) below, which is equivalent to Eq. (9),

\begin{array}{l} | S 〉 & = & \sum_{n = 1}^{N} \int d^{2} k_{‖} d k_{0} d z | k_{‖}, k_{0}, z_{n} 〉 H (k_{‖}, k_{0}) \\ \times 〈 - 2 \sqrt{k_{0}^{2} - {(k_{‖} / 2)}^{2}} | z - z_{n} 〉 w_{r} {(z - z_{n})}^{- 1} 〈 k_{‖}, z | η 〉 . \end{array}

The detected signal |S〉 can be seen as the result of two operators P and Z acting on the sample susceptibility |η〉, so that |S〉 = PZ|η〉. As is described by Eq. (13), the P operator takes the data and remaps them on to a different grid,

P = \int d^{2} k_{‖} d k_{0} | k_{‖}, k_{0} 〉 H (k_{‖}, k_{0}) 〈 k_{‖}, - 2 \sqrt{k_{0}^{2} - {(k_{‖} / 2)}^{2}} | .

The Z operator takes data from |η〉 and applies a depth-dependent weighting of w_r(z − z_n)⁻¹ in magnitude in the z dimension, that is

Z = \sum_{n = 1}^{N} \int d z | z - z_{n}, z_{n} 〉 w_{r} {(z - z_{n})}^{- 1} 〈 z | .

To reconstruct |η〉, we need to derive P⁺ and Z⁺, the pseudoinverses of the two operators P and Z, respectively. The reconstructed |η〉 then can be expressed as |η〉 = Z⁺P⁺ |S〉.

The pseudoinverse of P can be derived using P⁺ = (P^*P)⁻¹P^*, which is of the same form as the inverse operator in ISAM,

P^{+} = \int d^{2} k_{‖} d β | k_{‖}, β 〉 H^{- 1} (k_{‖}, \frac{1}{2} \sqrt{β^{2} + {k_{‖}}^{2}}) 〈 k_{‖}, \frac{1}{2} \sqrt{β^{2} + {k_{‖}}^{2}} | .

Using similar method, the exact inverse of operator Z can also be derived, that is

Z^{- 1} = \sum_{n = 1}^{N} \int d z | z 〉 w_{r} (z - z_{n}) 〈 z - z_{n}, z_{n} | .

When noise is considered, the exact inverse does not give the optimum results. For different noise models, the corresponding regularization can be applied by replacing w_r(z − z_n) with its function, so that

Z^{+} = \sum_{n = 1}^{N} \int d z | z 〉 R [w_{r} (z - z_{n})] 〈 z - z_{n}, z_{n} |,

to achieve the least square minimum norm solution.

Thus the general regularized solution to the inverse problem is

\begin{array}{l} 〈 r | η 〉 & = & \sum_{n = 1}^{N} R [w_{r} (z - z_{n})] \int d^{2} k_{‖} 〈 r_{‖} | k_{‖} 〉 \int d β 〈 z | - β 〉 \\ \times H^{- 1} (k_{‖}, \frac{1}{2} \sqrt{β^{2} + {k_{‖}}^{2}}) 〈 k_{‖}, \frac{1}{2} \sqrt{β^{2} + {k_{‖}}^{2}}, z_{n} | S 〉 . \end{array}

Here specifically for the multiplicative noise, Tikhonov regularization promises the least square minimum norm solution of the problem. The regularizer of Tikhonov regularization is

R_{Tik} [w_{r} (z - z_{n})] = \frac{w_{r} {(z - z_{n})}^{- 1}}{w_{r} {(z - z_{n})}^{- 2} + λ},

where λ is the regularization parameter. Thus, the Tikhonov regularized solution of the inverse problem is then

\begin{array}{l} 〈 r | η 〉 & = & \sum_{n = 1}^{N} \frac{1}{1 + λ w_{r} {(z - z_{n})}^{2}} [w_{r} (z - z_{n}) \int d^{2} k_{‖} 〈 r_{‖} | k_{‖} 〉 \int d β 〈 z | - β 〉 \\ \times H^{- 1} (k_{‖}, \frac{1}{2} \sqrt{β^{2} + {k_{‖}}^{2}}) 〈 k_{‖}, \frac{1}{2} \sqrt{β^{2} + {k_{‖}}^{2}}, z_{n} | S 〉], \end{array}

where the factor H⁻¹ depends on k_‖ and β.

2.3. Implementation

Equation (20) can be implemented based on a conventional ISAM code. Table 1 presents a step-by-step list to aid the implementation. To start with, multiple sets of OCT data should be acquired and pre-processed to compensate for dispersion and to remove the complex conjugate image. The later part of Eq. (20) in the square brackets can be seen as the ISAM reconstruction of each channel, corresponding to steps 2 through 6 in Table 1. Optional aberration compensation using computational adaptive optics (CAO) can be performed at this stage (not included in Eq. (20), and details will be covered in Section 2.5). The MISAM reconstruction is then the weighted sum of the ISAM reconstructions, with the weighting being 1/ [1 +λw_r(z − z_n)²], corresponding to step 7 and 8 in Table 1.

Table 1. A step-by-step guide for implementing MISAM.

View Table

2.4. Simulation

A computer simulation demonstrates the MISAM algorithm in Fig. 3. The sample used in the simulation consists of a series of point-scatterers on the x–z plane with 5 μm separation in depth between neighboring point-scatterers. An OCT system with an NA of 0.75, refractive index of 1, and wavelength of 800 nm to 1000 nm was used, resulting in a transverse resolution of 0.38 μm and a DOF of 1.02 μm in air.

Fig. 3 Panels (a)–(c) are three channels of OCT data, each focused at −20 μm, −60 μm, and −100 μm respectively, (the contrast has been lowered to show the signal far from focus). Panel (d) is the multifocal OCT reconstruction based on (a)–(c). Panels (e)–(g) are three channels of ISAM reconstruction based on (a)–(c) respectively. Panel (h) is the MISAM reconstruction based on (e)–(g). The MISAM image shows spatially invariant resolution and superior SNR.

Download Full Size | PDF

The results in Figs. 3(a)–3(c) show how the spreading in the OCT image acquired with focal plane at different depth becomes more severe at large distances from focus. Figure 3(d) shows the stitched image multifocal OCT synthetized from data in Figs. 3(a)–3(c). Although Fig. 3(d) sees an immediate improvement over Figs. 3(a)–3(c), the transverse resolution achieved by multifocal OCT still varies spatially.

The results in Figs. 3(e)–3(g) shows the corresponding ISAM reconstruction from the OCT data in Figs. 3(a)–3(c). Note that significant noise can be seen at lower part of Fig. 3(e) and upper part of Fig. 3(g). The regularization works as a spatial weighting to scale-down the noisy signal far from focus at each channel before the coherent combination. As a result, the noisy part contributes very little to the final reconstruction. The results in Fig. 3(h) show the regularized inverse solution works effectively in correcting the defocus and coherently combining the three datasets to achieve a reconstruction with spatially invariant resolution. Compared with the multifocal OCT image at Fig. 3(d), resolution in Fig. 3(h) again shows a significant improvement.

With the same parameters as above, Fig. 4 shows plots of theoretical transverse resolution and signal-to-noise ratio of MISAM compared with multifocal OCT and single channel ISAM. In Fig. 4(a), MISAM shows superior transverse resolution, especially at larger distance from each focus, compared to multifocal OCT. Moreover, in Fig. 4(b), MISAM shows overall superior SNR characteristics, compared to single-focal ISAM. In this case, MISAM provides 200% to 700% improvement in SNR far from focus.

Fig. 4 (a) Theoretical transverse resolution for three channels of OCT, multifocal OCT and MISAM, based on a three-aperture multifocal OCT system with focal plane at −20 μm, −60 μm, and −100 μm. (b) Theoretical SNR of single-focal ISAM and MISAM, based on the same multifocal OCT system.

Download Full Size | PDF

2.5. Experimental results

Experimental data have been taken and the regularized inverse for MISAM has been applied to demonstrate the MISAM reconstruction. Three datasets were acquired separately using a single channel SD-OCT setup. The focal plane depths of the scanning beams were adjusted before each scan. As shown in Fig. 5, three datasets (a)–(c) were acquired with focal plane at −1200 μm, −1000 μm, and −800 μm, respectively. The OCT system has a NA of 0.1, and the optical spectrum recorded by the spectrometer had central wavelength of 1330 nm and FWHM bandwidth of 150 nm. Thus, the setup has a transverse resolution of ∼ 7 μm and a Rayleigh range of ∼ 50 μm in air. The sample consists of a collection of sub-resolution sized TiO₂ particles embedded in a semitransparent silicone matrix, which has a refractive index of ∼ 1.4.

Fig. 5 Panels (a)–(c) are three channels of OCT data, each focused at −1200 μm, −1000 μm, and −800 μm respectively. Panel (d) is the multifocal OCT reconstruction based on (a)–(c). Panels (e)–(g) are three channels of ISAM reconstruction based on (a)–(c) respectively. Panel (h) is the MISAM reconstruction based on (e)–(g). Three outsets on the right side show the en face planes of Multifocal OCT, MISAM and MISAM with computed adaptive optics (CAO) at −300 μm.

Download Full Size | PDF

Figures 5(a)–5(c) are the OCT reconstructions of each set of signals. It can be observed that the reconstructions become defocused far away from their focal plane; thus, each set provides ∼ 200 μm usable effective DOF. After applying the image stitching method over Figs. 5(a)–5(c), Fig. 5(d) is the multifocal OCT reconstruction of the datasets. It provides an effective DOF of ∼ 800 μm, from −600 μm to −1400 μm.

The MISAM algorithm was applied to the three channels of OCT data. Because the datasets were sequentially acquired, they were not phase coherent, thus the absolute value of each dataset were taken before the regularized combination were applied.

Because inexpensively priced lens sets were used for the imaging system, moderate aberration, mainly spherical, were observed in the experimental data, which degrades the imaging quality. In order to eliminate the effect, we have incorporated computational adaptive optics [31, 32] to the datasets. A phase plate H(k_‖) = exp[−ikΦ_h(−2πz_f k_‖/k₀)], where Φ_h represents the weighted sum of Zernike polynomials, was multiplied to 〈k_‖, z|η〉 for each z. The weighting parameters for the Zernike polynomials were searched by a script for each z to optimize for the maximum pixel intensity in the en face plane.

By comparing Fig. 5(d) with Fig. 5(h), it can be observed that after regularized combination and CAO, the reconstruction in h) shows not only higher resolution especially far from focus, but also high SNR throughout all depths. A set of en face planes at the same depth form the multifocal OCT reconstruction, MISAM reconstruction without CAO and MISAM reconstruction with CAO are shown on the right side. It can be observed that MISAM corrects for the spreading seen in multifocal OCT images. Based on the MISAM image, CAO further corrects the aberration of the system to achieve more ideal point spread functions.

Figure 6 shows a comparison of the three-dimensional rendering of dual-channel multifocal OCT on the left and dual-channel MISAM on the right using only datasets 1 and 3, and also two en face plots at the middle depth between the two focal planes. In addition to the obvious advantage of MISAM at the top and bottom of the reconstruction, this comparison shows that MISAM is able to provide higher quality reconstruction with fewer number of foci than multifocal OCT.

Fig. 6 Panels (a) and (d) show the three-dimensional volume rendering of size of 500 μm by 500 μm by 1800 μm in depth of multifocal OCT and MISAM constructed with channel 1 and 3. Panels (b) and (c) show the en-face plot of multifocal OCT and MISAM about 200 μm from both focal planes of channel 1 and 3.

Download Full Size | PDF

3. Conclusion

We have demonstrated multifocal ISAM to achieve improvements in SNR and image quality over OCT alone. In this paper, the method was demonstrated with three focal planes, though the technique can be readily extended to accommodate arbitrary numbers of foci.

Data for MISAM can be acquired either in parallel using a multifocal OCT setup or in series using a conventional OCT setup. For serial acquisition, the OCT scan is repeated multiple times with a variety of different focal depths. Serial acquisition trades longer acquisition time (proportional to number of foci desired) and possible sample drift for lower experimental complexity.

With the aid of computational adaptive optics, inexpensive lenses can be used for data acquisition without severely affecting the imaging quality. This trades computational cost for lower equipment cost. The computational cost for MISAM is proportional to the number of foci desired. Reconstruction at a similar scale to Fig. 5 can be easily performed using a current mainstream notebook computer (Thinkpad X220 with Intel i7-2620m CPU) within one minute. Furthermore, it is possible to implement the algorithm on a GPU using parallel programming language, such as CUDA, to achieve real time reconstruction [33, 34].

MISAM is based on the assumption of the first Born approximation, as is OCT. Under such approximation, single scattering is assumed and higher-order scattering is neglected. In many applications of OCT, this proves to be a good approximation when dealing with many kinds of samples. So for any sample imaged with OCT, MISAM can be applied to produce images with high resolution, large effective DOF, and high SNR. Reconstruction artifacts may appear if the sample is strongly scattering. Although the effective DOF is improved by using MISAM, the effective DOF may still be bounded by penetration depth for samples with low transparency at the scanning wavelength.

We have shown that MISAM can realize the benefits of both single-focus ISAM and multifocal OCT by producing a reconstruction with spatially invariant resolution and theoretically high signal-to-noise over a large range. MISAM provides a physically meaningful method to combine multiple ISAM data from the same sample with different foci to scale up the effective DOF of the image.

Acknowledgments

This research was supported in part by a grant from the National Institutes of Health (NIBIB, R01 EB012479).

References and links

1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254, 1178–1181 (1991). [CrossRef] [PubMed]

2. S. A. Boppart, B. E. Bouma, C. Pitris, G. J. Tearney, J. G. Fujimoto, and M. E. Brezinski, “Forward-imaging instruments for optical coherence tomography,” Opt. Lett. 22, 1618–1620 (1997). [CrossRef]

3. S. A. Boppart, B. E. Bouma, C. Pitris, J. F. Southern, M. E. Brezinski, and J. G. Fujimoto, “In vivo cellular optical coherence tomography imaging,” Nat. Med. 4, 861–865 (1998). [CrossRef] [PubMed]

4. J. G. Fujimoto, C. Pitris, S. A. Boppart, and M. E. Brezinski, “Optical coherence tomography: An emerging technology for biomedical imaging and optical biopsy,” Neoplasia 2, 9–25 (2000). [CrossRef] [PubMed]

5. J. F. de Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearney, and B. E. Bouma, “Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coherence tomography,” Opt. Lett. 28, 2067–2069 (2003). [CrossRef] [PubMed]

6. M. Choma, M. Sarunic, C. Yang, and J. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11, 2183–2189 (2003). [CrossRef] [PubMed]

7. E. A. Swanson, J. A. Izatt, M. R. Hee, D. Huang, C. P. Lin, J. S. Schuman, C. A. Puliafito, and J. G. Fujimoto, “In vivo retinal imaging by optical coherence tomography,” Opt. Lett. 18, 1864–1866 (1993). [CrossRef] [PubMed]

8. J. A. Izatt, M. R. Hee, E. A. Swanson, C. P. Lin, D. Huang, J. S. Schuman, C. A. Puliafito, and J. G. Fujimoto, “Micrometer-scale resolution imaging of the anterior eye in-vivo with optical coherence tomography,” Arch. Ophthalmol. 112, 1584–1589 (1994). [CrossRef] [PubMed]

9. J. S. Schuman, M. R. Hee, C. A. Puliafito, C. Wong, T. Pedutkloizman, C. P. Lin, E. Hertzmark, J. A. Izatt, E. A. Swanson, and J. G. Fujimoto, “Quantification of nerve fiber layer thickness in normal and glaucomatous eyes using optical coherence tomography,” Arch. Ophthalmol. 113, 586–596 (1995). [CrossRef] [PubMed]

10. M. Wojtkowski, A. Kowalczyk, R. Leitgeb, and A. F. Fercher, “Full range complex spectral optical coherence tomography technique in eye imaging,” Opt. Lett. 27, 1415–1417 (2002). [CrossRef]

11. V. Guedes, J. S. Schuman, E. Hertzmark, G. Wollstein, A. Correnti, R. Mancini, D. Lederer, S. Voskanian, L. Velazquez, H. M. Pakter, T. Pedut-Kloizman, J. G. Fujimoto, and C. Mattox, “Optical coherence tomography measurement of macular and nerve fiber layer thickness in normal and glaucomatous human eyes,” Ophthalmology 110, 177–189 (2003). [CrossRef] [PubMed]

12. M. Wojtkowski, V. Srinivasan, J. G. Fujimoto, T. Ko, J. S. Schuman, A. Kowalczyk, and J. S. Duker, “Three-dimensional retinal imaging with high-speed ultrahigh-resolution optical coherence tomography,” Ophthalmology 112, 1734–1746 (2005). [CrossRef] [PubMed]

13. J. Welzel, “Optical coherence tomography in dermatology: a review,” Skin Res. Technol. 7, 1–9 (2001). [CrossRef] [PubMed]

14. M. C. Pierce, J. Strasswimmer, B. H. Park, B. Cense, and J. F. de Boer, “Advances in optical coherence tomography imaging for dermatology,” J. Invest. Dermatol. 123, 458–463 (2004). [CrossRef] [PubMed]

15. T. Gambichler, G. Moussa, M. Sand, D. Sand, P. Altmeyer, and K. Hoffmann, “Applications of optical coherence tomography in dermatology,” J. Dermatol. Sci. 40, 85–94 (2005). [CrossRef] [PubMed]

16. S. A. Boppart, G. J. Tearney, B. E. Bouma, J. F. Southern, M. E. Brezinski, and J. G. Fujimoto, “Noninvasive assessment of the developing Xenopus cardiovascular system using optical coherence tomography,” Proc. Natl. Acad. Sci. U. S. A. 94, 4256–4261 (1997). [CrossRef] [PubMed]

17. R. Leitgeb, M. Villiger, A. Bachmann, L. Steinmann, and T. Lasser, “Extended focus depth for Fourier domain optical coherence microscopy,” Opt. Lett. 31, 2450–2452 (2006). [CrossRef] [PubMed]

18. C. Blatter, B. Grajciar, C. M. Eigenwillig, W. Wieser, B. R. Biedermann, R. Huber, and R. A. Leitgeb, “Extended focus high-speed swept source oct with self-reconstructive illumination,” Opt. Express 19, 12141–12155 (2011). [CrossRef] [PubMed]

19. J. Holmes, S. Hattersley, N. Stone, F. Bazant-Hegemark, and H. Barr, “Multi-channel Fourier domain OCT system with superior lateral resolution for biomedical applications,” in “Biomedical Optics (BiOS) 2008” (International Society for Optics and Photonics, 2008), 68470O. [CrossRef]

20. J. Holmes, “Theory and applications of multi-beam OCT,” in “1st Canterbury Workshop and School in Optical Coherence Tomography and Adaptive Optics” (International Society for Optics and Photonics, 2008), 713908.

21. J. Holmes and S. Hattersley, “Image blending and speckle noise reduction in multi-beam OCT,” in “SPIE BiOS: Biomedical Optics” (International Society for Optics and Photonics, 2009), 71681N.

22. T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Inverse scattering for optical coherence tomography,” J. Opt. Soc. Am. A 23, 1027–1037 (2006). [CrossRef]

23. T. S. Ralston, D. L. Marks, S. A. Boppart, and P. S. Carney, “Inverse scattering for high-resolution interferometric microscopy,” Opt. Lett. 31, 3585–3587 (2006). [CrossRef] [PubMed]

24. T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Interferometric synthetic aperture microscopy,” Nat. Phys. 3, 129–134 (2007). [CrossRef]

25. B. J. Davis, S. C. Schlachter, D. L. Marks, T. S. Ralston, S. A. Boppart, and P. S. Carney, “Nonparaxial vector-field modeling of optical coherence tomography and interferometric synthetic aperture microscopy,” J. Opt. Soc. Am. A 24, 2527–2542 (2007). [CrossRef]

26. B. J. Davis, T. S. Ralston, D. L. Marks, S. A. Boppart, and P. S. Carney, “Autocorrelation artifacts in optical coherence tomography and interferometric synthetic aperture microscopy,” Opt. Lett. 32, 1441–1443 (2007). [CrossRef] [PubMed]

27. B. J. Davis, D. L. Marks, T. S. Ralston, P. S. Carney, and S. A. Boppart, “Interferometric synthetic aperture microscopy: Computed imaging for scanned coherent microscopy,” Sensors 8, 3903–3931 (2008). [CrossRef] [PubMed]

28. P. S. Carney, S. G. Adie, and S. A. Boppart, “Interferometric synthetic aperture microscopy,” in Emerging Imaging Technologies in Medicine (CRC Press, 2012), pp. 293–301.

29. T. S. Ralston, S. G. Adie, D. L. Marks, S. A. Boppart, and P. S. Carney, “Cross-validation of interferometric synthetic aperture microscopy and optical coherence tomography,” Opt. Lett. 35, 1683–1685 (2010). [CrossRef] [PubMed]

30. V. X. Yang, N. Munce, J. Pekar, M. L. Gordon, S. Lo, N. E. Marcon, B. C. Wilson, and I. A. Vitkin, “Micromachined array tip for multifocus fiber-based optical coherence tomography,” Opt. Lett. 29, 1754–1756 (2004). [CrossRef] [PubMed]

31. S. G. Adie, B. W. Graf, A. Ahmad, P. S. Carney, and S. A. Boppart, “Computational adaptive optics for broadband optical interferometric tomography of biological tissue,” Proc. Natl. Acad. Sci. U. S. A. 109, 7175–7180 (2012). [CrossRef] [PubMed]

32. S. G. Adie, B. W. Graf, A. Ahmad, B. Darbarsyah, S. A. Boppart, and P. S. Carney, “The impact of aberrations on object reconstruction with interferometric synthetic aperture microscopy,” in “SPIE BiOS” (International Society for Optics and Photonics, 2011), 78891O.

33. T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Real-time interferometric synthetic aperture microscopy,” Opt. Express 16, 2555–2569 (2008). [CrossRef] [PubMed]

34. A. Ahmad, N. D. Shemonski, S. G. Adie, H.-S. Kim, W.-M. W. Hwu, P. S. Carney, and S. A. Boppart, “Real-time in vivo computed optical interferometric tomography,” Nat. Photon. 7, 444–448 (2013). [CrossRef]

1	Acquire OCT data with various focal plane depth {z₁, z₂,..., z_N} and preprocess, 〈r_‖, k₀, z_n\|S〉 = 〈r_‖, k₀\|S_n〉.
2	For each dataset, take the 2-D Fourier transform with respect to transverse coordinates, $〈 k_{‖}, k_{0} \| {\tilde{S}}_{n} 〉 = \frac{1}{2 π} \int d^{2} r_{‖} exp (i k_{‖} \cdot r_{‖}) 〈 r_{‖}, k_{0} \| S_{n} 〉$ .
3	For each dataset, apply the inverse filter (optional), 〈k_‖, k₀\|S̃′_n〉 = H⁻¹ (k_‖, k₀) 〈k_‖, k₀\|S̃_n〉.
4	For each dataset, perform the resampling, $〈 k_{‖}, - β \| {\tilde{S}}_{n}^{″} 〉 = {〈 k_{‖}, k_{0} \| {\tilde{S}}_{n}^{'} 〉 \|}_{k_{0} = \frac{1}{2} \sqrt{β^{2} + {k_{‖}}^{2}}}$ .
5	For each dataset, take the 3-D inverse Fourier transform with respect to k_‖ and −β, 〈r_‖, z\|S″_n〉 = (2π)^−3/2 ∫ d²k_‖dβ exp(−ik_‖ · r_‖ + iβz) 〈k_‖, −β\|S̃″_n〉.
6	For each dataset, multiply with the relative beam radius function, 〈r_‖, z\|η_n〉 = w_r(z − z_n) 〈r_‖, z\|S″_n〉.
7	For each dataset, multiply with the weighting, $〈 r_{‖}, z \| η_{n}^{'} 〉 = \frac{1}{1 + λ w_{r} {(z - z_{n})}^{2}} 〈 r_{‖}, z \| η_{n} 〉$ .
8	Take the sum of all datasets, $〈 r_{‖}, z \| η_{MISAM} 〉 = \sum_{n = 1}^{N} 〈 r_{‖}, z \| η_{n}^{'} 〉$ .

Multifocal interferometric synthetic aperture microscopy

Abstract

1. Introduction

2. Multifocal interferometric synthetic aperture microscopy

2.1. Forward model

2.2. Regularized inverse

2.3. Implementation

2.4. Simulation

2.5. Experimental results

3. Conclusion

Acknowledgments

References and links

Cited By

Figures (6)

Tables (1)

Equations (20)

Optics Express