Robust motion tracking based on adaptive speckle decorrelation analysis of OCT signal

: Speckle decorrelation analysis of optical coherence tomography (OCT) signal has been used in motion tracking. In our previous study, we demonstrated that cross-correlation coefficient (XCC) between Ascans had an explicit functional dependency on the magnitude of lateral displacement ( δx ). In this study, we evaluated the sensitivity of speckle motion tracking using the derivative of function XCC( δx ) on variable δx . We demonstrated the magnitude of the derivative can be maximized. In other words, the sensitivity of OCT speckle tracking can be optimized by using signals with appropriate amount of decorrelation for XCC calculation. Based on this finding, we developed an adaptive speckle decorrelation analysis strategy to achieve motion tracking with optimized sensitivity. Briefly, we used subsequently acquired Ascans and Ascans obtained with larger time intervals to obtain multiple values of XCC and chose the XCC value that maximized motion tracking sensitivity for displacement calculation. Instantaneous motion speed can be calculated by dividing the obtained displacement with time interval between Ascans involved in XCC calculation. We implemented the above-described algorithm in real-time using graphic processing unit (GPU) and demonstrated its effectiveness in reconstructing distortion-free OCT images using data obtained from a manually scanned OCT probe. The adaptive speckle tracking method was validated in manually scanned OCT imaging, on phantom as well as in vivo skin tissue.

In our previous study, we have demonstrated that Pearson cross-correlation coefficient (XCC) between slightly displaced OCT Ascans had an explicit functional dependency on the magnitude of lateral motion [25]. Based on this finding, we developed a method for motion tracking by calculating XCC between subsequently acquired Ascans and converting XCC value to displacement. It was noted that there existed an optimized degree of decorrelation motion tracking. When the decorrelation between adjacent Ascans was very small (XCC close to 1) or very big (XCC close to 0), the result of motion tracking was dominated by noise. As a result, it was challenging to track motion speed that could vary in a large range using XCC between Ascans acquired with the same time interval. Therefore, there is an unmet need to achieve improved robustness in motion tracking based on OCT speckle analysis.
In this study, we evaluated the sensitivity (S) of speckle motion tracking using the derivative of function XCC (δx) on variable δx, where δx indicates lateral displacement. We demonstrated that the magnitude of the derivative could be maximized and the optimized motion tracking robustness could be achieved using signals with appropriate amount of decorrelation for XCC calculation. Based on this finding, we developed an adaptive speckle decorrelation analysis strategy for robust motion tracking. For adaptive speckle tracking, we calculated a set of XCC values using Ascan pairs obtained with different intervals and selected the XCC that maximized S for motion tracking. Instantaneous motion speed can be calculated by dividing the obtained displacement with time interval between Ascans involved in XCC calculation. We also implemented the algorithm for adaptive speckle decorrelation analysis in real-time using graphic processing unit (GPU) and demonstrated its effectiveness in reconstructing distortion-free OCT image using data obtained from a manually scanned OCT probe. Imaging experiments were conducted on phantom as well as in vivo skin tissue.

Principle
In this study, we use a Cartesian coordinate system (x, y, z) to describe the 3D space. Light beam propagates along z direction. x and y are lateral coordinates orthogonal to z axis. In general, motion in 3D space is characterized by a vector with three independent components, δx, δy and δz. We assume the sample is static and the motion is attributed to light beam scanning in x direction. Therefore, the only nonzero component in the motion is δx, and δy = δz = 0.
In our previous work, we theoretically derived and experimentally validated that the Pearson correlation coefficient XCC (ρ) between displaced Ascans has an explicit functional dependency on spatial displacement, as shown in Eq. (1) where ω 0 indicates the beam waist of the probing beam [25]. Equation (1) suggests δx can be obtained from speckle decorrelation analysis.
An estimation of XCC ( ) between Ascan I t (z) and I t+Δt (z) is calculated using Eq. (2).
In Eq. (2), < > indicates taking the mean value of a signal; I t (z) and I t+Δt (z) are the intensity of Ascans acquired at time t and t + Δt; m t = <I t (z)> and m t+Δt = <I t+Δt (z)> are the mean values for I t (z) and I t+Δt (z), respectively; σ t and σ t+Δt are the standard deviations for I t (z) and I t+Δt (z), respectively.
With XCC calculated using Eq. (2), an estimated value of lateral displacement ( ) can be obtained as below.

