Optical elastography, as defined here, is the field of study that aims to use coherent light techniques to evaluate the mechanical behavior of materials. Of immediate interest is the use of laser speckle methods for investigating the mechanical behavior of biological tissues. The use of laser speckle techniques for the mechanical characterization of materials is well known in the nondestructive evaluation community.1 2 3 4 One particular application is to infer strain by monitoring the motion of the speckle pattern that results from coherently illuminating a stressed object. Typically a reference image of the speckle pattern is acquired before deformation of the object. Motion (with respect to this reference image) of subsequent speckle patterns, which occurs when the object is stressed, is used to infer the resulting strain. A problem experienced in using this technique for measurements of hydrated tissues is the rapid decorrelation of the speckle patterns.5 Thus, application of speckle techniques to assessment of strain in biological tissues often relies upon rapid sampling of the speckle patterns and the use of processing algorithms that are aimed at inferring strain rates rather than absolute strains.
The goal of optical elastography is to noninvasively, or minimally invasively, quantify meaningful mechanical constants of tissues in a manner that provides clinically relevant information. It is well known that many disease processes, such as tumors of the breast and prostate, manifest themselves as stiff, hard nodules relative to the surrounding tissue. This feature frequently allows for their detection through manual palpation. However, soft tissue palpation is not only qualitative, but also highly subjective. Furthermore, it provides information on a relatively large spatial scale. Thus the ability to detect small tumors by palpation is limited. Manual palpation is also limited to those areas of the body which are accessible to touch. Optical elastography offers the potential for increased spatial resolution (i.e., smaller lesions may be detected) and better strain resolution (low-contrast elastic modulus distributions may be visualized) than other elastographic methods, such as vibration amplitude sonoelastography,6 vibration phase gradient sonoelastography,6 elastography,7 spectral tissue strain measurement,8 and various methods employing magnetic resonance imaging data.9 Strain resolution on the order of single microstrain has been evaluated in biological tissue optically.10 However, this increased spatial and strain resolution is gained at the cost of decreased probing depth as coherent optical methods are limited to the outer few millimeters of tissue. Furthermore, most optical elastography methods are limited in that only relatively small areas or volumes of tissues may be probed at any one time. Nevertheless, optical methods can still be useful in the early detection of neoplastic changes because many of these early changes occur in the mucosa and submucosa of the affected organs. Subsurface skin tumors present themselves as objects with distinct mechanical properties relative to the surrounding normal tissue. The displacement of fibrillar papillary dermis by the softer, cellular mass of a growing melanoma is one such example of this. Optical elastographic techniques may provide a means by which to probe these masses to determine their state of progression and thereby help to determine a proper means of disease management. Other skin afflictions, such as psoriasis and icthyosis, also present as localized tissue areas with distinct mechanical properties that can be delineated optically.
The conventional strategy behind most one- and two-dimensional speckle techniques that are aimed at measuring surface displacements is to compute the correlation between a reference and displaced (sample) speckle patterns. This approach is straightforward and has been used successfully in processing data from speckle strain gauges.11 12 13 This approach, however, is subject to errors and certain limitations for biomedical diagnostics. For example, if the test is subject to vibrations, there will always be some question as to the validity of the reference exposure. Furthermore, if the speckle shifts are substantial, the speckle pattern will decorrelate from the reference and a new reference must be chosen. This can lead to compounding of errors.
However, since the goal of these conventional approaches is simply to quantify a lateral shift in a “noisy” signal, a number of other data collection and data analysis options present themselves. One data collection scheme that has proven to be very useful for evaluating strains in biological tissues begins with the fundamental concept of the laser speckle strain gauge as described by Yamaguchi.11 12 13 The configuration is appropriate for either objective or subjective laser speckle.
The scheme is based upon observing translating laser speckle with a linear array charge coupled device (CCD) camera through an observation angle, θo, as the specimen is sequentially illuminated through two equal, but opposite, illumination angles, ±θs by a collimated laser beam (Figure 1). Alternatively, a two detector, single beam configuration is also possible.13 14 15 Due to the fast decorrelation of the observed speckle patterns as a result of random movements within the tissues (water, blood, etc.), the switching of the illumination angle and subsequent triggering of the CCD array must be rapid enough to acquire at least several records with minimal decorrelation between subsequent exposures. Typically, the switching must be on the order of 50 Hz.10 16 This frequency is beyond the range of most mechanical shutters so electro-optical devices, such as ferroelectric crystals (FLC), in combination with polarizing beam splitters have been used in the past.10 16 In this case, the FLC acts as a binary, switchable half wave plate and the beam is thus transmitted through or reflected by the polarizing beam splitter, depending upon the polarization of the light exiting the FLC.
Using a physical optics approach, Yamaguchi11 has shown that for an object undergoing strain the speckle motion observed through θo for illumination angle θs is given byax is an in-plane motion, az is an out-of-plane motion, εxx is the linear strain in the plane of the detector and laser beams (ultimately, the desired term), Ωy is a rotation about the axis perpendicular to the measurement plane, Ls is the radius of the illuminating wave front (i.e., the source distance), and Lo is the observation distance. By using the illustrated configuration where θo=0°, by using collimated beams (Ls→∞), and by subtracting the speckle motions as observed from the two equal, but opposite, illumination angles, a relation describing the differential speckle motion, δA, can be derived from Eq. (1) 2) would have taken the form az was made negligible to the strain term. This can be accomplished through a proper choice of L0 and θo.
The goal of this strain measurement concept is to determine the shift in a speckle pattern resulting from an applied load. Towards this end, the one-dimensional records are stacked into what is termed a “stacked-speckle history.” 14 15 Stacked speckle histories are time series of one-dimensional views of the speckle patterns combined in a spatio-temporal array such that the spatial dimension (camera pixel) is along the abscissa and the temporal axis is the ordinate. In the configuration shown in Figure 1, two stacked speckle histories are generated, one for each θs. Figure 2 is a gray scale display of one such stacked speckle history taken from an object undergoing a slow linear strain. The figure shows a sequence of 200 one-dimensional speckle patterns stacked one atop another. The desired information, that is, the shift of the speckle pattern, is reflected in the tilt of the structure.
Processing of these stacked speckle data has been accomplished using the transform method.10 15 The method has been adequately described in the literature and therefore will be only summarized here. In the transform method of processing stacked speckle histories, the speckle histories are visually inspected and matching sets of 10–30 records from each history that display minimal decorrelation, or for some other experiment-driven reason, are selected for analysis. These records are transformed into the frequency domain using a fast Fourier transform (FFT) algorithm in the spatial direction and an autoregressive spectral estimator (modified covariance, 3–5 poles) in the temporal direction. Alternatively a two-dimensional FFT may be employed. The result of this operation is a pair of “images” in the frequency domain consisting of a bright band of energy oriented perpendicular to the tilt in the stacked speckle histories. The slopes, m1 and m2, respectively, have the dimensions of length/time and quantify the time rate of speckle pattern shift, δx˙. Thus, by taking the time derivative of Eq. (2), and rearranging to isolate the desired strain term, we get a simple expression for directly estimating the time rate of in-plane strain, ε˙xx (± confidence intervals)N−2° of freedom at a probability level of α. Absolute strains can be determined by an integration over the time course of the experiment.
There are numerous alternative approaches that can be employed for determining speckle movement. Generally, these approaches can be categorized as either nonparametric or parametric approaches. The fundamental difference between the two types of approaches is that nonparametric estimators make no a priori assumptions regarding the nature of the speckle shift. An example of a nonparametric scheme is the correlation approach.17 In this approach, the cross correlation of sequential records is calculated and the shift of the correlation peak determines the amount of speckle pattern shift. Examples of parametric approaches include: a minimum mean square error estimator, a maximum likelihood estimator, and a maximum entropy estimator. Each of these processing schemes is discussed in the following sections.
Nonparametric Speckle Shift Estimators
The traditional approach to processing speckle data makes use of the cross correlation between successive records. Regardless of the dimensionality of the data the concept is the same; we therefore illustrate with data in one dimension. Two sequential speckle records are modeled asδx. The cross correlation between these two records is given by τ=δx. 18 This is the basic concept; we seek the lag (τ) for which the correlation is maximum. It is usually implemented with the FFT19 through use of the (auto) correlation theorem F(F−1).
One of the fundamental assumptions in this (and other) approaches is that the two speckle records are simple translates of one another. With analogy to one of the basic tools of the theory of propagation through atmospheric turbulence, Taylor’s frozen turbulence hypothesis,20 this may be called the “frozen speckle” model.
Parametric Speckle Shift Estimators
Parametric estimators do incorporate a priori knowledge of the experiment. The objective, of course, in using a such model-based estimators is to gain a degree of sensitivity. By making use of information that is known about the conditions of the measurement, nonphysical outcomes can be eliminated. While this often results in greater resolution, this increased resolution sometimes comes with a price. Typically nonparametric approaches are more robust. Because the parametric approaches implicitly eliminate estimates that are “nonphysical,” when data are encountered that are less than ideal, these estimators sometimes break down.
In the material that follows, the relative merits of a couple of specific parametric approaches to estimation of speckle pattern shift or movement are explored. The first is a minimum mean square error (MMSE) technique that, for small amounts of speckle motion, requires no search for an optimum solution. Additionally, when the noise component in the data is Gaussian, it is shown that this approach leads to the maximum likelihood estimate of the speckle pattern shift. Another approach makes use of maximum entropy concepts. This method, like the maximum likelihood technique, is shown to have the ability to estimate subpixel motion with temporal resolution that is much finer than the nonparametric approaches.
Minimum Mean Square Error Estimator
For this discussion of parametric approaches a frozen speckle model is adopted. It is assumed that, over a time on the order of a couple of sequential exposures of the camera, the structure of the speckle pattern is fixed. The only change with time is its lateral motion (Figure 3). Thus the speckle motion can be modeled asδx is small compared to a pixel, Eq. (8) can be approximated as 4) δx that minimizes the error is then determined δx that brings these two records into registration is thus sought. Equation (11) may be solved numerically by making use of a gradient search algorithm.21 If, however, we make use of the approximation in Eq. (9) (small speckle motions), then differentiation with respect to δx and rearranging yields the formula 22 to the derivative δx, is the time rate at which the speckle pattern shifts; units are pixels/record.
Although the means by which the estimate for δx was arrived at was quite specific, this estimation approach is more general than it would seem. For instance, the term approximating the first central difference for the estimate of the temporal derivative arose because the speckle records on either side of the record of interest where chosen for inspection. We could just as easily have included additional records and weighted their contributions appropriately to estimate higher order approximations to the derivative. For example, instead of using the weights23 13) would be 22 This improvement, however, comes at the expense of reduced temporal resolution. Nevertheless, by making use of these higher order approximations, the processing can be tailored to the demands of the experiment.
Up to this point, no specific assumptions about the statistics of the speckle pattern or associated noise have been made. The only a priori knowledge that was introduced is that the speckle shift is small with respect to the pixel size. If the assumption is now made that the measured speckle signal comprises a deterministic speckle signal(!) plus noise, the measured signal can be modeled asjth record, we get 24) is commonly referred to as the likelihood function.24 To choose the δx that maximizes this likelihood we set 12). It is straightforward to show that this is an unbiased estimate of the speckle pattern shift and that the variance of the estimate attains the Crame´r–Rao lower bound.25
Maximum Entropy Estimator
Another parametric approach to calculating speckle motion relies on the concept of entropy. As before, it is assumed that, over a time on the order of a couple of sequential exposures of the camera, the speckle pattern is fixed. Thus the speckle motion is modeled asj+1 is that the latter is shifted by the amount δx. The procedure for estimating the shift parameter is to shift a pair of records in opposite directions and form the quotient of these shifted records 24 in this ratio is then calculated 29), we thus have 28), the shift parameter that most nearly shifts the numerator and denominator of Eq. (27) into registration is chosen. This is the shift that produces a maximally flat quotient. In the event that the distribution, p, is perfectly flat, the associated entropy is simply log(N). A simple gradient search algorithm is all that is required to accomplish our objective.
From the proceeding discussion, it would appear that the original speckle pattern must be interpolated to determine the shifted pattern. While this is true, it can be accomplished very simply by use of the Fourier shift theorem.26 Specifically, making use of the fact that27) is generated by δx that makes this ratio the “flattest.” In the actual implementation, the innermost transforms are computed only once, so that at each iteration of the gradient search, we need only compute 26) is only approximate, there are occasional noisy spikes in the ratio calculated as in Eq. (27). To avoid the possibility of the gradient search algorithm getting trapped in a local maximum, a median filter is applied to the ratio prior to normalization and computation of the entropy. This ensures that the entropy surface is convex. The gradient algorithm can also be aided if one has some a priori knowledge of the amount of speckle shift. In this case the search can be restricted to a limited region. Of course, the objective can also be reached by searching for the minimum of the negative entropy.
Finally, we note that the MMSE solution (Eq. 11) can also be made more efficient using shift theorem and a gradient search algorithm. In this case we need not restrict our interest to small speckle motions.
In evaluating the relative merits of these various estimators there are several criteria. One, of course, is the sensitivity, i.e., the ability to estimate very small speckle motions. On the other hand, the technique should be capable of dealing with very large speckle motions as well. Another evaluation criterion is robustness. This is the ability of the algorithm to cope with other than ideal data.
The first example is for relatively large speckle motions. In this case, a simple nonparametric approach such as correlation works quite well. An example of speckle history for large motion is shown in Figure 5. The speckle motions predicted using a correlation approach are shown in Figure 6. For this calculation, a central differencing technique was used; records to either side of the record of interest were cross correlated. Of course, the speckle shift amount per record is half the total, so that the minimum speckle shift is one half pixel. Speckle shifts estimated using the maximum likelihood estimator (assuming a small speckle shift) are illustrated in Figure 7. This result shows some differences with those derived using correlation processing. On the other hand, the maximum likelihood estimate for arbitrary speckle shift (using a gradient search algorithm to determine the shift), Figure 8, is in excellent agreement with the correlation results. The correlation between these two results is in excess of 0.99. Obviously the small-shift requirement is violated by these data. These results simply show that estimating speckle shifts for large motions is a relatively straightforward problem.
Now an example displaying small speckle motions (Figure 9) is examined. Correlation and max likelihood processing produce the results shown, respectively, in Figures 10 and 11. Clearly, the speckle motion in these data is below the resolution of the correlation processing algorithm (it cannot perceive motions of less than a half pixel), while the max likelihood estimator performs quite well. Processing results using the max entropy estimator are shown in Figure 12. These results are similar to the max likelihood estimate, although noisier and of slightly smaller magnitude.
In the discussion so far it has been assumed implicitly that objective (nonimaged) speckle is being considered. In this case a particular portion of the speckle record cannot be identified with a particular spot on the object. Most of these techniques, however, are appropriate for subjective speckle as well. The example previously discussed was subjective speckle. However, the experiment was carefully contrived so that the specimen (a cortical bone sample) was uniformly illuminated, and uniformly strained. Success of this effort is confirmed by the plots of speckle motion shown in Figure 13. These two traces represent estimates of speckle motion made with the max likelihood estimator using the first 128 columns of the speckle history and the last 128 columns. Although there are some minor differences, which we attribute to the fact that these are two distinct speckle realizations, the speckle motions are very similar (r=0.97). With this concept in mind, it is easy to see that these algorithms can be used for a full two-dimensional sequence of images by separately processing a two-dimensional speckle history for each of the image dimensions as in Figure 14.
Some illustrative parametric processing algorithms have been discussed. Using the basic approach, one could imagine a variety of others. For example, the conventional processing approach searches for the lag for which the cross correlation is maximum. Rather than using the FFT, one might search directly for the shift that maximizes the cross correlation. Specifically, we could seek the speckle pattern shift, δx, that maximizes the cross correlation coefficientδx, an approximate closed form expression for the shift may be obtained readily. Otherwise, a straightforward gradient search algorithm could be used to determine a solution.
For any of the differential speckle measurement techniques, whether it be single detector and two laser beams or single laser beam and two detectors, one arrives at a pair of stacked speckle histories. These histories display the time rate of movement of the individual speckles whether they are subjective or objective. The relationship between these speckle motions and the strain undergone by the test specimen is dictated by the actual measurement configuration. As an example, consider the configuration shown in Figure 1. If δx1 and δx2 represent the speckle motions derived from the two speckle histories (one for laser beam 1 and the other for laser beam 2), then the object strain in the plane of the detector and beams is given byLo is the effective object difference and θs is the illumination angle. For objective speckle, the effective object distance is the physical distance between the object and the detector focal plane. In the case of subjective speckle, it is the misfocus distance. For the complementary configuration of one normally incident laser beam and two cameras, the corresponding relationship is θo is the observation angle. Note that these are the same results as Eqs. (2 3 4).
This work was funded through Research Grant Nos. BES-9807497 and BES-0086719 from the National Science Foundation.
Portions of this manuscript will also appear in S. J. Kirkpatrick and D. D. Duncan, “Optical assessment of tissue mechanics” in Handbook of Optical Biomedical Diagnostics Imaging, V. V. Tuchin, Ed., SPIE, Bellingham, WA (in press). Address all correspondence to Sean J. Kirkpatrick. Tel: 503-216-5246; Fax: 503-216-2422; E-mail: email@example.com