Blur resolved OCT: full-range interferometric synthetic aperture microscopy through dispersion encoding

Jonathan H. Mason; Mike E. Davies; Pierre O. Bagnaninchi

doi:10.1364/OE.379216

1. Introduction

Optical coherence tomography (OCT) offers high resolution non-invasive imaging of tissues and other semi-transparent materials [1–4]. Through the reflection interference between a reference and sample arm, the structure of scatterers along depth are encoded. In spectral-domain OCT (SD-OCT) [5], this interferometry signal is diffracted onto a detector array, from which the one-dimensional structure (A-scan) can be reconstructed through an inverse fast discrete Fourier transform (IFFT). The three-dimensional structure of the specimen can then be formed by raster scanning the sample and combining the resulting profiles.

One deficit of this simplistic scheme is that A-scans are not independent, due to the widening of the beam away from the focal plane of the lens in the sample arm. With this, objects appear blurred in the out-of-focus region of the image, leading to a non-uniform resolution with depth. With interferometric synthetic aperture microscopy (ISAM), Ralston et al. [6–8] showed that this effect may be actively compensated for by resampling the spatial frequency components and filtering, which effectively combines information from neighboring A-scans, and leads to a uniform resolution with depth.

Another potential deficit of SD-OCT is due to the measurements at the spectrometer being real intensities. Therefore, its Fourier transform will be conjugate symmetric, effectively halving the available depth range. In practice, one often ensures the absence of objects in the negative optical delay region, and then ignore the superfluous mirror image after applying an IFFT. There are several hardware approaches to utilize the entire range, such as placing a phase modulator in the sample arm [9,10], offsetting the scanning mirror pivot [11], or measuring the quadrature component of the interferometry signal [12], although these solutions increase system complexity and reduce scanning rate due to requiring several measurements [13].

It is also possible to differentiate the conjugate terms, by introducing a dispersion discrepancy between sample and reference arms. This is well approximated as a non-linear phase lag, which acts in an opposite direction for the mirrored complement. Therefore, after compensating for dispersion in one direction, the mirror component becomes ‘doubly dispersed’ leading to a blurring and distinction from the desired sharpened signal. In dispersion encoded full-range (DEFR) OCT [13], one takes a greedy optimization approach to resolving the object, by iteratively removing the blurred mirror associated with the highest magnitude component. This uses the implicit approximation that A-scans are independent. There have been several extensions to this, including removing several components on each iteration [14,15], and removing autocorrelation artefacts also [16]. It was recently shown that DEFR may indeed allow faithful reconstruction even under subsampling regimes [17]. Interestingly, there are strong parallels between these approaches and radio frequency interference suppression in synthetic aperture radar (SAR) [18].

In order to perform full-range imaging from real spectral measurements, one must accept an inherent sampling deficit from inferring as many complex parameters as real samples, which is the case in DEFR methods. As commonly employed in the field of compressed sensing [19,20], one can introduce a sparsity constraint that allows the faithful reconstruction of sparse signals from few measurements. It has been demonstrated that images from OCT are typically sparse in some domain, and compressed sensing restoration methods have been proven successful [13,21–23].

It was recently shown in [23] that the ISAM resampling can be utilized in a model-based iterative reconstruction (MBIR) in the half-range setting under sub-sampling, with sparsity promoting regularization as is used in compressed sensing. In this work, we extend MBIR to the dispersion encoded full-range setting. We present an accelerated MBIR algorithm, along with an enhanced variation, and evaluate it with synthetic and real full-range measurements of several complex samples.

1.1 Novel contributions

We demonstrate, for the first time, full-range dispersion encoded ISAM. As well as showing how this may be achieved through a naive combination of two existing algorithms (DEFR+ISAM), we propose a novel MBIR optimization approach as a solution. Unlike the greedy DEFR methods, this has the potential to exploit the shared information between A-scans, such as multidimensional sparsity in the ISAM refocussed space. We provide analysis from quantitative simulation and real data, where we show the feasibility of reconstructing refocussed full-range measurements with computational methods, and observe a significant performance gain of MBIR over alternative approaches.

2. Background

In this section, we describe the dispersion encoded full-range ISAM model exploited by our MBIR method. Here, we wish to reconstruct the complex susceptibility, $\boldsymbol {\eta }\in \mathbb {C}^{N}$, from the real spectrometer measurements, $\boldsymbol {s}\in \mathbb {R}^{N}$, where $N = N_xN_yN_z$ with $N_x$ and $N_y$ the number of lateral measurements in $x$ and $y$, and $N_z$ the number of axial measurements (also the resolution of the spectrometer). Since we wish to infer $N$ complex numbers with $2N$ unknowns from only $N$ measurements, this represents a clear sampling deficit, which we attempt to overcome by exploiting sparsity.

