1. Introduction

The determination of elastic and viscoelastic properties of soft tissues is of fundamental interest in the detection and diagnosis of many diseases, since these properties depend in a very significant manner on the pathological state of the tissue. As is well known from manual palpation, many pathological tissue regions have a different strain response to an imposed stress than does the surrounding healthy tissue of the same type. Indeed, this is the exact reason why manual self examination is effective in locating breast and testicular lesions.

Elastography in general is a relatively new method for the in vivo quantitative imaging of strain and elastic modulus distributions in tissue [1]. The modality is based on the estimation of strain due to tissue compression or expansion [2] and traditionally employs either ultrasonic or magnetic resonance signals to carry the information and to form the images, known as elastograms. Both ultrasonic and magnetic resonance elastography have been used to quantify the elastic behavior of breast tissue and breast lesions [3–9] as well as prostate tissue for the detection of prostatic cancer [10, 11].

Optical elastography employs optical frequency radiation to detect and image the strain response of a tissue subjected to mechanical or acoustic forces. Various implementations of optical elastography have been used to evaluate the mechanical behavior of ex vivo cortical bone [12], chicken skin [13], human and rat arteries [14], elastin [15], and porcine skin [16]. Elastographic methods have also been combined with optical coherence tomography (OCT). Schmitt [17] was the first to implement such a system and investigated the microscopic deformation of gelatin tissue phantoms, pork meat, and human skin as a function of depth under a compressive load. OCT-based elastography has also been applied to the investigation of arterial wall biomechanics [18, 19], atherosclerotic plaques [20, 21], and engineered and developing tissues [22]. All of these OCT-based elastography applications have relied upon 2-dimensional cross-correlation algorithms to quantify the speckle shifts and thereby estimate the motions and velocities in the OCT image sequences. Chan et al., [18], however, employed a priori knowledge regarding the velocity fields in the pulsating arterial wall to add robustness to the measurements while avoiding image smoothing that may reduce the information content of the final elastograms.

In OCT-based elastography, however, cross-correlation approaches for estimating the shift in the speckle patterns are of limited usefulness both when the deformations are either very small [23,24] or very large between frames. For example, when the speckle motion between successive frames (images) is only a fraction of a pixel, cross-correlation approaches lack the resolution to robustly track the shift. Such a case is experienced when the random motion in hydrated tissues leads to rapid speckle decorrelation, thus dictating a high enough sample rate to ensure that the speckles remain correlated at least through several image acquisitions. Unless the loading rate is rapid with respect to the acquisition rate, the speckle motion will be small between frames. Small speckle motions are also observed when the tissue is exposed to small loads, or is very stiff.

At the other extreme, without the a priori knowledge of [18], the basic assumption for these correlation-based approaches is that the speckle in the OCT image produced from the coherence-volume is stable, except for a lateral translation, between successive B-scans (i.e., a frozen speckle model [23]). This assumption limits these techniques to motions between the B-scans to somewhat less than the spatial coherence length of the light source used. However, relatively large tissue movement rapidly decorrelates the speckle from one frame to another, limiting the large-motion applicability of correlation based approaches. In addition, the speckle motion in correlation based algorithms is typically found by searching for the lag at which the cross-correlation between two moving windows over the OCT images is a maximum. This is a time consuming procedure, making the real time, on-line measurement difficult.

Herein we present two alternative numerical approaches to speckle tracking, and ultimately the generation of elastograms, in OCT-based elastography; one is particularly useful for small speckle motions and the other for large and/or fast speckle motions. The experimental conditions, of course, dictate which approach is more appropriate for a given situation. We demonstrate our elastography approaches through a series of experiments during which 2-layered polyvinyl alcohol tissue phantoms are subjected to a dynamic compressive load while being imaged with a spectral-domain OCT system.

2. Theory

2.1 Small speckle motions