#247422
Due to the nonlinear relationship between ρ(δx) and δx (Eq. (1)) and between  and (Eq. (3)), the derivative of ρ(δx) on δx is not a constant. As a result, the same infinitesimal change in lateral displacement (dδx) will result in different change in XCC (dρ), for different value of δx. In practice, the estimation of ρ suffers from inaccuracy, and the same amount of inaccuracy in estimating ρ will lead to different amount of inaccuracy in estimating displacement δx. Therefore, it is critical to select an optimized operating point for decorrelation analysis. This is illustrated in Fig. 1 where the functional relationship between ρ and δx is plotted according to Eq. (1) with ω 0 = 17μm. Consider three values of XCC (ρ 1 , ρ 2 and ρ 3 in Fig. 1) and assume the same amount of uncertainty ( ± 0.035) in XCC measurement (red bars in vertical direction). Clearly, the same uncertainty in XCC will lead to different amount of inaccuracy in displacement estimation (blue bars in horizontal direction with different length) for different operating point in speckle decorrelation analysis, due to the nonlinear relationship between ρ and δx. We define the sensitivity (S) of speckle motion tracking as the absolute value of the derivative of ρ(δx) on δx (Eq. (4)) [26]. Figure 2 show how S varies as δx. Clearly, there exists a maximum value of S (S max ). Analytically, the value of δx that maximizes S can be found when dS/dδx = 0. This relationship is shown explicitly in Eq. (5). Clearly, S is maximized when δx max = ω 0 /√2 or ρ max = e -1/2 .
In this study, we propose to use the XCC value that approximately equals ρ max and maximizes S for robust motion tracking. In other words, we select XCC value most sensitive to change in displacement δx for speckle decorrelation analysis [26]. This is because a large S is desirable for speckle tracking. With a large S, change in displacement δx can result in measurable change in XCC. Otherwise, the measurement of XCC is overwhelmed by random noise [27]. For example, if two Ascans with small displacement (δx<< ω 0 ) or large displacement (δx>>ω 0 ) are chosen for decorrelation analysis, S is approximately 0, leading to a compromised sensitivity in speckle motion tracking.
Theoretically, there may exist an optimal ω 0 for the proposed method. However, the beam waist ω 0 is determined by the optics of imaging probe. Once the probe is designed and fabricated, ω 0 remains the same. Therefore, it is challenging to practically optimize ω 0 in manually scanned OCT instruments.
In our robust speckle tracking, we will calculate XCC values using multiple pairs of Ascans: I t and I t+Δt, I t and I t+2Δt , I t and I t+3Δt , …, I t and I t+NΔt , where Δt is the time interval between the acquisition of adjacent OCT Ascans. Using N different pairs of Ascans, N different values of XCC are obtained: ρ Δt , ρ 2Δt , ρ 3Δt ,… and ρ NΔt . With sufficiently small Δt and sufficiently large N, we will be able to find one XCC close enough to ρ max . We can thus choose the value of XCC closest to ρ max to calculate displacement for robust motion tracking.

OCT system
We used a spectral domain OCT (SD OCT) engine at 1.3μm for this study. More details of this system can be found in our previous work [23]. The system uses a superluminescent diode as light source (SLD1325 Thorlabs, 100 nm bandwidth, corresponding to a 7.4 μm axial resolution). Interferometric signal is detected by a CMOS InGaAs camera (SUI1024LDH2, Goodrich). A frame grabber (PCIe-1433, National Instrument) takes the interferometric signal from the camera. Results in 4.1 were obtained based on a Michelson interferometer with reference and sample arms. In the sample arm of the Michelson interferometer, the light beam is scanned by a galvanometer and a 3X scanning lens (SLM04, Thorlabs) is used in front of the sample. The e 1 lateral resolution of the system is 17 μm according to experimental characterization and therefore ω 0 = 17 μm. This configuration allowed us to accurately control the transverse motion and thus validate the proposed adaptive speckle tracking method. Data analysis was performed with Matlab. Signal to noise ratio (SNR) for images used in speckle decorrelation analysis in Section 4.1 was 43 dB, according to the following definition: SNR = 10log 10 (I mean 2 /σ 2 ). Here I mean indicates mean signal amplitude and σ 2 indicates noise variance. Results in 4.2 were obtained from a common path interferometer and signal was processed in real-time with GPU in a host computer (Dell Precision T7600) based on software developed in C + + (Microsoft Visual Studio, 2012). In the common path interferometer, a single mode fiber with flat tip was used as probe arm. Reference light derived from the fiber tip and shared the same probe path as sample light. Such a system acquired OCT signal from a simple, robust single mode fiber probe. The e 1 lateral resolution ω 0 as 43 μm for the depth range used in speckle tracking according to experimental characterization.

