Analysis of speckle photographs by subtracting phase functions of digital Fourier transforms Part 1

Abstract

This paper presents a method for measuring displacements and strain in digital speckle photography that is an alternative to currently used correlation techniques. The method is analogous to heterodyne speckle photogrammetry wherein optical Fourier transforms are taken of individually recorded specklegrams and combined in a heterodyne interferometer where an electronic phase meter measured the phase differences between the two transforms. Here, digital photographs are recorded and Fourier transformed so that their phase functions can be subtracted and fitted to a linear function of the transform coordinates. The effect of different recording and processing parameters is investigated. It is found that incoherent speckles give better results that those formed by coherent laser light. In addition, image correlation is used to process an identical data set so that comparison of the two methods can be made.

Introduction

Laser speckle photography1 arose as an alternative to holographic interferometry. When an object is illuminated with laser light, the speckles that form in its image move as if attached to the object itself. If an object moves between two exposures of a photograph, the resultant doubling of the speckle pattern can be observed via an optical Fourier transform created by illuminating a small region of the photograph with a narrow, converging laser beam. The doubling of speckles in the photograph gives rise to linear fringes in the transform plane where the beam comes to focus, and because of the similarity to Young’s experiment, they are often called Young’s fringes. The fringes are normal to the direction of the speckle displacement, and their spacing is inversely proportional to its magnitude.

Whereas speckle photographs, or specklegrams, can measure object displacements, they have limitations for measuring strain, which is defined as the change in displacement between two object points divided by their separation. Resistive strain gages can measure strains down to 10-5 over gage lengths as short as one to two mm, and to duplicate this with speckle photography is very difficult. This is especially so in regions of the object where the displacement is so small that no fringes are observable in the transform plane and in regions where the displacement is so large that the fringes are too narrow to observe. To overcome these problems, a technique was developed2-6 for heterodyne readout of halo fringes.

In heterodyne speckle photogrammetry, two photographic glass plates are used to record separate specklegrams before and after a stress is applied to the object. These two plates are placed side by side in an interferometer on a common translation stage and aligned so that the same region on each can be illuminated with each of two small, mutually coherent, converging laser beams. The transmitted beams and scattered halos are combined after equal propagation by a set of mirrors and a beamsplitter. Adjustments are provided so that the halo fringes can be minimized, and the two plates aligned to eliminate relative rotation. The interferometer is provided with the capability of shifting the optical frequency of one beam relative to the other to cause the halo fringes to move and generate sinusoidal irradiance fluctuations. These are detected by an array of photodiodes located in the transform plane, processed by a phase meter, and the phases recorded. As the pair of plates is moved, any translation of one speckle pattern relative to the other causes a change in the relative phases in the fluctuating halo fringe pattern. These changes are recorded and used to calculate the relative displacement of the speckle pattern, which, divided by the amount by which the pair of plates is translated, gives the object strain.

Strain measurement by this technology requires that the object be photographed along the surface normal by means of a tele-centric lens system, i.e., one that has no spherical perspective. The high accuracy available from electronic phase meters, 0.1 degree out of 360, makes it possible to measure strain down to a theoretical level of 10 microstrain (i.e., 10-5) from recordings made with lens systems whose f/numbers are as high as f/10. Furthermore, multiple recordings allow cumulative strain to be measured up to several percent. Although this technique was demonstrated, it has drawbacks. The photographic plates require chemical development and drying before they can be analyzed. The setup and alignment of the interferometer is quite complex and requires considerable space in a darkened laboratory, and the alignment of the photographs is quite critical.

They must be placed so that the same regions are at least approximately illuminated in order to obtain halo fringes at all. Finally, the procedure of recording the phases, moving the pairs of plates, and computing the strains is time consuming. For these reasons, it is desirable to reconsider this process from the point of view of modern digital photography and digital image analysis to see which aspects can be improved and which cannot.

Digital photography has been used for specklegram analysis for quite some time,8,9 however, speckle displacements have been analyzed by image correlation. (See www.lavision.de, and Aramis at www.trilion.com, for example. ) The purpose of this paper is to present and investigate an alternative method of specklegram analysis based upon the subtraction of the phases of digital Fourier transforms in a manner analogous to heterodyne speckle photogrammetry. It begins with a description of the process, followed by the mathematical analysis and experimental study. Several parameters were investigated – recording lens f/number, 8-bit versus 12-bit digitization, and incoherent versus coherent speckle patterns. Finally, a comparison is made to an image correlation measurement of displacement of the best translation data set.

