1. Introduction

Most of the current optical coherence tomography (OCT) systems are based on frequency-domain (FD) detection of the interference spectrum. A drawback of FD-OCT is the inability to retain the complete signal phase since only the real part of a generally complex spectrum is detected. Thus, with standard FD-OCT systems only half of the available depth range can be used for imaging. Otherwise complex conjugate mirror terms overlap with sample structures.

To overcome this problem, several full range techniques using multiple measurements of the interference spectrum have been proposed. A large variety of approaches has been individually tested so far, e.g., stepping phase with piezo-mounted reference mirror [1–4] or using electro-optic phase modulators [5,6], instantaneous detection with 3×3 fibre couplers [7–9], frequency shifting [10–13], or polarization diversity [14]. Recently, techniques that enforce a sequential phase shift between consecutive A-lines have been evaluated [15–22]; here, phase instabilities due to subject motion can introduce imaging artifacts or cause reduced mirror term suppression.

Dispersion mismatch between reference and sample arm fields in the interferometer leads to a wavelength dependent phase shift of the spectral interference fringes. Numerical dispersion compensation [23–25] serves to compensate this phase shift. It was noted that dispersion causes a blurring of complex artifacts [26] and artifact suppression can be increased when dispersion compensation is combined with complex FD-OCT [20].

Recently, a method for full range reconstruction from a single measurement of the interference spectrum using a standard FD-OCT system has been proposed [27]. This iterative reconstruction method uses the dispersive spreading of mirror terms to suppress complex conjugate terms in individual depth scans. The technique is inherently phase stable but computational complexity is higher. Later it was demonstrated [28] that increasing the level of dispersion results in better in vivo images and that cancelation of multiple signal components will increase convergence speed as compared to the basic scheme [27]. After one DEFR iteration, conjugate artifacts were already found to be reduced on in vitro samples [29]. Recently DEFR has also been employed for spectroscopic measurements in single and multilayered non-scattering phantoms [30].

In this manuscript we develop and extensively evaluate a fast version of the iterative DEFR scheme. The algorithm is explained in detail after briefly reviewing a discrete signal model for FD-OCT. The proposed fast DEFR algorithm detects multiple signal components in each iteration. A crucial step here is the determination of the detection threshold in order to obtain stable reconstruction and to avoid destruction of tomogram data. We show that reconstruction quality can be retained with the correct choice of the threshold. Our novel implementation requires only two fast Fourier transforms (FFTs) per iteration, further reducing computational complexity. The basic behavior of the algorithm is demonstrated via simulation. The dependency on the amount of dispersion mismatch is illustrated experimentally using spectrometer based FD-OCT. The depth dependency is evaluated and the influence of limiting effects due to finite spectrometer resolution on DEFR reconstruction quality is discussed. We evaluate the algorithm on in vivo retinal tomograms acquired with a 800 nm FD-OCT setup using various levels of dispersion and different algorithm parameters. This allows us to draw conclusions on how to choose parameters for best imaging results with DEFR. Finally we demonstrate the applicability of DEFR to data acquired with a 1060 nm system.

2. Background

2.1. Physical signal model

Using a simple physical model with constant bulk dispersion mismatch between reference and sample arm beams [25, 27], the spectral interference signal S̃(ω) = S(ω) - S̅(ω) (with S̅(ω) denoting the background signal and ω denoting optical frequency) can be expressed as

Here, the nth reflective layer in the sample arm has optical delay τ_n and intensity I_n(ω). Assuming reflective sites positioned at positive and negative delays τ_n the desired full range signal S̃_d(ω) generally overlaps with the conjugate mirror components S̃_m(ω) identified in Eq. (1).

Common bulk dispersion results in the frequency-dependent phase term ϕ(ω). Using dispersion coefficients ã_i the phase can be written as

From Eq. (1) it follows that the dispersive phase occurs with opposite sign in S̃_d(ω) and S̃_m(ω). This phase difference allows to iteratively identify and separate the mirror components from the full range signal on data from a single depth-scan [27].

To obtain the interference signal the background signal S̅(ω) is subtracted from the detected spectrum S(ω). Here, S̅(ω) is either determined by separate measurements with blocked sample and reference beams respectively or estimated from tomogram data using the assumption of a small and heterogeneous sample arm signal. Generally the spectrum is sampled at non-uniformly distributed frequency points ω_p, i.e., with pixel index p ∈ {1,…,N} of the CCD linescan camera employed in the spectrometer. It is convenient to resample S̃(ω) on an equidistant frequency grid ω_u in order to use fast Fourier transform (FFT) for algorithm implementation.

2.2. Digital signal and problem formulation

Next we formulate a digital signal model employing vector notation [27] that adequately represents the physical model summarized in the preceding paragraph. While finite resolution effects of the spectrometer are not addressed, the influence of a sampling function that deviates from an ideal Dirac pulse will be discussed in Section 5. We assume that N samples of the detected spectrum S(ω) are measured. The corresponding digital signals are summarized in Table 1.

Table 1. Vector notation used for algorithm formulation.

View Table | View all tables in this article

We define the discrete Fourier transform (DFT) matrix Φ := [ψ ₁, ψ ₂,…, ψ_N] ∈ ℂ^N×N and the N×N dispersive measurement matrix Φ= diag(ϕ). The digital signal model corresponding to Eq. (1) and augmented by measurement noise w ∈ ℝ^N then reads

