1. Introduction

Optical coherence elastography (OCE) is an imaging technique that maps mechanical parameters (e.g., strain, stress, transverse-wave velocity, and resonance frequency) and mechanical properties of tissue (e.g., elasticity and viscoelasticity) into micro-scale images, termed micro-elastograms [1–6]. Phase-sensitive compression OCE is one of the most prominent variants of OCE, providing rapid image acquisition with micro-scale resolution over millimeter-scale fields-of-view. This approach has shown promise in a number of applications, particularly in oncology [7–11] and mechanobiology [12–14]. In addition, as the compressive load is co-axially aligned with the optical axis, phase-sensitive compression OCE is readily implemented in compact imaging probes [15–19], which holds promise for translation of OCE to pre-clinical and clinical applications. Phase-sensitive compression OCE comprises three main steps: measuring sample deformation (i.e., displacement) using phase-sensitive optical coherence tomography (OCT), calculating local strain as the spatial derivative of displacement, and relating strain to elasticity via a mechanical model. Thus, a key factor in determining the image quality of compression OCE is the method used to estimate local strain.

A main challenge in compression OCE is obtaining accurate strain estimation in the presence of noise in the measured displacement. In an optimally configured OCE system, displacement sensitivity is determined by shot noise [20] and, also, as has recently been shown, speckle brightness [21]. The most straightforward approach to estimating strain is the finite difference method [22], i.e., the calculation of the change in displacement with depth between two discrete depth locations. However, this results in noisy strain as the noise in measured displacement is amplified in the estimated strain [23,24]. To alleviate this, several more sophisticated strain estimation methods have been proposed [5,24]. In particular, in the last decade, linear regression methods have been extensively utilized [19,24–28]. For instance, in an early demonstration of vibrational OCE [25], which is effectively a variant of compression OCE, ordinary least squares (OLS) was used to estimate the slope of vibration amplitude with depth. In this approach, a range of consecutive vibration or displacement measurements in depth are used to estimate strain, often referred to as the fitting range (typically, 5-10 times larger than the OCT axial resolution). As OCT signal-to-noise ratio (SNR) is one of the main determinants of the accuracy and precision of the measured displacement at a given depth [20], weighted least squares (WLS) was then proposed, with larger weights given to pixels corresponding to higher OCT SNR [24]. However, phase wrapping (occurring when the phase difference exceeds 2π) imposes a limit on the maximum measurable displacement and, by extension, the maximum measurable local strain, requiring the additional step of phase unwrapping to extend the strain dynamic range. Phase unwrapping techniques are inherently sensitive to noise, which can cause artifacts in strain micro-elastograms. Alternatively, a variant of the finite difference method implemented in the complex plane, referred to as the vector method [29,30], has been proposed to estimate local strain by calculating the phase gradient (i.e., the gradient in phase difference calculated between consecutive pixels in the axial direction). As this phase gradient is typically smaller than 2π and the intermediate steps for calculating phase gradients are implemented in the complex plane, this method avoids phase unwrapping. Importantly, the vector method performs vector summation, i.e., a type of spatial averaging to increase phase stability by summing a set of complex OCT signals, providing improved image quality over the finite difference and OLS methods [29].

Whilst both WLS and the vector method have been demonstrated in the majority of recent compression OCE studies [5,13,24,28–30], a systematic comparison of these two methods has yet to be performed. Consequently, it is difficult to determine which method is optimal for a given application. A main challenge is that, until now, each method has been demonstrated on different samples and using specific parameters that are difficult to compare between the methods, e.g., the fitting range and size of the phase averaging window are main determinants of spatial resolution in WLS and the vector method, respectively, but it is not clear how they compare. To address this, we introduce a framework to quantitatively compare strain imaging metrics (e.g., strain sensitivity [24,31,32], strain SNR [24,31,33], strain resolution [34–36], and strain accuracy [21,32]) calculated on the same data acquired from phantoms. In addition, we demonstrate the utility of the framework by comparing the performance of a new phase unwrapping algorithm for WLS, fast phase unwrapping (FPU), to both the existing WLS method and the vector method. The existing implementation of WLS applies a phase unwrapping algorithm that corrects phase discontinuities in depth by subtracting an integer multiple of 2π from the weighted phase difference [37], termed weighted phase unwrapping (WPU). Whilst WPU accounts for the dependence of phase sensitivity on OCT SNR, it is sensitive to phase noise in regions of low OCT SNR, resulting in phase unwrapping errors that manifest as prominent, streak-like artifacts in strain micro-elastograms [15]. This issue is addressed by FPU, which identifies phase discontinuities by calculating phase gradients using Laplace operators that can be calculated using Fourier transforms to improve computational efficiency [38]. The three-dimensional (3-D) version of the algorithm was initially developed and implemented in a graphics processing unit for quantitative retinal blood flow measurements in Doppler OCT [39,40]. Here, we propose iterative-based 3-D FPU combined with Volkov symmetrization [38,41] to fulfill periodic boundary conditions of Fourier transforms so that the accuracy of the approach is optimum for OCE. We then perform a comparison between the existing WLS method, termed WLS_WPU, the vector method, and the new implementation of WLS, termed WLS_FPU. Our results show that whilst these three methods provide similar strain sensitivity, strain SNR, and strain resolution, WLS_FPU provides improved dynamic range of accurate strain estimation with an error of ∼4% that is ×11 and ×15 times smaller than the error measured for WLS_WPU and the vector method, respectively, for an expected strain of 2.5 mɛ. We also provide preliminary analysis of sub-resolution contrast in compression OCE, i.e., the capability to detect changes in strain occurring within the strain axial resolution and, also, how this varies for the different strain estimation methods. Furthermore, we present strain micro-elastograms of mouse skeletal muscle and human breast tissue. We demonstrate that strain imaging artifacts caused by phase unwrapping errors in WLS_WPU are avoided in the other two methods and that the higher dynamic range of accurate strain estimation provided by WLS_FPU translates to improved contrast in strain micro-elastograms.

2. Theory: strain estimation methods in phase-sensitive compression OCE

Consistent with the broader field of elastography, in compression OCE, the sample is commonly assumed to be linearly elastic and mechanically isotropic. Under these assumptions, the elasticity of the sample can be represented by Young’s modulus, E, which is defined as the ratio between the normal uniaxial stress, ${\sigma _{zz}}$, and the normal uniaxial strain, ${\varepsilon _{zz}}$,

