1.

Introduction

Elasticity imaging in optical coherence tomography (OCT) has attracted much attention over the last 15 years following the seminal publication by Schmitt,¹ in which speckle-tracking by correlation processing was proposed to reconstruct displacements and then strains in quasistatically deformed tissues. Since then, a number of dynamic methods have also been studied, requiring auxiliary excitation of vibrations or shear-wave fields and evaluation of their characteristics by OCT (e.g., see Refs. 23.–4).

By analogy with medical ultrasound, where elastographic imaging of tissues deformed by the ultrasound probe itself is implemented in several platforms, many elastography-related works in OCT have also been used its probe to produce quasistatic tissue deformation.^5–8 Following the initial idea of Schmitt,¹ many of these studies are based on the initial reconstruction of the displacement field in the inspected region of the tissue, where approximately uniform uniaxial stress is created by the OCT probe acting as a piston. Then, the reconstructed displacement field is differentiated to evaluate local strains that are inversely proportional to the local tissue stiffness. Independent of a particular technique used for the displacement reconstruction [including both phase-resolved methods and approaches based on speckle tracking using various digital image correlation (DIC) techniques], all such methods can be considered as displacement-based.

For phase-resolved OCT methods of determining displacements and strains, discussion of their ultimate possibilities supported by numerical simulations and physical experiments has been recently reported.⁹ The influence of image decorrelation induced by displacements and strains themselves, as well as errors introduced by additive noises, were discussed. Here, we present similar analysis of displacement and strain estimations based on correlational speckle tracking in OCT. We consider simplified analytical arguments and then present rigorous numerical results based on a recently developed and validated model¹⁰ of speckle-pattern formation in OCT and its evolution, accounting for strain-induced speckle “blinking” and splitting/desplitting due to the motion of scatterers. Such “decorrelation noise”¹¹ is much more pronounced in OCT than in conventionally discussed processing of photographic images of deformed textures in mechanical engineering problems of tracking displacements and strains.¹² This significantly complicates direct adaptation of DIC techniques to OCT, as evidenced the fact that speckle-tracking approaches have not lead to creation of practical and robust elastographic methods in OCT despite over 15 years of significant efforts. The reported results reveal ultimate accuracy and intrinsic requirements (allowable strain ranges, required spectral bandwidths, and so on) for implementing correlation-tracking methods in OCT. In particular, we revisit the widely accepted notion that for sufficiently small strains which only weakly perturb the OCT speckle pattern, speckle tracking should be always feasible.

2.

Speckle Tracking Limitations Due to Speckle Blinking–Analytical Arguments

In the OCT methods developed for reconstructing the displacement field using cross-correlation processing,^1,5,6,7,8 the information is extracted from the variations in the intensity of each image pixel produced by the displacements of scatterers in the strained tissue. In contrast to deformed photographic images processed in mechanical engineering applications,¹² coherently formed OCT scans exhibit an additional effect—the relative motion of subresolution scatterers causes speckle-intensity variations known as speckle “blinking” and “boiling.” Similar effects have been known for a long time in various problems related to coherently formed images (see e.g., Refs. 13 and 14). In the context of the discussed correlational speckle tracking in OCT images, such effects are especially important since they can significantly mask the “useful” intensity variations in OCT scans produced by collective displacements of scatterers. To enable sufficiently exact reconstruction of these displacements via correlational speckle tracking, the displacement-produced variations in the intensity pattern of OCT images should significantly exceed the masking effect of strain-induced speckle blinking. The relative contributions of the two effects can be estimated using the following simple analytical arguments.

For clarity, consider the axial displacements caused by the axial strain $s$ (assumed approximately uniform in the studied region). We place the origin of the coordinate system at the interface between the tissue and the surface of OCT probe that produces the tissue strain, and assume that the $z$-axis is directed into the tissue bulk. Then at depth $z$, the axial displacement relative to the OCT probe surface is $d\approx zs$, so that for widely used common-path OCT schemes, the upper boundary of OCT scans is motionless. We also assume that the OCT signal intensity distribution as a function of depth is $I(z)$. Then, the variations $\delta I_{\mathrm{geom}}$ in the intensity for a given depth $z$ (where the tissue constituents experience a displacement of $d=zs$) can be estimated as

In Eq. (1), we use a conventional finite difference approximation of the derivative of a function exhibiting spatial variations around some mean values. The characteristic amplitude of intensity variations between minimal and maximal values is $I_0$ and the characteristic distance between the neighboring maxima is the characteristic speckle size $D$. Since the intensity varies between minimal and maximal values by a characteristic value $\sim I_0$ over the spatial scale $\sim D/2$, the derivative can be estimated as $\sim I_0/(D/2)$.

We also note here that compressional elastography (especially for ex vivo tissue samples) is realized by placing the sample onto a motionless substrate, with the sample surface compressed by a transparent rigid plate, such that maximum displacement now occurs at the surface. In such a case, it makes sense to place the origin of the coordinate system at the motionless bottom. Correspondingly, instead of the maximum displacement in the bottom part of OCT scans [as assumed in Eq. (1)], larger displacements should occur in the upper part of the scan. Mathematically, the form of Eq. (1) would then remain largely the same (but the distance would be referenced from the motionless bottom), so we do not repeat the mathematics for that case here (bearing in mind that the only difference is whether larger displacements occur near the bottom or the top of the scan).

Next, a similar estimate as in Eq. (1) is necessary for the masking effect due to speckle blinking caused by the relative motion of subresolution scatterers inside the sample volume, the axial extent of which is determined by the coherence length $L_c$ of the probing optical field. Following Ref. 15, we argue that for a pair of such subresolution scatterers separated in the axial direction by $\mathrm{\Delta }z\le L_c$, the difference in the relative phase $\mathrm{\Delta }\phi$ of the scattered signals equals $\mathrm{\Delta }\phi =2k_0\mathrm{\Delta }z$, where $k_0$ is the central optical wavenumber. The axial strain $s$ then produces the variation in the relative phase equal to $\delta [\mathrm{\Delta }\phi (s)]=2k_0\mathrm{\Delta }zs$. For a pair of scatterers with identical scattering strengths (corresponding to intensities $I_1=I_2=I_0$), the phase-sensitive contribution to the resulting speckle intensity is determined by the interference term, $2\sqrt{I_1{I}_2}\cos [\mathrm{\Delta }\phi (s)]$. Consequently, for the initial phasing ensuring the maximum sensitivity of the interference term to small variations $\delta [\mathrm{\Delta }\phi (s)]$, the corresponding variation in the speckle intensity (strain-induced blinking) is given by

In Eq. (2), the average distance $\mathrm{\Delta }z$ between the scatterers within the sampling volume is approximated by $\mathrm{\Delta }{z}_{\mathrm{aver}}\sim L_c/2$. Note that formally the estimate given by Eq. (2) is obtained for initial mutual phasing $2k\mathrm{\Delta }z=\pi /2+\pi m$ of the scattered fields, corresponding to maximum sensitivity with respect to variations in $\mathrm{\Delta }z$. However, the slope of a sinusoidal function decreases fairly slowly towards its extremum points, such that Eq. (2) can reasonably approximate a significant portion of all speckles formed by randomly distributed scatterers.

To estimate characteristic strains that can be considered “small” or “large” for speckle blinking effects described by Eq. (2), one can argue that large strains evidently correspond to relative phase variation $\delta [\mathrm{\Delta }\phi (s)]$ in the order of $\pi$ radians and greater. The corresponding characteristic strain $s_{\text{blink}}^{*}$ producing phase variations $\delta [\mathrm{\Delta }\phi (s)]=2k_0\mathrm{\Delta }z{s}_{\mathrm{blink}}^{*}˜\pi$ and causing pronounced blinking under the reasonable assumption $\mathrm{\Delta }z\sim L_c/2$ is given by