Procedure for digital Fourier processing

The procedure begins by photographing the object with a digital camera. Unless a telecentric lens system is used, any translation of the object toward or away from the object will result in a magnification change, and this will give an apparent strain that must be considered in any subsequent analysis. The optical axis of the lens should be aligned along the surface normal and the camera should have black and white pixels with square pixel format. Twelve-bit image digitization should offer increased accuracy; however, shot noise may be expected to restrict the useful digitization to eight bits. Photographs are captured before and after object perturbation and designated as A and B. For simple displacement analysis, an entire camera image may be used; however, for strain analysis it will be necessary to divide the camera images into segments for separate processing.

The next step is to compute the digital Fourier transforms FAnm and FBmn for each photograph or segment thereof. A digital Fourier transform has the advantage over an optical Fourier transform that it is possible to calculate the phase of the pixels in the transform. The phase values will range randomly from -n to n; however, a speckle displacement will generate a linear phase change across the transform plane whose slope is proportional to the displacement and whose gradient is in the direction of displacement as described by the shift theorem of Fourier transforms. The phase function can be obtained by subtracting the phases of the two transforms; however, simple subtraction will exhibit what may be called random wrapping. This occurs when the value of a pixel phase before displacement is made to exceed -n or n when the additional phase is added to it. For example, if the phase of a pixel is n – A and the phase change is +2A, the resulting phase will be -n+A, and subtracting the two values will give 2A – 2n. Wrapping the phase difference from the range of -2n to 2n into the range of -n to n will remove these random wrapping effects.

Figure 1. 1a. Subtraction of two analytically created random 8-bit phase patterns with a linear phase difference between them. 1b. The same data as 1a after wrapping into 8-bit pixels. Note the wrapping removes the effect of random wrapping  shown in 1a.

After it is unwrapped, the next step is to fit the unwrapped phase difference to a linear function of the spatial frequencies, rax and ray, and these slopes correspond to the x and y translations of the speckles. If strain analysis is being performed, the image will be divided into segments and the slope values from neighboring segments can be used to calculate the average strains experienced between the segments as described below.

Mathematical Analysis

The digital Fourier transform used for this analysis is described by the following equation: where:

f(m,n) is the discrete function whose transform is being calculated, m and n are the pixel indices in the x and y directions, M and N are the number of pixels in the x and y directions of the image or segment thereof, j and k are the horizontal and vertical indices of the coordinates of the transform, and F(j,k) is the Fourier transform of f(m,n) with respect to the variables j and k.

Consistent with the usage in television, we take the y axis, represented by n, to be positive downward. We may calculate the effect of a displacement of one pixel in the x or y directions, dpx and dpy, for the function f(m,n) by substituting m-1 for m or n-1 for n in Eq.(1). This will generate the respective phase functionsin the transform plane, where

The slopes of these phase functions, sj1 and sk1 per incremental values of j or k, are

If dx and dy are the actual displacements of the image pattern and sk and sj are the corresponding measured phase changes per pixel of their transforms, then the corresponding fractions of pixel displacements, dx/dpx and dy/dpy, will equal the corresponding fractional changes in slope of the transform function, sk/sj1 and sk/sk1.

Substituting from Eqs.(3) gives  where the units of sj and sk are radians per pixel in the Fourier transform plane. Equations (5) allow calculation of the object displacement in terms of the pixel spacing, the transform plane slope, and the total number of pixels in the direction considered. The products, dpxM and dpyN, are the physical sizes in x and y of the camera array or the segments of the array. If the object is magnified relative to its image on the camera, then both sides of equations 5 must be multiplied by that magnification to obtain the object displacement.

Strain analysis requires measuring the change in displacement between two points on a surface and dividing that by the distance between the two points. Because, in general, surface strain is expressed as a 2×2 matrix, we need to measure the fractional displacements of four points on the surface. Let these four points be the centers of four neighboring segments which we identify by the subscripts shown below.

11	12
21	22

The average relative displacements of these segments may be defined as:

Given the fact that the M by N pixel segments are separated in x by M pixels and in y by N pixels, the x strain, y strain, and shear are defined as:

When Eqs. (3) and (4) are substituted into Eqs. 6a-6d, and the results substituted into Eqs. 7a-7c, it is seen that the factors of pxM and pyN cancel from the strain calculations. Equations 7a-7c may be rewritten in terms of the measured Fourier transform plane slopes as

