Selection of convolution kernel in non-uniform fast Fourier transform for Fourier domain optical coherence tomography

Kenny K.H. Chan; Shuo Tang

doi:10.1364/OE.19.026891

1. Introduction

Fourier domain optical coherence tomography (FD-OCT) is an imaging modality that provides cross-sectional images with micrometer resolution. Interference fringes measured in the wavenumber k domain is Fourier transformed to axial depth z domain in the reconstruction. The Fourier transform (FT) can be achieved by the computationally efficient fast Fourier transform (FFT) algorithm if the spectral data is uniformly sampled. Unfortunately, this is not the case for most FD-OCT systems. One type of FD-OCT, spectral domain OCT (SD-OCT) is typically designed to capture the spectral components linearly in wavelength λ with a spectrometer. The inverse relationship between wavenumber and wavelength, $k = 2 π / λ$ , causes nonlinearity in the sampling and therefore resampling of the spectral data is needed to achieve uniform spacing in k. The second type, swept source OCT (SS-OCT), time-encodes wavenumber by rapidly tuning a narrowband source through a broad optical bandwidth non-linearly, which also requires a resampling. The image reconstruction algorithm however is interchangeable between the two kinds of FD-OCT.

Various solutions have been developed to address the non-linear k sampling. Hardware approaches such as linear-in-wavenumber spectrometer for SD-OCT [1], linear-k swept laser [2] and k-triggering [3,4] for SS-OCT were successfully implemented but not commonly used due to the added complexity and cost. Usually software techniques based on linear interpolation [5] and cubic spline interpolation [6,7] is used to resample the data. Although these interpolation algorithms offer acceptable processing speed, they generally exhibit large reconstruction error [8–12]. An FD-OCT system employing non-uniform discrete Fourier transforms (NDFT) can directly compute the Fourier transform with matrix multiplication on the unevenly sampled data to produce the exact reconstruction [12]. The slow processing speed of NDFT, however, restricts its use to non-real time applications.

The use of non-uniform fast Fourier transform (NFFT or NUFFT) for FD-OCT processing was recently demonstrated and it was shown to offer better performance than traditional interpolation methods [8–11,13]. NUFFT was formally presented by Dutt and Rokhlin [14] in 1991. NUFFT is an algorithm with reduced computational complexity for approximating the NDFT [14,15]. Its high efficiency has led to its wide acceptance in the field of magnetic resonance imaging (MRI) and is known as gridding for 2D Fourier transform in MRI reconstruction [16–18].

The NUFFT algorithm includes a convolution with a windowing kernel, an FFT on the oversampled data, and a deconvolution. The accuracy of the reconstruction is affected by the choice of kernel function, the parameters related to each kernel, the non-linearity in the sampling and the measured spectral data [17,19]. D. Hillmann et al. reported the use of Kaiser-Bessel (KB) kernel in NUFFT for FD-OCT, but did not study the effect of the kernel design parameters [8]. Several other groups implemented NUFFT for FD-OCT independently using the Gaussian kernel [9,10,13] and KB kernel [11] with optimized design parameters for the respective kernel function. Currently, there is no consensus for the optimum choice of the kernel function for FD-OCT. This manuscript, to our knowledge, will be the first attempt to compare the different kernels in NUFFT for FD-OCT applications. The performance of the NUFFT algorithm will be compared by simulated interferograms and experimental data from a mirror reflector. Finally the NUFFT method is applied to high speed in-vivo imaging of a human finger.

Section 2 describes the NUFFT processing method and the different kernels under consideration. Section 3 provides an error comparison with the simulation of an SD-OCT system. Section 4 presents the experimental accuracy and processing time and demonstrates the proposed method with in-vivo imaging.

2. The NUFFT processing method

2.1 Gridding based NUFFT

In FD-OCT, the axial reflectivity profile of a sample in z domain is obtained by a Fourier transform of the interference fringes sampled in the k domain. However, most systems acquire spectra non-uniformly sampled in wavenumber. In such cases, the axial reflectivity can be directly calculated by using NDFT [12],

f_{N D F T} (z_{m}) = \frac{1}{N} \sum_{n = 0}^{N - 1} F (k_{n}) \exp (- i \frac{2 π}{Δ K} m (k_{n} - k_{o})), m = 0, \dots, N - 1

Here f_NDFT(z_m) is the axial reflectivity at the axial depth location z_m, F(k_n) is the interference signal sampled at non-uniform k spacing, k_o is the wavenumber of the first sample, ΔK is the wavenumber range, and N is the number of sample points. The NDFT computation has a complexity of O(N²). Therefore, it is often desirable to resample the data and use FFT with a complexity O(Nlog₂N) for faster processing speed. This requires that the intensity values to be resampled to points evenly spaced in k. Traditionally linear interpolation [5] and cubic spline interpolation [6,7] are used, but recently more accurate reconstruction with comparable speed by NUFFT with Gaussian [9,10,13] and KB [8,11] kernels have been reported.

