Fast dispersion encoded full range optical coherence tomography for retinal imaging at 800 nm and 1060 nm

The dispersion mismatch between sample and reference arm in frequency-domain optical coherence tomography (OCT) can be used to iteratively suppress complex conjugate artifacts and thereby increase the imaging range. In this paper, we propose a fast dispersion encoded full range (DEFR) algorithm that detects multiple signal components per iteration. The influence of different dispersion levels on the reconstruction quality is analyzed experimentally using a multilayered scattering phantom and in vivo retinal tomograms at 800 nm. Best results have been achieved with 30 mm SF11, with neglectable resolution decrease due to finite resolution of the spectrometer. Our fast DEFR algorithm achieves an average suppression ratio of 55 dB and typically converges within 5 to 10 iterations. The processing time on non-dedicated hardware was 5 to 10 seconds for tomograms with 512 depth scans and 4096 sampling points per depth scan. Application of DEFR to the more challenging 1060 nm wavelength region is also demonstrated by introducing an additional optical fibre in the sample arm. © 2010 Optical Society of America OCIS codes: (110.4500) Optical coherence tomography; (170.4500) Optical coherence tomography; (300.6300) Fourier transforms; (100.3010) Image reconstruction techniques; (100.3020) Image reconstruction-restoration; (100.5070) Phase retrieval. References and links 1. M. Wojtkowski, A. Kowalczyk, R. Leitgeb, and A. F. Fercher, “Full range complex spectral optical coherence tomography technique in eye imaging,” Opt. Lett. 27, 1415–1417 (2002). http://www.opticsinfobase. org/ol/abstract.cfm?URI=ol-27-16-1415 2. R. A. Leitgeb, C. K. Hitzenberger, A. F. Fercher, and T. Bajraszewski, “Phase-shifting algorithm to achieve highspeed long-depth-range probing by frequency-domain optical coherence tomography,” Opt. Lett. 28, 2201–2203 (2003). http://ol.osa.org/abstract.cfm?URI=ol-28-22-2201 3. P. Targowski, M. Wojtkowski, A. Kowalczyk, T. Bajraszewski, M. Szkulmowski, and W. Gorczynska, “Complex spectral OCT in human eye imaging in vivo,” Opt. Commun. 229, 79–84 (2004). http://dx.doi.org/ doi:10.1016/j.optcom.2003.10.041 #122272 $15.00 USD Received 5 Jan 2010; revised 9 Feb 2010; accepted 12 Feb 2010; published 24 Feb 2010 (C) 2010 OSA 1 March 2010 / Vol. 18, No. 5 / OPTICS EXPRESS 4898 4. P. Targowski, W. Gorczynska, M. Szkulmowski, M.Wojtkowski, and A. Kowalczyk, “Improved complex spectral domain OCT for in vivo eye imaging,” Opt. Commun. 249, 357–362 (2005). http://dx.doi.org/doi: 10.1016/j.optcom.2005.01.016 5. J. Zhang, J. S. Nelson, and Z. P. Chen, “Removal of a mirror image and enhancement of the signal-to-noise ratio in Fourier-domain optical coherence tomography using an electro-optic phase modulator,” Opt. Lett. 30, 147–149 (2005). http://www.opticsinfobase.org/ol/abstract.cfm?URI=ol-30-2-147 6. E. Götzinger, M. Pircher, R. A. Leitgeb, and C. K. Hitzenberger, “High speed full range complex spectral domain optical coherence tomography,” Opt. Express 13, 583–594 (2005). http://www.opticsinfobase.org/ oe/abstract.cfm?URI=oe-13-2-583 7. M. A. Choma, C. Yang, and J. A. Izatt, “Instantaneous quadrature low-coherence interferometry with 3×3 fiber-optic couplers,” Opt. Lett. 28, 2162–2164 (2003). http://ol.osa.org/abstract.cfm?URI= ol-28-22-2162 8. M. V. Sarunic, M. A. Choma, C. H. Yang, and J. A. Izatt, “Instantaneous complex conjugate resolved spectral domain and swept-source OCT using 3x3 fiber couplers,” Opt. Express 13, 957–967 (2005). http://www. opticsinfobase.org/oe/abstract.cfm?URI=oe-13-3-957 9. M. V. Sarunic, B. E. Applegate, and J. A. Izatt, “Real-time quadrature projection complex conjugate resolved Fourier domain optical coherence tomography,” Opt. Lett. 31, 2426–2428 (2006). http://www. opticsinfobase.org/ol/abstract.cfm?URI=ol-31-16-2426 10. S. Yun, G. Tearney, J. de Boer, and B. Bouma, “Removing the depth-degeneracy in optical frequency domain imaging with frequency shifting,” Opt. Express 12, 4822–4828 (2004). http://www.opticsinfobase. org/oe/abstract.cfm?URI=oe-12-20-4822 11. A. M. Davis, M. A. Choma, and J. A. Izatt, “Heterodyne swept-source optical coherence tomography for complete complex conjugate ambiguity removal,” J. Biomed. Opt. 10, 064 005–6 (2005). http://link.aip. org/link/?JBO/10/064005/1 12. A. H. Bachmann, R. A. Leitgeb, and T. Lasser, “Heterodyne Fourier domain optical coherence tomography for full range probing with high axial resolution,” Opt. Express 14, 1487–1496 (2006). http://www. opticsinfobase.org/oe/abstract.cfm?URI=oe-14-4-1487 13. A. H. Bachmann, R. Michaely, T. Lasser, and R. A. Leitgeb, “Dual beam heterodyne Fourier domain optical coherence tomography,” Opt. Express 15, 9254–9266 (2007). http://www.opticsexpress.org/ abstract.cfm?URI=oe-15-15-9254 14. B. J. Vakoc, S. H. Yun, G. J. Tearney, and B. E. Bouma, “Elimination of depth degeneracy in optical frequencydomain imaging through polarization-based optical demodulation,” Opt. Lett. 31, 362–364 (2006). http:// www.opticsinfobase.org/ol/abstract.cfm?URI=ol-31-3-362 15. Y. Yasuno, S. Makita, T. Endo, G. Aoki, M. Itoh, and T. Yatagai, “Simultaneous B-M-mode scanning method for real-time full-range Fourier domain optical coherence tomography,” Appl. Opt. 45, 1861–1865 (2006). http: //ao.osa.org/abstract.cfm?URI=ao-45-8-1861 16. R. K. Wang, “In vivo full range complex Fourier domain optical coherence tomography,” Applied Physics Letters 90, 054 103 (2007). http://link.aip.org/link/?APL/90/054103/1 17. R. A. Leitgeb, R. Michaely, T. Lasser, and S. C. Sekhar, “Complex ambiguity-free Fourier domain optical coherence tomography through transverse scanning,” Opt. Lett. 32, 3453–3455 (2007). http://ol.osa.org/ abstract.cfm?URI=ol-32-23-3453 18. B. Baumann, M. Pircher, E. Götzinger, and C. K. Hitzenberger, “Full range complex spectral domain optical coherence tomography without additional phase shifters,” Opt. Express 15, 13 375–13 387 (2007). http:// www.opticsexpress.org/abstract.cfm?URI=oe-15-20-13375 19. Y. K. Tao, M. Zhao, and J. A. Izatt, “High-speed complex conjugate resolved retinal spectral domain optical coherence tomography using sinusoidal phase modulation,” Opt. Lett. 32, 2918–2920 (2007). http://ol. osa.org/abstract.cfm?URI=ol-32-20-2918 20. S. Makita, T. Fabritius, and Y. Yasuno, “Full-range, high-speed, high-resolution 1-μm spectral-domain optical coherence tomography using BM-scan for volumetric imaging of the human posterior eye,” Opt. Express 16, 8406–8420 (2008). http://www.opticsexpress.org/abstract.cfm?URI=oe-16-12-8406 21. M. Szkulmowski, I. Grulkowski, D. Szlag, A. Szkulmowska, A. Kowalczyk, and M. Wojtkowski, “Flow velocity estimation by complex ambiguity free joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express 17, 14 281–14 297 (2009). http://www.opticsexpress.org/abstract.cfm?URI=cfm?URI=ol-32-23-3453 18. B. Baumann, M. Pircher, E. Götzinger, and C. K. Hitzenberger, “Full range complex spectral domain optical coherence tomography without additional phase shifters,” Opt. Express 15, 13 375–13 387 (2007). http:// www.opticsexpress.org/abstract.cfm?URI=oe-15-20-13375 19. Y. K. Tao, M. Zhao, and J. A. Izatt, “High-speed complex conjugate resolved retinal spectral domain optical coherence tomography using sinusoidal phase modulation,” Opt. Lett. 32, 2918–2920 (2007). http://ol. osa.org/abstract.cfm?URI=ol-32-20-2918 20. S. Makita, T. Fabritius, and Y. Yasuno, “Full-range, high-speed, high-resolution 1-μm spectral-domain optical coherence tomography using BM-scan for volumetric imaging of the human posterior eye,” Opt. Express 16, 8406–8420 (2008). http://www.opticsexpress.org/abstract.cfm?URI=oe-16-12-8406 21. M. Szkulmowski, I. Grulkowski, D. Szlag, A. Szkulmowska, A. Kowalczyk, and M. Wojtkowski, “Flow velocity estimation by complex ambiguity free joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express 17, 14 281–14 297 (2009). http://www.opticsexpress.org/abstract.cfm?URI= oe-17-16-14281 22. K. Wang, Z. Ding, Y. Zeng, J. Meng, and M. Chen, “Sinusoidal B-M method based spectral domain optical coherence tomography for the elimination of complex-conjugate artifact,” Opt. Express 17, 16 820–16 833 (2009). http://www.opticsexpress.org/abstract.cfm?URI=oe-17-19-16820 23. B. Cense, N. A. Nassif, T. C. Chen, M. C. Pierce, S. H. Yun, B. H. Park, B. E. Bouma, G. J. Tearney, and J. F. de Boer, “Ultrahigh-resolution high-speed retinal imaging using spectral-domain optical coherence tomography,” Opt. Express 12, 2435–2447 (2004). http://www.opticsinfobase.org/oe/abstract.cfm?URI= oe-12-11-2435 24. D. L. Marks, A. L. Oldenburg, J. J. Reynolds, and S. A. Boppart, “Autofocus algorithm for dispersion correc#122272 $15.00 USD Received 5 Jan 2010; revised 9 Feb 2010; accepted 12 Feb 2010; published 24 Feb 2010 (C) 2010 OSA 1 March 2010 / Vol. 18, No. 5 / OPTICS EXPRESS 4899 tion in optical coherence tomography,” Appl. Opt. 42, 3038–3046 (2003). http://www.opticsinfobase. org/ao/abstract.cfm?URI=ao-42-16-3038 25. M. Wojtkowski, V. J. Srinivasan, T. H. Ko, J. G. Fujimoto, A. Kowalczyk, and J. S. Duker, “Ultrahighresolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12, 2404–2422 (2004). http://www.opticsinfobase.org/oe/abstract.cfm? URI=oe-12-11-2404 26. K. E. OHara and M. Hacker, “Method to suppress artifacts in frequency-domain optical coherence tomograghy,” US7330270 (2008). 27. B. Hofer, B. Považay, B. Hermann, A. Unterhuber, G. Matz, and W. Drexler, “Dispersion encoded full range frequency domain optical coherence tomography,” Opt. Express 17, 7–24 (2009). http://www. opticsexpress.org/abstract.cfm?URI=oe-17-1-7 28. B. Hofer, B. Povazay, B. Hermann, A. Unterhuber, G. Matz, and W. Drexler, “Dispersion encoded full range frequency domain OCT,” in Proc. SPIE, J. G. Fujimoto, J. A. Izatt, and V. V. Tuchin, eds., 7168, 71 681I (2009). http://link.aip.org/link/?PSI/7168/71681I/1 29. S. Witte, M. Baclayon, E. J. Peterman, R. F. T


