Strain estimation in phase-sensitive optical coherence elastography

Brendan F. Kennedy; Sze Howe Koh; Robert A. McLaughlin; Kelsey M. Kennedy; Peter R. T. Munro; David D. Sampson

doi:10.1364/BOE.3.001865

1. Introduction

Disease often causes alteration of the elastic properties of tissue [1,2]. The emerging technique of optical coherence elastography (OCE) is aimed at differentiating tissues by imaging their microscopic elastic properties [3–6]. In OCE, a tissue sample is subjected to a load (stress) and the resulting local displacement is measured using optical coherence tomography (OCT) [7]. The spatial derivative of the displacement is the strain [1,2], which is indicative of the sample’s elastic properties. The strain is displayed in a two- or three-dimensional image, known as an elastogram [2]. The limits on resolution and penetration depth of displacement measurements in OCE are set by OCT, the underlying imaging modality and, thus, lie in the range 1-20 μm and 1-3 mm, respectively.

Early work in OCE used speckle-tracking algorithms to estimate displacement [3–6]. Speckle tracking is complex and, thus, difficult to implement in real time, and has restrictive upper and lower limits on displacement, set by speckle decorrelation [4,8–10] and spatial sampling in the underlying OCT system, respectively. In phase-sensitive OCE, proposed in part to address these issues [11–16], the axial component of displacement is calculated from the phase difference between A-scans originating from within a single B-scan [11] or from successive B-scans [12]. The accuracy and precision of the measured displacement are limited by the signal-to-noise ratio (SNR) of the OCT signal [17], and can be corrupted by phase wrapping [18]. Since the strain is a differential measure, any errors are amplified in its calculation [19]. To date, there has been no published framework that relates errors in the displacement measurement to the calculated strain, making it difficult to optimize OCE, assess elastogram image quality and quantitatively compare OCE techniques.

In this paper, we present a theoretical framework for strain imaging parameters in phase-sensitive OCE based on a statistical analysis of the measured displacement. We focus on two previously proposed strain estimation methods: finite difference (FD) [3] and ordinary least squares (OLS) [20,21]. We define expressions for strain sensitivity, SNR and dynamic range. In this analysis, we treat each measurement as a realization of an underlying random variable. We compare the performance of both methods using theoretical predictions and experimental measurements, and demonstrate the significant sensitivity advantage of OLS strain estimation. We propose a new modified least-squares strain estimation method for OCE, based on weighted least squares (WLS). In this method, a weight is assigned to each displacement measurement, based on the known OCT SNR at that location. We present elastograms generated using each method and demonstrate an improvement in the strain sensitivity of WLS strain estimation compared to both FD and OLS strain estimation. Finally, we also demonstrate an improvement in strain imaging parameters using Gaussian smoothing.

The paper is structured as follows: Section 2 describes the measurement of strain with OCE and provides a statistical analysis of both FD and OLS strain estimation. In Section 3, this framework is used to define strain SNR, sensitivity and dynamic range, and the new WLS method for strain estimation is introduced. Section 4 provides details of the experimental setup and methods used to validate the theoretical framework. Section 5 presents elastograms recorded from phantoms and freshly excised porcine airway. Section 6 presents discussion of the results and, finally, Section 7 presents the main conclusions of the work.

2. Phase-sensitive OCE strain

Young’s modulus, Y, is often used to characterize a material’s elastic properties. It is defined as the ratio of stress, the force per unit area, to strain, the resulting fractional change in length of the sample, and is given by [1]

Y = \frac{σ}{ε_{b}} = \frac{\frac{F}{A}}{\frac{Δ l}{l_{0}}},

where σ is the stress,

ε_{b}

is the bulk strain, F is the applied force, A is the surface area, $Δ l$ is the change in length, and

l_{0}

is the original length of the sample. The objective of OCE is to estimate a sample’s depth-resolved elastic properties from OCT-derived displacement data. The strain can be used as a surrogate for elasticity, if the stress introduced throughout the sample is assumed to be uniform [1,2]. However,

ε_{b}

, defined in Eq. (1), is not depth-resolved; it refers to a single measurement of the bulk strain present in a sample. In OCE, the local strain,

ε_{l}

, is a more suitable definition [1–3,20]:

ε_{l} = \frac{Δ d}{Δ z},

where

Δ d

is the change in displacement measured over an axial depth range

Δ z

, which defines the axial resolution of the local strain measurement. Assuming uni-axial loading, if

Δ z

is the same as the total length of the sample,

l_{0}

, then

ε_{l}

and

ε_{b}

are equivalent. In general, strain measured in elastography refers to local rather than bulk strain [1,2]. We will adopt this generalization and, therefore, refer to local strain as simply ε.

The goal of strain estimation in OCE is to obtain the strain defined in Eq. (2) from a set of measured displacements. Two major factors in determining the strain measurement accuracy are the accuracy of the displacement measured from the OCT data and the accuracy of the method used to calculate strain from the measured displacement.

2.1. Displacement measurement

A single complex OCT A-scan containing N data points may be represented by the set

I^{(j)} = {(z_{i}, a_{i}^{(j)} \exp (ι φ_{i}^{(j)})) | z_{i} = (i - 1) δ z + z_{1} \in R, a_{i}^{(j)} \in R^{+}, φ_{i}^{(j)} \in [- π, π)}_{i = 1}^{N},

where 𝚤 is the imaginary unit,

z_{i}

are the depths to which the complex OCT measurements

a_{i}^{(j)} \exp (ι φ_{i}^{(j)})

correspond,

δ z

is the separation between successive

z_{i}

, and

z_{1}

is the first depth in the A-scan. In this paper, the superscript j is used to denote a particular A-scan and could be, for example, one of many making up a B-scan, or two A-scans acquired in the same lateral position in the sample in successive B-scans. The phase,

φ_{i}^{(j)}

, is in general random, although a map of the displacement of scattering centers, occurring between the acquisition of A-scans

I^{(2)}

