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

TL;DR

This paper introduces a computational approach for full-range interferometric synthetic aperture microscopy that doubles OCT depth range and improves resolution using dispersion encoding and a model-based iterative reconstruction method.

Contribution

The paper presents a novel MBIR method incorporating dispersion encoding and NUFFT for fast, stable full-range ISAM imaging without hardware changes.

Findings

Successfully doubles OCT depth range

Enhances resolution away from focal plane

Outperforms existing methods in complex sample imaging

Abstract

We present a computational method for full-range interferometric synthetic aperture microscopy (ISAM) under dispersion encoding. With this, one can effectively double the depth range of optical coherence tomography (OCT), whilst dramatically enhancing the spatial resolution away from the focal plane. To this end, we propose a model-based iterative reconstruction (MBIR) method, where ISAM is directly considered in an optimization approach, and we make the discovery that sparsity promoting regularization effectively recovers the full-range signal. Within this work, we adopt an optimal nonuniform discrete fast Fourier transform (NUFFT) implementation of ISAM, which is both fast and numerically stable throughout iterations. We validate our method with several complex samples, scanned with a commercial SD-OCT system with no hardware modification. With this, we both demonstrate full-range ISAM imaging, and significantly outperform combinations of existing methods.

Tables

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

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

Equations

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

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

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

$${\mathbf{s}}_c=\mathbf{K}\mathbf{\eta },$$

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

$$\mathbf{s}=\frac{1}{2}\left[{\mathbf{s}}_c+{\bar{\mathbf{s}}}_c^{\ast }\right],$$

$$e^{j\mathbf{\varphi }}=\exp \left(j\sum _{i=1}^{N_p}{a}_i(\mathbf{k}-k_0{)}^i\right),$$

$${\mathbf{s}}_d=\frac{1}{2}\left[{\mathbf{s}}_c\odot e^{j\mathbf{\varphi }}+{\bar{\mathbf{s}}}_c^{\ast }\odot e^{-j\mathbf{\varphi }}\right],$$

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

$$\hat{\mathbf{z}}=\underset{\mathbf{z}}{\mathrm{argmin}}\Vert 2\Re \{{\mathcal{F}}_{\updownarrow }(\mathbf{z})\odot e^{j\mathbf{\varphi }}\}-{\mathbf{s}}_d{\Vert }_2^2,$$

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

$$\tilde{\mathbf{\eta }}=\underset{\mathbf{\eta }}{\mathrm{argmin}}\frac{1}{2}\Vert \Re (\hat{\mathbf{K}}\mathbf{\eta })-{\mathbf{s}}_d{\Vert }_2^2+\lambda {\mathbf{w}}^T|\mathbf{\eta }|,$$

$$\hat{\mathbf{K}}\equiv \mathrm{diag}(e^{j\mathbf{\varphi }})\mathbf{K}.$$

$${\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}}.$$

$$r_r^k=\frac{\Vert {\mathbf{z}}^k-{\mathbf{\eta }}^k\Vert }{\max \left\{\Vert {\mathbf{g}}^k\Vert ,\Vert {\mathbf{z}}^k-{\mathbf{\eta }}^k+{\mathbf{g}}^k\Vert \right\}+ϵ},$$

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

$${\mathbf{K}}^H\mathbf{s}={\mathcal{N}}^{\prime }({\mathcal{F}}_{\leftrightarrow }(\mathbf{s})).$$

Full text

1 Introduction

Optical coherence tomography (OCT) offers high resolution non-invasive imaging of tissues and other semi-transparent materials Huang et al., (1991); Fujimoto et al., (2000); Tomlins and Wang, (2005); Liu et al., (2017). Through the reflection interference between a reference and sample arm, the structure of scatterers along depth are encoded. In spectral-domain OCT (SD-OCT) Wojtkowski et al., (2004), 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. Ralston et al., 2006a ; Ralston et al., 2006b ; Ralston et al., (2007) 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 Götzinger et al., (2005); Kim et al., (2010), offsetting the scanning mirror pivot An and Wang, (2007), or measuring the quadrature component of the interferometry signal Mao et al., (2008), although these solutions increase system complexity and reduce scanning rate due to requiring several measurements Hofer et al., (2009).

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 Hofer et al., (2009), 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 Witte et al., (2009); Hofer et al., (2010), and removing autocorrelation artefacts also Köttig et al., (2012). It was recently shown that DEFR may indeed allow faithful reconstruction even under subsampling regimes Yi and Sun, (2018). Interestingly, there are strong parallels between these approaches and radio frequency interference suppression in synthetic aperture radar (SAR) Kelly and Davies, (2013).

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 Candes et al., (2006); Donoho, (2006), 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 Hofer et al., (2009); Liu and Kang, (2010); Mohan et al., (2010); Mason et al., (2019).