Here, t ∈ ℂ^N represents the full range time domain data which is assumed to be a sparse vector, i.e., only K ≪ N elements are non-zero. The non-zero elements correspond to the desired signal components in the order of reflecting sites from the physical model above (we note that the Fourier transform of the intensities I_n(ω) is usually very narrow).

DEFR reconstruction aims to recover the sparse vector t from the measured interference signal f given the dispersive phase ϕ [27]. To this end, an ℓ ₁ optimization problem using the dispersive basis V = ΦΨ can be phrased from Eq. (3) as follows

Thus, DEFR reconstruction searches for a sparse representation of the measurement vector f in the dictionary {v _i} consisting of the column vectors of the dispersive basis V (cf [31]).

3. Algorithm

3.1. Fast DEFR reconstruction

In [27] a greedy matching pursuit (MP) algorithm [31] was adapted in order to iteratively find locally optimum solutions of Eq. (4). The algorithm requires specification of the maximum number I ≥ K of iterations and of the fraction ε of the energy of f that may remain in the residual, which usually depends on the noise level and modeling errors (cf. [31]).

We propose to improve the algorithm of [27] by including detection of multiple signal components in each iteration. Thus the convergence speed can be rapidly increased, depending on the amount of dispersion inserted in the system. This also dramatically reduces the computational load required for DEFR reconstruction while retaining reconstruction quality as will be verified experimentally. The key element of the fast algorithm is a multi peak detector. This additional step requires a predetermined threshold γ(ϕ,b) ∈ [0,1] and a stability parameter δ ∈ [0,1].

FAST MP ALGORITHM FOR DEFR FD-OCT SIGNAL RECONSTRUCTION

3.2. Implementation

The block diagram of the fast DEFR algorithm is shown in Fig. 1. The iterative procedure is initialized by loading the residual r with the interference spectrum f, and clearing the depth scan estimate t̂. Element-wise multiplication of r with the conjugate dispersive phase vector ϕ ^* and transformation to the spatial domain via inverse FFT equals the standard processing with numeric dispersion compensation used in FD-OCT. Regarding the DEFR scheme, the resulting intermediate spatial signal c _i = (c _i,1…c_i,N)^T provides all N normalized inner products c_i,n = ⟨v_n, r _i-1)⟩/∥v_n∥, n = 1,…,N in Eq. (5). The strongest complex signal component c_i,ni. is detected at the position n_i of the element of |c _i| with maximum value as with the basic algorithm [27].

Further signal components are found by searching for elements of c _i with absolute value higher than contributions due to the blurred mirror artifacts [cf. Eq. (6)]. Thus if the elements c_i,n, n = 1,…,N normalized by the magnitude of the strongest signal component |c_i,ni| are higher than the threshold 1 - δγ, the multi peak detector (MPD) decides to retain these elements in c̃_i whereas all other elements of the detected signal component signal c̃_i are set to zero, since they are likely to be corrupted by artifact terms of the stronger elements.

 

Fig. 1. Block diagram of the fast dispersion encoded full range (DEFR) algorithm. The iterative procedure is indicated by iteration index i. The multi peak detector is denoted as MPD; further notations: measured interference signal f, dispersive phase term ϕ, fast Fourier transform (FFT), inverse FFT (IFFT), residual signal r, intermediate spatial signal c; and finally the complex full range tomogram line t̂ ∈ ℂ_N.

Download Full Size | PDF

The threshold is precalculated for a given stability parameter δ ∈ [0,1], the threshold factor γ(ϕ, b) depends on dispersion and spectral shape, i.e., corresponding to b, and is determined as follows:

The frequency domain representation ψ _N/4 of a unit-pulse at position N/4 is used as a test signal in Eq. (10). The ratio of the maximum magnitude of the dispersed and the dispersion compensated test signal provides a measure for the disparity of desired signal components and mirror artifacts. High dispersion provides a larger spreading and thus suppression of mirror artifacts; correspondingly, the denominator in Eq. (10) will increase, resulting in a lower threshold.

The detected signal components c̃_i are added to the output signal t̂ [cf Eq. (7)]. Update of the residual r _i according to Eq. (8) is performed by directly subtracting the frequency domain representation of the detected signal components along with their corresponding mirror artifacts. To this end c̃_i is transformed via FFT, dispersion is applied, the corresponding mirror signals are generated by forwarding only the real part of the dispersed spectrum of c̃_i and finally a factor 2 is included to preserve signal energy. The remaining new residual r _i+1 has high energy signal parts and their conjugate terms removed, lower energy parts become visible and can be detected by continuing the iterative scheme. The algorithm stops when either the maximum number of iterations I is reached or the residual contains only noise, i.e., ∥r _i∥² < ε∥f∥². As a last step the remaining final residual may be added to the output signal. This optional operation will be considered later when only a small number of iterations is performed, i.e., the final residual contains significant signal energy.

3.3. Simulation

To demonstrate the algorithms’ principle we explore a simple simulated sample consisting of two reflective sites, e.g. originating from a glass plate. Whereas this simulation does not include limiting effects of the spectrometer or detection electronics it allows visualization of the iterative behavior and influence of the stability parameter.

Consider N = 2048 sample points, discrete frequency index $u=[-\frac{N}{2},\dots ,\frac{N}{2}-1]$ and dispersion coefficients â ₂ = 1, â ₃ = 7 · 10^-5. Dispersive phase is then determined as $\varphi \left(\omega \right)=\frac{1}{N-1}\left({\hat{a}}_2{u}^2+{\hat{a}}_3{u}^3\right)$ and the resulting dispersive broadening corresponds to a dispersion mismatch of ~ 20 mm SF11. With two pulses, g ₁ = 1.7955exp(jπ/5) positioned at N/4 and g ₂ = 0.17955 exp(jπ23/19) at 3N/4 the detected interference spectrum equals [cf. Eq. (3)]

