Quantitative Mapping of Strains and Young Modulus Based on Phase-Sensitive OCT

Abstract

In this chapter we consider mapping of local strains and tissue elasticity in optical coherence tomography (OCT) based on analysis of phase-sensitive OCT scans. Conventional structural OCT scans correspond to spatially resolved mapping of the backscattering intensity of the probing optical beam. Deeper analysis of such sequentially acquired multiple OCT scans can be used to extract additional information about motion of scatterers in the examined region. Such detailed analysis of OCT scans has already resulted in creation of OCT-based visualization of blood microcirculation, which has been implemented in several commercially available devices, especially for ophthalmic applications. Another functional extension of OCT emerging in recent years is the OCT-based elastography, i.e., mapping of local strains and elastic properties in the imaged region. Here, we describe the main principles of local strain mapping in phase-sensitive OCT with a special focus on the recently proposed efficient vector method of estimation of interframe phase-variation gradients. The initially performed mapping of local strains is then used for realization of quantitative compressional elastography, i.e., mapping of the Young modulus and obtaining stress-strain dependences for the studied samples. The discussed principles are illustrated by simulated and experimental examples of elastographic OCT-based visualization. The presented elastographic principles are rather general and can be used in a wide area of biomedical and technical applications.

1. Introduction

Scans in optical coherence tomography (OCT) strongly resemble those obtained by ultrasound scanners. Conventional OCT images correspond to spatially resolved visualization of the backscattering intensity for optical waves similarly to visualization of backscattering intensity of ultrasonic waves. Deeper analysis of such images, especially applied to sequentially acquired multiple scans, opens possibilities to extract a rich additional information about motion of scatterers in the examined region. The development of such functional extensions in OCT imaging was in many aspects stimulated by the analogous trends in ultrasound.

In this context, probably the most well-known extension of ultrasound-based imaging is visualization of flows (first of all, imaging of blood flows in medical ultrasound angiography [1]). In OCT, generically similar principles in the 2000s were also used to realize OCT-based visualization of blood microcirculation. Since the typical imaging depth in OCT is ∼1–2 mm and lateral field of view ∼several millimeters, OCT-based angiography visualizes the microcirculation on a smaller scale in comparison with ultrasound but correspondingly with a higher resolution typical of OCT, where the typical resolution is 5–15 μm. Now the angiographic modality in OCT is implemented in several commercially available devices, especially for ophthalmology [2] but also for other medical applications, e.g., in oncology [3, 4, 5].

Another functional extension of OCT, the development of which was also inspired by analogous trends in ultrasonic imaging, is the optical coherence elastography (OCE), i.e., mapping of deformations (including local strains) and elastic properties (first of all, the Young modulus) in the imaged region. This direction in OCT development was triggered in 1998 by the seminal publication by Schmitt [6]. In that paper, Schmitt considered the possibility of transferring to OCT the ideas proposed in medical ultrasound in the very beginning of 1990s [7]. This approach was rather successfully developed in subsequent years [8] and since ∼2000 was realized in several commercially available ultrasound platforms.

In OCT, however, the development of elastographic mode passed with 10–15 years delay. Despite the evident similarity between ultrasound and OCT scans, the transferring to OCT, the elastographic principles successfully realized in ultrasound appeared to be rather challenging. Sufficiently successful realizations of OCT-based elastography were demonstrated only during the last ∼5 years. Similarly to OCT-based microangiography, elastography in OCT is mostly focused on biomedical applications (e.g., [9, 10, 11, 12]), but certainly similar principles attract interest for engineering applications, e.g., for testing polymers [13].

In what follows, we briefly overview the main trends in the development of OCE, including measurement of strains and quantitative OCT-based mapping of the Young modulus with focus of the so-called compressional phase-sensitive OCE that became one of the most active directions in the development of OCE.

2. Basic principles used in OCE

The basic principles used for realization of OCE generically are rather similar to those used in ultrasound imaging for elastographic purposes. The primary goal of elastography is to estimate the shear modulus of the studied material (at least in the relative sense without absolute quantification). In biomedical applications this interest is explained by the fact that all soft biological tissues have very similar value of the bulk modulus that varies very insignificantly in various states of the tissue. In contrast, the shear modulus for the same state of the tissue may exhibit much stronger variability, up to several times and even orders of magnitude. Therefore, it is the observation of shear modulus variability which is especially interesting for studying structural changes in soft materials. The materials that are called “soft” at the intuitive level formally belong to the class of “nearly incompressible” or “water-like” materials, for which their Poisson’s ratio $\nu$ is very close to the upper physically allowable limit $\nu \to 0.5$, i.e., to the value typical of liquids. The condition $\nu \to 0.5$ also means that the shear modulus $G$ of such a material is much smaller than its bulk modulus $K$: $G/K\ll 1$. When a rod, made of such a material with free boundaries in lateral direction, is subjected to axial loading that causes its axial strain ${\epsilon }_{\mathit{zz}}$, the soft material experiences lateral strains ${\epsilon }_{\mathit{xx}}={\epsilon }_{\mathit{yy}}\approx -0.5{\epsilon }_{\mathit{zz}}$. Consequently, the expansion in one direction is compensated by contraction in the orthogonal directions, such that the volume of the so-loaded sample remains almost invariable. This explains why materials with $\nu \to 0.5$ are called “nearly incompressible.”

Another remarkable feature of such materials with $\nu \to 0.5$ is a rather specific relationship between their Young modulus $G$ and shear modulus $E=2\left(1+\nu \right)G\approx 3G$. Therefore, instead of direct measurements of modulus $G$ (via observation of propagation of shear or surface waves, which is used in wave variants of both ultrasonic elastography [14] and OCE [15]), it is possible to use reaction of the tissue under longitudinal uniaxial stress to evaluate the Young modulus $E=3G$.

The idea to use quasistatic uniaxial stress for estimation of the Young modulus was proposed for ultrasound [7] and was transferred to OCT in paper [6]. In practice, the strain, which can be characterized as fairly close to uniform and uniaxial, is created in the vicinity of a piston pressed onto a tissue that is not stuck at the piston-tissue interface and can fairly freely slide laterally. This strain is created by compressional loading, so that the idea of measuring the Young modulus in such a configuration is called compressional (or compression) elastography.