and

I^{(1)}

can be found, without the loss of generality, from the phase difference between complex OCT signals as

D = {(z_{i}, d_{i}) | d_{i} = \frac{λ (ϕ_{i}^{(2)} - ϕ_{i}^{(1)})}{4 π n} = \frac{λ Δ ϕ_{i}}{4 π n}}_{i = 1}^{N},

where λ is the mean source wavelength and n is the refractive index of the material. Displacements deriving from a phase change of greater than

\pm π

are ambiguous. This is a limitation referred to as phase wrapping [18] and sets the upper limit of measurable phase difference, without the need for complicated phase unwrapping algorithms.

The minimum measurable phase difference is limited by the standard deviation of the random variable for which $Δ φ_{i}$ is a realization [17]. Under the assumption that the SNR of the detected OCT signal intensity $S N R_{O C T_{i}}$ satisfies, $S N R_{O C T_{i}}$ $> > 1$ , it can be shown that the phase sensitivity, $σ_{Δ φ_{i}}$ , is given by [17]:

σ_{Δ φ_{i}} = \frac{1}{\sqrt{S N R_{O C T_{i}}}} .

For example, considering a typical

S N R_{O C T}

of 20 dB, combining Eqs. (4) and (5) gives a displacement sensitivity of 4.7 nm, assuming λ = 800 nm. It should be noted that if _{$S N R_{O C T_{i}} > > 1$}, the phase difference is unbiased; thus, the expected value of the displacement random variable corresponds to the true sample displacement [22].

2.2. Strain estimation

In OCE, two methods have been proposed to estimate strain from the measured displacement map: FD strain estimation [3] and OLS strain estimation [20,21]. To determine the accuracy and precision of these methods, we define expressions for their mean and variance.

2.2.1. Finite difference strain estimation

Strain estimation using the FD method was previously proposed in the related field of ultrasound elastography [23], and later applied to OCE [3]. FD strain estimation, $ε_{i}^{f}$ , is defined as

ε_{i}^{f} = \frac{d_{i + m - 1} - d_{i}}{z_{i + m - 1} - z_{i}} = \frac{Δ d}{Δ z} .

where

Δ z = m δ z

. As the strain is obtained from the linear combination of two random variables, its variance can be defined as [24]:

σ_{Ε_{i}^{f}}^{2} = \frac{σ_{D_{i}}^{2} + σ_{D_{i + m - 1}}^{2} - 2 c o v (D_{i}, D_{i + m - 1})}{Δ z^{2}},

where

Ε_{i}^{f}

,

D_{i}

and

D_{i + m - 1}

are random variables of which

ε_{i}^{f}

,

d_{i}

and

d_{i + m - 1}

are realizations. Assuming the measured displacements are from non-overlapping positions, they can be assumed to be statistically independent [24] and the covariance term in Eq. (7) can be set to zero. The expression is further simplified assuming the displacement measurements are homoscedastic, i.e., the variance is the same for

D_{i}

and

D_{i + m - 1}

, and is written simply as

σ_{D}^{2}

:

σ_{Ε_{i}^{f}}^{2} = \frac{2 σ_{D}^{2}}{Δ z^{2}} .

From Section 2.1, the displacement estimation is unbiased. Therefore, FD strain estimation is also unbiased; the mean of the random variable

Ε_{i}^{f}

corresponds to the true underlying strain:

μ_{Ε_{i}^{f}} = ε .

2.2.2. Ordinary least squares strain estimation

An alternative method to estimate strain is ordinary least squares (OLS). In contrast to the FD method, in the OLS method use is made of more than two data points to obtain greater accuracy. In ultrasound elastography, the OLS method has been shown to provide more accurate strain estimation than the FD method [19]. We aim here to provide such a comparison in OCE for the first time. By restricting z to the interval $z_{i} \leq z \leq z_{i} + Δ z$ , we can rewrite Eq. (2), implicitly defining ε:

d = ε (z - z_{i - 1}) + c .

Given a set of depth values and associated displacement measurements, d, within the interval defined above, the strain is determined by minimizing the sum of squared differences, R, between the measured displacements and the displacements given by Eq. (10), i.e., by minimizing the following expression:

R = {\sum_{j = i}^{i + m - 1} [d_{j} - ε (z_{j} - z_{i - 1}) - c]}^{2} .

A complete derivation of OLS is provided elsewhere [25]. Following this derivation, the strain can be estimated using the OLS method according to

ε_{i}^{o} = \frac{(\sum_{j = i}^{i + m - 1} 1) (\sum_{j = i}^{i + m - 1} (z_{j} - z_{i - 1}) d_{j}) - (\sum_{j = i}^{i + m - 1} (z_{j} - z_{i - 1})) (\sum_{j = i}^{i + m - 1} d_{j})}{(\sum_{j = i}^{i + m - 1} 1) (\sum_{j = i}^{i + m - 1} {(z_{j} - z_{i - 1})}^{2}) - {(\sum_{j = i}^{i + m - 1} (z_{j} - z_{i - 1}))}^{2}},

where m is the number of data points used in the strain calculation. Using this expression, we obtain an estimation of strain at each data point, i, within an A-scan. Equation (12) can be simplified by reducing the number of summation operators such that

ε_{i}^{o} = \sum_{j = i}^{i + m - 1} (\frac{κ_{0} (z_{j} - z_{i - 1}) - κ_{1}}{κ_{0} κ_{2} - κ_{1}^{2}}) d_{j},

where

κ_{x} = \sum_{j = i}^{i + m - 1} {(z_{j} - z_{i - 1})}^{x}, x = 0, 1, 2.

Substituting

z_{j} - z_{i - 1} = (j - i + 1) δ z

in Eqs. (13) and (14) allows

ε_{i}^{o}

to be expressed as

ε_{i}^{o} = \sum_{j = i}^{i + m - 1} (\frac{κ_{0} (j - i + 1) δ z - κ_{1}}{κ_{0} κ_{2} - κ_{1}^{2}}) d_{j} .