The ISAM model tells us that the 3D Fourier transform of the object, given as

(1)$${\boldsymbol{H}}(\boldsymbol{q}_x,\boldsymbol{q}_y,\boldsymbol{\beta}) = \mathcal{F}_\mathrm{3D}(\boldsymbol{\eta}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z})),$$

where $\boldsymbol {q}_x$, $\boldsymbol {q}_y$ and $\boldsymbol {\beta }$ are transverse and axial spatial frequencies respectively, is related to the transverse Fourier transform of the complex interferometry signal

(2)$$\boldsymbol{S}(\boldsymbol{q}_x,\boldsymbol{q}_y,\boldsymbol{k}) = \mathcal{F}_\leftrightarrow(\boldsymbol{s}_c),$$

where $\boldsymbol {s}_c$ is the complex interferometry signal, of which we directly measure its real part; we use the notation $\mathcal {F}_\leftrightarrow (\cdot )$ for the Fourier transform in the transverse spatial dimensions, and $\boldsymbol {k}$ is the wavenumber.

The ISAM relationship can then be expressed as [24]

(3)$$\boldsymbol{S}(\boldsymbol{q}_x,\boldsymbol{q}_y,\boldsymbol{k}) = \boldsymbol{B}(\boldsymbol{q}_x,\boldsymbol{q}_y,\boldsymbol{k})\odot{\boldsymbol{H}}\left(\boldsymbol{q}_x,\boldsymbol{q}_y,\boldsymbol{\beta}\right)|_{\boldsymbol{\beta} = -\sqrt{4k^2-\boldsymbol{q}_x^2-\boldsymbol{q}_y^2}},$$

where $\boldsymbol {B}(\boldsymbol {q}_x,\boldsymbol {q}_y,k)$ is a filter to account for the lens and intensity drop off [25] and a frequency warping effect is seen through the relation $\beta = -\sqrt {4k^2-\boldsymbol {q}_x^2-\boldsymbol {q}_y^2}$, $\odot$ represents element-wise multiplication. This resampling is known as the Stolt mapping, and is used in the fields of seismology and SAR [26,27].

Since the resampling through interpolation, filtering and Fourier transform are all linear operations, we can express this by [13]

(4)$$\boldsymbol{s}_c = \boldsymbol{K}\boldsymbol{\eta},$$

where $\boldsymbol {K}$ is a matrix representing the mapping from susceptibility image to the complex spectroscopic signal.

We only directly measure the real part of this signal [3], which we can express as

(5)$$\boldsymbol{s} = \Re(\boldsymbol{K}\boldsymbol{\eta}),$$

where $\Re (\cdot )$ selects the real part. This is equivalent to

(6)$$\boldsymbol{s} = \frac{1}{2}\left[\boldsymbol{s}_c+\bar{\boldsymbol{s}}_c^*\right],$$

where $\bar {\boldsymbol {s}}_c^*$ is the complex spectrum from the conjugate component we wish to suppress. In practice, there will also be background and autocorrelation signals on top of $s$ measured at the spectrometer. Throughout this work, we assume the background can be easily removed and the autocorrelation components are small. Since autocorrelation artifacts have been shown to be influenced in both ISAM [28] and DEFR [16], their treatment may be useful in highly scattering specimens.

In either the half-range of full-range setting [5,13], when a dispersion discrepancy between sample and reference arms exists, this may be accurately modelled through a non-linear phase term given as

(7)$$e^{j\boldsymbol{\phi}} = \exp\left(j\sum_{i=1}^{N_p} a_i(\boldsymbol{k}-k_0)^i\right),$$

where $k_0$ is the central frequency component, $a_i$ are the polynomial coefficients, and $N_p$ is the order. In practice, $N_p = 3$ is usually sufficient to capture significant dispersion effects and $\boldsymbol {a}$ may be found through an autofocus algorithm [5]. With this, Eq. (6) becomes

(8)$$\boldsymbol{s}_d = \frac{1}{2}\left[\boldsymbol{s}_c\odot e^{j\boldsymbol{\phi}}+\bar{\boldsymbol{s}}_c^*\odot e^{-j\boldsymbol{\phi}}\right],$$