NUFFT is a fast algorithm that approximates NDFT [14,15,19,20]. Its steps are summarized in the following:

1) Convolve the non-uniformly sampled k domain data with a kernel function.
2) Resample the result onto evenly spaced locations on an oversampled grid.
3) Compute the Fourier transform using an FFT on the oversampled data.
4) Apply a deconvolution in the z domain by a division of the FT of the kernel function (also known as apodization correction or deapodization).

The speed and accuracy of the NUFFT algorithm can be adjusted by modifying the oversampling ratio and the kernel width. The choice of kernel is also a determining factor in reconstruction accuracy. O’Sullivan noted that the optimal kernel is the infinite length sinc function [21], but kernel functions are typically evaluated only for a finite length for computational practicality.

The first step in the NUFFT algorithm is the convolution of the kernel function to the acquired interference fringes. The convolution of the kernel function C(k) with F(k) gives the intermediate function F_r(k) on the oversampled grid that can be defined in discrete form as,

F_{r} (k_{l}) = \sum_{n = 0}^{N - 1} F (k_{n}) C (k_{l} - k_{n}), l = 0, \dots, N_{r} - 1

where N_r=RN is the number of points on a grid with oversampling ratio R. Note that Eq. (2) completes both convolution and resampling steps of the NUFFT algorithm. Figure 1 illustrates the convolution and resampling process. The kernel width W defines the number of grid points adjacent to the original data point that the kernel is evaluated. The kernel width will be further examined in section 2.3.

Fig. 1 Illustration of the resampling of data with Gaussian kernel into equally spaced grid with a width of 6 [9]. The circles are the original unevenly sampled data and the vertical dashed lines are the new uniform grid. A Gaussian function is convolved with each original data point, spreading its power over a few adjacent grid points as shown by the crosses. The new evenly sampled data value is the summation of the values of all the crosses on each grid line.

Download Full Size | PDF

Once the data points are resampled to a uniform grid, the Fourier transform of Eq. (1) can then be computed using standard FFT algorithm on the oversampled grid with N_r points.

f_{r} (z_{m}) \approx \frac{1}{N_{r}} \sum_{l = 0}^{N_{r} - 1} F_{r} (k_{l}) e^{- i 2 π \frac{m l}{N_{r}}}, m = 0, 1, \dots, N_{r} - 1

Once f_r(z_m) has been calculated, the actual profile f(z_m) can be determined by a deconvolution in k space or alternatively with a simple division by the Fourier transform of C(k) in z space,

f_{N U F F T} (z_{m}) = \frac{f {}_{r}{(z_{m})}}{c (z_{m})}

where c(z_m) is the Fourier transform of the kernel function C(k). The resulting vector f(z_m) will have N_r points, which is larger than the original N points input because of the oversampling. This vector is truncated to N/2 points to remove extra points due to oversampling that contain artifacts [9].

2.2 Sources of error in NUFFT

The main source of error in the NUFFT algorithm is due to aliasing, which is more easily visualized in the Fourier domain. An overview of the NUFFT algorithm is illustrated in Fig. 2 in both domains. Convolution with the kernel function in spectral domain results in a multiplication in the Fourier domain with the kernel response. Resampling with spacing 1/N_r introduces identical copies of the spectrum in the Fourier domain at every N_r interval, which leads to aliasing for broadband signal such as those in FD-OCT [19]. The deconvolution recovers the roll-off caused by the kernel function, but also amplifies the aliasing signal which leads to reconstruction errors. Extra data points due to oversampling are discarded after the deconvolution.

Fig. 2 The overview of the NUFFT algorithm in the spectral domain (k space) and the spatial domain (z space) (color online).

Download Full Size | PDF

A measure of the reconstruction accuracy can be obtained by determining the absolute difference between the NUFFT and NDFT reconstructions as shown in Eq. (5). In the Fourier domain, convolution and deconvolution are multiplication and division respectively. By further recognizing that the aliasing signal is a summation of repeated spectrums, the NUFFT can be rewritten as shown in Eq. (6) where f_NDFT is the exact FT of the interference fringes, N_r is the repetition interval for the repeating spectrum and i is any integer [19]. The structural information of the sample is represented in the data for i=0 and the aliasing signal is located at i≠0, and thus the error equation reduces to Eq. (7) [17,19].

e (z_{m}) = | f_{N U F F T} (z_{m}) - f_{N D F T} (z_{m}) |

= | \frac{1}{c (z_{m})} \sum_{i = - \infty}^{\infty} [f_{N D F T} (z_{m + i N_{r}}) c (z_{m + i N_{r}})] - f_{N D F T} (z_{m}) |