Equation (14) can now be simplified using mathematical induction to remove the summation operators:

κ_{0} = \sum_{j = i}^{i + m - 1} {((j - i + 1) δ z)}^{0} = m,

κ_{1} = \sum_{j = i}^{i + m - 1} {((j - i + 1) δ z)}^{1} = \frac{m (m + 1)}{2} δ z,

κ_{2} = \sum_{j = i}^{i + m - 1} {((j - i + 1) δ z)}^{2} = \frac{m (m + 1) (2 m + 1)}{6} δ z^{2} .

Substituting Eqs. (16)-(18) in Eq. (15) and simplifying, we obtain

ε_{i}^{o} = \sum_{j = i}^{i + m - 1} (\frac{6 (2 (j - i) - m + 1)}{δ z (m^{3} - m)}) d_{j} .

To allow expressions for the mean and variance of

ε_{i}^{o}

to be derived, we assume that the set of displacements has the same variance, i.e., homoscedasticity, and use the scaling and summing properties of independent Gaussian random variables. The mean of

Ε_{i}^{o}

, the random variable of which

ε_{i}^{o}

is a realization, can be expressed as

μ_{Ε_{i}^{o}} = \sum_{j = i}^{i + m - 1} (\frac{6 (2 (j - i) - m + 1)}{δ z (m^{3} - m)}) μ_{D_{j}} .

As the random variable

D_{j}

is unbiased (see Section 2.1), its expected value,

μ_{D_{j}}

, equals the true displacement. Replacing

z - z_{i - 1}

with

(j - i + 1) δ z

in Eq. (10), and then substituting this into Eq. (20), and simplifying using mathematical induction, we obtain

μ_{Ε_{i}^{o}} = \sum_{j = i}^{i + m - 1} (\frac{6 (2 (j - i) - m + 1)}{δ z (m^{3} - m)}) (ε (j - i + 1) δ z + c) = ε .

Equation (21) confirms that the mean of

Ε_{i}^{o}

equals the true underlying strain. Similarly, the variance of

Ε_{i}^{o}

may be defined as

σ_{Ε_{i}^{o}}^{2} = {\sum_{j = i}^{i + m - 1} (\frac{6 (2 (j - i) - m + 1)}{δ z (m^{3} - m)})}^{2} σ_{D}^{2} = \frac{12 σ_{D}^{2}}{δ z^{2} (m^{3} - m)} \approx \frac{12 σ_{D}^{2}}{Δ z^{2} m} .

3. Strain imaging parameters

Having derived the mean and variance of the strain estimated using both FD and OLS methods, we now define the strain sensitivity, SNR and dynamic range in terms of these parameters. The strain sensitivity is the minimum detectable strain in an elastogram. Analogously to the definition of phase difference sensitivity in Eq. (5), the strain sensitivity for FD and OLS strain estimation, _{$S_{Ε_{i}^{f / o}}$}, is defined as the standard deviation of the strain [26]:

S_{Ε_{i}^{f / o}} = σ_{Ε_{i}^{f / o}} .

In ultrasound elastography, the strain SNR,

S N R_{Ε_{i}^{f / o}}

, is defined as the ratio of the mean to the standard deviation of the strain [26]. We choose the same definition for strain SNR in OCE:

S N R_{Ε_{i}^{f / o}} = \frac{μ_{Ε_{i}^{f / o}}}{σ_{Ε_{i}^{f / o}}} .

Strain dynamic range,

D R_{Ε_{i}^{f / o}}

, is the ratio of the largest to the smallest measurable strain [26]. Maximum strain in OCE occurs when the maximum measureable change in displacement, Δd_max, occurs over one strain axial resolution element. To avoid the ambiguity caused by phase wrapping, the upper limit of displacement considered in this paper corresponds to a phase difference of π radians. Minimum strain is defined by _{$σ_{Ε_{i}^{f / o}}$}, allowing the dynamic range to be defined as

D R_{Ε_{i}^{f / o}} = \frac{(\frac{Δ d_{max}}{Δ z})}{σ_{Ε_{i}^{f / o}}} = \frac{λ}{4 n σ_{Ε_{i}^{f / o}} Δ z} .

3.1. Strain imaging parameters using the FD method

Substituting the mean and variance of the strain estimated using the FD method, defined in Section 2.2.1, into the definitions of the strain imaging parameters describes the performance of FD strain estimation. The FD strain sensitivity, $S_{Ε_{i}^{f}}$ _, SNR, $S N R_{Ε_{i}^{f}}$ , and dynamic range, $D R_{Ε_{i}^{f}}$ , are found to be

S_{Ε_{i}^{f}} = \frac{\sqrt{2} σ_{D}}{Δ z};

S N R_{Ε_{i}^{f}} = \frac{μ_{Ε_{i}^{f}} Δ z}{\sqrt{2} σ_{D}}; and

D R_{Ε_{i}^{f}} = \frac{λ}{\sqrt{32} σ_{D} n} .

3.2. Strain imaging parameters using the OLS method

Similarly, substituting the mean and variance of the strain estimated using OLS, defined in Section 2.2.2, in the expressions for the strain imaging parameters defined above allows us to assess the performance of this method. Using this approach, the OLS strain sensitivity, $S_{Ε_{i}^{o}}$ , SNR, $S N R_{Ε_{i}^{o}}$ , dynamic range, $D R_{Ε_{i}^{o}}$ , are defined as

S_{Ε_{i}^{o}} = \frac{\sqrt{12} σ_{D}}{Δ z \sqrt{m}};

S N R_{Ε_{i}^{o}} = \frac{μ_{Ε_{i}^{o}} Δ z \sqrt{m}}{\sqrt{12} σ_{D}}; and

D R_{Ε_{i}^{o}} = \frac{λ \sqrt{m}}{\sqrt{192} σ_{D} n} .

3.3. Theoretical predictions