where $\boldsymbol {s}_d$ represents the real measurements as in Eq. (6) with the inclusion of dispersion, and the phase shift has opposite effect on each of the conjugate components [13]. We highlight that the derivation of Eq. (8) as in [13] uses a thin sample approximation, whereby dispersion between the layers of the sample is not considered.

In standard half-range imaging, dispersion compensated reconstruction can be performed by

(9)$$\mathcal{F}^{-1}_\updownarrow(\boldsymbol{s}_d\odot e^{-j\boldsymbol{\phi}}) = \frac{1}{2}\left[\mathcal{F}^{-1}_\updownarrow(\boldsymbol{s}_c)+\mathcal{F}^{-1}_\updownarrow(\bar{\boldsymbol{s}}_c^*\odot e^{-2j\boldsymbol{\phi}})\right],$$

where we use the notation $\mathcal {F}^{-1}_\updownarrow$ for an IFFT in the axial dimension. If the object of interest only occupies the positive delay area, then $\mathcal {F}^{-1}_\updownarrow (\bar {\boldsymbol {s}}_c^*\odot e^{-2j\boldsymbol {\phi }})$ will not overlap with the desired signal, and can be easily ignored. When there is an overlap, and given sufficient dispersion encoding, then DEFR [13] can approximately remove the unwanted blurred artefacts in an iterative manner.

2.1 DEFR+ISAM

One approach we offer to achieve full-range ISAM under dispersion encoding, is through combining the two existing algorithms as presented in [13,24].

Firstly, DEFR attempts to solve the following optimization problem based on the dispersion encoding from Eq. (9) in a greedy fashion

(10)$$\hat{\boldsymbol{z}} = \mathop{\arg\!\min}\limits_{\boldsymbol{z}} \|2\Re\{\mathcal{F}_\updownarrow(\boldsymbol{z})\odot e^{j\boldsymbol{\phi}}\}-\boldsymbol{s}_d\|_2^2,$$

The iterative method in [13] works like Matching Pursuit (MP) [29,30], by updating $\boldsymbol {z}$ one component at a time, giving a solution that is sparse.

From here, one can then estimate the susceptibility by applying the ISAM back-projection operator from Eq. (4) as

(11)$$\boldsymbol{\eta} = \boldsymbol{K}^H(\mathcal{F}_\updownarrow(\hat{\boldsymbol{z}})+\boldsymbol{r}),$$

where $\boldsymbol {r}$ is the residual term from the approximate solution of Eq. (10) as defined in [13].

Since DEFR operates on each A-scan independently, it is unable to exploit multidimensional sparsity in the image domain. This is especially relevant when applied before ISAM, as one expects significant portions of the image to be out of focus, which will decrease the observed sparsity. Secondly, the application of Eq. (4) calculates the back-projection, which is only a crude approximation of the inverse of the system model in $\boldsymbol {K}$.

3. Method

To overcome the potential deficits of the simple two-step approach suggested in Section 2.1, we propose that full-range DEFR–ISAM can be achieved by solving of the following optimization problem

(12)$$\tilde{\boldsymbol{\eta}} = \mathop{\arg\!\min}\limits_{\boldsymbol{\eta}}\frac{1}{2}\|\Re(\hat{\boldsymbol{K}}\boldsymbol{\eta})-\boldsymbol{s}_d\|_2^2 + \lambda \boldsymbol{w}^T|\boldsymbol{\eta}|,$$

where $\boldsymbol {w}^T|\boldsymbol {\eta }|$ is a weighted $\ell _1$ norm with the weighting vector $\boldsymbol {w}\in \mathbb {R}^N$. The unweighted $\ell _1$ penalty function — usually denoted as $\|\boldsymbol {\eta }\|_1$ — more commonly employed in compressed sensing can be achieved through choosing $\boldsymbol {w} = \boldsymbol {1}$, although choosing a nonuniform $\boldsymbol {w}$ can account for the increased sparsity with depth typical in OCT. We also adopt the compact notation in Eq. (12) for the dispersion corrected ISAM model

(13)$$\hat{\boldsymbol{K}} \equiv \mathrm{diag}(e^{j\boldsymbol{\phi}})\boldsymbol{K}.$$

The $\ell _1$ penalty term in Eq. (12) assumes the specimen is spatially sparse after refocussing and the conjugate artifacts have been removed. However, this method and algorithm are also easily extendable for any convex regularization function, such as total variation (TV) and wavelet sparsity [23], which may be more appropriate for specimens with different spatial structure.