The key point in realization of compressional OCE is, therefore, estimation of axial strains in the so-compressed material by analyzing a series of OCT scans acquired during the material compression. Comparison of such scans can be used to reconstruct axial displacements $U\left(z\right)$ of scatterers, and the local axial strains can be estimated by evaluating spatial gradients of the displacement field:

$${\epsilon }_{\mathit{zz}}=\frac{\partial U}{\partial z}$$E1

For tracking the displacements of scatterers, conventionally correlation principles have been discussed and fairly successfully realized in various applications, including medical ultrasound elastography [8] and engineering problems, where this processing is applied to sequences of photographic images of a deformed surface [16]. In paper [6] similar principles were supposed to be transferred to the analysis of OCT images. However, the attempts to directly transfer the correlational principles of estimating strains by tracking displacements using consequently acquired OCT scans were not very successful [17] in the sense that the correlationally reconstructed displacement fields were rather noisy. In view of this, their numerical differentiation required for reconstructing local strains did not give satisfactory results.

The reason of this was that unlike photographic images, OCT scans are characterized by a peculiar speckle structure originated in OCT scans from the interference of optical waves scattered from sub-resolution scatterers. This speckle pattern is rather sensitive to deformation of the imaged material, because straining produces relative displacements of the sub-resolution scatterers (for which a quarter-wavelength mutual displacement in the axial direction changes the character of interference of the scattered waves from constructive to destructive and vice versa). In view of this, for fairly moderate strains ∼ a few percent and even less, speckles in OCT images may demonstrate pronounced “boiling” and “blinking” resulting in strong decorrelation of the compared OCT scans. The intuitively attractive idea to use the correlational speckle tracking for “sufficiently small strains,” for which the abovementioned decorrelation could be avoided, did not help too much. The reason is that in fact both masking distortions producing the decorrelation and the variations in the speckle structure that are used for the tracking purposes do appear simultaneously and are proportional to the same order of strain [18]. Thus, even for small strains, the decorrelation of OCT images usually does not allow for sufficiently precise correlational speckle tracking allowing for performing subsequent differentiation and finding strains via Eq. (1). The possibilities of correlational speckle tracking in OCT imaging of deformed tissues were analyzed in detail [19]. The conclusion was that by the abovementioned reasons, the correlational approach to strain reconstruction in OCT can be practically feasible for super-broadband spectrum of the illuminating source (like used in [20]), for which the deformation-induced speckle blinking/boiling is strongly suppressed. Otherwise the correlation approach is operable for only approximate strain estimations, i.e., with averaging over large portions of the entire OCT scan. By this reason in what follows, we will focus on another realization of compressional OCE, based on phase-resolved measurement. This approach has proven to be rather promising and is especially actively developed in recent years.

In OCT the phase of the backscattered signal is naturally available and can be readily used to track the displacements of scatterers using the well-known relationship between the variation $\mathrm{\Phi }$ in the backscattered-wave phase and the axial scatterer displacement $U$:

$$U=\frac{{\lambda }_0\mathrm{\Phi }}{4\mathit{\pi n}},$$E2

where ${\lambda }_0$ is the optical wavelength in vacuum and $n$ is refractive index of the material. It can be shown that the phase of the OCT signal can be more tolerant to strain-induced decorrelation [21]. Consequently, even for “typical” OCT systems (i.e., without the need to ensure a super-broadband spectrum for reduction of deformation-induced decorrelation), phase measurements related to the displacements of scatterers can be made much more reliably than the correlational speckle tracking. This is an important advantage of phase-sensitive approaches to estimation of strains.

Thus, estimation of the axial gradient of phase variations $\mathrm{\Phi }\left(z\right)$ makes it possible to estimate local strains via Eq. (1). It is important to point out that for pixelated OCT images, both the displacements and distances are naturally measured in pixels (the physical values of both quantities being dependent on the refractive index). Therefore, the gradients of the phase variations, also calculated in pixels, give correct values of strain without the necessity to know the refractive index, which is a positive feature of phase-sensitive OCT-based strain measurements.

Another important point is that, for unambiguous relation between the observed interframe phase variation $\mathrm{\Phi }\left(z\right)$ and displacements $U\left(z\right)$ of scatterers in the observation point, the condition $U\left(z\right)<\lambda /4$ should be fulfilled. For larger displacements the phase wrapping occurs because of periodicity of the dependence of parameters of a wave on its phase. In view of this, the interframe displacement can be directly extracted from the interframe phase variation only with an uncertainty to the unknown integer number of wave periods. To exclude this ambiguity in estimations of the displacements, the conventional approach to realization of unambiguous estimation is to ensure sufficiently small interframe displacements that do not exceed $\pm \lambda /4$ [22].

Alternatively, in the cases when the displacement of scatterers is caused by the material straining, the displacements may gradually increase over the OCT scan in a wide range from essentially sub-wavelength (without phase wrapping) to super-wavelength values (with phase wrapping). Although multiple phase wrappings may occur over the entire imaged depth, in order to reconstruct a continuous function $U\left(z\right)$ even for super-wavelength displacements, the gradual increase in the phase variation makes it possible to apply conventional phase unwrapping procedures by adding $2\pi$ rad. at every depth where the phase variation exhibits a $2\pi$ jump. However, the unwrapping procedure is error-prone, so that because of measurement noises always present in real measurements (especially in areas with weaker scattering and reduced signal-to-noise ratio (SNR)), the unwrapped phase may be rather noisy. Consequently, numerical differentiation of such a noisy function $\mathrm{\Phi }\left(z\right)$ (or equivalently $U\left(z\right)$) can be made with a reasonable accuracy only with averaging over large portions of the entire imaged depth.