In Fig. 1 , we present theoretical predictions obtained using the expressions derived in Sections 3.1 and 3.2 for FD and OLS strain estimation. The OCT SNR and source mean wavelength were assumed to be 20 dB and 835 nm, respectively. The variation in strain sensitivity with strain axial resolution, presented in Fig. 1(a), demonstrates the increasing sensitivity advantage of OLS over FD strain estimation as the resolution degrades. An improvement in sensitivity of ~8 dB is achieved using OLS estimation at the lowest simulated axial resolution of Δz = 120 μm.

Fig. 1 Theoretical predictions of strain estimation using FD (dashed red) and OLS (solid blue) methods: (a) sensitivity; and (b) SNR versus strain axial resolution. (c) SNR versus change in displacement, Δd, in one resolution element.

Download Full Size | PDF

In Fig. 1(b), the SNRs of the strain using FD and OLS estimation are compared. The strain was chosen to be $1200 \times 10^{- 6}$ , corresponding to the maximum change in displacement, Δd_max, over the lowest strain axial resolution considered, i.e., 120 μm. As discussed in Sections 3.1 and 3.2, the expected value of strain using both methods is unbiased and, therefore, corresponds to the true strain. As the SNR is the ratio of mean strain to strain sensitivity, the difference between the FD and OLS methods is determined by the strain sensitivity. This is confirmed in Fig. 1(b), where the SNR of the OLS estimate is also improved by ~8 dB for Δz = 120 μm. Figures 1(a) and 1(b) highlight a fundamental trade-off in strain estimation in OCE, namely, improved strain sensitivity and SNR, using both methods, at the expense of axial resolution.

In Fig. 1(c), the SNR of FD and OLS strain estimations versus the change in displacement, Δd, measured in one strain axial resolution element, is presented. The resolution was fixed at 120 μm, corresponding to the lowest value shown in Figs. 1(a) and 1(b). For each Δd, the SNR of the OLS estimate is higher than that of the FD estimate. Figure 1(c) highlights the importance of maximising Δd in each strain resolution element.

3.4. Weighted least squares strain estimation

According to the Gauss-Markov theorem [25], OLS provides the best linear unbiased strain estimation if the errors have the same variance (i.e., assuming homoscedasticity). However, in phase-sensitive OCE, the variance of each displacement measurement is dependent on the OCT SNR (Eq. (5)). To address this and to improve the strain estimation, we adopt a weighted least-squares (WLS) method and assign a weight to each displacement measurement based on the known variance of that measurement. A strain estimation method based on WLS has been proposed in ultrasound elastography [27], however, we report its first use in OCE. The derivation of WLS strain, $ε_{i}^{w}$ , follows the same procedure as for $ε_{i}^{o}$ , enabling the following expression to be defined:

ε_{i}^{w} = \frac{(\sum_{j = i}^{i + m - 1} w_{j}) (\sum_{j = i}^{i + m - 1} w_{j} (z_{j} - z_{i - 1}) d_{j}) - (\sum_{j = i}^{i + m - 1} w_{j} (z_{j} - z_{i - 1})) (\sum_{j = i}^{i + m - 1} w_{j} d_{j})}{(\sum_{j = i}^{i + m - 1} w_{j}) (\sum_{j = i}^{i + m - 1} w_{j} {(z_{j} - z_{i - 1})}^{2}) - {(\sum_{j = i}^{i + m - 1} w_{j} (z_{j} - z_{i - 1}))}^{2}},

where the weighting variable,

w_{j}

, is chosen to be

w_{j} = \frac{1}{σ_{D_{j}}^{2}} .

In Eq. (33), _{$σ_{D_{j}}^{2}$} is determined by combining Eqs. (4) and (5). Analogously to

ε_{i}^{o}

, the strain axial resolution,

Δ z

, and the change in displacement over this range, Δd, are determined by the number of data points, m, over which

ε_{i}^{w}

is calculated. As this method incorporates more information about the measured displacement, it has the potential to provide more accurate strain estimations than either the FD or OLS methods. Note that because

w_{j}

depends on

S N R_{O C T_{j}}

, which is non-uniform and not known a priori, we are unable to derive for WLS strain estimation simple analytical expressions for strain imaging parameters from Eq. (32). Instead, we demonstrate its superiority experimentally in Section 5.

4. Experimental setup and procedure

A schematic diagram of the SD-OCE system is presented in Fig. 2(a) . The system employs a superluminescent diode source with mean wavelength of 835 nm and 55.4 nm bandwidth, illuminating the sample with 10 mW of optical power. The theoretical lateral resolution is 15 μm and the sensitivity was measured to be 108 dB for an exposure time of 98 μs. The spectrometer consists of a diffraction grating with 1800 lines/mm, a lens with a focal length of 50 mm and a line-scan camera, consisting of 4096 pixels. The sensitivity roll-off was measured to be 4 dB/mm. The A-scan acquisition time was 100 μs. Lateral scanning over a 4 mm range was performed using a galvanometer mirror, driven by a sawtooth waveform. For each elastogram, 1000 A-scans were recorded in an acquisition time of 0.1 s.

Fig. 2 (a) Schematic diagram of the OCE system; RB: rigid boundary; RA: ring actuator; SM: scanning mirror; BS: beam splitter; RM: reference mirror; SLD: superluminescent diode; (b) Saw-tooth pattern of lateral scanning beam synchronized with motion of ring actuator. (c) The measured phase difference with no lateral scanning and (d) with lateral scanning.

Download Full Size | PDF

Samples were mechanically loaded using a piezoelectric ring actuator described previously [20] that performed imaging and compression from the same side of the sample. Samples were preloaded by compression against a rigid upper boundary to ensure full contact with the sample surface. A 5 Hz square wave synchronized with the lateral scanning was applied to the actuator, as illustrated in Fig. 2(b), to produce a displacement amplitude of 0.15 μm, corresponding to the upper limit set by phase wrapping. The resulting sample displacement was calculated from the phase difference between A-scans at the same lateral position acquired in successive B-scans. We note that by using this procedure the phase sensitivity is not affected by the lateral oversampling in the x-direction, as is the case when the phase difference is measured between A-scans acquired within one B-scan [17]. However, when performing 3-D OCE measurements, scanning in the y-direction degrades phase sensitivity, and oversampling in this axis is required to minimize this reduction in phase sensitivity.