Problems with the form of Eq. (12) have been extremely well studied in the field of compressed sensing [19,20], in which many algorithms for its solution have been developed. In this work, we opt for the fast iterative thresholding shrinkage algorithm (FISTA) [31], which is an accelerated gradient descent method with soft-thresholding to minimize the $\ell _1$ function. FISTA applied to the objective function in Eq. (12) is given in Algorithm 1, with the weighted soft-thresholding operator defined as

(14)$$\mathcal{T}_{\lambda w_i}(u_i) = \frac{u_i\max\left(|u_i|-\frac{\lambda w_i}{N},0\right)}{\max\left(|u_i|-\frac{\lambda w_i}{N},0\right)+\frac{\lambda w_i}{N}}.$$

3.1 Parameters

There are two elements in Algorithm 1 that must be appropriately chosen: the regularization constant $\lambda$, and the termination condition.

For terminating the method we adopt the relative residual stopping condition from [32] defined through the value

(15)$$r_r^k = \frac{\|\boldsymbol{z}^k-\boldsymbol{\eta}^k\|}{\max\left\{\|\boldsymbol{g}^k\|,\|\boldsymbol{z}^k-\boldsymbol{\eta}^k+\boldsymbol{g}^k\|\right\}+\epsilon},$$

where $\boldsymbol {g}^k = \frac {1}{N}\hat {\boldsymbol {K}}^H\left (\Re (\hat {\boldsymbol {K}}\boldsymbol {\eta }^{k})-\boldsymbol {s}\right )$, and $\epsilon$ is a small constant to ensure a non-zero denominator. With this, one terminates the iterations once $\|r_r^k\|\;<\;tol$, where $tol$ is the desired tolerance for convergence (we use $1\times 10^{-3}$ in our testing).

The regularization parameter $\lambda$, in combination with the optional weighting vector $\boldsymbol {w}$, will control the level of sparsity and hence conjugate artifact suppression in the reconstruction. It will also have an impact on the convergence, with larger values of $\lambda$ typically requiring fewer iterations to reach a given $tol$. From our testing on a range of different samples, we find $[0.1,0.8]$ to be a range providing good performance.

The introduction of the weighting term $\boldsymbol {w}$ gives the option to use non-uniform regularization throughout the reconstructed image. From numerical and empirical testing, we have found an increase in $\boldsymbol {w}$ with depth produces better images. An intuition for this is that there are typically fewer photons collected from deep in the specimen due to scattering in the top layers, hence the local signal–to–noise ratio is reduced, so one should give less weight to the data fidelity term and more to the sparsity prior. We find that using a simple linearly increasing $\boldsymbol {w}$ with depth from 0.5 to 1, provides a large performance gain in our quantitative experiment presented in Section 4.3. However, we suggest that other weighting schemes may be superior for different lenses or focal plane positions.

As an optional final step in Algorithm 1, one can add the final residual onto the solution after termination of the MBIR. This is similar to the inclusion of residual in DEFR [13], and is also used in radar imaging [33]. The advantage of this is to retain non-sparse low intensity but informative signal, which may otherwise be suppressed by the regularization, but is too small to contribute significant conjugate artifacts. For applications such as optical coherence elastography [34], retaining the low level speckle structure is critical for displacement tracking. Another advantage of including the residual term is that more aggressive regularization can be used within the MBIR optimization, which results in faster convergence and hence faster reconstruction speeds.

3.2 Efficient and robust ISAM through non-uniform FFT

The proposed approach requires repeated realization of the ISAM model within the optimization algorithm and therefore it is necessary to have an accurate yet efficient implementation of the ISAM operator. In this work, this is realized through the non-uniform FFT (NUFFT) algorithm [35].

We rewrite the ISAM operator as

(16)$$\boldsymbol{K}\cdot \equiv N_r\mathcal{F}_\leftrightarrow^{-1}(\mathrm{diag}(\boldsymbol{b})\mathcal{N}(\cdot)),$$

with the matrix $\boldsymbol {K}$ as in Eq. (4), $\mathcal {N}(\cdot )$ is the NUFFT operator and $\boldsymbol {b}$ is vector representation of the filter $\boldsymbol {B}(\boldsymbol {q},\boldsymbol {k})$ in Eq. (3). We will henceforth treat the unfiltered solution in this work, whereby we exclude $\boldsymbol {b}$, which has been shown to have minimal qualitative effect [6,24].

For a standard ISAM as in [24], this is equivalent to back-projection as

(17)$$\boldsymbol{K}^H\boldsymbol{s} = \mathcal{N}'(\mathcal{F}_\leftrightarrow(\boldsymbol{s})).$$