= | \frac{1}{c (z_{m})} \sum_{i \neq 0} f_{N D F T} (z_{m + i N_{r}}) c (z_{m + i N_{r}}) |

As evident in the above equation, the reconstruction error is dependent on c(z_m), the FT of the kernel function. Readers should note that the error is also dependent on the amplitude of f_NDFT(z_m) and therefore it is dependent on the wavenumber sampling pattern and the signal acquired [17,19]. Since the FT of the interference fringes is multiplied with the FT of the kernel function, larger errors will be expected near depth locations with stronger signals and reflectivities, such as interfaces within the sample.

To better represent the error over the full range of the axial scan, the relative L₂ error known also as normalized RMS error is usually used to evaluate the accuracy of the reconstruction [14,19]. The L₂ error gives the root mean square error over the full range of the axial scan. It is defined as [14,19]

e_{2} = \frac{{({\sum_{m} | f_{N U F F T} (z_{m}) - f_{N D F T} (z_{m}) |}^{2})}^{\frac{1}{2}}}{{(\sum_{m} {| f_{N D F T} (z_{m}) |}^{2})}^{\frac{1}{2}}}

In OCT, images are commonly displayed in dB scale. By taking the difference of amplitude (in dB scale) between NUFFT and NDFT reconstructed images on each pixel, the average absolute error of an N point axial scan can also be calculated, which is defined as

e_{a b s} = \frac{1}{N} \sum_{m} | 20 \log_{10} [f_{N U F F T} (z_{m})] - 20 \log_{10} [f_{N D F T} (z_{m})] |

2.3 Kernel width and oversampling rate

The FT of kernel function c(z_m) is dependent on the kernel width W and the oversampling ratio R. These two parameters are fundamental to all NUFFT implementation regardless of the kernel function used. The kernel width W defines the number of grid points adjacent to the original data point that the kernel is accounted for in the NUFFT evaluation. An infinite kernel width would produce the most accurate result, but the value of W is often set to a small finite value for computational efficiency. The use of a finite W value introduces a truncation error [22] because the tails of the kernel are not used. Limiting the kernel function to a width W is equivalent to multiplying a rectangular function to the kernel function in the k domain. The FT of the rectangular function is a sinc function and it contributes to the side lobes as illustrated in Fig. 3 for positive m. The side lobes within the aliasing regions will combine with the region of interest (ROI) at $m / N = 0 to 0.5$ . Increasing the kernel width decreases the amplitudes of the side lobes as illustrated in Fig. 3, which reduces the amount of unwanted signal that is aliased.

Fig. 3 The effect on the kernel response by changing the kernel width W (left) and oversampling ratio R (right) of the Kasier-Bessel kernel function.

Download Full Size | PDF

The oversampling ratio R determines the number of data points N_r in sampling the convolution. In the Fourier domain, increasing R only decreases the amplitudes of the side lobes slightly but allows the main lobe to have a wider width as illustrated in Fig. 3. This decreases the amount of roll-off and increases the signal to noise ratio in the region of interest, where the aliasing signal is considered to be noise in this context. Although the deconvolution or apodization correction can compensate for the roll-off of the original signal, it also amplifies the aliased signal. Therefore it is desirable to minimize the roll-off by increasing the oversampling ratio.

Although the accuracy of the NUFFT algorithm is improved by using a large kernel width and oversampling rate, the processing speed is reduced owing to the extra computation. Designers often choose the smallest values that satisfy the accuracy requirement of the system. The computation complexity of NUFFT will be discussed in section 2.5 following the discussion of the kernel functions.

2.4 Kernel functions

Aside from the kernel width and oversampling ratio, the kernel function C(κ) is central to the NUFFT algorithm. It is used to interpolate and resample the data points onto the uniform oversampled grid. The shape of the function will affect both the roll-off and side lobe amplitudes, and therefore dictates the amount of reconstruction error. Gaussian [9,10,13] and KB [8,11] kernels have recently been reported for use in FD-OCT reconstruction and cosine based kernels have been studied in other fields [16]. Note that there has been no report in the literature on the comparison of the different kernel functions in FD-OCT. The functions under consideration are illustrated in Fig. 4 and are defined as:

Fig. 4 (left) Plot of kernel functions with R = 2, W = 3 and (right) Fourier transform of the kernel function a) Gaussian b) Kaiser-Bessel c) Two term cosine d) Three term cosine.

Download Full Size | PDF