In Fig. 2, the phase difference, measured using a stationary glass slide in the sample arm, is presented for each of 100 B-scans, both without, Fig. 2(c), and with, Fig. 2(d), lateral scanning. In both cases, the OCT SNR was measured to be 50 dB. In the non-scanning case, $σ_{Δ φ}$ was measured to be 4 mrad, within 1 mrad of the value given by Eq. (5), and corresponding to a displacement sensitivity of 0.2 nm. In the scanning case, $σ_{Δ φ}$ was measured to be 25 mrad, corresponding to a displacement sensitivity of 1.2 nm. This is similar to the phase sensitivity reported previously in OCT systems with galvanometer-based lateral scanning [28]. The increased phase noise is due to galvanometer jitter between B-scans and is not related to oversampling, as no scanning is performed in the y-direction.

For the elastograms presented in Section 5.2, the variance of the measured displacement was reduced by averaging 100 B-scans at each square wave half-cycle prior to calculation of the phase difference. The averaging was performed in the complex plane, before extracting the phase, in order to avoid the systematic error caused by the 2π phase modulus [18,29].

Both phase and strain sensitivity are determined by the OCT SNR (Eq. (5)), which is strongly affected by speckle and signal attenuation. To improve the strain sensitivity, Gaussian smoothing (GS) was applied by convolving the elastogram with a Gaussian kernel with full-width at half maximum (FWHM) of 20 μm × 20 μm.

4.1. Phantom fabrication

Tissue-mimicking phantoms with approximate dimensions (length × width × height) of 10 mm × 10 mm × 2 mm were fabricated using room-temperature vulcanizing (RTV) silicone as the bulk medium and titanium dioxide (TiO₂) as optical scatterers. Three phantom designs were fabricated. Phantom 1 was designed to be optically and mechanically homogeneous and relatively soft. It was fabricated from soft silicone (Elastosil® P7676, Wacker, Germany) and contained 1.5 mg/ml of TiO₂ particles. Its elastic modulus was measured using an Instron materials testing system to be 25 kPa. Phantom 2 was designed to be a soft bulk medium containing hard inclusions. It was fabricated using the same silicone used in Phantom 1. The inclusions were fabricated using a hard silicone (Elastosil® RT 601, Wacker, Germany), with measured elastic modulus of 4 MPa. Inclusions of measured average diameter of 0.4 mm were cut from a block by hand using a scalpel under a dissecting microscope. The inclusions contained 3.5 mg/ml TiO₂ particles. They were embedded in the soft surrounding silicone using a two-stage curing process. A 1.5 mm-thick layer of soft silicone was cured in a mold. The inclusions were placed on the surface of the soft layer using a hypodermic needle. A second layer of soft silicone, from the same batch as the first layer, was poured on top of the inclusions and allowed to cure, encasing them in the soft medium of total thickness 2 mm. Phantom 3 was designed to be a bilayer phantom with hard material of thickness 0.6 mm on top of soft material of thickness 1.4 mm. The top layer was fabricated using RT 601 silicone and contained 1.5 mg/ml TiO₂ scatterers. The bottom layer was fabricated using P7676 silicone and contained 3.5 mg/ml TiO₂. The elastic moduli of the hard and soft layers were measured to be 4 MPa and 25 kPa, respectively.

5. Experimental results

5.1. Comparison between theoretical predictions and experiment

In this section, experimental results are presented and compared to theory. Experiments were performed using Phantom 1 with uniform elastic properties so that the strain introduced was constant for all positions in the sample. To maximize phase sensitivity, no lateral scanning was performed. The theoretical and measured strain sensitivity, mean and SNR are presented in Fig. 3 . The measured parameters were calculated from experimental data for FD (red dots), OLS (blue dots), WLS (green dots) and Gaussian smoothed weighted least squares (GS-WLS, black dots) strain estimation methods. For each method, the mean and standard deviation of 100 calculations from consecutive axial positions within the sample were obtained, for strain resolutions from zero to 120 μm. The strain sensitivity and SNR were calculated from the experimental data using Eqs. (23) and (26). The theoretical predictions of the strain sensitivity and SNR are shown for the FD (solid red) and OLS (solid blue) methods. To best compare experiment and theory, the standard deviation of the experimentally measured displacement was substituted for $σ_{D_{i}}$ in the expressions for $Ε_{i}^{f}$ and $Ε_{i}^{o}$ defined in Sections 3.1 and 3.2. In Fig. 3, close agreement is seen between experiment and theory for both FD and OLS strain estimation methods. As expected, for both strain sensitivity (Fig. 3(a)) and strain SNR (Fig. 3(c)), the WLS method provides better performance than FD (by ~12 dB) and OLS estimates (by ~4 dB). The resultsobtained using GS-WLS strain estimation are also presented in Figs. 3(a)-3(c). The smoothing results in a ~3 dB improvement in strain sensitivity and SNR in comparison to WLS estimation. The mean strain calculated for each method is presented in Fig. 3(b). For strain resolutions >40 μm, the mean strain converges to −38 dB for all methods. This value was, therefore, used in the simulations. Below ~40 μm, the mean strain measured with all methods varies by ~3 dB, suggesting that below this threshold the measured strain is unreliable.

Fig. 3 (a) Strain sensitivity; (b) mean strain; and (c) strain SNR of Phantom 1. Experimental results are presented for FD (red dots), OLS (blue dots), WLS (green dots), and GS-WLS (black dots) strain estimation. Theoretical results are presented for FD (red line) and OLS (blue line) strain estimation.

Download Full Size | PDF

5.2. Elastograms