The approach for small speckle motion estimation is essentially a 2-dimensional implementation of the one-dimensional maximum-likelihood estimator that we previously described [23, 24]. The problem is to assess the motion from frame to frame of the features in a sequence of x-z OCT images, where the x-direction is a lateral dimension along the surface and the z-direction is the depth dimension. We denote the intensity or gray level variations in this k ^th image as g(x,z,k), and wish to determine the time progression of the lateral shift of features in this image. Towards this end, we inspect this shift based on the pair of images prior to and after the k ^th image. For this image, we define a mean square difference between a portion of the pair of adjacent images;

The summation is over a small neighborhood, say of dimension 3×3, of the images.

The objective is to determine the local velocity, f̄=(f_x, f_z) that will bring the features in this neighborhood pair into registration. Note that the units of velocity are implicit, i.e., pixels/record. Obviously, the shift estimates are specified for the neighborhood chosen by the summation in Eq. (5), for example, 3×3. As detailed previously [23, 24], expanding each image in terms of a two-dimensional Taylor series results in the expression

where

In the above, g _x and g _z are respectively the average x- and z-gradients at the k ^th record, and g _t is the central difference approximation to the temporal gradient. Minimization of this error by taking the partial derivatives with respect to the velocity components results in a pair of simultaneous equations,

where 〈…〉 represents the local averaging operation. These equations are easily inverted to yield estimates of the velocity, f̄. Note that each element in Eq. (7) is a full size N×M array. Solution of these linear simultaneous equations produces velocity estimates for each pixel of each plane of the image cube.

More formally, Eq. (6a) is expressed as

where the window function, w, denotes the local averaging operation. In its simplest form, this window function is a simple boxcar that evenly weights the residuals within the region of interest.

An alternative derivation of this equation is based on the concept of optical flow. This concept relies on the assumption that the brightness, g(x,z,t), is conserved (only its spatial distribution changes with time) so that the total time derivative, dg/dt, is zero, i.e.,

where f=[f_x f_z]^T is the optical flow. Note that this equation (commonly known as the brightness change constraint equation - BCCE [25]) characterizes a single point on the image. We presume that the image gradient, ∇g, and the partial derivative with respect to time, g_t, are known, i.e., measured. We have a single equation [Eq. (9) represents a single constraint] and two unknowns, the velocities f_x & f_z. This is referred to as the aperture problem of motion estimation. Minimizing the residual within the region specified by the function, w (the aperture), leads to the previous result, Eq. (8).

Once the velocities (f_x, f_z) are estimated they can be integrated over time and mapped on a pixel-by-pixel basis to create an image whose gray values encode displacements in either the x or z direction, or the vector sum via the Pythagorean theorem. Normalizing these displacement maps by the initial dimensions of the sample yields elastograms that encode strains in the limit of small deformations (i.e., infinitesimal strains) [26].

Visualization of regions that display a different strain response is enhanced through the use of a neighborhood operator in the form of a discrete convolution filter [27]. Pixel operations as described above are not particularly suitable for discriminating regions of interest in images because the gray value of each pixel is determined with no consideration of the gray values of the neighboring pixels. In contrast, neighborhood operations, such as convolution filters, analyze the spatial relations of the gray values. Thus, in effect, this operator converts the gray scale elastograms into feature-based elastograms [27]. The convolution kernel (mask) was adaptively generated directly from the gray values of a small background region of the strain-encoded elastogram. The elastogram was then convolved with this kernel in 2 dimensions to emphasize the local magnitude of the strains.

2.2 Large speckle motions