Assuming that each pixel within the OCT field-of-view is mechanically isolated, E estimated at each pixel, i.e., local elasticity, is independent throughout the micro-elastogram. In early demonstrations of compression OCE, the assumption of uniform stress was commonly made [1,5,42], whereby strain is inversely proportional to elasticity. More recently, techniques have been developed to estimate local stress, providing a more direct estimate of elasticity [43–45]. In both cases, developing accurate strain estimation methods is a crucial determinant of image quality in compression OCE. In this section, we describe in detail the two most prominent methods, WLS and the vector method.

2.1 Weighted least squares regression

2.1.1 Displacement measurement

In phase-sensitive compression OCE, the OCT phasor difference B-scan, $b(x,z)$, can be calculated using the Kasai estimator [46,47] from complex OCT B-scans of the sample in unloaded, ${a_U}(x,z)$, and loaded, ${a_L}(x,z)$ states,

Combining Eqs. (3) and (4), the displacement sensitivity, ${S_{{u_z}}}$, is

However, phase difference is a circular variable and wraps when the phase difference exceeds 2π, creating ambiguity in the measured displacement. Thus, to increase the dynamic range of displacement and, therefore, strain estimated using WLS, phase unwrapping is required.

2.1.2 Weighted phase unwrapping (WPU) algorithm

The WPU algorithm has been described in detail in a previous study [37]. Here, we provide a brief summary. The algorithm was initially proposed for a common-path phase-sensitive compression OCE configuration [37], where the back surface of the imaging window in contact with the sample is used as a reference reflector. Furthermore, the imaging window is affixed to an annular PZT actuator and acts as a compression plate. In this configuration, the phase difference at the surface of the sample is zero, as the window and sample surface have the same displacement. For a sample under uniform compression, as the resulting magnitude of phase difference increases with depth, i.e., as the change in displacement between the imaging window and the sample increases, the direction of required phase unwrapping is known. In addition, as highlighted in Eq. (4), phase difference sensitivity is dependent on OCT SNR, which allows the phase difference at each pixel to be weighted by the corresponding OCT SNR.

The weighted phase differences at the first n_z pixels in depth are assumed to be free of wrapping events, where the mean of weighted phase difference in these pixels, µ_{Δϕ, z}, is used as a reference to determine if a wrapping discontinuity has occurred. Then, for subsequent depths, starting from n_z+1, each pixel is unwrapped by subtracting an integer multiple of 2π such that the difference between each weighted phase difference and µ_{Δϕ, z} calculated from the preceding n_z pixels is minimized. Similar to the axial phase unwrapping, in the subsequent lateral phase unwrapping, a sliding window with a size of (2m+1) × (2m+1) in the xy plane is specified to calculate the mean of weighted phase difference in the en face plane, µ_{Δϕ, xy}, used as a reference to determine if lateral phase unwrapping is needed. Each pixel in the en face plane is laterally unwrapped by subtracting an integer multiple of 2π to minimize the difference between the axially unwrapped phase difference at each pixel and µ_{Δϕ, xy}.

2.1.3 Fast phase unwrapping algorithm

Here, we propose to use 3-D FPU algorithm as an alternative phase unwrapping algorithm [39] for use in phase-sensitive compression OCE. The algorithm was initially developed for unwrapping blood velocity profiles in retinal vessels imaged with Doppler OCT [39,40] and a complete derivation of the method can be found there. The method can be regarded as one of the transport of intensity equation solvers described in detail in [38,48]. Compared to the original FPU method, we introduce two modifications here. The first is the use of iterative error correction as described in [38], and the second is the use of Volkov symmetrization [41] in the FPU algorithm. FPU calculates the estimated unwrapped phase, $\Delta {\psi _{est}}(\mathbf{r})$, by solving the Poisson equation,

Figure 1 presents an illustrative diagram of the modified FPU algorithm using experimental data acquired from a silicone inclusion phantom. After the calculations of $\Delta {\psi _{est}}(\mathbf{r})$ from the enlarged data, we extract one part (represented by a solid black frame in Fig. 1) for the subsequent steps of the FPU algorithm. As it is common that $\Delta {\psi _{est}}(\mathbf{r})$ may not be fully unwrapped from the first estimation, the phase unwrapping errors can be minimized by iterating Eq. (9) based on the condition that $\Delta \phi (\mathbf{r})$ and $\Delta {\psi _{est}}(\mathbf{r})$ should be the same in the complex plane, as highlighted in Eq. (7). Therefore, we define the error function for the estimation of the unwrapped phase difference, $Q(\mathbf{r})$, as

 

Fig. 1. Illustrative diagram of the FPU algorithm for phase-sensitive compression OCE on a silicone inclusion phantom. “Sym.” is the Volkov symmetrization that precedes Eq. (9). “Cut.” is the operation of cutting off data that was marked as a solid black frame on enlarged symmetrized data. “Angle” is the operator to extract phase difference from complex functions.

Download Full Size | PDF

As the number of iterations in Eq. (11) increases, $Q_{est}^{iter}(\mathbf{r})$ is reduced towards zero. When $Q_{est}^{iter + 1}(\mathbf{r})$ is equal to zero, the iterations stop and ${\Delta }\psi _{est}^{iter + 1}(\mathbf{r})$ represents the final estimation of unwrapped phase difference. In practical implementation, the FPU algorithm faces two main problems when using $Q(\mathbf{r})$ as the stop criterion. The first is that $\Delta {\psi _{est}}(\mathbf{r})$ can have different phase offset equal to an integer multiple of 2π, leading to an ambiguity of applying Eq. (10). The second is variation in noise that still exists after a number of iterations. Thus, we instead compute spatial variance of $Q_{est}^{iter}(\mathbf{r})$ as the stop criterion since it is insensitive to constant value difference between $\Delta \phi (\mathbf{r})$ and $\Delta {\psi _{est}}(\mathbf{r})$. Then, the iteration only stops when the variance of $Q_{est}^{iter}(\mathbf{r})$ reaches a threshold accuracy δ, which is obtained by measuring the noise floor of (wrapped) phase difference in the deep regions of experimental OCT scans.

2.1.4 WLS strain estimation

After disambiguating the relationship between phase difference and displacement using phase unwrapping, strain can be estimated using WLS. WLS has been shown to improve strain sensitivity over OLS by accounting for spatially varying OCT SNR at each pixe [24] over the fitting range Δz, which determines the axial resolution of strain estimation [24,34]. Specifically, WLS assigns a weight, ${w_i}$, to each displacement measurement, ${u_{{z_i}}}$, where ${w_i}$ is defined as the inverse variance of each displacement measurement,

A complete derivation of the analytical expression for WLS strain estimation is provided in [24]. Here, we directly provide the expression,

2.2 Vector method