In this section, elastograms and the corresponding structural OCT images are presented for Phantoms 2 and 3, and for an excised porcine airway sample. For each image, the scale bars correspond to physical depth, assuming a refractive index of 1.4, the strain axial resolution is 75 μm, and the lateral resolution is 15 μm. Imaging and compression were introduced from the top in all elastograms. Each sample was tilted with respect to the beam axis to minimize spectral reflection from the air-sample interface. To maintain uniform compression, the ring actuator compression plate and the upper rigid boundary were tilted by the same angle in both x and y planes. In each elastogram, the strain SNR is plotted. Each pixel in an elastogram corresponds to the local strain divided by the noise floor, which was determined from the standard deviation calculated in an area (50 μm × 50 μm) with uniform intensity in the OCT structural image.

5.2.1. Phantom 2: Soft with hard inclusions

In Fig. 4(a) , the OCT structural image of a central portion of Phantom 2 is presented. A rectangular-shaped inclusion is visible commencing at ~200 μm below the surface. In Figures. 4(b)-4(e), the corresponding elastograms are presented for (b) FD, (c) OLS, (d) WLS and (e) GS-WLS strain estimation methods. Using FD strain estimation (Fig. 4(b)), it is difficult to distinguish the inclusion due to the low strain SNR. The result obtained with OLS, Fig. 4(c), is significantly improved, as expected, and the inclusion is clearly distinguishable. The strain SNR is further improved using WLS, as evident in Fig. 4(d).Large striations or bands are evident in the elastograms presented in Figs. 4(b)-4(d). The origin of this banding is discussed in Section 6. Its effect can be reduced using Gaussian smoothing, shown in Fig. 4(e). The tradeoff inherent in Gaussian smoothing is the improvement in SNR at the expense of spatial resolution. In Fig. 4(e), the axial resolution is marginally reduced from 75 μm to 78 μm, but the lateral resolution is degraded from 15 μm to 25 μm. Figures 4(f)-4(j) further illustrate the performance of the different strain estimation methods, showing the variation of OCT and strain SNR over the same lateral range as the images, at the depth indicated by the blue arrow in Fig. 4(a). The high spatial frequency noise in FD, OLS and WLS strain estimation is visible in Figs. 4(g)-4(i), as is its reduction in Fig. 4(j).

Fig. 4 Phantom 2: (a) OCT structural image; and (b) FD; (c) OLS; (d) WLS; and (e) GS-WLS elastograms. In (f)-(j), corresponding lateral traces are shown for the depth indicated by the blue arrow in (a).

Download Full Size | PDF

The low OCT SNR in Fig. 4(a) below the hard inclusion is a shadow artifact caused by signal attenuation in the highly scattering inclusion. High strain is observed in this location in all elastograms. From a mechanical standpoint, this is expected due to the increased stress caused by the inclusion. However, as the OCT SNR decreases, the probability distribution of the measured phase difference tends toward uniform [22]. This reduction in the mean phase difference causes the measured strain to increase below the inclusion. This is an artifact in strain estimation and is visible in all elastograms in Fig. 4. It could potentially be removed using additional thresholding to set a lower limit on OCT SNR at which strain can be measured reliably.

5.2.2. Phantom 3: Hard/soft bilayer phantom

In Fig. 5(a) , the OCT structural image of Phantom 3 shows the layer boundary at a depth of ~600 μm. In Fig. 5(b), the corresponding elastogram, generated using GS-WLS, shows that the layers can also be distinguished based on strain. The low strain SNR in the hard first layer suggests that it is displaced predominantly as a bulk, with most of the stress being dissipated in the underlying soft layer, as expected [20,30]. The banding effect described in Section 5.2.1 is also visible in Fig. 5(b), and will be discussed in Section 6. It can be seen in Fig. 5(b) that the strain in the second layer is non-uniform. This is due to two factors. Firstly, the silicone, like many soft tissues, is incompressible and therefore has a Poisson’s ratio of ~0.5. In this case, axial deformation detected using phase-sensitive OCE, is accompanied by a corresponding lateral deformation. In the softer second layer, near the boundary with the stiffer first layer, the lateral deformation of the sample is restricted by adhesion to the stiffer layer. This causes lower strain in the softer layer near the boundary, resulting in a higher apparent Young’s modulus. Secondly, the lower axial resolution in the elastogram reduces the sharpness of interfaces and contributes to the strain gradient in the second layer in Fig. 5(b).

Fig. 5 (a) OCT structural image; and (b) GS-WLS elastogram of Phantom 3. (c) OCT structural image; and (d) GS-WLS elastogram of a thin section of excised porcine airway with layers as labeled in (c).

Download Full Size | PDF

5.2.3. Porcine airway

To demonstrate an elastogram of tissue, we performed measurements on fresh ex vivo porcine airway harvested from a healthy animal. The airway was sectioned into a flat 10 mm × 10 mm × 2 mm sample. The sample was positioned such that the inner wall of the airway was in contact with the actuator compression plate. An OCT image of the airway is presented in Fig. 5(c). Three distinct layers are visible and readily attributable to the expected morphology. The highly scattering, most superficial layer is the mucosa; the middle low-scattering layer is the cartilage; and the deepest layer is the adventitia. The layers visible in Fig. 5(c) correspond well with those reported in previous OCT airway studies [31,32]. In Fig. 5(d), the GS-WLS elastogram is presented (generated in the same manner as for Figs. 4(e) and 5(b)). The three layers are apparent in the elastogram. Both the innermost mucosal layer and the outermost adventitial layer have relatively high strain, in comparison to the stiff cartilage, which has relatively low strain. The high stiffness of the cartilage is expected due to its composition, largely of collagen fibers [33].

6. Discussion