where w is zero-mean white Gaussian noise with variance 0.025 corresponding to a noise floor of ~-65 dB and b is a Tukey window. A good approximation of the laser spectrum from the experiments is achieved by choosing 0.89 as ratio of cosine tapered to constant sections for b. The pulse positions are conjugate, thus mirror terms are maximally overlapping with desired signal components.

For DEFR reconstruction the maximum number of iterations was set to I = 12 and stopping parameter was chosen as ε = 0. The resulting spatial domain signals for δ = 0.75 are shown in Fig. 2, e.g. Fig. 2(a) depicts the magnitude of the initial residual r ₀ = f after standard processing step. According to Eq. (3) with Eq. (5) and Eq. (11) we find that in the first iteration c ₁ = Ψ _H(d + Φ ^* Φ ^* d ^* + Φ ^* w), i.e., c ₁ is a superposition of the desired signal d = d ₁ + d ₂, the suppressed and spread conjugate terms and the dispersed noise. With white noise the application of the phase term Φ ^* will not change the appearance of the noise floor. However, autocorrelation artifacts and remaining DC terms in real world OCT signals will become broadened and suppressed. With large dispersion the maximum element c _1,n1 with magnitude ~0 dB at n ₁ = 512 is barely disturbed by overlapping mirrored portions of the second pulse [cf. green signal in Fig. 2(a)]. It can be extracted from the residual, and r can be cleaned from the corresponding strongest conjugate artifact [27].

The multi peak detector accelerates this concept by extracting all signal components above the threshold [indicated by dashed line Fig. 2(a)] within a single iteration. It is based on the assumption that contributions of conjugate mirror terms in |c _i| are all well below the absolute maximum signal multiplied by the predetermined threshold [-10.8 dB as determined by Eq. (6) with Eq. (10)], thus not only is the maximum signal extracted but also further high energy signal parts, without the need to individually correlate with the dispersive basis V. According to Eq. (7) the detected signal components are added to the reconstructed signal t̂ [blue line in Fig. 2(b)]. After update of the residual Eq. (8) the smaller second peak at position n ₁ = 1536 bunches out and can be extracted in the second iteration since the disturbing mirror artifacts of the first peak have already been canceled. The iterative procedure continues by adding weaker signal parts to the output and after 12 iterations the dispersion compensated final residual [black line in Fig. 2(a)] is also added [cf. Eq. (9)].

Without dispersion compensation and assuming ideal complex detection the original signal is given by t̃ = Ψ ^H̃ f with f̃ = d ₁ + d ₂ + w [cf. red signal in Fig. 2(b)]. Some noise is rejected in the reconstructed signal since the algorithm is only able to reconstruct the noise parts which can be expressed by V. The values t̃₅₁₂ and t̃₁₅₃₆ are shown in Table 2 and compared to the values obtained with DEFR reconstruction or just dispersion compensation, i.e., c ₁, at the two peak positions which are exactly conjugate. Successively the iterative scheme allows reconstruction of the original complex signal values with absolute error smaller than 0.5%. Furthermore DEFR reconstruction error is about ten times smaller compared to the values extracted after standard processing.

 

Fig. 2. Simulation of 2 reflective sites positioned at n ₁ = 512 and n ₂ = 1536. (a) Solid green line: initial residual r ₀ = f after dispersion compensation and IFFT; dashed magenta line: threshold level of multi peak detector for first iteration, δ = 0.75; solid black line: residual after 12 iterations. (b) Solid red line: original signal t̃ without dispersion and complex detection; solid blue line: DEFR reconstruction of complex full range signal t̂ (Media 1).

Download Full Size | PDF

Table 2. Complex signal values and average absolute errors for simulated signals at two different positions.

View Table | View all tables in this article

Convergence behavior and signals for all 12 iterations and different stability parameters δ = 0.5, 0.75 and 0.95 can be viewed online (Media1); results are summarized in Table 3. It follows that with increasing threshold the number of iterations required to reach the noise-floor in the reconstructed signal decreases. However the reconstruction stability is reduced which is expressed by a residual level above the noise floor. As a side effect of the iterative reconstruction scheme it was found that the noise floor was reduced with average noise reduction of about 3 dB. DEFR scrambles the noise floor when it tries to reconstruct it using the dispersive basis V, i.e., noise components from successive iterations add up incoherently which results in a reduction of the final noise-floor.

Table 3. Simulation results for different values of stability parameter.

View Table | View all tables in this article

4. Methods

4.1. Experimental setups

Two fibre-based OCT systems operating in the wavelength regions of 800 nm and 1060 nm respectively have been used in this study [27, 32, 33]. For in vivo imaging the systems were interfaced to a modified fundus camera patient module (OCT-2, Carl Zeiss Meditec, CA, USA) with a collimator, two galvanometric mirrors in close pair configuration and variable focusing optics in front of the subject’s eye optimized for the respective wavelength region. Maximum sensitivity close to ZD was 96 dB for 800 nm and 97 dB for 1060 nm.

