1.

Introduction

Following the promising extraction of various stiffness metrics from ultrasound scans of biological tissues (e.g., elastographic visualization of shear-modulus distribution), development of various elastographic techniques has been initiated in other areas of biomedical imaging (see review¹), including optical coherence tomography (OCT).² The increasing interest in using OCT for visualization and quantification of shear/Young moduli distribution stems from the fact that other biophysical tissue parameters, such as optical scattering properties or ultrasound velocity, usually demonstrate a lower range of variations.^2,3 Furthermore, compared with parameters typical of medical ultrasound and most other imaging techniques,¹ characteristic OCT spatial scales are an order of magnitude smaller (typically 1 to 2 mm in depth and several millimeters in the lateral direction, with spatial resolutions of 5 to $20\mu \mathrm{m}$). This higher resolution (albeit with smaller imaging depths) is advantageous for elastographic applications, for example, for obtaining better accuracy of detecting boundaries of tumors (known to possess higher stiffness^3,4) and in moving toward “optical biopsy” to better discriminate pathological and normal tissues ex-vivo^5,6 and in vivo.⁷ Another emerging OCT elastography area is ophthalmology,⁸ where OCT imaging capabilities are well suited to the spatial scales of eye structures; in fact, multimodality OCT (spatial, angiographic, polarimetric) is now widely accepted in ocular care as a standard diagnostic tool.^9–11

Starting from seminal works by Schmitt² by analogy with strain mapping problems typical in mechanical engineering applications^12,13 and medical ultrasound,¹ where correlational processing of deformed speckle patterns is successfully applied, different variants of correlational processing methods have been considered in OCT-based elastography.^14–19 However, the adaptation of digital image correlation methods to OCT has not been very successful, primarily due to the interferometric nature of speckles in coherent OCT images and the resultant additional effects of pronounced speckle decorrelation that are not common in photographic images (see, e.g., Refs. 20 and 21). In view of this, phase-resolved measurements have been attracting ever increasing attention in OCT elastography, where the phase variations between compared A-scans or B-scans intrinsically exhibit better tolerance to decorrelation related to deformation-induced speckle blinking and boiling caused by mutual motion of scatterers located within the coherence volume.²² The measured phase variations can then be used for calculating the axial displacements of the scatterers^23–26 in elastographic mapping or for Doppler shifts’ estimation in angiographic imaging.^22,27,28 However, in OCT elastography, one is often interested in local strains rather than total displacements of scatterers. This way, using the determined differences in strains produced by the OCT probe pressed onto the inspected biological tissue, the spatial distribution of the tissue stiffness can be visualized. This is the basic principle of compressional optical coherence elastography (OCE).

Usually, the phase data are somewhat corrupted by intrinsic measurement noises of the OCT scanner and decorrelation noises caused by displacements of scatterers in the deformed tissue. In particular, in compressional OCE, the compression-induced strains cause inhomogeneous displacements that depend on depth. The resultant strain-produced displacements within a given OCT A- or B-scan may vary from subwavelength values to those exceeding not only wavelength, but pixel scale as well (Fig. 1). Conventional phase variation measurement methods compare phases from pixels with the same coordinates in the reference and deformed scans.^23,24,26 This is reasonable in regions of fairly small displacements, where the compared pixels contain the same scatterers in the referenced and strained scans, so that the phase variations are directly related to the scatterer displacements. However, in scan regions corresponding to suprapixel displacements, the same location in the two compared scans may correspond to different scatterers, thus the phases of the two OCT signals can be unrelated, their difference essentially yielding a random value. Therefore, “regular” phase variations in regions of suprapixel displacements become strongly masked by such displacement-produced decorrelation effects, as illustrated in Fig. 1.

Fig. 1

Schematic elucidation of the origin of the displacement-induced “decorrelation noise” effects arising in phase measurements when the same pixels in the original and deformed scans are directly compared. Upper panel schematically represents pixelated images of two neighboring localized scatterers in the original OCT scan, and the lower panel shows their displaced positions due to pixel-scale displacement of the scatterers (at a certain depth corresponding to the one-pixel shift to the right). It is clear that direct comparison of pixels with the same coordinate in the reference and deformed scans (marked by the vertical arrow) is contaminated by cross-talk between neighboring scatterers with unrelated initial phases. Consequently, the determined phase difference may yield a random value even in the absence of any other sources of measurement noises (the actual new position of scatterer 2 now corresponds to the neighboring pixel marked by the inclined arrow).

The straightforward way of reducing such displacement-induced decorrelation noises is to choose sufficiently small strains (say, below ${10}^{-3}$) that produce limited subpixel displacements over the entirety of compared scans. However, even this may not improve the strain-estimation accuracy: strains are proportional to phase-variation gradients and for small strains, the small phase gradients become significantly corrupted by the measurement noises of the system. The latter can be reduced by applying auxiliary periodic small-amplitude actuators to enable averaging of compared OCT scans, but this is probably impractical for most in vivo OCE applications based on hand-held-probe operations. Therefore, one faces a dilemma—on the one hand, needing to ensure small displacements to avoid different scatterer cross talk in original-deformed scan comparisons (as per Fig. 1); but on the other hand, these resultant small phase differences may be overwhelmed by measurement noises, thus reducing the accuracy of phase-gradient estimations.