Our approach to measure the large speckle motion is based on the Doppler effect induced by the tissue motion, called tissue Doppler optical coherence elastography (tDOCE) [28]. Because the basic principle of OCT is based on the interference between the reference light and the backscattered light from within a sample, it is extremely sensitive to any particle movement in the sample which causes a Doppler shift in the resulting interference signal. This phenomenon was first utilized in time-domain OCT (TdOCT) [29] and later in frequency-domain (FdOCT) systems [30] to accurately measure fluid flow, including blood flow. Clearly, because of light scattering from tissue, any tissue movement also induces a Doppler shift in the interference signals. This tissue Doppler effect, to our knowledge, has not yet been exploited in the OCT community. Our method to measure the large movement of the speckle was developed based on the measurement of the Doppler frequency shift induced by instantaneous tissue motion. Consequently the displacement, strain rate and strain imaging of local tissue can be quantified accurately in real time. Because FdOCT has inherent advantages over TdOCT in terms of imaging speed, detection sensitivity and phase stability [31], we will describe tDOCE used in this study through the use of a FdOCT system based on the spectral domain implementation. However it should be appreciated that the theory applies equally to TdOCT systems.

In spectral domain OCT, the interference signal between the reference light and the scattered light from within a sample is spectrally resolved by a linear array detector, usually a line scan CCD camera. The light intensity incident on each element of the linear array detector is proportional to the spectral density I_d(υ) of the combined reference and sample light, which can be expressed as:

where υ is the frequency of the light source. Range information is determined from the propagation time τ of the light backscattered by the scatters within the sample, R(τ) is the normalized intensity backscattered from the scatter at timeτ, and S(υ) is the spectral density of the light source. The 3^rd term in the brackets represents the self interference, R_self of the backscattered light between different scatterers in the sample. We assume the reflection of the reference arm to be unity for clarity. The 4^th term contained within the brackets is the spectral interference between the sample and the reference and is the term that contributes to the OCT signal. Thus, we disregard the first 3 terms within the brackets and concentrate on the final term yielding,

The depth information regarding the local scatters in the sample is obtained by Fourier-transforming the spectral interference signal above:

Note that Eq. (12) is a complex function where the amplitudes are used to generate the traditional OCT structural images. The phase term Φ(z) is generally random along the depth z, but fixed in position for the (relatively) static scatters. A translation of the scatterer at position z by a distance Δd(z) during the time interval Δt between two successive A-scans will induce a change in the measured phase of the reflected light given by

where n is the refractive index of the sample and k is the wavenumber of the light source (k=2π/λ=2πυ/c). Calculating this phase difference at each depth z yields depth-resolved measurements of both the magnitude and direction of the axial (parallel to the imaging beam) displacement of the tissue at the time t [28],

Similar to the phase-resolved optical Doppler tomography where the concern is the flow velocity [29, 30], the depth-resolved tissue velocity v(z,t) in the beam direction can be obtained by

To ensure correlation between the phase measurements of successive A-scans, the transverse displacement of the imaging beam between A-scans must be small enough relative to the probe beam size to avoid rapid spatial speckle decorrelation. This constraint can be met by effectively dense sampling in the transverse direction. After the depth-resolved instantaneous displacement and velocity are obtained, the strain rate map $\dot{\epsilon }$(z,t) can be generated and color coded to represent the localized elastic properties of tissue,

where z ₀ is the initial depth of the sample before the tissue movement, i.e., before the compression and the dot over the character implies the time derivative. So far, we have derived the depth-resolved instantaneous displacement (Eq. 14) and strain (Eq. 15) of the sample at the time t. The total displacement and strain over a time period T can therefore be obtained by integration of the instantaneous displacement and strain over the elapsed time, respectively. As a consequence, the depth-resolved displacement d(z) and strain ε(z) maps of the sample over the time duration T can be written as:

Unlike the optical coherence elastography (OCE) approaches previously developed where maintaining correlated speckle between successive OCT B-scans is required [17, 18, 22], the current approach only requires that the speckle patterns between the successive A-scans are correlated and the calculation of tissue motion and strain maps from the phase information is straightforward and rapid. With our A-scan rate of 29.2kHz (see below), tDOCE presents an opportunity to map the elastic properties of in vivo tissue in real time.

3. Methods

3.1 OCT system