The 800 nm setup was based on a 90/10 fibre optic beam-splitter (Ipitek Inc., CA, USA) and used a Ti:Sapphire laser with 140 nm FWHM bandwidth (Integral OCT, Femtolasers GmbH, Austria). In the detection arm a refractive spectrometer was used that utilized a 2048 pixel silicon CCD-camera (AVIIVA M2 CL2014-BAO, Atmel, CA, USA) and imaging was performed at a line rate of 20 kHz. Sample and reference arm fibres (SM 650, Fibrecore Ltd., UK) were looped through polarization control paddles, polarization was adjusted for maximal fringe visibility. A gold-coated hollow corner cube reflector mounted on a motorized translation-stage was used in the reference arm. A reflective neutral-density (ND) filter-wheel allowed to adjust reference arm power.

To obtain test signals for algorithm evaluation at 800 nm the fundus camera was replaced by a scanning microscope employing an objective with working distance 10 mm and numerical aperture ~0.1 was used. The optics of the scanning microscope have been pre-compensated by using a BK-7 prism pair in the reference arm. Dispersion mismatch was achieved by inserting SF-11 prism pairs of different length in the reference arm. The resulting thickness of dispersive material was either 5, 10, 20, 40, 60, 80 or 150 mm, respectively. For zero-delay (ZD) adjustment additional air path was introduced after the collimator at the scanning microscope entrance. To measure depth dependent signal parameters a plane mirror and a ND filter were used in the sample arm. Power was adjusted equally for reference and sample signal to about 1/4 of the CCD-saturation level which resulted in maximum fringe contrast.

The scanning microscope was also used to image a multi-layered phantom including scattering spheres. The scattering phantom was composed of a stack of 8 coverslips (each ~ 200μm thick) that were separated by oil/water suspension (~ 15μm). The surface of the phantom was tilted by about 8° to reject specular reflections from detection. Power at the objective was 2.1 mW and maximum sensitivity close to ZD was 94 dB. At the entrance of the spectrometer a bandwidth of 128 nm FWHM remained. The same sample position was then imaged with different levels of SF-11 in the reference arm to allow direct comparison of the algorithms’ performance for effective cancelation of mirror terms from scattering objects.

For in vivo retinal imaging the power at the cornea was 950 μW. Typical 19° (5.5 mm) wide horizontal scans at the fovea and 38° (11 mm) wide-field scans spanning through fovea and optic nerve head were conducted with the different amounts of SF-11 from 5 mm to 150 mm. The dispersion introduced by the vitreous of the subjects eye was pre-compensated by additionally placing a 25 mm deep glass-cuvette filled with water in the reference arm.

The 1060 nm setup was based on a 80/20 beam-splitter (Ipitek Inc., CA, USA) and used an amplified spontaneous emission source with 72 nm FWHM bandwidth (NP-Photonics Inc., Arizona, USA). In the detection arm a spectrometer with all reflective components in a Czerny-Turner geometry was used. The spectrometer utilized a 1024 pixel InGaAs camera (SU-LDH 1024, SUI-Goodrich, NJ, USA) and imaging was performed at a line rate of 47 kHz. For in vivo retinal imaging the power at the cornea was 1.9 mW otherwise imaging protocols as described above have been performed. To achieve a sufficient amount of dispersive broadening 160 mm SF-11 was introduced in the reference arm. As alternative method to obtain a substantial dispersion mismatch at 1060 nm the fundus camera was connected via an additional 2.1m patch-cord. For ZD adjustment additional air path corresponding to the optic length of the patch-cord was introduced after the collimator in the reference arm. Thus, the free-space length in the reference arm was about 3 m. With good alignment reference arm power remained constant, also when adjusting the position of the reference mirror, which was mounted on a translation stage.

4.2. Signal processing

We briefly summarize the additional signal processing tasks that have been performed as described previously [27].

As a first step the background signal b was estimated from the measured raw signal s ∈ ℝ^N. To this end detected signals from a whole tomogram were averaged, i.e., $\mathbf{b}=\frac{1}{N_x}{\sum }_{x=1}^{N_x}{\mathbf{s}}_x,$ with x being a spatial index that determines the transversal position of an A-scan on the tomogram consisting of N_x lines. When measuring a plane mirror the background signal was determined as b = b _s + b _r - b _d, with b _s, b _r and b _d being the background signals when blocking sample arm, reference arm and both interferometer arms respectively.

Autocorrelation terms and fixed pattern noise are removed from s via background subtraction, i.e., s̃ = s - b ⟨b, s⟩/ ∥b∥². Incorporation of the normalized inner product accounts for energy fluctuations of the background signal over successive tomogram lines, which are probably caused by timing issues within the CCD-electronics of the employed camera.

The spectrometer frequency sampling points ω_p have been calibrated independently using a free-space interferometer [27]. The resulting discrete mapping function $\tilde{g}\left(u\right)\in [-\frac{N}{2},\frac{N}{2}-1]$ was used for linear resampling of the signals s̃ to uniform discrete frequencies u = [-N/2.…N/2 - 1]. Up-sampling by a factor 2 and sinc-interpolation filter with length 2N was used prior to the resampling to prevent additional aliasing components close to end of depth (EOD) and to reduce resampling errors of the fast linear procedure. The resampled signals can be decimated again to obtain the interference signals f with N samples, alternatively in this study the DEFR algorithm was executed on the upsampled signals consisting of 2N samples which further reduced partial aliasing components when imaging near EOD.

The processing steps outlined above have been implemented in Matlab and all images in the results section are based on this. For real-time view a simplified processing scheme was implemented in the LabView based data acquisition software. Thus background subtraction did not account for energy fluctuations and signals have not been upsampled for real-time view. Processing time for real-time view tomograms with 512×2048 pixels was ~0.5 seconds for standard processing and ~1.8 for DEFR with 5 iterations (Intel 2GHz dual core processor, 8 GB RAM).

