|Home | About | Journals | Submit | Contact Us | Français|
We present a novel algorithm for reconstructing interferograms acquired in optical frequency domain imaging (OFDI). The algorithm was developed specifically for processing in field programmable gate arrays (FPGAs) and featured the use of a finite-impulse-response (FIR) filter implementation of B-spline interpolation for efficiently re-sampling k-space. When implemented in FPGAs, the algorithm allowed for real-time processing of interferograms acquired with a high-speed OFDI system at 54 kHz and a sampling rate of 100 MS/s.
Optical frequency domain imaging (OFDI), also known as swept-source optical coherence tomography (OCT), is a noninvasive imaging modality that has undergone rapid development in recent years –. An interferometric method, OFDI measures depth-resolved reflectance of infrared light with micron-scale spatial resolution. With its high imaging speed, OFDI allows for 3-D image volumes to be acquired across large surface areas in vivo. As such, it has tremendous potential as a screening tool for identifying pathologies in cardiovascular, retinal, and gastrointestinal tissues –.
OFDI images are acquired in the frequency domain, with interferogram samples acquired successively as the peak laser wavelength is rapidly swept in time –. Frequency domain detection gives OFDI systems an inherent sensitivity advantage over their time domain counterparts –. By means of a discrete Fourier transform (DFT), interferograms are transformed from the frequency domain (k-space) to the axial depth domain. Since most wavelength-swept lasers demonstrated to date exhibit nonlinear trajectories in k-space, direct application of the Fourier transform to interferograms acquired with temporally linear sampling results in a loss of axial resolution, especially for frequencies near the Nyquist limit. Although a nonlinearly sampled clock that is synchronized to the laser sweep rate can be used to control the digital acquisition and achieve linear k-space sampling , this approach is potentially problematic at high (>10 MS/s) digitization rates and is incompatible with OFDI systems incorporating frequency shifters. An alternative is to apply interpolation in the data processing algorithm. While this approach imposes additional computational overhead, we demonstrate in this letter that a novel algorithm, optimized for field programmable gate array (FPGA) implementation, can yield real-time processing performance at a sampling speed of 100 MS/s. This algorithm differs significantly from that utilized for real-time processing of signals acquired in the time-domain, which did not involve a DFT or k-space re-sampling . The algorithm’s performance was evaluated with numerical simulations and with FPGAs that were integrated into a high-speed OFDI system.
The “sinc” interpolation method involves convolution of a sampled function with the normalized periodic-sinc kernel, also known as the Dirichlet kernel. A method motivated by the Nyquist-Shannon sampling theorem, periodic-sinc interpolation has been used for reconstruction of OFDI images with high accuracy. In the study of Yun et al. , it consisted of the following steps, in which the FFT was employed to perform the convolution operation.
The parameter M determines the factor by which the interferogram is over-sampled prior to k-space resampling. Qualitatively, it has been observed that little improvement is obtained for M > 8, and that excellent image quality is achieved for M = 2. A disadvantage of the sinc method is that it is computationally expensive: it requires two FFT operations, one real transform of length Nin and one complex transform of length Nin × M, in addition to the linear interpolation in step iv). Real-time processing of images with the sinc method and M ≥ 2 is beyond the capabilities of most commercially-available CPUs and field programmable gate arrays (FPGAs), however.
B-splines provide a continuous representation of discrete samples. As such, they are used in a wide variety of image processing applications such as interpolation, rotation, and edge detection. For this OFDI processing application, we consider a cubic B-spline approximation f3(t) as an empirical model of the continuous-time interferogram,f(t). With this notation, t denotes acquisition time t′ normalized by the time interval separating adjacent samples Δt, so that t = t′/Δt. This approximation can be represented as
where c(n) is an l2 sequence of reals called the B-spline coefficients, and β3(t) is a centered cubic B-spline –. Interpolation of the interferogram in the frequency domain is performed by applying (1) to a set of normalized time points such that the interferogram is uniformly sampled in k-space. This set can be obtained by inverting k(t), the mapping from normalized time to k-space corresponding to the wavelength-swept source. With F3(z), C(z), B3(z) and denoting the z-transforms of f3(n), C(n), and β3(n) respectively, (1) can be written as
For a given interferogram, the coefficients C(z) can be recovered with so-called direct IIR cubic B-spline filter [B3(z)]−1.
Vrcelj et al. suggested performing interpolation with filter coefficients derived from a short finite impulse response (FIR) filter that approximates the IIR cubic B-spline filter . This FIR approach to B-spline interpolation is highly desirable from the standpoint of interpolating OFDI interferograms, because sampled data points can be processed sequentially in real-time as the laser is swept in time. By comparison, sinc interpolation and other methods such as cubic-spline interpolation are non-sequential, necessitating acquisition of the entire interferogram before interpolation is performed. The two computational steps involved with FIR B-spline interpolation were as follows.
The results of  motivated the use of a cubic B-spline filter and a Kaiser window. They also motivated the choice of 1.76 for the Kaiser parameter which controlled the tradeoff between the main-lobe width and the side-lobe area for the Kaiser window.
The FIR B-spline interpolation method was included into a complete OFDI processing algorithm that included background subtraction, demodulation, and postprocessing. Background subtraction was required to remove components of the interferogram that did not derive from the sample. Demodulation was required in order to compensate for the effects of frequency-shifters . The effect of the frequency shifters was to induce a time-varying phase in the interferogram signal. The resulting frequency shift, Δν, was chosen to be exactly one quarter of the sampling rate, vsamp. Conceptually, digital demodulation was performed by multiplying the input signal separately by cos(Δν · n) for the real part of the output and sin(Δν · n) for the imaginary part. In practice, because, Δν = vsamp/4, the interferogram was multiplied by the sequences (1,0,−1,0,…) and(0,1,0,−1,…) to generate the real and imaginary parts, respectively. Digital demodulation mapped a single frequency to two frequencies, the first being lower than the original by vsamp/4, and the second higher by the same amount. The latter frequencies were removed with a symmetric, 40th-order low-pass FIR filter. This filter was designed to remove frequencies with magnitude greater than vsamp/4, preserving the range (−vsamp/4, vsamp/4). In summary, the following steps comprised the complete B-spline reconstruction algorithm optimized for sequential processing, with an input interferogram of length Nin as shown schematically in Fig. 1.
A standard step in OFDI image reconstruction is the multiplication of the interferogram with a window function prior to performing the Fourier transform in order to reduce side-lobe artifacts. In this B-spline reconstruction algorithm, the weights of a Hamming window were included into the interpolation weights as multiplicative factors (step b) of the FIR B-spline interpolation method], so that a separate windowing step was not required. The choice of a Hamming window was based on continuity with previous studies; a different window type could also have been chosen to achieve different axial point-spread function characteristics.
Whether or not the proposed reconstruction algorithm is more computationally efficient than others depends largely on the hardware platform on which it is implemented. It is expected to perform particularly well in contexts in which processing is performed on the interferogram in real-time as the laser is swept in time. Indeed, B-spline interpolation imposes only modest requirements on FPGAs: for each interpolation point, a five point convolution and a weighted sum of at most four terms is computed, relative to a weighted sum of two terms in the case of linear interpolation. The increase in memory requirements required by B-spline interpolation derives only from a 1.5-fold increase in the number of stored weights and, trivially, three unique values corresponding to the symmetric FIR filter. Both B-spline and linear interpolation are therefore well suited to real-time processing. In comparison, the sinc interpolation algorithm involves two FFTs, which can only be performed after an entire interferogram is acquired.
Interferograms from reflectances at different axial positions were simulated numerically, and corresponding point-spread-functions (PSFs) were obtained by application of the B-spline reconstruction algorithm. The quadratic and cubic chirp coefficients that determined the nonlinearity of k-space sampling were measured from the wavelength-swept laser that was incorporated into the OFDI system. Nonlinear k-space sampling was modeled as k(t) = c0 + c1(t − Nin/2) + c2(t − Nin/2)2 + c3(t − Nin/2)3. The ratios c2/c1 and c3/c1, which determine the strength of quadratic and cubic coefficients, were −6.00 and −3.12, respectively. Random white noise was applied additively. The ratio of the interferogram amplitude squared to the variance of the random white noise was 38 dB, which is a signal power typical of that received by tissue. The resulting PSFs were compared with ones that were obtained from three variants of the B-spline reconstruction algorithm. Two variants employed linear and sinc interpolation, respectively; a third employed no interpolation. With these three variants, all other processing steps were identical to those of the B-spline reconstruction algorithm. The importance of PSFs is that, speckle notwithstanding, their profiles completely determine the OFDI image that is obtained from a particular sample . All numerical simulations were performed in Matlab (Mathworks).
For interferogram frequencies with magnitudes approaching the limit of 0.25 (expressed as a fraction of the sampling frequency), the performance of B-spline interpolation was found to be significantly higher than that of linear interpolation and very similar to that of sinc interpolation. Indeed, the PSF corresponding to an interferogram frequency of 0.2 obtained with linear interpolation had a pedestal that was approximately 15 dB above the noise floor [Fig. 2(a)]. This pedestal, which was not apparent with B-spline interpolation or with sinc interpolation, is a direct result of the discrepancy between the interferogram and its piecewise-linear approximation that was employed in linear interpolation. For an interferogram frequency of 0.1, the pedestal was not apparent with all three interpolation methods [Fig. 2(b)]. With lower interferogram frequencies, interpolation inaccuracies tended to be smaller because there was less signal variation across consecutive samples. For both simulated interferograms, PSFs obtained without interpolation were significantly broader than those obtained with interpolation. The signal-to-noise (SNR) ratio was found to be insensitive to the interpolation method: for both interferogram frequencies, it varied by less than 0.2% (in dB) with linear, B-spline, and sinc interpolation methods.
The B-spline reconstruction algorithm was integrated into two FPGAs (Red River, Virtex-2), each of which corresponded to one polarization channel. Each FPGA contained an identical copy of the algorithm, and it received data from a 12-bitanalog-to-digital converter (ADC) that was clocked at MS/s. Operations in the FPGA were performed with 18-bit precision, with the exception of the magnitude-squared operation, which was performed with 36-bit precision and subsequently delivered as output with 32 bits. Each ADC-FPGA pair was contained within a single-board computer (Curtiss-Wright, Champ-AVIII). A custom program written for the single-board computers channeled the outputs of the FPGAs simultaneously to a single RAID hard drive array (Apple Xserve) via separate fiber channels, and to a user-interface computer (Apple, G5; 8GB RAM) via gigabit Ethernet. The RAID array received all data delivered by both single board computers; its write-speeds accommodated continuous data storage. The coefficients (c0,…c3) of the mapping k(t) were measured by a) applying the Hilbert transform to an interferogram corresponding to reflection from a single depth; b) unwrapping the phase of the result from part a); c) fitting a third-order polynomial to the result of step b). Prior to the start of data acquisition, parameters were uploaded to the FPGA from the user-interface computer. These parameters consisted of background signals acquired from the ADCs in the absence of an imaging sample, the B-spline interpolation parameters, and a flag indicating whether to output raw or processed data. The optical elements of the OCT system, including the laser and the interferometer, were developed for high-speed cardiovascular imaging . The average power of the laser was 50 mW; its tuning repetition rate was 54 kHz. For each A-line, the input vector had a length of Nin = 1780 wavelength samples, and the FPGA output consisted of samples Nout = 1024. Two postprocessing steps were performed on the user-interface computer: signal powers from the two FPGAs were added, and the output was logarithmically transformed using a look-up table. The interface computer displayed a subset of the processed frames immediately after FPGA processing; in that sense, real-time display was achieved. The display rate of approximately 5 Hz (with 1080 A-lines per frame) could be significantly improved with software optimization, however.
For imaging of human skin (finger tip) in vivo, a glass window angled at 20° with respect to the incident beam was placed beneath the focusing lens to provide mechanical stabilization and to reduce specular reflections. In a typical image reconstructed in real-time by the two FPGAs, the boundary between the stratum corneum and the papillary dermis and that between glass and the skin surface were sharply delineated [Fig. 3(a)]. The spiral structure of the sweat ducts within the stratum corneum was clearly resolved. Overall, the image quality was good, with no prominent image artifacts apparent. The FPGA-processed image compares favorably with those reconstructed offline in software [Fig. 3(b) and (c)], and therefore it demonstrates the feasibility of accurate real-time processing at sampling frequencies of 100 MS/s with the B-spline reconstruction algorithm.
The advantages and disadvantages of the proposed OFDI reconstruction algorithm compared with existing algorithms depend on the context in which it is implemented. As demonstrated in this study, its accuracy is significantly higher than algorithms with linear interpolation; in general however, higher accuracy can be expected from algorithms with sinc interpolation. The relative performances of different interpolation methods can be expected to depend on the extent to which the trajectory of the laser deviates from linearity in k-space, which in turn depends on the optical configuration of the laser. The reconstruction speeds achieved with the FPGAs utilized in this study were comparable to those that could be achieved with optimized software running on current desktop CPUs. Looking forward, it is likely that FPGA processing of OFDI data will be most valuable for real-time processing in cases where the tuning repetition rate of the source is faster than that of the system utilized in this study , .
A novel OFDI reconstruction algorithm was proposed, and its accuracy was demonstrated with numerical simulations and an FPGA implementation integrated into a high-speed OFDI system. By using an FIR implementation of B-spline interpolation, reconstruction accuracy was shown to be significantly higher than linear interpolation and very similar to sinc interpolation. With interpolation and demodulation performed sequentially in time, the reconstruction algorithm was particularly well suited to FPGA processing. It allowed for continuous, real-time processing in conjunction with sampling frequencies of 100 MS/s. From a clinical standpoint, the significance of this algorithm lies in its ability to make OFDI images available more quickly to clinicians after acquisition, which in turn can lead to faster diagnostics.
The authors are grateful to B. Vrcelj (Caltech DSP Group) and to the anonymous reviewers for their critical feedback.
Adrien E. Desjardins, Harvard Biophysics Program, Boston, MA 02115 USA.
Benjamin J. Vakoc, Wellman Center for Photomedicine, Massachusetts General Hospital, Boston, MA 02114 USA.
Melissa J. Suter, Wellman Center for Photomedicine, Massachusetts General Hospital, Boston, MA 02114 USA.
Seok-Hyun Yun, Wellman Center for Photomedicine, Massachusetts General Hospital, Boston, MA 02114 USA.
Guillermo J. Tearney, Wellman Center for Photomedicine, Massachusetts General Hospital, Boston, MA 02114 USA.
Brett E. Bouma, Wellman Center for Photomedicine, Massachusetts General Hospital, Boston, MA 02114 USA.