Temporally and spatially adaptive Doppler analysis for robust handheld optical coherence elastography.

Optical coherence elastography (OCE), a functional extension of optical coherence tomography (OCT), can be used to characterize the mechanical properties of biological tissue. A handheld fiber-optic OCE instrument will allow the clinician to conveniently interrogate the localized mechanical properties of in vivo tissue, leading to better informed clinical decision making. During handheld OCE characterization, the handheld probe is used to compress the sample and the displacement of the sample is quantified by analyzing the OCT signals acquired. However, the motion within the sample inevitably varies in time due to varying hand motion. Moreover, the motion speed depends on spatial location due to the sample deformation. Hence, there is a need for a robust motion tracking method for manual OCE measurement. In this study, we investigate a temporally and spatially adaptive Doppler analysis method. The method described here strategically chooses the time interval (δt) between signals involved in Doppler analysis to track the motion speed v(z,t) that varies temporally and spatially in a deformed sample volume under manual compression. Enabled by temporally and spatially adaptive Doppler analysis, we report the first demonstration of real-time manual OCE characterization of in vivo tissue to the best of our knowledge.


Introduction
Optical coherence tomography (OCT) allows structural and functional imaging of biological tissue with high resolution and high speed [1]. The imaging capability of OCT can be integrated into handheld instruments using fiber optic components [2][3][4]. A compact, lightweight handheld OCT probe allows a clinician to interrogate tissue characteristics at different anatomical locations [5,6]. Therefore, handheld OCT imaging instrument is attractive for delineating tu breast or pros to reveal mo characteristics extracted. On such as breast normal tissue been used in elastography mechanical pr frequency, el mechanical ex imaging and spatial resolut Despite gr depth resolve surrogate for remains const the sample (ε sample). Ther derivative of smaller for a pathophysiolo displacement has a larger excision with generated thro depth, and th boundary be characterizatio informed clin tissue. Fig palpation of the tissue with great convenience and flexibility. Moreover, fiber optic OCE instruments can be integrated into a needle device, further delivering the capability of mechanical characterization to tissue that is deeply embedded. However, the major challenge for manual OCE characterization of tissue is the unpredictable and unstable motion within the tissue. The deformation of the sample under the known pattern of mechanical excitation can be reliably tracked by analyzing OCT signal. In conventional compression OCE measurement, the sample has a well-defined geometry and undergoes quasi-static compression. Alternatively, the mechanical excitation can be impulse or sinusoidal function [22,23]. In previous studies of 1D or 3D OCE for the mapping of mechanical properties, an actuator was often used to deform the sample. However, the need to use an actuator for mechanical loading limits the development of a light-weight and compact OCE probe. For a handheld, manually actuated OCE instrument, it is challenging to impose mechanical excitations that are quasi-static, impulse or sinusoidal. The manual loading process often generates a motion speed that varies with time. The quality of in vivo OCE signal is also affected by involuntary motion from the subject and from the user who holds the probe. In addition, the sample deforms under compression, implying spatial variation of motion characteristics. Our software approach of adaptive Doppler analysis enables 1D OCE measurement based on a manual loading process, and allows the acquisition of high quality in vivo OCE signal from a handheld probe.
Motion tracking in OCE can be achieved through Doppler analysis or speckle decorrelation analysis. Speckle analysis has a smaller dynamic range and is more appropriate to track motion with larger magnitude [24,25]. In this study, Doppler analysis is used to quantify the axial motion speed and displacement. A simple and effective method for temporally and spatially adaptive Doppler analysis is investigated here. The adaptive Doppler analysis method strategically chooses the time interval (δt) between signals involved in Doppler analysis, to track the motion speed v(z,t) that varies temporally in a manual compression process and spatially in a deformed sample volume. The method is validated in an OCE system with a handheld single fiber probe and real-time signal processing software based on graphic processing units (GPU). To achieve robust motion tracking, we calculate high density (HD) Doppler phase shift that is most unlikely to have phase wrapping artifact and average the HD Doppler signal to estimate the speed of axial motion from which we derive a time interval to achieve a large yet artifact free Doppler phase shift. The premise of this method is that (1) directional motion affects larger scale characteristics of the Doppler signal and can be estimated through averaging; (2) noise characteristics in estimated Doppler phase shift are independent of the time interval δt while the signal due to directional motion does. Enabled by high signal acquisition and processing speed, we perform an online estimation of the motion speed, select an optimal δt adaptively, and perform robust motion tracking for OCE measurement.
The manuscript describes the first handheld fiber optic instrument that allows in vivo realtime OCE characterization, to the best of our knowledge. The manuscript is organized as follows. First, we introduce the principle of the adaptive Doppler analysis method. Afterwards, the imaging system and data acquisition are described. We then show results obtained from phantom experiments and in vivo tissue characterization, to demonstrate the effectiveness of the adaptive Doppler analysis for motion tracking in a dynamic manual loading process.
It is worth mentioning that we do not intend to extract quantitative mechanical properties of the sample through OCE measurement in this study. Instead, we acquired one dimensional data to reveal depth resolved sample displacement. By observing the variation of the displacement, particularly the slope, mechanical contrast of the sample can be revealed. The spatial variation of stress, 3D nature of displacement, as well as viscoelastic behavior of the sample are not considered. Moreover, OCE is referred to the application of OCT technology for mechanical characterization. In this manuscript, OCE signal indicates depth resolved displacement of the sample under compression, although the strain (derivative of displacement) is more directed related to the mechanical properties of the sample.