For the ratio $L_c/{\lambda }_0\sim 10$, quite typical of OCT scanners (where ${\lambda }_0=2\pi /k$ is the optical wavelength), one obtains the estimate $s_{\mathrm{blink}}^{*}\sim {\lambda }_0/2L_c\sim 0.05$ (we recall that this ratio is determined by the OCT spectral width and central wavelength: for a Gaussian source, $L_c/{\lambda }_0\sim 0.4{\lambda }_0/\mathrm{\Delta }\lambda$;¹⁶ alternatively, the ratio $\mathrm{\Delta }k/k_0\approx \mathrm{\Delta }\lambda /{\lambda }_0$ can also be used for such estimates).

For significantly smaller strains $s<s_{\mathrm{blink}}^{*}$ (e.g., $s\sim 0.01$), the speckle blinking affects the image fairly weakly (i.e., cross-correlation between the compared images remains high), and the variations $\delta I_{\mathrm{blink}}$ due to strain-induced blinking are still proportional to the strain $s$ as given by Eq. (2). An important difference in Eq. (2) from Eq. (1) for $\delta I_{\mathrm{geom}}$ is that the variation $\delta I_{\mathrm{blink}}$ does not depend on the depth $z$ of the scatterers. Correspondingly, the ratio of these variations is indeed independent of strain [$\delta I_{\mathrm{geom}}/\delta I_{\mathrm{blink}}\ne f(s)$], but depends on the depth $z$, i.e., $\delta I_{\mathrm{geom}}/\delta I_{\mathrm{blink}}\propto z$. This means that for sufficiently large $z$, the “useful” intensity variations $\delta I_{\mathrm{geom}}$ of geometric origin can exceed the masking variations $\delta I_{\mathrm{blink}}$ due to blinking. Conversely, in more superficial tissue regions of smaller depth $z$, the blinking-induced variations $\delta I_{\mathrm{blink}}$ can significantly mask the geometric variations $\delta I_{\mathrm{geom}}$ that are used to detect displacements. We emphasize that this conclusion holds even for arbitrarily small strains, despite the fact that in this regime, speckle blinking seems insignificant in the sense that the compared images remain highly correlated. Nevertheless, for insufficiently large depth $z$, speckle blinking $\delta I_{\mathrm{blink}}$ masks $\delta I_{\mathrm{geom}}(z)$ and renders speckle tracking problematic even at arbitrarily small strains. This conclusion contradicts the common expectation that sufficiently small strains can assure the feasibility of speckle tracking in OCT images for reconstructing displacements.

Let us now use Eqs. (1) and (2) to estimate the minimal depth $z_{\min }$ beyond which “useful” geometrical variations may dominate over the blinking effect $\delta I_{\mathrm{geom}}/\delta I_{\mathrm{blink}}\ge 1$. Considering the quantity $D$ is Eq. (1), we note that the minimal speckle size (satisfying the Kotelnikov–Nyquist criterion¹⁷ for pixelated images) is two pixels. Then, Eqs. (1) and (2) yield

Assuming again $L_c/{\lambda }_0=10$, one obtains $z_{\min }\ge 40\pi \sim 100$ (px). Since the pixel size is $\sim L_c$ (5 to $10\mu \mathrm{m}$ in a standard OCT system), the scale $z_{\min }$ may reach hundreds of $\mu \mathrm{m}$, which is likely unacceptable in practice. Note that for $z_{\min }$ given by Eq. (3), we still have $\delta I_{\mathrm{geom}}/\delta I_{\mathrm{blink}}\sim 1$, such that a significant dominating influence of $\delta I_{\mathrm{geom}}$ requires even greater depths. Equation (4) is thus a reasonable estimate of the depth scale $z_{\min }$ required for OCT speckle tracking to accurately determine the displacement field.

Another important point is that for elastographic purposes, reconstruction of displacement field is only an intermediate step towards estimating local strains. To this end, one needs to ensure a sufficiently exact measurement of the difference in displacements over a certain scale $Z=z_2-z_1$. Therefore, this scale $Z$ plays a similar role as the depth $z_{\min }$ discussed previously. Under the adopted approximation that the displacement is proportional to the depth $d(z)\approx zs$, a similar proportionality holds for the difference $d_2-d_1\approx (z_2-z_1)s=Zs$, where $Z$ is the depth scale over which the estimated local strain is averaged. In fact, Eq. (4) then gives the estimate of the minimal scale (i.e., spatial resolution of the strain field) $Z_{\min }$ over which the strain can be reasonably estimated by means of speckle tracking in the presence of masking influence of speckle blinking. The above performed estimate shows that for OCT scanners with typical bandwidths $\sim 10\%$ (in terms of the ratio $\mathrm{\Delta }\lambda /{\lambda }_0$ or equivalently $\mathrm{\Delta }k/k_0$), this resolution may be degraded to unacceptable scales of hundreds of micrometers, because speckle blinking strongly limits the accuracy of displacement reconstruction. Note that these estimates may require modifications for ultrawide bandwidth OCT sources that are occasionally present in high-end research systems. For such very broadband OCT scanners with $L_c/{\lambda }_0\sim 1$ (e.g., Ref. 18), the conditions for performing the displacement-field reconstruction are more favorable because the effect of speckle blinking is strongly suppressed.

In the context of the obtained conclusions, recall that after Eq. (1) we discussed two variants of compressional elastography with maximal displacements either near the bottom, or near the sample surface. If the tissue compression is produced by the OCT probe itself, sufficiently large displacements occur at depth. There, even if the “useful” displacement-induced intensity variations exceed the masking variations related to speckle blinking, the OCT signal may have already decayed strongly, such that inevitable system noise may further degrade the efficiency of correlational speckle tracking. In such situations, an OCE scheme where maximal displacements occur at the sample surface (fixed reference arm with a separate compression plate) may yield more favorable conditions over a larger section of an OCT scan. In the latter case, the distance $z$ in Eq. (4) should be counted from the bottom.

The aforementioned simple arguments help explain why the apparently simple idea¹ of correlational speckle tracking in intensity-only OCT scans has not yielded satisfactory results in experimental configurations utilizing “standard” optical sources (e.g., ${\lambda }_0\sim 1300\mathrm{nm}$, bandwidth $\mathrm{\Delta }\lambda \sim 100\mathrm{nm}$) with $L_c/{\lambda }_0\gg 1$. We now supplement the aforementioned arguments with rigorous numerical simulations of realistic OCT image formation, to demonstrate key differences between DIC methods in mechanical engineering problems and in their application to coherently formed OCT images.

3.

Model of Image Formation Based on Close-to-Reality Procedures Used in Spectral-Domain Optical Coherence Tomography

The simulation results that follow are based on a model of OCT image formation that is described in detail previously¹⁰ and is briefly summarized here. To generate the speckle patterns, the model reproduces the main steps of image formation in real spectral-domain OCT systems, in which the scans are constructed by superposing a set of optical spectral components (typically several hundreds) ballistically scattered from tissue. These harmonics are equidistantly spaced within the bandwidth of the light source used in the OCT scanner. The model allows one to compare resulting image features for different parameters of the scanner/tissue and to simulate the speckle structures of pixelated OCT images. Of direct relevance to elastography, one can study the evolution of the speckle structure under the influence of tissue deformation/motion for various optical properties of the strained medium; the number of optical harmonics and total spectral width of the OCT-scanner source can also be varied. Displacements of the point-like scattering particles in the model can vary in a wide desired range, from subpixel and subwavelength values to those greater than the pixel size. In the simplest approximation of a uniform amplitude and phase distribution over the optical beam cross-section, the expression for a single A-scan of the OCT image can be written as¹⁰