Here, we provide a brief description of the vector method. A more detailed description can be found elsewhere [30]. In the vector method, strain is estimated from the finite difference between two OCT phasor differences from consecutive depths within $b(x,z)$. We refer to this finite difference as the OCT phasor gradient, $c(x,z)$. This method removes the need for phase unwrapping as the computation is performed in the complex plane, and the phase angle of $c(x,z)$, $\Delta {\phi ^{\prime}}(x,{z_i}) = [\Delta \phi (x,{z_i}) - \Delta \phi (x,{z_{i + 1}})]$, referred to as the phase gradient, is typically smaller than 2π. From Eq. (3), as ${u_z}(x,z)$ is proportional to $\Delta \phi (x,z)$, it follows that since ${\varepsilon _{zz}} = \varDelta {u_z}/\varDelta z$, then ${\varepsilon _{zz}}(x,z)$ is proportional to the $\Delta {\phi ^{\prime}}(x,z)$,

3. Experimental setup and procedure

3.1 Phantom fabrication, tissue preparation and histology

The bulk material of the inclusion phantom (Phantom 1) is a 3 mm thick cylinder with a diameter of 10 mm, which was made from Elastosil P7676 (Wacker, Germany) with Young’s modulus of ∼19 kPa, as measured by a compression testing device. To minimize the out-of-plane deformation which could reduce the accuracy of strain estimation in OCE [49], the inclusion was designed to be a strip with dimensions of ∼8 × 0.6 × 0.5 mm (xyz), and was made from Elastosil RT601 (Wacker, Germany) with measured Young’s modulus of ∼327 kPa. To provide optical contrast, titanium dioxide (TiO₂) particles (refractive index = 2.3, Sigma-Aldrich, Germany) with a mean diameter of 1 µm were mixed in the bulk and inclusion silicone materials (both of refractive index = 1.4) at concentrations of 2 mg/ml and 10 mg/ml, respectively. Separately, to investigate the accuracy of strain estimation methods, a mechanically homogeneous phantom (Phantom 2) with a thickness of ∼1 mm and a diameter of 10 mm was made from Elastosil P7676 (Wacker, Germany). A 2 mg/ml concentration of TiO₂ particles was incorporated in the silicone mixture for Phantom 2.

In addition, we performed OCE experiments on data acquired from freshly excised mouse skeletal muscle and human breast tissue. An extensor digitorum longus (EDL) muscle with dimensions of ∼10 × 3 × 1 mm was excised from a 10-month-old male wild-type (WT) mouse sacrificed via cervical dislocation under anesthetic, as per institutional ethics requirements. Similar to a previous OCE study [50], we encapsulated the EDL muscle in 3-D methacrylated gelatin (GelMA [51]), a clear photopolymerizable hydrogel (∼15 × 15 × 3 mm), to keep the tissue hydrated during scan acquisition, and more importantly, to reduce image artifacts induced by non-uniform compression due to uneven surface topology [52]. The small size of the EDL muscle made it challenging to obtain meaningful co-registration between histology and micro-elastograms. When we attempted this, we observed freeze artifacts that caused staining issues in tissue histology. Thus, the corresponding histology is not provided here. In addition, we scanned a human breast tissue specimen excised during a mastectomy surgery. Informed consent was obtained from the patient and ethics was approved by the Sir Charles Gairdner and Osborne Park Health Care Group Human Research Ethics Committee (HREC No: 2007-152). This specimen had approximate dimensions of 4 × 4 × 2 cm and was scanned within an hour of excision. After data acquisition, the specimen was marked by ink for orientation, dissected, placed in two cassettes, and fixed in 10% neutral-buffered formalin for at least one day before histology processing was performed. Paraffin embedding was performed prior to sectioning and staining with hematoxylin and eosin, following the standard histology protocols at the hospital.

3.2 Phase-sensitive compression OCE system

The compression OCE system used has been detailed in previous studies [8,32]. Hence, a brief summary is provided here. Compression OCE measurements were performed using a fiber-based spectral-domain OCT system (Telesto 220; Thorlabs Inc., USA) with a central wavelength of 1300 nm and a 3 dB bandwidth of 170 nm. The measured full-width-at-half-maximum (FWHM) OCT axial resolution is 4.8 µm (in air). For imaging tissue-mimicking phantoms and the encapsulated muscle tissue, a scan lens (LSM03; Thorlabs Inc., USA) with a measured lateral resolution of 7.2 µm (FWHM in air) and a maximum field-of-view of 10 mm was used. A scan lens (LSM04; Thorlabs Inc., USA) with a measured lateral resolution of 13 µm (FWHM in air) and a maximum field-of-view of 16 mm was used to image the relatively larger breast specimen for cancer margin assessment. An imaging window (Edmund Optics Inc., USA) with a diameter of 75 mm was bonded with wax to a ring actuator (Piezomechanik GmbH, Germany) with an internal aperture of 65 mm, enabling imaging and compression to be introduced from the same side of the sample. The OCT system was implemented in a common-path configuration where the bottom side of the imaging window that contacts the sample during scanning provides the reference reflection. Silicone oil was applied between the sample and the imaging window to reduce surface friction. For biological tissue scans, a compliant silicone layer (made from Elastosil P7676) with a thickness of ∼500 µm was placed between the imaging window and the underlying tissue sample to improve the uniformity of compressive loading throughout the sample [43,52]. To ensure uniform contact, a pre-strain (i.e., a bulk strain prior to actuation) of 20% was imparted to the layer-sample system via a motorized vertical translation stage (MLJ150; Thorlabs Inc., USA). At this pre-strain, the ring actuator was driven collinearly with the imaging beam and synchronized with the OCT B-scan rate, delivering a maximum stroke of 10 µm. The A-scan acquisition time was 20 µs. A single B-scan consisting of 1000 A-scans was recorded in an acquisition time of 20 ms.

To characterize the strain sensitivity, strain SNR, and strain resolution, we utilized methods defined previously [31,32,34]. Specifically, in this study, 1,000 A-scans per B-scan, and 1,000 B-scan pairs (unloaded and loaded) per C-scan were acquired from Phantom 1 over a 2.5 mm field-of-view with a lateral sampling density of 2.5 µm. To enable a more accurate characterization of strain accuracy, denser lateral sampling was performed to include more OCT realizations; 2,000 A-scans per B-scan were acquired over a 2 mm lateral extent (1 µm per pixel in x) and a total of 256 B-scans were acquired at the same spatial location in Phantom 2. Similarly, to reduce the effect of optical noise on strain accuracy measurement, temporal averaging of 16 phasor difference B-scans was performed, which resulted in a total of 8 averaged phasor difference B-scans. The middle B-scan (the 4^th B-scan) was used to characterize the accuracy of the strain estimation methods, described in Section 4.2. Scans of the encapsulated EDL muscle were acquired by taking 1,000 A-scans per B-scan, and 1,000 B-scan pairs (unloaded and loaded) per C-scan over a 10 mm field-of-view at a lateral sampling density of 10 µm. Additionally, for the breast tissue scans, 808 A-scans per B-scan and 808 B-scan pairs (unloaded and loaded) per C-scan were acquired over a 16 mm field-of-view at a lateral sampling density of 19.8 µm. We note that the OCT lateral resolution was sub-sampled for measurements on biological tissues to allow for rapid data acquisition by trading off resolution to reduce the impact of temporal factors (e.g., dehydration [53] or viscoelastic creep [54]) that could alter the mechanical response of tissue in the process of data acquisition. The imaging specifications used in this study are summarized in Table 1.