In what follows we consider an efficient recently developed approach to strain-reconstruction phase-sensitive OCE based on local estimation of the gradient $\mathit{dU}\left(z\right)/\mathit{dz}$ without preliminary reconstruction of the displacement function in the depth ranges where the estimates of the total displacement $U\left(z\right)$ may be ambiguous because of phase wrapping. Then some examples of the application of the phase-sensitive strain visualization for mapping both fairly rapidly varying strains and very slow varying ones (with special optimizations of the measurement procedures) will be considered. Finally, we will discuss a realization of quasistatic compressional OCE with application of reference translucent layers with pre-calibrated stiffness for obtaining quantitative estimates of the Young modulus.

3. Local estimates of strains in phase-sensitive OCE using the “vector method”

Significant progress in strain mapping has been achieved in OCE in recent years due to transition to the use of phase-resolved OCT data. Quite a detailed discussion of axial strain estimation based on phase-resolved data for compared deformed and reference OCT scans was presented in [23]. In that paper, the least-square method (including the improved version with amplitude weighting) was considered to estimate local gradients of function $\mathrm{\Phi }\left(z\right)$ using averaging over a window with a size significantly smaller than the entire image size. In [24, 25] another procedure for finding phase gradients was proposed. That approach was called “vector method,” because it operates with complex-valued OCT signals (i.e., with signals characterized by amplitude and phase). Such signals can be considered as vectors in the complex plane, which explains why the method is called “vector.” The estimated phase gradient is singled out at the very last stage of the signal processing.

Comparison with the least-square fitting of the $\mathrm{\Phi }\left(z\right)$ slope (even with amplitude weighting to suppress noisy small-amplitude pixels) for the same processing window size demonstrated superior robustness of the vector method with respect to strain-induced speckle-decorrelation noise and other measurement noises [24]. One of the advantages of the vector method is due to the fact that the amplitude weighting is also intrinsically made in the vector method (since the signal amplitude determines the absolute value of the corresponding vector) and, furthermore, especially strong phase errors exceeding $\pi /2$ rad. are naturally suppressed in this method even for strong signals. The increased tolerance of the method to decorrelation noises makes this method operable under elevated interframe strains (up to ∼${10}^{-2}$). Consequently, for comparing other noises, the possibility to operate with elevated strains with larger phase gradients corresponds to effectively higher signal-to-noise ratio (SNR). This makes the vector method especially suitable for visualization of aperiodic strains (with magnitudes in the range ${10}^{-3}-{10}^{-2}$), for which enhancement of SNR via conventional periodic averaging is not possible.

Consider the main steps of signal processing in the vector method. Let the complex-valued signal in each pixel $\left(m,j\right)$ in the reference scan be written as $a_1\left(m,j\right)=A_1\left(m,j\right)exp\left[i\cdot {\varphi }_1\left(m,j\right)\right]$ and $a_2\left(m,j\right)=A_2\left(m,j\right)exp\left[i\cdot {\varphi }_2\left(m,j\right)\right]$ correspond to the deformed scan. In each scan the signal amplitudes A_1,2 and phases ϕ_1,2 are random because of random positions and scattering strength of scatterers. However, the strain-induced variations in the phase are caused by fairly regular displacements of scatterers and can be found by singling out the phase in the product:

$$a_2\left(m,j\right)a_1^{\ast }\left(m,j\right)\equiv b\left(m,j\right)=B\left(m,j\right)exp\left[i\cdot \mathrm{\Phi }\left(m,j\right)\right]$$E3

Here, the asterisk denotes complex conjugation, $B\left(m,j\right)=A_2\left(m,j\right)A_1\left(m,j\right)$ and $\mathrm{\Phi }\left(m,j\right)\equiv {\varphi }_2\left(m,j\right)-{\varphi }_1\left(m,j\right)$. Phase variation $\mathrm{\Phi }\left(m,j\right)$ for pixel $\left(m,j\right)$ is related to the axial displacements $U\left(m,j\right)$ of the scatterers with coordinates close to pixel $\left(j,m\right)$ via Eq. (2).

Usually the axial strain $\partial U/\partial z$ is evaluated by finding axial gradient of the discreet phase variation $\mathrm{\Phi }\left(m,j\right)$ with averaging within a processing window $N_x\times N_z$ pixels in size (e.g., using the least-square method [23]). Alternatively, the averaging procedures can be performed with complex-valued quantities (3) considered as vectors in the complex plane. Then the phase gradients can be singled out only at the very final stage.

In [24] the vector approach was considered for laterally weakly inhomogeneous phase variations (i.e., for nearly uniaxial straining of the material along $z$ axis). In such a case, the complex-valued quantities $b\left(m,j\right)$ (which may exhibit phase fluctuations due to various measurement noises and strain-induced decorrelation) can be laterally averaged within a chosen processing window to obtain an array $\overline{b\left(j\right)}\equiv B\left(j\right)exp\left[i\cdot \mathrm{\Phi }\left(j\right)\right]$ with more regular phase $\mathrm{\Phi }\left(j\right)$:

$$\overline{b\left(j\right)}={\sum }_{m=1}^{N_x}b\left(m,j\right)={\sum }_{m=1}^{N_x}{A}_2\left(m,j\right)A_1\left(m,j\right)exp\left\{i\cdot \left[{\varphi }_2\left(m,j\right)-{\varphi }_1\left(m,j\right)\right]\right\}$$E4

Here $m=1..N_x$ is the horizontal index of the vertical columns in the processing window, and $j=1..N_z$ is the index of the horizontal rows (see schematic Figure 1). For the averaged vector $\overline{b(j})$, contributions of noisy small-amplitude pixels become significantly suppressed, such that they do not strongly distort the phases $\mathrm{\Phi }\left(j\right)$ of the resultant vectors $\overline{b(j})$. Furthermore, Figure 1, where Eq. (4) is represented as the summation of vectors, illustrates that the strongest phase errors give minimal distortions of the averaged phase $\mathrm{\Phi }\left(j\right)$. (See the second panel in Figure 1, where it is shown that the individual noisy vectors $b\left(m,j\right)$ with almost opposite directions with respect to the direction of the averaged vector $\overline{b\left(j\right)}$ weakly affect the orientation of $\overline{b\left(j\right)}$.) For nonzero strain, the so-constructed phase variation $\mathrm{\Phi }\left(j\right)$ depends on the discrete (pixelated) vertical coordinate $j$. In principle, similarly to [23], the vertical phase gradient then can be found by conventional least-square fitting of $\mathrm{\Phi }\left(j\right)$ within the vertical size $N_z$ of the processing window as was used in [26].

