1. Introduction {#sec1}
===============

Compared with other biomedical imaging technologies, optical coherence tomography (OCT) has the advantages of a relatively high resolution but a relatively low set-up cost, making it a very successful imaging modality to migrate from research to industry products \[[@r1], [@r2], [@r3]\]. To date, OCT technology has been through three generations: time-domain OCT (TD-OCT), spectral-domain OCT (SD-OCT), and swept-source OCT (SS-OCT) \[[@r4], [@r5]\]. Resolution of an OCT system has been improved from 10 ∼ 20*μ*m to several *μ*m. Data acquisition speed has also increased from 2 A-lines per second in 1991 to 480,000 A-lines per second in 2011 \[[@r6]\]. The applications of OCT also extend from ophthalmology to other areas such as endoscopy, coronary atherosclerosis, dermatology, dental, pulmonary medicine, gynecology, and urology \[[@r2],[@r7]\]. Among the three kinds of OCT systems, SD-OCT and SS-OCT have the advantage of a higher signal-to-noise ratio (SNR) and a higher sensitivity at a large depth position compared to TD-OCT \[[@r8], [@r9]\].

In an SD-OCT system, the measurements are collected uniformly in wavelength *λ*. This is because the measurement device, a spectrometer, operates by tuning to different wavelengths. On the other hand, in reconstructing the image, we use a Fourier transform of the measurement with respect to the wavenumber *k*, where *k* = 2*π*/*λ*. Therefore, the measurements are in fact nonuniform samples in *k*-space \[[@r10]\]. Customarily, the SD-OCT measurements undergo a resampling process to generate a uniformly-spaced data in *k*-space, where the resampling is achieved through linear, or cubic, interpolations \[[@r10], [@r11]\]. However, these methods need to face various challenges, such as improving system sensitivity falloff along the depth direction, increasing system local SNR, or suppressing sidelobe artifacts \[[@r12], [@r13]\].

An alternative to the re-sampling process is recently developed, which makes use of the non-uniform discrete Fourier transform (NUDFT) for SD-OCT signal reconstruction \[[@r12], [@r13], [@r14]\]. NUDFT is a modified version of the discrete Fourier transform, aiming to compute DFT directly on those signal points that are not uniformly sampled \[[@r15], [@r16]\]. It has been shown that with NUDFT, the reconstructed signal in SD-OCT presents less sensitivity falloff at the deeper position \[[@r13]\]. However, this gain comes at the expense of a slower computation speed. Fast versions of NUDFT are developed with GPU programming \[[@r12]\], but they tend to sacrifice some reconstruction quality. Other approaches to deal with the nonuniformly sampling issue include the development of a more sophisticated linear-in-wavenumber (linear-*k*) spectrometer \[[@r17]\] for measurement collection.

In this paper, we present a purely computational alternative by using inverse imaging \[[@r18]\]. We model the SD-OCT system with a set of linear equations, and by solving an inverse problem we can also directly reconstruct a cross-section of the sample instead of estimating each Aline independently. Because of the ill-posedness of the inverse problem, we need to impose a regularization term in the optimization function. In this work, we use the total variation (TV) norm in the regularization \[[@r19], [@r20], [@r21]\], which is known to enhance the edges in the images and reduce the background noise \[[@r22], [@r23]\].

The paper is organized as follows. In Section 2, we develop the model of the SD-OCT system. Then, in Section 3, we discuss the TV method used for signal reconstruction. In Section 4, experimental results are presented to demonstrate the feasibility and advantages of our technique. Concluding remarks are then given in Section 5.

2. SD-OCT system modeling {#sec2}
=========================

[Fig. 1](#g001){ref-type="fig"} shows an SD-OCT system. A low coherence Superluminescent Diode (SLD) is used as an illumination source. Most of the source light is split equally through the coupler into a reference arm and a sample arm. Then, the reflected beams from a reference mirror and the beams reflected from the sample are coupled again into a detecting arm for signal detection. In the detecting arm, a spectrometer consisting of a grating and a line-scan CCD is used for collecting measurements. During a detector exposure, each pixel of the CCD collects a measurement corresponding to one source wavelength. With all pixels of the CCD, the data for one A-line along a sample's axial direction are collected. To obtain all A-line data of a sample, two galvano mirrors, labeled GM1 and GM2, are used to scan through the sample's transversal plane.

![The SD-OCT system diagram.](boe-3-4-741-g001){#g001}

Let *Ĩ*(*k*) be the system measurement as a function of the wavenumber *k*, and *G*(*k*) the source power spectrum. The SD-OCT system measurements for one A-line is \[[@r10]\] $$\widetilde{I}(k) = G(k)\left\{ {p_{r}^{2} + 2{\int_{0}^{\infty}{p_{r}p_{s}(z)\text{cos}(2a_{s}kz)\text{d}z + {\int_{0}^{\infty}{\int_{0}^{\infty}{p_{s}(z)p_{s}(z^{\prime})\text{e}^{\text{j}2ka_{s}(z - z^{\prime})}\text{d}z\text{d}z^{\prime}}}}}}} \right\},$$where *p~r~* is the reflective ratio of the reference mirror (normally we take *p~r~* = 1), *p~s~*(*z*) is the reflective ratio of the sample varying with depth value *z*, *a~s~* is the sample refractive index (also assumed to be unity), and $\text{j} = \sqrt{- 1}$. [Eq. (1)](#e01){ref-type="disp-formula"} consists of three terms. The first term corresponds to the reflected light from the reference mirror, the second term records the interference between the reflected lights from the mirror and the sample, while the third term is the cross-interference between the reflected lights from different sample layers. Of the three, the second term is the most important, and from it we aim to reconstruct the sample signal *p~s~*(*z*).

Conventionally, the signal *p~s~*(*z*) is further assumed to be symmetric with respect to *z* = 0, and therefore it can be estimated from measurements *Ĩ*(*k*) using the fast Fourier transform (FFT) \[[@r13]\]. On the other hand, we note that if the source spectrum *G*(*k*) can be estimated independently, we can substract $G(k)p_{r}^{2}$ (since *p~r~* = 1) from [Eq. (1)](#e01){ref-type="disp-formula"} and represent the last term as an error *e*(*k*), which is usually small compared with the other two terms. Thus, the SD-OCT measurement is simplified to $$I(k) = 2G(k){\int_{0}^{\infty}{p_{s}(z)\text{cos}(2kz)\text{d}z + e(k).}}$$In the discrete domain, this measurement can be written as $$I(k_{m}) = 2G(k_{m})\sum\limits_{n = 0}^{N - 1}{p_{s}(z_{n})\text{cos}(2k_{m}z_{n}) + e(k_{m}),}$$where we have discretized to *N* sections in the axial direction, and *m* ∈ {0*,*1*,\...,M* -- 1} denotes the *m*th pixel of the line-scan CCD camera. The value of *N* is not required to be the same as *M*, which is different from the FFT or NUDFT methods. However, the largest *N* should be limited by the largest detectable depth, which is determined by the spectrometer's spectral resolution. Writing [Eq. (3)](#e03){ref-type="disp-formula"} into a matrix formulation, we have $$\begin{bmatrix}
{I(k_{0})} \\
{I(k_{1})} \\
 \vdots \\
{I(k_{M - 1})} \\
\end{bmatrix} = \begin{bmatrix}
H_{0,0} & H_{0,1} & \cdots & H_{0,N - 1} \\
H_{1,0} & H_{1,1} & \cdots & H_{1,N - 1} \\
 \vdots & \vdots & \ddots & \vdots \\
H_{M - 1,0} & H_{M - 1,1} & \cdots & H_{M - 1,N - 1} \\
\end{bmatrix}\begin{bmatrix}
{p_{s}(z_{0})} \\
{p_{s}(z_{1})} \\
 \vdots \\
{p_{s}(z_{N - 1})} \\
\end{bmatrix} + \begin{bmatrix}
{e(k_{0})} \\
{e(k_{1})} \\
 \vdots \\
{e(k_{M - 1})} \\
\end{bmatrix},$$where *H~m~*~,~*~n~* = 2*G*(*k~m~*) cos (2*k~m~z~n~*).

To write the above compactly, we use **y** = \[*I*(*k*~0~) *I*(*k*~1~) \... *I*(*k~M~*~−1~)\]*^T^*, **x** = \[*p~s~*(*z*~0~) *p~s~*(*z*~1~) \... *p~s~*(*z~N~*~−1~)\]*^T^*, and **e** = \[*e*(*k*~0~) *e*(*k*~1~) ··· *e*(*k~M~*~−1~)\]*^T^* to represent the SD-OCT measurement vector, one A-line sample signal, and the detection error vector, respectively. Matrix H is defined as the SD-OCT system impulse response matrix with *H~m,n~* in the (*m*,*n*)th position. Using these notations, we can write the measurement process as **y** = H**x** + **e**. This formulation represents the basic idea for the FFT and NUDFT methods as well. Furthermore, the SD-OCT measurements for *L* number of A-lines can be represented as $$\text{Y} = \text{HX} + {\text{E},}$$where matrices Y and E are of size *M* × *L*, and X is of size *N* × *L*. Each column of these quantities represents the measurement for one A-line.

3. Signal reconstruction using TV regularization {#sec3}
================================================

From [Eq. (5)](#e05){ref-type="disp-formula"}, our task is to reconstruct X given the observation Y and the system H. This is an inverse problem. We solve it by minimizing the Frobenius norm \[[@r24]\] (which is equivalent to the *ℓ*~2~ norm for a vector) for signal fidelity together with a total variation regularization, i.e. $$\text{X}_{\text{min}} = \text{arg}\,\underset{\text{X}}{\text{min}}\left\| {\text{HX} - \text{Y}} \right\|_{F}^{2} + \alpha\left\| \text{X} \right\|_{\text{TV}},$$where *α* \> 0 is a penalty parameter, and \|\|X\|\|~TV~ is the total variation of X defined as $$\left\| \text{X} \right\|_{\text{TV}} = \sum\limits_{n = 0}^{N - 1}{\sum\limits_{l = 0}^{L - 1}{\sqrt{\left| {{({\nabla_{n}\text{X}})}_{n,l})} \right|^{2} + \left| {{({\nabla_{l}\text{X}})}_{n,l})} \right|^{2}},}}$$with (∇*~n~*X)*~n~*~,~*~l~* = X*~n~*~+1,~*~l~* − X*~n~*~,~*~l~* if *n* \< *N* − 1 and is zero when *n* = *N* − 1, and (∇*~l~* X)*~n~*~,~*~l~* = X*~n~*~,~*~l~*~+1~ − X*~n~*~,~*~l~* if *l* \< *L* − 1 and is zero when *l* = *L* − 1.

Several approaches exist to solve this minimization problem. In our work, we make use of the augmented Lagrangian technique, which has shown to be a fast converging TV method \[[@r25]\], and with results comparable to other TV minimization techniques \[[@r26]\]. We incorporate an additional intermediate matrix U such that we solve for $${({\text{X}_{\text{min}},\text{U}_{\text{min}}})} = \text{arg}\,\underset{\text{X},\text{U}}{\text{min}}\left\| {\text{HX} - \text{Y}} \right\|_{F}^{2} + \alpha_{1}\left\| {\text{X} - \text{U}} \right\|_{F}^{2} + \alpha_{2}\left\| \text{U} \right\|_{\text{TV}},$$where *α*~1~ \> 0 and *α*~2~ \> 0 are two predefined positive penalty parameters. We then use an alternating minimization algorithm where each iteration consists of the following two steps: (using the superscript (*i*) to denote the *i*th step) $$\left\{ {\begin{array}{l}
{\text{X}^{(i)} = \text{arg}\,\underset{\text{X}}{\text{min}}\left\| {\text{HX} - \text{Y}} \right\|_{F}^{2} + \alpha_{1}\left\| {\text{X} - \text{U}^{(i - 1)}} \right\|_{F}^{2}} \\
{\text{U}^{(i)} = \text{arg}\,\underset{\text{U}}{\text{min}}\alpha_{1}\left\| {\text{X}^{(i)} - \text{U}} \right\|_{F}^{2} + \alpha_{2}\left\| \text{U} \right\|_{\text{TV}}} \\
\end{array}.} \right.$$Conjugate gradient (CG) and Chambolle's projection algorithms are used to search for X^(*i*)^ and U^(*i*)^, respectively. Details about these two algorithms can be found in \[[@r26], [@r27], [@r28], [@r29]\]. Note that the usage of the matrix U does not change the minimization problem in [Eq. (6)](#e06){ref-type="disp-formula"}, but is there to improve the convergence. In addition, we use the difference between X^(*i*)^ and U^(*i*)^ as a stop criteria for the iterative algorithm. The parameters *α*~1~ and *α*~2~ are used to balance the measurement error and the TV regularization penalty in [Eq. (8)](#e08){ref-type="disp-formula"}. Large *α*~1~ and *α*~2~ values are for the case of a high measurement error. However, if they are set at very big values, the images will be significantly smoothed in the reconstruction process.

4. Experimental results {#sec4}
=======================

Our SD-OCT system is equipped with a SLD source with the center wavelength of 850 nm, a bandwidth of 50 nm, and a Gaussian-shaped power spectrum. In air, the system has a depth resolution of 6.38*μ*m, and a maximum detectable depth of 2.42 mm. The mirror in the reference arm is silver-coated with a two-inch diameter. System measurements are collected using a customized spectrometer consisting of a one-inch volume-phase holographic (VPH) grating and a line-scan camera. As depicted in [Fig. 1](#g001){ref-type="fig"}, the camera and the galvano mirrors are controlled by a computer.

In the first experiment, 800 A-lines with 2048 samples per A-line are collected for a cross-sectional image of a metal plate sample. This metal plate represents defective industry elements with air cavities inside. [Fig. 2(a)](#g002){ref-type="fig"} presents these raw measurements. To reconstruct the sample signal, first we need to estimate the source power spectrum *G*(*k~m~*) from the raw data. We note that the raw measurements, denoted *Ĩ~l~*(*k~m~*), resemble a Gaussian-shaped low-frequency signal modulated by high-frequency components. The former originates from the SLD source power spectrum *G*(*k~m~*), while the high-frequency components are the interference signals represented by the second and third terms in [Eq. (1)](#e01){ref-type="disp-formula"}, together with detector noise. These can be suppressed by a lowpass filter when our objective is to estimate *G*(*k~m~*).

![(a) Raw measurements *Ĩ~l~* (*k~m~*), where each A-line is represented by a different color. (b) The averaged Fourier transform of *Ĩ~l~*(*k~m~*). (c) The estimated source power spectrum. (d) *I*(*k~m~*) after adjustment.](boe-3-4-741-g002){#g002}

To reduce the detector noise, we first average the 800 A-lines we acquired. Because the Fourier transform of a Gaussian function is also Gaussian, the cut-off frequencies of the low-pass filter are chosen as the first pair of zero points in the averaged transform result. We then obtain the estimated source power spectrum, *G*~est~(*k~m~*), from the inverse Fourier transform of the filtered signal. [Fig. 2(b) and (c)](#g002){ref-type="fig"} present the averaged transformation and the spectrum estimation, respectively. Subtracting the latter from the system raw measurements, we have the 800 "adjusted" A-line measurements *I*(*k~m~*) as plotted in [Fig. 2(d)](#g002){ref-type="fig"}. For each A-line, *I*(*k~m~*) can now be deemed a linear combination of *p~s~*(*z~n~*) using [Eq. (3)](#e03){ref-type="disp-formula"}.

Next, we note that *I*(*k~m~*) also has an envelope of a Gaussian function due to *G*(*k~m~*). However, the envelope is not symmetric with respect to zero, which seems to contradict the symmetry of the cosine function in H of [Eq. (5)](#e05){ref-type="disp-formula"}. In fact, this is caused by several factors. Chief among them is the interference between different sample layers *p~s~*(*z~n~*), represented by the third term of [Eq. (1)](#e01){ref-type="disp-formula"}. We can re-write this term as $$\left( {\int_{0}^{\infty}{p_{s}(z)\text{e}^{\text{j}2ka_{s}z}\text{d}z}} \right)\left( {\int_{0}^{\infty}{p_{s}(z^{\prime})\text{e}^{- \text{j}2ka_{s}z^{\prime}}\text{d}z^{\prime}}} \right) = \left| {\int_{0}^{\infty}{p_{s}(z)\text{e}^{\text{j}2ka_{s}z}\text{d}z}} \right|^{2},$$where the values of *p~s~*(*z*) are real. This nonnegative term shifts the measurement envelope towards the positive side with respect to zero. Other factors causing the asymmetric envelope include the detector noise, quantization error in data acquisition, and detector nonlinear responsivity. All of these limit the reconstruction performance of the TV regularization.

At this point, we compare three different methods to reconstruct the sample signal. The first is the conventional FFT method where the measurements for each A-line are resampled before the transformation. The second uses the NUDFT, where the reconstructed signal for one A-line is obtained by multiplying the measurement vector **y** with a matrix Ĥ, with elements *Ĥ~m~*~,~*~n~* = exp (−*jk~m~z~n~*) \[[@r12], [@r13], [@r14]\]. The third uses the TV regularization detailed in the previous section. We choose *α*~1~ = 20 and *α*~2~ = 100. Totally 80 iterations are used for the algorithm defined in [Eq. (9)](#e09){ref-type="disp-formula"}. The cross-sectional sample reconstructions using these three methods are presented in [Fig. 3](#g003){ref-type="fig"}. The horizontal and vertical directions are the A-line (*l*) and the index of the depth (*n*), respectively.

![Reconstruction using (a) the FFT, (b) the NUDFT, and (c) the TV regularization.](boe-3-4-741-g003){#g003}

When we compare the NUDFT with the FFT, the former improves the reconstruction quality, especially in the large depth area marked by the white dash lines. However, with TV regularization, the reconstructed image preserves the sample signal and presents much less noise residue, delivering a better SNR than the other two existing methods. Note that SNR here refers to the difference between the desired sample signal *p~s~*(*z*) and any remaining error in the reconstruction. Because $G(k_{m})p_{r}^{2}$ or the signal from the reference mirror has been subtracted from the system raw measurement *Ĩ*(*k~m~*), this SNR is different from the value obtained by comparing the reconstruction using a mirror sample and the signal by simply detecting the reflection from the reference arm \[[@r9]\]. To have a better visualization of these results, we plot the signals for the A-line at *l* = 150 in [Fig. 4(a)](#g004){ref-type="fig"}, and at the depth position *n* = 425 in (b). The smoothness of the results corresponding to the TV regularization suggests that it minimizes fluctuations in the signal.

![One line of the reconstructed signals using FFT, NUDFT, and TV along (a) the axial direction and (b) the transversal directions.](boe-3-4-741-g004){#g004}

Furthermore, the TV method discussed in this paper focuses on improving OCT system reconstruction performance. The gain demonstrated here comes at the expense of computation speed. In our Matlab code for the FFT, NUDFT, and TV methods, the computation time using the TV method to reconstruct an image consisting of 1000 A-lines with 2048 measurements per A-line is slightly shorter than 35 minutes. For the same image, it takes about half a minute using either the FFT or NUDFT method. We observe roughly the same order of computation time in the following experiments.

In a second experiment, we investigate the effect of the TV regularization to an important SD-OCT system property, namely, sensitivity falloff. We obtain the data for a cross-section of a silver plate and reconstruct using the three methods discussed above. We capture 1000 A-line, each with 2048 measurements. The parameters *α*~1~ and *α*~2~ for the TV regularization are still *α*~1~ = 20 and *α*~2~ = 100. We use 100 iterations to obtain the reconstruction result. [Fig. 5](#g005){ref-type="fig"} presents the reconstructed signals and their zoomed-in views.

![Reconstructions using (a) the FFT, (c) the NUDFT, and (e) the TV regularization. Zoomed-in views of the reconstructions using (b) the FFT, (d) the NUDFT, and (f) the TV regularization.](boe-3-4-741-g005){#g005}

A general trend in [Fig. 5](#g005){ref-type="fig"} is that as the depth value *z~n~* increases, as evident from the slanted lines approaching the left, errors such as the vertical line segment indicated by the arrows in [Fig. 5(b)](#g005){ref-type="fig"} appear in the reconstruction. This is because the wavenumber *k* is inversely proportional to the wavelength *λ*, and therefore a uniformly-spaced *λ* corresponds to fewer measurements in the region where *k* is large. As a result, more reconstruction error exists in a larger depth. The situation is considerably better with NUDFT and TV. If we compare [Fig. 5(c) with (e)](#g005){ref-type="fig"}, we see that the reconstruction using the TV regularization presents much less noise residue than using NUDFT, consistent with the earlier experiment.

Now to study sensitivity falloff, we make use of [Fig. 5(b), (d), and (f)](#g005){ref-type="fig"}, and pick 3 A-lines where the two slanted lines appear in various depths. The lines we selected are at *l* = 9 (for deep in the sample), *l* = 601 (intermediate), and *l* = 930 (shallow). Notice that in these images, we masked out the places far from the two slanted lines. For example, for *l* = 9, we only selected *n* from about 660 to 760; for *l* = 601, roughly for 460 \< *n* \< 560; and for *l* = 930, around 350 \< *n* \< 450. The 3 A-lines for the different methods are presented in [Fig. 6](#g006){ref-type="fig"}. As shown in (c) for the TV regularization, two peaks corresponding to the two layers in the sample can be observed clearly in the reconstruction, whether for deep, intermediate, or shallow regions. In contrast, for FFT and NUDFT, the two peaks corresponding to the intermediate and deep regions are difficult to be identified. We locate the value of *n* for each peak in [Fig. 6(c)](#g006){ref-type="fig"}. Then *p~s~*(*z~n~*) at these locations using the three methods are summarized in [Table 1](#t001){ref-type="table"}. Note that, *n* = 384, 497, and 704 are for the first sample layer, while the other three are for the second layer. As expected, generally *p~s~*(*z~n~*) becomes smaller as *n* increases or the depth is enlarged. For example, for the first layer, as *n* increases from 384 to 704, the normalized *p~s~*(*z~n~*) decreases from −5.60 dB to −21.68 dB in the FFT method, from 0 dB to −16.23 dB in the NUDFT method, and from −1.68 dB to −12.91 dB in the TV method. However, the difference between the minimum and maximum *p~s~*(*z~n~*) varies with the three schemes. The smallest is the 12.91 dB in the TV method, which is 3.32 dB smaller than that for NUDFT (16.23 dB) and 7.71 dB smaller than that for FFT (20.62 dB).

![Zoom-in reconstructed signals at *l* = 9, 601, and 930 using (a) the FFT, (b) the NUDFT, (c) the TV regularization.](boe-3-4-741-g006){#g006}

###### 

Normalized peak values (dB) of the reconstructed signal in [Fig. 6](#g006){ref-type="fig"}.

                                FFT      NUDFT    TV
  ----------------- ----------- -------- -------- --------
  *l* = 930         *n* = 384   −5.60    0        −1.68
  *n* = 410         −3.54       −2.33    0        
  *l* = 601         *n* = 497   −14.59   −3.87    −6.12
  *n* = 525         −1.06       −1.37    −0.60    
  *l* = 9           *n* = 704   −21.68   −16.23   −12.91
  *n* = 733         −9.50       −0.52    −3.05    
  min -- max (dB)   −20.62      −16.23   −12.91   

In [Fig. 7](#g007){ref-type="fig"}, we present the reconstruction results obtained using the three methods for an orange flesh and a pearl sample. The former is located at a small depth position, while the pearl sample has a large depth. We still use 100 iterations in the TV method. The parameters used in the method are *α*~1~ = 20 and *α*~2~ = 100. Comparing the three reconstructions for the orange flesh sample, the TV regularization shows the best reconstruction quality. In the reconstruction using the FFT method for the pearl sample, once again we observe the artifact at a large depth region, which does not exist in the other two reconstruction schemes. For the pearl sample we also observe speckle noise. This is suppressed with TV regularization, because TV acts as a prior to enhance the edges and reduce random noise.

![Reconstructions for an orange flesh and a pearl using (a--b) the FFT, (c--d) the NUDFT, and (e--f) the TV regularization.](boe-3-4-741-g007){#g007}

In the last experiment, we use a fingernail sample to test the three methods. Reconstruction results are presented in [Fig. 8](#g008){ref-type="fig"}. Once again, we observe the TV method improves the reconstruction quality, especially in the large depth area pointed by the green arrow. With all of these studies, we can therefore conclude that TV regularization has less falloff, and shows better edge performance. It is a purely computational approach to improve the SD-OCT system performance. It has the advantage of a lower cost because there is no need for additional hardware.

![Reconstructions for a fingernail sample using (a--b) the FFT, (c--d) the NUDFT, and (e--f) the TV regularization.](boe-3-4-741-g008){#g008}

5. Conclusions {#sec5}
==============

In this work, we use a linear system of equations to model the SD-OCT measurements, and reconstruct the image from nonuniformly sampled signal via solving an inverse imaging problem. With a TV constraint as the regularization, we show that this method provides a better noise suppression and delivers a better contrast compared to the conventional FFT-based and NUDFT-based methods. Furthermore, we can reduce the system's sensitivity falloff using the TV regularization. In addition, by modifying the model parameters, the inverse imaging method can also be used for other OCT systems besides SD-OCT.

This work was supported in part by the University Research Committee of the University of Hong Kong under Project 10208648. The help of Mr. Ray Zhu in providing the biological samples is gratefully acknowledged.