In any imaging technique, a framework characterizing performance, defining parameters such as spatial resolution, sensitivity, SNR, and dynamic range, enables assessment and optimization of image quality. Such a framework has played an important role in the closely related technique of ultrasound elastography [19,24,26]. In this paper, we have presented the first such framework for OCE in considering the performance of three strain estimation methods: two previously proposed, based on FD and OLS strain estimation, and a method that we have proposed for the first time in OCE, WLS strain estimation. Importantly, the application of these methods in OCE is distinct to their application in ultrasound elastography because of the different noise characteristics and processing steps used in each case. For example, in ultrasound elastography, the variance of the measured displacements is determined by parameters such as the transmitter bandwidth and the correlation coefficient [26], whereas in phase-sensitive OCE, statistical optics determines this variance [22].

In the elastograms presented, a banding effect is visible, i.e., the modulation between high and low strain in adjacent axial regions of elastograms. We demonstrated that Gaussian smoothing reduces this effect, at the expense of degradation of the spatial resolution. Future directions will include the investigation of more sophisticated filters to mitigate this degradation, e.g., edge-preserving bilateral filters [34]. These banding artifacts are similar in appearance to those reported in ultrasound elastography [2]. Their origin has been attributed to cyclic bias error in displacement measurements and correlated noise patterns arising from the use of oversampling. This is the first time such artifacts have been reported in OCE, and further work is required to confirm their origin and develop methods for their mitigation. We anticipate that these artifacts are related to the change in the probability distribution of the phase difference measurements when the condition $S N R_{O C T} > > 1$ is not satisfied. Under these conditions, Eq. (5) is not valid and the phase difference distribution is uniform rather than Gaussian [22]. This introduces a systematic bias error in the measured displacement, as the mean displacement, $μ_{D_{i}}$ , is zero rather than the true displacement. This bias influences the estimated strain and may explain the observed banding effect. Low SNR largely results from attenuation as light propagates through a sample, which could explain why the banding effect is more prominent at greater depths in the elastograms presented in Figs. 4 and 5.

7. Conclusions

We have presented a theoretical framework for strain estimation in OCE based on FD, OLS and WLS methods. We have defined, for the first time, expressions for strain sensitivity, SNR and dynamic range. For FD and OLS strain estimation, we demonstrated both close agreement between theory and experiment, and the advantage of OLS strain estimation. Our experimental results indicate that strain estimation using WLS provides more accurate results than using FD (by ~12 dB) or OLS (by ~4 dB). We demonstrated that Gaussian smoothing conveyed an additional improvement in strain sensitivity at the expense of lateral resolution. The framework presented provides a quantitative means to relate phase sensitivity to elastogram image quality in phase-sensitive OCE. We expect that this framework will allow quantitative comparisons to be made between different strain estimation methods, as well as different OCE experimental procedures, greatly facilitating the future development of OCE.

Acknowledgments

This project was funded in part by the National Breast Cancer Foundation, Australia. BFK acknowledges funding received from the Raine Medical Research Foundation; KMK from the Scholarship for International Research Fees and University International Stipend, UWA; RAM from the Cancer Council Western Australia; and PRTM from the Australian Research Council Discovery Early Career Research Award.

References and links

1. J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imaging elastic properties of biological tissues,” Annu. Rev. Biomed. Eng. 5(1), 57–78 (2003). [CrossRef] [PubMed]

2. J. Ophir, S. K. Alam, B. Garra, F. Kallel, E. Konofagou, T. Krouskop, and T. Varghese, “Elastography: ultrasonic estimation and imaging of the elastic properties of tissues,” Proc. Inst. Mech. Eng. H 213(3), 203–233 (1999). [CrossRef] [PubMed]

3. J. M. Schmitt, “OCT elastography: imaging microscopic deformation and strain of tissue,” Opt. Express 3(6), 199–211 (1998). [CrossRef] [PubMed]

4. R. C. Chan, A. H. Chau, W. C. Karl, S. Nadkarni, A. S. Khalil, N. Iftimia, M. Shishkov, G. J. Tearney, M. R. Kaazempur-Mofrad, and B. E. Bouma, “OCT-based arterial elastography: robust estimation exploiting tissue biomechanics,” Opt. Express 12(19), 4558–4572 (2004). [CrossRef] [PubMed]

5. J. Rogowska, N. A. Patel, J. G. Fujimoto, and M. E. Brezinski, “Optical coherence tomographic elastography technique for measuring deformation and strain of atherosclerotic tissues,” Heart 90(5), 556–562 (2004). [CrossRef] [PubMed]

6. H. J. Ko, W. Tan, R. Stack, and S. A. Boppart, “Optical coherence elastography of engineered and developing tissue,” Tissue Eng. 12(1), 63–73 (2006). [CrossRef] [PubMed]

7. W. Drexler and J. G. Fujimoto, Optical Coherence Tomography: Technology and Applications (Springer-Verlag, Berlin, 2008).

8. B. F. Kennedy, T. R. Hillman, A. Curatolo, and D. D. Sampson, “Speckle reduction in optical coherence tomography by strain compounding,” Opt. Lett. 35(14), 2445–2447 (2010). [CrossRef] [PubMed]

9. B. F. Kennedy, A. Curatolo, T. R. Hillman, C. M. Saunders, and D. D. Sampson, “Speckle reduction in optical coherence tomography images using tissue viscoelasticity,” J. Biomed. Opt. 16(2), 020506 (2011). [CrossRef] [PubMed]

10. T. R. Hillman, A. Curatolo, B. F. Kennedy, and D. D. Sampson, “Detection of multiple scattering in optical coherence tomography by speckle correlation of angle-dependent B-scans,” Opt. Lett. 35(12), 1998–2000 (2010). [CrossRef] [PubMed]

11. R. K. Wang, Z. H. Ma, and S. J. Kirkpatrick, “Tissue Doppler optical coherence elastography for real time strain rate and strain mapping of soft tissue,” Appl. Phys. Lett. 89(14), 144103 (2006). [CrossRef]

12. R. K. Wang, S. Kirkpatrick, and M. Hinds, “Phase-sensitive optical coherence elastography for mapping tissue microstrains in real time,” Appl. Phys. Lett. 90(16), 164105 (2007). [CrossRef]