Introduction
Most of the current optical coherence tomography (OCT) systems are based on frequencydomain (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.
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][24][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 recon-struction 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.

Physical signal model
Using a simple physical model with constant bulk dispersion mismatch between reference and sample arm beams [25,27], the spectral interference signalS(ω) = S(ω) −S(ω) (withS(ω) denoting the background signal and ω denoting optical frequency) can be expressed as . (1) 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 componentsS 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 inS d (ω) andS 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 signalS(ω) 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 nonuniformly 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 resampleS(ω) on an equidistant frequency grid ω u in order to use fast Fourier transform (FFT) for algorithm implementation.

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.

Description
Physical model Digital signal . The digital signal model corresponding to Eq. (1) and augmented by measurement noise w ∈ R N then reads Here, t ∈ C 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 1 optimization problem using the dispersive basis V = Φ Φ Φ Ψ Ψ Ψ can be phrased from Eq.
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]).

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].
2. Determine the dictionary vector v n which is maximally collinear with the residual, 3. Detect further signal components with sparse representation in the dictionary {v i }, multi peak detector (MPD): 4. Update the estimate of the corresponding signal components and the residual: 5. Increment i. If i < I and r i 2 > ε f 2 , then continue iterating with Step 2; otherwise proceed with Step 6.
6. Obtain the final estimate for the desired full range complex time domain FD-OCT signal t and optionally add the residual:t

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 estimatet. 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 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,n i | are higher than the threshold 1−δ γ, the multi peak detector (MPD) decides to retain these elements inc i , whereas all other elements of the detected signal component signalc i are set to zero, since they are likely to be corrupted by artifact terms of the stronger elements. 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 componentsc i are added to the output signalt [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 endc i is transformed via FFT, dispersion is applied, the corresponding mirror signals are generated by forwarding only the real part of the dispersed spectrum ofc 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 2 < ε f 2 . 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.

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 = − N 2 ,..., N 2 − 1 and dispersion coefficientsâ 2 = 1,â 3 = 7 · 10 −5 . Dispersive phase is then determined as φ (ω u ) = 1 N−1 (â 2 u 2 +â 3 u 3 ) and the resulting dispersive broadening corresponds to a dispersion mismatch of ∼ 20 mm SF11. With two pulses, g 1 = 1.7955 exp( j π/5) positioned at N/4 and 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 0 = f after standard processing step. According to Eq. (3) with Eq. (5) and Eq. (11) we find that in the first iteration , c 1 is a superposition of the desired signal d = d 1 + d 2 , 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,n 1 with magnitude ∼0 dB at n 1 = 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 signalt [blue line in Fig. 2(b)]. After update of the residual Eq. (8) the smaller second peak at position n 1 = 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 byt = Ψ Ψ Ψ Hf withf = d 1 + d 2 + 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 valuest 512 andt 1536 are shown in Table 2 and compared to the values obtained with DEFR reconstruction or just dispersion compensation, i.e., c 1 , 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. 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.

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

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 ∈ R N . To this end detected signals from a whole tomogram were averaged, i.e., b = 1 N x ∑ N x x=1 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 2 . 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 functiong(u) ∈ − N 2 , N 2 − 1 was used for linear resampling of the signalss 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).

Dispersion estimation
Signals from a plane mirror allowed to directly access the dispersion parameters by phase analysis. To this end a phase signal ϕ ϕ ϕ ∈ R N was extracted from the interference signal f by unwrapping arg f + j H {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 1 N−1 (â 0 +â 1 u +â 2 u 2 +â 3 u 3 ) with u = − N 2 ,..., N 2 − 1 ,â 0 being a phase offset andâ 1 corresponding to OPD τ was fitted to the concave part of the phase data. Finally the dispersive phase was constructed as φ (ω u ) = 1 N−1 (â 2 u 2 +â 3 u 3 ). 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.

Results and discussion
The magnitude of the dispersion compensated signal c 1 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. 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.  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.

Alias positions
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    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).
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  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 . 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.
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 inhomogeneous 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 depthscan 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. 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 O (2IN(1 + log N)), whereas standard processing with numeric dispersion compensation obeys O (N(1 + log N)). 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.

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