Table 1. Phase-sensitive compression OCE scanning parameters and FPU parameters used in this study.

View Table | View all tables in this article

To unwrap the phase difference using the WPU algorithm, n_z = 7 axial pixels were utilized to calculate the phase unwrapping reference along the axial direction (µ_{Δϕ, z}). After axial phase unwrapping, a window of 9 × 9 lateral pixels (m = 4) was specified to calculate the phase unwrapping reference along the lateral direction (µ_{Δϕ, xy}). To determine the optimal number of iterations in FPU, a threshold accuracy of δ was measured from the standard deviation of phase difference calculated within a spatial window in the xz plane (typically 20 × 30 pixels) at a depth of 1.4 mm in OCT scans. The measured δ and the corresponding number of iterations for different datasets are also summarized in Table 1. We note that different values of δ in Table 1 correspond to varying sample deformation that increases strain-induced decorrelation and amplifies the spatial variation of phase difference in the strain accuracy experiment. Based on Eq. (3), the axial displacement at each pixel was then calculated from the unwrapped phase difference using the WPU and FPU algorithms, respectively, prior to WLS strain estimation. In terms of the computational cost, using the homogenous OCT data as an example, it requires 146.43 ± 0.38 s, 142.39 ± 0.86 s, and 22.62 ± 0.43 s for reconstructing 3-D strain from phasor difference using WLS_WPU, WLS_FPU, and the vector method, respectively, in MATLAB (R2016a) with a workstation that has two Intel Xeon E5-2690 8-core CPUs, 192 GB of RAM, and a 1 TB solid-state drive.

3.3 Comparison of strain estimation methods

Figure 2 provides a comparison of strain micro-elastograms generated using WLS_WPU and the vector method from data acquired on Phantom 1, using processing parameters previously reported for both methods [30,37]. For WLS_WPU, a 1-D axial fitting range of 100 µm, corresponding to 30 axial pixels, was implemented to estimate strain. For the vector method, a pre-filtering window with a size of 3×3 pixels (xz) was applied to OCT phasor differences, followed by vector summation over 16 pixels of normalized phasor gradients in x and a subsequent vector summation over 16 pixels of normalized, laterally-averaged phasor gradients in z. In addition, the phasor gradient at each pixel was calculated between two pre-filtered phasor differences at a step of 2 pixels (i.e., n = 2).

 

Fig. 2. Comparison of strain micro-elastograms generated with WLS_WPU and the vector method using typical parameters described in previous studies for both techniques. (a) OCT B-scan (xz plane), and the corresponding strain micro-elastograms generated using (b) WLS_WPU with a 1-D axial fitting range of 100 µm and (c) the vector method with a pre-filtered window size of 3×3 pixels (xz) and a main processing window size of 16×16 pixels (xz). The OCT SNR plotted as a function of (d) lateral position in x and (f) axial position in z, respectively. The strain estimated using WLS_WPU (blue) and the vector method (purple) as a function of (e) lateral position in x and (g) axial position in z, respectively.

Download Full Size | PDF

Figure 2(a) shows the OCT B-scan of Phantom 1, and the corresponding strain micro-elastograms obtained with WLS_WPU and the vector method are presented in Figs. 2(b) and 2(c), respectively, where the difference in the strain texture is pronounced due to the different spatial averaging techniques used in both methods. Figures 2(d)–2(g) show the variation of OCT SNR and strain over the lateral and axial ranges indicated by dashed black lines in Figs. 2(a)–2(c), respectively. Compared to WLS_WPU (Fig. 2(b)), the lateral boundaries of the inclusion are more blurred in the strain micro-elastogram obtained using the vector method (Fig. 2(c)) due to the incorporation of lateral averaging. Consequently, in Fig. 2(e), the strain estimated by WLS_WPU along the x-axis is noisier than that of the vector method. Similarly, as the axial averaging window size of WLS_WPU is larger than that of the vector method, the axial boundary of the inclusion is more blurred in the case of WLS_WPU in Fig. 2(b), and the corresponding strain along the z-axis shown in Fig. 2(g) is relatively smoother, compared to the case of the vector method.

From Fig. 2, it can be seen that differences in the averaging techniques and window size used for different strain estimation methods make it challenging to analyze and compare performance. To enable a fair comparison, for the results presented in Section 4, we modified both strain estimation methods and, importantly, aimed to preserve the advantages of each method, by restricting the processing to two layers of spatial averaging: (1) pre-filtering OCT phasor difference using a 3-D vector summation window (2×2×7 pixels, xyz), which has been typically incorporated in the vector method in previous studies [30,45], and (2) estimating strain over a 100 µm axial range of consecutive displacement and phasor gradient (n = 1 herein) data points for WLS-based methods and the vector method, respectively. To facilitate a numerical comparison of imaging performance between the different strain estimation methods, we define the following imaging metrics: strain sensitivity, strain SNR, strain resolution, and strain accuracy, as summarized in Table 2.

Table 2. Imaging metrics used to compare strain estimation methods

View Table | View all tables in this article

4. Experimental results

4.1 Strain sensitivity, strain SNR, and strain resolution