Figure 1.

Schematic of evaluation of the axial phase gradient in the vector approach for laterally nearly homogeneous strains. In the vector diagrams, the complex-valued summands $b\left(m,j\right)$ and $c\left(j\right)$ are shown as vectors in the complex plane. All intermediate transformations are performed with complex-valued signals, and the sought phase gradient is singled out at the last stage.

However, by analogy with the initial vector averaging, the next step of finding the vertical gradient of function $\mathrm{\Phi }\left(j\right)$ can also be done without explicit extraction of the phases $\mathrm{\Phi }\left(j\right)$. Namely, one can define complex-valued quantities $c\left(j\right)$ that contain vertical increments $\mathrm{\Delta }\Phi \left(j\right)=\mathrm{\Phi }\left(j+1\right)-\mathrm{\Phi }\left(j\right)$ of the horizontally averaged interframe phase variations $\mathrm{\Phi }\left(j\right)$ (see the third panel in Figure 1):

$$c\left(j\right)=\overline{b\left(j+1\right)\cdot }\overline{b\left(j\right)}{}^{\ast }\equiv B\left(j+1\right)B\left(j\right)exp\left\{i\cdot \left[\mathrm{\Phi }\left(j+1\right)-\mathrm{\Phi }\left(j\right)\right]\right\}$$E5

The quantity $c\left(j\right)$ can also be used in the normalized form $c_{\mathit{norm}}\left(j\right)=c\left(j\right)/{\left|c\left(j\right)\right|}^{\alpha }$. It was verified that good results can be obtained using normalization with $\alpha =1$

$$c_{\mathit{norm}}\left(j\right)=exp\left[i\cdot \mathrm{\Delta }\Phi \left(j\right)\right],$$E6

which corresponds to retaining information about the phase increment only.

The complex-valued quantities $c\left(j\right)$ (or $c_{\mathit{norm}}\left(j\right)$) can also be considered as vectors and vertically averaged in the vector sense. If the vertical size of the processing window is smaller than the characteristic vertical scale of the strain inhomogeneity, the phase-variation increments $\mathrm{\Delta }\Phi \left(j\right)$ are nearly identical (but in reality may be distorted by decorrelation and other noises). Vector averaging over the vertical size of the processing window then gives a complex quantity with a much more stable phase $\overline{\mathrm{\Delta }\Phi }$:

$$\overline{c}\equiv \overline{C}exp\left[i\cdot \overline{\mathrm{\Delta }\Phi }\right]={\sum }_{j=1}^{N_z-1}c\left(j\right)/{\left|c\left(j\right)\right|}^{\alpha }.$$E7

Summations in Eqs. (4) and (7) actually correspond to obtaining of averaged real and imaginary parts of the complex-valued signals $b\left(j,m\right)$ and $c\left(j\right)$, with subsequent singling out the resultant phase instead of direct averaging of individual phases for $b\left(j,m\right)$ and $c\left(j\right)$. As is clear from Figure 1, the vector summation makes the phases of the averaged vectors especially tolerant to small-amplitude erroneous vectors and signals with especially strong phase errors $\sim \pi$ rad which may occur for individual summands $b\left(j,m\right)$ or $c\left(j\right)$.

For the vertical inter-pixel distance $d_p$ in the imaged tissue, the so-found vertical phase increment $\overline{\mathrm{\Delta }\Phi }$ is proportional to the sought strain $\partial U/\partial z=\gamma \overline{\mathrm{\Delta }\Phi }/d_p=\overline{\mathrm{\Delta }\Phi }\left({\lambda }_0/4\pi d_p^0\right)$ and does not depend on the refractive index. Indeed, according to Eq. (2), the coefficient $\gamma ={\lambda }_0/4\mathit{\pi n}$, and the inter-pixel distance in the material is $d_p=d_p^0/n$ (since for a material with refractive index n, the same phase of the probing optical wave is accumulated along $n$ times smaller distance than in vacuum).

The above-considered averaging procedures are well applicable for fairly horizontal plane-parallel phase-variation isolines. However, in many case of practical interest, noticeable lateral inhomogeneities of the strain distribution may occur, so that the interframe isophase lines may be noticeably inclined within the processing window. Consequently, the directions of vectors $b\left(j,m\right)$ in Eq. (4) may strongly differ as a function of horizontal index. Thus, the straightforward lateral averaging like in Eq. (4) instead of improvement may significantly distort the visualized strain distribution.

To suppress this negative effect of averaging and retain advantages of the vector method, the following improvements can be proposed [25]. At the first stage, the complex-valued interframe signal $b\left(m,j\right)$ is averaged over small regions (having lateral and horizontal sizes $N_{\mathit{\text{small}}}\sim 2-3$ pixels only (see Figure 2)). This initial averaging does not yet noticeably affect the resolution but yields less noisy interframe matrix $b\left(m,j\right)\to \tilde{b}\left(m,j\right)$ by suppressing distorting contributions of weakest amplitude and most noisy pixels. At the next step, one creates a matrix containing vertical increments of interframe phase variations (similarly to Eq. (5) but without horizontal averaging over the processing window):

Figure 2.

Schematic of evaluation of the axial phase gradient in the vector approach adapted for laterally inhomogeneous strains with non-horizontal phase-variation isolines. All intermediate transformations are performed with complex-valued signals, and the sought phase gradient is singled out at the last stage.

$$d\left(m,j\right)=\tilde{b}\left(m,j+1\right){\tilde{b}}^{{}^{\ast }}\left(m,j\right)$$E8