4.3. Dispersion estimation

Signals from a plane mirror allowed to directly access the dispersion parameters by phase analysis. To this end a phase signal φ ∈ ℝ^N was extracted from the interference signal f by unwrapping arg {f + j𝓗{f}}. Depending on OPD and the amount of dispersion mismatch φ has concave and convex parts. According to Eq. (1) and Eq. (2) dispersion parameters were estimated from concave parts of the phase signal. Thus the reference mirror was first positioned such that most part of the phase signal obeys a concave function and the OPD between reference and sample mirror was negative. Then a third order polynomial $\frac{1}{N-1}\left({\hat{a}}_0+{\hat{a}}_{1}u+{\hat{a}}_2{u}^2+{\hat{a}}_3{u}^3\right)\mathrm{with}u=[-\frac{N}{2},\dots ,\frac{N}{2}-1],$ being a phase offset and â ₁ corresponding to OPD τ was fitted to the concave part of the phase data. Finally the dispersive phase was constructed as $\varphi \left({\omega }_u\right)=\frac{1}{N-1}\left({\hat{a}}_u{u}^2+{\hat{a}}_3{u}^3\right)$. Table 4 shows the extracted dispersion coefficients dependent on the amount of dispersion mismatch introduced for up-sampled data (N = 4096). For measurements of the scattering phantom and in vivo imaging the extracted coefficients have been used as an initial guess and have been optimized using a sequential quadratic programming algorithm and information entropy as a sharpness metric. The resulting coefficients have also been included in Table 4.

Table 4. Dispersion coefficients employed in the different experiments.

View Table | View all tables in this article

5. Results and discussion

The magnitude of the dispersion compensated signal c ₁ from a plane mirror positioned at -100 μm is shown in Fig. 3 for various levels of dispersion (SF11 of length d). For the 800 nm spectrometer the end of depth (EOD) was found to be 1373 μm in air, correspondingly the full depth range was 2746 μm. After dispersion compensation the peak at -100 μm was clearly visible for all levels d. The conjugate artifact terms appear broadened and suppressed as wide arcs with maximum 18 to 28 dB below the signal peaks. Signals have been normalized by their maximum value to allow easier comparison of different dispersion levels. For 5 mm SF11 the absolute width of the blurred mirror term was about EOD/4, for 10 mm ~EOD/2, 20 mm ~EOD, 40 mm ~2EOD and further increasing for 60 mm, 80 mm and 150 mm. The maximum magnitude of the artifact was decreasing by about 10 dB from 5 mm to 80 mm, however it was increasing for 150 mm SF11. To visualize this behavior the inset of Fig. 3 depicts the dispersion diversity DD defined as ratio of maximum conjugate signal part to maximum signal peak (blue curve). To simulate DD dependent on length d of the dispersive material we use:

which is shown as green curve and is strictly decreasing as opposed to the measured points. The reason for this was the increasing signal loss for larger dispersion due to the finite resolution of the spectrometer. This effect is similar to the signal decay observed for increasing spatial frequency because a high amount of dispersion mismatch results in a large local phase change on parts of the interference spectrum corresponding to a high instantaneous frequency. The proposed threshold (magenta curve) as determined by Eq. (10) better follows the measured points since it penalizes high dispersion levels by using the ratio of maximum conjugate signal part to dispersed signal peak. For a composite signal which was the superposition of individual signals from the plane mirror positioned between -1200 μm to -100 μm with 100 μm steps, dispersion diversity becomes drastically smaller which was the reason for inclusion of the stability factor δ in the multi-peak detector. Also, any uncertainty or error in the employed dispersion parameters would cause residual dispersion and reduced dispersion diversity. In addition to undesired image blurring due to the uncompensated dispersion part, the effect on the reconstruction quality of the algorithm will be similar as with a reduced dispersion mismatch.

 

Fig. 3. Signal of single reflection after dispersion compensation for various levels of dispersion. The signal peak is located at -100 μm and the overlapping broad arc corresponds to the double dispersed conjugate mirror artifact. Inset: Dispersion diversity obtained with different amounts of SF11 [cf. Eq. (12)].

Download Full Size | PDF

Spectral envelopes for 3 different positions of the mirror (-600 μm, 200 μm and 1000 μm) and 4 different levels of dispersion (no dispersion mismatch, SF11 5 mm, 40 mm and 150 mm) are illustrated in Fig. 4. With balanced dispersion the spectral envelope does not much depend on the mirror position [Fig. 4(a)]. A dispersion mismatch of 5 mm SF11 [Fig. 4(b)] resulted in a small change of the envelope shape as compared to Fig. 4(a) at high frequencies, depth dependency was still moderate. When increasing dispersion to 40 mm SF11 [Fig. 4(c)] spectral shape becomes depth dependant. This effect is even more prominent with dispersion mismatch due to 150 mm SF11 and in addition the spectral width was visibly reduced [Fig. 4(d)]. These results demonstrate the signal deterioration expected by increasing dispersion mismatch due to the finite spectral resolution of the spectrometer.