The OCE system used in this study includes two main parts, the spectral domain OCT system and the loading rig that was used to apply a dynamic force to compress the tissue. The OCT system has been described in detail previously [32]. It uses a super-luminescent diode with a central wavelength of 842nm and FWHM bandwidth of 50 nm that yields a measured axial resolution of ~8 µm in air. After passing through an optical isolator, the light was coupled into a fiber-based Michelson interferometer. The reference light was delivered to a double-pass grating based rapid scanning system [33]. However, unlike TdOCT systems, the reference mirror in the scanning system was kept stationary. This double pass grating-based system was used to efficiently compensate for the second-order dispersion in the system. The sample light was coupled into a probe, consisting of an x-scanner and the optics to deliver the probing light onto, and collect the light backscattered from, the sample. The lateral imaging resolution was approximately 20 µm as governed by the focusing optics. The light returning from the reference and sample arms was recombined and sent to a custom designed high speed spectrometer, consisting of a 30 mm focal length collimator, a 1200 lines/mm diffraction grating and an achromatic focusing lens with 150 mm focal length. The focused light spectrum impinged on a line scan CCD camera consisting of 2048, 10 µm×10 µm pixels with 10 bit digital depth, and capable of a 29.2 kHz line rate. Polarization controllers were used in the reference, sampling and detection arms in order to maximize the interference fringe contrast at the detector. The spectrometer has a designed spectral resolution of 0.055 nm, resulting in an imaging depth of approximately 3.2 mm in air. The signal sensitivity of 95 dB was measured at z=0.5 mm and dropped to 80 dB at z=2.0 mm when the camera integration time was set at 34.1 µs. The B-scan OCT images consisting of 500 A-scans, corresponding to 2.5 mm, were obtained by lateral scanning of the probe controlled by the x-scanner. The A-scan spatial interval was 5 µm ensuring that the speckle in the interferograms was correlated between successive A-scans, because the probing beam spot was 20 µm (transverse resolution). The imaging rate was set at 20 frames per second (fps) with τ=100 µs for this investigation, although 58 fps is achievable using this system [32].

The loading fixture consisted of a calibrated load cell mounted on a precision motorized translation stage that was used to compress the sample. Both the signals from the load cell representing the force applied and the translation stage indicating the actual displacement of the sample were acquired by an analog-to-digital converter and stored in the computer for later correlation with the OCT measurements. A transparent glass slide (with a thickness of ~1.5 mm) was glued to a fixed stand that had a clear aperture of 1 cm diameter in the middle that allowed the OCT probe beam to pass through freely. The tissue phantom was sandwiched between the glass slide and the load cell. The computer-controlled translation stage moved the load cell and applied the necessary dynamic compression to the sample against the fixed glass slide.

3.2 Experiments

The tissue phantoms used in this study were made of 20% w/w polyvinyl alcohol (PVA) that was subjected to a series of freezing and thawing cycles. Repeated freezing and thawing increases the dynamic mechanical stiffness and also increases the reduced scattering coefficient of the material [34]. We created a bi-layer phantom, with the top layer being subjected to 7 freezing-thawing cycles and the bottom layer only 5 cycles. Based on prior mechanical testing, the expected ratio of dynamic stiffness was on the order of 1.5:1 between the top and bottom layers, respectively.

For the small strain case, the bi-layer phantom was compressed using a triangular waveform with a peak-to-peak amplitude of 5 µm at a rate of 4 µm/s. Total compressive stress was on the order of 500 Nm^-2 as determined by dividing the output of a calibrated load cell in series with the phantom by the cross-sectional area of the phantom.

For the large deformation case, an identical phantom was dynamically compressed using a triangular waveform with a peak-to-peak amplitude of 400 µm at a rate of 400 µm/s repeatedly for 8 cycles. The phantom was allowed to rest for ~5 seconds, during which time the OCT system continuously acquired the velocity and strain maps of the phantom at a rate of 20 frames/s during the loading, unloading and rest periods.

