1.

Introduction

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.⁵ 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,⁶ vibration phase gradient sonoelastography,⁶ elastography,⁷ spectral tissue strain measurement,⁸ and various methods employing magnetic resonance imaging data.⁹ Strain resolution on the order of single microstrain has been evaluated in biological tissue optically.¹⁰ 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.

Figure 1

Optical arrangement for collecting speckle data.

Using a physical optics approach, Yamaguchi¹¹ has shown that for an object undergoing strain the speckle motion observed through θ_o for illumination angle θ_s is given by

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.

Figure 2

Stacked speckle history of a sample undergoing a linear strain. Time is given by the y axis and space is given by the x axis.

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, m₁ and m₂, 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)

2.

Methods

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.¹⁷ 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.

2.2.

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.

2.2.1.

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

Eq. (8)

$$g_{j+1}({\mathit{x}}_1)={\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}}-\delta x)\text{,}$$

where the subscript i denotes the pixel (spatial dimension) and the subscript j represents the record (temporal dimension). Assuming that the shift δx is small compared to a pixel, Eq. (8) can be approximated as

Eq. (9)

$$g_{j+1}({\mathit{x}}_{\mathit{i}})\approx {\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})+\delta {\mathit{xg}}_{\mathit{j}}^{\prime }({\mathit{x}}_{\mathit{i}})\text{.}$$

This is simply the first two terms of the Taylor series expansion for g. To introduce a degree of symmetry into the problem, the two speckle records on either side of the record of interest are inspected (Figure 4)

Eq. (10)

$$[g_{j+1}({\mathit{x}}_{\mathit{i}}),g_{j-1}({\mathit{x}}_{\mathit{i}})]\text{.}$$

The δx that minimizes the error is then determined

Eq. (11)

$${\mathit{\epsilon }}_{\mathit{j}}^2=\sum _{i=1}^N{[g_{j+1}({\mathit{x}}_{\mathit{i}}+\delta x)-g_{j-1}({\mathit{x}}_{\mathit{i}}-\delta x)]}^2\text{,}$$

where the sum is over all pixels in the array. The δ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.²¹ 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

Eq. (12)

$$\delta {\mathit{x}}_{\mathit{j}}=\frac{-{\sum }_{j=1}^N[g_{j+1}({\mathit{x}}_{\mathit{i}})-g_{j-1}({\mathit{x}}_{\mathit{i}})][g_{j+1}^{\prime }({\mathit{x}}_{\mathit{i}})+g_{j-1}^{\prime }({\mathit{x}}_{\mathit{i}})]}{{\sum }_{j\mathrm{=}1}^N{[g_{j+1}^{\prime }({\mathit{x}}_{\mathit{i}})+g_{j-1}^{\prime }({\mathit{x}}_{\mathit{i}})]}^2}\text{.}$$

The term in the first square bracket in the numerator is simply the first central difference approximation²² to the derivative

Eq. (13)

$$\frac{\mathit{\partial }{\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})}{\mathit{\partial }{\mathit{t}}_{\mathit{j}}}\approx \frac{g_{j+1}({\mathit{x}}_{\mathit{i}})-g_{j-1}({\mathit{x}}_{\mathit{i}})}{2}\text{.}$$

The spatial derivatives may be approximated similarly

Eq. (14)