Fig 5 depicts depth dependancy of signal parameters; signal decay is shown in Fig. 5(a), conjugate artifact suppression in Fig. 5(b) and axial resolution in Fig. 5(c). For dispersion mismatch up to 20 mm SF11 signal decay and axial resolution are almost identical within a depth range of ±1200 μm to values obtained with no dispersion. Only close to EOD deviations become observable, also indicated by outliers due to alias positions. At ±1000μm a signal loss of ~-8 dB could be quantified. At the same positions axial resolution was ~3μm and resolution at zero delay position was 2.85μm in air. With 40 mm SF11 a 1 dB reduction of achievable intensity and a decrease of ~0.1μm in axial resolution was found. With increasing dispersion the maximum intensity was further reduced and axial resolution decreased as well. This is a direct consequence of the spectral contraction observed in Fig. 4. DEFR algorithm was applied on the signals with δ = 0.5, I = 10 and the final residual was added. The suppression ratio was defined as magnitude of the signal peak divided by the magnitude of the signal at the conjugate peak position according to [17]. Shape of suppression ratio was related to signal decay shape and within ±1200 μm was higher than 40 dB. Average values for suppression ratio are 59.3 dB for 5 mm SF11, 59.3 dB (10 mm), 58.6 dB (20 mm), 55.3 dB (40 mm), 55.8 dB (60 mm), 56.3 dB (80 mm) and 55.4 dB (150 mm) respectively. The observed decrease of suppression ratio around ZD is caused because the current implementation of the DEFR algorithm only provides local optimum solutions of Eq. (4).

 

Fig. 4. Depth dependency of spectral envelope for various levels of dispersion and different mirror positions. (a) No dispersion mismatch; (b) SF11 5 mm; (c) SF11 40 mm; (d) SF11 150 mm. Dotted line: mirror at -600 μm; solid line: mirror at 200 μm; dashed line: mirror at 1000 μm.

Download Full Size | PDF

Further results for algorithm evaluation are illustrated in Fig. 6. Individual signals from the plane mirror positioned between -1200 μm to -100 μm with 100 μm steps have been combined in a composite signal [black line in Fig. 6(a) for 40 mm SF11]. Conjugate artifacts add up with highest parts just 14 dB below the maximum signal [cf. also inset of Fig. 3]. After 1 iteration (red line) with δ = 0.5 the artifacts are already reduced to about -30 dB and 2 iterations (cyan line) result in further decrease of ~6 dB. After 5 iterations (green line) most artifact terms are already rejected below -40 dB and after 30 iterations (yellow line) the algorithm almost converged. Conjugate artifact suppression (CAS) was calculated as ratio of mean value of residual terms in the positive range to maximum signal in negative range. For algorithm test the positive range composite signal (mirror between 100 μm to 1200 μm) was also used and CAS values presented have been averaged for positive and negative range composite signal. Fig. 6(b) depicts the dependancy of CAS on stability paramater δ and the number of iterations I for 40 mm SF11. With δ = 0 in each iteration only the maximum signal component is retained, which corresponds to the slow convergence curve of the basic scheme presented in [27]. When only 1 iteration is applied the achievable maximum CAS was about 35 dB when choosing δ = 0.82. This specific setting has also been investigated in [29]. The algorithm then operates in a region with nearly instable reconstruction properties, as can be seen by the increasing and oscillating CAS values for δ > 0.82. Choosing δ between 0.5 and 0.75 allows stable signal reconstruction, high CAS values ~50dB and fast convergence within about 5 iterations. Similar behavior was also observed for the other amounts of SF11 tested. For illustration of the dependency on dispersion Fig. 6(c) depicts CAS after 1 iteration in response to different stability parameters δ. With increasing δ up to the point of critical stable reconstruction (δ ≈ 0.8) CAS becomes larger for all dispersion levels. Highest suppression is achieved with 80 mm SF11 since dispersion diversity was also maximal for 80 mm SF11 [cf. inset of Fig. 3]. Weaker suppression is obtained with 150mm SF11 due to the limiting effects of the spectrometer. With DEFR operating with stable reconstruction at δ= 0.5 fastest convergence within 5 iterations was found for 40 mm SF1 [cf. Fig. 6(d)]. From the convergence curves Fig. 6(d) it was observed that increasing dispersion levels result in a higher residual artifact error. This was caused by the increasing depth dependancy of the spectral envelopes [cf. Fig. 4], currently not included in the signal model Eq. (3). The increased residual level found with 20 mm SF11 resulted from aliased terms close to EOD.

 

Fig. 5. (a) Signal loss, (b) artifact suppression ratio of DEFR algorithm, (c) resolution with different amounts of SF11 and data from a plane mirror. Parameters for DEFR: δ = 0.5, I = 10, final residual added.

Download Full Size | PDF

Tomograms from the multilayered phantom with scattering spheres are shown in Fig. 7. For the DEFR algorithm, 7 iterations and δ = 0.5 were chosen. The final residual was not added thereby saving one transformation step, as there was no visible difference when adding the residual since residual level was already below the noise floor. With no dispersion mismatch between sample and reference arm, true signal components and overlapping conjugate artifact terms cannot be distinguished [Fig. 7(a)]. With 5 mm SF11 artifact terms appear blurred after dispersion compensation [Fig. 7(b)] and further blurring resulted from increased dispersion mismatch by 80 mm SF11 which heavily disturbs the image [Fig. 7(c)]. DEFR is able to reconstruct the full range tomogram, zero delay position in the middle of the tomogram is indicated by a magenta arrow. Slight residual artifact terms close to ZD which are visible with 5 mm SF11 [Fig. 7(d)] vanish with 80 mm SF11 [Fig. 7(e)], however 80 mm exhibits an increased noise floor since the limiting effects of the spectrometer are currently not included in the signal model and the algorithm only achieves a local optimum solution of Eq. (4). Vertical zooms for different dispersion levels are depicted in Fig. 7(f), demonstrating the decreasing resolution and signal intensity for increasing dispersion. A dispersion level of 20 mm to 40 mm SF11 was found to be a good compromise for imaging of scattering structures. Residual artifacts were barely visible while resolution and sensitivity were still retained. Processing time was 0.7 seconds for dispersion compensation and 8.3 seconds for 7 iterations of DEFR (512 depth scans and 4096 spectral sampling points per depth scan).

 

