1. Introduction

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 [1] or circumferential [10] 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 [11].

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 [14]. 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 [14]. 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 [18]; 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 [14] 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.

2. Theory

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 f_A, as shown in Eq. (1).

 

Fig. 1 Schematic of manual-scanned OCT imaging with time-varying scanning speed.

Download Full Size | PDF

According to Nyquist theorem, the sampling rate, R, has to be larger than twice the highest spatial frequency of the specimen (F_n): R = 1/∆x = f_A/v>2F_n. Therefore, the scanning speed has to be smaller than v_m shown in Eq. (2).

Equation (2) also implies that a scanning velocity smaller than v_m 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 I_x,y(z) is the intensity of an A-scan at (x,y). I_x,y(z) is calculated by taking the square of the amplitude of the A-scan. Denote the complex valued OCT signal as S_x,y (z); then I_x,y(z) = S_x,y(z)S*_x,y(z). Similarly, I_x+Δ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 I_x,y(z) and I_x+Δx,y+Δy(z + Δz).

As we assume the scanning is in x direction, ∆y = ∆z = 0, I_x+Δx,y+Δy(z + Δz) becomes I_x+Δx,y(z) and Eq. (3) becomes:

If we assume the specimen has a homogeneous distribution of scatterers with a uniform scattering strength [19], e.g., the speckle is fully developed, the following relationship exist: <I_x,y(z)> = <I_x+Δx,y(z) > = I₀; <I_x,y(z)²> = <I_x+Δx,y(z)²> = I_RMS². 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 [12] and [20]. 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), |•|² is the square of the amplitude of a complex value. Signal S_x,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 S_x+Δx,y(z) can be expressed as

It is worth mentioning that ∫ indicates integration over (-∞, + ∞) in the expressions of S_x,y(z), S_x+Δx,y(z) and in the following derivations.

Plugging the expression of S_x,y(z) and S_x+Δ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, a₀ 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 [18]. As a result, P(x,y,z) can be written explicitly as P(x,y,z) = P_x(x) P_y(y) P_z(z) and therefore we have:

Lateral PSF P_x(x) can be expressed as:

In Eq. (6), w₀ is the Gaussian beam waist of probing beam [21]. 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].

Plugging Eq. (5) into Eq. (4), we have:

Using the expression of P_x(x) shown as Eq. (6), ρ can be re-written as:

We are able to calculate the integration of Gaussian function over (-∞, + ∞):

Therefore,

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 i_f and i_l 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: x₁, x₂, ..., x_N. Although ∆x_i, the interval between x_i and x_{i + 1}, is not a constant for different A-scans due to non-constant scanning speed, we could extract ∆x_i using ρ_i, the XCC between I_i(z) and I_i+1(z) (A-scans obtained at spatial coordinate x_i and x_{i + 1}). As a result, we were able to estimate ∆x_total, the displacement (in x direction) between the first and the last A-scan in a frame by summing up ∆x_i: ∆x_total = ∑∆x_i. With a pre-set sampling interval ∆x_s, we could calculate the number of A-scans required for this particular frame of data by dividing ∆x_total with ∆x_s: M = ∆x_total/∆x_s. Afterwards, we performed interpolation to obtain A-scan data at spatial points 0, ∆x_s, 2∆x_s,...,M∆x_s to obtain A-scans that were evenly distributed in x dimension.

 

Fig. 2 Flow chart for the scanning speed correction using cross-correlation coefficient.

Download Full Size | PDF

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 [24]. 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 w₀. From our calculation, w₀ 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 w₀ to be 6μm for this system when using Eq. (7) to correlate ρ and ∆x. With the 12μm lateral resolution, F_n = (1/12)(μm⁻¹); 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 [26]. 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 [27].

 

Fig. 3 (a) OCT system diagram; (b) principle of common path interferometer based on single mode fiber; (c) single mode fiber probe integrated with needle.

Download Full Size | PDF

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 [14]. 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 w₀ 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.

 

Fig. 4 Block diagram of signal processing for the OCT system with real-time scanning speed variance correction.

Download Full Size | PDF

To obtain an accurate value of w₀ in our software for our single mode fiber probe, we have first used an estimation of w₀ 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 w₀ that equals (wL)/(M∆x_s). We varied the value of ∆x_s and obtained a consistent value of w₀. All images in section 5.2 were acquired based on w₀ = 18.5μm which was constant for different values of ∆x_s.

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 i^th 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:

 

Fig. 5 (a) Image obtained in calibration experiment (∆x = 0.96μm); (b) relationship between ρ and ∆x obtained by calculating XCC between adjacent A-scans from different B-scans (red circles: experimental; black, dashed line: theoretical); (c) relationship between ρ and ∆x obtained by calculating XCC using A-scans with different offsets from the same B-scan (red, solid line: experimental; black, dashed line: theoretical); (d) the ratio between standard deviation and mean of ρ_i, at different sampling intervals; (e) ratio between σ_Δxtotal and Δx_total.

Download Full Size | PDF

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 (δx_i). The results are shown in Fig. 5(d). With each obtained ρ_i, we used Eq. (8) to calculate a corresponding displacement value Δx_i and then assess the variation of Δx_i, σ_xi². Assume that the displacement between each pair of the A-scans follows the same statistics and the probe travels for a given distance Δx_total. In this case, M, the number of A-scans acquired is Δx_total/δx_i. Based on these assumptions, σ_xtotal², the variance of the estimated displacement approximately equals Mσ_xi². In Fig. 5(e), we show the ratio between σ_Δxtotal and Δx_total, for different values of Δx_total and δx_i.

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 δx_i) 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. Results

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.

 

Fig. 6 (a) Displacement of probing beam versus time (upper); ∆x as a function of time (lower); (b) pseudo B-scan obtained from sinusoidal scanning pattern; (c) upper inset: XCC calculated from adjacent A-scans; lower inset: ∆x calculated from XCC (red) and ground truth ∆x calculated from the driving voltage (black); (d) pseudo B-scan with artifact induced by non-constant scanning speed; (e) B-scan after non-constant scanning speed correction.

Download Full Size | PDF

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 ∆x_s 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.

 

Fig. 7 OCT images obtained from manual scan: (a) before scanning speed correction; (b) after scanning speed correction. Red arrows in Fig. 7(a) indicates areas with motion artifacts.

Download Full Size | PDF

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 ∆x_s 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 ∆x_s, 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 ∆x_s a which is parameter in our software.

 

Fig. 8 Images of IR viewing card with sampling interval ∆x_s equal to 1μm (a), 2μm (b), and 4μm (c).

Download Full Size | PDF

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 (M_i) shown in Fig. 9(c). Afterwards, we extracted zero-crossing points of M_i 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).

 

Fig. 9 (a) photo of quality resolution chart; (b)OCT image obtained from manual scan with scanning speed correction; (c) Blue curve: Mean OCT signal of difference A-scans in Fig. 9(a); red circles: zero-crossing points of the blue curve; (d) OCT image obtained from manual scan without scanning speed correction

Download Full Size | PDF

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).

 

Fig. 10 Manually scanned OCT image of human skin from finger tip (a) and palm (b).

Download Full Size | PDF

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.

 

Fig. 11 Image of onion cells obtained from manual scanning.

Download Full Size | PDF

6. Discussion

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 w₀, 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 w₀. If w₀ 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 [28]. 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 [29], 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.