1) Gaussian [9,10,13,15] $C (κ) = e^{\frac{- κ^{2}}{4 τ}}$
2) Kaiser-Bessel [11,17] $C (κ) = \frac{1}{W} I_{o} {β {[1 - {(2 κ / W)}^{2}]}^{\frac{1}{2}}}$
where I_o() is the zeroth-order modified Bessel function of the first kind.
3) Two term cosine (cos2) [16] $C (κ) = α + (1 - α) \cos (\frac{2 π}{W} κ)$
4) Three term cosine (cos3) [16] $C (κ) = α + β \cos (\frac{2 π}{W} κ) + (1 - α - β) \cos (\frac{4 π}{W} κ)$
Here
$κ = | \frac{k_{n} - k_{l}}{δ k} | \leq \frac{W}{2}, δ k = \frac{k_{\max} - k_{\min}}{N_{r}}$

The variables α, β, τ are determined by a function that is dependent on the oversampling ratio and kernel width to optimize accuracy for the approximation to the NDFT. An equation for the optimal value for τ of the Gaussian kernel is [15,19],

τ = \frac{W}{2 N^{2}} \frac{π}{R (R - 0.5)}

and the optimal β value of the KB kernel is given as [11,17]:

β = π {[\frac{W^{2}}{R^{2}} {(R - \frac{1}{2})}^{2} - 0.8]}^{\frac{1}{2}}

No closed form equations have been reported for α, β of the cosine based functions, but Jackson et al. determined the optimum values for a range of kernel widths as summarized in Table 1 [16]. They were chosen to give the lowest average aliasing energy over the image. Jackson et al. did not report values for two term cosine with width ≥ 4 and for three term cosine of width ≥ 6 because the first side lobe alias into the ROI. In these cases, the value for the next closest W was used in this study.

Table 1. The parameters for cosine functions with an oversampling ratio of 2

View Table | View all tables in this article

In addition, to facilitate deconvolution post FFT in the NUFFT algorithm, a solution for the Fourier transform for the Gaussian kernel is given as [15,19]

c (z_{m}) = {(2 τ)}^{\frac{1}{2}} \exp (- z_{m}^{2} τ), m = 0, \dots, \frac{N}{2} - 1

and for the KB kernel it is defined as [11,17]

c (z_{m}) = \frac{\sin ({[{(m π W / N_{r})}^{2} - β^{2}]}^{\frac{1}{2}})}{{[{(m π W / N_{r})}^{2} - β^{2}]}^{\frac{1}{2}}}

Although no closed form equations are given for the cosine kernels, they can be calculated by simply taking the Fourier transform.

From the shapes of the kernel functions displayed in Fig. 4 (left), one might notice that they are very similar but their Fourier transforms in Fig. 4 (right) are significantly different. The ideal kernel function should have a flat main lobe in the ROI and zero amplitudes in the side lobes. By examining Fig. 4, one can see that the kernels have non-ideal responses. The amplitude in the ROI rolls-off and aliasing signals are amplified by the side lobes in the spectral response. Aside from selecting the oversampling ratio and kernel width, it is also important to choose a kernel function that has good response in the Fourier Domain which minimizes the error.

2.5 Computational complexity and memory requirements

In order to choose the optimum parameters and kernel function, their effects on processing time must be also considered. The total complexity of the NUFFT is dependent on kernel width and the oversampling ratio and is given as [19]

O [p W N + R N \log_{2} (R N) + N]

The terms in the above notation represent the convolution-resampling steps, oversampled FFT and the apodization correction respectively. It can be seen that the oversampling ratio R has a larger effect on the complexity than the kernel width. Therefore it is desirable to use the lowest oversampling ratio together with the appropriate kernel width to achieve the desired accuracy [11,17]. FFT algorithms are often more efficient at calculating vectors with a power of two number of elements [23]. To obtain best efficiency, designers should therefore choose a minimum oversampling ratio R>1 that will give a power of two sized oversampled array based on the hardware. For most SD-OCT systems using a camera with power of two number of pixels, the optimal oversampling rate would be two. Other OCT systems with non power of two inputs, such as SS-OCT, should use the smallest fractional oversampling ratio that results in a power of two sized oversampled array. As determined by the hardware used in this study, this manuscript will focus on using an input array of 1024 elements with the optimal oversampling ratio of two. Finally, the kernel calculation constant p in Eq. (19) is dependent on the kernel function and the implementation of the calculation.

For a static sampling pattern that are known a priori, such as in SD-OCT, the values of C(κ) on the grid can be evaluated before the start of the reconstruction to save computation. For this pre-computed mode, the complexity of NUFFT is not dependent on the kernel function and p is constant for all kernel functions. For sampling patterns that are not known a priori or for high speed SS-OCT system with frequency jitters and low repeatability of frequency scanning [24], pre-computation of the kernel function will compromise the reconstruction accuracy. The kernel function should be calculated on-the-fly and p varies for the different kernels. The computation of a Bessel function for the KB kernel is considerably larger than the computation of exponentials and cosines. Therefore the efficiency of on-the-fly mode NUFFT is highly dependent on the choice of kernel functions.