4. Experiments

4.1 Materials and measurements

All measurements were acquired using a Wasatch Photonics 800nm SD-OCT system, with 2048 spectrometer elements. The system’s imaging arm included a length of optical fibre, sufficient to introduce a dispersion discrepancy against the reference mirror. We recorded 1024 A-scans over a 2 mm lateral distance using its standard protocol, and extracted the raw spectrometer data for processing. In each case the focal point was adjusted, by eye, to lie within the sample, and at the zero delay position.

Preprocessing from the raw data consisted of background subtraction, obtained by averaging across A-scans, and non-linear calibration from detector element to wavenumber, according to parameters from the manufacturer.

The samples used were as follows:

1. Beaded gel: TiO$_2$ micro-beads suspended in 2% agarose gel, at a concentration of 1 mg/ml. The powdered TiO$_2$ ($<$5 $\mu$m, Sigma-Aldrich) was dispersed throughout the gel before curing, through combination of stirring, pipette agitation and sonication.
2. Tape: a roll of GiftWrap Scotch adhesive tape.
3. Cucumber: a slice of cucumber flesh sectioned and blotted to remove excess moisture.

The reconstruction was performed retrospectively on a commodity PC with an Intel i7-8700 CPU, 16 GB of RAM, and an NVIDIA Titan Xp GPU. All software was written in and run with Matlab.

4.2 Methods under test

The various methods and their implementation that we analyze in this experimental section are as follows:

• Direct IFFT: dispersion compensation is applied to the measurements, with the polynomial coefficients $a_2$ and $a_3$ in Eq. (7) found through the autofocus method of [5], followed by an IFFT. We then used the same dispersion parameters in each of the following methods.
• DEFR: the method as described in [13], with the residual included as it resulted in a superior quantitative performance. We ran the algorithm for 500 iterations to ensure convergence.
• ISAM: we follow dispersion compensation with ISAM reconstruction as described in [24], through the NUFFT adjoint operator in Eq. (17), although without the cropping step to reduce the data to half-range.
• DEFR+ISAM: the combination of DEFR followed by full-range ISAM. Although this method is not described in the literature, it is a naive approach that we use as a point of comparison for our proposed MBIR.
• MBIR: a simple implementation of Algorithm 1, with no weighting ($\boldsymbol {w} = \boldsymbol {1}$), and the residual step included. In every case we use $\lambda =0.5$, which provided a good balance between numerical accuracy and visual clarity, and use the termination condition $tol=1\times 10^{-3}$.
• MBIR+: a more advanced implementation of Algorithm 1 with the weighting term $\boldsymbol {w}$ linearly increasing with depth from 0.5 to 1, and the residual term also included.

4.3 Quantitative test

To objectively assess the performance of our approach relative to the alternatives described in Sec. 4.2, we performed analysis from pseudo-full-range data, against a ground truth obtained from real half-range measurements. In generating these measurements, we adjusted the focal plane to approximately lie halfway through the positive delay range, and ensured the negative delay was absent of scattering media. This process is illustrated in Fig. 1, with each step detailed as follows:

1. Preprocessing: the raw measurements are processed to remove background, resample from detector element to discrete frequency (k-space), a Hilbert transform is taken to ensure zero negative delay, and finally dispersion compensation is applied.
2. ISAM mapping: the processed measurements are reconstructed through fully sampled ISAM. This image represents a ground truth against which subsampled reconstructions can be assessed. The ground truth is also scaled by a factor 0.5, to match the magnitude change from taking real part during measurement synthesis.
3. IFFT: an IFFT is applied to the processed data to map into image space. Due to the Hilbert transform, all values in the negative delay region will be zero.
4. Measurement synthesis: the negative delay region of the IFFT image is truncated to place a virtual zero-delay at the dashed yellow line. The pseudo-full-range measurements are then generated by taking the axial FFT, applying the dispersion factor, and finally taking the real part.
5. Direct mapping: directly applying an IFFT to the under-sampled measurements produces conjugate artifacts and lateral blurring in out–of–focus regions. We include this to demonstrate the magnitude of the artifacts, and the relative performance of the computational approaches.
6. Reconstruction: we produce reconstructions with our proposed and tested methods described in Sec. 4.2 and evaluate their quantitative accuracy against the ground truth.

Fig. 1. Flow diagram showing the generation of synthetic full-range data from real measurements with ground truth. Descriptions of each processing step are detailed in the enumerated list in Sec. 4.3.

Download Full Size | PDF