To avoid confusion in the notations instead of array $c\left(j\right)$, here we introduced the matrix of complex-value quantities $d\left(m,j\right)$. The so-constructed matrix $d\left(m,j\right)=D\left(m,j\right)exp\left[i\mathrm{\Phi }\left(m,j\right)\right]$ contains the sought vertical inter-pixel increments $\mathrm{\Phi }\left(m,j\right)$ of the interframe phase variations. The vertical increments $\mathrm{\Phi }\left(m,j\right)$ of the interframe phase variations should already be nearly identical even if the isolines of the interframe phase variations are inclined within the processing window. Therefore, the quantities $d\left(m,j\right)$ can be efficiently averaged in the horizontal direction even for noticeably inclined interframe phase-variation isolines. Thus, the next steps are performing the vector averaging of $d\left(m,j\right)$ over the processing window with the size $N_x\times N_z$ pixels. For this procedure, the order of averaging over indices $m$ and $j$ is not essential, so that we first average over index $m$ and then over $j$.

Note that, by analogy with the previously considered normalization of vectors $c_{\mathit{norm}}\left(j\right)$ containing vertical phase-variation increments, the utilization of normalized quantity $d\left(m,j\right)$ can be useful, so that with normalization the averaging over index $m$ takes the form

$$\overline{d\left(j\right)}={\sum }_{j=1}^{N_x-1}d\left(m,j\right)/{\left|d(m,j)\right|}^{\alpha }$$E9

Here, notation $\overline{d\left(j\right)}$ is introduced to denote averaging over the horizontal dimension of the processing window. By analogy with Eq. (7), it was verified that the simplest choice $\alpha =1$ yields quite good results.

Similarly, averaging over the vertical index $j=1..\left(N_z-1\right)$ yields better results with additional normalization:

$$\overline{\overline{d}}={\sum }_{j=1}^{N_z-1}\overline{d\left(j\right)}/{\overline{d\left(j\right)}}^{\beta }.$$E10

Here, normalization exponent $\beta =1$ can also be recommended. Equation (10) represents a complex-valued quantity $\overline{\overline{d}}=\overline{D}exp\left[i\overline{\mathrm{\Delta }\Phi }\right]$ (averaged over the processing window), in which $\overline{\mathrm{\Delta }\Phi }$ is the sought vertical inter-pixel phase variation for the current position of the processing window. Phase increment $\overline{\mathrm{\Delta }\Phi }$ is directly analogous to the angle estimated in Eq. (7) and is also similarly linked to the sought axial strain. Geometric interpretation of Eqs. (8)–(10) is illustrated in Figure 2.

To illustrate the applicability of the vector approach in the most clear form, we first consider a numerically simulated example in which the reference and deformed OCT scans are simulated using model [27]. The simulation parameters for Figure 3 correspond to a typical OCT system, for which the central wavelength is 1300 nm, spectral width 90 nm, and A-scan of 256 pixels corresponds to the imaging depth 2 mm in air. The inhomogeneity of scattering amplitudes in the simulated structural image (Figure 3a) are chosen similar to typical experimental scans (compared with Figure 6). In the simulations initially 1024 scatterers were randomly distributed over each A-scan, and then the initial positions of the scatterers were displaced according to the assumed strain distribution shown in Figure 3b. Figure 3c shows the color map for interframe phase variations corresponding to the strain distribution in Figure 3b showing pronounced lateral inhomogeneity and tilting of isophase lines.

Figure 3.

Simulation of OCT-based mapping of laterally inhomogeneous strain field based on model [27]. (a) Is the simulated structural image; (b) shows the strain distribution adopted in the model; (c) color map corresponding to the interframe phase variations shown in (b); (d) is the reconstructed strain map using the vector method with straightforward horizontal averaging showing good results in laterally fairly homogeneous areas but prone to errors in the regions of inclined isophase lines; (e) is the reconstructed strain map based on the modified vector method adapted to averaging in regions of inclined isophase lines. The processing window is 16 × 16 pixels in size and the preliminary averaging area for panel (e) 2 × 2 pixels. In this example only strain-induced decorrelation noise is present.

The reconstructed strain obtained using the above-described vector approach is shown in Figure 3c and d. The processing window size is 16 × 16 pixels. The noisy areas in Figure 3d correspond to the regions of strongly inclined isophase lines, for which the straightforward lateral averaging Eq. (4) worsens the quality of elastographic mapping. Figure 3e shows much better strain-reconstruction quality for the modified vector method corresponding to Eqs. (8)–(10). We note that for Figure 3, only strain-induced decorrelation noise is taken into account. Figure 4 illustrates the vector method tolerance to other measurement noises. The latter were simulated by adding to each pixel of the image random complex-valued numbers with Gaussian distribution in order to obtain a preselected ratio between the average intensity of the OCT image and the added noise. Figure 4 demonstrates that the vector approach ensures fairly satisfactory reconstructions of strain maps down to average SNR∼1 (i.e., 0 dB) over the image area.

Figure 4.

Illustration of the vector method tolerance to measurement noises (other than strain-induced decorrelation noises) for the same simulated data as in Figure 3 but with additional noise (in the form of Gaussian random complex values added to each pixel). Panels (a) and (c) show the interframe phase-variation maps in the presence of noises for SNR = 6 dB and SNR = 0 dB, respectively. Panels (b) and (d) show the corresponding strain maps reconstructed using the vector method adapted for inclined phase-variation isolines.

4. Examples of the vector method application for strain mapping by elastographic processing real OCT scans