Another consideration is the memory requirement to store the WN number of constants when the kernel function is pre-computed. With increasing interest in high-resolution images (in terms of number of pixels N), high memory demands for these system might make on-the-fly computation preferable, especially those on embedded systems with limited memory such as field-programmable gate array [25] and digital signal processing units [26]. Even with traditional computer systems, large array of kernel constants that cannot fit entirely on the cache of processor are stored in the memory of higher latency, which could decrease system performance. These problems should also be considered when choosing a kernel function for NUFFT.

3. Reconstruction error based on simulated data

To assess the accuracy of using different kernel functions, axial scans at different depths are simulated from measured non-linearity of an SD-OCT system. Interferograms of a prefect reflector are simulated at 17 positions along the imaging depth with 1024 data points. The data set are processed using NUFFT with different kernels while varying the kernel width and oversampling rate. The L₂ error is shown in Fig. 5 and the average absolute error is shown in Fig. 6 .

Fig. 5 Relative L₂ error with varying oversampling ratio and kernel width for different kernel functions. a) Gaussian, b) Kaiser-Bessel, c) Two term cosine, d) Three term cosine.

Download Full Size | PDF

Fig. 6 Average absolute error with varying oversampling ratio and kernel width for different kernel functions. a) Gaussian, b) Kaiser-Bessel, c) Two term cosine, d) Three term cosine.

Download Full Size | PDF

These plots serve as a reference for the amount of error expected. The L₂ error is defined in Eq. (8) as the normalized root mean square error of each A-line. The average absolute error is defined in Eq. (9) and shows the image error in dB scale. Readers should note that dB scale is not sensitive when displaying small errors. Note that increasing the kernel width provides a larger improvement on the accuracy than the oversampling ratio. Compared to other kernels, the KB kernel exhibits the lowest error for W≥3, regardless of the oversampling ratio.

4. Experiment

The NUFFT kernels are evaluated with experimental data from a custom-built SD-OCT system [9]. The source is a superluminescent diode (Superlum) with a center wavelength of 845 nm and a bandwidth of 45 nm. The spectrometer consists of an achromatic collimation lens, a transmission grating (Wasatch Photonics) with 1200 lines/mm and a set of air-spaced lenses with an effective focal length of 100 mm. The linescan CCD camera (Atmel) has 1024 pixels with 14×14 μm² pixel size. The maximum line rate of the camera is 53 kHz and the data transfer is handled by a PCI camera link frame grabber (National Instrument) in base configuration with 32MB of onboard memory. The spectrometer was designed to realize a source limited axial resolution of 7 μm and an imaging depth of 1.73 mm.

4.1 Experimental reconstruction error and computational speed

To measure the reconstruction errors using different kernels, A-lines were acquired from a mirror reflector in the sample arm at 17 positions along the imaging depth. The camera exposure time is 20 μs for each A-line. The A-lines are processed with NUFFT with different kernels at an oversampling ratio of 2 and varying kernel width. Recall that the computational complexity is the same for any kernel function when they are pre-computed. Therefore we only expect a difference in processing time when the kernels are computed on-the-fly. Table 2 summarizes the processing time and A-scan rate of NUFFT and NDFT. The processing time and reconstruction error of NUFFT with different kernels are presented in Fig. 7 . The data is measured on a Dell Vostro 420 with an Intel Q9400 Core 2 Quad (2.66 GHz) and 3 GB of memory. The acquisition and processing program is written in VC + + and is compiled with Intel C++ compiler. The algorithm converts raw data to an image which includes the NUFFT, logarithmic scale calculation, contrast and brightness adjustment as well as on screen display. The Fourier transform of data is calculated using FFTW [23] and the algorithm is accelerated by processing with all four cores available in the computer using OpenMP [27]. Due to the high processing speed of NUFFT, memory transfer latency of the PCI interface can become a bottleneck in the system. This can be improved by using faster data link and batch data transfer. In our experiment, the 32MB onboard memory of the frame grabber is used to temporarily store two frames of data (2×512 A-lines). After the capture of two complete frames the data are transmitted in a batch to the RAM which are then processed and displayed sequentially. This dataflow structure reduces the transfer overhead associated with immediately transferring each A-line to the RAM after capture and allows the system to achieve the maximum camera line rate of 53 kHz,

Table 2. NUFFT computation time (μs) for a 1024 pixels A-line and A-scan rate (lines/s)

View Table | View all tables in this article