Note that a difference in the strain response of the two layers in the direction of the applied load is only observable because the layered phantom is arranged in series with the loading stage (Maxwell model). Had the two layers been arranged in parallel with one another (i.e., Voigt or Kelvin model), a differential strain response in the direction of the load between the two regions should not be observed, as a parallel arrangement demands that the two regions deform identically for a uniformly applied load [26]. If the compression tests are performed in an un-confined manner, however, differences in the strain response in the direction orthogonal to the loading direction may still be observed in the Voigt arrangement, reflecting a difference in the Poisson’s ratio of the two mechanically disparate regions.

4. Results

4.1 Small deformations

A spectral OCT video of the bi-layered tissue phantom being compressed is shown in Fig. 1. Compression was in the direction of the arrow. The difference in the scattering properties of the two layers is seen in the different intensities between the two layers. This disparity in scattering, however, has no influence on the generation of the elastograms. Indeed, the constraints on the image quality and resolution are minimal for elastography. All that is required is that the speckles are resolved. No effort (e.g., image averaging) was made to reduce the speckle in the OCT images and the speckles are quite obvious.

 

Fig. 1. (2.25 MB) OCT video sequence of the bi-layered tissue phantom under compression showing small deformations (9.0 MB version).

Download Full Size | PDF

We are concerned with deformations in the direction of the applied force (z-direction) because this is the only direction in which we could make knowledgeable predictions of how the bi-layered phantom should deform. The analysis was performed over the first 75 frames of the video shown in Fig. 1. Note that the large rigid-body motions in the x-direction observable in the video sequence (Fig. 1) prevented a reasonable assessment of the strains in the x-direction because it was not possible to determine which motions should be attributed to strain and which should be attributed to rigid body motion. Rigid body motion was not found to be a significant source of speckle motion in the z-direction. The video sequence of Fig. 2 displays the cumulative strain-encoded elastogram as the load on the phantom is increased. Using the infinitesimal strain assumption, the integrated deformations, ∫f_zdt, were normalized by the thickness of the sample to arrive at cumulative strain. The two regions of the phantom are readily seen in the strain-encoded video sequence (Fig. 2), with the stiffer, upper region displaying a smaller strain response than the less-stiff lower region. The interfacial region between the two layers is of intermediate stiffness and is a result of the manufacturing process when warm, liquefied PVA is pored onto the solidified layer that has already gone through 2 freezing/thawing cycles. The mean cumulative strain in the stiffer (top) region was ~170 µε and that in the less stiff bottom region was on the order of 280 µε giving a ratio of strains between the two areas of 1.6, which is in good agreement with our expected ratio of approximately 1.5 (see above).

 

Fig. 2. (1.7 MB) Movie of strain development in the bi-layered tissue phantom as the phantom was compressed.

Download Full Size | PDF

 

Fig. 3. Neighborhood operation on the elastogram of cumulative strain. Figure 3(a). (left) is the elastogram of the total cumulative strain as determined by the pixel-by-pixel speckle motion estimator. Figure 3(b). (center) is an enlarged picture of the 40×40 pixel convolution kernel used in the neighborhood operation. The kernel was taken directly from gray-scale values within the small box outlined in Figure 3(a). Figure 3(c). (right) is the final featurebased elastogram encoded to display the relative cumulative strains in the different layers of the tissue phantom. The highest strain in the less-stiff region was normalized to unity. The mechanical distinction between the two layers is evident and the interfacial region is quite visible (greenish-blue).

Download Full Size | PDF