Figure 3(a) shows an experimental OCT B-scan acquired from Phantom 1, where the top surface of the embedded stiff inclusion can be observed at a depth of 380 µm. The corresponding strain micro-elastograms generated using the vector method, WLS_WPU, and WLS_FPU are presented in Figs. 3(b)–3(d), respectively. We note that there are variations in strain different between the homogeneous regions of Phantom 1, i.e., at the top left and top right of Figs. 3(b)–3(d). This is most likely caused by the non-uniform deformation in the sample due to sample imperfections (e.g., uneven surface topology [55]) arising from the phantom fabrication process. The quality of the strain micro-elastogram obtained using WLS_FPU (Fig. 3(d)) is similar to WLS_WPU (Fig. 3(c)). However, the texture of the strain micro-elastogram obtained using the vector method in Fig. 3(b) is distinct to that of the WLS-based methods, where the strain at each pixel estimated by the vector method appears to be less stretched along the z-axis (i.e., the fitting direction), despite the same axial range for strain estimation being used. This is attributed to the difference between weighted (i.e., WLS) and unweighted averaging (i.e., the vector method) in strain estimation. Because of weighted averaging, the linear regression in WLS tends to fit towards those displacement pixels with high OCT SNR (i.e., large weights) over the axial range. In contrast, the second layer of vector summation in the vector method is performed on the normalized phasor gradients, which means that each pixel over the axial range contributes equally to the estimated strain. This effect caused by weighting can also result in distinct features in strain micro-elastograms generated using WLS and the vector method, respectively. For example, Fig. 3(e) shows the magnified OCT image of a feature, a clump of TiO₂ particles, highlighted by the dashed blue box in Fig. 3(a). This feature extends over an axial range of 30 µm, i.e., it is smaller than the axial resolution of the OCE system (FWHM, ∼71 µm [34]). The corresponding magnified region in the strain micro-elastograms obtained using the vector method, WLS_WPU, and WLS_FPU are shown in Figs. 3(f)–3(h), respectively. Noticeably, using the WLS-based methods, the feature is identified as a region of low strain, indicating a stiff “inclusion”. However, this feature is not present in the micro-elastogram obtained with the vector method. Figures 3(i) and 3(j) present the OCT SNR and displacement depth profiles from a single A-scan through the middle of the bright feature. Importantly, the depth range where the slope of the displacement is close to zero (i.e., zero strain), marked by dashed black lines in Fig. 3(j), corresponds to the location of the highest OCT SNR, indicating that, in this case, the WLS-based methods correctly detect the presence of the stiff feature, whilst the vector method does not. The locations of the dashed lines in Figs. 3(i) and 3(j) represent the top and bottom boundaries of the aggregate determined from the abrupt changes of OCT SNR ranging from 30 dB to 41 dB and 37 dB to 26 dB, respectively. This result demonstrates a coupling between optical and mechanical contrast brought about by the weighting used in WLS-based methods. For a stiff and highly scattering feature that is smaller than the fitting range (×3 smaller in this case), the fitting can be biased towards the bright region, resulting in the detection of the feature whilst the absence of this weighting precludes the detection of the feature using the vector method. It should be noted, however, that although the feature is detected using the WLS-based methods, as it extends over an axial length that is lower than the fitting length, it is not resolved, as can be observed by its elongation over a larger axial range than is the case in the OCT image. This phenomenon will be further discussed in Section 5.

 

Fig. 3. Characterizations of strain sensitivity and strain SNR using the vector method, WLS_WPU and WLS_FPU using a B-scan acquired from Phantom 1. (a) OCT image. The corresponding strain micro-elastograms obtained using (b) the vector method, (c) WLS_WPU and (d) WLS_FPU. (e) The magnified region-of-interest where a TiO₂ aggregate is resolved, marked by the dashed blue box in (a). The corresponding magnified strain obtained using (f) the vector method, (g) WLS_WPU and (h) WLS_FPU. The depth profiles of (i) OCT SNR and (j) displacement across the middle of the TiO₂ aggregate. The dashed black lines in both (i) and (j) indicate the location of the aggregate.

Download Full Size | PDF

To quantify the comparative performance of the three strain estimation methods, we used the imaging metrics defined in Table 2. We measured the strain sensitivity (S_mε), mean strain (µ_mε), and strain SNR (SNR_mε) from five different regions (each ∼375 × 50 µm in area, 3,000 pixels in total) in Phantom 1, indicated by the green boxes in Figs. 3(a)–3(d). The imaging metrics corresponding to the regions are summarized in Table 3. Overall, in the regions of high OCT SNR (R1-R4), the measured strain sensitivity is similar (<5% difference), whereas in the region of low OCT SNR (R5), the strain sensitivity of the WLS-based methods is slightly superior (∼12% improvement) to that of the vector method. The improved strain sensitivity of WLS is likely due to the advantage of weighted averaging over unweighted averaging in regions of low OCT SNR. Similar to strain sensitivity, the mean strain is also comparable between each method, within 5%. As such, the measured strain SNR, in general, is close between the three different methods except in region R5 where the strain SNR of both WLS_WPU and WLS_FPU is slightly higher than that of the vector method due to the higher strain sensitivity. In addition, we reported the resolution (mean ± standard deviation) of strain micro-elastograms for all three methods, as shown in Table 3. As indicated by the dashed black lines in Figs. 3(b)–3(d), the axial resolution for each method was measured from the middle A-scan across the inclusion axial boundary, and the lateral resolution was measured from the strain step response over the inclusion lateral boundary at a depth location of ∼134 µm from the inclusion surface. The mean and standard deviation were computed from three resolution measurements. Although the texture in the strain micro-elastograms gives the appearance of different resolution, as shown in Figs. 3(b)–3(d), the measured lateral and axial strain feature resolutions for three different methods are close: 181 µm and 182 µm, respectively. Overall, both lateral and axial strain resolution for the three different methods are similar: 181 µm and 182 µm, respectively. As the same number of pixels were used to estimate strain from Phantom 1 and the same OCT system was used, it is expected that the strain resolution is similar for the three different methods.

Table 3. Strain sensitivity, mean strain, strain SNR, and strain resolution (mean ± standard deviation) for the three strain estimation methods.

View Table | View all tables in this article

4.2 Strain accuracy