Here, $A(q)$ is the complex amplitude in $q$th pixel (the pixel numbers $q=1\dots N$ are counted from the tissue surface); $z_q$ is the physical axial coordinate of the center of $q$th pixel; $H$ is the total depth of the image; $A_j$ is the amplitude of scattering from $j$th localized scatterer located within the optical beam, and $z_j$ is the axial coordinate (depth) of the scatterer ($0\le z_j\le H$). The wavenumbers of individual spectral components are given by $k_n=k_0+\delta kn$, where $k_0$ is the central wavelength, $n=1\dots \pm (N-1)/2$ and the distance $\delta k$ between the spectral components ensuring the desired imaging depth $H$ is $\delta k=\pi /H$. The total spectrum width is then $\mathrm{\Delta }k=\delta k·(N-1)$ for $N$ spectral components. The factor $S_S(k_n)$ in Eq. (5) describes the spectral shape of the source (often Gaussian in practice), such that the number of effectively contributing spectral components can be somewhat smaller than $N$ (the latter corresponding to an ideally flat rectangular “top-hat” spectrum). If necessary, exponential decay of the optical beam $\propto \exp (-\mu z)$ with the decay parameter $\mu$, additive noise representing the noises of the receiving array, and the desired transversal amplitude and phase profiles of the optical beam can all be introduced; however, these are not necessary for the present consideration.

A single A-scan as per Eq. (5) is an elementary part of the simulated image. Two-dimensional (2-D) B-scans (and similarly three-dimensional scans) can be formed by simulating adjacent A-scans with arbitrary degrees of separation/overlap (and desired transversal profile of the optical beam if necessary). The so-formed images contain complete complex-valued signals (i.e., both phase and amplitude), which allows one to clearly demonstrate the difference in processing intensity-only and complex-valued OCT data.

4.

Numerical Simulations of Displacement Field Reconstruction

The correlational speckle tracking implies that the correlation window taken from the reference OCT image is moved and cross-correlated with the similarly sized regions of the deformed/displaced tissue image in order to maximize the cross-correlation coefficient. The so-found displacement of the moving correlation window is understood to represent the displacement of the tissue constituents. For discrete (pixelated) images, the corresponding procedure performed using the correlation window $m_1\times m_2$ pixels in size can be represented as maximization of the following coefficient defined in the ($x$, $z$) plane

Here, the correlation window $S$ of $m_1\times m_2$ elements from the reference image centered at a point ($x,z$) is moved over the search area of $n\times k$ pixels in size and is correlated via Eq. (6) with similarly sized areas $F$ from the deformed image; the quantities ${\mu }^S$ and ${\mu }^F$ are the mean values over the correlated areas in the reference and deformed images, respectively. For identical images, $C_{x,z}(n,k)$ reaches its maximum value of unity for $n=k=0$. If the compared image is distorted, the coordinates $n*$ and $k*$ are found for which coefficient in Eq. (6) reaches a maximum within a certain search region around the initial position $(x,z)$ of the window in the reference image. The so-found coordinates $(n*,k*)$ correspond to the new position of the same group of scatterers located within the subset $S$ in the reference image. In the discrete form of Eq. (6), the correlation window can be moved with a minimal step of one pixel, which in many cases is too crude and better (subpixel) resolution of the correlation procedures is indispensable. Correspondingly, different methods of ensuring a subpixel resolution have been developed in DIC literature.¹⁹

One such approach is the use of an assumed shape of the auto-correlation function for the processed speckle patterns (correlation curve fitting method¹²). After performing rough maximization with an accuracy of an integer numbers of pixels, one can estimate the interpixel position of the correlation-curve maximum for an assumed correlation-function shape (e.g., Gaussian or locally parabolic). The efficiency of this procedure, however, strongly depends on the similarity between the chosen and actual correlation-function profile, so that the attained accuracy is not evident a priori and the error may be a significant fraction of a pixel. Thus, more efficient methods have been developed, in which assumptions about the correlation-function form are not required. These achieve maximum similarity between the reference and deformed images by applying various procedures that compensate for the distortions produced by displacements and local strains in the deformed images. In engineering applications of DIC methods, such iterative procedures based on the Newton–Raphson method demonstrate the possibility of obtaining essentially subpixel accuracy of displacement measurements up to $\sim {10}^{-3}$ pixel, and strain estimates on the order of ${10}^{-3}$.¹²

For correlational speckle tracking in OCT, it is often implicitly assumed that DIC methods should be as efficient as in mechanical engineering problems,^6–8 in accord with the aforementioned estimate. However, despite the use of the same term “speckles” in mechanical engineering (where mostly photographic images of natural or artificial textures are cross-correlated), in OCT speckles are formed due to coherent superposition of spectral components from a relatively narrow-band optical source and scattered by groups of subresolution scatterers. Due to this interferrometric (coherence) nature, speckles in OCT experience strain-induced blinking that consequently causes image decorrelation, which is not typical for photographic images. Another, somewhat more technical feature is that speckle sizes in photographic images are often significantly larger than one pixel, whereas OCT images are usually obtained with maximal possible axial resolution determined by the total spectral width of the optical source. This leads to significantly more pronounced pixelated structure of OCT images, for which interferrometrically formed speckles can be more localized (down to 1 to 2 pixels). Consequently, in addition to strain-induced speckle blinking, this feature of OCT images results in a much stronger influence of the image discrete structure on decorrelation, even for subpixel displacements. We point out that this displacement-related decorrelation can be caused either by pure translational shifts of scatterers (without their relative motion) or can be caused by strain-induced displacement (and in such a case accompanied by the aforementioned speckle blinking). In what follows, we show that these effects strongly complicate the use of DIC methods for speckle tracking in OCT, especially if amplitude-only (intensity) information is used.

To illustrate these statements, Fig. 1 schematically shows a pair of A-scans where a point-like scatterer is initially located in the middle of a pixel in the reference scan [panel (a)], and then is displaced by exactly 1/2 pixel [panel (b)]. In an A-scan, a single point-like (subresolution) scatterer is visualized as one or several bright pixels, the amplitudes of which in Fig. 1 are shown as rectangles of different height. Figure 1(a) shows that the scatterer in the reference scan is well localized mostly in one pixel, and after a half-a-pixel shift becomes spread over two neighboring pixels with significantly lower amplitudes (and eventually leaks even more distantly). Correlational displacement tracking with a sliding window requires that the shifted and split images should be shifted back in subpixel steps and cross-correlated with initial images to maximize the correlation coefficient. According to the theorem on spectral functions with shifted arguments,^20,21 any pixelated one-dimensional image can be shifted by an arbitrary value $\mathrm{\Delta }z$ (including fractional subpixel values), through the introduction of an additional phase shift of $\exp (i{k}_n\mathrm{\Delta }z)$ in each of its discrete Fourier component. Introducing such a compensating phase factor for spectral components of the shifted image and adjusting the value of the shift parameter $\mathrm{\Delta }z$, one can maximize the cross-correlation between the reference and the back-shifted image with any chosen accuracy. We emphasize that for this operation to work properly, however, complete complex-valued A-scans given by Eq. (5) are needed (i.e., both signal amplitude and phase, or equivalently two quadrature components). If a complete complex-valued signal is available for the “split” image of a scatterer [Fig. 1(b)], then such a phase shift applied to each complex-valued spectral component allows one to shift such a “split” profile back by 0.5 pixel and correctly recover the original A-scan, where the scatterer image is predominantly localized in one pixel [Fig. 1(c)]. Thus, perfect unity correlation can be recovered in the absence of other noises and the desired shift (half-pixel in this example) can be found with a high accuracy. Conversely, if a similar operation is applied to the intensity-only scan of the shifted scatterer, then the back-shifted split profile does not recover the original form [Fig. 1(d), instead of Fig. 1(a)] and the cross-correlation is significantly reduced resulting in much greater ambiguity of displacement tracking. The reason for this is clear from Eq. (5) representing the complex-valued A-scan as a result of wave interference, for which both phases and amplitudes of spectral components are important. Indeed, if we take the absolute values $|A(q)|$ of the pixels in an A-scan, its Fourier spectrum does not retain the same phases/amplitude of the spectral components as for the original complex-valued $A(q)$ and correspondingly does not give the same result of the back-shift procedure.