Fig. 7 Computation time–error curves of the NUFFT (R = 2) in pre-computed mode (a, c) and in on-the-fly mode (b, d) for kernel width W = 2, 3, 4, 5, 6. The top row shows the L₂ error and the bottom row shows the average absolute error. Increasing the kernel width increases the processing time as expected. The arrow points to the most efficient kernel function: Kaiser-Bessel for the pre-computed mode and the Gaussian for the on-the-fly mode.

Download Full Size | PDF

The NUFFT (R=2, W=3) in the pre-computed mode has a processing time of 10.51 μs for a single A-line. This results in a theoretical reconstruction speed at over 95k A-lines per second. This is two orders of magnitude faster than the NDFT method. The processing speed for on-the-fly mode, however, depends primarily on the kernel functions used. The KB NUFFT (R=2, W=3) has a long computation time of 39.72 μs for one A-line due to the evaluation of the Bessel function. This limits the usefulness of the KB kernel for on-the-fly computation. Conversely, the calculation of the Gaussian function is based on a single exponential evaluation per point, and therefore much less computational expensive than the KB kernel. Although less accurate than the KB kernel at the same width, Gaussian based NUFFT for R=2, W=5 attains an error that is comparable to a KB based NUFFT with R=2, W=3, which makes it a suitable alternative when the kernel function cannot be pre-computed due to system limitations. The Gaussian kernel is often used in multidimensional MRI due to this speed advantage [15,28].

The relationship between processing time and accuracy is shown in Fig. 7. If the processing time and accuracy is weighted equally, the point closest to the origin will represent the most efficient kernel. It can be seen in Fig. 7 that the KB function is the most efficient in the pre-computed mode and the Gaussian function is the best for the on-the-fly mode.

4.2 Reconstruction error in high-speed in-vivo imaging of a human finger

The performance of the reconstruction methods are demonstrated on in-vivo imaging of a human finger. The images are taken with an integration time of 20 μs, which limited the A-scan rate to 50k A-line/s. The raw data is processed first with NUFFT with different kernels at oversampling ratio R=2 and kernel width W=3. Due to the high accuracy of NUFFT and the use of log scale display common to FD-OCT images, there are only subtle observable differences in the images in Fig. 8 .

Fig. 8 Images of the palmar surface of a finger reconstructed using NUFFT with different kernel functions and their differences from NDFT reconstructed image. From top to bottom the rows are Gaussian, Kaiser-Bessel, two term cosine, and three term cosine, respectively. The field of view is 2.5 mm wide by 1.7 mm deep. The scale bars are in dB.

Download Full Size | PDF

The errors can be better visualized by comparing them to an image generated by the NDFT method which is considered as the true image without processing error. The absolute error which is the difference between the NUFFT and NDFT processed images is displayed beside the respective images in Fig. 8. The horizontal lines in the error images represent the peaks of the side-lobes that are aliased, which are located at the same locations for all A-lines. Readers should also note that large errors are found near the regions of strong signals, which is predicted by Eq. (7). This creates artifacts that are especially evident in images reconstructed using the cosine kernels. Those images appear to show higher signal strength and a thicker epidermis layer than suggested by the NDFT reconstructed image. These artifacts could obscure critical structural information.

To further illustrate the change when the kernel width is increased, the same data set are also processed by NUFFT with W = 5 and the results are shown in Fig. 8. Improvements can be seen with all the kernels. Note that reconstruction with Gaussian NUFFT at W = 5 displays an error that is comparable to the KB NUFFT at W = 3, therefore making the Gaussian kernel a suitable choice when on-the-fly computation is needed.

5. Conclusion

Processing FD-OCT data with NUFFT is a recently developed method. To date, NUFFT based on Kaiser-Bessel and Gaussian kernels have been shown independently by several groups to be a highly efficient processing method for FD-OCT. However, the choice of the kernel function, the oversampling ratio and kernel width all affects both the speed and accuracy of the reconstruction, which makes it difficult for designers to choose these parameters. In this paper, the effects of changing these design parameters of NUFFT are investigated systematically on both simulated and experimental data. In addition, the use of NUFFT in the pre-computed kernel mode and on-the-fly mode are examined.

The kernel function can be pre-computed if the sampling non-linearity is static and known a priori. In this mode of NUFFT the processing time is based only on the oversampling ratio and the kernel width, but the accuracy depends on an addition parameter – the kernel function. The most efficient processing method would employ the kernel that produces the lowest reconstruction error for a given R and W. This is found to be the Kaiser-Bessel function in this study.

The Kaiser-Bessel function, however, is difficult to compute. When sampling non-linearity has low repeatability, not known a priori or for memory limited systems, the kernel is computed on-the-fly within the NUFFT algorithm. For these types of systems the Gaussian kernel, requiring only simple exponential evaluations, offers better efficiency than the Kaiser-Bessel kernel.

