|Home | About | Journals | Submit | Contact Us | Français|
An interferometric synthetic aperture microscopy (ISAM) system design with real-time 2D cross-sectional processing is described in detail. The system can acquire, process, and display the ISAM reconstructed images at frame rates of 2.25 frames per second for 512 × 1024 pixel images. This system provides quantitatively meaningful structural information from previously indistinguishable scattering intensities and provides proof of feasibility for future real-time ISAM systems.
Interferometric synthetic aperture microscopy  (ISAM) is an imaging modality using computed imaging and synthetic aperture techniques [1–5] to extend the capabilities of instrumentation derived from optical coherence tomography [6–9] (OCT) and optical coherence microscopy [10–12] (OCM). These computed imaging techniques are obtained from physics-based scattering models for OCT and formulated via a solution of the inverse scattering problem. One of the principal demonstrated advantages of ISAM over OCT is the ability to resolve features in the sample that are outside of the confocal region. Ultimately, a more quantitatively accurate and faithful representation of the sample structure is computed.
ISAM requires more computation than OCT and reconstructions have previously been computed offline. Currently, OCT systems provide real-time visualization, an ability that is useful for immediate feedback in time critical situations or for monitoring transient dynamics [13–16]. A real-time ISAM system that could provide spatially-invariant resolution has obvious benefits for the biomedical optics community. One real-time capable system produces an extended focus for OCT scans using a Bessel beam . Other numerical methods for OCT have been used to numerically move the focus to areas of interest via a diffraction kernel [18, 19].
ISAM has previously been described and demonstrated with processing rates of roughly one frame per minute on a modern personal computer [1,2]. In this work, an implementation of a real-time ISAM algorithm is described. The critical functionality of the algorithm is identified where optimizations have been made. In addition, adequate compromises to a prototype algorithm are able to further improve the speed. These improvements are described in detail and the negative effects of the compromises are discussed.
The optimized system can acquire, process, and display the cross-sectional ISAM reconstructed image at frame rates of 2.25 frames per second for 512 × 1024 pixel images on a personal computer with two 3.0 GHz Intel Xeon processors. As previously demonstrated, this system provides quantitatively meaningful structural information from previously indistinguishable scattering intensities and now does so in real-time. The limits to speed now rely on the parallelizability of the processing hardware.
In this paper, the theoretical expression from the signal to the reconstructed object is derived for a prototype algorithm. Then, the prototype algorithm is described in a step-by-step computational procedure. Next, the real-time algorithm is derived where each of the steps is explicitly detailed. Finally, an implementation is described, and experimental results are presented.
An object has a susceptibility represented by η (r||, z), where r|| is the transverse position and z is the longitudinal position. The object may be imaged with a variety of OCT system designs to acquire complex field measurements in time-domain OCT (TD-OCT) S(r||, tω) or spectral-domain OCT (SD-OCT) (r||,ω), which take the arguments transverse position of the beam r||, and time tω or frequency ω. It should be noted here that only real intensities are measured from which the complex representations can be computed unambiguously [20, 21]. To represent the correction of dispersion, a change of variables is introduced from ω to k, the uniformly spaced wavenumber, and from tω to t, time as if waves were propagated in a dispersionless environment. After correcting for dispersion, the TD-OCT signal S (r||, t) and the SD-OCT signal (r||, k) are related by a Fourier transform with respect to k. Similarly, the 3-D Fourier transform of S (r||, t) is S͌(q, k), where q is the transverse wavevector. The analytical coordinates and signals are outlined in Table 1 and Table 2, respectively.
The algorithm for ISAM is derived as in references 2 and 3. Explicitly, a relation is derived that associates the Fourier transform of the object η with the Fourier transform of the signal S. The 2-D Fourier transform of the spectral-domain OCT signal (r||, k), is given by the expression
where η͌ is the 3-D Fourier transform of η, k is the wavenumber, , z0 is the fixed position of the beam focus, A(k) is the square root of the power spectral density, α = π/NA, and NA is the numerical aperture of the beam. The corresponding Tikhonov regularized solution, η͌+, takes the form
β is the longitudinal frequency coordinates of the object, and λ is the regularization constant. The term z0 will be zero when the origin of the z coordinates are at the focal plane. Implementation of this is described in more detail in Step 4 of the prototype algorithm. Thus, rearranging the terms, the pseudo inverse can be rewritten as
The resampling step in Eq. (4) corrects the depth-dependent defocus and is crucial for the performance of the algorithm. The filter B͌(q, k) can be highly singular which introduces noise, hence the need for regularization. Furthermore, applying the filter provides a quantitative, but not a significant qualitative change to the data. In place of Eq. (4), it is sensible to compute the unfiltered solution,
The prototype algorithm was designed in Matlab, where all the calculations were done in double precision, 64-bit. Figure 2 shows a data flow diagram for the prototype ISAM algorithm for 2D imaging where in the paraxial limit the 2D coordinates r|| are separable in each transverse dimension. The discrete coordinates and signals are described in Tables 3 and and4,4, respectively. The digitized spectra Srω is a function of the M discrete positions r in the plane perpendicular to the beam axis and the N discrete frequencies ω. To account for dispersion , a non-uniform interpolation is needed to map the sampled frequencies ω to the sampled wavenumbers k. Similarly, the ISAM solution requires a non-uniform interpolation mapping wavenumbers k to longitudinal frequency coordinates of the object β. The cubic B-spline  is chosen as the non-uniform resampling interpolator; although, a windowed sinc interpolator, such as the prolate-spheroidal interpolator , may be chosen to select specific convergence characteristics. Unfortunately, the frequency response for many non-uniform interpolation procedures drops in performance for frequencies higher than half the Nyquist rate. To mitigate this effect, each spectrum is interpolated by performing a fast Fourier transform (FFT), padding with N zeros, and performing an inverse FFT (IFFT) . The interpolated spectra are stored in a buffer with 2N rows and M columns. Thus, the new interpolated, digitized spectra Srω′ is a function of the sampled positions r and the new sampled frequencies ω′. Similarly, Srk is interpolated to twice the Nyquist rate to get Srk′ as a function of r and the new uniformly sampled wavenumbers k′.
The non-uniformly sampled spectrum in spectral-domain OCT can be corrected in a fashion similar to material dispersion correction  by resampling the Fourier spectrum from ω′ to k. A desired reindexing array in is based primarily on calibrated, second- and third-order correction factors α2 and α3, respectively.
where n is an integer between 0 and N-1, 2N is the number of samples of the interpolated spectrum Srω′, and ωctr is the calculated centroid of the spectrum on a scale from 0 to 1. α2 and α3 are set through a calibration routine. A mirror models a perfect reflector with a minimized width of the point spread function (PSF). Thus, α2 and α3 are adjusted such that the imaged PSF is indeed minimized.
Figure 2 shows a data flow chart of the prototype algorithm which has been broken down into the eight steps listed below.
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 Figs. 1(b) and 1(c), 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. Figure 1(b) 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 Fig. 3.
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.
Coherence imaging measurements are made using a spectral-domain OCT system [27,28], Fig. 4. The system can be minimally reconfigured to also collect time-domain (TD) data. A femtosecond laser (Kapteyn-Murnane Laboratories, Boulder, CO) provides the broadband illumination. The bandwidth of the source is 100 nm, with a center wavelength of 800 nm. The field fluctuates too rapidly to be detected directly, thus the optical signal must be measured with an interferometer. The illuminating source is divided by a 50:50 fiber-optic coupler (splitter) for interference measurements. The 90:10 coupler is used only for time-domain OCT reference measurements, which are not used when imaging in the spectral domain. The associated reference diodes are used to remove the source intensity in TD-OCT. The reference arm contains a delay mirror and the sample arm contains a lens (objective) to focus the Gaussian beam into the sample. In the sample arm, the objective has a focal length of 12 mm, a spot size of 5.6 μm, a confocal parameter (depth-of-focus) of 239 μm, and a NA of 0.05. Since the ISAM algorithm is relatively robust to focal plane alignment, the focal plane can be placed at the center of the image either by a rough measurement or with a calibration tissue phantom.
The spectral detector collects light from the sample and reference arms. In the spectrometer, the light is collimated with a 100 mm focal length achromatic lens and dispersed from a blazed gold diffraction grating (53004BK02-460R, Spectra-Physics), with 830.3 grooves per millimeter and a blaze angle of 19.70 degrees for a blaze wavelength of 828 nm The dispersed optical spectrum is focused using a pair of achromatic lenses each having a focal length of 300 mm. The focused light is incident upon a line-scan camera (P2-2x-02K40, Dalsa.) which contains a 2048-element charge-coupled device (CCD) linear array with 10 × 10 μm detection elements. The spectral resolution provides 2 mm of range depth per axial scan. The camera operates at a readout rate of over 29 kHz, thus one axial scan can be captured during an exposure interval of about 34 μs.
The tissue phantom was placed on a stage and the 2D image data were acquired by scanning the incident beam in the transverse plane using a pair of computer-controlled galvanometer-scanning mirrors. The frame rate depends on the number of axial scans acquired per image or volume. In our experiment, a series of spectrally resolved images (500 × 2048 pixels) was each acquired in 17 ms. A frame capture card (PCI-1428, National Instruments, Inc.) in a computer collects camera data and writes it to the RAM. The frame capture card receives an external trigger from a galvanometer controller card (PCI-6711, National Instruments, Inc.) which activates the frame acquisition. The computer CPUs are dual Intel Xeon processors operating at 3.0 GHz each.
A tissue phantom consisting of TiO2 scatterers suspended in silicone with an average diameter of 1 μm is imaged. Figure 5 shows a transverse scan of 0.4 mm where the real-time cross-sectional OCT image is on the left and the cross-sectional ISAM reconstruction is on the right. The scattering effects of the beam profile in this tissue phantom are evident and shown in Fig. 5 on the left. The previously indistinguishable scatterers are easily identified in the reconstruction, shown in Fig. 5 on the right. The reconstruction better represents the tissue phantom and thus, provides better visualization and a potential for more accurate tissue classification and disease diagnosis. Note that the image boundaries on the left and right, near the top of the phantom, still have some blurring. The solution to the inverse scattering problem is underdetermined for those particular scatterers in these regions because those scatterers were not measured by the full aperture of the beam. The real-time ISAM system provides quantitatively meaningful structural information from previously indistinguishable scattering intensities and provides proof of feasibility for future real-time systems. The acquisition speed has improved from 45 seconds per image to under 450 milliseconds per image, a factor of over one-hundred. Real-time translation of the sample is shown in an online movie (movie1.avi, 6.1 MB).
When imaging 3D volumes, the real-time 2D ISAM processed measurements can be used as a partial solution to the 3D reconstruction. By computationally rotating the data in the transverse direction and re-running the algorithm in the slow scan direction, the full 3D reconstruction can be computed. Figure 6 shows an en face area (0.4 mm × 0.4 mm) of excised human breast tumor tissue imaged with OCT (left) and ISAM processed with the real-time algorithm (right), following post-processing with the same algorithm in the slow scanning direction. The en face slice shown is 450 μm above the focus, such that the cellular and nuclear features in the OCT image are not apparent. Notice that the corresponding features become apparent in the ISAM processed data. A representative histological slice of the tissue is shown in Fig. 7 for comparison. It should be noted that at increasing depths, the probability of multiple scattering increases. Multiple scattering obviates the fundamental assumption of single scattering in ISAM, as it does for OCT. The real-time 2D cross-sectional acquisition movie (movie2.avi, 6.1 MB) and a representative en face movie (movie3.avi, 15 MB) are available online.
Through modifications of a reconstruction algorithm and hardware upgrades, we have demonstrated the feasibility for real-time imaging of the inverse scattering solution. There are a number of modifications made to the algorithm which compromise certain aspects of the processing for speed. The largest drawback appears to be the degradation of the frequency response for the cubic B-spline interpolator, which decreases the SNR for high frequency interpolations. As described in this paper, a variety of interpolators may be used to gauge the tradeoffs in speed versus performance.
The system described achieves spatially invariant resolution for all depths in a cross-sectional scan. Such a system has a potential to provide important cellular scattering information when used to image biological tissues, particularly for real-time, ISAM image-guided surgical interventions. Furthermore, existing OCT systems may be modified to produce similar results with only minor modifications to the processing algorithm. The limits to processing speed now depend largely on the parallelizability of the processing hardware. The parallel nature of this ISAM algorithm suggests that specialized computing hardware such as Graphics Processing Units (GPU) may offer further speed advantages.
This work was supported in part by the National Institutes of Health (Roadmap Initiative/NIBIB, 1 R21 EB005321, and NIBIB, 1 R01 EB005221, S.A.B.), the National Science Foundation (CAREER Award, BES 03-47747, and BES 05-19920, and BES 06-19257, S.A.B., CAREER Award, 0239265, P.S.C.), and the Grainger Foundation (S.A.B.). Human tissue was acquired under Institutional Review Board protocols approved by the University of Illinois at Urbana-Champaign and Carle Foundation Hospital. Additional information can be found at http://www.biophotonics.uiuc.edu.
OCIS codes: (100.3190) Inverse problems; (170.1650) Coherence imaging; (170.4500) Optical coherence tomography; (180.3170) Interference microscopy.