Fig. 1

Difference between the back-shifted one-dimensional (1-D) images of a scatterer based on complex-valued and intensity-only optical coherence tomography (OCT) A-scans. (a) 1-D image of a single scatterer with coordinate corresponding to the middle of a pixel; (b) the scatterer shifted by 0.5 pixel with the image split between two pixels; (c) back-shifted image with correct compensation of splitting using complex-valued signal; (d) back-shifted image with incorrect compensation of splitting using intensity-only signal [resulting in significant decorrelation with respect to the initial image (a)].

Back-shift procedures that do not accurately take into account displacement-produced speckle splitting/desplitting do not correctly restore the pattern of the back-shifted OCT images (as in methods using amplitude-only information in the OCT signal, or the aforementioned approaches based on autocorrelation-curve fitting,²² or other DIC methods based on up-scaling/interpolation¹²). Consequently, in application to OCT images, the residual displacement-induced (i.e., unrelated to speckle blinking caused by relative motion of scatters) decorrelation even in the absence of other masking noises strongly impedes correlational tracking of speckles and causes strong uncertainty in strains that should be extracted from the reconstructed displacements. It is clear that the discussed speckle splitting/desplitting effect should be the largest for initial maximal localization of individual speckles [which can be as small as a single pixel for a rectangular, top-hat-like spectral shape $S_S(k_n)$ in Eq. (5)].

Smoother and narrower $S_S(k_n)$ profiles reduce the influence of speckle splitting/desplitting, but they also lead to spreading of a point-like scatterer image over several pixels by increasing the coherence length $L_c$, thus worsening the axial resolution and leading to stronger strain-induced speckle blinking; this, additionally complicates correlational speckle tracking. Indeed, according to Eq. (3), the value $s_{\mathrm{blink}}^{*}$ decreases for larger $L_c$, and consequently less strain is needed to produce blinking, which increases the “decorrelation noise.” In contrast to decorrelation of speckles due to purely geometrical displacement-induced splitting, the strain-induced speckle blinking of interference origin is due to the relative motion of the subresolution scatterers. In contrast to the previously discussed splitting/desplitting, speckle blinking cannot be compensated in principle, because the compensating phase shifts would be different for different subresolution scatterers with unknown initial positions. Such strain-induced decorrelation is specific to the coherent nature of OCT speckles and is absent for DIC processing of photographic images. Thus, speckle blinking becomes the main reason that limits the accuracy of correlational tracking of displacements in OCT compared with mechanical engineering problems.

To clearly demonstrate the importance of speckle-splitting and possibility of its compensation for complex-valued OCT scans, Fig. 2 shows several modeled 2-D correlational maps for reference and deformed OCT images.¹⁰ By analogy with similar study⁹ of phase-sensitive displacement-tracking methods, in our simulations, the tissue containing randomly distributed scatterers is subjected to uniaxial straining with the axial strain $s$, such that the displacement as a function of depth is a known linear function $d(z)=sz$. Thus, the axial coordinate $z_j$ of $j$th scatterer in Eq. (5) is changed by the value $\mathrm{\Delta }{z}_j=s{z}_j$. For the strain $s=0.003$ corresponding to Fig. 2, the previously discussed effects of splitting/desplitting of speckles are especially strong for scatterers located near the depth $z\sim 165$ pixels, where the shift value is $\sim 0.5$ pixel. Figure 2(a–1) shows a cross-correlation map found via shifting the correlation window by an integer number of pixels. The results correspondingly show a pronounced minimum (down to $\sim 0.8$ to 0.85) around $z\sim 165$ pixels. This correlation reduction is strongly dominated by speckle splitting/desplitting rather than speckle blinking that is still insignificant for $s=0.003$. In the considered case, the displacements are subpixel over almost the entire scan, and conssequently such displacement tracking in steps of integer number of pixels is definitely unacceptable for reconstructing $d(z)$. Figure 2(a–2) shows the difference between the actual dependence $d(z)=sz$ and the result of the correlational displacement-tracking performed in one-pixel steps. Consequently, the tracking error in such a case is as large as $\pm 0.5$ pixel.

Fig. 2

(a–1),(b–1), and (c–1) Simulated cross-correlation maps and (a–2), (b–2), and (c–2) corresponding speckle-tracking error based on the reference and uniaxially strained images. The speckle-tracking error as a function of depth in the right column is given for a fixed lateral position of the correlation window in the middle of the scan (i.e., not laterally averaged over the entire scan width). Strain value is $s=0.003\ll s_{\mathrm{blink}}^{*}$; the corresponding strain-induced displacement as a function of depth is $d^{\mathrm{strain}}(z)=0.003·z$. The correlation-window size is $16\times 16$ pixels, and the spectral width ratio $\mathrm{\Delta }k/k_0\sim 1/8$ (typical of many OCT scanners). Displacements close to 0.5 px cause maximal speckle-splitting at $z\sim 165\mathrm{px}$. Panels (a-1) and (a–2) are obtained via maximizing the cross-correlation by moving the correlation window in integer-pixel steps, and demonstrate strong tracking errors (up to $\pm 0.5\mathrm{px}$) and strong decorrelation near the depth $z\sim 165\mathrm{px}$. Panels (b–1) and (b–2) are obtained for complex-valued data using spectral sub pixel back-shift procedure which correctly compensates speckle splitting over the entire scan area, thus ensuring uniform correlation map and minimal tracking errors. Panels (c–1) and (c–2) correspond to the same subpixel tracking procedure, but applied to amplitude-only data with lost phase (label A), as well as to amplitude-only data with pseudophase found via Hilbert transform (label H). Compared to full complex-valued signals (b–1 and b–2), these give significantly larger residual decorrelation and tracking errors due to incorrectly compensated speckle splitting. Note that the color bars are not identical for all three correlation maps, to emphasize the spatial inhomogeneity of the residual decorrelation; furthermore, the vertical scale in plot (a–2) is 10 times larger than in (b–2) and (c–2).

The next row in Figs. 2(b–1) and 2(b–2) demonstrates the correlation map and the displacement-tracking error found with a subpixel search and accurate compensation of speckle splitting/desplitting. This was ensured by introducing a compensating phase factor $\exp (i{k}_n\mathrm{\Delta }z)$ for each complex-valued spectral component of the image contained within the correlation window in the deformed image. Then the value of $\mathrm{\Delta }z$ was found to maximize the cross-correlation coefficient, which corresponds to finding the correlation-window displacement. Since $\mathrm{\Delta }z$ is a continuous rather than discrete quantity, this maximization can be found with a desirable subpixel accuracy (for example, using a conventional bisection method ensuring an accuracy of $1/2^n$ pixels in $n$ steps²³). The resultant correlation map of Fig. 2(b–1) demonstrates rather high and uniform values over the entire scan cross-correlation level (higher than 0.98) including the region around the depth $z\sim 165$ pixels, where the displacements are close to 0.5 pixel and the effect of speckle splitting is especially strong. However, the use of complex-valued information renders it possible to efficiently compensate speckle splitting induced by translational displacements. The difference between the tracked displacements and the exact function $d(z)=sz$ assumed in the simulations is much smaller than in Fig. 2(a–2), $\sim {10}^{-2}$ pixel as shown Fig. 2(b–2). However, this error remains nonzero because of the residual image decorrelation due to strain-induced speckle blinking even in the absence of other noises. This illustrates that in contrast to the displacement-induced decorrelation, the blinking-induced decorrelation cannot be compensated even for full complex-valued signal (since random initial positions of all scatterers cannot be known for real-life scans).