Alternatively, a moderate increase in strain (say, from $\sim {10}^{-3}$ up to $\sim {10}^{-2}$, for which pronounced strain-induced speckle “boiling” does not yet occur) may strongly reduce the influence of additive measurement noises (e.g., those of the photodetector). However, stronger strain-induced displacements will cause strong cross-talk of the displaced scatterers, as shown in Fig. 1, thus masking displacement-induced decorrelations. Another larger-strain complication is the effects of phase wrapping with increasing measurement depth, thus making the determined phase difference—displacement relationship ambiguous and necessitating error-prone unwrapping procedures.

To reconcile such contradicting requirements while retaining the advantages of using larger strains up to $\sim {10}^{-2}$ and avoiding the abovementioned complications, we develop a method in which the suprapixel displacements in deformed scans are tracked to ensure comparison of scattered fields corresponding to the same scatterers even if they move to the neighboring pixel and their displacement significantly exceeds the wavelength. We demonstrate that despite multiple phase wrapping and strong decorrelation noise arising if using the straightforward pixel-to-pixel comparison, the proposed phase gradient method yields significantly better accuracy and better tolerance to measurement noises compared to the limiting small-strain measurements. We emphasize that the proposed way of using larger strains (up to the onset of significant speckle blinking) appears attractive for practical compressional elastography based on hand-held probe operation in probable clinical conditions. Note that the application of periodic straining and use of time-averaging for improving SNR may also be possible but is arguably more challenging.

2.

Phase Comparison for Supra-Pixel Displacements Between Native and Deformed OCT Scans

To demonstrate phase comparison for complex-valued OCT scans suitable for suprapixel displacements of scatterers, we use our model²⁹ according to which complex-valued amplitudes of pixels in a one-dimensional (1-D) A-scan formed by spectral-domain OCT can be written as follows:

To demonstrate the phase variations caused by the axial displacement of scatterers, let us first consider a single subwavelength scatterer with initial coordinate $z_j$ that is displaced by $\mathrm{\Delta }z$. For simplicity and brevity, we assume that the total number of spectral harmonics $N+1$ is odd [i.e., in Eq. (1), $N$ is supposed even without loss of generality], whereas the spectrum shape $S(n)$ is rectangular and symmetrical around the center. Then, performing summation in Eq. (1), the 1-D signals $A^{(1)}(q)$ and $A^{(2)}(q)$ in the initial and $\mathrm{\Delta }z$-displaced pixelated scans are

The particular representation of Eq. (4) for the factor $T(g,N)$ is suitable for a rectangular spectrum shape $S(n)$, whereas for a smoother and more realistic (and likely effectively narrower) spectral shape, the sign may remain constant or may flip at a distance greater than one pixel. In any case, the phase change by a constant value does not affect the phase gradient that is of primary interest for estimating strain in elastographic applications. Therefore, even if for $\pm 1$ pixel, the function $T(g,N)$ changes its sign, this fact can easily be verified by comparing the phases of the reference scan and the same scan shifted by one pixel (in image processing), so that if necessary, this effect can be readily taken into account. In other words, it can be said that around a single scatterer, there is a kind of “frozen phase-variation hat” described by the real-valued function $T[g(\mathrm{\Delta }z=0),N]$ for the initial position and $T[g(\mathrm{\Delta }z\ne 0),N]$ for the displaced position, so that the phase difference between any pixels within these “hats” is described by the same phase factor $\exp (i2k_0\mathrm{\Delta }z)$.

Therefore, as schematically illustrated in Fig. 1, comparison of phases between the $q$’th pixel in the reference image can be made using not only the same pixel $q$ in the second scan but also pixels $q\pm 1$ or $q\pm 2$ (that better correspond to the position of the maximum of the “hat” in the displaced images for tracking of different methods can be used as discussed below). This pixel-scale tracking of the new positions of the scatterers can strongly reduce uncontrollable errors due to various noises (both decorrelation ones and noises from other sources). Indeed, it is clear from Fig. 1 that the errors related to leakage (cross-talk) of the fields from neighboring scatterers can be significantly reduced if instead of pixel $q$ in the second image, one considers the neighboring pixel corresponding to the displaced position of the tracked scatterer, say, ($q\pm 1$)’th pixel. This means that even for suprapixel displacements, the influence of the displacement-induced decorrelation noise can be reduced down to the much lower level typical of subpixel displacements. At the same time, phase gradients (and, therefore, strains) can easily be made an order of magnitude greater than for the commonly used small subpixel (and even subwavelength) displacements, so that one can strongly reduce the strain-estimation errors related to additive measurement noises (e.g., thermal noises in the photodetector). Thus, somewhat paradoxically, despite an unknown number of wrapped phase periods for suprapixel displacements, one can accurately measure strains by estimating strain gradients over relatively small depths intervals (say, $\sim 10$ to 20 pixels) along A-scans with very high tolerance to noises, as demonstrated in Ref. 30. The above-mentioned property of the “frozen phase-variation hat” around the maximum of the point-spread function $T(g,N)$ enables one to extend the range of measurable strains up to $\sim {10}^{-2}$, relative to maximal strains of $\sim {10}^{-4}$ to ${10}^{-3}$ that limit conventional pixel-to-pixel phase comparison methods. Such an extension of the strain range to $\sim {10}^{-2}$ is sufficient to encompass the entire strain range over which OCT scans maintain sufficient similarity at a pixel-scale level; for strains $>{10}^{-2}$, OCT images become strongly corrupted by speckle blinking and “boiling” caused by sufficiently strong mutual displacements of scatterers located within the same pixel; detailed discussions of similarity/difference between displacement- and strain-induced decorrelation of OCT scans can be found in Refs. 17, 2021.–22, 29, 30.

