There are a number of operations hindering the performance of this prototype ISAM algorithm which have been overcome. Using C++ instead of Matlab allows for a number of speed improvements. The 64-bit operations such as the FFT and interpolations can be replaced with 32-bit operations with a small, but tolerable, increase of quantization noise. A speed advantage is gained by replacing the ‘for’ loops within Matlab scripts by vectorized Intel SSE (Streaming SIMD Extentions) commands and/or compiled ‘for’ loops. Time-expensive, built-in Matlab interpolation and FFT functions are replaced with hardware optimized functions as found in the Intel Math Kernel Library (MKL). An FFT of the real spectrum is implemented using the customized real-to-complex FFT in the MKL.
The focal plane is fixed such that t0=0 by imposing a restriction that the focus be placed at the center of the image. Therefore, the complex analytic signal does not need to be padded with any zeros, and thus the computation time of the FFT is reduced because the FFT is always a power of two. While the prototype functions utilized dynamically-allocated memory, the real-time program allocates memory prior to imaging time. Similarly, a table of resampling coefficients for the dispersion compensation and the ISAM reconstruction are pre-calculated and stored in memory. The prototype algorithm is interpolated to twice the Nyquist rate to mitigate the high frequency roll-off of the non-uniform interpolator. Although the interpolation has good frequency response, it is not necessary, especially when speed is an issue. The phase stabilization procedures, which might be needed for long acquisition periods, are not necessary for 2D imaging since this SD-OCT system provides good phase stability over the short duration of the scan, about 17 ms. A good estimate of maximum drift for maintaining phase stability is no more than a half of a wavelength for the length of the synthetic aperture, the depth dependent beam width.
The key computation for ISAM is the resampling in the Fourier space. The new design is an efficient routine which requires two interpolations of the columns, one one-dimensional (1D) real-to-complex (R2C) fast Fourier transform (FFT) of the columns, and two 2D FFTs. The FFT routine from the Intel Math Kernel Library (MKL) was used for all 1D and 2D transforms. The 1D 2048-point R2C FFT is calculated for every column of the 512 × 2048 real-valued data, while the 2D FFT and 2D IFFT are calculated for the 512 × 1024 complex values.
The reindexing array in and the corresponding interpolation table is pre-computed and stored in random access memory (RAM) for repeated use. The coefficients of the cubic B-spline are pre-computed. The same type of calculations could be simply done for a table of linear spline coefficients, but the cubic B-spline is described below for completeness. The integer parts of the indices used for calculation of the coefficients to the cubic B-spline are given by
The fractional coefficients are given by
Similarly, the Fourier resampling indices of ISAM need to be determined and pre-calculated. To determine the ISAM resampling indices, the constants which specify the axial and transverse bandwidths of the object are determined. These are based on system parameters and are βmin, βmax, qmin, and qmax, respectively. By approximating the longitudinal bandwidth parameters of the object, one can avoid computationally costly zero-padding for calculation of the resampled solution. In this case, βmin and βmax are well approximated with the range of the optical bandwidth, qmin is set to zero, and qmax is set to the maximum transverse scanning frequency. There may be a small loss of spectral information across the βmin and βmax boundary as illustrated by , but this has only a marginal effect on the signal-to-noise ratio. The maximum and minimum wavenumbers are calculated for the region of interest,
A range of values for q[m] and β[n] is created in the 2-D Fourier space and the resampling grid kq[m, n] is pre-calculated. Notice here that M and N′ = N/2 are assigned according to number of rows and columns in the complex analytic sampled Fourier data.
The kq[m, n] grid is used to calculate a table of values for the cubic B-spline coefficients. The values are calculated as shown above except each column of resampling parameters is different and therefore must also be saved. shows the plot of the kq[m, n] grid where the curved lines represent the constant values of β[n]. This 2D grid is used to calculate another table of cubic B-spline coefficients.
Similar to the spline coefficient calculations shown above, the reindexing array kq[m, n] and the corresponding cubic B-spline interpolation coefficient table is pre-computed and stored in random access memory (RAM) for repeated use. The integer part of the index used for calculation of the in the cubic B-spline is given by
The fractional coefficients are given by
Using the pre-calculated tables, a flow diagram of the real-time algorithm is shown in .
Computational flow chart for memory allocation for each step of the inverse scattering algorithm.
Here Srω[m, n] is the raw interferometric data captured from the camera and has M columns and N rows. In this implementation, M=512 columns and N=2048 rows.
- Step 1. The pre-calculated table is used to perform the interpolation as follows.
for all integers 0 ≤ n < N and 0 ≤ m < M.
- Step 2. The real-to-complex 1-D FFT routine from the Intel Math Kernel Library (MKL) is used on all the columns.
The real-to-complex FFT will compute N/2 complex values. The new number of rows of the complex data is given by N′ = N/2.
- Step 3. The contribution of the noise from the average spectral intensity on the detector is removed by setting Srt[m, n] equal to zero at the t=0 plane. Also, Srt[m, n] is circularly shifted by half such that the focus will be the new t′=0 plane.
- Step 4. The complex 2-D inverse FFT (IFFT) of the complex analytic signal Srt′[m, n] is calculated
- Step 5. The pre-calculated table is used to perform the cubic B-spline interpolation as follows.
where the calculated cubic B-spline coefficients are from the lookup table.
- Step 6. The complex 2-D FFT of the Fourier transformed object ηqβ [m, n] is calculated
- Step 7. ηrz′[m, n] is circularly shifted such that the focus is in the middle of the image,
then the magnitude |ηrz [m, n]| is displayed.