The samples used in the quantitative study were the beaded gel, tape and cucumber, as described in Sec. 4.1, although with the focal plane positioned at half the positive delay range. Reconstructions are shown in Fig. 2, Fig. 3 and Fig. 4, with results of the quantitative analysis presented in Table 1. We quantify the root mean squared error (RMSE) between the ground truth and raw reconstructed data. On top of this, we use the peak signal–to–noise ratio (PSNR) and structural similarity index (SSIM) [36] calculated on the 16-bit images displayed in Fig. 2, Fig. 3 and Fig. 4 after taking the logarithm and scaling the data. We use these additional metrics as they correlate better with human perception of image quality [36], and are in the format in which the data is likely to be interpreted.

Fig. 2. Reconstructions from synthetic beaded gel data with ground truth and measurements produced as shown in Fig. 1. The methods and their implementation are detailed in Sec. 4.2.

Download Full Size | PDF

Table 1. Quantitative results from synthetic full-range data (best method in bold).

View Table

Fig. 3. Reconstructions from synthetic Scotch tape data with ground truth and measurements produced as shown in Fig. 1. The methods and their implementation are detailed in Sec. 4.2.

Download Full Size | PDF

Fig. 4. Reconstructions from synthetic cucumber data with ground truth and measurements produced as shown in Fig. 1. The methods and their implementation are detailed in Sec. 4.2.

Download Full Size | PDF

From the quantitative results in Table 1, the performance of our proposed approaches are compelling. Firstly, the simple two-step DEFR+ISAM is effective, and significantly outperforms each of its constituent algorithms, with a significant reduction in error over DEFR for the three samples, but also tends to introduce lateral artifacts in areas from which conjugate components are removed. On top of this, we achieve a 17–24% reduction in RMSE between DEFR+ISAM and our simple MBIR implementation and a 26–34% reduction against our MBIR+, which adequately justifies our new algorithm. Although the PSNR and SSIM values between MBIR and DEFR+ISAM are reasonably similar, again our MBIR+ shows a significant gain over both of these, which highlights the benefit of the weighting term.

By comparing the reconstructed images in Fig. 2, Fig. 3 and Fig. 4, the performance and action of the methods can be seen: DEFR effectively mitigates the conjugate artifacts that are present in the IFFT and ISAM; ISAM brings the areas into focus away from the focal plane that are blurred in IFFT and DEFR; and DEFR+ISAM achieves the combined benefit of its two constituent steps. As well as also having this combined effect, MBIR+ produces images very close to the ground truth, and significantly outperforms any existing model or combination thereof.

The advantages from MBIR approaches over DEFR+ISAM suggests that being able to exploit the multidimensional sparsity present in the focussed image after ISAM resampling is valuable.

4.4 Real data validation

Whilst the results in Section 4.3 allow us to objectively analyze reconstructions against a ground truth, and between various methods, we also validated this with real full-range measurements. In this section, we apply the same methods to extend the depth range from measurements directly from a commercial OCT system. Since we no longer have a ground truth, this analysis will necessarily be qualitative, with our observations qualified by the results in Sec. 4.3.

Firstly, we have shown reconstructions from the similiar beaded gel and cucumber samples as used in Sec. 4.3. The samples are shown in Fig. 5, Fig. 6 and Fig. 7 respectively, but with the focal plane positioned at the zero-delay position and taking the fully-sampled spectroscopic data.

Fig. 5. Reconstructions of beaded gel sample with real measurements. The various methods and their implementation are detailed in Sec. 4.2.

Download Full Size | PDF

Fig. 6. Reconstructions of Scotch tape sample with real measurements. The various methods and their implementation are detailed in Sec. 4.2.

Download Full Size | PDF

Fig. 7. Reconstructions of cucumber with real measurements. The various methods and their implementation are detailed in Sec. 4.2.

Download Full Size | PDF

Several observations can be made from the full-range reconstruction of the beaded gel and tape images in Fig. 5 and Fig. 6 respectively. Similarly to the images in Sec 4.3, DEFR is effective at suppressing conjugate components, ISAM is effective in refocussing blurred structures away from the focal plane, and DEFR+ISAM or MBIR combine both of these benefits. Unlike the DEFR+ISAM, where there is visible lateral blurring around some of the structures, especially in the inset images, both MBIR and MBIR+ produce sharper looking images. Another feature from DEFR+ISAM, especially of the beaded gel, is an apparent double lobe effect above and below the focal plane. Since this is not present in any other image, it is likely an artifact, but one that does not seem to affect MBIR. We suggest it is an effect of subtracting the conjugate components, which will only be approximated by the implicit DEFR dispersion model, and can hence lead to cumulative errors, especially in areas around the focal plane with greatest intensity. Between the MBIR+ and MBIR, the former does have visibly reduced residual artifacts, whilst preserving all structural information.