Returning to Eqs. (2)–(4) obtained for a rectangular spectrum with the odd total number $N+1$ of optical spectral components symmetrically located around $k_0$, one can easily see that for even total number $N$ of spectral components (for which there is no central component $k_0$), the same trends hold. That is, the function $T(g,N)$ remains real, and the only difference is that the central wavenumber in the “rapid” phase factor $\exp (i2k_0\mathrm{\Delta }z)$ becomes slightly renormalized, $k_0\to k_0+\pi /2H$, because of the absence of a “genuine” central component. For smooth spectral shapes with decreased spectral amplitudes at the spectrum wings, $S(\pm N/2)\to 0$, the difference between odd and even numbers of harmonics becomes even much weaker and can be ignored. Eventual slight asymmetries of $S(n)$ result in only a slight shift of the effective central wavelength (wavenumber $k_0$) and in a somewhat modified shape of the real-valued $T(g,N)$. Therefore, the considered property of “frozen phase-variation hat” is quite general.

In the following sections, we illustrate these phase variation properties by simulated and experimental examples, and demonstrate how these can be efficiently used for improving the phase-gradient estimate accuracy. In particular, the mentioned possibility of comparing phases in different pixels will be used for making compensating shifts (in combination with additional optimizations considered below); note that the entire developed phase processing methodology is shown in a summary flowchart in Fig. 9. However, before focusing on these optimizations, we show that the effect of the “frozen phase-variation hat” should be taken into account with care when estimating strains in the vicinity of bright highly scattering image features, near which artifacts in the form of false stiff inclusions may appear in OCT elastograms.

3.

Apparent Stiff Inclusion Artifacts near Bright Scatterers due to Amplitude-Phase Interplay in Optical Coherence Tomography Phase Gradient Elastography

Although phase measurements in OCT are usually considered as essentially decoupled from amplitude measurements, and amplitude approaches are contrasted to phase ones, here, we demonstrate some interdependence of the two. Specifically, the above-discussed properties of phase variations in the vicinity of displaced scatterers with significantly enhanced brightness can result in the appearance of “stiff inclusions,” i.e., artifacts in the form of apparent stiff regions in elastographic images obtained by comparing complex-valued OCT scans in compressional elastography. The physical origin of such artifacts is schematically elucidated in Fig. 2, where a fragment of an A-scan containing a bright scatterer near weakly scattering neighbors is shown. Even if the displacement is produced by a fairly uniform strain, the amplitudes corresponding to strong scatterers can “leak” to neighboring pixels and significantly exceed the reflectivity amplitudes of weaker scatterers located in those pixels. Consequently, the apparent phase variation in the vicinity of such a localized bright scatterer with initial position $z_0$ is dominated by its field [see Figs. 2(a) and 2(b)]. Correspondingly, the phase variation in this group of pixels is almost identical and equal $2k_{0}d(z_0)$, because it is dominated by the displacement $d(z_0)$ of the bright scatterer rather than by the actual displacements of the local weaker scatterers. For the assumed uniform strain $s=\mathrm{const}$. corresponding to the linear dependence of displacements of scatterers on their coordinates, $d(z)=s·z$, the phase variations of the scattered fields should also represent a linear function $2k_{0}s·z$, as shown by the dashed line in Fig. 2(c), with the constant gradient (derivative) as shown by the dashed line in Fig. 2(d). However, the masking field of the strong scatterer with nearly constant phase variation $2k_{0}d(z_0)\approx \mathrm{const}$. around its position $z=z_0$ forms a kind of plateau in the phase-variation dependence [solid line in Fig. 2(c)]. Correspondingly, the phase-variation gradient in the vicinity of such strong scatterers remains very weak [see solid line in Fig. 2(d)], erroneously suggesting the presence of very stiff inclusions.

Fig. 2

Elucidation of how the group of high-amplitude pixels around a strong scatterer can result in the appearance of an apparent stiff ($\sim \text{undeformable}$) inclusion, even if the actual strain is uniform (and the displacement of scatterers with initial coordinate is a linear function). (a) and (b) the initial and shifted positions of such a strong scatterer and several weaker scatterers in its vicinity are shown. Dashed line with constant slope in (c) shows the unwrapped phase shift expected for uniform strain; solid line shows how the group of high-amplitude pixels around the strong scatterer with phase variation leads to a plateau in the observed phase dependence. (d) The derivative (gradient) of the phase dependences in (c): dashed line for the case of uniform strain without a strong scatterer, and solid line in the presence of a localized strong scatterer is shown. To summarize, interpixel cross talk near a strong scatterer masks the fields of weaker scatterers nearby, thus forming a flat plateau in the depth dependence of local phase variation. Correspondingly, the phase gradient (and thus, the apparent strain) in the region of this plateau is nearly zero, as shown in (d), suggesting the presence of a nonexistent stiff inclusion.