Further, Fig. 2(c) demonstrates significant worsening of the correlational tracking accuracy if the phase information contained in complex-valued signal is lost, although the used procedures of subpixel tracking are the same as for Fig. 2(b) obtained for the complete complex-valued signal. The correlation map shown in Fig. 2(c–2) is obtained for the scans in which amplitude-only information is retained (left half of the plot marked by letter A), as well as for the same intensity-only scans in which the pseudophase information is reconstructed via Hilbert transform¹⁷ [right half of the correlation map Fig. 2(c–1) marked by letter H]. The corresponding tracking errors (also marked by letters A and H, respectively) are shown in Fig. 2(c–2). Both curves demonstrate noticeably lower ($~6$ times) tracking accuracy and lower cross-correlation in the correlation maps than for the complete complex-valued signal, which is explained by incorrect compensations of speckle splitting/desplitting effects because of lost phase information in the amplitude-only signal. We emphasize that the pseudophase found via Hilbert transform is not the phase of the actual complex-valued signal, such that the back-shift procedure applied to such pseudocomplex-field data does not ensure correct compensation of the speckle-splitting effect and does not give better accuracy of the maximization of cross-correlation and the displacement tracking compared to the intensity-only data [Figs. 2(c–1) and 2(c–2)].

4.1.

Accuracy of Displacement Tracking for Complex-Valued Optical Coherence Tomography Scans Limited by Strain-Induced Speckle Blinking

In view of the aforementioned arguments that show limitations of analyzing amplitude-only scans, we now focus on complex-valued signal processing to demonstrate the influence of strain-induced speckle blinking in more detail, because speckle blinking is the main factor responsible for insufficient accuracy of DIC methods in OCT. Simulated speckle patterns are $256\times 256$ pixels in size. In the middle depth, a stiffer layer with a thickness of 50 pixels and twice the stiffness is placed. The initial pattern is subjected to uniaxial straining, such that over the majority of the simulated image, uniaxial strain is close to the average value 0.5%; within the stiffer layer, the local strain is half, i.e., 0.25%.

Figure 3 shows cross-correlation maps (left column), the reconstructed depth-dependences $d(z)$ of vertical displacements (middle), and their derivatives corresponding to estimated local strains (right column). The tracking of displacements was made by maximizing the correlation coefficient using a sliding correlation window ($16\times 16$ pixels in size) and the above-described spectral back-shift procedure applied to simulated complex-valued scans for $\mathrm{\Delta }k/k_0=1/16$, $\mathrm{\Delta }k/k_0=1/8$, and $\mathrm{\Delta }k/k_0=1/2$. The first two values are close to typical parameters of OCT scanners, for which typical central wavelengths fall in the range of 0.8 to $1.3\mu \mathrm{m}$, and typical spectral widths of $\sim 80$ to 120 nm. Ratio $\mathrm{\Delta }k/k_0=1/2$ corresponds to a super-broadband light source, for example similar to that used in Ref. 18.

Fig. 3

Simulated cross-correlation maps and correlational speckle tracking results from complex-valued scans obtained for uniaxially strained sample for three different spectral parameters of the source. The image size is $256\times 256$ pixels. The middle imbedded depth layer (50 pixels in thickness) is two times stiffer than the surrounding media. The average vertical strain is 0.5%, and consequently the displacements of the scatterers range from zero at the surface (imitating the rigid OCT probe) to $\sim 1.28$ pixels at the bottom. The number of scatterers in the presented examples is ${10}^4$ for each A-scan and correlation-window size is $16\times 16$ pixels in all cases. Speckle tracking utilizes complex-valued data using the subpixel spectral back-shift procedure as per Fig. 2. Row (a) is for the source spectral width ratio $\mathrm{\Delta }k/k_0=1/16$; (b) is for $\mathrm{\Delta }k/k_0=1/8$; and (c) is for super-broadband source with $\mathrm{\Delta }k/k=1/2$. Correspondingly, the effect of speckle blinking is much weaker for correlation plot (c–1) than for (a–1) and (b–1). (a–2), (b–2), and (c–2) show the reconstructed displacements, with the thin dashed lines showing the difference in the slopes between the twice-stiffer layer and surrounding tissue. (a–3), (b–3), and (c–3) The derivative of the displacement (local strain) found via simplest one-pixel finite differences in the depth direction for a fixed horizontal coordinate of the correlation-window center (blue lines); red lines are the averaged values of strain found for 16 adjacent cuts in the vertical direction. (c–3) For super-broadband case, the reduced strain within the stiffer layer is fairly well-visible for the averaged (red) curve. All plots correspond to zero additive noise except for the inset in (c–3) that is calculated in the presence of $-20\mathrm{dB}$ noise.

Correlation maps in the left column of Fig. 3 demonstrate that the spectral method efficiently compensates translational motion of scatterers and makes the residual decorrelation independent of the depth-dependent displacements. Since in Fig. 3 the additive noise (i.e., unrelated to the displacement- and strain induced decorrelation) is absent and the applied processing compensates the displacement-induced decorrelation, the residual reduction in the correlation coefficient is entirely determined by local strains. In view of this, the stiffer and less-deformed layer in the middle demonstrates smaller strain-induced decorrelation, which is especially clear in Fig. 3(a–1). The fluctuations of the decorrelation level for the regions with uniform stiffness are related to random distances between subresolution scatterers in different pixels, such that the degree of speckle blinking is randomly variable. In the absence of other noise sources, the inaccuracy of the reconstructed displacements is entirely due to the strain-induced decorrelation noise. In agreement with Eq. (3), Fig. 3 shows that under the same strain, the decorrelation noise is lower for more broadband sources characterized by smaller coherence length $L_c$, for which the characteristic strain $s_{\mathrm{blink}}^{*}$ of speckle blinking is larger.

The reconstructed displacements shown in the middle column of Fig. 3 agree well with the general trend of $d(z)$ dependence corresponding to uniaxial compression of the three-layer sample. However, for typical bandwidth sources of common OCT imagers, with $\mathrm{\Delta }k/k=1/16$ and $\mathrm{\Delta }k/k=1/8$, the variation in the slope of $d(z)$ within the stiffer layer is not so clearly visible as for the wider-broadband source with $\mathrm{\Delta }k/k=1/2$. On the other hand, for the most narrow-bandwidth case of $\mathrm{\Delta }k/k=1/16$ with the strongest speckle blinking, the stiffer layer most clearly manifests itself via lower decorrelation in the correlation map, which can be used for elastographic mapping in OCT.^24,25

The plots in the right column of Fig. 3 are shown in connection with the commonly-held opinion that mapping of local strains in OCT can be made on a micrometer scale, because $\sim 5$ to $10\mu \mathrm{m}$ resolution is typical of structural OCT scans. This statement may in fact be incorrect. Indeed, the local strain corresponds to the derivative of the displacement $d(z)$, so that for a pixelated dependence $d(z)$ extracted from structural OCT scans, the derivative found with the maximal resolution of one pixel corresponds to the difference between the neighboring discrete values $d(n)-d(n-1)$, where the depth is measured in pixels. The results of such straightforward differentiation of the reconstructed discrete functions $d(z)$ with the maximum one-pixel resolution are shown in panels (a–3), (b–3), and (c–3) of Fig. 3. The thin lines in the derivative plots are found along a single A-scan (note that the correlational-tracking procedure already performs averaging over the window size while estimating the displacements of the correlation-window center). Consequently, essentially without worsening the lateral resolution, the determined derivatives can be averaged for 16 adjacent A-scans corresponding to the correlation-window width. The results of such averaging are shown in the right column of Fig. 3 by thick lines. These curves are much less noisy, although the twice-lower-strain layer within the stiffer region is sufficiently clearly seen only for the most broadband source with $\mathrm{\Delta }k/k_0=1/2$. For “conventional” width OCT sources (top two rows of Fig. 3), even larger spatial averaging (several tens of adjacent A-scans or more) may be needed to reliably detect the presence of the stiffer layer using correlational speckle tracking. The minimal averaging over the correlation window size (in the considered $\text{case}=16\times 16$ pixels), which is intrinsic to the correlational-tracking procedure, may be insufficient for conventionally used spectral widths OCT sources. Furthermore, in real scans, additive noises are always present. The inset in Fig. 3(c–3) corresponding to 20-dB signal-to-noise ratio for the super-broadband source with $\mathrm{\Delta }k/k_0=1/2$, demonstrates noticeable worsening of differentiation quality even with smoothing over the correlation-window size.

