Hand-held OCT systems that offer physicians greater freedom to access imaging sites of interest could be useful for many clinical applications. In this study, by incorporating the theoretical speckle model into the decorrelation function, we have explicitly correlated the cross-correlation coefficient to the lateral displacement between adjacent A-scans. We used this model to develop and study a freehand-scanning OCT system capable of real-time scanning speed correction and distortion-free imaging—for the first time to the best our knowledge. To validate our model and the system, we performed a series of calibration experiments. Experimental results show that our method can extract lateral scanning distance. In addition, using the manually scanned hand-held OCT system, we obtained OCT images from various samples by freehand manual scanning, including images obtained from human in vivo.
© 2012 OSA
Optical coherence tomography (OCT) is a high resolution optical imaging modality widely used in biological and medical fields [1,2]. For many clinical or intraoperative applications, a hand-held OCT system could be particularly useful; it would offer physicians greater freedom to access imaging sites of interest [3–10]. In a hand-held OCT system, it is desirable to have a robust and lightweight probe which can image detailed anatomical structures with a large field-of-view.
In conventional OCT systems, a mechanical scanner steers the OCT probe beam to perform lateral scans. Sequentially acquired A-scans are assembled according to a pre-defined raster  or circumferential  scanning patterns to form two dimensional (2D) or three dimensional (3D) images. Scanners used for OCT include galvanometer-mounted mirrors, piezoelectric transducers (PZT) and microelectromechanical systems (MEMS). Galvanometers have a high linearity and accuracy; however, they are usually bulky and heavy, especially in the case of 3D imaging which requires two galvanometers to perform 2D transverse scans. PZT scanners are smaller than galvanometers and therefore are more suitable for hand-held probes. However, they require a high driving voltage, which is a safety concern. MEMS scanners are smaller but relatively expensive, and they require relatively high voltage .
On the other hand, OCT scans can also be performed manually, similar to manually-scanned ultrasound imaging systems [12,13]. A manually-scanned OCT probe without any mechanical scanner to steer the beam could be much simpler, cost-effective, and easy to use during intraoperative settings . It has been shown that a simple 1D, hand-held OCT probe integrated with standard surgical instruments can be used for 2D OCT imaging and depth ranging during surgery [8,15]. When surgeons manually scan the OCT probe integrated with a surgical tool across the target transversally, the time-varying A-scans can be acquired sequentially and can be used to form pseudo B-scan images. Due to the non-constant scanning velocity of the surgeon’s hand, such a pseudo B-scan results in a non-uniform spatial sampling rate in lateral dimension. Such artifact varies widely between surgeons depending on the stability and dexterity of their hands.
Researchers in the ultrasound community have developed various methods in the last decade to correct the artifact induced by the non-constant scanning velocity in manual scanning, and ultrasound imaging systems have benefited from the use of manually scanned probes. In addition, methods including position tracking and speckle decorrelation have recently been adopted by the OCT community [9,14,16,17]. The speckle decorrelation algorithm is particularly interesting and was demonstrated a few years ago by A. Ahmad et al. in OCT systems for the first time . Compared to a video position tracking system, the speckle decorrelation technique may achieve better accuracy because the dimension of OCT speckle is in the order of micrometers ; which is sufficient for high-resolution OCT with a micrometer-resolution. Speckle decorrelation algorithm is attractive also because it does not require extra hardware components and is easy to implement.
In this study, we incorporated the theoretical speckle model into the decorrelation function to explicitly correlate the cross-correlation coefficient (XCC) to the lateral displacement between adjacent A-scans. We performed a series of experimental calibrations to validate our model and to show that lateral displacement between adjacent A-scans can be extracted quantitatively based on XCC. With the displacement extracted, we were able to correct the artifact induced by the non-constant scanning velocity. To test the method, we built and demonstrated a freehand-scanning OCT system capable of real-time scanning speed correction—for the first time to the best our knowledge. Our system consists of a simple hand-held probe, a spectral-domain OCT engine, and software for image reconstruction and scanning speed correction. We integrated a single-mode fiber with a needle to serve as our probe. The spectral domain OCT engine uses a line scan CCD camera for spectral interferogram acquisition. We also developed our high speed software for real-time signal processing based on a general-purpose graphic processing unit (GPGPU). To demonstrate the system, using a simple 22 gauge needle integrated with a single-mode fiber probe, we obtained OCT images from various samples by freehand manual scanning, including images obtained from human in vivo.
One significant difference between this work and Ref  is that our method can directly calculate lateral displacement from the value of cross-correlation coefficient based on the speckle model we derived. Moreover, in this work, we have developed real-time scanning speed correction algorithm, because the scanning speed correction algorithm proposed in the manuscript could be easily parallelized and implemented using GPGPU.
In this manuscript, we use a Cartesian coordinate system (x, y, z) to describe the 3D space. z indicates the axial direction; x is the lateral direction or the direction of the manual scan. For simplicity, we assume the motion of the OCT needle probe is limited to x-z plane and the specimen is static.
2.1. Manually scanned OCT imaging
As shown in Fig. 1 , when a simple hand-held OCT probe without mechanical scanner is scanned manually in x direction, the displacement between adjacent A-scans, ∆x, is a function of the instantaneous scanning velocity v and the A-scan acquisition rate fA, as shown in Eq. (1).
According to Nyquist theorem, the sampling rate, R, has to be larger than twice the highest spatial frequency of the specimen (Fn): R = 1/∆x = fA/v>2Fn. Therefore, the scanning speed has to be smaller than vm shown in Eq. (2).
Equation (2) also implies that a scanning velocity smaller than vm would lead to oversampling and information redundancy. Under the oversampling condition, there is correlation between adjacent A-scans. The degree of correlation can be measured by Pearson cross-correlation coefficient (XCC) shown as Eq. (3).
In Eq. (3), < > indicates to take the mean value of a signal. Here Ix,y(z) is the intensity of an A-scan at (x,y). Ix,y(z) is calculated by taking the square of the amplitude of the A-scan. Denote the complex valued OCT signal as Sx,y (z); then Ix,y(z) = Sx,y(z)S*x,y(z). Similarly, Ix+Δx,y+Δy(z + Δz) is the intensity of A-scan that is displaced by (Δx,Δy,Δz). σIx,y(z) and σIx+Δx,y+Δy(z+Δz) are the square roots of variance for Ix,y(z) and Ix+Δx,y+Δy(z + Δz).
As we assume the scanning is in x direction, ∆y = ∆z = 0, Ix+Δx,y+Δy(z + Δz) becomes Ix+Δx,y(z) and Eq. (3) becomes:
If we assume the specimen has a homogeneous distribution of scatterers with a uniform scattering strength , e.g., the speckle is fully developed, the following relationship exist: <Ix,y(z)> = <Ix+Δx,y(z) > = I0; <Ix,y(z)2> = <Ix+Δx,y(z)2> = IRMS2. Therefore, we have:
Based on the above relationships, we can simplify Eq. (3) to:
Similar to ultrasound images with fully developed speckle, with the moment theorem for the jointly zero mean, Gaussian random variables, and assuming that the real and imaginary parts of S are uncorrelated, we have [12,13,20],
It is worth mentioning that to derive Eq. (4), we only utilized statistical properties of the random variable involved; such statistical properties are exactly the same for OCT and for ultrasound imaging in  and . On the other hand, such statistical property is not related to the physical mechanisms of OCT image formation; and therefore the Eq. (4) is applicable to OCT signal.
In the Eq. (4), |•|2 is the square of the amplitude of a complex value. Signal Sx,y(z) is determined by the physics of OCT image formation mechanism and can be expressed as the convolution of scattering distribution function a(x,y,z) with system's 3D point spread function (PSF) P(x,y,z)
Similarly, OCT signal Sx+Δx,y(z) can be expressed as
It is worth mentioning that ∫ indicates integration over (-∞, + ∞) in the expressions of Sx,y(z), Sx+Δx,y(z) and in the following derivations.
Plugging the expression of Sx,y(z) and Sx+Δx,y(z) into Eq. (4) and utilizing the fact that OCT system's PSF is not random, we have:
Assuming that the speckle is fully developed and thus scatterers in different spatial location are described by identical but independent random variables, we have the following relationship:
In the above equation, a0 is a constant representing the scattering strength. Using the sifting property of delta function, we have
In OCT, the axial PSF P(z) and the lateral PSF P(x,y) are separable because axial and lateral PSFs are governed by different physical principles: axial PSF is determined by the temporal coherence of the light source while lateral PSF is determined by the imaging optics in the sample arm. Furthermore, in Gaussian optics model, P(x,y) is the product of PSFs in x and y dimensions . As a result, P(x,y,z) can be written explicitly as P(x,y,z) = Px(x) Py(y) Pz(z) and therefore we have:
Lateral PSF Px(x) can be expressed as:
In Eq. (6), w0 is the Gaussian beam waist of probing beam . It is worth mentioning that Gaussian beam waist in this definition is the distance from the beam axis where the intensity of OCT signal drops to 1/e. This PSF expressed in Eq. (6) is valid and consistent with literature [22,23].
Using the expression of Px(x) shown as Eq. (6), ρ can be re-written as:
We are able to calculate the integration of Gaussian function over (-∞, + ∞):
Equation (7) shows that the value of ρ is merely determined by ∆x for fully developed speckle; therefore, we can calculate the cross-correlation coefficient ρ between adjacent A-scans and use the value of ρ to derive the time-varying ∆x as:
It is worth mentioning that ∆x with opposite signs can lead to the same value of ρ, according to Eq. (7). In other words, the scanning direction cannot be determined by calculating XCC.
For digitized sample points in A-scans, ρ can be calculated with Eq. (9).
In Eq. (9), i is the index of pixel in an A-scan and j is the index of A-scan. Segmentation of signal between if and il is selected to calculate ρ.
Although Eq. (3) and Eq. (9) implies that XCC can be either positive or negative, it is unlikely that XCC of two adjacent A-scans has very small or negative value when lateral dimension is highly over-sampled. As demonstrated previously, speckle pattern is formed by convolving random scattering field with OCT system's PSF. Due to the finite dimension of PSF, adjacent A-scans are correlated as long as lateral displacement is small compared to the lateral width of PSF, no matter how sample field is like. However, it is true that XCC can become very small or negative. To deal with this problem, in our software implementation, we took the absolute value of XCC for displacement estimation so that the term inside of logarithm operator is never negative. In addition, we applied thresholding to the value of XCC because small XCC indicates decorrelation between A-scans and thus is not reliable for the displacement assessment. It is also worth mentioning that when XCC turns negative, it does not represent the change of scanning direction.
2.2. Scanning speed correction based on speckle decorrelation
Equation (8) indicates that lateral interval between A-scans can be extracted from the XCC and therefore we could use XCC to correct artifact induced by non-constant scanning speed. The flow chart of the scanning speed correction for one frame of data is shown in Fig. 2 . In our OCT system, the interferometric spectra are detected by a line scan CCD camera and the data is transferred to the host computer through a frame grabber as frames. One frame consists of N spectra acquired at lateral locations: x1, x2, ..., xN. Although ∆xi, the interval between xi and xi + 1, is not a constant for different A-scans due to non-constant scanning speed, we could extract ∆xi using ρi, the XCC between Ii(z) and Ii+1(z) (A-scans obtained at spatial coordinate xi and xi + 1). As a result, we were able to estimate ∆xtotal, the displacement (in x direction) between the first and the last A-scan in a frame by summing up ∆xi: ∆xtotal = ∑∆xi. With a pre-set sampling interval ∆xs, we could calculate the number of A-scans required for this particular frame of data by dividing ∆xtotal with ∆xs: M = ∆xtotal/∆xs. Afterwards, we performed interpolation to obtain A-scan data at spatial points 0, ∆xs, 2∆xs,...,M∆xs to obtain A-scans that were evenly distributed in x dimension.
3. OCT system and software implementation
The basic properties of the OCT system used for our calibration experiments have been reported in our previous work . Briefly, it was a spectral domain OCT system based on GPGPU processing with a 70kHz A-scan rate. In that set-up, an achromatic doublet lens at the sample arm was used to focus the incident light beam and collect back-scattered photons. We used the focal length of imaging objective, focal length of the collimator, and mode field diameter of the fiber to calculate w0. From our calculation, w0 equaled 6μm and this indicated a 12μm 1/e lateral resolution of our OCT. This was further verified by using our OCT system to image 1951 USAF resolution target. The obtained en face OCT image clearly showed that our OCT system can resolve the 5th element in group 6 which corresponds to a 10μm FWHM lateral resolution and therefore 12μm 1/e lateral resolution. Therefore, we assumed w0 to be 6μm for this system when using Eq. (7) to correlate ρ and ∆x. With the 12μm lateral resolution, Fn = (1/12)(μm−1); therefore the largest scanning speed that satisfies the requirement of Nyquist sampling is about 420mm/s, as implied by Eq. (2). For our calibration experiments, a high-speed galvanometer was used to perform lateral scanning.
Our freehand-scanning OCT was a modified version of this system. As shown in Fig. 3(a) , we used a superluminescent diode centered at 840nm with a full width half maximum bandwidth of 55nm (Superlum, S840-B-I-20) as a broadband source. The interferometric signals are dispersed spectrally by a spectrometer that consists of a collimator, diffraction grating (Wasatch Photonics, HD 1800 l/mm @ 840nm), and two identical achromatic doublet lenses (AC508-250-B, f = 250mm) for focusing. The spectra are detected by a line scan CCD camera (e2v, AVIIVA EM4, maximum line rate 70kHz). Spectral data is transferred through a high-speed frame grabber (Matrox Solios eV-CLF) into our host computer (Dell Precision R5500 Rack Workstation, 6-Core Intel® Xeon® Processor X5690 3.46GHz, 12GB RAM, 64bit Windows 7 operating system). For high-speed signal processing, we implemented massive computation on a GPGPU (Nvidia GeForce GTX 480) with 480 cores, with each core operating at 1.4 GHz. To build a simple, lightweight, and small probe which can have arbitrary length, we adopted a common path (CP) configuration for our interferometer [25,26]. In the CP interferometer, the reference and sample light shares the same probe arm which is simply a single mode fiber with its tip cleaved in right angle, as in Fig. 3(b). The reference light comes from the Fresnel reflection at the fiber tip. The sample and reference light is routed by a 50/50 fiber optic coupler to the spectrometer. To protect the fragile fiber tip, we integrated the fiber probe with a 22 gauge needle attached to a syringe so that it can be easily hand-held, shown as Fig. 3(c). The system performance was characterized in our previous work . Despite the fact that Gaussian beam diverges, the lateral resolution of our bare fiber probe is higher than 25μm from results obtained from both experimental measurement and ZEMAX simulation .
The signal processing procedure of our system is briefly summarized in Fig. 4 . A data frame containing N spectra was acquired and transferred to GPU memory. Afterwards, we re-sampled the spectral data from wavelength (λ) space to wavenumber (k) space using cubic spline interpolation and then performed fast Fourier transformation (FFT) to obtain A-scans. With the obtained A-scans, we calculated the XCC between adjacent A-scans using the OCT signal intensity and re-distributed the A-scans using algorithm shown in Fig. 2 to achieve uniform spatial sampling. Before calculating ρ, we processed each A-scan with a moving average filter that averages three adjacent pixels in axial direction and subtracted the output of the filter from each A-scan to reduce low spatial frequency components in the A-scan and therefore increase the sensitivity of cross correlation calculation for lateral motion estimate . Due to the high A-scan rate of our OCT system, lateral dimension was highly over-sampled during the manual scan; therefore simple nearest neighbor interpolation could achieve satisfactory result in re-distributing A-scans and achieve uniform lateral sampling. To match the dynamic range of the OCT data and display device, we applied a truncated log transform to the OCT signal before transferring the signal back to the host computer for display. Procedures shown as red blocks in Fig. 4 were implemented with GPU. Moreover, we implemented multi-thread programming so that data acquisition was in parallel with processing; therefore, we were able to acquire spectral data detected by the CCD continuously as long as we kept the data acquisition rate slightly lower than processing rate. Our software is able to process over 62,000 A-scans every second. Although slightly less than the maximum data acquisition rate of our camera (70k line rate), this speed allows a maximum lateral scanning speed to be approximately 620mm/s assuming w0 is about 20μm for the single mode fiber probe and this lateral scanning speed is significantly larger than moderate manual scanning speed (several millimeters per second). With the software optimization, we will be able to further increase the processing speed in future studies.
To obtain an accurate value of w0 in our software for our single mode fiber probe, we have first used an estimation of w0 based on the experimentally measured lateral resolution of our CP OCT system from our previous work [26,27] and indicate this value as w. Afterwards, we manually scanned a highly scattering phantom for 1cm (L = 1cm) and acquired a certain number of A-scans that were uniformly distributed. We performed such scanning for 10 times and calculated the average A-scan number in all the measurements which is indicated by M. According to Eq. (8), we were able to obtain a better estimation of w0 that equals (wL)/(M∆xs). We varied the value of ∆xs and obtained a consistent value of w0. All images in section 5.2 were acquired based on w0 = 18.5μm which was constant for different values of ∆xs.
4. Calibration of the relationship between cross-correlation coefficient and lateral displacement
In our calibration experiments, we scanned a phantom consisting of 9 layers of cellophane tape to verify the relationship between displacement ∆x and ρ (XCC), shown as Eq. (7) and (8). We used a galvanometer to perform lateral scans with known scanning speeds. We applied a periodical sawtooth voltage V from a function generator to the galvanometer and synchronized V with the acquisition of a frame of data which contained N A-scans (N = 1000). For a 100% duty cycle sawtooth driving voltage, ∆x, lateral interval between adjacent A-scans stays constant because the driving voltage increases linearly during signal acquisition. Therefore, we could calculate the displacement between adjacent A-scans directly from the amplitude of the sawtooth function: ∆x = γV/N. Here γ is a coefficient that relates the driving voltage (V) applied to galvanometer and the probing beam displacement (D) at the focal plane of the imaging lens: γ = D/V. γ was measured to be 1.925mm/V in the OCT setup for our calibration experiments. As a result, by applying different V, we could achieve different scanning speeds and thus different ∆x. We acquired B-scans at various scanning speeds. One example of the image obtained is shown 5(a) which contain 1000 A-scans, with a sampling interval ∆x equal to 0.96μm.
For B-scan acquired at a certain spatial sampling interval determined by the driving voltage, we calculated ρi, XCC between the ith and (i + 1)th A-scan, using pixels within the range as indicated by the double-head red arrow in Fig. 5(a) . Afterwards, with all the ρi obtained, we took the ensemble average of XCC:
Using the above processing, we obtained a value of ρ from each B-scan corresponding to a certain ∆x. The result is shown as red circles in Fig. 5(b). As a comparison, we also plotted the theoretical relationship between ρ and ∆x shown in Eq. (7) as a black, dashed curve in Fig. 5(b). By calculating ensemble average of XCC using A-scans with different offsets from the same B-scan, we obtained decorrelation curve as shown in Fig. 5(c). Similarly, theoretical relationship between ρ and ∆x is shown as black, solid curve in Fig. 5(c). The consistency between experimental results and the analytical model described by Eq. (7) and (8) implies that we can use ρ to quantitatively extract the lateral sampling interval and thus correct non-constant scanning speed.
XCC, as the correlation of two random variables, is inherently a random variable. It is very important to evaluate the statistics of XCC to assess the accuracy in displacement estimation. To evaluate the statistics of XCC, we calculated the standard deviation and the mean of ρi from B-scans acquired experimentally at different spatial sampling intervals (δxi). The results are shown in Fig. 5(d). With each obtained ρi, we used Eq. (8) to calculate a corresponding displacement value Δxi and then assess the variation of Δxi, σxi2. Assume that the displacement between each pair of the A-scans follows the same statistics and the probe travels for a given distance Δxtotal. In this case, M, the number of A-scans acquired is Δxtotal/δxi. Based on these assumptions, σxtotal2, the variance of the estimated displacement approximately equals Mσxi2. In Fig. 5(e), we show the ratio between σΔxtotal and Δxtotal, for different values of Δxtotal and δxi.
As shown in Fig. 5(b) and especially Fig. 5(c), when interval between A-scans is small (<5μm), the measured XCC values are highly consistent with results from theoretical calculation. However, with larger interval (>5μm) between A-scans, the measured XCC values become slightly different from the values calculated using Eq. (7). Moreover, Fig. 5(d) and 5(e) show that a smaller sampling interval (smaller δxi) and larger lateral distance travelled would result in a higher accuracy of displacement estimation. In other words, a large number of sampling points for a given displacement can provide a better displacement estimation, due to the inherent statistics of XCC. Other than the random nature of XCC, there are several reasons for displacement calculated using XCC and actual displacement to be different. First, OCT signal suffers from optical and electrical noise, such as shot noise, excess noise, thermal noise, etc; in addition, OCT signal decorrelates overtime due to random environmental disturbance. Second, the waist of probing beam varies at different lateral displacement due to lens abbreviation. Third, parts of A-scans with low signal intensity produced high correlation due to signal absence, so that the measured XCC was higher than the theoretical values for large displacements in Fig. 5(b) and 5(c). However, those factors might be negligible as compared to the inherent random noise of XCC, because OCT is a high speed and high sensitivity imaging modality.
Results in Fig. 5(d) and 5(e) implies that errors in displacement estimation are minimized with small sampling interval which can vary over time and oversampling in x dimension is necessary for accurate scanning speed correction. Comprehensive evaluation of the error in distance measurement using XCC will be our future work and is beyond of the scope of this manuscript.
5.1 Quantitative lateral sampling interval extraction
We used the same OCT setup for this experiment as in the calibration experiments. To demonstrate that we can use XCC for the quantitative displacement extraction and thus quantitatively correct artifacts from non-uniform scanning speed, we applied a sinusoidal driving voltage (f = 28Hz, Vpp = 1.5V) to the galvanometer and scanned the multilayered tape phantom. Knowing voltage applied to galvanometer and the value of γ from our previous measurement, we were able to calculate the displacement of the probe beam with regards to its neutral position when the voltage applied to the galvanometer was 0. The calculated displacement of the probe beam at different time is shown in the upper inset of Fig. 6(a) . Instantaneous ∆x was calculated by taking the absolute value of the difference between the displacements of adjacent sampling points, as shown in the lower inset of Fig. 6(a). In this experiment, 5000 A-scans were acquired sequentially. By stacking the A-scans, we obtained the pseudo B-scan shown in Fig. 6(b). When the driving voltage/displacement reaches their extreme points, the interval between adjacent A-scans became smallest, which can be clearly seen in Fig. 6(a). With a small sampling interval, data is redundant, as in areas indicated by the arrows in Fig. 6(b). We calculated the XCC between the adjacent A-scans in Fig. 6(b). XCC as a function of time is shown in the upper inset of Fig. 6(c) which was processed with a low pass filter for noise reduction and normalized to the maximum value. We further calculated ∆x using Eq. (8) with XCC obtained in the upper inset of Fig. 6(c) and show the result as the red curve in the lower inset of Fig. 6(c). In the lower inset of Fig. 6(c), we also plot ∆x calculated from the known driving voltage applied to the galvanometer, as the black curve. The consistency between the red and black curve verified our assumption that ∆x could be extracted from XCC quantitatively. There are several reasons for the red and black curve in the lower inset of Fig. 6(c) to be slightly different, as discussed in the previous section. However, the inherence statistics of XCC plays the most significant role in resulting errors. As shown previously in Fig. 5(d), if the sample is laterally homogenous and the displacement between the adjacent A-scans is small, for example less than 1μm, the errors in displacement due to the inherent randomness of XCC are small. Therefore, the errors might be mostly due to other random noise in OCT signal. With a larger sampling interval, errors come from the inherent statistics of XCC rather than other noise in OCT signal. As shown in the lower inset of Fig. 6(c), difference between estimated and actual interval is smaller when interval between A-scans is smaller.
To validate the scanning speed correction algorithm, we took A-scans between 8ms and 13ms when the scanning velocity did not change its direction; thus, we didn’t have to consider the ambiguity of scanning direction. Simply stacking A-scans acquired, we obtained Fig. 6(d) which exhibits an oversampling artifact on the left side of the image, as indicated by the arrow. To remove such artifact, we set sample interval ∆xs to be 5μm and performed nearest neighbor interpolation as described in section 2.2. The resultant image is shown in Fig. 6(e) in which the oversampling artifact is removed.
To demonstrate the effectiveness of our method more clearly, we manually scanned a phantom consisting of multiple layers of tape and saved all the A-scans obtained. Stacking all the A-scans, we obtained Fig. 7(a) which suffers from motion artifacts, as indicated by red arrows. After correcting A-scans using the method proposed in this paper, we obtained Fig. 7(b) which is free of distortion due to non-constant scanning speed. It is worth mentioning that scale bars in Fig. 7 are only applicable to axial dimension.
5.2. Images obtained from manually scanned OCT probe with real-time correction
Using our scanning speed corrected, hand-held OCT system, we manually scanned our single mode fiber probe across the surface of an infrared (IR) viewing card by moving the probe with a freehand. In our real-time scanning speed correction software, we set the spatial sampling interval ∆xs to be 1μm, 2μm, and 4μm and show the corresponding images obtained in Fig. 8(a) , 8(b), and 8(c). The plastic covering film and the underneath fluorescence materials of the IR card can be clearly seen from Fig. 8. With different spatial sampling interval ∆xs, the same physical length is represented by different numbers of A-scans. As the lateral axis of Fig. 8 is A-scan index, the scale of porous structure of the fluorescence materials decreases from Fig. 8(a) to 8(c) due to the increasing sampling interval. Results in Fig. 8 verify that we were able to achieve uniform spatial sampling interval during manual scan and the sampling interval is explicitly determined through ∆xs a which is parameter in our software.
To further evaluate the accuracy of our scanning speed correction method, we imaged a quality resolution chart with 1 line per mm from Edmund Optics through manual scan, as shown in Fig. 9(a) . The red arrow in Fig. 9(a) indicates the scanning direction. Figure 9(b) is the image obtained from the software with real-time correction capability. Periodical structure is clearly shown in Fig. 9(b). To quantitatively evaluate the accuracy of our re-sampling algorithm, we averaged signal amplitude of each A-scan in Fig. 9(b), performed mean subtraction and obtained the blue curve (Mi) shown in Fig. 9(c). Afterwards, we extracted zero-crossing points of Mi to detect the edge of each line, as indicated by red circles in Fig. 9(c). Then we calculated widths of the lines, their mean and standard deviation (STD). The ratio between STD and mean of the width was 0.025, which indicates that our method effectively removed artifacts induced by non-constant manual scanning speed. In comparison, in Fig. 9(d), we show the image obtained from manual scan, but without correction using cross-correlation. Artifacts due to non-uniform scanning speed are clearly visible in Fig. 9(d).
We have also manually scanned the skin of a healthy volunteer using our hand-held OCT probe with 2μm digital sampling interval. To perform manual scan, one of the author held the probe almost perpendicular to the sample surface and moved the probe laterally. Images obtained from the index finger tip and the palm are shown in Fig. 10(a) and 10(b), respectively. White scale bars in Fig. 10 represent 100μm and arrows indicate sweat duct. Epidermis and dermis layer can be visualized in Fig. 10. As light can penetrate deeper in the palm skin, the subcutaneous layer can also be seen in Fig. 10(b).
To further demonstrate the effectiveness of our method on heterogeneous samples, we performed manual scan on an onion sample and show the obtained image as Fig. 11 in which hexagon shaped onion cells can be observed.
Equation (7) forms the mathematical foundation of this work. In deriving Eq. (7), it assumes that the speckle is fully developed. Therefore, to validate Eq. (7) experimentally, we used several different models to test our method. We used a multi-layered phantom without much heterogeneity in lateral dimension for our calibration experiments. However, most OCT images using real specimens have partially developed speckle instead of fully developed speckle. If the sample is heterogeneous, the correlation coefficient between adjacent A-scans not only depends on the lateral interval, but also depends on sample structure. Moreover, when the probe scans across a boundary within the sample, due to the abrupt change in OCT signal, a new A-scan will be attached to the data set regardless of the lateral displacement between the current and previous A-scans. As a result, heterogeneous sample can cause inaccuracy in lateral motion correction of our method. However, for highly scattering samples such as skin, it is usually reasonable to assume that areas corresponding to sample boundary take only a few pixels and therefore do not contribute significantly in the calculation of XCC. As a result, Eq. (7) is valid for highly scattering specimens when a significant portion of the specimen contains homogeneous scatterers, although speckle does not develop fully in most biological specimens. This was verified in the experiment using a quality resolution chart with abrupt changes in lateral features as a sample. We further tested and verified the method using in situ tissue imaging.
As indicated by Eq. (7), to quantitatively estimate ∆x from XCC, it requires that we know w0, the Gaussian beam waist of probing beam which can be calculated or experimentally measured. As a result, the calibrating decorrelation curve shown in Fig. 5(b) is only valid for an OCT system with a certain w0. If w0 used in the calculation for lateral interval between adjacent A-scans is different from the actual beam size, the image reconstructed from our algorithm will be different from the “true” image by a scaling factor in the lateral dimension. However, under such circumstance, uniform sampling can still be achieved to obtain an image that is easy for human to comprehend. In fact, the size of the imaging beam from the probe changes as the beam propagates, and the lateral PSF of OCT system depends on the image depth. Therefore, the speckle decorrelation curve has depth dependency as well. Considering the overall effect, the lateral resolution defined by the Gaussian beam waist is always slightly different from the decorrelation length of OCT signal. Moreover, to reduce the effect of the diverging beam, we took only part of an A-scan to calculate XCC when implementing our software. As a result, the statistics of different pixels within the segment of an A-scan does not vary significantly.
In this work, our assumption is that the manual scan is one dimensional in x axis and that there is no axial motion from the scanning probe, which is not exactly true in a realistic scenario. For example, human hands suffer from physiological tremor and this makes the probe to move randomly in both lateral and axial directions. However, experimental results have shown that the tremor during retinal microsurgery has low frequency motion (<1Hz) with amplitude in the order of 100μm . As a result, with high data acquisition rate, adjacent A-scans usually do not have offset in axial direction for more than a few pixels. Moreover, a cross-correlation maximization-based shift correction algorithm was recently proposed to suppress artifact due to axial motion , which might be helpful to improve the performance of our image acquisition algorithm in the future. In our future study, we are going to conduct a more comprehensive theoretical and experimental study on motion tracking using OCT speckle texture analysis, by considering different types of motion including axial translation, lateral translation and rotation.
Another drawback of our method is that it is not able to determine the direction of scanning; therefore a correct reconstruction of sample structure requires that the scan is in a single direction. This might results in ambiguity if the surgeon's hand moves back and forth around a region of interest. The problem can be solved by combining OCT scan with video microscope technology. The spatial location of OCT tool can be tracked in videos and therefore the scanning direction can be obtained with slightly lower time accuracy.
The research reported in this paper was supported in part by NIH BRP grants 1R01 EB 007969, R21 1R21NS063131-01A1, NIH/NIE R011R01EY021540-01A1, and in part by fellowship support from the ARCS Foundation.
References and links
1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef] [PubMed]
2. A. M. Zysk, F. T. Nguyen, A. L. Oldenburg, D. L. Marks, and S. A. Boppart, “Optical coherence tomography: a review of clinical development from bench to bedside,” J. Biomed. Opt. 12(5), 051403–051421 (2007). [CrossRef] [PubMed]
3. S. A. Boppart, B. E. Bouma, C. Pitris, G. J. Tearney, J. G. Fujimoto, and M. E. Brezinski, “Forward-imaging instruments for optical coherence tomography,” Opt. Lett. 22(21), 1618–1620 (1997). [CrossRef] [PubMed]
5. J. Wu, M. Conry, C. Gu, F. Wang, Z. Yaqoob, and C. Yang, “Paired-angle-rotation scanning optical coherence tomography forward-imaging probe,” Opt. Lett. 31(9), 1265–1267 (2006). [CrossRef] [PubMed]
6. S. Han, M. V. Sarunic, J. Wu, M. Humayun, and C. Yang, “Handheld forward-imaging needle endoscope for ophthalmic optical coherence tomography inspection,” J. Biomed. Opt. 13(2), 020505 (2008). [CrossRef] [PubMed]
7. J. Han, M. Balicki, K. Zhang, X. Liu, J. Handa, R. Taylor, and J. U. Kang, “Common-path Fourier-domain optical coherence tomography with a fiber optic probe integrated Into a surgical needle,” Proceedings of CLEO Conference (2009).
8. M. Balicki, J. Han, I. Iordachita, P. Gehlbach, J. Handa, J. U. Kang, and R. Taylor, “Single fiber optical coherence tomography microsurgical instruments for computer and robot-assisted retinal surgery,” Proceedings of the MICCAI Conference, London, 108–115 (2009).
10. L. Huo, J. Xi, Y. Wu, and X. Li, “Forward-viewing resonant fiber-optic scanning endoscope of appropriate scanning speed for 3D OCT imaging,” Opt. Express 18(14), 14375–14384 (2010). [CrossRef] [PubMed]
11. W. G. Jung, J. Zhang, L. Wang, P. Wilder-Smith, Z. P. Chen, D. T. McCormick, and N. C. Tien, “Three-dimensional optical coherence tomography employing a 2-axis microelectromechanical scanning mirror,” IEEE J. Sel. Top. Quantum Electron. 11(4), 806–810 (2005). [CrossRef]
12. J.-F. Chen, J. B. Fowlkes, P. L. Carson, and J. M. Rubin, “Determination of scan-plane motion using speckle decorrelation: theoretical considerations and initial test,” Int. J. Imaging Syst. Technol. 8(1), 38–44 (1997). [CrossRef]
14. A. Ahmad, S. G. Adie, E. J. Chaney, U. Sharma, and S. A. Boppart, “Cross-correlation-based image acquisition technique for manually-scanned optical coherence tomography,” Opt. Express 17(10), 8125–8136 (2009). [CrossRef] [PubMed]
15. K. Zhang, W. Wang, J. Han, and J. U. Kang, “A surface topology and motion compensation system for microsurgery guidance and intervention based on common-path optical coherence tomography,” IEEE Trans. Biomed. Eng. 56(9), 2318–2321 (2009). [CrossRef] [PubMed]
17. B. Lau, R. A. McLaughlin, A. Curatolo, R. W. Kirk, D. K. Gerstmann, and D. D. Sampson, “Imaging true 3D endoscopic anatomy by incorporating magnetic tracking with optical coherence tomography: proof-of-principle for airways,” Opt. Express 18(26), 27173–27180 (2010). [CrossRef] [PubMed]
18. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). [CrossRef]
19. J. W. Goodman, Statistical Optics (Wiley, 1985).
20. R. F. Wagner, M. F. Insana, and D. G. Brown, “Statistical properties of radio-frequency and envelope-detected signals with applications to medical ultrasound,” J. Opt. Soc. Am. A 4(5), 910–922 (1987). [CrossRef] [PubMed]
21. A. Yariv, Optical Electronics in Modern Communications (Oxford U. Press, 1991).
24. 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]
26. X. Liu and J. U. Kang, “Progress toward inexpensive endoscopic high-resolution common-path OCT,” Proc. SPIE 7559, 755902, 755902-11 (2010). [CrossRef]
27. J. U. Kang, J. Han, X. Liu, K. Zhang, C. Song, and P. Gehlbach, “Endoscopic functional Fourier domain common path optical coherence tomography for microsurgery,” IEEE J. Sel. Top. Quantum Electron. 16(4), 781–792 (2010). [CrossRef]
28. S. Sinch and C. Riviere, “Physiological tremor amplitude during retinal microsurgery,” Proc. 28th IEEE Northeast Bioeng. Conf, 171–172 (2002).
29. J. Lee, V. Srinivasan, H. Radhakrishnan, and D. A. Boas, “Motion correction for phase-resolved dynamic optical coherence tomography imaging of rodent cerebral cortex,” Opt. Express 19(22), 21258–21270 (2011). [CrossRef] [PubMed]