To illustrate operability of the vector method for elastographic processing of real phase-sensitive OCT scans, we demonstrate some examples of mapping thermally induced interframe strains which were produced in samples of biopolymers (collagenous samples of cartilage and eye cornea) by pulse-periodic irradiation of the sample by an infrared laser operating at a wavelength of 1.56 μm that is efficiently absorbed in the water that is present in the tissue. The heating pulses had duration ∼ a few seconds and the temperature reached ∼50–60°C in the corneal tissue and ∼60–80°C in cartilaginous samples. In more detail the experimental conditions are described in papers [28, 29, 30]. The used custom-made OCT device had parameters close to those assumed for the simulated examples in Figures 3 and 4. A typical experimental configuration is shown in Figure 5. The measurements were made in contact mode such that the OCT-probe surface contacted the studied sample though an intermediate layer of translucent silicone with pre-calibrated Young modulus. Certainly, measurements in the noncontact mode when the studied samples had free boundaries were also possible. However, during the irradiation of the water-saturated biopolymers, the silicone layer played a useful auxiliary role to protect the sample from drying during the heating. Furthermore, the main destination of such a layer was to play the role of compliant sensor to estimate the pressure exerted by the deformed silicone onto the studied sample (similarly to the discussion in [31]). OCT-based monitoring of strains in the pre-calibrated reference silicone and studied sample during mechanical compression of the silicone-sample sandwich was used to obtain stress-strain curves for the studied samples and estimate their Young modulus.

Figure 5.

Schematic of OCT imaging setup: 1 is the OCT scanner (λ = 1.3 μm); 2 is the soft silicone layer; 3 is the studied biopolymer sample (e.g., samples of cartilaginous or corneal tissues); 4 is the source of the heating infrared irradiation (erbium fiber laser operating at λ = 1.56 μm).

Figure 6 gives an example of experimentally reconstructed maps of interframe (“instantaneous”) strain obtained using the experimental configuration shown in Figure 5 for a rabbit cornea layer during one of the heating pulses. Structural image in Figure 6a demonstrates that the irradiated region exhibits visually appreciable local expansion after a series of heating pulses, and Figure 6b shows the interframe phase-variation map with pronounced inhomogeneity in the lateral direction with regions of very steep isophase lines. Similarly to simulated Figure 3, the strain map Figure 6c demonstrates the interframe-strain map found using straightforward horizontal averaging within the processing window (with sizes ∼1/20 of the entire image sizes), and Figure 6d is obtained using the method variant adapted for processing regions of steep isophase lines. Figure 6c and d clearly demonstrates that the irradiation-induced expansion of the heated water-saturated sample causes straining with the opposite sign (contraction) of the surrounding silicone layers. Figure 6e for one of the heating pulses demonstrates a pseudo-3D waterfall image showing the complex spatiotemporal evolution of the interframe-strain profile corresponding to the depth marked by the dashed line in the B scans in Figure 6d. We emphasize that the robustness of the vector method to measurement noises makes it possible to obtain rather clear quantitative strain maps for aperiodic strain evolution when periodic averaging for enhancement of SNR is impossible.

Figure 6.

Experimental demonstrations of OCT-based mapping of interframe strains. (a) is a structural OCT image of a biopolymer layer I (rabbit cornea sample) placed between translucent silicone layers (II); (b) is a typical color map of interframe phase difference, where the heated region is characterized by pronounced lateral inhomogeneity; (c) and (d) are the strain maps obtained with straightforward horizontal averaging over the processing window and using the vector method adapted for processing tilted isophase lines, respectively; (e) is a waterfall image showing the time evolution of the interframe-strain profile at the depth labeled by the dashed line in panel (d).

It can be said that in the above examples, strains were “instantaneous,” corresponding to the time interval between the neighboring scans (for the examples in Figure 6, this interval was 50 ms). However, in other biomedical and technical applications, it may be interesting to monitor cumulative strains over larger time intervals, during which the studied sample may be subjected to an external action (like laser irradiation in the above examples) or may exhibit other structural variations, e.g., chemical curing of polymers, influence of drying or impregnation by a liquid, etc.

The developed technique readily makes it possible to perform the summation of interframe strains in order to find the resultant cumulative strain. The issue of finding cumulative strains is discussed in more detail in [1] since this can be done in somewhat different ways, the result of which do not strongly differ for fairly small strains, but this difference may become appreciable for larger strains (∼10%). It can be said that the choice of the method of strain cumulation depends on the particular problem. For example, straightforward summation of interframe strains on larger interval may give a biased estimate of the total change in the sample thickness but may be preferable if the difference in strain over the imaged area should be used for estimating the differences in the elastic modulus in different image regions. Referring to [32] for details, for illustration purposes we will use simple straightforward summation of interframe strain. Figure 7a shows an example of the so-found evolution of lateral profile of cumulative strain (in the same experimental configuration as for evolution of interframe strains along a particular depth shown in Figure 6e). Figure 7b shows a representative graph for the time dependence of cumulative strain for a given lateral position. The profile demonstrates that during the heating, the tissue experiences internal structural changes with changes in the functional behavior upon reaching certain threshold temperatures (see peculiar inflection points in the dependence in panel Figure 7b). These changes are related to heating-induced active generation of pores as discussed in more detail in [30].

Figure 7.

Experimental demonstrations of OCT-based mapping of cumulative strains in the same experimental configuration as for Figure 6. (a) is a waterfall image showing the time evolution of the cumulative-strain profile at the depth labeled by the dashed line in Figure 6d; (b) is a representative graph showing the cumulative strain as a function of time for a particular lateral coordinate and depth.

In examples shown in Figures 6 and 7, we showed the results for only 2D mapping of evolving strain that evolved fairly rapidly, so that their 3D mapping via comparison between entire 3D sets of OCT data was impossible for a conventional scanning OCT system with a moderate rate of obtaining A-scans of ∼20–80 kHz. Utilization of a much faster scanning system may be also problematic for sufficiently long-time monitoring because of huge data flows and total amounts of OCT data that would require special means for signal acquisition and storage, as well as supercomputing performance for real-time processing. However, there are many processes of high interest, for which strains evolve much slower (e.g., drying of a sample, gradual curing of a polymer, various osmotic phenomena [33], etc.) To monitor such slow processes with acceptable time resolution, much larger time steps ∼1 sec and even greater may be sufficient. Furthermore, it can be shown that from the viewpoint of enhancing SNR in problems of monitoring of slow processes, temporal rarefaction of OCT data may be even more advantageous than acquisition with the maximal rate (see details of the corresponding discussion in [34]).