In this section, we investigate the accuracy of strain estimation for the different methods over a range of strains imparted to Phantom 2, a homogeneous silicone phantom. Figure 4(a) shows the OCT image of Phantom 2. To determine strain accuracy, we require a method to measure the strain introduced to the phantom independently of the three strain estimation methods under analysis. To achieve this, we defined the expected strain using a method previously implemented in quantitative micro-elastography (QME) [43] where strain is calculated as the displacement introduced to the phantom using a PZT actuator (as described in Section 3.2) divided by the unloaded phantom thickness, measured from the OCT intensity image. In a common-path system, the displacement is calculated from the phase difference at the phantom boundary distal from the imaging plate, i.e., at the back plate, importantly, for the specific case of Phantom 2, it was possible to image through the entire phantom. This allowed us to calculate strain without using one of the strain estimation methods typically used in OCE. Also, as described below, this method provided a strain accuracy far superior to that achievable using the strain estimation methods, providing a “ground truth” strain measurement. It is important to note that this method is not generally applicable in OCE, as it is only suited to samples whose thickness can be determined directly from an OCT image, which is not typically the case in turbid samples. Another critical aspect is that the boundaries should be well lubricated to ensure uniform axial strain throughout the sample. Using this approach, phase-sensitive detection was used to measure the displacement of the phantom at the distal back plate, with the number of wrapping events manually counted and added to the phase difference measured at the interface, avoiding any inaccuracies introduced by phase unwrapping. The phantom thickness was measured from the distance between the first pixel in the OCT image (Fig. 4(a)) and a pixel of high OCT SNR (∼58 dB) at the sample-plate interface. The error in the strain measured using this method is the combination of errors in the displacement measurement at the interface and errors in the phantom thickness measurement. For the case of a specular reflection, it can be assumed that the phase difference error measured from any specular reflection (e.g., a bright interface) is inversely proportional to OCT SNR [20,21,56] and, thus, the displacement error is 3.7 nm for a corresponding OCT SNR of 58 dB. For the phantom thickness measurement, the maximum error of the measured peak intensity is determined by the axial pixel size, corresponding to ∼2.4 µm in the phantom (assuming a bulk refractive index of 1.4). Applying the law of error propagation, the error of expected strain is ∼7.2 µɛ, which is more than an order of magnitude higher strain sensitivity than that reported in Table 3.

 

Fig. 4. Accuracy of strain estimation using the vector method, WLS_WPU and WLS_FPU, respectively. (a) OCT B-scan of Phantom 2, where a spatial window (marked by a white box) with a size of ∼1000 × 250 µm is used to quantify strain sensitivity (S_mɛ), strain accuracy (µ_mɛ) and strain SNR (SNR_mɛ). Magnified views of the window in the strain micro-elastograms at 2.5 mɛ obtained using (b) the vector method, (c) WLS_WPU and (d) WLS_FPU, respectively. Scale bars in (b)-(d) represent 100 µm. The black arrow in (b) highlights a line artifact corresponding to positive strain in the case of the vector method. (e) S_mɛ, (f) µ_mɛ and (g) SNR_mɛ as a function of expected strain. The black line in (f) represents the ideal strain measurements where OCE strain matches expected strain.

Download Full Size | PDF

S_mɛ, µ_mɛ and SNR_mɛ at different expected strain were measured from the ∼1000 × 250 µm white box indicated in Fig. 4(a) in the strain micro-elastograms of Phantom 2 obtained using the three different methods. Figures 4(b)–4(d) show the strain micro-elastograms at 2.5 mɛ for the vector method, WLS_WPU and WLS_FPU, respectively. Figure 4(e) shows the measured strain sensitivity of the different strain estimation methods plotted as a function of expected strain measured. Within the low strain regime (i.e., <1.1 mɛ), the strain sensitivity of all three methods is similar, but degrades slowly as the expected strain increases. The degradation of strain sensitivity is attributed to strain-induced phase decorrelation, which invalidates the assumption that the unloaded and loaded scans are fully correlated [31]. Within the strain regime 1.4-2.5 mɛ, the strain sensitivity of all three methods degrades more rapidly, particularly for the vector method. This is likely due to the positive strain artifacts (i.e., vertical bright lines indicated by a black arrow in Fig. 4(b)), which start to appear in the strain micro-elastogram obtained using the vector method at ∼1.4 mɛ. Compared to WLS-based methods, these positive strain artifacts occur at lower expected strain and introduce spatial variation to the strain micro-elastogram of a mechanically homogeneous phantom, decreasing the precision (i.e., strain sensitivity) of strain micro-elastograms generated using the vector method. The formation of this artifact results from the accumulated errors in the unweighted vector summation over a set of normalized phasor gradients (i.e., only phase angles retained) including some noisy phase signals, which are likely to be inaccurate in regions of dark speckle [21]. As both strong and weak signals are normalized to unity, these weak and inaccurate phase values reduce the accuracy of the averaged phase gradient after vector summation.

Figure 4(f) shows the mean strain estimated using the different methods plotted against the expected strain, overlaid with a solid black line that indicates the ideal result where expected strain and estimated strain are equal. Up to the expected strain of 1.4 mɛ, the strain estimated using all three methods is similar and shows good agreement with the ground truth. Above 1.4 mɛ, the vector method significantly underestimates the expected strain. For instance, at 2.1 mɛ, the error of strain estimation using the vector method is approximately 5 and 7 times larger than WLS_WPU and WLS_FPU, respectively. Analogously to the degradation of strain sensitivity using the vector method, the accuracy is also reduced due to a significant underestimation of expected strain resulting from the formation of positive strain line artifacts. Above 2.1 mɛ, WLS_WPU also significantly underestimates the strain. At 2.5 mɛ, the vector method and WLS_WPU deviate by ∼42% and ∼60% from the expected value, respectively. In contrast, at this point, the error of strain estimation using WLS_FPU is ∼4%, which is nearly 11 and 15 times smaller than WLS_WPU and the vector method, respectively. The underestimation of strain using WLS_WPU is attributed to increased strain-induced decorrelation resulting in decreased phase difference sensitivity, causing errors in the calculation of the axial phase unwrapping reference in WPU. We observed that, as strain increases above 2.5 mɛ, a high spatial frequency noise due to phase decorrelation occurs at local depth locations, resulting in zero strain (i.e., a flattened unwrapped phase difference-depth profile) in these locations. This suggests that the unwrapped phase difference using WPU no longer reflects the true displacement in the sample, as indicated in Eq. (3). This ultimately leads to the reduction of strain accuracy using WLS_WPU, setting an upper limit to its dynamic range of accurate strain estimation. Compared to WPU, by performing Volkov symmetrization in 3-D space, FPU incorporates more data points (i.e., eight duplicates of the whole data volume), providing a more accurate estimation of unwrapped phase difference in the presence of noise. In addition, FPU is implemented in an iterative manner in which phase unwrapping errors are further reduced by iterating Q(r) (i.e., the phase unwrapping error function) until its spatial variance meets δ (i.e., the threshold accuracy), leading to a larger dynamic range of accurate displacement estimation and, by extension, strain estimated by WLS_FPU, in comparison to the other two methods. More importantly, as highlighted in Fig. 4(f), the extended strain dynamic range of WLS_FPU could potentially lead to improved strain contrast. For instance, for expected strains of 1.4 mɛ and 2.5 mɛ, the strain contrast (i.e., strain ratio) between these two points measured by WLS_FPU is ×2 larger than that of the vector method, suggesting the suitability of WLS_FPU for imaging strain in mechanically heterogeneous samples. Figure 4(g) presents the measured strain SNR for the three different methods, with WLS_FPU providing the highest strain SNR when the expected strain is above 2.1 mɛ.

