Numerical study on spectral domain optical coherence tomography spectral calibration and re-sampling importance

A spectral calibration technique, a data processing method and the importance of calibration and re-sampling methods for the spectral domain optical coherence tomography system were numerically studied, targeted to optical coherence tomography (OCT) signal processing implementation under graphics processing unit (GPU) architecture. Accurately, assigning the wavelength to each pixel of the detector is of paramount importance to obtain high quality images and increase signal to noise ratio (SNR). High quality imaging can be achieved by proper calibration methods, here performed by phase calibration and interpolation. SNR was assessed employing two approaches, single spectrum moving window averaging and consecutive spectra data averaging, to investigate the optimized method and factor for background noise reduction. It was demonstrated that the consecutive spectra averaging had better SNR performance.


Introduction
Optical coherence tomography (OCT) is a noninvasive, contact-free and cost-effective imaging modality providing high resolution cross-sectional images of the probed sample [1,2]. It is a diagnostic imaging technique that employs white light sources with extremely short coherence lengths, making submicron depth resolution attainable. OCT is based on low-coherence interferometry, measuring the interference between the backscattered light from the sample and the light traveling through a reference arm. From a resolution and depth penetration standpoint, OCT complements conventional imaging modalities such that it fills the gap between confocal microscopy and ultrasound imaging. OCT provides depth resolution of typically around 1 μm to 15 μm and penetration depth of around 2 mm-3 mm in biological tissue [3].
In spectral domain optical coherence tomography (SD-OCT), the detected light intensity is a function of the wavelength, I(λ), and the sample's reflectivity depth profile, or A-line, is obtained by applying a Fourier transform over the corresponding modulated optical spectrum. Since there is a nonlinear relationship between the wavelength (λ) and wevenumber (k), and the sampled data is often almost evenly spaced in λ-space (spectrometer detection), it will be then not evenly spaced in wavenumber space (k-space).
Therefore, in order to achieve the highest image quality, spectral calibration is required, and the data has to be re-sampled in k-space in a way that is equally-spaced prior to applying the Fourier transform. Otherwise, applying the fast Fourier transform (FFT) to the unevenly spaced data would yield a depth-dependent increase in axial resolution and consequently blurred images. Numerous approaches have been proposed for re-sampling of the spectral data. These approaches fall into two categories: hardware and software solutions.
As an optics hardware solution, SD-OCT linear-in-wavenumber spectrometers were developed to obtain real time imaging and avert time-consuming calculations for the re-sampling process [4,5]. The detection arm was modified by using a customized prism placed directly after the diffraction grating to evenly disperse the spectrum in the optical frequency or optical wavenumber. These proposed hardware approaches eliminate the need for numerical re-sampling, thus reducing computing time. Another approach uses an external k-trigger electronic board for swept source OCT to acquire the data uniformly sampled in k-space [6]. However, these methods add additional cost and the latter also increases the complexity of the acquisition system.
Software solutions are widely used to make the data uniform in k-space. Numerous approaches have been proposed in order to achieve the best performance, for both OCT signal quality and computational efficiency. For spectral calibration, one approach is to measure the interference spectra at two different positions [7,8]. This eliminates the impact of dispersion mismatch and thus accurately assigns the wavelength to each pixel in the detector. For data processing techniques, the simplest way is to perform a basic linear interpolation method [9,10]. However, this method usually yields poor accuracy results (errors at greater depth). In other approach, the interpolation method was carried out by applying a discrete Fourier transform (DFT), zero padding, inverse DFT and then re-sampling of the resultant data at regular intervals [11,12]. The most popular re-sampling method offering the best cost-performance among the available interpolation methods is the cubic B-spline interpolation [13][14][15]. This method is more accurate; however, is computationally expensive.
This study was focused on the spectral calibration and data processing techniques employed for processing of the SD-OCT spectrometer data targeted to short term implementation in graphics processing unit (GPU) architecture. In this study, white noise was added to the signal, and the signal to noise ratio (SNR) was assessed using two approaches of numerical background noise minimization. The two considered approaches were moving the window over spectrum averaging and consecutive spectra data averaging.

Theory
A simple SD-OCT configuration is represented in Fig. 1. A broad spectrum optical source (BOS) is required. BOS light is split by a broadband beam splitter (BS) into reference and sample paths. Sample path back-scattered light interferes with reference path back-reflected light at the BS. The resultant interference beam is finally dispersed by a diffraction grating (DG) over a detector (e.g. a charge coupled device linear camera known as LCCD). The lens inside the spectrometer focuses a narrow spectral window over each pixel of the LCCD. The total detected interference signal at each pixel of the LCCD detector I(λ) includes the incoherent superposition of reference and sample electromagnetic fields, and the reflection information along depth of the object encoded in the coherent superposition of reference and sample fields, which can be expressed as follows [7]: where a r and a s are the reflection coefficients of the reference mirror and sample, respectively, S(λ) is the spectral shape of the light source, and g(λ) is the dispersion mismatch between the two optical paths. The first two terms can be removed by proper subtraction of the background intensities. The desired information is thus encoded in the last interference term that includes the existence of the dispersion mismatch between the two optical paths: Spectral calibration is performed by a phase fitting method. First, the phase of (2) at the optical path difference (OPD) z 1 is extracted using the Hilbert transform. The same procedure is repeated for the optical path difference z 2 . The difference between two phases (ф) is calculated. The wavenumber (k) is obtained by ф/z, ф=(ф 2 -ф 1 ) and z=(z 2 -z 1 ). The subtraction of two phases taken from two different optical path lengths eliminates the impact of the dispersion mismatch since g(λ) is the same as the optical path length changes.
The vector of N, pixel number, as a function of k is fitted by a polynomial. The attained polynomial coefficients are known as the calibration coefficients, and the corresponding curve may then be used to determine the interpolation points prior to the inverse Fourier transform.
The following section presents the numerical simulations for these procedures for a typical SD-OCT system.