Experimental and simulated (based on model²⁹) demonstrations of this stiff-inclusion artifact are presented in Fig. 3. In the experiments, we used a custom-built spectral-domain OCT system operating at a central wavelength of $1.3\mu \mathrm{m}$, with axial resolution $\sim 8\mu \mathrm{m}$ in air, focus beam diameter $\sim 15\mu \mathrm{m}$, 20-kHz rate of obtaining spectral fringes, with phase stability better than $\sim 0.01\text{rad}$ (i.e., sensitive to $<1\mathrm{nm}$ displacements of scatterers, as quantified in well-controlled vibration isolation experiments with piezoactuators). However, in the examples discussed here, the measurements were made without advanced vibration isolation typical of biomedical labs, but rather under typical “clinical-like” conditions with manual deformation of samples. Thus, minimal measurable strains (usually on order $\sim {10}^{-5}$) were influenced by basic factors, such as ambient vibrations and photodetector noises, rather than being limited by the ultimate phase stability of the scanner.

Fig. 3

Demonstration of the artifacts in the form of apparent stiff inclusions (indicated by arrows in the strain maps) appearing near bright localized scatterers due to the “frozen phase-gradient” effect. The three-layer sample has layers I = transparent and stiff, II = transparent and soft, and III = scattering and medium stiffness. (a) Measured structural OCT scan of a three-layer phantom with strongly increased brightness of the interfaces. (b) Strain map obtained by fitting the gradient of phase difference between the reference and deformed scans. The “vector” method used here for gradient fitting is described in Sec. 5. (c) Simulated structural image based on the model²¹ in which the distribution of brightness of the scatterers is chosen to imitate the real OCT scan in (a), and the contrast in the stiffness of the layers also approximates the real sample. (d) Corresponding strain map calculated from the simulated image in (c). Note the similarity in the experimental (b) and simulated (d) strain maps, including the appearance of artifactual stiff layers near the bright interface between deformable layers II and III and localized stiff inclusions between stiff layer I and deformable layer II.

B-scan rate was 20 frames per second and the OCT probe displacements were manually controlled with characteristic interframe strains in the range between $10{}^{-4}$ and $10{}^{-2}$, which made it possible to use a single pair of B-scans for estimation of phase-gradients using an unconventional and robust method that is discussed in Sec. 5 (see also the summary flowchart in Fig. 9).

Figure 3(a) shows an experimental 2-D B-mode structural scan of a three-layer phantom pressed from above by the OCT probe; Fig. 3(b) shows the reconstructed strain map obtained by comparing a pair of initial/deformed structural scans. The phantom consists of a rather transparent and very stiff upper plastic layer (I) pressed by the output window of the OCT probe, the much softer and also relatively transparent layer (II), and a scattering layer (III) of intermediate stiffness [$3\times$ stiffer than layer (II)]. It is clear from the structural scan of Fig. 3(a) that layer (III) is much brighter that other layers owing to its stronger scattering. The bright interface between layers (II) and (III) thus creates the above-mentioned “frozen phase-gradient” zone. This yields an artifact in the form of an apparent stiff and thin interface between layers (II) and (III) in the calculated strain map shown in Fig. 3(b). Figures 3(c) and 3(d) show the corresponding simulation results using our model that describes the speckle-pattern evolution in OCT images of deformed tissues.²⁹ Equation (1) of the present paper with $N=256$ spectral components was used to generate B-scans containing 256 A-scans with input experimental system parameters, to compare directly with real OCT images of the left panels of Fig. 3; the optical and mechanical phantom properties used in the simulation also correspond to the experiment. The strain map thus obtained in Fig. 3(d) displays similar features to those derived from the experimental measurements of Fig. 3(b), importantly including the abovementioned effect of localized false stiff inclusion in the vicinity of bright scatterers [at the boundary between layers (II) and (III)]. Note that although the layer boundaries in the model are flat, the reconstructed strain map exhibits some undulations similar to the experimental results [Fig. 3(b)]. The reason for this is due to the fluctuating density of randomly distributed scatterers in the simulated scans. Although we did not try to imitate the finer features of the real structural scan (such as existence of aggregates of scatterers surrounded by more transparent regions) and only roughly simulated the gradual decay of the signal with increasing depth, the similarity of the main features of strain maps in Figs. 3(b) and 3(d), including the appearance of the false stiff layer near the II-III interface, as well as artifactual stiff “teeth” penetrating from the stiff layer (I) into the more deformable layer (II) near scatterers with especially bright contrast, is encouraging.

Figure 3 clearly demonstrates that the “frozen phase-gradient” effect near bright scattering boundaries can indeed be rather pronounced and should be taken into account when interpreting the elastographic strain maps obtained via phase-gradient analysis. Note that normalization of the OCT-image brightness does not help to eliminate this effect—even in the normalized form, the optical field “leaking” from bright scatterers to the neighboring pixels representing weak scatterers still dominates, so that the “frozen phase-gradient” effects persist even in the amplitude-normalized images. This artifact near localized bright areas may thus be easily misinterpreted as the existence of localized stiff inclusions (e.g., microcalcification zones in tissues, and so on), so that OCT strain maps in such circumstances should be interpreted with caution.

The magnitude of the “frozen phase-variation” effect depends on such features as the optical spectrum shape (this determines the degree of localization of the point-spread function for localized scatterers) and signal processing specifics (size of the processing window that determines the area of smoothing in the strain map, the particular method of finding the phase-variation gradient, and so on). Some examples of how the gradient-finding procedures affect the quality of the resultant strain maps are given in Sec. 5.

4.

Strain-Map Improvements with Pixel-Scale Shift Corrections