To enhance the visual discrimination between the two layers of the phantom in the final elastogram, we applied a neighborhood operator in the form of a convolution filter to the raw elastogram that displayed the total cumulative strain response with the goal of highlighting the specific features (in this case the strain response) of the gray-scale elastogram. Figures 3(a)–3(c) displays the results of this procedure. Figure 3(a) is the total cumulative strain-encoded elastogram which was subjected to the discrete convolution operation with the kernel shown in Fig. 3(b). This kernel was taken directly from the gray scale values in the 40×40 pixel region of the elastogram indicated by the box in Fig. 3(a). While no rigorous effort was made to optimize the size and shape of the kernel, the final form of the kernel was determined adaptively by visually inspecting the resultant feature-based elastogram and modifying the dimensions of the kernel until a visually optimal feature-based elastogram [Fig. 3(c)] was generated. The strain values in the feature-based elastogram were normalized so that the highest strain in the image was given the value of unity. The two regions of the tissue phantom are clearly discriminated in this normalized, feature-based elastogram [Fig. 3(c)].

4.2 Large deformations

In the demonstration of our approach to measure large tissue deformations, we used the same bi-layer phantom as in our small deformation case. The bi-layer phantom was first pre-loaded with a force 0.3 kg, corresponding to a stress of ~2.94 N/cm², and then allowed enough time to relax the internal stress. The thickness of the phantom at this time (i.e. t=0) was measured as ~2 mm with a precision caliper. It was then compressed repeatedly as described in Section 3.2 by a triangular waveform for 8 cycles (16 s) and then allowed to rest for 5 s, during which time the OCT system continuously read out in real time the velocity maps of the instantaneous displacement, and strain rate maps of the phantom at an imaging rate of 20 fps during the loading, unloading and rest periods.

 

Fig. 4. Video sequences of the real time OCT structural image (a, 3.6 MB), velocity map (b, 2.7 MB), and strain rate map (c, 2.3 MB) of the bi-layer phantom during the dynamic loading, respectively. The OCT image was coded as gray-scale from 20 (black) to 50dB (White), the velocity was color coded with dark blue representing -120 µm/s (minimum) and purple-red representing +120 µm/s (maximum), and the strain rate map was color coded with dark blue representing -0.25 s^-1 and purple red representing +0.25 s^-1. The physical size of the images is 1.2×2.5 mm.

Download Full Size | PDF

Figure 4 shows the real time movies for the OCT image (a), velocity map (b), and strain rate map (c) when the sample was being loaded. Here the physical size of the images shown was 1.2×2.5 mm (depth x width). In generating these movies, a threshold of about 15 dB above the noise floor in the OCT image was used to calculate the phase changes s Δφ(z,t) between successive A-scans, which meant that the velocities were set to zero if the intensity of the OCT signal fell below this threshold value. A moving average window of 4×4 pixels, corresponding to a physical size of 20×20 µm, was used to smooth the phase map from which the velocity, instantaneous displacement and strain rate maps were obtained. As was expected, the values of the velocity map changed sign when the loading direction was reversed (red corresponds to loading, blue corresponds to unloading). More importantly, the absolute values in the top layer are less than those in the bottom layer, indicating that the top layer of the phantom is stiffer than that of the bottom layer, which was expected. From this real time velocity map, the instantaneous displacement, strain rate and strain maps were obtained according to Eq. (14) and Eq. (16), respectively. The instantaneous displacement and strain movies are not shown here because they are visually similar to the velocity map.

To further quantify the elastic properties of the bi-layer phantom used in this study, we collapsed each B-scan velocity image and averaged them into a single A-scan velocity profile. This collapsed A scan velocity profile was then plotted as a function of time as the load was applied to the sample. Figure 5(a) gives such a time varying velocity map over a time duration of 21 s, including the 8 loading and unloading cycles followed by a 5 s resting period. It is clearly seen that the tissue velocity alternates its sign during the loading and unloading periods, and then it approaches zero when at rest. In addition, the velocity in the top layer is seen to be slower than that in the bottom layer, which is expected because the top layer is stiffer than the bottom layer. By applying Eqs. (16), (17), and (18) to Fig. 5(a), the results can be obtained and color coded to display the time varying strain rate [Fig. 5(b)], displacement [Fig. 5(c)] and strain [Fig. 5(d)] maps, respectively. The units of the colormaps used in Figs. 5(a), 5(b), 5(c), and 5(d) were µm/s, s^-1, µm, and percentage (%), respectively.

 