Fig. 6. Evaluation of DEFR scheme with composite signal. (a) Reconstructed composite signal for 40 mm SF11 and δ = 0.5, (b) conjugate artifact suppression (CAS) for 40 mm SF11, (c) CAS after 1 iteration dependent on stability parameter δ, (d) convergence curves with δ = 0.5 for different levels of dispersion.

Download Full Size | PDF

Applicability for in vivo imaging with different algorithm parameters was investigated with results illustrated in Fig. 8 for 20 mm SF11. The left upper row tomogram in Fig. 8 was obtained with standard processing (dispersion compensation), processing time was 0.7 seconds (512 depth scans and 4096 spectral sampling points per depth scan). The final residual was added in all tomograms and with only 1 iteration (1.6 seconds, upper row in Fig. 8) this causes a disturbing residual noise floor even when choosing a critical stable threshold [29]. When operating in the instable regime (δ= 0.9 and 1) black vertical stripes appear and parts of the tomogram become destroyed. Choosing a stable threshold δ= 0.5 (lower row in Fig. 8) allows successive decreasing of the artifact terms with each iteration. After 6 iterations (7.2 seconds) the tomogram is free of artifacts, after 8 iterations (9.5 seconds) the noise floor is reduced and the algorithm has converged. Each iteration costs about 1.2 seconds computational time (on oversampled data) and adds about 6 dB of dynamic range. When not adding the final residual for real time preview, acceptable image quality was achieved after 4 iterations (8 FFTs) since already all the high energy signal components have been reconstructed corresponding to a dynamic range of ~24 dB. Total dynamic range shown (as indicated by the colorbar in Fig. 7) was 40 dB in all in vivo images and identical to the images of the scattering phantom .

 

Fig. 7. Demonstration of fast DEFR algorithm on data from scattering phantom. (a) Dispersion balanced; (b,c) dispersion due to 5 mm and 80 mm SF11 respectively has been numerically compensated; (d,e) DEFR reconstruction after 7 iterations with stability parameter δ=0.5 for 5 mm and 80 mm SF11 respectively; each tomogram is 2048×512 pixels corresponding to 2746 μm×445 μm in air. (f) Vertically zoomed tomogram portions for balanced dispersion, 5 mm, 40 mm, 80 mm and 150 mm SF11 respectively; zoomed region is indicated by magenta colored frame for 80 mm SF11.

Download Full Size | PDF

Similar results have been found for different levels of dispersion and a qualitative comparison is depicted in Fig. 9. The same subject was imaged several times at approximately the same position starting with dispersion matched measurement and aiming for best quality tomograms. Slight tilt between measurements and subjects eyestrain influenced signal strength which makes it difficult to draw direct conclusions as opposed to the measurements of the scattering phantom with identical setting, however, the global trend of reduced signal and resolution for dispersion levels above 40 mm SF11 could be verified. After 1 iteration (upper row in Fig. 9) and choosing a close to critical threshold δ= 0.7 artifact terms are reduced with increasing dispersion which is consistent with Fig. 6(c). Much better results have been achieved with 5 iterations and a stable threshold δ = 0.5 (lower row in Fig. 9). Choosing a dispersion mismatch of 20 mm to 40 mm SF11 (corresponding to a dispersion diversity of ~25 dB with a peak spreading of EOD to 2 EOD) and a stability parameter of δ= 0.5 to 0.7 allows high quality full range image reconstruction within 4 to 10 iterations while still retaining axial resolution and signal strength.

 

Fig. 8. Behavior of fast DEFR algorithm under various parameters for in vivo retinal imaging at 800 nm with dispersion mismatch due to 20 mm SF11. Each tomogram is 2048×512 pixels corresponding to 2746 μm×5.5 mm in air.

Download Full Size | PDF

 

Fig. 9. Behavior of fast DEFR algorithm for in vivo retinal imaging at 800 nm with dispersion mismatch due to different amounts of SF11. Each tomogram is 2048×512 pixels corresponding to 2746 μm×5.5 mm in air. (a) DEFR with 1 iteration and stability parameter δ=0.7; (b) DEFR with 5 iterations and δ=0.5.

Download Full Size | PDF

 

Fig. 10. Application of fast DEFR algorithm for wide field in vivo retinal imaging at 800 nm. Each tomogram is 2048×2048 pixels corresponding to 2746 μm×11 mm in air.

Download Full Size | PDF

With DEFR, phase stable reconstruction can be achieved since no additional measurements are required as DEFR operates on data from individual single depth-scans. Thus DEFR is well suited for wide-field retinal imaging, without any constraints on transversal sampling. This is demonstrated in Fig. 10(middle) for a tomogram spanning 19° (11 mm) directly through the fovea and optic nerve head. When imaging deep structures at the optic nerve head full range imaging is particularly useful since with half depth range structures are easily overlapping with artifact terms, as can be seen in the left tomogram of Fig. 10 for matched dispersion (processing time 2.8 seconds, slightly different spatial location). DEFR was performed with 5 iterations and by choosing δ=0.5, no final residual was added and processing time was 19.5 seconds on oversampled data (4096 samples per depth scan, 2048 depth scans). For comparison the right tomogram in Fig. 10 was obtained after numeric dispersion compensation [same raw data as Fig. 10(middle)] and is heavily distorted by blurred conjugate mirror artifacts. The inhomo-geneous appearance of the noise floor with DEFR was caused by non-adaptive operation of the algorithm on individual depth-scans, i.e. depth-scans with smaller signal energy obey faster convergence. Future improvements in implementation of the algorithm could include depth-scan adaptive choice of residual energy ε and iteration number I to achieve a homogeneous convergence level.