OCE’s high robustness to additive noises was previously discussed for elastographic mapping,¹⁸ based on phase-gradient analysis with corrections of pixel-scale shifts and interpixel phase comparison. We now show further useful utilization of this technique for improved phase-gradient estimation. In particular, this pixel-scale compensation may help in handling the displacement-induced decorrelation effects, recognized as one of the main sources of error in phase-gradient estimates.³⁰ These may be important in compressional OCT elastography, for example, in the use of an intermediate layer of a reference material with calibrated stiffness to estimate the stiffness of the inspected tissue. We propose that this deformable layer can be placed between the inspected biological tissue and the rigid surface OCT probe, producing the tissue compression, as shown in Fig. 4. This approach is somewhat similar to the “optical palpation” method^31,32 and the use of a “compliant sensor” for quantification of tissue stiffness.⁵ In these works, the variation in the entire thickness of the transparent calibrated layer is directly measured to estimate its strain, in order to quantify the underlying tissue stiffness or indirectly estimate its spatial distribution (underlying biological tissue is not directly visible on OCT). However, allocating most of the OCT scan depth for directly visualizing the tissue is desirable, which could be enabled by a thinner reference layer ($\sim 100\mu \mathrm{m}$); but for thin layers, boundary effects (sticking at the contact with the rigid OCT probe) can strongly distort its apparent stiffness, so their use is impractical. Alternatively, here, the visualized biological tissue occupies the majority of the OCT scan, while the weakly scattering but thick reference layer is only partially visible in the upper part of the scan; this significantly reduces the artifactual effect of possible layer sticking to the OCT probe (Fig. 4). In this configuration, the strain map would represent the strain values in both the inspected biological tissue and the partially visible reference layer placed above the tissue. Due to the latter’s deformability (that can have a millimeter-scale thickness that is mostly located above the visualized region), the scatterers in the visualized region can be displaced by a supra-wavelength and even suprapixel scale between the initial and deformed positions. In such a case, the situation corresponding to the schematic of Fig. 1 can occur. Consequently, conventional comparison of phases of pixels with the same coordinates in the reference and deformed scans may become altogether impossible and/or give very noisy results due to cross-talk of the scattered fields between neighboring pixels. However, judicious integer-pixel correction of this displacement can significantly improve the quality of the phase comparison and correspondingly improve the accuracy of phase-gradient estimation.

Fig. 4

Schematic illustration how a calibrated reference (sensor) intermediate layer with known stiffness properties (for quantification of tissue stiffness values in compressional OCE) can produce pixel-scale shift of the entire imaged tissue region. Besides this shift as a whole due to compression of the reference layer (shown by three arrows), scatterers in the strained tissue experience additional depth-dependent displacements, as shown by the larger thick arrow. This introduces strong displacement-induced decorrelation preventing direct pixel-to-pixel comparison of the two scans; however, this problem can be handled by introducing compensating integer-pixel shift.

We will discuss the application of a calibrated reference layer elsewhere, and here, we demonstrate the usefulness of this pixel-scale shift compensation procedure in Fig. 5, using the same three-layer phantom examined in Fig. 3(a) but choosing another experimental run. In these experiments, to eliminate undesirable strong scattering at the interface between the upper stiff layer (I) and the glass window of the OCT probe, glycerol was used as immersing fluid, so that it could be squeezed out of the contact area upon compression, resulting in a translational pixel-scale displacement of the scatterers over the entire OCT scan. In Fig. 3, a strong pixel-scale interframe shift was absent, but for Fig. 5, we intentionally chose the case where the contact between the OCT-probe surface was especially loose; thus, the squeezing of the immersing liquid was stronger and resulted in $>\mathrm{one}\text{-}\mathrm{pixel}$ initial shift of the sample relative to the OCT-probe surface. This way we imitate translational shift due to squeezing of the reference layer shown in Fig. 4. Furthermore, compression of the deformable layers (II) and (III) produced additional pixel-scale displacements of scatterers in the lower part of the scan. Such significant overall displacements result in noticeable interscan decorrelation, necessitating some method(s) of compensation in the strain estimation routine.

Fig. 5

Estimating the phase-difference between reference and deformed OCT images with and without correcting pixel-scale shift [same phantom sample as in Fig. 3(a)]. (a) The phase-difference map found by straightforward comparison of pixels with the same coordinates in the reference and deformed images when the latter experiences about a one-pixel translational shift (because of squeezing out of the immersing glycerol between the OCT probe and sample surface). (b) The same phase-difference map after correcting pixel-scale shift (both the initial one-pixel translational displacement of the entire scan is introduced together with cumulative strain-induced one-pixel shift at the depths below $\sim 150\mathrm{px}$). (c) and (d) the corresponding strain maps without and with correcting shift, respectively, are shown. The application of the correcting shift appears especially useful in the region of the weakly scattering and stronger deformed intermediate layer, where the improvement is noticeable. Inset in (a) shows zoomed ($7\times 7\text{pixels}$ in size) view of the interframe shift $\sim 2\mathrm{px}$ and shape distortion of one of bright speckles in the initial structural image (the dynamic range of the color palette in the inset is 10 dB and conventional interpixel interpolation is used).

