Open Access
Add to BibSonomyAbstract
Presented herein is a performance analysis of a maximum likelihood estimator for calculating small speckle motions. Such estimators are important in a variety of speckle techniques used in non-destructive evaluation. The analysis characterizes the performance (bias and RMS deviation) of the estimator as a function of the signal-to-noise ratio. This SNR parameter is a convenient surrogate for decorrelation of sequential speckle patterns such as are seen in biological tissues. Although the particular estimator is predicated on speckle motions that are a small fraction of a pixel, accurate performance is demonstrated for instantaneous motions of up to ±0.8 pixel/record. Beyond this point, more conventional approaches have been shown to perform well.
©2002 Optical Society of America
1. Introduction
Various algorithms have been presented in the past for estimating the shift in a laser speckle pattern originating from an object (e.g., biological tissue) undergoing a dynamic strain. These algorithms have included the transform method of processing for speckle strain rate measurements [1,2], a maximum likelihood speckle shift estimator, and a maximum entropy estimator [3,4]. In addition, to these parametric estimators, there is the traditional approach to processing speckle data that makes use of the cross correlation between successive speckle records [3,4]. Experimental factors dictate which speckle tracking algorithm is best suited for a given situation. For example, for speckle motions on the order of fractions of a pixel between successive records, the cross correlation approach lacks sufficient resolution to perceive these small shifts. Whereas both the maximum likelihood approach and the maximum entropy approach provide good performance even when the speckle shift is less than a half pixel [3]. For larger speckle shifts, on the order of a pixel or more between successive records, conventional cross correlation approaches offer good performance. The ability to resolve small speckle motions is of importance for high-resolution non-destructive evaluation applications.
Robustness of a speckle shift estimator refers to the quality of its performance in the presence of noise or speckle decorrelation [5]. The purpose of this paper is to evaluate the robustness of a maximum likelihood estimator in the presence of varying amounts of Gaussian white noise. The robustness is characterized in terms of the RMS variation and bias error in estimating small speckle shifts for decreasing signal-to-noise (SNR) ratios in simulated speckle data. The SNR parameter is a convenient surrogate for decorrelation of sequential speckle patterns such as are seen in biological tissues [6].
Conceptual Experiment
The simulated speckle patterns used herein represent data that are acquired using a laser speckle strain gauge [2,7,8–10]. For this speckle strain gauge, the test specimen is sequentially illuminated through two equal, but opposite illumination angles, ±θs, by a coherent optical source. The resulting speckle pattern from each illumination angle is recorded with a linear array CCD camera. These speckle pattern can be either subjective (i.e., imaged) or objective (non-imaged) [11]. As the specimen is dynamically stressed, each successive record from one of the illumination angles is placed as a sequential row in an array so as to generate two (one from each illumination angle), two-dimensional arrays with the units of space (pixels) on the abscissa and time (record number) on the ordinate. As the speckle patterns shift due to the increasing (or decreasing) strain in the sample, the arrays display a ‘tilt’ as the speckles translate in space. We refer to these arrays as ‘stacked speckle histories’ [2]. The crux of the problem is to calculate the shift in the speckle patterns, β, as a function of record number. From this information, the strain, ε, in the sample can be calculated directly as [1–4]
where Lo is the observation distance.
2. Theory
A parametric ‘frozen speckle’ model, much like Taylor’s frozen turbulence hypothesis [11], is adopted. It is assumed that over a time on the order of a couple sequential exposures of the camera, the structure of the speckle pattern is fixed. The only change with time is its lateral
motion (Fig. 1). Thus, the speckle motion can be modeled as
where g represents a pixel gray level, the subscript i denotes the pixel (spatial dimension) and the subscript j signifies the record (temporal dimension). If the shift, β, is small compared to a pixel, Eq. (2) can be approximated as
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 can be employed,
The objective then is to determine the β that minimizes the error,
where the summation is over all pixels in the two records. In other words, the proper β brings these two records into registration. If use is made of the approximation in Eq. (3), then differentiation with respect to β and rearrangement yields the formula
The term in the first square bracket in the numerator is proportional to the first central difference approximation [12] to the derivative; the spatial derivatives (primed variables) may be approximated similarly. Note that the shift parameter, β, is the time rate at which the speckle pattern shifts; units are pixels/record. This estimator was derived in terms of the shift that minimized the mean-square error (MMSE) between adjacent records. Under certain mild constraints, it is also the maximum likelihood (ML) estimate [13] of the shift [3,4].
Although the means of arriving at the estimate for β was quite specific, the estimation approach is more general than it appears. For instance, the term representing the 1st central difference estimate of the temporal derivative arose because the two speckle records adjacent to the record of interest were used. Additional records, appropriately weighted, could be used for higher order approximations to the derivative [3,4]. 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.
Central difference operators as discussed above, are based on the requirement that they be exact for polynomials of specific degree [14]. For example, the operator [-1/2, 0, +1/2] yields an estimate of the derivative that is exact for a polynomial of degree 2 or less; [1/12, -3/12, 0, 3/12, -1/12] is exact for polynomials of degree 4, etc. There are many other possibilities for derivative operators. One that is particularly useful is based on a MMSE approximation to the transfer function of the derivative [12]. Typically, these operators have irrational coefficients such as in the following:
3. Analysis model
Performance of the ML estimator was analyzed using a series of synthetic speckle histories [2]. For a proper speckle pattern as recorded by a pixel array, one must ensure that the pixels sample the speckle pattern adequately. Specifically, to meet the Nyquist criterion, the speckles must be at least twice the pixel dimension. To meet this requirement, a long sequence of samples from an exponential probability distribution (as would be the case for a perfectly polarized speckle pattern [11]), was generated by calculating
where F denotes the Fourier transform and x is an under-filled vector of samples from a distribution that is uniformly distributed on(-π,π). A typical sequence is shown in Fig. 2.
Next, the speckle history was simulated by sub-sampling this sequence by a factor of 10. For the first record the sequence (10, 20, 30, …) was used. For the second record the sequence (10+shift, 20+shift, 30+shift, …) was used, etc. For example, if “shift” is equal to one, the simulated speckle pattern shifts by 1/10 at each record. The nominal shift pattern used is illustrated in Fig. 3, and the resulting speckle history in Fig. 4.
In the analysis, the previously described ML algorithm was used to estimate the shift in a series of simulated speckle histories. To evaluate the susceptibility of the estimator to noise, varying amounts of Gaussian white noise were added to a sequence of noiseless speckle histories. For each case the performance of the estimator was quantified in terms of the bias of the estimate and the RMS variation about the estimate.
Bias error performance was quantified using a sigmoid function;
where α is a slope parameter and SNRo is a threshold level. The actual function fit to the data was
where b is a bias level. Similarly, the RMS performance was assessed by fitting the function
where l is the low-SNR asymptote and h is the high-SNR asymptote. Each of these characterizations was in terms of the steady state strain rate performance. Discontinuities in the instantaneous strain rate were avoided.
4. Results
Performance of a variety of spatial derivative operators was evaluated first. Figures 5 through 8 show estimates of instantaneous speckle shifts for an arbitrarily chosen SNR of 22dB.
Fig. 7. Instantaneous speckle shifts using spatial derivative operator [3, -32, 168, -672, 0, 672, -168, 32, -3]/840 (4th order central difference)
Fig. 8. Instantaneous speckle shifts using spatial derivative operator [0.12019, - 0.74038, 0, 0.74038, - 0.12019] (MMSE optimized derivative operator [12])
This series of figures shows that the algorithm displays much lower bias error using the optimized MMSE derivative operator [12]. As such, the remainder of the simulations used this derivative approximation for calculating the spatial derivatives. The central difference, however, was retained for estimating the temporal derivative.
Shown in Figs.. 9 and 10 are the bias and RMS errors parameterized on SNR. For this portion of the study, a sequence of 20 unique noiseless speckle histories was generated and varying degrees of Gaussian white noise added to each. These figures illustrate the qualitative behavior of all the variations on the fundamental algorithm. The bias for high SNRs approaches a small (generally non-zero) error and at low SNRs approaches an error of 100%. This latter limit makes sense on physical grounds. For large amounts of additive noise, no speckle motion can be perceived; under these conditions the algorithm estimates zero motion (100% error). Similarly, for high SNRs the RMS error approaches some small, non-zero value.
Figure 11 illustrates the effect of increased speckle motion. In this case all conditions were the same as in the previous study except the speckle motion was ±0.2 pixel. Functional dependencies of the bias and RMS error are very similar to the previous case with the exception of increased bias (overestimate for high SNRs). This effect is seen more clearly in Fig. 12 for a speckle motion of ±0.3 pixel/record. Here the overestimate for high SNRs is approximately 4%.
As seen in these last two cases, violating the assumption of small speckle motions leads to overestimates for high SNRs. One way of coping with this problem might be to include additional terms in the Taylor series expansion of the shifted speckle pattern. Unfortunately this approach does not lead to a convenient closed-form estimator. A more straightforward method is to iteratively estimate the speckle shift. The approach is as follows: 1) use the estimator as previously described; 2) use this estimate to remove the shift from the speckle history; 3) re-apply the speckle shift estimator to this residual. The total shift is simply the sum of the two separate estimates. Figures 13 through 17 illustrate this idea. Figure 13 is a typical noiseless speckle history for an instantaneous speckle shift of ±0.3 pixel. Based on the initial estimate of this shift one can interpolate an aligned speckle history as shown in Fig. 14. The dark triangular areas at the left hand side reflect a lack of knowledge of the portion of these records that was shifted into the visible frame. A number of different interpolation routines were tried for this alignment step; a spline-based scheme performed the best. Fig. 15 shows an estimate of the instantaneous speckle shift using the iterative refinement scheme discussed above. The effect is a sharply lower bias level at the expense of slight ringing when the instantaneous speckle motion changes discontinuously. Behavior of the iterative refinement approach is summarized in Figs. 16 and 17. Aside from the lower bias and slightly increased RMS error, the algorithm shows a decreased susceptibility to noise through a lower threshold, SNRo.
Fig 16. Bias performance of iterative refinement algorithm for instantaneous speckle motion of ±0.3 pixel, using temporal derivative kernel [-1, 0, 1]/2 . Note lower SNR threshold.
An alternative method of coping with larger speckle motions is to use a higher order approximation to the central difference for the temporal derivative as discussed in the theory section. In the following results, the temporal derivative kernel [1, -8, 0, 8, -1] was used. Figure 18 shows typical estimates using this kernel for an instantaneous speckle shift of ± 0.3 pixel. Results for various amounts of speckle motion with and without iterative refinement are summarized in Figs. 19 and 20. For all of these simulations the MMSE operator [12] was used for the spatial derivative and the 2nd order central difference operator was used for the temporal derivative. As seen in Fig. 19, the non-iterative algorithm excels for speckle motions of less than ±0.5 pixel/record and the iterative approach for larger speckle motions. The iterative approach displays a bias error that increases only very slowly with increasing speckle motion. Figure 20 shows that both the non-iterative and the iterative approaches display decreasing amounts of RMS noise with increased speckle motion.
Fig. 17. RMS noise performance of iterative refinement algorithm for instantaneous speckle motion of ±0.3 pixel, using temporal derivative kernel [-1, 0, 1]/2
Fig. 18 Instantaneous speckle shift estimates for a ±0.3 pixel shift using temporal derivative kernel [1, - 8, 0, 8, - 1]
5. Discussion and Conclusions
Herein we have demonstrated the performance of a maximum likelihood estimator for determining small speckle motions. The technique estimated instantaneous speckle motion based on records prior to and after the record of interest; cumulative motion was obtained by integrating these motion estimates. Critical parameters that affected the performance of the estimator were the speckle motion per record and the signal to noise ratio of the speckle pattern. Thus, performance was quantified in terms of the bias error and RMS variation about the mean as a function of signal-to-noise-ratio for a variety of speckle motions less than 1 pixel/record. Varying SNRs were contrived by adding different amounts of Gaussian white noise to a series of unique noiseless speckle histories. An exponential distribution (as would be the case for fully polarized speckle) was used to generate these noiseless speckle histories.
Bias as a function of SNR displayed a sigmoidal behavior. For high SNRs the bias approached a small (generally non-zero) error and at low SNRs approached an error of 100%. RMS variation too, displayed a sigmoidal behavior on SNR. For high SNRs it approached some small, non-zero value and gradually approached a higher-level asymptote for low SNRs.
Based on these simulations, the ML estimator was shown to produce accurate predictions for instantaneous speckle motions of approximately ±0.1 pixel/record up to ±0.8 pixel/record. For very small speckle motions the ML estimator slightly underestimated the motion. As actual speckle motion increased, so too did the bias error. Absolute bias in the estimate of the instantaneous motion estimate was under 0.5% for motions as large as 0.4 pixel per record. For speckle motions larger than approximately ±0.5 pixel/record, an iterative refinement approach was shown to be useful. This variation on the basic algorithm produced lower bias errors at the expense of slightly higher noise in the estimate. This noise would be of secondary importance, however, if the feature of interest were the cumulative strain. In this case, variations about the bias level would have little impact. For speckle motions greater than ±0.8 pixel/record, other approaches, e.g., correlation methods [15] are more appropriate.
The crux of making the ML estimator perform accurately was found to be the choice of the derivative kernels for the spatial and temporal derivatives. We found that the 3rd order MMSE optimized derivative kernel [12] worked well for the spatial derivative while a 2nd order central difference operator performed well for the temporal derivative. Use of the MMSE optimized derivative kernel for the temporal derivative displayed no performance gain. This is because speckle motion generally displays a polynomial temporal behavior and the central difference operators, as discussed earlier, are exact for particular polynomial orders. On the other hand, the spatial structure of the speckle pattern is more complex and a better derivative kernel that a central difference is required.
We note that for the case of non-polarized speckle, the noiseless speckle intensity variations display Rayleigh rather than exponential statistics. Tests of the ML algorithm for such Rayleigh distributed speckle patterns displayed very similar results to those cited here. Further, we note that these simulations employed speckle patterns that possessed full floating point precision. Simulations were also performed on speckle patterns that were quantized to 8 bits. No appreciable performance differences were observed.
Acknowledgments
This work was funded in part by NSF Grants #0196172, #0201841 and NIH Grant #1R24EB000224-03.
References and Links
1. D. D. Duncan, F. F. Mark, and L. W. Hunter, “A new speckle technique for noncontact measurement of small creep rates,” Opt. Eng. 31, 1583–1589 (1992). [CrossRef]
2. D. D. Duncan, S. J. Kirkpatrick, F. F. Mark, and L. W. Hunter, “Transform method of processing for speckle strain-rate measurements,” Appl. Opt. 33, 5177–5186 (1994). [CrossRef] [PubMed]
3. D. D. Duncan and S. J. Kirkpatrick, “Processing algorithms for tracking speckle shifts in optical elastography of biological tissues,” J Biomed Opt. 6, 418–426 (2001). [CrossRef] [PubMed]
4. S. J. Kirkpatrick and D. D. Duncan, “Optical Assessment of Tissue Mechanics” in Handbook of Optical Biomedical Diagnostics, V. V. Tuchin, ed. (SPIE Press, Bellingham, WA. 2002).
5. E. E. Konofagou, T. Varghese, J. Ophir, and S. K. Alam, “Power spectral strain estimators in elastography,” Ultrasound Med. Biol. 27, 1115–1129 (1999). [CrossRef]
6. A. Oulamara, G. Tribillon, and J. Duvernoy, “Biological activity measurement on botanical specimen surfaces using a temporal decorrelation effect of laser speckle,” J. Mod. Opt. 36, 165–179 (1989). [CrossRef]
7. I. Yamaguchi, “Speckle displacement and decorrelation in the diffraction and image fields for small object deformation,” Opt. Acta. 28, 1359–1376 (1981). [CrossRef]
8. S. J. Kirkpatrick and D. D. Duncan, “Noncontact microstrain measurements in orthodontic wires,” J. Biomed. Mater. Res. 29, 1437–1442 (1995). [CrossRef] [PubMed]
9. S. J. Kirkpatrick and B. W. Brooks, “Micromechanical behavior of cortical bone as inferred from laser speckle data,” J. Biomed. Mater. Res. 39, 373–379 (1998). [CrossRef] [PubMed]
10. S. J. Kirkpatrick and M. J. Cipolla, “High resolution imaged laser speckle strain gauge for vascular applications,” J Biomed Opt. 5, 62–71 (2000). [CrossRef] [PubMed]
11. J. W. Goodman, “Statistical properties of laser speckles” in Topics in Applied Physics, Vol. 9: Laser Speckle and Related Phenomena, Second Enlarged Edition, J. C. Dainty ed. (Springer Verlag, New York1984).
12. B. Jähne, Practical Handbook on Image Processing for Scientific Applications (CRC Press, Boca Raton1997).
13. R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification, 2nd edition (John Wiley & Sons, Inc. New York2001).
14. C. F. Gerald, Applied Numerical Analysis (Addison-Wesley Publishing Company, Reading1973).
15. T. Takemori, K. Fujita, and I. Yamaguchi, “Resolution improvement in speckle displacement and strain sensor by correlation interpolation,” Laser Interferometry IV: Computer Aided Interferometry, R. J. Pryputniewicz, ed., Proc. SPIE 1553, 137–148 (1990).