4.3 Comparison of the performance of strain estimation methods in biological tissue

To compare the performance of the three strain estimation techniques on biological tissue, we first imaged a mouse muscle, as such tissue has previously been shown to provide large variations in strain across muscle fibers [57]. Figure 5 presents the results of the encapsulated EDL muscle using the different strain estimation methods. An en face OCT image of the tissue at a depth of ∼340 µm below the tissue surface is presented in Fig. 5(a). Figures 5(b)–5(d) demonstrate the corresponding strain micro-elastograms. To better visualize the difference in image quality between the different methods, we selected a local region (∼3 × 1 mm) that contains a characteristic feature showing a muscle striation close to the highly scattering tendon, which is highlighted by a white arrow in the magnified OCT image shown in Fig. 5(e). The corresponding magnified strain micro-elastograms are shown in Figs. 5(f)–5(h). Similar to a previous study [57], striation patterns manifest in both the OCT image and the strain micro-elastograms. Noticeably, the striation is much more clearly delineated by WLS_FPU in Fig. 5(h), in comparison to the other two methods, especially for the case of the vector method where the striation is hardly visible (Fig. 5(f)). The improved contrast of muscle striation provided by WLS_FPU, corresponding to larger compressive strains in the micro-elastogram of WLS_FPU (Fig. 5(h)), likely results from the extended dynamic range of WLS_FPU (Fig. 4(f)).

 

Fig. 5. Excised mouse EDL muscle encapsulated in the hydrogel. En face (a) OCT image at a depth of ∼340 µm below the tissue surface. The corresponding strain micro-elastograms obtained using (b) the vector method, (c) WLS_WPU and (d) WLS_FPU, respectively. Magnified regions-of-interest (marked by dashed blue boxes in (a)-(d)) for (e) OCT, (f) the vector method, (g) WLS_WPU and (h) WLS_FPU. White arrows indicate a location where additional strain contrast of a muscle striation can be observed in (h) using WLS_FPU. En face unwrapped phase difference images obtained using (i) WPU and (j) FPU, respectively. Unwrapped phase difference B-scans obtained using (k) WPU and (l) FPU, respectively. Cross-sections in the B-scan plane and en face plane are represented by dashed green and red boxes. The dashed black arrows in (i) and (k) indicate a streak-like phase unwrapping artifact caused by WPU, corresponding to an artifact in the strain micro-elastogram generated using WLS_WPU in (c) and (g).

Download Full Size | PDF

Figures 5(i) and 5(j) show the en face unwrapped phase difference images obtained using the WPU and FPU algorithms, respectively. Interestingly, distinctive imaging artifacts (highlighted by a dashed black arrow in Fig. 5(i)) that propagate from the tissue boundary can be observed in the unwrapped phase difference image generated using WPU, which result in streak-like imaging artifacts appearing in the corresponding strain micro-elastogram obtained with WLS_WPU (highlighted by a dashed black arrow in Fig. 5(c)). As the tissue is encapsulated within clear hydrogel with correspondingly low OCT SNR, the phase difference sensitivity in the hydrogel is low. As such, the axial phase unwrapping reference calculated from the preceding n_z (i.e., n_z = 7 herein) axial pixels is inaccurate. As highlighted in the unwrapped phase difference B-scan obtained using WPU in Fig. 5(k), this error is then propagated in the axial phase unwrapping of the subsequent axial pixels, followed by the lateral phase unwrapping [15]. In contrast, owing to the advantage of using an enlarged data volume through Volkov symmetrization and 150 iterations for the estimation of unwrapped phase difference, such an artifact is not present in the unwrapped phase difference images in the same plane generated using FPU (Figs. 5(j) and 5(l)) and the corresponding strain micro-elastograms (Figs. 5(d) and 5(h)). Similarly, since the vector method avoids the need for phase unwrapping, the strain micro-elastogram is free of this imaging artifact in Figs. 5(b) and 5(f), as expected.

In addition, we analyzed the performance of each strain estimation method on excised human breast tissue. Figure 6(a) shows an en face OCT image at a depth of 100 µm from the tissue surface. Analogously to Fig. 5, we selected a ∼4 × 4 mm region-of-interest (highlighted by a dashed blue box in Fig. 6(a)) in OCT, for better visualization of differences in micro-elastogram quality, and the magnified OCT image of this region is presented in Fig. 6(b). The co-registered histology and the corresponding strain micro-elastograms of the vector method, WLS_WPU and WLS_FPU are presented in Fig. 6(c)–6(f), respectively. In contrast to the OCT image in Fig. 6(b) where the optical contrast between the tumor and benign tissue is limited, the tumor is clearly delineated by distinctive mechanical contrast provided by all three strain estimation methods in Figs. 6(d)–6(f), showing a close correspondence with the co-registered histology in Fig. 6(c). The strain estimated by each method presents a heterogeneous pattern in regions corresponding to tumor, which is consistent with previous OCE studies [7–9]. Whilst the improvement in contrast provided by WLS_FPU is lower here than the case of the mouse muscle shown in Fig. 5, additional contrast in the tumor region is still visible in the WLS_FPU strain micro-elastogram (Fig. 6(f)), highlighted by white arrows, particularly with respect to the vector method (Fig. 6(d)). The relatively low improvement in contrast observed in Fig. 6(f) is likely attributed to lower strain induced in this much thicker breast tissue (×7 thicker than the mouse muscle demonstrated in Fig. 5) under compression. As such, the advantage of larger dynamic range of accurate strain estimation using WLS_FPU is less pronounced. Similar to Fig. 5(c), in Fig. 6(e), streak-like artifacts in the WLS_WPU strain micro-elastogram can be seen, which are attributed to phase unwrapping errors in WPU.

 

Fig. 6. Excised human breast tissue containing invasive tumor. (a) En face OCT image and (b) the magnified region-of-interest marked by the dashed blue box in (a) (∼4 × 4 mm in area). (c) The corresponding co-registered histology. L: Lobules, T: Tumor. The corresponding zoomed-in strain micro-elastograms for (d) the vector method, (e) WLS_WPU and (f) WLS_FPU. White arrows showing additional strain contrast provided by WLS_FPU.

Download Full Size | PDF

5. Discussion