The manifestation of this decorrelation effect is illustrated in Fig. 5(a), where the direct same-pixel comparison phase difference map between the reference and deformed images is shown, using the conventional “rainbow” phase periodicity color map. The decorrelation noise in the phase-variation map is very pronounced and strongly masks the regular strain-induced phase variations. However, these errors in the phase-variation estimation can be significantly reduced by introducing compensating shifts with a proper integer number of pixels. Figure 5(b) shows the corresponding phase-difference map obtained with application of such a correcting one-pixel shift of the entire deformed image (a possible way of finding this correcting shift is discussed at the end of Sec. 5). It is clear that the resultant phase difference after correction becomes much less noisy. Certainly, in both Figs. 5(a) and 5(b), the phase difference is wrapped, so that its value over most of the scan is determined with an uncertainty proportional to an integer number of periods. However, this uncertainty does not affect the estimation of strain, because an unknown number of wrapped periods does not prevent correct estimation of the phase gradient distribution. Figures 5(c) and 5(d) show the corresponding strain maps without [panel (c)] and with [panel (d)] initial translational one-pixel displacement compensation of the entire scan, combined with a cumulative strain-induced correcting shift applied for depths greater than $\sim 150\mathrm{px}$. It is evident that both the phase variation map [panel (b)] and strain map [panel (d)] become significantly less noisy due to the used pixel-scale compensation of the displacements.

This difference is especially clear in the vicinity of the weakly scattering deformable middle layer, where the strain is $>1\%$ and decorrelation noise related to the mutual motion of scatterers in the same pixel is especially strong. However, application of the appropriate compensating shift in Fig. 5(d) made it possible to obtain the strain map with quality comparable to the experimental Fig. 3(b) (where was no initial suprapixel interframe shift causing strong displacement-induced interframe decorrelation). We emphasize that the experimental strain maps in Figs. 3 and 5 obtained by the proposed phase-resolved method demonstrate fairly good results even for strains $>1\%$ that cause noticeable speckle blinking and markedly reduce the interframe correlation coefficient (even after applying the compensating shift it remains $\sim 0.3$ to 0.4 or less). By contrast, processing based on digital image correlation is much less tolerant to speckle blinking and would be hardly helpful for determining strains in such conditions.^21,30 The applied “vector” method of determining the phase gradient within the processing window ($16\times 16\text{pixels}$ for Figs. 3 and 5), as well as its comparison with other methods of gradient estimation, is discussed in greater detail and summarized in Sec. 5 (see specifically the summary flowchart in Fig. 9).

5.

Comparison of Phase-Gradient Estimation Methods for Finding Strains and Strain-Induced Supra-Wavelength Shifts

An important issue of the proposed optimization of phase-variation estimation is evaluation of the strain-induced displacements to enable the needed compensation shifts. The conventional way of displacement measurement in phase-resolved OCT is to directly unwrap the phase-variation maps in order to find the cumulative phase variation and the corresponding cumulative supra-wavelength displacements. In principle, differentiation of the so-found displacement with respect to the depth coordinate can be used for finding strains, although this approach often yields very noisy strain results, because conventional unwrapping procedures are error-prone and subsequent differentiation further enhances the inaccuracies in the displacement reconstruction.

An alternative way to determine the cumulative strain-induced displacements, originally proposed in Ref. 30, is based on the direct estimation of strains via determining phase gradients within a chosen analysis window starting from the regions closest to the compressing surface (we assume that this is the OCT-probe surface), where the strain-induced displacements are minimal (subpixel) for all realistic compression levels, thus their compensation is not necessary. Note that the estimation of phase gradients automatically removes the need for the error-prone spatial differentiation step that plagues conventional methodologies. Then with increasing depth, the contribution of the previously found strains in the upper part of the scan is summed to determine the cumulative displacement at the current depth. If this value exceeds $1/2\mathrm{px}$, $3/2\mathrm{px}$, $5/2\mathrm{px}$, and so on, we then add a one-pixel vertical compensating shift at the depth in question in the deformed image. As discussed in detail in Ref. 30 and illustrated by numerous simulations (results not shown), this method allows one to obtain very smooth displacement fields even in low signal-to-noise environments. Indeed, the noise tolerance of this cumulative approach is significantly higher than that of conventional phase-unwrapping procedures.

An important issue in the proposed strain-determination methodology is the choice of the actual method of determining the phase-variation gradient. This can be accomplished, for example, using conventional procedures of least-squares fitting²⁴ of the depth-dependent phase-variation within a selected size of the processing window. Before fitting the slope of the axial (vertical) dependence of the phase variation, for every horizontal row of pixels within the processing window, it is useful to determine the average phase variation at that depth, for example, using the Kasai estimator^33,34 or a somewhat different but analogous procedure described below in this section. The least-square procedure²⁴ of fitting the vertical phase-variation to determine the phase-variation gradient [Fig. 6(a)] implies that within the vertical size of the processing window, the fitted phase is not wrapped. Since the window height is usually fairly small (e.g., $\sim 1$ to 30 px) compared to the entire scan of the order of hundreds of pixels, this constraint of no phase wrapping inside the window is much less restrictive than for no phase wrapping over the entire scan depth (typically an order of magnitude greater vertical extent). However, for larger strains of the order of $\sim {10}^{-2}$, the necessity of phase unwrapping even inside the processing window may already be required [Fig. 6(a)], which can introduce additional errors.

Fig. 6