$$\frac{\mathit{\partial }{g}_{j+1}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\approx \frac{g_{j+1}(x_{i+1})-g_{j+1}(x_{i-1})}{2}\text{,}$$

Eq. (14)

$$\frac{\mathit{\partial }{g}_{j-1}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\approx \frac{g_{j-1}(x_{i+1})-{\mathit{g}}_{j-1}(x_{i-1})}{2}\text{.}$$

Note that the shift parameter, δx, is the time rate at which the speckle pattern shifts; units are pixels/record.

Figure 3

Illustration of shifting frozen speckle pattern.

Figure 4

Central differencing technique.

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 weights

Eq. (15)

$$\frac{1}{2}\left[-1,0+1\right],$$

we could use²³

Eq. (16)

$$\frac{1}{12}[1,-8,0,8,-1]\text{.}$$

In this case, the formulation for the mean square error would be

Eq. (17)

$${\in }_j^2={\displaystyle \sum _{i=1}^N[g_{j+2}(x_i+2\delta x)-8g_{j+1}(x_i+\delta x)}+8g_{j-1}(x_i-\delta x)-g_{j-2}(x_i-2\delta x){]}^2,$$

the temporal derivative term of Eq. (13) would be

Eq. (18)

$$[-g_{j-2}({\mathit{x}}_{\mathit{i}})+8g_{j-1}({\mathit{x}}_{\mathit{i}})-8g_{j+1}({\mathit{x}}_{\mathit{i}})+g_{j+2}({\mathit{x}}_{\mathit{i}})]\approx 12\frac{\mathit{\partial }{\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})}{\mathit{\partial }{\mathit{t}}_{\mathit{j}}}$$

and the term involving the spatial derivatives would be

Eq. (19)

$$[-2g_{j-2}^{\prime }({\mathit{x}}_{\mathit{i}})+8g_{j-1}^{\prime }({\mathit{x}}_{\mathit{i}})+8g_{j+1}^{\prime }({\mathit{x}}_{\mathit{i}})-2g_{j+2}^{\prime }({\mathit{x}}_{\mathit{i}})]\text{.}$$

Note that this is simply the weighted average

Eq. (20)

$$12\{\frac{4}{3}\left[\frac{g_{j-1}^{\prime }+g_{j+1}^{\prime }}{2}\right]-\frac{1}{3}\left[\frac{g_{j-2}^{\prime }+g_{j+2}^{\prime }}{2}\right]\}$$

of the terms on either side of the record of interest. Further, these higher order operators for the derivative can be used to estimate the spatial derivatives

Eq. (21)

$$\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\hfill & \approx \hfill & \frac{1}{12}[-{\mathit{g}}_{\mathit{j}}({\mathit{x}}_{i-2})+8{\mathit{g}}_{\mathit{j}}({\mathit{x}}_{i-1})-8{\mathit{g}}_{\mathit{j}}({\mathit{x}}_{i+1})\hfill \\ \hfill & \hfill & +{\mathit{g}}_{\mathit{j}}({\mathit{x}}_{i+2})]\text{.}\hfill \end{array}$$

These higher order approximations to the derivative have somewhat better noise characteristics.²² 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 as

Eq. (22)

$${\mathit{d}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})={\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})+{\mathit{n}}_{\mathit{i}}({\mathit{x}}_{\mathit{i}})\text{,}$$

where n is a zero-mean noise signal. If we calculate the central difference about the jth record, we get

Eq. (23)

$$\begin{array}{c}{\mathit{d}}_{j+1}({\mathit{x}}_{\mathit{i}}+\delta x)-{\mathit{d}}_{j-1}({\mathit{x}}_{\mathit{i}}-\delta x)\hfill \\ =g_{j+1}({\mathit{x}}_{\mathit{i}}+\delta x)-g_{j-1}({\mathit{x}}_{\mathit{i}}-\delta x)+{\mathit{n}}_{j+1}({\mathit{x}}_{\mathit{i}}+\delta x)\hfill \\ \hspace{1em}-{\mathit{n}}_{j-1}({\mathit{x}}_{\mathit{i}}-\delta x)\text{.}\hfill \end{array}$$

It is assumed that the sequential speckle patterns are constant mean. Further it is assumed that the noise is statistically independent Gaussian with zero mean and constant variance. As a result the noise probability density function can be written as

Eq. (24)

$$\begin{array}{l}p(d_j|g_j)=p(n_j)\\ =C \text{exp}\left\{\frac{-1}{2{\sigma }^2}{\displaystyle \sum _{i=1}^N[g_{j+1}(x_i+\delta x)]}-g_{j-1}(x_i-\delta x){]}^2\right\},\end{array}$$

where C is a constant. Equation (24) is commonly referred to as the likelihood function.²⁴ To choose the δx that maximizes this likelihood we set

Eq. (25)