From the cucumber tissue images in Fig. 7, the gain in structural clarity of MBIR+ over DEFR+ISAM can be seen in the inset images. The cell boundaries have greater definition, and some fine structures that are distorted in the other cases, are visible with MBIR. This further confirms the advantage of our proposed method.

In our tested samples, there are no significant autocorrelation artifacts around the zero delay visible, but we suggest these are likely to increase with samples of higher scattering intensity.

5. Conclusions

We have demonstrated full-range ISAM through dispersion encoding, and presented a MBIR algorithm for its implementation. While a naive DEFR+ISAM is an alternative way to achieve this, it does not fully exploit the ISAM signal model or multidimensional sparsity, resulting in structural degradation and observed artifacts. In contrast, MBIR produces images of enhanced structural clarity and significantly lower errors. Within this, we have adopted an efficient NUFFT implementation of ISAM, which is numerically stable through many iterations. Ongoing and future work includes reducing the computational time of MBIR, actively modelling autocorrelation artifacts, extending the method to work efficiently with 3D datasets, and to explore compelling biomedical applications.

Funding

Engineering and Physical Sciences Research Council (EP/P031250/1); H2020 European Research Council (ERC-ADG-2015-694888); Royal Society (Wolfson Research Merit Award).

Acknowledgments

The authors sincerely thank Graham Anderson from the University of Edinburgh, for assistance creating the beaded gel phantom. We acknowledge NVIDIA Corporation for kindly donating the Titan Xp GPU. This work has made use of the resources provided by the Edinburgh Compute and Data Facility (ECDF).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

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

2. 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(1-2), 9–25 (2000). [CrossRef]

3. P. H. Tomlins and R. K. Wang, “Theory, developments and applications of optical coherence tomography,” J. Phys. D: Appl. Phys. 38(15), 2519–2535 (2005). [CrossRef]

4. Y.-Z. Liu, F. A. South, Y. Xu, P. S. Carney, and S. A. Boppart, “Computational optical coherence tomography,” Biomed. Opt. Express 8(3), 1549–1574 (2017). [CrossRef]

5. M. Wojtkowski, V. J. Srinivasan, T. H. Ko, J. G. Fujimoto, A. Kowalczyk, and J. S. Duker, “Ultrahigh-resolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12(11), 2404–2422 (2004). [CrossRef]

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

7. 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(5), 1027 (2006). [CrossRef]

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

9. E. Götzinger, M. Pircher, R. A. Leitgeb, and C. K. Hitzenberger, “High speed full range complex spectral domain optical coherence tomography,” Opt. Express 13(2), 583–594 (2005). [CrossRef]

10. D. Y. Kim, J. S. Werner, and R. J. Zawadzki, “Comparison of phase-shifting techniques for in vivo full-range, high-speed Fourier-domain optical coherence tomography,” J. Biomed. Opt. 15(5), 056011 (2010). [CrossRef]

11. L. An and R. K. Wang, “Use of a scanner to modulate spatial interferograms for in vivo full-range Fourier-domain optical coherence tomography,” Opt. Lett. 32(23), 3423–3425 (2007). [CrossRef]

12. Y. Mao, S. Sherif, C. Flueraru, and S. Chang, “3x3 Mach-Zehnder interferometer with unbalanced differential detection for full-range swept-source optical coherence tomography,” Appl. Opt. 47(12), 2004–2010 (2008). [CrossRef]

13. B. Hofer, B. Považay, B. Hermann, A. Unterhuber, G. Matz, and W. Drexler, “Dispersion encoded full range frequency domain optical coherence tomography,” Opt. Express 17(1), 7–24 (2009). [CrossRef]

14. S. Witte, M. Baclayon, E. J. G. Peterman, R. F. G. Toonen, H. D. Mansvelder, and M. L. Groot, “Single-shot two-dimensional full-range optical coherence tomography achieved by dispersion control,” Opt. Express 17(14), 11335–11349 (2009). [CrossRef]

15. B. Hofer, B. Považay, A. Unterhuber, L. Wang, B. Hermann, S. Rey, G. Matz, and W. Drexler, “Fast dispersion encoded full range optical coherence tomography for retinal imaging at 800 nm and 1060 nm,” Opt. Express 18(5), 4898–4919 (2010). [CrossRef]