Schematic illustration of different methods for estimating the vertical gradient of phase-variation within the processing window. (a) The least-squares estimation of the slope of the horizontally averaged within the window phase-variation function, where index enumerates the horizontal rows within the window (finding requires verification of eventual phase wrapping within the processing window). (b) The stepwise summation of interpixel vertical differences in the phase variation for determining the total phase variation over the processing-window size in the vertical direction (obviates the necessity of phase unwrapping at each step even if the phase is wrapped over the entire processing-window size). (c) The vector method based on stepwise summation of complex-valued interpixel increments that are represented as vectors in the complex-valued plane (also obviates the necessity of unwrapping at each step even if the phase is wrapped over the entire processing-window size). It is clear from the geometrical interpretation that the vector method is especially robust to large errors in phase variation close to $\pi \mathrm{rad}$, because the rotation of the incremental vectors by $\pi \mathrm{rad}$ (as illustrated by the dashed vector) does not affect the direction of the resultant summed vector.

Alternatively, another computationally efficient and fairly noise tolerant procedure that does not require phase unwrapping over the window extent may be advantageous. In this procedure, for each horizontal row within the processing window, we also initially determine the averaged phase difference between the reference and deformed B-scans. Then, after performing this horizontal averaging for each depth along the window centerline, one obtains a vertical array $\mathrm{\Phi }(j)$ of phase differences between the reference and deformed scans, where index $j=1\dots N_z$ enumerates the horizontal rows within the window [of $N_z$ pixels in the vertical (depth) direction]. The total phase variation $\mathrm{\Phi }(N_z)-\mathrm{\Phi }(1)$ required for estimating its gradient over the vertical extent of the window is

Fig. 7

Strain maps for the above described three-layer sample found for the same processing window size $16\times 16\text{pixels}$ using conventional least squares fitting of the phase-variation gradient for [panel (a), see Fig. 6(a)], direct averaging of the inter-pixel phase-variation increments [panel (b), see Fig. 6(b)], and vector averaging via Eqs. (8) and (9) [panel (c), see Fig. 6(c)].

To further improve the phase gradient finding procedure, we propose another approach [Fig. 6(c)], which also completely obviates the necessity of unwrapping, is computationally efficient, and produces smoother strain maps, thus avoiding the pronounced patchy reconstructions in the noisy regions. This method uses full complex-valued OCT signal amplitudes in the reference scan, $a_1(j,m)=A_1(j,m)\exp [i·{\phi }_1(j,m)]$, and in the deformed one, $a_2(j,m)=A_2(j,m)\exp [i·{\phi }_2(j,m)]$. For every position of the processing window, we first find the following complex-valued quantity containing the horizontally averaged phase difference $\mathrm{\Phi }(j)$ between the horizontal rows with numbers $j$ in the deformed and reference scans:

The horizontally averaged vertical phase increments $\mathrm{\Delta }\mathrm{\Phi }(j)$ appearing in Eqs. (7) and (8) can be processed (for example, via the least-square method²⁴) to find the window-averaged vertical phase gradient. However, we propose a somewhat different method that can give even better results as demonstrated below. It is based on processing the entire complex-valued quantities $c(j)$ or $c_{\mathrm{norm}}(j)$. We verified that the normalized phase factor [Eq. (8)], without any additional amplitude weighting, gave the best results. This observation makes sense – at the first stage of horizontal averaging with amplitudes via Eq. (6), contributions of weakest and most noisy pixels become significantly suppressed. Subsequent amplitude weighting of these already-averaged quantities $\overrightarrow{c}(j)$ effectively decreases the number of information-bearing signals available for phase-gradient estimation over the vertical size of the window. The resultant detriment from this decrease in the number of averaged signals dominates, so that otherwise useful amplitude weighting does not yield an improvement of the elastographic image. This is true both for conventional least-square processing of the phases $\mathrm{\Delta }\mathrm{\Phi }(j)$ and for the alternative processing of the complex-valued signals $c(j)$ [Eqs. (7) and (8)].

In the proposed method, the complex-valued quantities $\overrightarrow{c}(j)$ [or ${\overrightarrow{c}}_{\text{norm}}(j)$] are considered as vectors in the complex-valued plane, for which the interpixel phase increment $\mathrm{\Delta }\mathrm{\Phi }(j)$ corresponds to the angle in the polar coordinate representation. In the ideal case, the interpixel phase variation increments $\mathrm{\Delta }\mathrm{\Phi }(j)$ should be identical within the window size, although in reality, they fluctuate due to decorrelation effects and other measurement noises or physical inhomogeneities. However, the sum of these elementary incremental vectors over the vertical size of the window yields a much more stable phase angle:

Figure 7 shows the resulting strain maps for the sample used in Figs. 3 and 5, obtained with three different phase-gradient-estimation methods illustrated in Fig. 6 ($16\times 16\text{pixels}$ processing window, $\mathrm{\Delta }j=2$). Figure 7(a) is obtained using the conventional method of least-squares fitting of the phase-variation²⁶ [cf Fig. 6(a)], in which the necessity to unwrap the phase should be checked and then the error-prone unwrapping performed if necessary. Figure 7(b) is based on the direct averaging of $\mathrm{\Delta }\mathrm{\Phi }(j)$ as per Fig. 6(b) and Eq. (5) [the errors in the resultant strain map are somewhat different, but the image quality is roughly the same as in Fig. 7(a)]. Finally, Fig. 7(c) is obtained via vector averaging [see Fig. 6(c) and Eq. (9)] in which normalized elementary vectors [Eq. (8)] are used. It is clear that the vector averaging procedure is preferable, yielding cleaner and smoother strain maps (without sacrificing the spatial resolution) even in regions where the least-square fitting and incremental phase-gradient estimation generate numerous patchy defects. Similar to the incremental phase-difference method, this approach obviates the necessity of phase unwrapping and is computationally efficient.