$$\frac{\mathit{\partial }}{\mathit{\partial }(\delta x)}\text{ln}p({\mathit{d}}_{\mathit{j}}|{\mathit{g}}_{\mathit{j}})=0\text{.}$$

Carrying out this calculation leads to the formula in Eq. (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.²⁵

2.2.2.

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 as

Eq. (26)

$$g_{j+1}({\mathit{x}}_{\mathit{i}})={\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}}-\delta x)\text{,}$$

where the index j denotes the record number, and the index i denotes the pixel number. This model suggests that the only difference between records j and j+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

Eq. (27)

$${\mathit{R}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})=\frac{g_{j-1}({\mathit{x}}_{\mathit{i}}-\delta x)}{g_{j+1}({\mathit{x}}_{\mathit{i}}+\delta x)}\text{.}$$

The entropy²⁴ in this ratio is then calculated

Eq. (28)

$${\mathit{H}}_{\mathit{j}}(\delta x)=-\sum _{i=1}^N{\mathit{p}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})\text{log}[{\mathit{p}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})]\text{,}$$

where

Eq. (29)

$${\mathit{p}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})=\frac{{\mathit{R}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})}{{\sum }_{i=1}^N{\mathit{R}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})}\text{.}$$

From the normalization of Eq. (29), we thus have

Eq. (30)

$$\sum _{i=1}^N{\mathit{p}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})=1\text{.}$$

By maximizing the entropy in Eq. (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.²⁶ Specifically, making use of the fact that

Eq. (31)

$${\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}}-\delta x)={\mathit{F}}^{-1}\{{\mathit{e}}^{-j2\mathrm{\pi \delta }{\mathit{xf}}_{\mathit{k}}}{\mathit{G}}_{\mathit{j}}({\mathit{f}}_{\mathit{k}})\}\text{,}$$

where

Eq. (32)

$${\mathit{G}}_{\mathit{j}}({\mathit{f}}_{\mathit{k}})=\mathit{F}\{{\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})\}\text{.}$$

With this formulation, the ratio of Eq. (27) is generated by

Eq. (33)

$${\mathit{R}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})=\frac{{\mathit{F}}^{-1}\{{\mathit{e}}^{-j2\mathrm{\pi \delta }{\mathit{xf}}_{\mathit{k}}}F\{g_{j-1}({\mathit{x}}_{\mathit{i}})\}\}}{{\mathit{F}}^{-1}\{{\mathit{e}}^{+j2\mathrm{\pi \delta }{\mathit{xf}}_{\mathit{k}}}F\{g_{j+1}({\mathit{x}}_{\mathit{i}})\}\}}\text{.}$$

The approach outlined above determines the δ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

Eq. (34)

$${\mathit{R}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})=\frac{{\mathit{F}}^{-1}\{{\mathit{e}}^{-j2\mathrm{\pi \delta }{\mathit{xf}}_{\mathit{k}}}{G}_{j-1}({\mathit{f}}_{\mathit{k}})\}}{{\mathit{F}}^{-1}\{{\mathit{e}}^{+j2\mathrm{\pi \delta }{\mathit{xf}}_{\mathit{k}}}{G}_{j+1}({\mathit{f}}_{\mathit{k}})\}}\text{.}$$

Since the model of Eq. (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.

3.

Results

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.

Figure 5

Speckle history displaying large speckle shifts.

Figure 6

Performance of correlation estimator for large speckle shifts.

Figure 7

Performance of closed-form max likelihood estimator for large speckle shifts.

Figure 8

Gradient search for max likelihood solution for large speckle shifts. Compare with results shown in Fig. 13.

Figure 13

Speckle motion estimates based on first and last 128 columns.

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.

Figure 9

Example of speckle history for small shifts.

Figure 10

Correlation estimator results for small speckle motions.

Figure 11

Max likelihood estimate of small speckle motions.

Figure 12

Max entropy estimate for small speckle motions.

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.

Figure 14

Generalization of processing approaches to higher-dimensional problem.

4.

Discussion

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

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 δx₁ and δx₂ 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 by