Dispersion encoded full range (DEFR) optical coherence tomography (OCT) has become highly attractive as it is a simple way to increase the measurement range of OCT systems. Full range OCT is especially favorable as it does not only increase the measurement range but also shifts the highest sensitivity into the center of the measurement range. While the early versions of DEFR were highly computational expensive, new versions reduce the number of necessary Fourier transforms. Recently it has been shown that a GPU based algorithm can perform DEFR with more than 20,000 A-lines per second. We present a new version of the DEFR algorithm that requires only one Fourier transform per A-scan and uses convolution in z-space instead of multiplication in k-space, therefore reducing the computational effort considerably. While dispersion encoding has so far only been used to suppress mirror artifacts, we show that, with dispersion encoding and only one more Fourier transform, autocorrelation terms can be removed likewise. Since very high values of dispersion reduce the effective measurement range in dispersion encoded OCT, we present an estimate for a sufficient amount of dispersion for a successful image recovery, which is depending on the thickness of the scattering layers. Furthermore, we demonstrate the usability of ZnSe as a new dispersive material with a very high dispersion and describe a simple method to extract the dispersive phase from the measurement of a single reflex of a glass surface. Using a standard consumer PC, an artifact-free recovery of 1000 – 2000 A-scans per second with 2048 depth values including autocorrelation removal was achieved. The dynamic range (sensitivity) is not reduced and the suppression ratio of mirror artifacts and autocorrelation signals is more than 50dB using ZnSe.
© 2012 Optical Society of America
Optical coherence tomography (OCT) is a noninvasive, contactless, in vivo optical imaging technique which provides cross-sectional or three-dimensional images of transparent to highly scattering samples, for example biological tissue, with a high axial resolution of 1 – 15μm. OCT has found its main application in ophthalmology, but many other areas in biomedical research [1, 2] and nondestructive testing are emerging. During the 20 years of development, the speed of OCT has increased dramatically from about 30 A-scans per second in the early work  to more than 20 million A-scans per second using FDML lasers . A breakthrough in speed was achieved with the swap from time domain OCT (TD OCT) to Fourier domain OCT (FD OCT) eliminating the need for a fast mechanical change of the reference arm length, with the additional benefit of an improved signal to noise ratio (SNR) [5–7]. In TD OCT, the axial resolution is limited by the spectral width of the light source (including other components of the optical system) and the measurement range is limited by the stroke of the mirror in the reference arm, which are mainly independent from each other. In FD OCT, on the contrary, axial resolution and measurement range depend on the width or resolution, respectively, of the acquired wavelength spectrum and therefore are not independent from each other. Because some applications demand both higher resolution and greater measurement range, the requirements for the FD OCT spectrometer system increase and cannot be fulfilled in all cases.
The ratio of depth range to resolution is limited in classical FD OCT by the number of spectral points . While in spectrometer based FD OCT systems, called spectral domain OCT (SD OCT), this number is limited by the pixel number of the line scan camera, in swept source OCT (SS OCT) this is limited by the ratio of the sampling frequency of the analog-to-digital converter (ADC) to the repetition frequency (neglecting the limited duty cycle) of the swept source. As both values cannot simply be increased there is a demand for other methods to increase this ratio. Furthermore, both kinds of OCT systems show a sensitivity roll off to the end of the measurement range [8, 9] having therefore a reduced sensitivity in the mostly interesting center of the image. Consequently, several methods to increase the measurement range have been designed. Full range OCT eliminates the complex conjugate artifacts to effectively double the measurement range around the zero delay position (ZDP), which is in most systems the position of highest sensitivity. The removal of the complex conjugate artifacts is often performed by acquiring a second spectrum at almost the same beam position and introducing a phase shift between both A-scans. This can be done by special interferometers providing quadrature signals [10–12] or by introducing a phase shift between successive A-scans . The latter can be done with a phase modulator in the reference arm [14, 15] or by changing the length of the sample arm. The length of the sample arm can easily be changed by pivoting the rotation of the deflecting scanning mirror (BM-scanning [16–19]). In these systems dense and, in the case of moving samples, fast sampling is needed for high efficiency of complex conjugate removal. While two adjacent scans allow the reduction of the mirror images, better results are achieved by taking up to five A-scans [20, 21]. Furthermore, high phase stability is essential for the algorithm. In contrast to the previously described systems, dispersion encoded full range (DEFR) OCT uses only one spectrum to calculate the full range depth signal [22–26]. While it has been demonstrated that DEFR removes the complex conjugate data very efficient in many samples, there is a more philosophical problem with the algorithm that from only N (real) spectral data points N complex valued depth values can be retrieved. As we will show later, the solution for this dilemma is that the algorithm performs only well if many of the depth values are zero or small.
The basic idea of DEFR OCT is to exploit a dispersion mismatch between the two interferometer arms, caused by a dispersive material in one arm, which leads to a broadening of signal peaks in z-space. This can be corrected by numerical dispersion compensation in k-space before the inverse Fourier transform, thereby restoring the true signal components and broadening their complex conjugate mirror artifacts even more. A peak detector reveals the true signal components and their mirror artifacts can be removed consequently. The previously described versions of the DEFR algorithm introduce one iterative step into the common OCT processing, which comprises two steps. In a first step the spectral data are corrected for detector background then linearized in k-space and spectrally shaped. In case of dispersion mismatch between the sample and reference arm a corresponding phase shift is introduced to restore the resolution [27, 28]. In the second step of the usual FD OCT processing the data are then Fourier transformed and displayed using only the amplitudes and logarithmic scaling. The DEFR algorithm is introduced iteratively in between these two processing steps by identifying the highest signals in the depth data and removing them from the spectral data. This takes two or even four Fourier transforms in each iteration step. While in the early versions only one component was removed in each step, what could result in as many steps as sample data points, the recent versions remove many components in one step. A high dispersion between sample and reference arm ensures a large amplitude ratio between the true signals and the complex conjugate artifacts thereby reducing the number of necessary iteration steps. The iteration can be stopped if all the complex conjugate artifacts are below the noise level. This grouping of data reduces the processing time by about two orders of magnitude to the early versions of the DEFR algorithm. Witte et al.  suggested for a system with high dispersion only one iteration, limiting the number of Fourier transforms to only three, achieving more than 40dB of mirror artifact suppression, which is sufficient for the described system with a dynamic range of 45dB. Recently the DEFR algorithm was implemented on a GPU using the parallel computing capabilities of these devices. By this, the processing time could be reduced by approximately two orders of magnitude again .
The main processing steps of the DEFR algorithm, proposed by Hofer et al. , are depicted in Fig. 1. The iterative DEFR algorithm requires two Fourier transforms for the subtraction of a single signal component. That is in each iteration, two Fourier transforms are needed to calculate from z-space to the complex conjugate spectrum in z-space again after application of a phase shift in k-space:
2. Theory and DEFR implementation
2.1. OCT signal with dispersion and autocorrelation terms
To understand the processing steps of our proposed DEFR algorithm, it is helpful to initially derive a detailed signal model of the OCT signal with dispersion and autocorrelation terms. We will discuss the contributions of these signals to the spectrum and exploit their different properties regarding the influence of the dispersion. This model is based on the general model derived in [23, 30], but includes the crucial terms introduced by dispersion. To evaluate the influence of an artificial dispersion mismatch, we consider a Michelson interferometer with dispersive material in the reference arm. The electric field of the incident polychromatic electromagnetic wave with the spectral amplitude of the electric field s(k) can be written as:31]: Eq. (7). The first component is not frequency dependent and is obtained using Eq. (7) and Eq. (9) as follows: Eq. (10) and thus causes no dispersion, and the dispersive phase term: Eq. (6)) can be rewritten as: Eq. (20) using the relation |x|2 = xx*, the different contributions to the power spectrum can be identified as follows: Equation (21), term (a) describes the DC component, determined by the power spectrum of the light source and the reflection and transmission coefficients. Equation (21), term (b) and Eq. (21), term (c) describe the cross-correlation components that means the ’true’ signal we want to retrieve from the spectrum. Equation (21), term (d) are autocorrelation signals, which are usually considered as artifacts in conventional half range FD OCT systems. As the reflection coefficients of the sample rS,l are usually very low compared to that one of the reference mirror rM, the autocorrelation terms are negligible in many samples. Nevertheless, they can have a disturbing impact on images of highly scattering samples as we will demonstrate later. Note that the autocorrelation signals are not corrupted by the artificial dispersion, which will later enable us to remove them from the spectrum. Numerical dispersion compensation is achieved by multiplying Eq. (21) with the compensating phase term e−iϕ(k): Eq. (22): Eq. (23) consists of the desired complex full range signal Eq. (23), term (b) as well as of its complex conjugate mirror components spread by convolution with the inverse Fourier transform of twice the inverse dispersive phase function ℱ−1(e−i2ϕ(k)) Eq. (23), term (c) and DC and autocorrelation components, Eq. (23), term (a) and Eq. (23), term (d), respectively, which are also spread by the numerical dispersion compensation. In order to remove autocorrelation artifacts from the spectrum, we need to consider both the power spectrum with and without numerical dispersion compensation, so we also need to calculate the inverse Fourier transform of the original, non dispersion compensated spectrum Eq. (21): Eq. (23) and Eq. (25), we can now develop an implementation of the DEFR algorithm that enables us to do both a direct subtraction of the complex conjugate mirror components of a true signal in z-space and to perform an effective autocorrelation suppression. Moreover, Eq. (23) and Eq. (25) mathematically prove that a subtraction of signal components directly in z-space using the corresponding dispersive phase functions equals the subtraction indirectly via k-space as proposed in the former implementations of the DEFR algorithm. The details of our approach will be explained in the following chapter using a block diagram and pseudocode.
2.2. Description of the simple DEFR algorithm without autocorrelation removal
Compared with the proposed DEFR algorithm from Hofer et al. [23, 26], we implement a very efficient method to remove the complex conjugate mirror components caused by each true signal. Our improved implementation of the DEFR algorithm without autocorrelation removal can be derived directly from Eq. (23) when neglecting DC and autocorrelation terms, Eq. (23), term (a) and Eq. (23), term (d), respectively. Although these components are broadened by the dispersion compensation, which leads to a convolution with the inverse Fourier transform of the inverse dispersive phase function ℱ−1(e−iϕ(k)) (z) in z-space, resulting in a lower amplitude according to Parseval’s theorem:Eq. (23), term (a) ( ) is removed via background subtraction, the other terms in Eq. (23), term (a) (scaling with the reflection coefficients of the sample) and Eq. (23), term (d) (autocorrelation terms) cannot be distinguished from the wanted cross-correlation signals Eq. (23), term (b). Nevertheless, we will start with a simple implementation of our proposed DEFR algorithm neglecting DC and autocorrelation terms coming from the sample.
Figure 2 shows the block diagram of our proposed DEFR algorithm without autocorrelation removal concentrating on the iteration procedure itself, which is our basic approach of subtracting the complex conjugate mirror components of a true signal directly in z-space. In other implementations of the DEFR algorithm the subtraction of the complex conjugate mirror components of a true signal is performed by transforming from z-space to the complex conjugate spectrum in z-space again after applying a phase shift k-space [23, 26, 29]. In Eq. (23), term (c) the contribution of the complex conjugate mirror components of a true signal found in Eq. (23), term (b) can be directly identified and thus removed from the spectrum by a simple subtraction. Step (1) comprises the calculation of the corresponding dispersive phase function via inverse Fourier transform. This calculation has to be done only once and the results can be saved as kind of calibration for each setup with a certain fixed dispersion. We will present a simple method of calculating the dispersive phase ϕ(k) from a single measurement in Sect. 3.2. In step (2) the inverse Fourier transform of the background corrected, in frequency resampled, dispersion compensated signal is calculated ( , note that the index 1 is used to differentiate from a second spectrum in our advanced implementation of the DEFR algorithm; i. . .iteration index). The output for the true signal components is denoted . In frequency domain OCT one inverse Fourier transform is anyway necessary for each A-scan so there is no additional (inverse) Fourier transform required for our proposed algorithm. The iterative process consists of step (3) and step (4). A peak detector finds the strongest signal component (highest amplitude) and provides its value and position in the spectrum in z-space , respectively. When comparing Eq. (23), term (b) and Eq. (23), term (c), we can distinguish between the true signal component and its complex conjugate mirror component spread by twice the inverse dispersive phase. According to Parseval’s theorem, the energy of a signal in z-space remains constant if it is multiplied with the inverse dispersive phase function e−iϕ(k) in k-space before applying the inverse Fourier transform as the magnitude of this function equals one (|e−iϕ(k)| = 1). As a narrow peak in Eq. (23), term (b) is broadened by twice the inverse dispersive phase in Eq. (23), term (c) and the signal energy in Eq. (23), term (b) and Eq. (23), term (c) is the same, the amplitude in Eq. (23), term (c) must be smaller by a factor if the peak is one channel wide in Eq. (23), term (b) and 2D channels wide in Eq. (23), term (c), assuming a nearly rectangular shape of the inverse Fourier transform of the dispersive phase function, which is realistic as we will show later (note that the energy of a signal is proportional to the square of the amplitude). The true signal component , found by the peak detector, is added to the output and both the true signal component and its corresponding complex conjugate mirror components in are subtracted (directly in z-space) using and . The iteration stops if a maximum number of iterations is reached or if the amplitude of the true peak falls below a predefined threshold, which results in only noise remaining in the spectrum (i.e. consider the complex conjugate mirror components of the true signal, if they are below the noise level their removal will not improve the image quality and the iteration can be stopped). The iteration process does not need any additional Fourier transforms, therefore the computational complexity is very low – one multiplication and one subtraction for each array element of in each iteration. Moreover, considering only that part of the dispersive phase function that significantly differs from zero (or whose amplitude, when multiplied with a strong peak, is significantly above the noise level) (for multiplication and subtraction), further reduces the number of computations. Optionally, the residual spectrum (remaining after the iteration was stopped) can be added to the output in order to maintain weak signals.
2.3. Description of the advanced DEFR algorithm with autocorrelation removal
After having described a simple implementation of the DEFR algorithm, we now want to present our advanced approach that effectively eliminates autocorrelation artifacts from the spectrum with only minor additional computational effort. The algorithm proceeds similar to that described in Sect. 2.2, but additionally exploits the properties of the autocorrelation signals in Eq. (23), term (d) and Eq. (25), term (d). Consequently, we can reconstruct both the complex full range signal and the autocorrelation signal. While other implementations of the DEFR algorithm do not consider autocorrelation signals or simply set the signals near the zero delay position (which contain in most cases DC and autocorrelation terms) to zero [23, 29], thus removing any true signal in this region likewise and causing artifacts, our algorithm effectively removes them from the spectrum without disturbing true signal components in this region. This is very beneficial as autocorrelation removal is essential in some samples with highly scattering layers as we will show later.
Figure 3 shows the block diagram of our advanced DEFR algorithm including autocorrelation detection and removal. Step (1) again comprises the calculation of the corresponding dispersive phase functions via inverse Fourier transform. Besides the already in Sect. 2.2 used phase function e−i2ϕ(kn), we now also need to calculate the inverse Fourier transform of the phase functions e±iϕ(kn) (compare Eq. (23), term (a) and Eq. (23), term (d) as well as Eq. (25), term (b) and Eq. (25), term (c)). In step (2) the inverse Fourier transform of both the dispersion compensated (upper path) and the non dispersion compensated (lower path) signal is calculated ( and respectively). To reduce the computational effort only one half of the spectrum for autocorrelation suppression is needed as this spectrum from a real valued signal is complex conjugate (the magnitude spectrum is symmetric around the zero delay position) and the other half does not contain any additional information. The outputs for the true signal components and autocorrelation components are denoted and , respectively (i. . .iteration index). The iterative process now consists of step (3) to step (5). A peak detector finds the strongest signal components in both spectra and provides their values , and positions , in the spectra in z-space respectively. In addition to the removal of complex conjugate mirror components of a true signal described in Sect. 2.2, the autocorrelation removal exploits the different properties of DC and autocorrelation signals in the two spectra and . When comparing Eq. (23), term (b) and Eq. (25), term (b) it is obvious that a true signal component found in is spread by the dispersion in Eq. (25), term (b) and results in a lower amplitude of the corresponding signal components in . Thus a real signal component can be distinguished from its corresponding contribution in . The distinction of autocorrelation signals from real signals is performed similar. By comparing Eq. (23), term (d) and Eq. (25), term (d) we can detect autocorrelation signals as their amplitude in is higher than the amplitude of the corresponding signal components in that were spread by the numerical dispersion compensation. The important step for autocorrelation suppression is to decide whether the strongest signal component in the dispersion compensated spectrum ( ) is higher than the strongest signal component in the non dispersion compensated spectrum ( ) or vice versa. In case of True, meaning , we can now subtract both the complex conjugate mirror components in and the corresponding signal components in using and . In case of False, that is the strongest signal component is found in the non dispersion compensated spectrum , this component is said to be an autocorrelation artifact. In this case we proceed similar and subtract the corresponding signal components in . In both cases the outputs are updated accordingly. Note that the subtraction of true signal components in and the subtraction of autocorrelation signal components in does neither affect any of the autocorrelation signals in nor any of the true signals in , so that our approach assures the correct reconstruction of the spectrum around the zero delay position even in highly scattering samples which is not the case in all other implementations of the DEFR algorithm. Although the autocorrelation removal requires one additional Fourier transform for the processing of each A-scan, the iterative procedure comprises only two multiplications and two subtractions for each array element of and in each iteration. As the inverse Fourier transform of the dispersive phase functions e±iϕ(kn) contains only half the number of values significantly greater than zero compared to the inverse Fourier transform of twice the inverse dispersive phase function e−i2ϕ(kn), the computational effort for autocorrelation removal can be further reduced by using an appropriate part of these functions. Again, the residual spectrum can be added to the output .
2.4. Implementation of the advanced DEFR algorithm with autocorrelation removal
We adopt the DEFR algorithm, which uses a greedy matching pursuit algorithm for ℓ1 optimization [23, 32], and change the method of subtracting complex conjugate mirror components. Furthermore, we add a method to remove autocorrelation artifacts by exploiting the information provided by the different spectra Eq. (23) and Eq. (25). The discrete inverse Fourier transforms of the dispersion compensated spectrum Eq. (23) and the non dispersion compensated spectrum Equation (25) as well as the necessary phase functions are calculated as follows:
Advanced DEFR Algorithm With Autocorrelation Removal
i = 0, , .
- Determine the signal components with the highest contribution to the spectra in , and :
- Compare and with predefined threshold value T:
if max then proceed with step 6, otherwise proceed with step 4.
- Update the outputs and subtract corresponding mirror artifacts:
- Check if maximum number of iterations M reached:
i = i + 1,
- if i ≤ M then proceed with step 2, otherwise proceed with step 6.
- (Optional) Add remaining spectrum to the output:
These six steps describe the implementation of the new DEFR algorithm. We found another correction helpful in some cases. Because a true signal component always causes a broadened complex conjugate mirror artifact in the spectrum, these two signals will overlap when the signal is near the zero delay position. Therefore, it is necessary to calculate the true signal component from the complex value in the spectrum using the following assumptions: A complex value in the spectrum is the sum of the true signal component and its complex conjugate mirror component spread by the dispersive phase function :Eq. (27) we can retrieve the true signal :
3. Materials and methods
3.1. Experimental OCT setup
The setup is slightly modified from a dual-band OCT system which enables simultaneous imaging in the 0.8μm and 1.3μm wavelength range . In brief, we use only the 800nm band of the dual-band OCT system, which is illuminated by a SuperK Versa Super Continuum Source (Koheras A/S, Denmark) light source. In the used short wavelength band the spectrometer has a central wavelength of λ0 = 791nm with a spectral width of Δλ = 166nm focused on the line scan camera. Compared to the original setup, the detector was replaced by a Basler sprint spL2048–70km line scan camera. The sensitivity was measured to be −93dB with a spectrometer roll-off of approximately 20dB over the entire depth scan range (dispersion balanced). The fiber coupled scanner head includes the Michelson interferometer and two galvanometer scanners for beam steering. Reference and sample arm of the interferometer are nearly dispersion balanced, leading to a nonlinear phase shift with a maximum mismatch less than 7rad. To introduce additional dispersion for DEFR imaging, a glass plate was placed in the reference arm of the interferometer. Several glass plates have been used starting from a single 6.35mm SF6 plate up to two SF6 plates resulting in total thickness of 12.70mm. For even higher dispersion a 6.25mm plate of ZnSe was inserted into the reference arm. ZnSe is transparent from 600nm on, has a refractive index of ≈ 2.5 at 800nm with dispersion approximately five times as high as SF6.
3.2. Dispersion measurement
There are different methods for estimating the dispersive phase. Using the Sellmeier Eq.:23], which has the disadvantage that only a limited number of coefficients can be considered and that it is not sure that the global maximum is found during the iterative optimization. Moreover, many spectra with different parameter settings must be calculated in order to have a reliable set of data used for optimization. We propose a method that is similar to that one proposed by Witte et al. and Hofer et al. [24, 29], but has the crucial advantages that it does not affect the spectrometer calibration and does not require measurements in different depths. If the function for resampling from detector pixel number to frequency is known, then this method is favorable as it uses only one measurement of a strong reflector and does not depend on the method for resampling. We can calculate the dispersive phase for a given amount of dispersive material in the reference arm of the interferometer by simply measuring the spectrum of a strong reflector like a mirror or a glass surface. The spectrum of a glass surface in z-space is shown in Fig. 4(a). Note that there is a strong DC component as the reflection coefficient of the glass surface is comparatively high (compare Eq. (25), term (a)). We filter the dispersion spread peak in the left half-plane (compare Fig. 4(b)) and calculate the corresponding signal in k-space via Fourier transform. The physical model for dispersion measurement can be derived from Eq. (21), term (b) and Eq. (25), term (b). The filtered signal in z-space: Eq. (31) in k-space is retrieved and unwrapped to remove 2π phase jumps: Fig. 4 for a 6.25mm ZnSe glass plate in the system described in Sect. 3.1.
The dispersive phases ϕ(k) for various amounts of dispersive material in the reference arm of the interferometer have been measured using the technique described in Sect. 3.2 and are depicted in Fig. 5(a). In order to suppress noise and artifacts, the average of 50 A-scans and a fourth order polynomial regression of ϕ(k) was used. While averaging of several A-scans improves the quality of the measured dispersive phases, therefore reducing the width of a numerically dispersion compensated peak in z-space, polynomial regression is not absolutely necessary but can remove existing detector inaccuracys. The group velocity dispersions at the center frequency λ0 = 791nm for SF6 and ZnSe, and , respectively, and the third order dispersions , and , respectively, cause the major part of the dispersion. While the ZnSe and SF6 glass plates have a comparable thickness of 6.25mm versus 6.35mm, ZnSe causes a dispersive phase that is more than five times as high as that one introduced by SF6 (compare GVD(ZnSe) and GVD(SF6)), which makes it favorable for systems with limited space or optical path length since it enables the correct reconstruction of images showing a much larger span of echoes. Figure 5(b) shows the inverse Fourier transforms of twice the inverse dispersive phase functions (ℱ−1 (e−i2ϕ(k)) corresponding to the dispersion spread complex conjugate mirror artifact of a true signal component (compare Eq. (23), term (b) and Eq. (23), term (c)). As already described in Sect. 2.2, according to Parseval’s theorem, a wider range (in z-space) of the mirror components due to a larger dispersion leads to a lower amplitude. A sufficient amount of dispersion, for instance 6.25mm ZnSe, leads to an amplitude ratio of a unity-peak compared to the amplitude of the complex conjugate mirror artifacts after numerical dispersion compensation of .
4. Results and discussion
4.1. Necessary amount of dispersion
The presented algorithm treats each signal independently and is not grouping signals as proposed by Hofer et al. [24, 26]. Therefore a high dispersion, leading to a large ratio between a true signal and its complex conjugate mirror artifact, does not accelerate the algorithm. Nevertheless, we found that the algorithm needs a minimum dispersion to sort the signals to the correct side. We noted that in images with only a few sharp reflexes a small amount of dispersion between the sample and the reference arm of the interferometer already leads to the expected images, while images with a larger span of echoes show large artifacts at low dispersion. False mapping can occur if neighboring signals on the false side that are distributed over a certain range add up to a signal higher than any (residual) component on the correct side. To get more insight in the amount of dispersion needed for a successful recovery of the image we made a simulation with synthetic data. We assumed a band of echoes with Rayleigh distributed amplitude  and arbitrary phase over a certain number N of Fourier components. On this data in one side of the full range spectrum in z-space we applied a quadratic phase shift of arbitrary amplitude, which is the largest part in most cases and applied the presented DEFR algorithm. A quadratic phase shift of an amplitude broadens a single peak to a FWHM of D Fourier components. The recovery of the depth signal was regarded as successful if no signal was found on the false side. We found that the error rate increases strongly with the ratio of the number of Fourier components N divided by the broadening D. While the error rate is not only a function of the quotient but is slightly higher at larger N, for a ratio of below 0.7 the error rate was below 5% and for a ratio below 0.5 it was much below 1%. This result, which should hold true if higher order terms contribute essentially to the dispersion, is in good agreement with the observations made during the recovery of different images at variable amount of dispersion. Note that this is only a statistical argument. As the amplitude of the dispersion broadened signal drops only with the square root of the dispersion D (see Eq. (26)), in rare cases a number of Fourier components with similar amplitude and appropriate phase may result in a false mapping of signals. We recommend to use as less dispersion as possible for reconstructing a scattering layer with a certain thickness, because a higher amount of dispersion reduces the measurement depth range due to partial aliasing of the broadened peaks.
4.2. High resolution full range imaging
In the following, we will present images of various samples using different amounts of dispersion and explain their special characteristics. All images are displayed using a dynamic range of 45dB and a Red Hot color scale. The zero delay position is marked with a white dashed line. The upper left image (a) shows the original image disturbed by the artificial dispersion. The upper right image (b) shows the dispersion compensated image, where dispersion compensation was performed by multiplication with the inverse dispersive phase function e−iϕ(k) in k-space prior to inverse Fourier transform. The lower left image (c) shows the image post-processed with the simple implementation of the DEFR algorithm presented in Sect. 2.2. Note that there is no difference in the results using our simple implementation of the DEFR algorithm and the implementation of Hofer et al.  as both implementations are mathematically equal. The lower right image (d) finally shows our advanced implementation of the DEFR algorithm with autocorrelation removal. The iteration for a single A-scan was finished if a maximum number of iterations was reached (here: 250 iterations) or if the remaining complex conjugate mirror components of the true signals were below the noise level and consequently not causing artifacts if not removed before the residual spectrum was added to the reconstructed A-scan. As the latter condition mostly ensues first, the processing speed benefits of a larger amount of dispersion (for very high dispersion the approach proposed by Witte et al. , using three Fourier transforms, is very efficient, but does not enable autocorrelation removal). If an exact signal recovery is necessary, our advanced DEFR algorithm is to the best of our knowledge the first implementation that provides full range OCT in combination with autocorrelation removal.
Figure 6 shows the images of a nanoparticle target used for resolution estimation  that was imaged using 6.35mm SF6. The axial scanning range is 3.89mm with in air and 2048 Pixel in z-direction. Figure 6(a) shows the original image that is symmetric around the zero delay position as the amplitude spectrum of the real valued signal from the detector is an even function. The effect of the artificial dispersion is clearly visible. Figure 6(b) shows that numerical dispersion compensation leads to a sharp reconstruction of the true image components and spreads the complex conjugate mirror components even more. The sharp structures in the image allow a correct reconstruction with little dispersion. Due to the tilted surface autocorrelation artifacts are not very severe (compare Fig. 6(c) and Fig. 6(d) – autocorrelation signals are hardly visible in Fig. 6(c) and removed in Fig. 6(d)). Figure 6 shows that both our simple and advanced DEFR algorithm successfully reconstruct images with sharp structures using only a small amount of dispersion, whereas images with a larger span of echoes require a larger amount of dispersion for correct reconstruction.
Figure 7 shows the images of a capillary originally used for Doppler measurements. The inner diameter is d = 300μm and the intralipid inside the capillary has a concentration c = 1% leading to a refractive index n = 1.34. The image does intentionally not cross the zero delay position as a false mapping of signals is not visible in this position. Note that the aspect ratio of the image is stretched in order to accentuate the errors due to a false signal mapping over the width of the homogeneous back scattering intralipid layer. We used this sample to experimentally investigate the relation between the amount of dispersion and the width of a correct reconstructed homogenous back scattering layer. The horizontal white dashed lines indicate the maximal width the intralipid layer may have for a proper reconstruction, which is around 200μm. This is in very good agreement with the estimations made in Sect. 4.1. The scattering layer with a width of 200μm has an optical length of 200μm × 1.34 = 268μm. According to , this results in an wide signal in z-space (number of Fourier components). Using the width of a dispersion spread peak D = 215 Pixel (12.70 mm SF6), the ratio is below 0.7, therefore implying an error ratio below 5%. As opposed to this, the scattering layer in the center of the image, having a width of around 300μm, leads to a ratio of , which is well above 0.7 and therefore indicates that there will be lots of errors in the image. These estimations hold true when comparing the number of errors outside and inside the vertical white dashed lines in the image.
Figure 8 shows the same sample used in Fig. 7, but imaged with a much higher dispersion using 6.25mm ZnSe. As can be seen from Fig. 8(c) and Fig. 8(d), the signal recovery is correct apart from some very little errors and the autocorrelation signals are not very disturbing. The horizontal white dashed lines indicate the maximal width of the properly reconstructed intralipid layer having a thickness of around 300μm. The 300μm wide scattering layer and a dispersion spread peak having a width of D = 525Pixel in z-space (6.25mm ZnSe) lead to a ratio of , which is below 0.5, therefore predicting an error ratio well below 1%. Again, the experimental results conform with the results of the simulation that there will be only few errors but deviate from a predicted width of a (for the most part) properly reconstructed layer of around 500μm. Namely, due to the increasing number of errors in the center of the capillary, we assume that a homogenous scattering layer with a width of around 300μm (geometrical width) is the limit for a correct signal recovery using 6.25mm ZnSe and our OCT system described above, although the simulation predicts a higher limit. Note that while the dispersive phase and consequently the width of the dispersion spread complex conjugate mirror components corresponding to 6.25mm ZnSe is twice as high compared to 12.70mm SF6, the maximal width of the correct reconstructed scattering layer is not twice as high. This is probably caused by the overlap of autocorrelation and mirror artifacts seen in Fig. 8(b), which is not the case in Fig. 7(b) and was not taken into account in the simulation of Sect. 4.1.
Figure 9 shows the images of an optical longpass filter made from RG850 Schott glass having a thickness of 1.8mm. DC and autocorrelation removal is essential as the two highly reflecting surfaces otherwise cause very strong artifacts near the zero delay position. Comparing Fig. 9(c) and Fig. 9(d) it is obvious how well our advanced DEFR algorithm works and that full range imaging of samples having strong DC and autocorrelation signals requires their removal. After numerical dispersion compensation, the front and back surface of the filter are a little shifted in z-direction, which is probably caused by a remaining linear component in the dispersive phase ϕ(k) due to a non-ideal estimation of the linear component φlin (k) during the phase measurement procedure described in Sect. 3.2 (compare Eq. 33).
Figure 10 shows the images of an eardrum covering a large area in z-direction and crossing the zero delay position. While autocorrelation artifacts are still visible in Fig. 10(c), our advanced DEFR algorithm removes them from the spectrum and reconstructs a very clear image without any artifacts. Even strongly disturbed areas are reconstructed properly. Comparing Fig. 10(c) and Fig. 10(d), it becomes apparent that autocorrelation removal is important in common samples and should not be neglected in full range OCT imaging. The widest scattering layer in the image has a width of around N = 140Pixel in z-space and shows no artifacts, which is accordance with the simulation in Sect. 4.1, where the corresponding ratio of is stating that there will be hardly any errors.
In this paper we presented an advanced processing algorithm for dispersion encoded full range OCT, which uses convolution in z-space instead of multiplikation in k-space, thereby signifi-cantly increasing the processing speed as only one Fourier transform per A-scan is necessary for image reconstruction. Besides the efficient complex conjugate artifact removal, autocorrelation artifacts can be removed likewise with only minor additional computational effort. Furthermore, we demonstrated the usability of ZnSE as a new dispersive material for DEFR OCT with a very high dispersion and less space demand. Apart from the stop criteria, the number of iterations M and the threshold value of noise T, no additional parameters are required for the execution of our DEFR algorithm enabling an automatic stable signal recovery. With increasing dispersion, the spatial frequency span on the detector increases thereby reducing the usable depth range as the Nyquist limit should be fulfilled on the total detector to avoid aliasing artifacts. Therefore, there is some optimal dispersion depending on the image. We presented an estimation for the necessary amount of dispersion depending on the width of the scattering layers in the image. This estimation shows that arbitrary images cannot be reconstructed correctly and thereby solves the philosophical question raised in the introduction that DEFR retrieves N complex values from only N real spectral data points. As the processing speed of our DEFR algorithm only slightly benefits from a larger dispersion (which reduces the amplitude of the complex conjugate artifacts and consequently stops the iteration earlier), our implementation is favorable if the amount of dispersion is limited. We used a standard consumer PC with an Intel Core i5-2400 CPU @ 3.10GHz (4 cores), 8GB RAM and LabVIEW 2010 (64-bit). The processing of 1024 A-scans upsampled to 4096 values in k-space, reconstructed using 6.25mm ZnSe and a maximum number of M = 250 iterations (which are very seldom reached as the mirror artifacts are well below noise before), requires around 1s including autocorrelation removal and 0.9s without, leading to around 1000 A-lines per second (images depicted in Fig. 10). A lower amount of dispersion can even accelerate the processing speed of our algorithm. For example, the inverse Fourier transforms of the dispersive phase functions for 6.35mm SF6 are narrower compared to 6.25mm ZnSe, thus reducing the number of array elements for each subtraction considerably. Using 6.35mm SF6 and M = 250, the reconstruction of 1024 A-scans requires around 0.5s leading to around 2000 A-lines per second (images depicted in Fig. 6). Therefore, using a GPU based processing, which increases the processing speed by a factor of more than 90 , around 100,000 – 200,000 A-lines per second could be reconstructed with our DEFR algorithm with autocorrelation removal. A possible range of applications is high resolution imaging of biological samples that require a large depth scan range, e.g. the human eardrum. Also a future extension to full range Doppler OCT is conceivable as only one measurement is required thereby eliminating the critical phase sensitivity.
The authors would like to thank Peter H Tomlins for providing the nanoparticle target. This project was supported by SAB (Sächsische Aufbaubank, project 14302/2475).
References and links
1. S. Marschall, B. Sander, M. Mogensen, T. Jørgensen, and P. Andersen, “Optical coherence tomography – current technology and applications in clinical and biomedical research,” Anal. Bioanal. Chem. 400, 2699–2720 (2011). [CrossRef] [PubMed]
2. J. Walther, M. Gaertner, P. Cimalla, A. Burkhardt, L. Kirsten, S. Meissner, and E. Koch, “Optical coherence tomography in biomedical research,” Anal. Bioanal. Chem. 400, 2721–2743 (2011). [CrossRef] [PubMed]
3. D. Huang, E. Swanson, C. Lin, J. Schuman, W. Stinson, W. Chang, M. Hee, T. Flotte, K. Gregory, and C. Puliafito et al. “Optical coherence tomography,” Science 254, 1178–1181 (1991). [CrossRef] [PubMed]
4. W. Wieser, B. R. Biedermann, T. Klein, C. M. Eigenwillig, and R. Huber, “Multi-megahertz oct: High quality 3d imaging at 20 million a-scans and 4.5 gvoxels per second,” Opt. Express 18, 14685–14704 (2010). [CrossRef] [PubMed]
7. J. F. de Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearney, and B. E. Bouma, “Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coherence tomography,” Opt. Lett. 28, 2067–2069 (2003). [CrossRef] [PubMed]
8. T. Bajraszewski, M. Wojtkowski, M. Szkulmowski, A. Szkulmowska, R. Huber, and A. Kowalczyk, “Improved spectral optical coherence tomography using optical frequency comb,” Opt. Express 16, 4163–4176 (2008). [CrossRef] [PubMed]
9. M. Hagen-Eggert, P. Koch, and G. Hüttmann, “Analysis of the signal fall-off in spectral domain optical coherence tomography systems,” in Optical Coherence Tomography and Coherence Domain Optical Methods in Biomedicine XVI, Joseph A. Izatt, James G. Fujimoto, and Valery V. Tuchin, eds., Proc. SPIE 8213, 82131K (2012).
10. Y. Mao, S. Sherif, C. Flueraru, and S. Chang, “3x3 mach-zehnder interferometer with unbalanced differential detection for full-range swept-source optical coherence tomography,” Appl. Opt. 47, 2004–2010 (2008). [CrossRef] [PubMed]
11. H.-C. Cheng, J.-F. Huang, and Y.-H. Hsieh, “Numerical analysis of one-shot full-range fd-oct system based on orthogonally polarized light,” Opt. Commun. 282, 3040–3045 (2009). [CrossRef]
12. P. Meemon, K.-S. Lee, and J. P. Rolland, “Doppler imaging with dual-detection full-range frequency domain optical coherence tomography,” Biomed. Opt. Express 1, 537–552 (2010). [CrossRef]
13. T. Wu, Z. Ding, C. Wang, and M. Chen, “Full-range swept source optical coherence tomography based on carrier frequency by transmissive dispersive optical delay line,” J. Biomed. Opt. 16, 126008 (2011). [CrossRef] [PubMed]
14. D. Y. Kim, J. S. Werner, and R. J. Zawadzki, “Comparison of phase-shifting techniques for in vivo full-range, high-speed fourier-domain optical coherence tomography,” J. Biomed. Opt. 15, 056011 (2010). [CrossRef] [PubMed]
16. F. Jaillon, S. Makita, M. Yabusaki, and Y. Yasuno, “Parabolic bm-scan technique for full range doppler spectral domain optical coherence tomography,” Opt. Express 18, 1358–1372 (2010). [CrossRef] [PubMed]
17. 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). [CrossRef] [PubMed]
18. 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). [CrossRef] [PubMed]
19. L. An and R. K. Wang, “Use of a scanner to modulate spatial interferograms for in vivo full-range fourier-domain optical coherence tomography,” Opt. Lett. 32, 3423–3425 (2007). [CrossRef] [PubMed]
20. C.-T. Wu, T.-T. Chi, Y.-W. Kiang, and C. C. Yang, “Computation time-saving mirror image suppression method in fourier-domain optical coherence tomography,” Opt. Express 20, 8270–8283 (2012). [CrossRef] [PubMed]
21. 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). [CrossRef]
22. B. Hermann, B. Hofer, C. Meier, and W. Drexler, “Spectroscopic measurements with dispersion encoded full range frequency domain optical coherence tomography in single- and multilayered non- scattering phantoms,” Opt. Express 17, 24162–24174 (2009). [CrossRef]
23. 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). [CrossRef] [PubMed]
24. B. Hofer, B. Považay, A. Unterhuber, L. Wang, B. Hermann, S. Rey, G. Matz, and W. Drexler, “Fast dispersion encoded full range optical coherence tomography for retinal imaging at 800 nm and 1060 nm,” Opt. Express 18, 4898–4919 (2010). [CrossRef] [PubMed]
25. L. Wang, B. Hofer, Y.-P. Chen, J. A. Guggenheim, W. Drexler, and B. Povazay, “Highly reproducible swept-source, dispersion-encoded full-range biometry and imaging of the mouse eye,” J. Biomed. Opt. 15, 046004 (2010). [CrossRef] [PubMed]
26. L. Wang, B. Hofer, J. A. Guggenheim, and B. Povazay, “Graphics processing unit-based dispersion encoded full-range frequency-domain optical coherence tomography,” J. Biomed. Opt. 17, 077007 (2012). [CrossRef] [PubMed]
27. M. Wojtkowski, V. Srinivasan, T. Ko, J. Fujimoto, A. Kowalczyk, and J. Duker, “Ultrahigh-resolution, high-speed, fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12, 2404–2422 (2004). [CrossRef] [PubMed]
28. B. Cense, N. Nassif, T. Chen, M. Pierce, S.-H. Yun, B. Park, B. Bouma, G. Tearney, and J. de Boer, “Ultrahigh-resolution high-speed retinal imaging using spectral-domain optical coherence tomography,” Opt. Express 12, 2435–2447 (2004). [CrossRef] [PubMed]
29. S. Witte, M. Baclayon, E. J. Peterman, R. F. Toonen, H. D. Mansvelder, and M. L. Groot, “Single-shot two-dimensional full-range optical coherence tomography achieved by dispersion control,” Opt. Express 17, 11335–11349 (2009). [CrossRef] [PubMed]
30. W. Drexler and J. Fujimoto, Optical Coherence Tomography: Technology and Applications (Springer, 2008). [CrossRef]
31. A. Fercher, C. Hitzenberger, M. Sticker, R. Zawadzki, B. Karamata, and T. Lasser, “Numerical dispersion compensation for partial coherence interferometry and optical coherence tomography,” Opt. Express 9, 610–615 (2001). [CrossRef] [PubMed]
32. M. Duarte, M. Davenport, M. Wakin, and R. Baraniuk, “Sparse signal detection from incoherent projections,” in Proc. Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP) 3, III305–308 (2006).
33. B. Karamata, K. Hassler, M. Laubscher, and T. Lasser, “Speckle statistics in optical coherence tomography,” J. Opt. Soc. Am. A 22, 593–596 (2005). [CrossRef]
34. P. Cimalla, J. Walther, M. Mehner, M. Cuevas, and E. Koch, “Simultaneous dual-band optical coherence tomography in the spectral domain for high resolution in vivo imaging,” Opt. Express 17, 19486–19500 (2009). [CrossRef] [PubMed]
35. P. D. Woolliams and P. H. Tomlins, “Estimating the resolution of a commercial optical coherence tomography system with limited spatial sampling,” Meas. Sci. Technol. 22, 065502 (2011). [CrossRef]