These numerical examples and the analytical estimates in Sec. 2 consistently show that the resolution attainable in the initial OCT images (1 pixel in the ultimate case) cannot be directly translated into corresponding resolution of the resultant elastographic maps. Furthermore, even if displacements could be found for individual pixels without averaging over a correlation window (as can be imagined for phase methods in which phase variations and displacements can in principle be for individual pixels), the resultant maps of local strains (i.e., derivative of tracked displacements) will have significantly lower resolution. In the Sec. 4.2, we will discuss this issue of the minimal scale of differentiation.

4.2.

Spatial Resolution and Accuracy of Strain Estimation Using Differentiation of the Retrieved Displacements

The accuracy of evaluating local strains via differentiation of the initially reconstructed displacement $d(z)$ by correlational speckle tracking depends on a number of factors, including the accuracy of initial reconstructed of $d(z)$ and the specific differentiation method [e.g., finite differences or spectral procedures, or other methods based on fitting local slopes of measured dependences $d(z)$²⁶]. In addition, the discussed type of the processed signal (complex-valued or intensity-only OCT scans) and the spectral width of the optical source, the correlation-window size and certain additive noises also determine the accuracy of the strain evaluation.

To better understand why the examples considered previously yielded rather inaccurate derivatives, although the additive noises were absent and did not corrupt the reconstructed displacements $d(z)$, recall that the derivative $d^{\prime }(z)$ can be approximated by the finite difference

Eq. (7)

$$d^{\prime }(z)\approx \frac{\mathrm{\Delta }d}{\mathrm{\Delta }z}=\frac{d_2-d_1}{z_2-z_1}.$$

The accuracy of such estimate is determined by both numerator and denominator. The pixelated nature of the image implies that the initial positions of the scatterers visualized within a given pixel can be known to within an uncertainty of $\pm 0.5$ pixel [attainable for unrealistic rectangular top-hat spectral shape $S_S(k_n)$]. For realistic spectral shapes (e.g., approximated by a Hanning or Gaussian function), this uncertainty becomes larger ($\pm 1$ pixel) and larger still if the spectrum is strongly apodized. This implies that for neighboring pixels, the denominator $\mathrm{\Delta }z=z_2-z_1$ gives 100% to 200% uncertainly of the derivative estimate. This uncertainty persists even in the hypothetical case error-free determination of the numerator $\mathrm{\Delta }d=d_2-d_1$. The role of uncertainties in $\mathrm{\Delta }z=z_2-z_1$ and $\mathrm{\Delta }d=d_2-d_1$ is illustrated schematically in Fig. 4.

Fig. 4

Roles of uncertainties in the initial position of scatterers $\mathrm{\Delta }z=z_2-z_1$ and tracked displacements $\mathrm{\Delta }d=d_2-d_1$ for estimation of the derivative $d^{\prime }(z)$, i.e., local strain. The slope of the actual derivative is shown by the thick solid line and the thinner dashed lines show the range of uncertainty.

The first conclusion based on the denominator properties in Eq. (7) is that the minimal scale for estimating the derivative will significantly exceed the intrinsic minimal uncertainty $2{\sigma }_z\ge 1$ pixels in the difference of initial positions $\mathrm{\Delta }z=z_2-z_1$ of subresolution scatterers. This intrinsic property of pixelated scans shows that the possibility of micrometer resolution of strain fields in OCT is not possible even for zero error of tracked displacements. In reality, the latter gives additional and often dominant contribution to the overall uncertainty. Thus, one can roughly estimate that in Eq. (7), the minimal intrinsic uncertainty $2{\sigma }_z\ge 1$ pixels determined by the finite resolution of the $z$ coordinate in discrete OCT images requires a minimal scale of

Eq. (8)

$$z_2-z_1=\mathrm{\Delta }{z}_{\mathrm{resol}}^{\min }\gg 2{\sigma }_z,$$

to ensure fairly accurate strain estimate [Eq. (7)]. For example, even for ultrawide (not apodized) rectangular spectrum (yielding minimal possible $2{\sigma }_z=1$), one needs $z_2-z_1=\mathrm{\Delta }{z}_{\mathrm{resol}}^{\min }\ge 10$ pixels (px) to ensure the derivative estimation error of $\sim 10\%$.

Another source of errors in the strain estimate is the inaccuracy $2{\sigma }_d$ of the numerator $\mathrm{\Delta }d=d_2-d_1$ related to the error ${\sigma }_d$ in the displacement measurement. Since within the estimation scale $\mathrm{\Delta }z$, the estimated strain is assumed constant, one can write that $\mathrm{\Delta }d=s\mathrm{\Delta }z$. Then by analogy with Eq. (8), one should require that the measured difference in the displacements should significantly exceed the measurement inaccuracy

Eq. (9)

$$\mathrm{\Delta }d=s\mathrm{\Delta }z\gg 2{\sigma }_d.$$

If by analogy with the previously discussed denominator, we require that the inaccuracy $2{\sigma }_d$ of the numerator $\mathrm{\Delta }d$ in Eq. (7) does not exceed 10%, i.e., $\mathrm{\Delta }d/2{\sigma }_d\equiv s\mathrm{\Delta }z/2{\sigma }_d\ge 10$. Then, this required accuracy of the derivative estimate determines another minimal scale $\mathrm{\Delta }{z}_d^{\min }$ related to the displacement-measurement error ${\sigma }_d$

Eq. (10)

$$\mathrm{\Delta }{z}_d^{\min }\ge 20{\sigma }_d(s)/s.$$

Generally speaking, the displacement-measurement error ${\sigma }_d$ can be a function of strain, ${\sigma }_d={\sigma }_d(s)$, because a larger strain produces larger displacements and stronger speckle blinking, such that the error ${\sigma }_d$ increases. It is important to note that at the initial stage of decorrelation of strained OCT images (before the onset of intense speckle blinking as shown in the subsequent examples), the displacement-measurement error is a linear-in-strain function, ${\sigma }_d(s)=\chi s$, so that the ratio ${\sigma }_d(s)/s=\chi \ne f(s)$ in Eq. (10) does not depend on strain. The proportionality coefficient $\chi$, however, depends on the spectral width of the OCT source, correlation-window size, and nature of the processed signal (full complex-valued or intensity-only). Numerical simulations readily allow one to estimate the proportionality coefficient $\chi$ and investigate its dependence on these parameters. In terms of $\chi$, Eq. (10) becomes

Eq. (11)

$$\mathrm{\Delta }{z}_d^{\min }\ge 20·\chi (\mathrm{px}).$$

This, together with Eq. (8) determines the minimal scales that ensure $\le 10\%$ errors in the denominator and numerator of strain estimation [Eq. (7)]. Both scales significantly exceed the ultimate (i.e., one pixel) resolution of conventional OCT structural images. Correspondingly, the attainable resolution in elastographic maps obtained via correlational speckle tracking is limited by the maximum of the values $\mathrm{\Delta }{z}_d^{\min }$ [Eq. (11)], $\mathrm{\Delta }{z}_{\min }^{\mathrm{res}}$ [Eq. (8)], and the correlation-window size $L_w$.