In such situations, even systems with moderate acquisition rate can be quite sufficient for realization of efficient monitoring of slow deformations. Figure 8 gives an example monitoring of fairly slow strains for drying cornea of a porcine eye (total duration of the record in 20 min with 1 s time step). The development of pronounced shrinking of the drying near-surface layer is clearly seen.

Figure 8.

Monitoring of slow-rate strains arising in the near-surface layer of porcine eye cornea due to its drying in open air. Panel (a) is a typical structural OCT scan of the near-surface cornea region; (b) is the time dependence of interframe strains found with averaging over the rectangle area shown in panel (a); (c) is an example of the spatially resolved map of cumulative strain after 20 min of drying; (d) is the time dependence of cumulative strain corresponding to interframe strain shown in (a).

Figure 9 shows a 3D example of monitoring of slow strain caused by mechanical relaxation of a cartilaginous sample, in which initially a local compression was produced by an inserted needle. In this example the strain was reconstructed via comparison of entire 3D sets of complex-valued OCT data (acquisition of one 3D data set required 1.6 s). The signal processing, as in the other previous examples, was made using the above-described vector method for estimating gradients of interframe phase variations. More detailed discussion of optimizations for the monitoring of slow strains can be found in [34].

Figure 9.

3D visualization of cumulative strains for a cartilaginous sample (porcine rib cartilage). Drying of the near-surface layer causes its gradual shrinking (negative sign of strain). The presented snapshots correspond a monitoring interval of 145 s. Panel (a) is the horizontal CT scan of cumulative-strain distribution corresponding to the plane marked by dashed lines in an isometric 3D view in panel (b). The localized area with positive strain sign is related to dilatation of the tissue area preliminary subjected to compression by inserting a needle.

5. Obtaining of quantitative stress-strain curves and estimation of Young modulus in phase-sensitive compressional OCE

Consider now possibilities of the developed approach to mapping strains for quantitative mapping of the Young modulus in the studied samples. As was discussed in paper [7] related to ultrasound-based elastography and pointed out in the Introduction to this chapter, if a rigid piston compresses a material that can fairly freely expand in the lateral directions, the stress distribution near the piston can be close to uniaxial. Consequently, the reaction to such compression is determined by the Young modulus of the material. This statement is the basic principle of compression (or compressional) quasistatic elastography [8]. Thus, the Young modulus can be estimated by measuring strain produced by a known stress applied to the sample by a compressing piston. It was also mentioned in paper [7] that a reference pre-calibrated layer overlaying the studied material can be used as a kind of compliant sensor to control the stress by measuring the strain within the reference layer. This idea has not found application in the ultrasound-based compressional elastography (where the compression approach is mostly used for visualization of relative distribution of stiffness); however, the application of a reference layer has proven to be useful in OCT-based quantitative elastography [31, 35]. In [35] special attention was paid to the effects of stiction between the compressing rigid piston (usually, the output window if the OCT probe) and the reference silicone layer. The stiction impedes free lateral expansion of the compressed layer. Consequently, its compressibility becomes lower (in other words, the distorted apparent Young modulus seems to be greater), because under constrained lateral expansion, the layer reaction to loading is determined by some mixture of the Young modulus and much greater bulk modulus (the proportion of their contributions being determined by the strength of stiction). If the possibility to observe inhomogeneity of the local strain inside the translucent reference layer is ensured, this makes it possible to detect the presence of stiction and take the necessary measures to reduce its distorting effect (e.g., to add a lubricant liquid between the compressing solid surface and the silicone layer). For a non-scattering, very transparent reference layer, its strain can be estimated only by observing its total thickness, so that the strain distribution inside the layer is impossible to detect optically. Besides, the accuracy of strain estimation via changes in the total thickness usually is lower, so that the usage of translucent layers is preferable in practice. Panels (a) and (b) in Figure 10 illustrate the application of a reference layer as a compliant sensor. The interframe phase-variation distribution in Figure 10b demonstrates evident inhomogeneity of the vertical phase-variation gradient, although the material of the layer is mechanically homogeneous. The layer region contacting with the stiff surface of the OCT probe looks more stiff because of the silicone-glass stiction. Figure 10c and d demonstrates the interframe-strain profiles within the mechanically compressed reference silicone layer and the underlying tissue in the experimental configuration shown in Figure 5 for two different degrees of silicone-glass stiction. In Figure 10c pronounced inhomogeneity of strain is visible, which looks as if the Young modulus of silicone varied ∼4 times over the layer thickness. In Figure 10d obtained in the same configuration by adding a lubricating liquid between the glass and silicone, the stiction is strongly reduced, so that the distribution of strain within the layer becomes rather homogeneous (as should be expected for the homogeneous silicone). Thus, one should pay due attention to controlling the stiction in order to avoid possible strong distortions in the apparent Young modulus of the reference layer.

Figure 10.

Elucidation of principle of OCE-based quantitative compressional elastography and potential stiction-related distortions. Panel (a) schematically shows the reference silicone layer over the studied tissue; (b) is a real example of interframe phase variation with inhomogeneous phase-variation gradient within the homogeneous silicone layer; (c) is an example of pronouncedly inhomogeneous phase variation as a function of depth because of noticeable stiction at the silicone-glass interface; (d) is a similar graph obtained in the same configuration after adding a drop of lubricating liquid at the silicone-glass interface for reduction of stiction.

If sufficiently low stiction at the interface with the glass is ensured, by measuring strain within pre-calibrated silicone layers, one can estimate stress applied to the studied underlying material. In this context another essential point is whether it is possible to neglect the dependence of the Young modulus of the silicone on the degree of silicone straining. In other words, what is the strain range within which silicone behaves as fairly linear material with stiffness independent of the applied stress? This is an important point, because unlike acoustics with typical strains below ${10}^{-4}-{10}^3$, in mechanical tests of polymers and biological samples, their strains may reach several percent and even 10%. In this context, Figure 11a demonstrates the results of a kind of “self-calibrating” tests allowing to verify linearity of mechanical behavior of silicones. In this test, sandwich structures composed of silicone layers with different Young moduli were used. The stiffer layer experienced smaller strain for the same stress remaining within an expectedly more linear region. If both compressed materials are linear, then for one strain plotted against another, one should expect a linear function. The curves shown in Figure 11a obtained for various combinations of softer/stiffer silicones are fairly linear with a rather broad strain ranging up to several tens of percent. This confirms that silicone is a good candidate for linear reference material in which strain is proportional to stress.