Increasing the kernel width W other than the oversampling ratio R has been shown to provide a larger improvement on the accuracy. Meanwhile, reducing the oversampling ratio other than the kernel width can reduce the computation complexity more efficiently. Therefore when speed and accuracy are both equally important, designers should choose the lowest oversampling ratio R>1 that will give a power of two sized oversampled array for efficient FFT calculation. The kernel width should then be chosen to achieve the desired reconstruction accuracy.

Although the NUFFT reconstruction errors with the studied kernel functions are only subtly detectable in the images by the human eyes, further reduction in reconstruction error can certainly aid in future application if feature and edge detection are integrated with the display of FD-OCT images.

Acknowledgments

This project is support by the Natural Sciences and Engineering Research Council of Canada and the Canada Foundation for Innovation.

References and links

1. Z. Hu and A. M. Rollins, “Fourier domain optical coherence tomography with a linear-in-wavenumber spectrometer,” Opt. Lett. 32(24), 3525–3527 (2007). [CrossRef] [PubMed]

2. C. M. Eigenwillig, B. R. Biedermann, G. Palte, and R. Huber, “K-space linear Fourier domain mode locked laser and applications for optical coherence tomography,” Opt. Express 16(12), 8916–8937 (2008). [CrossRef] [PubMed]

3. D. C. Adler, Y. Chen, R. Huber, J. Schmitt, J. Connolly, and J. G. Fujimoto, “Three-dimensional endomicroscopy using optical coherence tomography,” Nat. Photonics 1(12), 709–716 (2007). [CrossRef]

4. J. Xi, L. Huo, J. Li, and X. Li, “Generic real-time uniform K-space sampling method for high-speed swept-source optical coherence tomography,” Opt. Express 18(9), 9511–9517 (2010). [CrossRef] [PubMed]

5. G. Liu, J. Zhang, L. Yu, T. Xie, and Z. Chen, “Real-time polarization-sensitive optical coherence tomography data processing with parallel computing,” Appl. Opt. 48(32), 6365–6370 (2009). [CrossRef] [PubMed]

6. M. A. Choma, M. V. Sarunic, C. Yang, and J. A. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11(18), 2183–2189 (2003). [CrossRef] [PubMed]

7. M. Wojtkowski, V. J. Srinivasan, T. H. Ko, J. G. Fujimoto, A. Kowalczyk, and J. S. Duker, “Ultrahigh-resolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12(11), 2404–2422 (2004). [CrossRef] [PubMed]

8. D. Hillmann, G. Huttmann, and P. Koch, “Using nonequispaced fast Fourier transformation to process optical coherence tomography signals,” Proc. SPIE 7372, 73720R (2009). [CrossRef]

9. K. K. H. Chan and S. Tang, “High-speed spectral domain optical coherence tomography using non-uniform fast Fourier transform,” Biomed. Opt. Express 1(5), 1309–1319 (2010). [CrossRef] [PubMed]

10. K. Zhang and J. U. Kang, “Graphics processing unit accelerated non-uniform fast Fourier transform for ultrahigh-speed, real-time Fourier-domain OCT,” Opt. Express 18(22), 23472–23487 (2010). [CrossRef] [PubMed]

11. S. Vergnole, D. Lévesque, and G. Lamouche, “Experimental validation of an optimized signal processing method to handle non-linearity in swept-source optical coherence tomography,” Opt. Express 18(10), 10446–10461 (2010). [CrossRef] [PubMed]

12. K. Wang, Z. Ding, T. Wu, C. Wang, J. Meng, M. Chen, and L. Xu, “Development of a non-uniform discrete Fourier transform based high speed spectral domain optical coherence tomography system,” Opt. Express 17(14), 12121–12131 (2009). [CrossRef] [PubMed]

13. K. Zhang and J. U. Kang, “Real-time intraoperative 4D full-range FD-OCT based on the dual graphics processing units architecture for microsurgery guidance,” Biomed. Opt. Express 2(4), 764–770 (2011). [CrossRef] [PubMed]

14. A. Dutt and V. Rokhlin, “Fast Fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. 14(6), 1368–1393 (1993). [CrossRef]

15. L. Greengard and J. Lee, “Accelerating the nonuniform Fast Fourier Transform,” SIAM Rev. 46(3), 443–454 (2004). [CrossRef]

16. J. I. Jackson, C. H. Meyer, D. G. Nishimura, and A. Macovski, “Selection of a convolution function for Fourier inversion using gridding [computerised tomography application],” IEEE Trans. Med. Imaging 10(3), 473–478 (1991). [CrossRef] [PubMed]

17. P. J. Beatty, D. G. Nishimura, and J. M. Pauly, “Rapid gridding reconstruction with a minimal oversampling ratio,” IEEE Trans. Med. Imaging 24(6), 799–808 (2005). [CrossRef] [PubMed]