To demonstrate the importance of the displacement measurement error, Fig. 5 shows the function ${\sigma }_d={\sigma }_d(s)$ simulated with the correlation-window of $16\times 16$ pixels and for two spectral widths $\mathrm{\Delta }k/k_0=1/16$ and $\mathrm{\Delta }k/k_0=1/4$ (for both amplitude-only and complex-valued scans). Additive noise is not included, such that these measurement errors are solely due to decorrelation effects. It is seen that for sufficiently small strains, the standard deviations grow linearly with increasing strain, ${\sigma }_d(s)\approx \chi s$. The slope $\chi$ for intensity-only scans is initially dominated by speckle splitting which is independent on the spectral parameters and yields ${\sigma }_d(s)$ that is significantly greater than for complex-valued scans (as argued previously). The initial independence of ${\sigma }_d(s)$ on $\mathrm{\Delta }k/k_0$ for intensity-only scans is intuitively clear, because the speckle splitting effect has purely geometrical origin due to collective motion of subresolution scatters in contrast to speckle blinking, which is due to their relative motion that changes the relative phasing of scattered waves and thus depends on the optical wavelength and bandwidth. For the used correlation-window size and the intensity-only data, ${\chi }_{\mathrm{only}}^{\mathrm{int}}\approx 7$ as seen from Fig. 5. According to Eq. (11), the corresponding minimal differentiation scale $\mathrm{\Delta }{z}_d^{\min }\ge 20$ $\chi =140\mathrm{px}$ which is usually unacceptably large, but can be significantly reduced using complex-valued signal as shown below. Figure 5 also shows that for larger strains (roughly $s>0.5\%$ for $\mathrm{\Delta }k/k_0=1/16$ and $s>2\%$ for $\mathrm{\Delta }k/k_0=1/4$), both intensity-only and full complex-valued signals give approximately the same measurement error ${\sigma }_d(s)$. Indeed, for such larger strains, the dominant contribution to the tracking errors is due to speckle blinking which cannot be compensated either for intensity-only or for complex-valued signal.

Fig. 5

Standard deviations of the reconstructed displacements ${\sigma }_d(s)$ found via correlational speckle tracking for intensity-only (curves 1) and full complex-valued OCT signals (curves 2) as a function of strain. Spectral width ratios of the source are: (a) $\mathrm{\Delta }k/k_0=1/16$ and (b) $\mathrm{\Delta }k/k_0=1/4$. Additive noise is absent. The correlation-window size is $16\times 16$ pixels, as per Figs. 2 and 3. For the used correlation-window size, the intensity-only curves (1) at sufficiently small strains have initial slope ${\chi }_{\mathrm{int}}\approx 7$ (having both a subscript and a superscript on the same symbol is confusing and unnecessary. I suggest we just use ${\chi }_{\mathrm{int}}$ and ${\chi }_{\mathrm{comp}}$, simpler and less confusing—here, in the text, in the Table caption, and in the caption of Fig. 6) independent of the source spectral properties, whereas for complex-valued data, the slopes are spectrum-dependent and significantly smaller. For larger-strain range, where speckle blinking becomes pronounced and dominate decorrelation, both complex-valued and intensity-only data give approximately the same tracking error. The thin dotted lines show initial slopes of the linear approximation of ${\sigma }_d(s)$ before the onset of pronounced blinking.

In contrast, Fig. 5 shows that at smaller strains, when speckle blinking is yet insignificant, the dominating effect of speckle splitting can be correctly compensated using the complex-valued scans, such that the residual errors are exclusively due to unavoidable speckle blinking. The latter essentially depends on the source spectral width and can be significantly reduced in very broadband systems. To demonstrate this in more detail, Fig. 6 shows the corresponding ${\sigma }_d(s)$ found for spectral width ratios $\mathrm{\Delta }k/k_0=1/16$, $\mathrm{\Delta }k/k_0=1/8$, and $\mathrm{\Delta }k/k_0=1/4$ using complex-valued scans with compensation of speckle splitting. For complex-valued signals at the initial stage of strain-induced speckle blinking, the proportionality coefficient scales as ${\chi }_{\mathrm{val}}^{\mathrm{comp}}\propto k_0/\mathrm{\Delta }k$ (see Fig. 6).

Fig. 6

Standard deviations ${\sigma }_d(s)$ of the displacements determined via correlational speckle tracking using complex-valued OCT scans for several spectral width ratios of the source, $\mathrm{\Delta }k/k_0=1/16$ (upper curve, ${\chi }_{\mathrm{comp}}\approx 4$), $\mathrm{\Delta }k/k_0=1/8$ (middle curve, ${\chi }_{\mathrm{comp}}\approx 2$) and $\mathrm{\Delta }k/k_0=1/4$ (lower curve, ${\chi }_{\mathrm{comp}}\approx 1$). The correlation-window size is once again $16\times 16$ pixels. The thin dotted lines show initial slopes that scale inversely with the bandwidth parameter $\mathrm{\Delta }k/k_0$.

In order to find the scaling of the parameter $\chi$ as a function of the correlation-window size $L_w$, it has been verified in simulations that for square windows $n\times n$ in size, the coefficient $\chi \propto 1/n$. Consequently, if one tries to decrease the correlation-window size in order to achieve better resolution (weaker smoothing) in the reconstructed displacements and strain maps, the actual result may be the opposite. Indeed, for a smaller window, greater measurement errors can result in an increase in the minimal scale required for determining strains, leading to even worse resolution of the elastographic map. Evidently, in the optimal case, the correlation-window size $L_w$ and the minimal scale of displacement differentiation given by Eqs. (8)–(11) should be the same. In this case, the window size corresponds to the best attainable elastographic resolution.

These derived scaling laws for the displacement-tracking error ${\sigma }_d$ render it possible to estimate the ultimate possibilities of strain reconstruction using correlational speckle tracking. Table 1 shows expected values of the minimal required differentiation scales depending on the signal type (complex-valued or intensity-only), spectral width $\mathrm{\Delta }k/k_0$ of the OCT system, and correlation-window size $L_w$. It is seen that the intensity-only data in general require rather large of minimal differentiation scales (on order of hundreds of pixels, which is unacceptable in practice) to extract strains from the initially reconstructed displacements. For complex-valued scans, the situation is better, although in contrast to correlational mapping of blood microcirculation,²⁷ small correlation-windows on order of a few pixels cannot be used for extracting elastographic maps (local strains) via speckle tracking. Table 1 shows that for complex-valued signals and sufficiently broadband sources, the optimal size of the correlation window can be reduced down to $~20$ to 30 pixels.

Table 1

Expected values of the minimal required differentiation scales Δzdmin for evaluating strain from the tracked displacements depending on the signal type (complex-valued or intensity-only), for different correlation-window size Lw and spectral width ratio Δk/k0 of the OCT system. The sizes are expressed in pixels, calculated for the displacement-tracking error σd corresponding to 10% accuracy in evaluating the numerator in Eq. (7) for strain estimation. Using the derived scaling laws χcomp∝k0/Δk, χint≠f(k0/Δk), and χcomp∝1/Lw (for a square window), the corresponding dependencies on other parameters can be readily obtained.

Correlation window size Lw	4×4	4×4	8×8	16×16	32×32	32×32
$\mathrm{\Delta }k/k_0$	1/16	1/8	1/16	1/16	1/16	1/8
$\mathrm{\Delta }{z}_d^{\min }$, px (intensity-only scans)	560	560	280	140	70	70
$\mathrm{\Delta }{z}_d^{\min }$, px (complex-valued scans)	320	160	160	80	40	20

4.3.

Minimal Detected Strain via Correlational Tracking of Displacements in the Presence of Noise

The previously estimated spatial resolution may or may not be acceptable in practice, but the estimates summarized in Table 1 are obtained for zero additive noises. In real OCT conditions, typical SNR values are 15 to 20 dB or less;⁹ the used model¹⁰ readily allows one to introduce various noise levels (a reasonable way is adding a random complex value with uniform phase and Gaussian amplitude distribution to each pixel). An example of noise influence on the strain-dependent displacement-tracking error ${\sigma }_d(s)$ is shown in Fig. 7 for the spectral width ratio $\mathrm{\Delta }k/k_0=1/8$ and SNR of 20 and 14 dB. The initial slope of the noise-free curve found using complex-valued scans is ${\chi }_{\mathrm{val}}^{\mathrm{comp}}\approx 2$; also in the presence of noise, ${\sigma }_d$ is seen to be initially independent of strain.