Simulations
All signal processing procedures were implemented under LabVIEW TM 2011 software.

Calibration
The source, S(λ), is simulated by adding two Gaussian functions: where A 1 and A 2 indicate the amplitudes, 1  and 2  are the central wavelengths, b 1 and b 2 are the standard deviations of the Gaussian functions, and λ(n) is the wavelength function going from λ min to λ max in steps of (λ min +λ max ) /N where N is the number of pixels at the LCCD.
Following (2), the photo-detected signal corresponding to an optical path difference z is simulated by modulation of S(λ) by a cosine function as follows: The simulated detected signal is demonstrated in Fig. 2, and the used parameters are shown in Table 1. The constant "const." was set to 0, corresponding to the signal after background removal. The difference of two phase signals acquired for optical path length differences z 1 =500 μm and z 2 =900 μm is computed. k is obtained by ф/z, z 2 -z 1 =400 μm. Figure 3 demonstrates the unwrapped phase obtained at z=900 μm. The vector of N as a function of k, N(k), is fitted by a polynomial order of 3. There is a non-linear relation between the wavenumber and pixel number confirming the necessity of the data interpolation prior to FFT.
After polynomial fitting, the calibration coefficients previously obtained are used to perform the interpolation, therefore obtaining evenly spaced k-space data. These points versus wavenumber k are demonstrated in Fig. 4.

Multiple reflections
Three reflecting surfaces along the sample arm are considered with total optical path differences z 1 , z 2 and z 3 . The resultant signal is implemented from (3) as where z i is the optical path length difference associated to the light reflected at the layer i contributing to the signal, and R i is an array of three elements corresponding to the reflectivity of each layer. In this simulation, the background term is considered as a constant number. Using (5), the photo-detected signal is simulated and demonstrated in Fig. 5.  Table 2.  The calibration coefficients are used to determine the interpolation points of the raw data to a regularly k-spaced corrected signal. Cubic B-spline interpolation is employed. 1024 interpolation points employing calibration coefficients are used for data interpolation. Figure 6 demonstrates non-interpolated and interpolated data spectra. The interpolated data is zero padded to smooth the interpolated signal. The number of zero points is chosen as a power of two since the FFT algorithm requires an input that has a power of two lengths to optimize the Fourier transform calculation. The depth profile is finally achieved by applying inverse FFT to the zero-padded data. The rescaling factor of 2π/k (where k is the difference between k min and k max obtained from the calibration stage) is required to coordinate the peaks representing each reflecting surface at the right position along the x-axis in the reflectivity depth profile. Note that all simulations here do not take into account the fringe wash-out with depth caused by the finite spectrometer resolution. Figure 7 demonstrates the reflectivity depth profile obtained from three reflecting surfaces at OPDs z=200 μm, 300 μm, and 400 μm.

Additive noise
The photo-detected signals of three reflecting surfaces contaminated by uniform white noise of Photonic Sensors 40 unity amplitude are simulated and shown in Fig. 9. Parameter values mentioned in Table 2 are used for this simulation.  Using the same calibration coefficients obtained from the calibration section, 1024 interpolation points are computed. The simulated interference data is then interpolated employing these new points. Figure 10 demonstrates a non-interpolated and an interpolated data spectrum where uniform white noise was added. After interpolation, the inverse FFT is applied to obtain the reflectivity depth profile. Figure 11 demonstrates the reflectivity depth profile attained for three reflecting surfaces at OPDs: z=200 μm,

Importance of calibration and re-sampling methods
Performing the Fourier transform of the collected interference signal requires uniform spacing in wavenumbers to avoid the deterioration of the SNR and also to prevent any degradation in the quality of the images. This fact is shown in  line is the output of the Fourier transform after spectral calibration and interpolation are applied. As it is demonstrated, the solid line is narrower and has higher amplitude proving the importance of the correct wavelength assignment to each pixel of the LCCD in the case of a spectrometer-based configuration.

SNR assessment
Some numerical measurements are performed for the SNR assessment of the simulated SD-OCT system. Uniform white noise with the amplitude of 1 is generated and added to the simulated spectrometer data (Fig. 9). Two approaches were studied to reduce the effect of the background noise and gain better SNR. One approach was to record several consecutive spectra and compute the average. The SNR was measured by 20log 10 (R s /R n ) where R s is the highest amplitude of the signal, and R n is the root mean square of the noise. The results up to 20 consecutive spectra records are shown in Fig. 13.  An ideal number of sets should be considered. It is clear that a big collection of spectra would enhance significantly the SNR, but the speed of SD-OCT would be potentially lost, since averaging over many signals will reduce the temporal response of the system. Figure 13 shows that 15 recordings can be set as an optimized value for this method as the SNR at this number reaches a plateau.
To overcome the potential temporal limitation of spectra averaging another possible approach is to consider a moving window with a specific size and calculate the mean as the window is moved over the whole spectrum. Windows with different sizes, up to 20, are considered. The resultant SNR is shown in Fig. 14 for the same amount of noise. The moving window averaging can be applied to the single spectrum, but might compromise the resolution associated to the point spread function of the SD-OCT system, since the averaging will smooth out sharp signal variations. The two techniques were compared. The results showed that the consecutive spectra averaging method provides a higher SNR than the moving window one. But the difference is not that relevant, so both methods can be potentially applied, preferring spectra averaging if good axial resolution is necessary and possibly using moving window averaging if the acquisition speed