In addition to phantoms, examples in Figs. 7 and 8 illustrate the advantage of the vector method of estimating the phase gradient and strain in tissue. Figure 8(a) reproduces the strain map from Ref. 30 obtained using the least-squares fitting of phase-variation gradient for a tissue sample excised during a gynecological procedure (the experimental studies of the excised samples were performed with obtained subjects’ consent and approved by the Ethics Committee of the Nizhny Novgorod State Medical Academy, Protocol No. 14 of December 10, 2013). The human gynecological tissue sample was $\sim 1\mathrm{mm}$ in thickness and $5\mathrm{mm}\times 7\mathrm{mm}$ laterally, was not fixed, and was examined within 2 h after excision. This specimen was surgically excised with biopsy forceps, such that the outer layers of the mucous tissue ended up on both sides (above and below) of the subsurface tissue with expectedly lower stiffness. Such a three-layer structure is clearly visible in the histology of Fig. 8(c) and in elastographic maps of Figs. 8(a) and 8(b), where the intermediate layer clearly exhibits greater deformability. In comparison with Fig. 8(a) obtained via least-square fitting of the phase variations, Fig. 8(b) now shows the improved strain map obtained with the vector method of phase-gradient estimation based on Eq. (9), in combination with normalization of elementary vectors according to Eq. (8) (it was verified that the use of normalization gave better results in this case).

Fig. 8

Strain maps obtained from processing a pair of compared OCT scans of human biological tissue excised during a gynecological operation. Panel (a) is from Ref. 22 where the least-square fitting was used for determining the strain gradient. Panel (b) is the result of applying the currently proposed vector method for phase gradient estimation, with the same size of the processing window $16\times 16\text{pixels}$ as in panel (a); note that the “patchy” strain errors are now significantly reduced. Panel (c) is also reproduced from Ref. 22, to illustrate that histology confirmed the presence of the three-layer structure revealed in the elastographic maps. [(a) and (c) reproduced with permission].

In the context of the pixel-scale compensation procedures described in Sec. 4, we note that summing the cumulative displacements of strains over depth (see also discussion in Ref. 18) is not able to account for possible translational (block or bulk) displacement of the entire scan that is not related to the strain field visible within the scan. Nevertheless, if such a translational displacement does occur (e.g., because of deformation of the reference layer that is not completely visualized or a translation shift due to another reason), it can be readily found by a conventional correlation search and thus compensated for. This can be rapid as it requires only an integer-pixel search and needs to be done only once at the beginning of processing (using a correlation window in the upper part of the scan, where the strain-produced displacements have not yet accumulated). The determined translational pixel-scale correction can thus be applied to the entire scan. After this, the additional depth-dependent strain-produced shift can be found using the described cumulative phase-gradient-estimation procedures. In fact, such a method of preliminary finding the “global” translational shift was used in generating the results of Fig. 5.

Summarizing the above-considered optimizations of phase-gradient estimation, the following flow chart in Fig. 9 shows the main steps including utilization of the proposed “vector” method and application of correcting shifts.

Fig. 9

Flow chart summarizing the main steps of the proposed optimizations of phase-gradient estimation including utilization of the proposed “vector” method and application of correcting shifts.

6.

Conclusion

The proposed optimizations of phase-gradient estimation between reference and deformed scans take into account some important but often neglected features of complex-valued pixelated images in low-coherence methods. Although it is true that phase and amplitude measurements are usually independent, the influence of amplitude of the OCT signal leaking to neighboring pixels may significantly affect the determined phase variations; this is particularly so in the vicinity of strongly reflecting image features, as illustrated in this study.

We emphasize that even for suprapixel displacements, it is possible to introduce appropriate compensating shifts that can diminish the role of the displacement-induced decorrelation noise, in contrast with direct pixel-to-pixel phase comparison methods. Such compensation can significantly enhance the accuracy and robustness of the estimated phase-gradients, both essential features for the success of the phase-resolved variants of compressional OCE. Consequently, here, we show that the proposed optimized methods of phase-gradient estimation can confer high tolerance of compressional OCE to various noise sources (see also Ref. 30) This is in part due to the possibility of using significantly increased strains and displacements, compared with conventional phase-resolved approaches that are usually limited to subpixel and even subwavelength displacements of scatterers.

Besides reducing the decorrelation noise by suitable integer-pixel displacement corrections, the proposed method of phase-gradient estimation eliminates the necessity of error-prone procedures of phase unwrapping that usually limit conventional approaches. Although our approach does not directly determine cumulative multiperiod phase variations via unwrapping, it does ensure reconstruction of not only local strains, but also cumulative strain-induced displacements on supra-wavelength and even suprapixel scales with high accuracies that are usually unattainable with conventional phase unwrapping. The proposed approach may be of special interest in situations for which it may be challenging to implement signal-to-noise enhancement by averaging strain maps via stable periodic deformation sources. Thus, this optimized phase-gradient-determination methodology may be especially relevant to compressional OCE using a hand-held OCT probe as the compression element, which is of direct interest for practical applications in clinic.

Finally, we also demonstrate the existence of the “frozen phase-variation” zones associated with bright scatterers that may yield artifactual stiff localized inclusions in the reconstructed strain maps. This phenomenon should be taken into account in interpretation of the resultant OCT elastograms.