Experimental validation of adaptive speckle decorrelation analysis for motion tracking
To validate Eq. (1) and (4), we performed Bscan imaging on a phantom made by mixing titanium dioxide with optical epoxy. Bscan images were obtained by steering the light beam in x direction with a galvanometer scanner (Thorlabs, GVS012, 99.9% linearity and a 400 μs step response time). The galvanometer was driven by sawtooth functions with different amplitude to acquire different Bscan images at an 89 Hz frame rate. For each Bscan, we excluded the first 150 Ascans (1600 μs after the Bscan started) and used the rest of Ascans for speckle decorrelation analysis, so that the only source of motion was constant speed lateral scanning of light beam in x direction. With identical line-scan interval (Δt = 10.9 μs) of the camera, the displacement between adjacent Ascans was proportional to voltage applied to the galvanometer.
The relationship shown in Eq. (1) assumes fully developed speckle with scatterers in different spatial location described by identical but independent random variables [25,29]. To validate the statistics of speckled signal, we calculated the probability density functions (PDF) of OCT signal (Bscan). To obtain the PDF, we first converted each OCT Bscan to unitless image by normalizing its magnitude with its mean value. Afterwards, we calculated the probability for normalized, unitless signal magnitudes. Results corresponding to Bscans obtained with different scanning speeds are shown as thin lines with different colors in Fig.  3(a). In Fig. 3(a), Rayleigh distribution (thick black curve) which is the standard PDF for fully developed speckle (for imaged normalized to its mean value) is also shown [28, 29]. Our results suggests PDFs obtained from experimental data are highly consistent with standard Rayleigh PDF and the model of fully developed speckle is valid for data obtained in our experiments.
Notably, speckle statistics of OCT signal particularly higher order statistics such as correlation depends on its lateral resolution which varies as imaging depth. Therefore, we used segment of Ascan signal in the vicinity (100 μm) of beam waist to calculate PDF and perform speckle tracking. The beam waist was located 150 μm under the sample surface and 300μm from the equal optical-path-length plane. The depth range was chosen because: first, the chosen Ascan segments had high SNR, due to focused light beam, small light attenuation, and negligible sensitivity roll-off; second, our model assumed a constant beam diameter and ω 0 remained approximately constant within the 100μm segment which was significantly smaller compared to the 700 μm Rayleigh range of the incident light beam, although ω 0 generally depends on depth. In addition, ω 0 = 17 μm is valid for light beam in vacuum (or air). The beam waist in our phantom (ω phantom ) was smaller than ω 0 and could be obtained by dividing ω 0 with phantom refractive index (n = 1.4). Therefore, δx max = ω phantom / 2 = ω 0 /   2n . Given ω 0 = 17 μm, the optimal displacement should be approximately 8.6 μm, as to be shown in our experimental results. To validate Eq. (1) using the same experimental data leading to Fig. 3(a), we calculated the cross-correlation coefficient between subsequently acquired Ascans, I t and I t+Δt , in Bscan image obtained with spatial sampling interval δx 0 . The resultant XCC is denoted as ρ t (δx 0 ). We averaged ρ t (δx 0 ) at different t to obtain ρ(δx 0 ). Similarly, ρ(δx 1 ), ρ(δx 2 ), ρ(δx 3 ) … were obtained using Bscan images with spatial sampling interval δx 1 , δx 2 , δx 3 …. δx 0 , δx 1 , δx 2 , δx 3 … were in a range from 1μm to 33 μm. The results are shown in Fig. 3(b) as blue circles. Error bars were obtained using the standard deviation of ensemble values of XCC. Figure 3(b) also shows the values of XCC at different displacements according to Eq. (1) (red curve with ω 0 = 17μm). In addition, we validated the relationship between sensitivity (S) in motion tracking and displacement (δx k ) shown in Eq. (4) using experimental data. For Bscan image obtained with spatial sampling interval δx k , we compared ρ(δx k ) and ρ(δx k + dδx) obtained from Ascan pairs with displacements δx k and δx k + dδx, and calculated the ratio between [ρ(δx k + dδx)ρ(δx k )] and dδx for the evaluation of S. The results are shown in Fig. 3(c) (circle: experimental data; red curve: theory). Figure 3  According to discussion in Section 2 and results shown in Fig. 3, it is ideal to select Ascans with appropriate displacement (δx max = ω 0 /   2n and appropriate degree of decorrelation (ρ max = e -1/2 ) for speckle tracking. To further validate this, we performed speckle decorrelation analysis using data from a Bscan with 0.37 μm lateral scanning interval. We calculated XCC between Ascans with small time interval (Δt, 2Δt, and 3Δt) and thus small lateral displacements (0.37 μm, 0.74 μm and 1.11 μm). We also calculated XCC between Ascans with medium time interval (21Δt, 22Δt, and 23Δt) and thus medium lateral displacements approximating δx max = 8 μm (7.73 μm, 8.09 μm and 8.47 μm); as well as XCC between Ascans with large time interval (45Δt, 46Δt, and 47Δt) and thus large lateral displacements (16.56 μm, 16.93 μm and 17.30 μm). In this analysis, Ascans with large displacement (17 μm) were obtained at a 480μs time interval. This time interval was small enough so that the decorrelation was induced by speckle rather than other random factors. XCC values obtained are shown in Fig. 4(a), 4(b) and 4(e), corresponding to small, appropriate and large degree of decorrelation. Using XCC values obtained, we extracted displacements using Eq. (3) and show results in Fig. 4(d), 4(e) and 4(f). Figure 4(d) shows that the estimated magnitude of motion is larger than its actual value if the XCC values are close to 1. With small displacement and thus XCC value close to 1, the decorrelation is dominated by random noise in OCT signal, resulting in an overestimated motion. On the other hand, motion is underestimated in Fig. 4(f) where Ascans with large displacement are used to calculate XCC. This is because the expected value of XCC is extremely small (Eq. (1), Fig. 1) when δx is large and the estimated displacement is determined by the non-zero noise floor in XCC measurement. In comparison, displacements obtained using XCC values close to s ρ max (ρ max = e -1/2 ) (Fig. 4(b) and 4(e)) are consistent with actual displacements. Results in Fig. 4  show that it is critical to use Ascans with appropriate degree of decorrelation for motion tracking. In practice, motion of sample or probe in OCT imaging can have speed varying in a large range and two Ascans acquired with a fixed time interval can have very small or very large displacement. According to results shown in Fig. 3 and Fig. 4, XCC between subsequently acquired Ascans can be close to 1 or close to 0, corresponding to a limited accuracy in motion tracking. Therefore, adaptive speckle decorrelation analysis that estimates displacement using Ascans with XCC close to ρ max can achieve more robust speckle tracking over a large range of motion speed.
To demonstrate the advantage of the adaptive motion tracking algorithm, we compared the results obtained from different speckle decorrelation analysis strategies for motion tracking, as shown in Fig. 5. We used OCT data in Bscans obtained with different lateral scanning speeds in our analysis and the horizontal axis in Fig. 5 indicates lateral scanning speeds.
For Bscan acquired with a certain lateral scanning speed, we calculated XCC values (averaged) using Ascan pairs with different time intervals (Δt, 2Δt, 4Δt, 8Δt, and 16Δt) and hence displacement between Ascan pairs shown as circles in different colors in Fig. 5. Notably, when Ascans with time interval kΔt were used for speckle tracking, the resultant XCC and the lateral displacement (δx(kΔt)) corresponded to time interval kΔt. Therefore, the displacement between subsequently obtained Ascans was calculated by dividing δx(kΔt) with k. On the other hand, we performed adaptive decorrelation analysis using the same set of experimental data. For each Bscan acquired with a certain lateral scanning speed, we calculated XCC values (averaged) using Ascan pairs with different time intervals (Δt, 2Δt, 3Δt, 4Δt…., 16Δt) and denoted the results as ρ Δt , ρ 2Δt , ρ 3Δt , ρ 4Δt , …, ρ 16Δt . We then chose the XCC value (ρ mΔt ) closest to ρ max = e -1/2 to calculate the displacement δx(mΔt) between two Ascans acquired with time interval mΔt and the displacement between adjacent Ascans which was δx(mΔt)/m. Results obtained from adaptive speckle decorrelation analysis are shown as black triangles in Fig. 5. To evaluate the performance of different speckle tracking methods, we also show the actual displacement between adjacent Ascans in Fig. 5 (solid black line) according to the known scanning speed and Ascan acquisition rate. As shown in Fig. 5, displacement obtained from XCC between Ascans with fixed time interval deviates from its actual value, no matter the time interval is small or large. Displacement estimated by calculating XCC between Ascans with small time intervals (Δt and 2Δt) is larger than its actual value for small speed. On the other hand, displacement obtained by calculating XCC between Ascans with large time intervals (4Δt, 8Δt and 16Δt) is smaller than its actual value for large speed. This is also consistent with results shown in Fig. 4(d) and 4(f). In comparison, adaptive speckle tracking provides a closer approximation to actual displacement over a large range of motion speed and outperforms speckle tracking based on XCC between Ascans with fixed time intervals.
From a different perspective, inaccuracy in displacement quantification for small or large decorrelation (XCC1 or XCC0) is due to the small value of motion tracking sensitivity, S. As shown in Fig. 3(c), when displacement is very small (XCC1) or value large (XCC0), S is approximately 0. With extremely small S, change in displacement does not lead to measurable change in XCC and therefore the estimated displacement appears saturated. To further validate that our adaptive speckle decorrelation analysis can track motion with nonconstant speed, we applied sinusoidal voltage with period T (T = 5.5 ms) to the galvanometer scanner and acquired OCT signal within a time span of 2T. We stacked sequentially acquired Ascans and show the resultant 2D data in Fig. 6(a). Due to the sinusoidal voltage applied to the galvanometer, the speed of lateral scanning was not constant. Therefore, distortion artifacts are clearly visible in Fig. 6(a). To obtain the magnitude of speed by analyzing acquired OCT signal, we performed adaptive speckle decorrelation analysis to derive the speed of motion and show the result as dashed black curve in Fig. 6(b). For comparison, we used Ascans acquired subsequently with time interval Δt for XCC calculation and motion tracking, with results shown the red dashed curve in Fig. 6(b). We also used Ascans pairs acquired with larger time interval (10Δt) for motion tracking, with results shown the green dashed curve in Fig. 6(b). The ground truth scanning speed is plotted as the solid black curve for comparison. Results of motion tracking shown in Fig. 6(b) correspond to two periods of the sinusoidal voltage and four speed peak can be observed because speckle decorrelation analysis extracts the magnitude of lateral motion that can have positive or negative values. Clearly, adaptive speckle decorrelation analysis provides a better approximation of the beam scanning speed.
Using scanning speeds obtained from different speckle decorrelation analysis strategies, we re-sampled Ascans shown in Fig. 6(a) to correct motion artifact due to nonconstant beam scanning speed. Reconstructed images are shown in Fig. 6(c), 6(d) and 6(c), based on adaptive speckle decorrelation analysis, speckle decorrelation analysis of Ascans with small time interval (Δt = 10.9 μs), and speckle decorrelation analysis of Ascans with large time interval (10Δt = 109 μs). We also used the actual sinusoidal beam scanning speed to correct image artifact and obtained Fig. 6(f). Compared to Fig. (d) and (e), Fig. 6(c) better resembles Fig. 6(f), suggesting that more robust motion tracking through adaptive speckle decorrelation analysis can provide better correction of image artifact.

Real-time distortion-free manual-scanning OCT imaging based on adaptive speckle decorrelation analysis
To demonstrate the effectiveness of adaptive speckle tracking in real-time imaging, we developed software to reconstruct distortion-free OCT images using signals obtained from a manually scanned probe. The algorithm is described in detail below. With M spectral interferograms acquired as a frame by the frame grabber, we subtracted reference spectrum from the measured interferograms, converted data from wavelength space to wavenumber space and performed fast Fourier transform (FFT) to obtain M Ascans. Afterwards, adaptive speckle decorrelation analysis illustrated in Fig. 7 was performed. For the i th Ascan (I i ) in the frame acquired at time t, we used Eq. (2) to calculate XCC values between I i and subsequent Ascans, I i + 1 , I i + 2 , I i + 3 , …, acquired at t + Δt, t + 2Δt, t + 3Δt, …, with Δt indicating the time interval between the acquisition of adjacent Ascans (Fig. 7(a)). Resultant XCC values, ρ i,Δt , ρ i,2Δt , ρ i,3Δt ,…, were expected to be different as illustrated in Fig. 7(b), according to our theoretical analysis and experimental results shown in Fig. 1 and Fig. 3(b). We then selected XCC that was closest to ρ max (dashed red line in Fig. 7(b)) for robust motion tracking, as indicated by the red arrow in Fig. 7(b). Denote the data point used for displacement calculation as ρ i,mΔt . The lateral displacement between the i th Ascan A i and the (i + m) th Ascan I i + m could be calculated using Eq. (3): δx i,m = ω 0 [ln(1/ρ i,m )] 1/2 . Assuming constant motion during the acquisition of I i , I i + 1 , I i + 2 , I i + 3 , I i + 4 …, I i+k , we were able to estimate the displacement between the i th Ascan and the subsequent Ascan: δx i = δx i,m /m. Procedures for XCC calculation and displacement estimation shown in Fig. 7(a) and 7(b) were performed for all M Ascans obtained in the frame, to obtain displacements between adjacent Ascans, δx 1 , δx 2 , …, δx M-1 . Due to the nonconstant beam steering speed in manual scanning, δx 1 , δx 2 , …, δx M-1 had different values, as shown in Fig. 7(c). For distortion-free OCT imaging, we established a grid to place Ascans with identical interval δx s . Signal values of the grid were obtained through interpolation using raw Ascan data and actual displacement values (δx 1 , δx 2 , …, δx M-1 ) derived from adaptive speckle decorrelation analysis. Given spatially oversampled OCT signal, we were able to obtain M decorr (M decorr <M) Ascans to form a distortion-free 2D image (I 0 ). We then attached I 0 to the end of the existing 2D distortion-free image and update the display. The above-described software was implemented in real-time, summarized as flow chart shown in Fig. 8(a). To estimate the motion speed between adjacent Ascans, we calculated 16 values of XCC using Ascans with different time intervals. Procedures indicated by red boxes were implemented using GPU. Signal processing speed was analyzed using NVidia Nsight embedded in Visual Studio. Ascans were processed on average at a speed of 120 kHz which is faster than the highest data acquisition speed of the line scan camera used in our OCT system (92 kHz). The most time consuming computation tasks in our algorithm are shown Fig. 8(b).  To validate the robustness of adaptive speckle tracking in real-time imaging, we made a phantom (photo shown in Fig. 9(a)) by covering 1 line/mm Ronchi Ruling Pattern with multiple layers of Scotch tapes and performed manual scanning using a single mode fiber OCT probe. Signal was processed in real-time using our home-made GPU software. Motion was limited in x dimension. To prevent axial motion (z dimension) from decorrelating the signal, we maintained the probe at the same height from the phantom surface during imaging. In addition, we constrained the motion with a barrier along x direction to exclude motion in y direction. We compared images reconstructed using XCC adaptive that was obtained through adaptive speckle decorrelation analyses, with images reconstructed using XCC between adjacent Ascans (XCC adj ). Notably, signal within depth range corresponding to Scotch tape layers was used to calculate XCC. Otherwise, the periodical lateral structure from the ruling pattern might affect the accuracy of speckle tracking. Images obtained using XCC adaptive and XCC adj are shown in Fig. 9(b) and 9(c). Clearly, Fig. 9(b) shows well defined periodical structure while residual distortion exists in Fig. 9(c). Distortion artifact can be seen more clearly within the dashed rectangular box in Fig. 9(c). To quantitatively compare the accuracy of motion tracking based on XCC adaptive and XCC adj , we averaged OCT signal within depth range indicated by the green arrow in Fig. 9(b). The result was subtracted with its mean value ( Fig. 9(d) and 9(e), blue curve) and low pass filtered ( Fig. 9(d) and 9(e), black, dashed curve). We then identified zero-crossing points corresponding to the edge of the bars in the ruling pattern, as indicated by red circles in Fig. 9(d) and 9(e). For clear visualization of zerocrossing points, Fig. 9(d) and 9(e) show results of edge detection within lateral range indicated by the red arrow in Fig. 9(b). Figure 9(d) shows a more uniform distribution of zerocrossing points compared to Fig. 9(e). The actual distance between adjacent zero-crossing points is the same for different bars in the periodical ruling pattern and a more uniform distribution of zero-crossing points suggests better performance in motion artifact removal. Therefore, results in Fig. 9(d) and 9(e) show that motion tracking based on XCC adaptive outperformed motion tracking based on XCC adj . We repeated the manual scanning experiments for four times, and performed motion tracking using XCC adaptive and XCC adj , respectively. To evaluate the quality of reconstructed image in terms of motion artifact, we performed zero-crossing detection on OCT signal and obtained the number of Ascans (N Ascan_bar ) between adjacent zero-crossing points. N Ascan_bar had different values for different pairs of zero-crossing points and the robustness of motion tracking was quantified using the ratio between the standard deviation and the mean of N Ascan_bar , denoted as R adaptive and R adj , for images reconstructed using XCC adaptive and XCC adj . R adaptive and R adj were averaged using data obtained from all the four experiments and turned out to be 0.11 and 0.27, indicating a significant improvement in motion artifact removal through adaptive speckle decorrelation analysis.  The advantage of adaptive speckle analysis is further illustrated in videos acquired in realtime using our GPU software. Visualization 1 and Visualization 2 were obtained by imaging the phantom shown in Fig. 9(a) with a manually scanned single mode fiber probe, based on motion tracking with XCC adaptive and XCC adj . Screen captures of the software are shown in Fig. 10(a) and 10(b). During the acquisition of Visualization 1 and Visualization 2, we first scanned the single mode fiber probe slowly and then with higher speed. Due to improved robustness through adaptive speckle decorrelation analysis, Visualization 1 can reveal the periodical structure of the resolution target in the reconstructed 2D image. In comparison, Visualization 2 obtained using speckle analysis between adjacent Ascans shows distorted 2D image due to limited dynamic range in motion tracking. We performed manual scanning using a single mode fiber probe on the skin of a healthy volunteer and obtained distortion-free images through adaptive speckle decorrelation analysis. Images were obtained from the skin of forearm ( Fig. 11(a) and Visualization 3), fingertip ( Fig. 11(b)), volar of hand (Fig. 11(c)) and dorsum of hand ( Fig. 11(d)). Scale bars in Fig. 11 indicate 250μm. Skin layers (stratum corneum, epidermis, dermis, subcutaneous adipose tissue) are clearly visible. Figure 11(b) and 11(c) also show sweat ducts as indicated by green arrows. In Fig. 11(d), blood vessels that cause image shadow can be observed, as indicated by white arrows.

Conclusion and discussion
In this study, we developed an adaptive speckle decorrelation analysis strategy. We calculated XCC using Ascans with different time intervals and selected the XCC value that maximizes motion tracking sensitivity for displacement calculation. This method was validated using experimental data. We implemented the algorithm for robust motion tracking in GPU software to reconstruct distortion-free image for a manually scanned OCT system. Images were obtained from phantom and human skin.
The motion tracking method described in this manuscript enables scanner-free OCT imaging based on a manually actuated miniature single mode fiber probe. Such an imaging configuration offers great flexibility to access pathological region of interest in clinical diagnosis and surgical guidance. It can achieve a large lateral field-of-view (FOV) that is independent of lateral resolution. This is a desirable when imaging region of interest has a large dimension, such as in OCT's application for tumor margin delineation. Conventionally, Bscan OCT image is obtained by steering collimated light beam with galvanometer in front of an imaging lens. In such a sample arm configuration, the lateral FOV limited by the dimension of imaging optics. To achieve a large lateral FOV at a given lateral resolution, the dimension of the imaging lens has to be large. To demonstrate this, we performed Bscan imaging on the phantom used in section 4.2 using conventional sample arm configuration with a 10X scanning lens (Thorlabs, LSM02, 13 μm 1/e 2 beam waist size). The resultant image is shown in Fig. 12 with a 3.4mm lateral FOV. The scale bar in Fig. 12 represents 500 μm. The surface of the phantom appears curved, because light beam had large aberration when passing through the edge of the scanning lens. In comparison, Fig. 9(b) and 9(c) are free of such image artifact due to beam aberration and have a significantly larger FOV compared to Fig. 12.   Fig. 12. Bscan image obtained from the phantom shown in Fig. 9(a), using a sample arm configuration based on galvanometer and scanning lens. Scale bars indicate 500μm.
Speckle tracking developed here provides quantification for the magnitude of motion. Although motion in 3D space is essentially a three dimensional vector, our speckle decorrelation analysis can provide satisfactory performance in tracking motion with limited degree of freedom. When OCT signal is decorrelated by motion with non-zero components in dimensions other than x direction, it is challenging to fully correct motion artifact through speckle tracking. However, 2D image consisting of re-sampled Ascans will appear undistorted. This allows morphological feature of the sample to be identified through visual inspection. Otherwise, multiplicative speckle noise randomly modulates the magnitude of OCT signal and morphological features are not discernable in simple Ascan display or 2D image obtained by directly stacking sequentially acquired Ascans. For the application of speckle tracking, it is critical to achieve a large dynamic range in motion speed measurement. Particularly, the adaptive speckle tracking in this paper allows more accurate tracking of small motion speed, compared to the method we developed before based XCC between adjacent Ascans. If Ascans acquired subsequently with time interval Δt are used in decorrelation analysis, the smallest speed that can be reliably measured is inversely proportional to Δt: v min = D min /Δt. Here, D min indicates the smallest displacement that can be reliably measured. As a result, it is challenging to track slow motion using an OCT system with high imaging speed (small Δt). In comparison, the adaptive speckle tracking methods calculates XCC using Ascans with larger time interval (kΔt where k is larger than 1) and can be sensitive to slow motion. On the other hand, the maximum motion speed that can be estimated accurately by the adaptive speckle tracking method is still determined by OCT system's data acquisition rate. If the scanning speed is so large that the displacement within Δt is significantly larger than the beam waist size ω 0 , subsequently acquired Ascans are completely decorrelated and calculating XCC using adjacent Ascans or Ascans with larger time interval cannot provide effective quantification of motion. However, current high-speed OCT system allows us to track a large scanning speed and does not pose a significant limitation in practice. Assume an OCT system has a 92 kHz Ascan rate (R) and the lateral interval between adjacent Ascans is δx max = 10 μm. The highest trackable speed is thus 0.46 m/s (v = δx max R), and is a very large scanning speed for the application of manual scanning OCT imaging.
Intuitively, there exists an optimal value of XCC between 0 and 1 for speckle decorrelation analysis. If Ascans involved in XCC calculation are completely decorrelated, the resultant XCC is small but does not carry any information about displacement. If Ascans involved in XCC calculation are almost identical, the resultant XCC value is dominated by additive noise in OCT signal rather than spatial displacement. In this study, we chose to select XCC value that maximized S defined in Eq. (4) to achieve robust motion tracking. We successfully demonstrated adaptive motion tracking using XCC that maximized. On the other hand, XCC measurement that suffers from uncertainty that can be quantified by σ xcc (standard deviation) [30]. An alternative strategy is to perform speckle tracking using XCC values that optimize the ratio between S and σ xcc : (η = S/σ xcc ). However, σ xcc is generally unknown for an arbitrary sample, therefore optimizing η is an untraceable problem.