Principle for adaptive Doppler analysis
In OCE characterization, the loaded sample deforms and has displacement dependent on spatial location (δL(z)). With local displacement δL(z) extracted by analyzing OCT signal, localized axial strain, the spatial derivative of the displacement (Eq. (1)), is calculated as the surrogate of sample stiffness. For sampled discrete OCT image, local strain can be estimated either through finite difference approach or least square estimation.
It is worth mentioning that the motion within a deformed sample under axial compression is generally 3D with axial and lateral components. However, Doppler phase analysis is only sensitive to axial motion. Therefore, motion in transverse plane is not considered in this study, as in most previous OCE studies based on Doppler analysis. Furthermore, our measurement geometry has cylindrical symmetry, and the light beam propagates along the axis of cylindrical symmetry. Hence the lateral displacement of an isotropic sample seen by the incident light beam is expected to be minimum.
To obtain 1D depth resolved OCE signal (δL(z)), Doppler phase shifts between OCT Ascans are calculated. Consider the OCT signal with complex value at the k th pixel at depth kδz of an A-scan (m th A-scan) and that at the k th pixel of another A-scan ((m + Δ k,m ) th A-scan).
The relationship between the estimated Doppler phase shift δ ,k m φ and the actual phase shift δφ k,m due to directional motion is shown in Eq. (4), where n k,m is the random phase noise deriving from various noise sources in OCT measurement (shot noise, thermal noise, excess noise, speckle noise, etc). On the other hand, N k,m is an integer and is non-zero when |δφ k,m |>π/2: N k,m = .
On the right hand side (RHS) of Eq. (5), the first term represents the actual displacement; the second term represents the phase wrapping artifact and the third term denotes the contribution from random phase noise. To improve the sensitivity, SNR and dynamic range for OCE characterization, it is desirable to have a smaller variance (Var(δ  k L -δL k )) for the estimated displacement. It is assumed that n k,m (m = 1, 2, 3, …, M) is Gaussian, and independent in different A-scans. The variance of n k,m is shown in Eq. (6), where SNR is the signal to noise ratio of the OCT signal and β is a constant [30-32]. With the above assumption and with N k,m ≡0, the variance in displacement tracking is given by Eq.
In Eq. (7), λ 0 depends on the OCT system used for the imaging study; M depends on the time period of the sample loading process; Var(n k,m ) is determined by the OCT system as well as the optical characteristics of the sample. Hence Δ k,m is the only parameter that can be varied to improve the displacement tracking, and noise in displacement tracking for a given compression process can be reduced with a larger value of Δ k,m .
On the other hand, for unbiased displacement tracking, it requires N k,m ≡0. Therefore, the time interval between A-scans used for Doppler phase calculation (δt = Δ k,m T 0 ) has to be sufficiently small, such that |δφ k,m | = |4πv k,m Δ k,m T 0 /λ 0 |≤ π/2. To prevention phase wrapping artifact in Doppler analysis, the following condition has to be satisfied.
Equation (8) implies the optimal choice of Δ k,m depends on the speed of the motion (v k,m ). For a handheld OCE instrument used to exert manual compression, the motion speed within the sample depends on the depth because the sample is axially deformed. The motion speed also varies as time due to the non-constant compression speed. Therefore, v k,m = v(kδz,mT 0 ) and there is a need to have a time interval adaptive to the spatial location and time (δt(z,t) = Δ(kδz,mT 0 )T 0 = Δ k,m T 0 where Δ k,m is an integer) for Doppler analysis. In other words, different values are chosen for Δ k,m at different depth (z = kδz) and at different time (t = mT 0 ) ( Fig. 2  |≤ 2 π to prevent phase wrapping from happening. For discrete OCT signal, δt = δt k,m = Δ k,m T 0 where the integer Δ k,m can be obtained by (9) with W>1. W is a coefficient that we introduce to scale the magnitude of phase shift δφ k,m with regard to π/2, to make sure that δφ k,m is large enough while free of phase wrapping artifact. Δ k,m is assigned to have a value of 1, if the result obtained by Eq. (9) is smaller than 1. Moreover, to calculate phase shift using A-scans within the same frame of OCT data, it requires Δ k,m to be smaller than M 0 /2. If the value calculated using Eq. (9) is larger than Notably, we choose the value of W to be larger than 1, such that the method is robust against phase wrapping when phase noise exists. As validated in previous studies including our recent work [30-32], the level of phase noise in the OCT imaging system is inversely proportional to the signal to noise ratio (SNR) of amplitude OCT signal. Consider a shot noise limited OCT system with noise level determined by the power of reference light. For such an OCT system, the phase noise is small for a sample that generates large amplitude OCT signal and the value of W can be close to 1. In comparison, the phase noise is large for a sample that generates small amplitude OCT signal and the value of W has to be sufficiently larger than 1. In this study, other than specifically mentioned, W = 2 for the calculation of adaptive time interval for Doppler analysis.
Using adaptively determined time interval for Doppler analysis (δt k,m = Δ k,m T 0 with Δ k,m obtained from Eq. (9) . By calculating the displacement within the entire compression process, M Doppler phases are averaged and the value of M is several orders of magnitude larger than 1. Therefore, the displacements obtained are orders of magnitude larger than the wavelength of the light source (Figs. 6,8,[9][10][11]. With depth resolved displacement, we then estimate the depth resolved strain of the loaded sample to evaluate its stiffness.
In summary, we have implemented the adaptive Doppler analysis illustrated in Fig. 3 in real-time through GPU accelerated parallel computation. The software acquires spectral interferograms frame by frame, performs fast Fourier transform on the spectral interferograms, calculates the HD Doppler phase shift to estimate the speed of axial motion, adaptively determines the optimal time interval for each frame of OCT data to perform Doppler analysis, and tracks depth resolved displacement for sample mechanical characterization.

Experime
The

Results
We used NVI software proc was much fa second). Mor OCT data was To demon manual comp in Fig. 4 Fig. 4(a).   with v motor = 0.1mm/s) are free of phase wrapping artifact for the same δt. Figure 5(a) suggests that the selection of time interval for Doppler analysis has to be adaptive to the mechanical excitation, i.e., the translation speed of the compressor during OCE characterization. In addition, the phase calculated using Eq. (3) also fluctuates randomly due to noise. Using  j δφ (δt) obtained, we assessed the random noise of the estimated Doppler phase σ φ (δt) = for phase obtained with different time intervals (δt = ΔT 0 ). σ φ (δt) are shown in Fig. 5(b) as blue (v motor = 0.2mm/s) and red (v motor = 0.1mm/s) curves. A peak can be observed in the blue signal in Fig. 5(b), because Doppler signal varies drastically when phase wrapping appears (blue signal in Fig. 5(a)). Other than the peak due to phase wrapping, the noise in Doppler phase estimation remains approximately constant for different values of δt. This is because random phase variation in OCT signal originates from factors (noise in OCT measurement and random environmental perturbations) that can be considered as temporally independent and identically distributed random variables, as indicated by Eq. (6). Therefore, the results of Doppler analysis are expected to have a similar level of noise, despite different time interval of δt. According to Eq. (7), a larger value of Δ k,m is desirable to achieve a reduced error in displacement tracking, because the phase noise does not increase with time ( Fig. 5(b)) while the phase shift due to directional motion increases with time ( Fig. 5(a)).
In addition, the displacement within the sample under OCE characterization varies as spatial location due to the sample deformation under compression. The deformation is quantified as axial strain (Eq. (1)) to reveal the mechanical properties of the sample. Therefore, Doppler analysis also has to be adaptive to the spatial location. To demonstrate this, one frame of OCT data acquired in the above described experiment (v motor = 0.2mm/s) was analyzed. We calculated Doppler phase shift between pixels at depth z = 1mm, as well as Doppler phase shift at a smaller depth (z = 0.5mm). The mean Doppler phase shifts for different δt are shown in Fig. 5(c) as blue (z = 1mm) and red (z = 0.5mm) curves. Similar to results shown in Fig. 5(a), the Doppler phase shift in Fig. 5(c) initially increases linearly with δt. For Doppler phase shift calculated for a larger depth (blue curve in Fig. 5(c) corresponding to z = 1mm), phase wrapping artifact arises when δφ approaches and exceeds π/2. In comparison, data obtained from a smaller depth (red curve in Fig. 5(c) with z = 0.5mm) are not affected by phase wrapping artifact for the same range of time interval. Therefore, it is necessary to select time interval adaptive to spatial location for Doppler analysis for robust Doppler analysis. Using the same set of OCT data obtained with v motor = 0.2mm/s, we also evaluated the random noise for Doppler phase shift for different depths. The results are shown in Fig. 5(d) as blue (z = 1mm) and red (z = 0.5mm) curves. Despite a peak observed in the blue curve due to the phase wrapping, the Doppler signals show a constant noise level for different values of δt. for Doppler an s. When the m ) to compress is because the larger than π/2 gest that the ti motion. In a 4(a), we trans acquired OCT g. 6(b). Notabl tween A-scans 00, 150. In Fi ed with a sma increases wit is larger for a   Fig. 6(b) ve Doppler an a handheld y compressed b were obtained ive (Eq. (9)), l e interval betw he displacemen 1 (adaptive tim displacement fi placement trac t does not inc o the phase wr clearly in Vis = T 0 (black sign displacement sy and free fro alysis is crucial maneuver is i nd 0.9mm/s), a Depth resolved own in Fig. 9 z and ε = δL mo s in Fig. 9(b)  We furthe analysis coul prepared (Sam cellophane ta manually, and the displacem 11(b)). In Fig  Sample 1 and mechanically depth due to (blue signal i with depth. In the tape laye obtained from reaching the layers did no analysis allow motion within location.  Fig. 11 terization base nical propertie con phantom w eld probe to c er analysis on 1(a)) and Sam nitude OCT sig antom (Sample ed signal in F hile the depth oppler analysis obtained from S pth in Fig. 11 ) is compresse alysis is show ( Fig. 12(b)) as the red curv or healthy and he displacemen nt for anothe arison, the disp With the assum ity of the disea within the dep he fingertip (Fi were obtained Fig. 12

Conclusion and discussion
In summary, we developed and validated a Doppler analysis method that is adaptive to time and spatial location, for robust manual OCE characterization based on a handheld instrument. Enabled by this adaptive Doppler tracking strategy, real-time, manual mechanical characterization of in vivo tissue was demonstrated for the first time to the best of our knowledge.
It is worth mentioning that our goal is to reveal local mechanical contrast of the tissue (Fig. 11(b), Fig. 12(c) and 12(f)). Absolute measurement of tissue mechanical properties in vivo is a much more challenging task. Quantification of tissue mechanical property requires to know the 3D spatial distribution of the stress and the strain. In addition, the measurement boundary condition determined by the probe as well as the sample structure/heterogeneity has to be known. Moreover, for biological tissue that is generally viscoelastic, the force and displacement generated in the compression process depends not only on the intrinsic properties of the tissue but also on the dynamic loading process. However, quantitative measurement of tissue mechanical properties is beyond the scope of this study. Doppler analysis only tracks motion in axial direction. Therefore, we simply used the slope of axial displacement to represent the magnitude of sample deformation. We also assumed the stress to have a uniform spatial distribution. Hence the displacement has a linear dependency on the depth within a mechanically homogeneous volume from which OCT signal is acquired. This is validated by our experimental data (Fig. 9). Currently, we neglect the viscoelastic behavior of the sample. Without force/stress quantification, the strain was used as a surrogate for the stiffness of the sample. Despite the above simplifications, our measurement has the capability to reveal local mechanical heterogeneity, as shown in Fig. 11(b), Fig. 12(c) and Fig. 12(f).
To improve the dynamic range of Doppler tracking, previous effort has been focused on tracking fast motion, through unwrapping algorithm or exciting the sample with a crawling wave [34]. However, when the Ascan acquisition rate was high and the manual compression speed was slow, it is essential to track small Doppler phase shift that might be overwhelmed by random phase noise. Therefore, our method improves the dynamic range of motion tracking mainly by using a longer time window (Δ k,m >1) to track slow motion. Consider the smallest measurable motion speed (v min ) to be equivalent to the noise magnitude in speed estimation. With Eqs. (6) and (7), we can estimate the minimal speed for the compression process, and v min takes its smallest value when Δ k,m takes its largest value (Δ k,m = M 0 /2).
Assuming constant SNR, we have: , min v is approximately 0.003mm/s and max v is approximately 120mm/s, which provides a sufficient dynamic range to track manual compression process.
The adaptive motion tracking method described here uses a longer time window (Δ k,m >1) to observe the change of signal by calculating the Doppler phase shift. This allows more accurate quantification of the motion but also reduces the temporal resolution of the measurement. However, with the assumption of quasi-static compression on elastic sample, we only consider the final accumulated displacement within the sample. Therefore, the compromised temporal resolution does not affect the characterization of local mechanical heterogeneity.
The potential benefit of using non-adjacent Ascans for Doppler analysis was noted by Adie et al [35]. However, in [35], there lacked detailed discussion and experimental validation for adaptive Doppler analysis. The authors of [35] either varied the data acquisition rate to track motion within the sample at different excitation frequencies, or chose sampling interval to be 1 or 2 depending on the known excitation frequency. The study described in this manuscript offers a practical solution for robust OCE characterization through a simple handheld device.
Compared to Eq. (3) that extracts Doppler phase using two-quadrant arctangent, the maximum phase free of phase wrapping artifact can be doubled by using four-quadrant arctangent function. Two-quadrant arctangent function was adopted in this study to simplify algorithm implementation for the analysis of noisy signal. Consider a complex OCT signal I(t). To obtain Doppler phase shift δφ within a time interval δt, Z(t) = I(t + δt)I*(t) = X + jY is calculated and δφ is estimated using the inverse of the tangent function:  δ φ = tan −1 (Y/X). The inverse of the tangent function can be obtained using either two-quadrant arctangent atan(Y/X), or four-quadrant arctangent (atan2(Y,X)) that calculates atan(Y/X) and uses the signs of both arguments to determine the quadrant of the resultant phase. Four-quadrant arctangent function is based on two-quadrant arctangent and hence has similar noise characteristics. When X is small and Y is significantly larger, the absolute value of atan(Y/X) is close to π/2 and can be overwhelmed by noise, because the small value (X) in the denominator effectively amplifies the noise. Therefore, when atan is used to extract δφ, two Ascans are chosen to generate a Doppler phase shift (absolute value) that is sufficiently smaller than π/2 to prevent phase wrapping artifact as well as amplified noise. When atan2 is used to achieve improved performance of Doppler analysis, the desirable Doppler phase shift (absolute value) has to be sufficiently larger than π/2 to prevent amplified noise and has to be sufficiently smaller than π to prevent phase wrapping artifact. Therefore, with Doppler phase shift calculated using four-quadrant arctangent function, the choice of time interval between Ascans depends on the noise of the amplitude OCT signal and is a non-trivial task.

Disclosures
The authors declare that there are no conflicts of interest related to this article.