Figure 11.

Panel (a) shows cumulative strain in one silicone layer plotted against strain in another layer. The numbers near the curves indicate the corresponding ratios of the Young moduli for the compared silicone types. The dependences in panel (a) remain fairly linear up to strains of ∼10%. Panel (b) shows another example of a pronounced nonlinear stress-strain curve, for which the strain in pre-calibrated silicone is recalculated in stress and the horizontal axis shows the strain in an excised sample of breast cancer. Note that the slope of this curve (corresponding to the Young modulus of the tissue) strongly varies from 400 kPa to 1400 kPa for the strain range in the tissue <2.5%.

Consequently, if the dependence of strain in a studied tissue plotted against strain in the reference silicone layer looks nonlinear, the reason for this is that the stress-strain dependence for this tissue is nonlinear. For pre-calibrated reference layers, this opens the possibility to obtain quantified stress-strain dependences as illustrated in Figure 11b. The example in Figure 11b demonstrates the importance of the tissue pre-straining: even for apparently insignificant pre-straining within 2.5%, the apparent Young modulus may vary several times because of possible nonlinearity of the studied sample.

Further, we demonstrate that the described OCT-based technique opens the possibility to observe spatially localized inhomogeneities of the material stiffness. The examples presented in Figure 12 are also obtained in the experimental configuration shown in Figure 5. Figure 12a and b shows structural images of the silicone-cornea-silicone sandwich before and after pulse-periodic infrared irradiation. The latter causes deformation (expansion) of the irradiated region such that after postirradiation cooling the residual expansion persists. It is slightly visible in the structural image Figure 12b and can be quantitatively estimated and clearly visualized as shown in Figure 12c. For the pre-compression used in the discussed experiment, the Young modulus of the corneal tissue is close to the Young modulus of the silicone (∼200 kPa). Consequently, before the irradiation the reconstructed distribution of the compressibility (inverse Young modulus) over the image found via comparison of spatial distribution of mechanically produced small straining ∼0.5% is fairly uniform and does not vary laterally in the cornea (see Figure 12d). In contrast, after laser heating, the compressibility of cornea demonstrates strong increase in the center of the heated region (see Figure 12e). This local increase in the compressibility of the heated region is attributed to the laser-induced appearance of crack-like microscopic pores that are not directly resolved in the OCT images, but due to enhanced deformability, they manifest themselves via reduction of the tissue Young modulus. Analysis of the post-heating local dilatation of the tissue and the complementary reduction in the Young modulus makes it possible to make quantitative conclusions about the volume content and averaged geometrical parameters of the pores. These results well agree with independent data of microscopic examination of the tissue (see details in [30]).

Figure 12.

Demonstration of complementary irradiation-induced residual strain and reduction in the Young modulus in the irradiated region of corneal sample. Panel (a) is a structural OCT image before irradiation, and (b) is a similar postirradiation structural image, in which slight expansion of the corneal tissue is distinguishable in the encircled zone. Panel (c) is an elastographic map of the postirradiation cumulative strain, where expansion of the heated region is clearly visible. Panel (d) is the elastographic map of spatially fairly uniform distribution of the inverse Young modulus (compressibility) before irradiation. Similar map (e) obtained after irradiation clearly demonstrates appearance of local increase in compressibility (reduction in the Young modulus) in the center of irradiation region.

6. Conclusions

The described elastographic OCT-based approach can be viewed as an optical counterpart of the elastographic approach [7] proposed in medical ultrasound. However, an important distinction is that in contrast to the initially proposed correlational tracking in ultrasound, in OCT the speckle tracking based on analysis of the scattered signal phase appeared to be more advantageous. The difficulties in realization of correlational speckle tracking in OCT arise because of strain-induced decorrelation of speckle patterns in OCT scans, so that super-broadband sources are required to reduce the strain-induced speckle blinking/boiling as discussed in detail in [19].

It can be shown that the phase-sensitive OCT-based speckle tracking is more tolerant to strains [21], so that phase-resolved tracking of axial displacements can be efficiently realized even using OCT systems with “typical” parameters. In problems of mapping strains, the difficulties in phase unwrapping for super-wavelength displacements of scatterers can be efficiently obviated by direct estimation of interframe phase-variation gradients within a chosen processing window. This can be done using the proposed robust “vector” method [24, 25]. The resolution of the resultant strain maps is mostly determined by the dimensions of the used processing window, the size of which should usually be at least 5–10 times greater (depending on the noise level) than the resolution scale in the initial structural OCT images to ensure sufficient SNR.

Additional ways of SNR enhancement can be based on application of stable periodic actuators producing strain and periodic averaging (e.g., see [11]; however, this solution cannot be used to determine the quasistatic Young modulus). Alternatively in the case of aperiodic strains, efficient averaging can be obtained by finding cumulative strains as discussed in [32], which is also a key point allowing for obtaining nonlinear stress-strain curves as demonstrated above. Additional improvements in terms of optimization of interframe intervals may be required in the case of long-term monitoring of slowly varying small strains [34]. In the context of biomedical applications, some examples related to biomechanics of such biopolymers such as cartilaginous and corneal tissues were given above (a more detailed discussion can be found in [28, 29, 30]). Other examples related to oncological applications can be found in [11, 12]. It is likely that monitoring of slow strains (as discussed in [34]) may be of special interest for technical applications related to studies of curing of polymers, etc.

Acknowledgments

The authors acknowledge the support of the Russian Foundation for Basic Research (Grant No. 18-42-520018 in part of the development of 3D data acquisition and Grant No. 18-32-20056 in part of the development of OCE-based reconstruction of nonlinear stress-strain dependences). Obtaining of the experimental data for cartilaginous and corneal samples was supported by the Russian Science Foundation (Grant No. 16-15-10274).