Fig. 5. Integrated time varying (a) velocity (µm/s), (b) strain rate (s^-1), (c) displacement (µm), and (d) strain (%) maps, respectively, of the tissue phantom subjected to 8 loading cycles and then followed by a 5 s rest period.

Download Full Size | PDF

Figs. 6(a) and 6(b) give the plots for the representative displacements and strains against the time at 2 different depths in the sample, respectively. The plotted depths were z=0.29 mm which was located in the top layer, and z=0.55 mm which was located in the bottom layer. These were compared with the actual displacement at z=0.58 mm and the force measured on line by the load cell and the translation stage, respectively, synchronized with the OCT measurements. Note that a linear relationship was assumed for the displacement measurement from the translation stage, thus the actual displacement measured at z=0.58 mm was scaled down from the 2 mm initial thickness of the sample. From Fig. 6(a), the time varying displacement of the sample during the loading and unloading cycles closely approximate the triangular waveform at both depths with the softer bottom layer displacing more and the stiff top layer displacing less. Correspondingly, the strain measured in the bottom layer is higher than in the top layer. From this plot, the averaged ratio of the strains of the bottom layer to the top layer was calculated as 1.65±0.10:1 that agrees well with the expected value of 1.5:1 (above).

 

Fig. 6. (a). Displacement and (b). strain profiles plotted against time compared with the synchronized separate measurements of actual displacement [top curve in (a)] and force [bottom curve in (b)] applied to the phantom, respectively. Rest of curves from bottom to top are the depth profiles at z=0.29 mm (blue) and 0.55 mm (red), respectively.

Download Full Size | PDF

5. Discussion

We have demonstrated two different approaches for quantifying displacements, strains, and strain rates using spectral OCT-based elastography. Both approaches are variants on speckle tracking. The first approach, which is suitable for small frame-to-frame speckle motions, is based on a maximum likelihood speckle shift estimator. In the 1-D case, we have shown previously that this approach excels over standard cross-correlation approaches when the motion between frames is below approximately 0.8 pixels [24]. In the case of larger motions, the algorithm yields large errors, which can be reduced through the use of an iterative scheme [24]. However, the iterative operations are somewhat computationally inefficient. Thus, we presented an alternative algorithm for use when the motion is greater than approximately 1 pixel per frame, tissue Doppler optical coherence elastography (tDOCE). tDOCE has the distinct advantage that it does not rely upon the frozen-speckle assumption (for successive B-scans) of previous OCE implementations, allowing for much faster and large motions to be quantified. Furthermore, the computation involved is quite straightforward and very rapid, allowing for the generation of real time elastograms. Both algorithms yielded similar results for the relative stiffnesses of the two layers of the tissue phantom and both sets of results agreed with our predicted stiffness ratio between the two layers based upon previous traditional mechanical testing of the PVA phantom material in our laboratory.

The choice of algorithms for the generation of elastograms is clearly dictated by the experimental conditions. Cases exhibiting very small speckle motions, such as are observed when the applied loads are small or when the material is very stiff, require an efficient algorithm that doesn’t yield unacceptable errors when the frame-to-frame shifts fall below approximately 0.8 pixels [24]. In this case, the 2-D maximum likelihood estimator would be a proper choice. For very fast and/or large inter-frame motions, the tissue Doppler approach is more appropriate. In situations where the motions are moderate, at least 0.5 pixels per frame and within the spatial coherence length of the source (i.e., when there is minimal spatial decorrelation due to large motions), when the frozen speckle assumption applies, and when computational efficiency is not a pressing issue, a cross-correlation algorithm is a very intuitive option.