Fig. 7

Standard deviations of the displacements determined via correlational speckle tracking for complex-valued OCT scans with and without additive noise. Spectral width ratio of the source is $\mathrm{\Delta }k/k_0=1/8$. The lower curve is for zero noise, middle one is for $\mathrm{SNR}=20\mathrm{dB}$, and upper curve is for $\mathrm{SNR}=14\mathrm{dB}$. The correlation-window size is $16\times 16$ pixels. In noise-free conditions (bottom curve), the error is proportional to strain (initially linearly, shown by thin dotted line). In the presence of noise (top two curves), the tracking error is initially independent of strain and is proportional to the noise level (i.e., $2\times$ larger for 6 dB decrease in SNR, dotted horizontal lines). For larger strains, the three error curves show similar increasing behavior with strain.

Figure 7 shows that up to $s*\sim 0.005$ for 20-dB SNR and up to $s*\sim 0.01$ for 14-dB SNR, the displacement-tracking errors are approximately strain-independent and are strongly dominated by the additive noise with the standard deviations ${\sigma }_d^{\mathrm{noise}}\approx 0.01$ and ${\sigma }_d^{\mathrm{noise}}\approx 0.02$, respectively. Thus, displacements smaller than ${\sigma }_d^{\mathrm{noise}}$ are not detectable. Since to estimate strain one should estimate the displacement difference between two points [Eq. (7)], the minimal detectable displacement difference is $\sim 2{\sigma }_d^{\mathrm{noise}}$. Next, since the correlational tracking averages displacements over the scale of the correlation-window size $\mathrm{\Delta }z\sim L_{\mathrm{win}}$, the corresponding minimal detectable strain with the best possible resolution defined by the correlation-window size can be estimated as

Eq. (12)

$$s_{\min }\sim \frac{2{\sigma }_d^{\mathrm{noise}}}{\mathrm{\Delta }z}=\frac{2{\sigma }_d^{\mathrm{noise}}}{L_{\mathrm{win}}}.$$

For SNR of 20 dB, the tracking error is ${\sigma }_d^{\text{noise}}=0.01$ px as seen from Fig. 7. Then for the used $L_{\mathrm{win}}=16\mathrm{px}$, Eq. (12) yields $s_{\min }\sim 0.02/16\sim {10}^{-3}$. This strain is still negligibly detectable and will have $\sim 100\%$ uncertainty. For constant resolution $\mathrm{\Delta }z\sim L_{\mathrm{win}}$ and within the strain range $s_{\min }\sim 0.02/16\sim {10}^{-3}$ (where the absolute displacement-tracking error is dominated by the noise ${\sigma }_d^{\mathrm{noise}}\approx \mathrm{const}$.), in the strain range $s_{\min }<s\le s*$, the uncertainty should be proportionally smaller than at the very detection threshold $s\approx s_{\min }$. Further, it is clear from Fig. 7 that for strains $s>s*$ (e.g., for $\mathrm{SNR}=20\mathrm{dB}$, $s>s*=0.005$), the displacement tracking error becomes dominated by the decorrelation noise rather than additive noise, such that in a certain strain range (for $\mathrm{SNR}=20\mathrm{dB}$, this is the range $0.005<s<0.015$), the displacement measurement error scales linearly with increasing strain, ${\sigma }_d\propto s$. Consequently, for a fixed chosen spatial scale [for example, defined by Eq. (10) corresponding to 10% uncertainty of strain estimation], the relative error in strain estimation remains constant in this strain range.

The linear dependence of the uncertainty ${\sigma }_d\propto s$ has an upper bound determined by the characteristic strain $s_{\mathrm{blink}}^{*}$; in other words, limited by the onset of intense speckle blinking [Eq. (3)]. Then for even larger strains $s>s_{\mathrm{blink}}^{*}$, the displacement-tracking error grows supralinearly (faster than ${\sigma }_d\propto s$), such that the relative error of strain estimation rapidly increases. Depending on the noise level and spectral width of the optical source, the range $s_{\min }<s<s_{\mathrm{blink}}^{*}$ of measurable strains via correlational speckle tacking can be rather narrow (in Fig. 7 plotted for $\mathrm{\Delta }k/k_0=1/8$, the transition to blinking occurs at $s=s_{\mathrm{blink}}^{*}\sim 0.015$), and reasonably accurate strain estimation may require rather large scales (see Table 1), especially for “standard” OCT sources. Such strong dependence on the speckle-decorrelation effects is not typical of DIC applications in mechanical engineering. This fact, coupled with the vulnerability to noises (which in OCT are often higher than in photographic images used in other DIC applications) account for the difficulties of correlational tracking of displacements for extracting local strains in elastographic OCT.

5.

Conclusions

The following main conclusions are summarized here for convenience. (1) The accuracy of OCT correlational speckle tracking depends on the nature of the signal to be analyzed: complex-valued images give much better accuracy compared to processing of intensity-only data. (2) For small strains, where speckle blinking is insignificant, tracking inaccuracy is dominated by the influence of speckle splitting-desplitting determined by the discrete character of the OCT images. Unlike speckle blinking, this effect is independent of the spectral bandwidth ratio $\mathrm{\Delta }k/k_0\approx \mathrm{\Delta }\lambda /{\lambda }_0$ (see Fig. 5) and can be compensated only using a full complex-valued OCT signal. (3) For typical spectral bandwidth OCT imagers, speckle blinking becomes pronounced for characteristic strains $s_{\mathrm{blink}}^{*}\sim {10}^{-2}$. In contrast to the displacement-induced speckle-splitting, the blinking effect cannot be compensated even using complex-valued data and becomes the main source of speckle-tracking inaccuracy even in the ideal case of no other noise sources. This effect is absent for classical DIC used in mechanical engineering applications. The use of ultrabroadband OCT scanners can reduce the influence of speckle blinking and render speckle tracking more accurate. (4) Quantitative predictions for ultimate resolution in elastographic OCT have been made, and indicate that using super-broadband OCT imagers and full complex-valued signal in the absence of other noises, strain can be estimated with 10% to 20% accuracy using correlation windows 20 to 30 pixels in size. Minimal detectable strain is determined by additive noises (unrelated to decorrelation noise) as discussed in Sec. 4.3. Therefore, the minimal detectable strain can strongly vary over the OCT scan area depending on the local SNR values (that usually decreases with increasing depth because of signal decay).

In elastographic mapping, the initially reconstructed displacements are differentiated for estimating local strains. The resultant accuracy of strain estimation and minimal detectable strain both strongly depend on the evaluation spatial scale. The latter scale limits the attainable resolution of strain mapping, and is considerably lower than the resolution of initial structural OCT images; even in the optimal case is limited by the correlation-window size that can hardly be reduced below 20 to 30 pixels (Table 1). Although the differentiation-related problems are common for correlational and phase methods of displacement tracking, the latter exhibits better tolerance to strain-induced speckle blinking and additive noises and thus seem more promising in OCT elastography (see analytical arguments in Ref. 15 and rigorous numerical results in Ref. 9). However, in contrast to phase methods that are sensitive primarily to axial displacements of scatterers, correlational speckle tracking can estimate lateral displacements as well; this is very important for reconstructing the full strain tensor and could be realized using super-broadband OCT scanners.

Finally, it is worthwhile to mention that in view of the error-prone differentiation procedure, it may be advantageous to exclude this intermediate stage of global reconstruction of the displacement field and then its differentiation. For example, direct measurement of shear- or Rayleigh wave velocities can be performed.²⁸ A possible alternative approach is considered in Refs. 24 and 29, and other options will be discussed elsewhere.