Application of the DEFR scheme for in vivo imaging at 1060 nm is depicted in Fig. 11. EOD was determined to be 3543 μm, thus full range was 7086μm in air, signal drop at 1750μm was ~10 dB and axial resolution ~10 μm. The wide-field tomogram after standard processing with dispersion mismatch of 160 mm SF11 exhibits strong distortions by mirror artifacts [cf. Fig. 11(a)]. Dispersive broadening due to SF11 at 1060 nm was much weaker than at 800 mm, dispersion diversity was 19 dB and absolute width of blurred conjugate artifact from a plane mirror was about EOD/2. Due to the weaker dispersion, residual artifact terms were still visible after DEFR with 10 iterations, δ = 0.5 [Fig. 11(b)]. Increased dispersion was obtained by introducing a 2.1 m patch-cord in the sample arm, dispersion diversity increased to 22 dB and artifact broadening was ~1.25 EOD. Thus residual artifacts could be suppressed further [Fig. 11(c)]. Processing time for 1024 depth scans with 2048 spectral sampling points per depth scan (oversampled data) was 0.6 seconds for dispersion compensation and 11.5 seconds for 10 iterations with DEFR.

 

Fig. 11. Application of fast DEFR algorithm for wide field in vivo retinal imaging at 1060 nm. (a) Dispersion due to 160 mm SF11 has been compensated numerically; (b) DEFR after 10 iterations and δ=0.5; (c) dispersion due to 2.1 m patch-cord, DEFR after 10 iterations and δ=0.5. Each tomogram is 1024×1024 pixels corresponding to 7086 μm×11 mm in air.

Download Full Size | PDF

A comparison of DEFR with standard processing for various imaging depths at 1060 nm with the 2.1 m patch-cord is finally shown in Fig. 12. Dispersion compensation is exhibited in the upper row tomograms of Fig. 12, the lower row tomograms show results after DEFR with 10 iterations and δ = 0.5. Alias terms are visible in the tomograms on the left side of Fig. 12. Figure 12 demonstrates that DEFR is also feasible at 1060 nm despite the smaller bandwidth (~70 nm) and reduced number of spectral sampling points (1024 pixel camera) as compared to the 800 nm system (130 nm, 2048 pixel). Although artifact terms are greatly reduced some reconstruction errors remain since the current algorithm implementation provides local optimum solutions of Eq. (4). Future improvements could target more sophisticated algorithmic steps (e.g. including position information of strong signal components and reiterating to exhibit global minimization, incorporation of additional parameters such as effective spectrum or inclusion of a suitable digital spectrometer model [34] to cope with limiting effects of the spectrometer).

The asymptotic computational complexity of the algorithm for a single depth scan is about 𝓞(2IN(1+logN)), whereas standard processing with numeric dispersion compensation obeys 𝓞(N(1 + logN)). Since the number of iterations required for algorithm convergence is only about 5 to 10 when using multiple peak detection, DEFR has an asymptotic computational complexity that is about 10 to 20 times larger than standard processing with numeric dispersion compensation.

 

Fig. 12. Comparison of fast DEFR algorithm with standard processing for retinal imaging at 1060 nm. Dispersion due to 2.1 m patch-cord; numerical dispersion compensation (upper row), DEFR with 10 iterations and δ=0.5 (lower row). Each tomogram is 1152×1024 pixels corresponding to 7972 μm×5.5 mm in air.

Download Full Size | PDF

6. Conclusion

Using detection of multiple signal components the fast dispersion encoded full range (DEFR) algorithm allows rapid iterative reconstruction of complex conjugate free FD-OCT signals on single depth-scans. Since DEFR does not require any constraints on transverse scanning patterns and works with individual depth scans it is inherently phase stable. It is based on the spreading of mirror terms due to bulk dispersion mismatch between sample and reference arm beams. After a single iteration of DEFR residual artifacts heavily disturb retinal tomo-grams even when choosing a low detection threshold where the algorithm is close to instability. Choosing a conservative threshold for detection, a few iterations with DEFR results in tomo-grams cleaned from mirror terms. At 800 nm best results regarding reconstruction quality and convergence speed have been achieved with dispersion levels of 20 mm to 40 mm SF11, corresponding to a spreading of individual signal components over the full depth range (by about 100 to 200% of EOD). Thereby the limiting effects due to the finite resolution of the 800 nm spectrometer was not dominant and dispersion diversity was already high (~25 dB). Thus after DEFR reconstruction, residual artifacts are small which resulted in good image quality. Average conjugate artifact suppression was about 55 dB, convergence was achieved within 5 to 10 iterations and processing times are in the range of 5 to 10 seconds per tomogram without using dedicated signal processing hardware. Sufficient dispersion mismatch can be simply introduced by a fibre-length difference as demonstrated at 1060 nm by inclusion of a 2.1 m patch-cord in the sample arm. In vivo full range retinal imaging using DEFR was successfully demonstrated at 800 nm and 1060 nm.