For a segment of 256 by 256 pixels, for example, the strain corresponding to a lateral shift of one pixel is (1/256) which equals 3906 microstrain. Measurement of strain to a level of 10-5 would require measurement of subpixel displacement to a level of 1/391 of the pixel spacing.

Experimental Study

An experimental study was carried out to determine the effectiveness of the process described above and investigate the effect of some experimental parameters. Images were obtained via a Prosilica EC650 monochrome TV camera with a 1/3 inch format sensor (4.7 mm horizontal by 3.5 mm vertical) with a pixel array of 659×493 elements on 7.4 |im centers. Image capture was done via two separate programs: the HoloFringe300K program that provided 8-bit images, 640×480 pixels, in a raw data format, and the Prosilica Viewer program that provided 12-bit images, 659×493, in tiff format. The object, mounted on a translation table, was a flat aluminum bar whose visible surface was whitened to provide a flat, diffuse reflection. The object was translated laterally by amounts that were read from the dial of the micrometer on the translation stage, single divisions of which corresponded to 10 |im each. A 25 mm lens was used on the camera which was set with its entrance pupil 355 mm from the object, which resulted in a demagnification of 12.9. The lens aperture was set to f/10, for which the characteristic speckle size should be 7.6 |im, approximately the size of the pixel spacing.

Recordings were made of the object at displacements of 0, -10, -20, -40, -80, -160, and -320 micrometers, each was transformed and the phase functions of all the transforms calculated. The phase of the transform at 0 micrometers was subtracted from the phase of each of the other the transforms and the results wrapped into eight bits obtain the wrapped phase differences. These operations were done using a program called DADiSP, which provides a number of useful functions in addition to the 2DFFT computation.These images were unwrapped using the method of calculated wrap regions from the HoloFringe300K electronic holography program7. Figure 2 shows, for the specklegrams at 0 and 320 micrometers, the wrapped and unwrapped phase difference. Note that the quality of the data is poorer at the left and right edges, which leads to errors in the unwrapping for those regions. To eliminate those errors, the unwrapped phase data was win- dowed to remove the left, right, top, and bottom 25% of the pixels. The corners of the window were located at, xmin=160, xmax=480, ymm=120, and ymax=360.

Figures 2. 2a shows the wrapped phase difference between the transforms of the specklegram recorded at 0 micrometers and the one recorded at 320 micrometers, and 2b shows the data of 2a unwrapped.

A linear function in x and y was then fitted to each data image for least square error by means of a program written in Liberty Basic. The resulting slope values in rax and ray, that is, sj and sk, were substituted into Eqs. 5a and 5b to obtain the x and y image displacements at the camera detector with M taken as 640 and N as 480. The image displacements were then multiplied by the magnification, 12.9, to get displacements of the object. Table 1 presents the measured displacements beside the displacements read from the micrometer dial.

x trans.m	x meas. m	y meas. m
-10.0	-10.0	2.57
-20.0	-23.7	7.68
-40.0	-49.3	8.62
-80.0	-90.6	8.06
-160	-166	7.03
-320	-324	6.85

Table 1. Tabulated values of measured displacement versus micrometer dial readings.

The measured displacements in the y direction should be zero, and it is not clear why they have a nearly constant value. It is possible that the translation stage had some vertical displacement associated with the first 20 |im of horizontal travel. The RMS error for the x displacements is 5.75 |im and 7.09 |im for the y displacements.

The role of lens aperture was investigated by recording another set of translations, 0 to +640 mm, with the lens set to f/2.8. This resulted in a noticeably different image, with a pixel histogram that looked much more like a Gaussian distribution. The results are tabulated in Table 2. The RMS error in the x measurement is 16.8 |m and for the y measurement is 2.33 |m. Clearly, the increase in lens aperture has not helped the measurement accuracy.

x trans. m	x meas. m	y meas. m
10.0	5.04	4.51
20.0	30.0	-1.86
40.0	57.5	-2.64
80.0	96.6	-3.02
160	181	-2.69
320	343	-2.37
640	667	-4.68

Table 2. Translations measured with the lens set to f2.8.