16. F. Köttig, P. Cimalla, M. Gärtner, and E. Koch, “An advanced algorithm for dispersion encoded full range frequency domain optical coherence tomography,” Opt. Express 20(22), 24925–24948 (2012). [CrossRef]

17. L. Yi and L. Sun, “Full-depth compressive sensing spectral-domain optical coherence tomography based on a compressive dispersion encoding method,” Appl. Opt. 57(31), 9316–9321 (2018). [CrossRef]

18. S. I. Kelly and M. E. Davies, “RFI suppression and sparse image formation for UWB SAR,” Int. Radar Symp. Proc. 2, 655–660 (2013).

19. E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52(2), 489–509 (2006). [CrossRef]

20. D. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]

21. X. Liu and J. U. Kang, “Compressive SD-OCT: the application of compressed sensing in spectral domain optical coherence tomography,” Opt. Express 18(21), 22010–22019 (2010). [CrossRef]

22. N. Mohan, I. Stojanovic, W. C. Karl, B. E. A. Saleh, and M. C. Teich, “Compressed sensing in optical coherence tomography,” Proc. SPIE 7570, 75700L (2010). [CrossRef]

23. J. H. Mason, Y. Reinwald, Y. Yang, S. Waters, A. E. Haj, and P. O. Bagnaninchi, “Model-based iterative reconstruction for spectral-domain optical coherence tomography,” Proc. SPIE 10883, 108830O (2019). [CrossRef]

24. T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Real-time interferometric synthetic aperture microscopy,” Opt. Express 16(4), 2555–2569 (2008). [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(9), 2527–2542 (2007). [CrossRef]

26. R. H. Stolt, “Migration by Fourier Transform,” Geophysics 43(1), 23–48 (1978). [CrossRef]

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(6), 3903–3931 (2008). [CrossRef]

28. 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(11), 1441–1443 (2007). [CrossRef]

29. S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Trans. Signal Process. 41(12), 3397–3415 (1993). [CrossRef]

30. M. F. Duarte, M. A. Davenport, M. B. Wakin, and R. G. Baraniuk, “Sparse Signal Detection from Incoherent Projections,” Int. Conf. Acoust. Spee. 3, 305–308 (2006). [CrossRef]

31. A. Beck and M. Teboulle, “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]

32. T. Goldstein, C. Studer, and R. G. Baraniuk, “A field guide to forward-backward splitting with a FASTA implementation,” CoRR abs/1411.3406 (2014).

33. S. I. Kelly, C. Du, G. Rilling, and M. E. Davies, “Advanced image formation and processing of partial synthetic aperture radar data,” IET Signal Process. 6(5), 511–520 (2012). [CrossRef]

34. B. F. Kennedy, P. Wijesinghe, and D. D. Sampson, “The emergence of optical elastography in biomedicine,” Nat. Photonics 11(4), 215–221 (2017). [CrossRef]

35. J. A. Fessler and B. P. Sutton, “Nonuniform fast Fourier transforms using min-max interpolation,” IEEE Trans. Signal Process. 51(2), 560–574 (2003). [CrossRef]

36. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. on Image Process. 13(4), 600–612 (2004). [CrossRef]

sample	metric	direct	ISAM	DEFR	DEFR+ISAM	MBIR	MBIR+
beaded gel	RMSE	5.23	3.02	4.27	0.797	0.648	0.548
	PSNR	8.16	8.65	15.5	18.7	19.1	21.3
	SSIM	0.147	0.341	0.199	0.371	0.370	0.595
Scotch tape	RMSE	4.61	2.68	3.76	1.05	0.874	0.777
	PSNR	9.11	9.70	12.8	15.1	14.8	16.0
	SSIM	0.178	0.365	0.168	0.314	0.349	0.411
cucumber	RMSE	0.860	0.493	0.704	0.388	0.295	0.258
	PSNR	16.7	18.0	19.2	22.4	23.1	25.1
	SSIM	0.105	0.315	0.147	0.332	0.369	0.517

Blur resolved OCT: full-range interferometric synthetic aperture microscopy through dispersion encoding

Abstract

1. Introduction

1.1 Novel contributions

2. Background

2.1 DEFR+ISAM

3. Method

3.1 Parameters

3.2 Efficient and robust ISAM through non-uniform FFT

4. Experiments

4.1 Materials and measurements

4.2 Methods under test

4.3 Quantitative test

4.4 Real data validation

5. Conclusions

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (7)

Tables (1)

Equations (17)

Optics Express