It was recently shown in Mason et al., (2019) 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_{x}N_{y}N_{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

[TABLE]

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

[TABLE]

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 Ralston et al., (2008)

[TABLE]

where $\boldsymbol{B}(\boldsymbol{q}_{x},\boldsymbol{q}_{y},k)$ is a filter to account for the lens and intensity drop off Davis et al., 2007b 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 Stolt, (1978); Davis et al., (2008).

Since the resampling through interpolation, filtering and Fourier transform are all linear operations, we can express this by Hofer et al., (2009)

[TABLE]

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 Tomlins and Wang, (2005), which we can express as

[TABLE]

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

[TABLE]

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 Davis et al., 2007a and DEFR Köttig et al., (2012), their treatment may be useful in highly scattering specimens.

In either the half-range of full-range setting Wojtkowski et al., (2004); Hofer et al., (2009), when a dispersion discrepancy between sample and reference arms exists, this may be accurately modelled through a non-linear phase term given as

[TABLE]

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 Wojtkowski et al., (2004). With this, Eq. (6) becomes

[TABLE]

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 Hofer et al., (2009). We highlight that the derivation of Eq. (8) as in Hofer et al., (2009) 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

[TABLE]

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 Hofer et al., (2009) 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 Hofer et al., (2009); Ralston et al., (2008).

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

[TABLE]

The iterative method in Hofer et al., (2009) works like Matching Pursuit (MP) Mallat and Zhang, (1993); Duarte et al., (2006), 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

[TABLE]

where $\boldsymbol{r}$ is the residual term from the approximate solution of Eq. (10) as defined in Hofer et al., (2009).

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

[TABLE]

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

[TABLE]

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 Mason et al., (2019), 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 Candes et al., (2006); Donoho, (2006), in which many algorithms for its solution have been developed. In this work, we opt for the fast iterative thresholding shrinkage algorithm (FISTA) Beck and Teboulle, (2009), 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

[TABLE]

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 Goldstein et al., (2014) defined through the value

[TABLE]

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 Alg. 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 Hofer et al., (2009), and is also used in radar imaging Kelly et al., (2012). 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 Kennedy et al., (2017), 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 Fessler and Sutton, (2003).

We rewrite the ISAM operator as

[TABLE]

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 Ralston et al., 2006a ; Ralston et al., (2008).

For a standard ISAM as in Ralston et al., (2008), this is equivalent to back-projection as

[TABLE]

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:

Beaded gel: TiO2 micro-beads suspended in 2% agarose gel, at a concentration of 1 mg/ml. The powdered TiO2 ($<$5 $\mu$m, Sigma-Aldrich) was dispersed throughout the gel before curing, through combination of stirring, pipette agitation and sonication. 2. 2.

Tape: a roll of GiftWrap Scotch adhesive tape. 3. 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 Wojtkowski et al., (2004), followed by an IFFT. We then used the same dispersion parameters in each of the following methods.

• •

DEFR: the method as described in Hofer et al., (2009), 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 Ralston et al., (2008), 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 Alg. 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 Alg. 1 with the weighting term $\boldsymbol{w}$ linearly increasing with depth from 0.5 to 1, and the residual term also included. We find the residual term allows for more aggressive regularization and hence used $\lambda=0.5$ throughout the simulated and real tests, although this resulted in visible suppression of structure in direct MBIR.

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:

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. 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. 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. 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. 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. 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.

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. 3, with results of the quantitative analysis presented in Tab. 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) Zhou Wang et al., (2004) calculated on the 16-bit images displayed in Fig. 2 and Fig. 3 after taking the logarithm and scaling the data. We use these additional metrics as they correlate better with human perception of image quality Zhou Wang et al., (2004), and are in the format in which the data is likely to be interpreted.

From the quantitative results in Tab. 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–340% 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 and Fig. 3, 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 analyse 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 same 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.

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 an 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

UK Engineering and Physical Sciences Research Council (EP/P031250/1); 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.