13. X. Liang, A. L. Oldenburg, V. Crecea, E. J. Chaney, and S. A. Boppart, “Optical micro-scale mapping of dynamic biomechanical tissue properties,” Opt. Express 16(15), 11052–11065 (2008). [CrossRef] [PubMed]

14. S. G. Adie, X. Liang, B. F. Kennedy, R. John, D. D. Sampson, and S. A. Boppart, “Spectroscopic optical coherence elastography,” Opt. Express 18(25), 25519–25534 (2010). [CrossRef] [PubMed]

15. X. Liang, S. G. Adie, R. John, and S. A. Boppart, “Dynamic spectral-domain optical coherence elastography for tissue characterization,” Opt. Express 18(13), 14183–14190 (2010). [CrossRef] [PubMed]

16. B. F. Kennedy, X. Liang, S. G. Adie, D. K. Gerstmann, B. C. Quirk, S. A. Boppart, and D. D. Sampson, “In vivo three-dimensional optical coherence elastography,” Opt. Express 19(7), 6623–6634 (2011). [CrossRef] [PubMed]

17. B. Park, M. C. Pierce, B. Cense, S. H. Yun, M. Mujat, G. Tearney, B. Bouma, and J. de Boer, “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 µm,” Opt. Express 13(11), 3931–3944 (2005). [CrossRef] [PubMed]

18. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley, New York, 1998).

19. F. Kallel and J. Ophir, “A least-squares strain estimator for elastography,” Ultrason. Imaging 19(3), 195–208 (1997). [PubMed]

20. B. F. Kennedy, T. R. Hillman, R. A. McLaughlin, B. C. Quirk, and D. D. Sampson, “In vivo dynamic optical coherence elastography using a ring actuator,” Opt. Express 17(24), 21762–21772 (2009). [CrossRef] [PubMed]

21. A. Grimwood, L. Garcia, J. Bamber, J. Holmes, P. Woolliams, P. Tomlins, and Q. A. Pankhurst, “Elastographic contrast generation in optical coherence tomography from a localized shear stress,” Phys. Med. Biol. 55(18), 5515–5528 (2010). [CrossRef] [PubMed]

22. J. W. Goodman, Statistical Optics (Wiley, New York, 1985).

23. J. Ophir, I. Céspedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastography: a quantitative method for imaging the elasticity of biological tissues,” Ultrason. Imaging 13(2), 111–134 (1991). [CrossRef] [PubMed]

24. I. Cespedes, M. Insana, and J. Ophir, “Theoretical bounds on strain estimation in elastography,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 42(5), 969–972 (1995). [CrossRef]

25. J. Fox, Applied Regression Analysis, Linear Models and Related Methods (Sage, Thousand Oaks, Calif., 1997).

26. T. Varghese and J. Ophir, “A theoretical framework for performance characterization of elastography: the strain filter,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 44(1), 164–172 (1997). [CrossRef] [PubMed]

27. J. I. Jackson and L. J. Thomas, “Ultrasound-based strain rate estimation of moving, fully developed speckle,” in 2001 IEEE Ultrasonic Symposium (IEEE, 2001), Vol. 2, pp. 1593–1596.

28. M. Pircher, B. Baumann, E. Götzinger, H. Sattmann, and C. K. Hitzenberger, “Phase contrast coherence microscopy based on transverse scanning,” Opt. Lett. 34(12), 1750–1752 (2009). [CrossRef] [PubMed]

29. A. Szkulmowska, M. Szkulmowski, A. Kowalczyk, and M. Wojtkowski, “Phase-resolved Doppler optical coherence tomography—limitations and improvements,” Opt. Lett. 33(13), 1425–1427 (2008). [CrossRef] [PubMed]

30. S. G. Adie, B. F. Kennedy, J. J. Armstrong, S. A. Alexandrov, and D. D. Sampson, “Audio frequency in vivo optical coherence elastography,” Phys. Med. Biol. 54(10), 3129–3139 (2009). [CrossRef] [PubMed]

31. G. J. Tearney, M. E. Brezinski, B. E. Bouma, S. A. Boppart, C. Pitris, J. F. Southern, and J. G. Fujimoto, “In vivo endoscopic optical biopsy with optical coherence tomography,” Science 276(5321), 2037–2039 (1997). [CrossRef] [PubMed]

32. S. Han, N. H. El-Abbadi, N. Hanna, U. Mahmood, R. Mina-Araghi, W. G. Jung, Z. Chen, H. Colt, and M. Brenner, “Evaluation of tracheal imaging by optical coherence tomography,” Respiration 72(5), 537–541 (2005). [CrossRef] [PubMed]

33. J. K. Rains, J. L. Bert, C. R. Roberts, and P. D. Paré, “Mechanical properties of human tracheal cartilage,” J. Appl. Physiol. 72(1), 219–225 (1992). [PubMed]

34. C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” in Sixth International Conference on Computer Vision, 1998 (IEEE, 1998), 839–846.

Strain estimation in phase-sensitive optical coherence elastography

Abstract

1. Introduction

2. Phase-sensitive OCE strain

2.1. Displacement measurement

2.2. Strain estimation

2.2.1. Finite difference strain estimation

2.2.2. Ordinary least squares strain estimation

3. Strain imaging parameters

3.1. Strain imaging parameters using the FD method

3.2. Strain imaging parameters using the OLS method

3.3. Theoretical predictions

3.4. Weighted least squares strain estimation

4. Experimental setup and procedure

4.1. Phantom fabrication

5. Experimental results

5.1. Comparison between theoretical predictions and experiment

5.2. Elastograms

5.2.1. Phantom 2: Soft with hard inclusions

5.2.2. Phantom 3: Hard/soft bilayer phantom

5.2.3. Porcine airway

6. Discussion

7. Conclusions

Acknowledgments

References and links

Cited By

Figures (5)

Equations (33)

Biomedical Optics Express