In this paper, we have presented the first framework that provides a comprehensive comparison between WLS and the vector method: the two most prevalent strain estimation methods. In addition, we proposed a new implementation of WLS, i.e., WLS_FPU, by implementing FPU in OCE for the first time. We demonstrated that WLS_FPU provides comparable strain sensitivity and resolution to the existing methods, but an increased dynamic range of accurate strain estimation with ∼4% error that is approximately 11 and 15 times lower than WLS_WPU and the vector method, respectively, for an expected strain of 2.5 mɛ. Through a qualitative comparison of tissue results, we have demonstrated that the improved dynamic range of accurate strain translates to a noticeable improvement in strain contrast of tissue provided by WLS_FPU. In addition, WLS_FPU removes phase unwrapping errors in WLS_WPU that result in imaging artifacts in strain micro-elastograms.

In Fig. 3, we demonstrated that a highly scattering, stiff feature with dimensions below the system resolution of the OCE system presents differently in the WLS-based methods and the vector method, respectively. In the example shown in Fig. 3, whilst the WLS fitting range is ×3 larger than the size of the bright feature, the corresponding strain of the feature is detected using WLS, demonstrating the detection of sub-resolution mechanical features in OCE for the first time. However, as the feature size is below the system resolution, the strain in the soft region adjacent to the feature is significantly underestimated by WLS (Figs. 3(g) and 3(h)). Furthermore, it can be inferred that, in the case that a soft or stiff feature is less scattering than the surrounding region, the weighting used in WLS may result in a feature being masked by the surrounding region. Whilst this effect has been demonstrated here to provide useful contrast, it should be highlighted that these scenarios are examples where the assumption used in WLS that a more accurate fitting is enabled by accounting for OCT SNR is violated. This coupling between OCT SNR and mechanical contrast suggest that it is prudent to analyze the OCT image before deciding on the optimal strain estimation technique to use for a given sample.

In Fig. 4(f), we demonstrated the effect of applied strain on the accuracy of strain estimation for three different methods. Compared to the WLS-based methods, in particular to WLS_FPU, the results suggest that the vector method is less tolerant to strain-induced decorrelation noise, where the “line” artifacts (Fig. 4(b)) corresponding to positive strain becoming dominant in the micro-elastogram at an expected strain of 2.1 mɛ, resulting in an underestimation of strain by ∼36%. As detailed in Section 4.2, the formation of this artifact is attributed to the accumulated errors of strain estimation arising from the weak and inaccurate phasor gradients within the axial range of unweighted vector summation. In addition, as the phasor gradient in the vector method is calculated from the finite difference of two adjacent phasor difference in depth, it is subject to errors caused by optical noise. As such, simply increasing axial range of the unweighted vector summation worsens the strain image quality. One solution is to assign a weight of OCT SNR or, ideally, speckle brightness [21] to the phasor gradient at each pixel over the axial range of vector summation, which should improve the sensitivity and accuracy of strain estimation using the vector method. Besides, it is straightforward to increase the size of the pre-filter window of the weighted averaging for improved strain sensitivity and strain accuracy at the expense of strain resolution. For the calculation of the phasor gradient, one could calculate the finite difference between two discrete depth locations at a step of n (n ≥ 2) pixels. This essentially improves the SNR of the phasor gradient since it increases n-fold the angle of the phasor gradient (i.e., the phase gradient) with a constant phase noise at each pixel. However, as it skips several pixels in depth, phase wrapping may occur within the step when the applied strain imparted to a sample is very high, making it difficult to choose an appropriate step size for the optimal strain SNR in mechanically heterogeneous biological tissues where phase wrapping could vary significantly. In addition, we note that both strain sensitivity and strain accuracy should be taken into account in OCE image interpretation. Whilst the mean strain measured in a local region could be close to the true value, the low strain sensitivity corresponding to larger spatial variation of strain measurements makes it challenging to differentiate a real mechanical feature from the noise. Thus, strain SNR (the ratio of mean strain to strain sensitivity), as illustrated in Table 2, can be used to assess the reliability of strain micro-elastograms. Separately, to trace the origin of the strain imaging artifact at various reconstruction steps, the images of OCT intensity, wrapped phase difference, unwrapped phase difference (if needed), and strain can be displayed at the same plane, and the artifacts can be identified from the comparison between these images, as shown in Fig. 5.

In Sections 4.1 and 4.2, we demonstrated the numerical comparisons of strain imaging metrics from one inclusion phantom and one homogeneous phantom. In general, to demonstrate the repeatability and reproducibility of each method, the experimental results of both intra-sample comparison and inter-sample comparison should be performed. However, previous studies suggested that, in addition to the contribution of strain estimation that we specifically focus on here, the strain micro-elastogram quality is determined by mechanical deformation [31–34]. As such, intra-sample variations (e.g., inconsistent frictional conditions [32,52,55]) and inter-sample variations (e.g., sample geometries and stiffness [33,34]) would obscure variations between strain estimation methods that we are targeting in this paper. A robust compression OCE setup and a framework to enable meaningful inter-sample comparison of strain imaging metrics are worth investigating in the future.

In this paper, we demonstrated that incorporating a robust phase unwrapping algorithm into WLS can provide improved strain dynamic range over the vector method, which leads to additional imaging contrast, as demonstrated by the results in Section 4.3. However, the necessity of phase unwrapping and the varying weights in WLS make it challenging to implement this method rapidly. A main advantage of the vector method is its computational efficiency as no phase unwrapping is required. This means that the computation can readily be accelerated using complex convolution rather than vector summation and the computation of phasor difference and phasor gradients can be parallelized. This suggests that, despite the improved dynamic range of accurate strain achievable using WLS_FPU, in applications requiring real-time visualization of strain, e.g., imaging strain dynamics induced by thermal heating [58–62], it may be preferable to implement the vector method. This highlights the need to take a considered approach in deciding which strain estimation method to use in a given application. We believe the analysis performed in this study can assist in making that decision.

6. Conclusion

We have presented a framework for analyzing strain estimation in phase-sensitive compression OCE. We have demonstrated the use of this framework by proposing a new strain estimation method, WLS_FPU, and more importantly, we have presented a quantitative comparison between this new method and the two most prevalent methods in compression OCE, i.e., WLS_WPU and the vector method. Our experimental results indicate that WLS_FPU provides an increased dynamic range of accurate strain estimation by reducing errors to ∼4% that is approximately 11 and 15 times smaller relative to WLS_WPU and the vector method, respectively, for an expected strain of 2.5 mɛ. Through a qualitative comparison of tissue results, we have shown that this increased dynamic range translates to higher mechanical contrast in strain micro-elastograms of heterogeneous tissue. We have also shown that imaging artifacts caused by phase unwrapping in WLS_WPU can be avoided using WLS_FPU. Furthermore, we have demonstrated the detection of sub-resolution mechanical features using WLS-based strain estimation. We believe that our framework provides a useful mechanism to compare strain estimation methods that could be adapted to compare different signal processing methods used in other OCE techniques such as transient OCE.