It was also of interest to learn if increasing the digitization to 12 bits would improve the measurement accuracy. Data was acquired via the Prosilica Viewer program for the same translations, the lens set to f/11 and processed in the same way. The results are presented in Table 3, and these show poorer results than with the 8-bit digitization; the RMS error for the x measurement is 11.9 |m and for the y measurement is 1.14 |m. Because the detectors are small, it is expected that their performance is limited by shot noise so that the difference between 8-bit and 12-bit digitization should not be significant. Why the RMS error for the x measurement is approximately twice that found with 8-bit digitization is not clear. It may be an artifact of the image capture programs which are different for the two cases.

x trans.m	x meas.m	y meas. m
10.0	16.6	1.42
20.0	316	0.494
40.0	56.4	-0.0219
80.0	95.8	-0.0142
160	172	-1.84
320	324	-1.45
640	527	-12.8

Table 3. Translations measured with the camera lens set to f/11 and the data digitized to 12 bits.

Live observation of the speckle patterns showed that, for continuous translation, the shifted patterns exhibit a periodic change of the speckles themselves. CCD detector arrays, with interline transfer, typically have active sensors that occupy only about 35% of the area of the array itself. Thus, the speckles integrated by the detectors vary periodically with the translation of the object. Also, the object is illuminated spherically, and this causes a shift of the field that actually passes through the entrance aperture of the lens so that the speckles decorrelate slightly with translation. To get a comparison, an incoherent speckle pattern was generated with the DADiSP program as an array of random numbers displayed as pixels. This was printed on a sheet of paper and cemented to the object surface. The recorded speckle pattern and its histogram of pixel values looked very much like those for laser speckles.

Images were recorded of this incoherent speckle pattern at the positions of the translation stage that were previously used. The camera lens was set to f/5.8, and only 8 bit images were captured based upon the evidence that 12 bit digitization did not improve the results. The results are shown in Table 4.

x trans. m	x meas. m	y meas. m
10.0	12.54	7.940
20.0	21.79	13.51
40.0	40.90	14.51
80.0	80.17	16.31
160	161.4	20.89
320	319.7	15.82
640	639.6	15.31

Table 4. Translations measured using an incoherent speckle pattern.

The results obtained using an incoherent speckle pattern are considerably more accurate than any of those obtained using laser speckles. The RMS error for the x measurements is 1.34 |im and for the y measurements is 15.3 |im. Observation of the images obtained by this method showed that the speckles translated with much less change in their pattern than was observed with the laser speckles, which is consistent with the improved measurement results. The large amount of measured displacement in the y direction; however, indicates that the stage has a significant amount of vertical displacement associated with its horizontal travel. For comparison, another stage was used made by New Focus, their model 8095, and the results are presented in Table 5.

These results show an RMS error for the x and y displacements of 1.97 |im and 3.22 |im respectively. Clearly, this stage moves with less vertical displacement that the previously used one. It is interesting to note that there is an almost linear increase in vertical displacement for the y displacement as a function of the quadratic increase in x displacement.

x trans. m	x meas. m	y meas. m
8.00	9.339	-1.043
16.0	16.27	-0.1983
32.0	34.14	1.619
64.0	66.12	3.152
128	131.6	3.351
256	257.2	4.044
472	473.3	4.941

Table 5. Displacement for incoherent speckles with the New Focus 8095 stage.

Strain was simulated by reorienting the translation stage so that the surface moved toward the camera. An apparent strain of 2817E-6 was generated between two recordings by moving the object toward the camera by 1 mm at a distance of 355 mm. The camera lens was set to f/5.8, and 8-bit digitization was used. The two recordings were divided into 12 segments of 160×160 pixels each, four horizontal by three vertical. Each segment in one recording was processed with its corresponding segment in the other recording to obtain the slope of the unwrapped phase difference for that segment pair. Sets of slope values for four adjacent segments were used in Eqs 8a and 8b to obtain six values of strain, three horizontal by 2 vertical. Table 6 presents the results.

Simulated	x strain	x strain	x strain
2.817E-03	2.829E-03	2.776E-03	2.845E-03
	2.792E-03	2.749E-03	2.802E-03
Simulated	y strain	y strain	y strain
2.817E-03	2.765E-03	2.815E-03	2.837E-03
	2.765E-03	2.754E-03	2.789E-03

Table 6. Strain measurements using digital FFT specklegram analysis.

All six results for x strain and all six for y strain should be equal to 2.817E-3. The average of the six x values is 2.799E-03 and 2.787E-03 for the y values which differ from the correct value by 18E-6 and 30E-6 respectively. The standard deviations are 32E-6 and 30E-6 respectively.