18. G. E. Sarty, R. Bennett, and R. W. Cox, “Direct reconstruction of non-Cartesian k-space data using a nonuniform fast Fourier transform,” Magn. Reson. Med. 45(5), 908–915 (2001). [CrossRef] [PubMed]

19. A. J. W. Duijndam and M. A. Schonewille, “Nonuniform fast Fourier transform,” Geophysics 64(2), 539–551 (1999). [CrossRef]

20. J. A. Fessler and B. P. Sutton, “Nonuniform fast Fourier transforms using min-max interpolation,” IEEE Trans. Signal Process. 51(2), 560–574 (2003). [CrossRef]

21. J. D. O’Sullivan, “A fast sinc function gridding algorithm for fourier inversion in computer tomography,” IEEE Trans. Med. Imaging 4(4), 200–207 (1985). [CrossRef] [PubMed]

22. D. Potts, G. Steidl, and M. Tasche, “Fast Fourier transforms for nonequispaced data: a tutorial,” in Modern Sampling Theory: Mathematics and Applications, J. J. Benedetto and P. Ferreira, eds. (Springer, 2001), Chap. 12, 249–274.

23. M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proc. IEEE 93(2), 216–231 (2005). [CrossRef]

24. B. Liu, E. Azimi, and M. E. Brezinski, “True logarithmic amplification of frequency clock in SS-OCT for calibration,” Biomed. Opt. Express 2(6), 1769–1777 (2011). [CrossRef] [PubMed]

25. T. E. Ustun, N. V. Iftimia, R. D. Ferguson, and D. X. Hammer, “Real-time processing for Fourier domain optical coherence tomography using a field programmable gate array,” Rev. Sci. Instrum. 79(11), 114301 (2008). [CrossRef] [PubMed]

26. A. W. Schaefer, J. J. Reynolds, D. L. Marks, and S. A. Boppart, “Real-time digital signal processing-based optical coherence tomography and Doppler optical coherence tomography,” IEEE Trans. Biomed. Eng. 51(1), 186–190 (2004). [CrossRef] [PubMed]

27. OpenMP Architecture Review Board, “The OpenMP API specification for parallel programming,” http://www.openmp.org/.

28. T. S. Sorensen, T. Schaeffter, K. O. Noe, and M. S. Hansen, “Accelerating the nonequispaced fast Fourier transform on commodity graphics hardware,” IEEE Trans. Med. Imaging 27(4), 538–547 (2008). [CrossRef] [PubMed]

Kernel Width	Two term cosine	Three term cosine
Kernel Width	α	α	β
2	0.5273	0.4149	0.4990
3	0.5076	0.3954	0.4977
4		0.3850	0.5000
5		0.3823	0.5000

Kernel function	Pre-computed mode (R=2, W=3)		On-the-fly mode (R=2, W=3)		On-the-fly mode (R=2, W=5)
Kernel function	Processing time	A-scan rate	Processing time	A-scan rate	Processing time	A-scan rate
Gaussian	10.51	95.14k	18.01	55.52k	23.77	40.27k
Kaiser-Bessel			39.72	25.17k	59.42	16.83k
Two term cos			20.01	49.96k	24.95	40.08k
Three term cos			25.56	39.13k	36.15	27.66k
NDFT	1470	0.6803k	-	-	-	-

Kernel Width	Two term cosine	Three term cosine
Kernel Width	α	α	β
2	0.5273	0.4149	0.4990
3	0.5076	0.3954	0.4977
4		0.3850	0.5000
5		0.3823	0.5000

Kernel function	Pre-computed mode (R=2, W=3)		On-the-fly mode (R=2, W=3)		On-the-fly mode (R=2, W=5)
Kernel function	Processing time	A-scan rate	Processing time	A-scan rate	Processing time	A-scan rate
Gaussian	10.51	95.14k	18.01	55.52k	23.77	40.27k
Kaiser-Bessel			39.72	25.17k	59.42	16.83k
Two term cos			20.01	49.96k	24.95	40.08k
Three term cos			25.56	39.13k	36.15	27.66k
NDFT	1470	0.6803k	-	-	-	-

Selection of convolution kernel in non-uniform fast Fourier transform for Fourier domain optical coherence tomography

Abstract

1. Introduction

2. The NUFFT processing method

2.1 Gridding based NUFFT

2.2 Sources of error in NUFFT

2.3 Kernel width and oversampling rate

2.4 Kernel functions

2.5 Computational complexity and memory requirements

3. Reconstruction error based on simulated data

4. Experiment

4.1 Experimental reconstruction error and computational speed

4.2 Reconstruction error in high-speed in-vivo imaging of a human finger

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (8)

Tables (2)

Equations (19)

Optics Express