- Email: [email protected]
Dynamic ESPI with subtraction–addition method for obtaining the phase
15 October 2002
Optics Communications 212 (2002) 35–43 www.elsevier.com/locate/optcom
Dynamic ESPI with subtraction–addition method for obtaining the phase Violeta Madjarova *, Satoru Toyooka, Rini Widiastuti, Hirofumi Kadono Department of Environmental Science and Human Engineering, Saitama University, 255 Shimo-Okubo, Saitama City, Japan Received 29 January 2002; received in revised form 2 May 2002; accepted 7 August 2002
Abstract Dynamic electronic speckle pattern interferometry (DESPI) was developed for in situ observations. The quantitative evaluation of the deformation field was performed through 2-D subtraction–addition method (SAM) for phase analysis. This method does not require additional phase modulation, which makes it applicable to studying dynamic events. The method utilizes the ratio of the subtraction and addition correlation fringe patterns. The bias intensity that has to be removed from the addition images was determined by temporal averaging of a sequence of speckle patterns. To reduce the error, Gaussian filter was applied to the correlation fringe patterns. The deformation field was evaluated with accuracy of k=10. Ó 2002 Published by Elsevier Science B.V. Keywords: Dynamic ESPI; Phase analysis; Gaussian filter; In situ observations; Plastic deformation
1. Introduction Optical methods can usually provide a whole field and in some cases a real-time picture of a given process, which makes them valuable means of measurement in different fields of science. For objects with optically rough surfaces, ESPI is one of the most appropriate method for high sensitivity measurements of out-of-plane and in-plane displacement components, strain, and vibrations. In addition, various image-processing techniques
*
Corresponding author. Fax: +81-48-853-3459. E-mail address: [email protected] (V. Madjarova).
could be applied, because the signal is obtained in digital form [1–4]. The deformation field in form of correlation fringes is obtained by comparing two images, one representing the reference object state and the other representing a deformed object state, either by subtracting or adding them. Conventionally in the ESPI method, the reference image is taken once, and the development of the deformation for a specified time can be received, by comparing the sequence of speckle patterns with the reference speckle pattern. In this manner, the technique could only be applied for studying deformation in a limited time, restricted by the deformation speed, CCD camera resolution, and de-correlation due to the speckle size [1]. Rather than keeping the reference image constant, in the
0030-4018/02/$ - see front matter Ó 2002 Published by Elsevier Science B.V. PII: S 0 0 3 0 - 4 0 1 8 ( 0 2 ) 0 1 9 0 9 - 0
36
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
method proposed by Suprapedi and Toyooka [5], the reference speckle pattern is renewed, and the subtraction interval is kept constant. This timesubtraction sequence method can allow the ESPI method to be successfully applied for studying dynamic events throughout their entire development. We call this time-sequence method dynamic ESPI (DESPI). We successfully applied DESPI to investigate dynamic events, such as plastic deformation of metals, in which various interesting phenomena were revealed [5–9]. These phenomena were studied qualitatively and semi-quantitatively, however, for the full understandings of the events quantitative evaluation of the phase is needed. In this paper, we aim to develop a quantitative method to obtain the phase and analyze the deformation field. There are numerous techniques for obtaining the phase that are employed presently. The essential of all these techniques is that additional modulation phase terms that can be controlled experimentally are introduced in time or space domain. Mainly three categories of methods exist: (1) fast Fourier transform (FFT) method with phase carrier proposed by Takeda [1–4,12–14], (2) temporal phase shifting (TPS) [1–4,11,15], and (3) spatial phase shifting [1–4]. In FFT method for correlation fringe analysis, a spatial carrier is introduced in either of the two speckle images. Provided that the frequency of the correlation fringes is lower than the carrier frequency it could be possible to extract the phase by FFT method. Thus, the spatial carrier limits the number of fringes, i.e., the maximum amount of deformation, which can be measured. Instead of applying a spatial carrier into the fringe pattern, Saldner et al. [14] showed that in the case of ESPI, FFT would be more successful if the carrier is applied to the speckle image itself, not to the correlation fringes. There are, though, certain requirements to be met: as the method requires the speckle size to be at least five times the size of the pixel, a camera with either high resolution or high sensitivity is needed. In addition, the use of large speckle sizes reduces the spatial resolution and as it is seen in the experimental results, places with denser fringes cannot be resolved. The FFT method can also be considered in the temporal
domain. Joenathan et al. [12] applied the temporal Fourier transform to obtain the phase. In this method, each pixel was considered independently and the Fourier transform of the temporal variation of the intensity at each pixel was performed. Some limitations and sources of error in the method are discussed in [13]. In TPS method, since the intensity is sampled temporally, the object should be stable during the acquisition of at least three frames. This makes the method difficult to apply in dynamic deformation measurements, where the object is continuously changing its state. The requirement for the object to be stationary might not be necessary if the following assumptions can be satisfied: (a) the applied load and the acquisition sampling rate are such that the time evolution of the deformation can be approximated locally into liner function; (b) the deformation rate does not produce a phase change larger than p (or jp=2j if a temporal carrier is introduced)at any position of the object between two successive frames. Then the phase can be obtained by combining five consecutive frames as proposed by Xavier Colonna de Lega and Pierre Jacquot [11]. The main drawback of the method is the assumption for almost linear deformation of the object. Moreover, if sudden change in the deformation state appears, the phase may vary faster than the value tolerable between five frames. The above-mentioned drawbacks can become significant limitations in real tension experiment, where the plastic deformation is a nonlinear phenomenon and sudden changes in the deformation state appears often. In addition, for very small deformation, the successive intensities values will be more susceptible to the noise due to CCD read out, vibrations of the object, digitization of pixels, and the deviation of the deformation from the linear trend. An alternative to TPS, in the case of dynamic events, could be the spatial phase shifting where the phase modulation is separated out in space (usually using a diffraction element) and recorded at the same time on one frame, containing several images with different phase. The main disadvantage of this technique is its susceptibility to rough environmental conditions and dependence on the CCD camera characteristics [2].
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
Reviewing the existing methods for dynamic events studies, and considering some of the restrictions of these methods, and the specifics of the events we are studying, we choose the method proposed by Yoshida et al. [10]. The reason for choosing this method is that it has no requirements for linearity in deformation or limitation in the maximum deformation (as long as the sampling frequency is satisfied), as is the case in [11]. In addition the method does not require any phase-shifting system, or spatial carrier [10]. In order to obtain the phase, both the correlation fringes from subtraction and addition images are used, and the phase is obtained from the ratio of these two patterns. We named this method subtraction–addition method (SAM). In the initially proposed method, however, there are some problems concerned with the bias intensity and unwrapping procedure. The bias intensity is received, simply by taking single beam illumination from both interfering beams and, adding these two values. The obtained bias intensity is applied for the rest of the experiment, assuming that it does not change. This can lead to large errors when carrying out experiments in which the object under study considerably alters its state. In addition, the proposed unwrapping procedure was applied only to the 1-D case for almost parallel fringes. Therefore, the successful application of this method is possible, if the method of obtaining the bias intensity allows the bias intensity value to be renewed, and if an appropriate 2-D phase unwrapping procedure is employed. In this paper, we present experimental results with an in-plane sensitive DESPI for in situ observation of dynamic events (plastic deformation of aluminum alloy). For obtaining the phase, we developed a 2-D phase unwrapping procedure for the SAM. To overcome the difficulty with the bias intensity, which was outlined above, we propose temporal averaging of sequence of speckle images and renewing the information about the bias intensity for every speckle image. This can simplify the optical set-up, because no additional arrangement, such as shutters or external phase modulator, is needed for the acquisition of the bias intensity.
37
2. Phase evaluation by SAM 2.1. Principle In general, interference speckle images obtained in ESPI can formally be represented by Iðx; y; tÞ ¼ I0 ðx; y; tÞ þ Im ðx; y; tÞ cos ðhðx; y; tÞ þ uðx; y; tÞÞ;
ð1Þ
where I0 ðx; y; tÞ and Im ðx; y; tÞ are the bias intensity and the modulation intensity, respectively. hðx; y; tÞ term is the random phase of the speckle field, which varies very rapidly in space domain. In these variables the time dependence is in general case omitted, since they usually vary significantly in the space domain but very slowly in the time domain. However, for long experiments, where the object alters significantly, this time dependence cannot be omitted. uðx; y; tÞ term is the phase introduced by the deformation of the object under study, which usually varies slowly in space and time domains. Typically in ESPI systems, to produce correlation fringes that correspond to the deformation of the object, the absolute difference between two images is taken. Addition images are not usually used, because of their very low visibility due to the bias intensity term. In the method proposed by Yoshida et al. [10], both the subtraction and addition images are used, and the phase is obtained from the ratio of these two images. The main idea of the method could be described from Eq. (2) to Eq. (4). The subtraction correlation fringe patterns for DESPI are represented by Isub ðx; yÞ ¼ Ii ðx; yÞ Iiþp ðx; yÞ Duðx; yÞ ¼ 2Im ðx; yÞ sin hðx; yÞ þ 2 Duðx; yÞ sin ; 2 i ¼ 1; 2; 3; . . . ; p ¼ 1; 2; 3; . . . ;
ð2Þ
where i indicates the frame number at the time ti , p indicates the subtraction interval, which is determined depending on the experimental conditions, and Duðx; yÞ ¼ uiþp ðx; yÞ ui ðx; yÞ is the phase change between two states of deformation. The
38
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
addition correlation fringe patterns for DESPI after removing the bias intensity are represented by Iadd ðx; yÞ ¼ Ii ðx; yÞ Iiþp ðx; yÞ I0 ðx; yÞ Duðx; yÞ ¼ 2Im ðx; yÞ cos hðx; yÞ þ 2 Duðx; yÞ cos ð3Þ : 2 The value of interest Du can be obtained by taking the arctangent function of the ratio of the subtraction image to the addition image as described by hIsub ðx; yÞi Du ¼ 2 arctan ; ð4Þ hIadd ðx; yÞi where h i denotes spatial averaging to remove the noise term, i.e., the sine and cosine terms containing the hðx; yÞ term.
Fig. 1. The time variation of the intensity at a given pixel. It demonstrates the sinusoidal variation of the intensity and also the slow temporal change in the bias and modulation intensities.
can be performed using a symmetrical time window around each speckle image. 2.3. Phase unwrapping
2.2. Elimination of bias intensity by temporal averaging over sequence of speckle images As was already mentioned in Section 1, the bias intensity also varies slowly in time. In the case when SAM is applied, using a constant value for the bias intensity, which is determined at the beginning of the experiment, will lead to considerable errors. In [10] additional arrangement in the optical set-up with mechanical shutters was employed to obtain the bias intensity. If the information about the bias intensity should be renewed, the experiment has to be interrupted to acquire the bias intensity again. This can be a considerable drawback in dynamic and high-speed event studies. Instead of using additional elements in the setup, we propose an alternative signal processing technique. In Eq. (1), it can easily be seen that, for continuously changing deformation the intensity is a cosine periodic function of time. To illustrate the temporal variation of the intensity, represented by Eq. (1), we show the intensity at a given pixel in Fig. 1, where the sinusoidal variation of the intensity can clearly be distinguished. This periodic temporal variation could allow determining the bias intensity by temporal averaging over a sequence of speckle images. The temporal averaging
The subtraction correlation fringes are obtained by subtracting two frames with a fixed subtraction interval, depending on the experimental conditions. After the bias intensity is calculated, it is subtracted from the addition image, and the addition correlation fringes are obtained. The next step is the calculation of the phase difference Du using Eq. (4). Unlike the other well-established phase calculation methods, the wrapped phase map is not obtained with the distinctive 2p phase jumps typical for the conventional phase-shifting methods. Instead the map is received as a modulus of the conventional wrapped phase map and the distinctive jumps disappear. The reason for this is the use of the modulus values of sinðDu=2Þ in Eq. (2) and cosðDu=2Þ in Eq. (3) and the sign in the wrapped phase values cannot be examined. In this paper, we propose a new 2-D phase unwrapping algorithm for the SAM. The only assumption we make is for monotonically increasing deformation. The main idea of this algorithm is shown in Fig. 2. For simplicity, 1-D linearly increasing noise-free signal is assumed in this graph in a limited interval of phase change from 0 to 4p. The two rectified sinusoidal functions represent the sine and cosine functions cor-
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
39
nðx; yÞ nðx;yÞ1 Duðx; yÞ ¼ 2p Du0 ðx; yÞ þ ð1Þ 2 ðn ¼ 1; 2; 3; . . .Þ;
ð5Þ
where ½
is the Gaussian symbol that stands for an integer number, n is an integer, Du is the unwrapped phase value, and Du0 is the wrapped phase value.
3. Experiments and discussions 3.1. DESPI system for in situ observation
Fig. 2. Illustration of the 1-D unwrapping procedure. The sine and cosine function are given in relative units.
responding to the intensity values of the subtraction Isub and addition Iadd , respectively. The wrapped phase Du calculated from the arctangent of the intensity ratio of the subtraction and addition fringes is shown as the triangle-shaped function, while it takes a saw tooth shape function in the conventional interferometry. Taking into account the assumption we made for the deformation, a solid straight line could be drawn to represents the unwrapped phase value. Although we make an assumption for the object deformation to be monotonically increasing, since the method has a sign ambiguity, there are many practical applications where the assumption is satisfied, such as the tensile experiments described in this study. The first step in the algorithm is to divide the interval into four areas, and assign an index number in each area. The assigned numbers in the areas correspond to the fringe order; therefore the area is indexed as Areas 1–4. The borders of two adjacent areas correspond to the position of the maximum or minimum of the wrapped phase values. These points are determined by scanning the positions of zero gradients from the obtained wrapped phase map expressed by Eq. (4). After correct index numbers have been assigned on the whole area, the phase values at the position ðx; yÞ are unwrapped following the rule:
The method described in Section 2 can be applied to in-plane as well as out-of-plane deformation measurement, because it was considered for the general ESPI. In this paper, we carried out a tensile experiment of metal object, where an inplane sensitive DESPI system was used. The experimental system is illustrated in Fig. 3. The light source was a second harmonic generation of a YAG laser with output power of 200 mW, and wavelength of 532 nm. The light emitted from the laser was expanded by an objective lens L to illuminate the specimen. Here, the upper half part of the expanded beam directly illuminates the object with an angle a, while the lower half part is reflected by the mirror M3 to illuminate the object
Fig. 3. The scheme of the optical set-up for in-plane deformation measurement. M1, M2, and M3 – mirrors, L – objective lens, CCD – camera, DL – data logger, PC1, PC2 – computers, and S – specimen.
40
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
with an angle of a. The propagation vectors of the two branches lie in the plane parallel to the z–y plane and the sensitivity vector lies in the x–y plane. Strictly speaking, the configuration of two divergent beams implies also the presence of sensitivity vector component that lies outside the x–y plane. Taking into account that the magnitude of the sensitivity vector that lies out of x–y plane is very small and that the displacement in the tensile experiment is mainly in-plane, we consider only the y-component of the in-plane displacement. The amount of deformation per unit fringe spacing is given by k=ð2 sin aÞ and was estimated to be 380 nm, with the incident angle a ¼ 45°. The CCD camera was positioned perpendicularly to the surface of the specimen and the interference speckle patterns are taken continuously with a constant acquisition rate. In the experiment, the camera acquisition rate was 2 frames per second (fps). The captured images are subtracted to get correlation fringes with a time interval of 4 s and are displayed on the monitor in real-time. Fifty-six thousand frames were acquired in total during the entire experiment. We performed a tensile experiment of a dumbbell-shaped plate of aluminum alloy A2017 with the effective length of 80 mm, width of 30 mm, and thickness of 5 mm. Before the tensile experiment, the specimens were annealed by exposing in atmosphere of 450 °C for 1 h and slowly cooled to the room temperature. The specimen was stretched on a tensile machine with a constant tensile speed, in a manner that the lower hole of the specimen was caught on a cylinder rod, which was fixed to a cramp grip, and the upper hole was caught on a cylinder, which was fixed on a moving crosshead. Several experiments were performed with different tensile speed of 0.5, 1, and 1.6 lm/s. One of the results from the experiment with the tensile speed of 0.5 lm/s is shown in Fig. 4. Fig. 4 presents the stress curve with the corresponding correlation fringes patterns obtained at different states of the deformation process, which are shown below the stress curve. The stress increases linearly in the elastic deformation state, followed by zigzag variations in the plastic state. The zigzag variation or serration is known as the phenomenon associated with strain aging called Portvin–Le Chartelier
Fig. 4. The stress curve in function of time with the corresponding fringe patterns. They demonstrate the complicated manner in which the deformation develops during the experiment.
effect [16]. The pictures presented demonstrate the development of the plastic deformation in rather complicated manner. Pictures 3 and 4 represent a stable band-like pattern at an angle of approximately 45° to the tensile axis. That pattern suddenly appeared many times during the experiment (when the stress curve changes sharply), and the position where it appeared moved over the specimen upwards or downwards. At the later stage of the plastic deformation state, there were moments at which the angle of the band pattern switches, i.e., the angle changes from approximately +45° to approximately )45° to the tensile axis (pictures 5, 6, and 7). 3.2. Phase analysis by SAM After the frame acquisition, the quantitative analysis based on SAM was performed to obtain the object phase. The procedure included: (1) calculating the bias intensity for every frame and subtracting it from the corresponding frame; (2) calculating the subtraction and addition images and applying filters for removing the speckle noise; (3) calculating the wrapped phase values using addition and subtraction images; and (4) applying the unwrapping procedure. There are two main sources of errors in the SAM with temporal evaluation of bias intensity
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
that reduce the accuracy of the measurement. The first source of errors is the inaccurate determination of bias intensity. Mainly, in real experiments, the accurate determination of the size of the temporal window that should include an integer number of cosine function periods in Eq. (1) is difficult to accomplish automatically, for the whole specimen and for the whole experiment, since it varies in time and space. Considering the fact that the averaging is performed over several intervals of a periodic function, the influence of the fraction of the period over the average value tends to zero with the increasing number of periods inside the window. In our experiment, the temporal averaging window was determined in a manner that it includes at least several periods of the cosine function. Another limitation in the size of the temporal window arises from the fact that the bias intensity is not a constant value. If the temporal window is too wide the bias intensity cannot be regarded as constant inside the window, and the average value received will not be its correct representation. Taking into account the tensile speed and the variation in the bias intensity, it was found that the optimum temporal window is 25 s, which includes on the average more than 6 periods of the cosine function in Eq. (1). The second source of errors comes from the inaccurate determination of the Areas 1–4, etc. Mainly the errors arise when the noise in the images is not filtered out successfully, which will bring inaccurate maximum and minimum positions in the wrapped phase. In order to get a
41
correct result on the phase unwrapping process, application of a filter to minimize the fluctuation due to the speckle noise on the obtained signal becomes inevitable. The noise on the obtained signal is associated with the presence of the high-frequency fluctuation caused by the speckle random phase h, which is uniformly distributed over all values between 0 and 2p. Taking into consideration that h and Du are statistically independent, and the spatial frequency of h is significantly higher than that of Du, a low-pass filter in the spatial domain can be used to minimize noise due to the presence of h. Performing an optimum low-pass filtering, the following approximation of hsinðh þ Du=2Þ sinðDu=2Þi
2 sinðDu=2Þ p
and 2 cosðDu=2Þ p can be achieved with significant accuracy, and consecutively phase value Du can be obtained. The brackets h i denote the spatial filtering operation. In this research, a Gaussian low-pass filter is applied to filter out the high-frequency term due to the random phase h. An important note considered here is that the filter must be adjusted depending on the spatial frequency of the obtained correlation fringes. Fig. 5 gives the results of addition, subtraction, and wrapped phase of a typical band-like pattern that makes angle of approximately 45° to the hcosðh þ Du=2Þ cosðDu=2Þi
Fig. 5. Addition, subtraction, and wrapped phase images with the 3-D plot of the deformation field representing the locally concentrated deformation field. The subtraction interval was 4 s, and temporal averaging window was 25 s.
42
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
Fig. 6. Addition, subtraction, and wrapped phase images with the 3-D plot of the deformation field representing the ‘switching’ of the position of the locally concentrated deformation field from approximately 45° to )45°.
tensile axis, which appeared several times during the experiment. The visibility of the addition correlation fringes is comparable to that of the subtraction correlation images. The 3-D plot of the deformation field after the unwrapping of the phases is also given in this figure. It can clearly be seen that the deformation is concentrated in a very narrow area, which we call locally concentrated strain band. The temporal behavior of this locally concentrated strain band was discussed in more detail in [5,17]. The ‘switching’ of the direction of the band from approximately +45° to )45° to the tensile axis can be seen in Fig. 6, where again addition, subtraction, and wrapped phase images are shown with the corresponding 3-D plot. The amount of deformation after the ‘switching’ is less, but it changes in a pulsating manner as it was discovered and discussed in [5,17]. The interval between the two states is approximately 130 s. The accuracy of the system was evaluated by estimating the standard deviation of the deformation for the area where almost no deformation was detected. It was on the average 50 nm, which corresponds to k=10.
4. Conclusions In this paper, we proposed a dynamic ESPI system for in situ observations and a 2-D subtraction–addition method for evaluating the phase in DESPI correlation images. The bias intensity in the addition speckle correlation images was de-
termined by temporal averaging over a sequence of speckle patterns. The phase map is not obtained with the distinctive 2p phase jumps, but as a modulus of the conventional wrapped phase map and the distinctive jumps disappear. Thus the main limitation of the method is the sign ambiguity and an a priory information for the object under study should be used for the correct unwrapping. However, there are many practical situations, where the proposed method can successfully be applied, as demonstrated in this study. In tensile experiments, the deformation direction is known a priori and is usually monotonically increasing. To perform the phase unwrapping, the folding lines of minimum and maximum of calculated wrapped phase distribution were utilized. Prior to the unwrapping procedure Gaussian filters were applied to reduce the speckle noise. In the present experiment the acquisition rate of the camera is 2 fps, but the system and the method can be applied with highspeed camera, because there is no limitation on the acquisition rate of the camera that is required by the method itself. It was demonstrated that DESPI with the proposed method for obtaining the phase could successfully be applied to dynamic events with an accuracy of approximately k=10. The main objective in developing the DESPI with SAM for obtaining the phase is to expand the possibilities of ESPI in studying dynamic events in rough environmental conditions. The method for receiving the bias intensity allows the application of the system in studying high-speed phenomena. The proposed method is expected to contribute to
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
quantitative analysis of deformation states in different materials, as it is not confined only to tensile experiments.
Acknowledgements This research was partly supported by a Grantin-Aid (13555193) for Science Research of the Ministry of Education, Science, Sport and Culture, Japan.
References [1] P.K. Rastogi, Digital Speckle Pattern Interferometry and Related Techniques, Wiley, Chichester, 2001. [2] D.W. Robinson, C.R. Raed (Eds.), Interferogram Analysis: Digital Fringe Pattern Measurement Techniques, Institute of Physics Publishing, Bristol, 1993.
43
[3] P.K. Rastogi (Ed.), Optical Measurement Techniques and Applications, Artech House, Boston, 1997. [4] P.K. Rastogi (Ed.), Photomechanics, Springer, Berlin, 1999. [5] Suprapedi, S. Toyooka, Opt. Rev. 4 (2) (1997) 284. [6] S. Toyooka, R. Widiastuti, Q. Zhang, H. Kato, Jpn. J. Appl. Phys. 40 (2001) 873. [7] S. Toyooka, X.L. Gong, Jpn. J. Appl. Phys. 34 (1995) L1666. [8] Suprapedi, S. Toyooka, Phys. Mesomech. 1 (1998) 51. [9] X.L. Gong, S. Toyooka, Exp. Mech. 39 (1999) 25. [10] S. Yoshida, Suprapedi, R. Widiastuti, E.T. Triastuti, A. Kusnowo, Opt. Lett. 20 (1995) 755. [11] X.C. de Lega, P. Jacquot, Appl. Opt. 35 (1994) 5115. [12] C. Joenathan, B. Franze, P. Haible, H.J. Tiziani, J. Mod. Opt. 44 (1998) 1975. [13] C. Joenathan, P. Haible, H. Tiziani, Appl. Opt. 38 (1999) 1169. [14] H.O. Saldner, N.E. Molin, K.A. Stetson, Appl. Opt. 35 (1996) 332. [15] J.M. Huntly, G.H. Kaufmann, D. Kerr, Appl. Opt. 38 (1999) 6556. [16] G.E. Dieter, Mechanical Metallurgy, McGraw-Hill, London, 1988. [17] S. Yoshida, S. Toyooka, J. Phys. 13 (2001) 6741.
Optics Communications 212 (2002) 35–43 www.elsevier.com/locate/optcom
Dynamic ESPI with subtraction–addition method for obtaining the phase Violeta Madjarova *, Satoru Toyooka, Rini Widiastuti, Hirofumi Kadono Department of Environmental Science and Human Engineering, Saitama University, 255 Shimo-Okubo, Saitama City, Japan Received 29 January 2002; received in revised form 2 May 2002; accepted 7 August 2002
Abstract Dynamic electronic speckle pattern interferometry (DESPI) was developed for in situ observations. The quantitative evaluation of the deformation field was performed through 2-D subtraction–addition method (SAM) for phase analysis. This method does not require additional phase modulation, which makes it applicable to studying dynamic events. The method utilizes the ratio of the subtraction and addition correlation fringe patterns. The bias intensity that has to be removed from the addition images was determined by temporal averaging of a sequence of speckle patterns. To reduce the error, Gaussian filter was applied to the correlation fringe patterns. The deformation field was evaluated with accuracy of k=10. Ó 2002 Published by Elsevier Science B.V. Keywords: Dynamic ESPI; Phase analysis; Gaussian filter; In situ observations; Plastic deformation
1. Introduction Optical methods can usually provide a whole field and in some cases a real-time picture of a given process, which makes them valuable means of measurement in different fields of science. For objects with optically rough surfaces, ESPI is one of the most appropriate method for high sensitivity measurements of out-of-plane and in-plane displacement components, strain, and vibrations. In addition, various image-processing techniques
*
Corresponding author. Fax: +81-48-853-3459. E-mail address: [email protected] (V. Madjarova).
could be applied, because the signal is obtained in digital form [1–4]. The deformation field in form of correlation fringes is obtained by comparing two images, one representing the reference object state and the other representing a deformed object state, either by subtracting or adding them. Conventionally in the ESPI method, the reference image is taken once, and the development of the deformation for a specified time can be received, by comparing the sequence of speckle patterns with the reference speckle pattern. In this manner, the technique could only be applied for studying deformation in a limited time, restricted by the deformation speed, CCD camera resolution, and de-correlation due to the speckle size [1]. Rather than keeping the reference image constant, in the
0030-4018/02/$ - see front matter Ó 2002 Published by Elsevier Science B.V. PII: S 0 0 3 0 - 4 0 1 8 ( 0 2 ) 0 1 9 0 9 - 0
36
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
method proposed by Suprapedi and Toyooka [5], the reference speckle pattern is renewed, and the subtraction interval is kept constant. This timesubtraction sequence method can allow the ESPI method to be successfully applied for studying dynamic events throughout their entire development. We call this time-sequence method dynamic ESPI (DESPI). We successfully applied DESPI to investigate dynamic events, such as plastic deformation of metals, in which various interesting phenomena were revealed [5–9]. These phenomena were studied qualitatively and semi-quantitatively, however, for the full understandings of the events quantitative evaluation of the phase is needed. In this paper, we aim to develop a quantitative method to obtain the phase and analyze the deformation field. There are numerous techniques for obtaining the phase that are employed presently. The essential of all these techniques is that additional modulation phase terms that can be controlled experimentally are introduced in time or space domain. Mainly three categories of methods exist: (1) fast Fourier transform (FFT) method with phase carrier proposed by Takeda [1–4,12–14], (2) temporal phase shifting (TPS) [1–4,11,15], and (3) spatial phase shifting [1–4]. In FFT method for correlation fringe analysis, a spatial carrier is introduced in either of the two speckle images. Provided that the frequency of the correlation fringes is lower than the carrier frequency it could be possible to extract the phase by FFT method. Thus, the spatial carrier limits the number of fringes, i.e., the maximum amount of deformation, which can be measured. Instead of applying a spatial carrier into the fringe pattern, Saldner et al. [14] showed that in the case of ESPI, FFT would be more successful if the carrier is applied to the speckle image itself, not to the correlation fringes. There are, though, certain requirements to be met: as the method requires the speckle size to be at least five times the size of the pixel, a camera with either high resolution or high sensitivity is needed. In addition, the use of large speckle sizes reduces the spatial resolution and as it is seen in the experimental results, places with denser fringes cannot be resolved. The FFT method can also be considered in the temporal
domain. Joenathan et al. [12] applied the temporal Fourier transform to obtain the phase. In this method, each pixel was considered independently and the Fourier transform of the temporal variation of the intensity at each pixel was performed. Some limitations and sources of error in the method are discussed in [13]. In TPS method, since the intensity is sampled temporally, the object should be stable during the acquisition of at least three frames. This makes the method difficult to apply in dynamic deformation measurements, where the object is continuously changing its state. The requirement for the object to be stationary might not be necessary if the following assumptions can be satisfied: (a) the applied load and the acquisition sampling rate are such that the time evolution of the deformation can be approximated locally into liner function; (b) the deformation rate does not produce a phase change larger than p (or jp=2j if a temporal carrier is introduced)at any position of the object between two successive frames. Then the phase can be obtained by combining five consecutive frames as proposed by Xavier Colonna de Lega and Pierre Jacquot [11]. The main drawback of the method is the assumption for almost linear deformation of the object. Moreover, if sudden change in the deformation state appears, the phase may vary faster than the value tolerable between five frames. The above-mentioned drawbacks can become significant limitations in real tension experiment, where the plastic deformation is a nonlinear phenomenon and sudden changes in the deformation state appears often. In addition, for very small deformation, the successive intensities values will be more susceptible to the noise due to CCD read out, vibrations of the object, digitization of pixels, and the deviation of the deformation from the linear trend. An alternative to TPS, in the case of dynamic events, could be the spatial phase shifting where the phase modulation is separated out in space (usually using a diffraction element) and recorded at the same time on one frame, containing several images with different phase. The main disadvantage of this technique is its susceptibility to rough environmental conditions and dependence on the CCD camera characteristics [2].
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
Reviewing the existing methods for dynamic events studies, and considering some of the restrictions of these methods, and the specifics of the events we are studying, we choose the method proposed by Yoshida et al. [10]. The reason for choosing this method is that it has no requirements for linearity in deformation or limitation in the maximum deformation (as long as the sampling frequency is satisfied), as is the case in [11]. In addition the method does not require any phase-shifting system, or spatial carrier [10]. In order to obtain the phase, both the correlation fringes from subtraction and addition images are used, and the phase is obtained from the ratio of these two patterns. We named this method subtraction–addition method (SAM). In the initially proposed method, however, there are some problems concerned with the bias intensity and unwrapping procedure. The bias intensity is received, simply by taking single beam illumination from both interfering beams and, adding these two values. The obtained bias intensity is applied for the rest of the experiment, assuming that it does not change. This can lead to large errors when carrying out experiments in which the object under study considerably alters its state. In addition, the proposed unwrapping procedure was applied only to the 1-D case for almost parallel fringes. Therefore, the successful application of this method is possible, if the method of obtaining the bias intensity allows the bias intensity value to be renewed, and if an appropriate 2-D phase unwrapping procedure is employed. In this paper, we present experimental results with an in-plane sensitive DESPI for in situ observation of dynamic events (plastic deformation of aluminum alloy). For obtaining the phase, we developed a 2-D phase unwrapping procedure for the SAM. To overcome the difficulty with the bias intensity, which was outlined above, we propose temporal averaging of sequence of speckle images and renewing the information about the bias intensity for every speckle image. This can simplify the optical set-up, because no additional arrangement, such as shutters or external phase modulator, is needed for the acquisition of the bias intensity.
37
2. Phase evaluation by SAM 2.1. Principle In general, interference speckle images obtained in ESPI can formally be represented by Iðx; y; tÞ ¼ I0 ðx; y; tÞ þ Im ðx; y; tÞ cos ðhðx; y; tÞ þ uðx; y; tÞÞ;
ð1Þ
where I0 ðx; y; tÞ and Im ðx; y; tÞ are the bias intensity and the modulation intensity, respectively. hðx; y; tÞ term is the random phase of the speckle field, which varies very rapidly in space domain. In these variables the time dependence is in general case omitted, since they usually vary significantly in the space domain but very slowly in the time domain. However, for long experiments, where the object alters significantly, this time dependence cannot be omitted. uðx; y; tÞ term is the phase introduced by the deformation of the object under study, which usually varies slowly in space and time domains. Typically in ESPI systems, to produce correlation fringes that correspond to the deformation of the object, the absolute difference between two images is taken. Addition images are not usually used, because of their very low visibility due to the bias intensity term. In the method proposed by Yoshida et al. [10], both the subtraction and addition images are used, and the phase is obtained from the ratio of these two images. The main idea of the method could be described from Eq. (2) to Eq. (4). The subtraction correlation fringe patterns for DESPI are represented by Isub ðx; yÞ ¼ Ii ðx; yÞ Iiþp ðx; yÞ Duðx; yÞ ¼ 2Im ðx; yÞ sin hðx; yÞ þ 2 Duðx; yÞ sin ; 2 i ¼ 1; 2; 3; . . . ; p ¼ 1; 2; 3; . . . ;
ð2Þ
where i indicates the frame number at the time ti , p indicates the subtraction interval, which is determined depending on the experimental conditions, and Duðx; yÞ ¼ uiþp ðx; yÞ ui ðx; yÞ is the phase change between two states of deformation. The
38
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
addition correlation fringe patterns for DESPI after removing the bias intensity are represented by Iadd ðx; yÞ ¼ Ii ðx; yÞ Iiþp ðx; yÞ I0 ðx; yÞ Duðx; yÞ ¼ 2Im ðx; yÞ cos hðx; yÞ þ 2 Duðx; yÞ cos ð3Þ : 2 The value of interest Du can be obtained by taking the arctangent function of the ratio of the subtraction image to the addition image as described by hIsub ðx; yÞi Du ¼ 2 arctan ; ð4Þ hIadd ðx; yÞi where h i denotes spatial averaging to remove the noise term, i.e., the sine and cosine terms containing the hðx; yÞ term.
Fig. 1. The time variation of the intensity at a given pixel. It demonstrates the sinusoidal variation of the intensity and also the slow temporal change in the bias and modulation intensities.
can be performed using a symmetrical time window around each speckle image. 2.3. Phase unwrapping
2.2. Elimination of bias intensity by temporal averaging over sequence of speckle images As was already mentioned in Section 1, the bias intensity also varies slowly in time. In the case when SAM is applied, using a constant value for the bias intensity, which is determined at the beginning of the experiment, will lead to considerable errors. In [10] additional arrangement in the optical set-up with mechanical shutters was employed to obtain the bias intensity. If the information about the bias intensity should be renewed, the experiment has to be interrupted to acquire the bias intensity again. This can be a considerable drawback in dynamic and high-speed event studies. Instead of using additional elements in the setup, we propose an alternative signal processing technique. In Eq. (1), it can easily be seen that, for continuously changing deformation the intensity is a cosine periodic function of time. To illustrate the temporal variation of the intensity, represented by Eq. (1), we show the intensity at a given pixel in Fig. 1, where the sinusoidal variation of the intensity can clearly be distinguished. This periodic temporal variation could allow determining the bias intensity by temporal averaging over a sequence of speckle images. The temporal averaging
The subtraction correlation fringes are obtained by subtracting two frames with a fixed subtraction interval, depending on the experimental conditions. After the bias intensity is calculated, it is subtracted from the addition image, and the addition correlation fringes are obtained. The next step is the calculation of the phase difference Du using Eq. (4). Unlike the other well-established phase calculation methods, the wrapped phase map is not obtained with the distinctive 2p phase jumps typical for the conventional phase-shifting methods. Instead the map is received as a modulus of the conventional wrapped phase map and the distinctive jumps disappear. The reason for this is the use of the modulus values of sinðDu=2Þ in Eq. (2) and cosðDu=2Þ in Eq. (3) and the sign in the wrapped phase values cannot be examined. In this paper, we propose a new 2-D phase unwrapping algorithm for the SAM. The only assumption we make is for monotonically increasing deformation. The main idea of this algorithm is shown in Fig. 2. For simplicity, 1-D linearly increasing noise-free signal is assumed in this graph in a limited interval of phase change from 0 to 4p. The two rectified sinusoidal functions represent the sine and cosine functions cor-
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
39
nðx; yÞ nðx;yÞ1 Duðx; yÞ ¼ 2p Du0 ðx; yÞ þ ð1Þ 2 ðn ¼ 1; 2; 3; . . .Þ;
ð5Þ
where ½
is the Gaussian symbol that stands for an integer number, n is an integer, Du is the unwrapped phase value, and Du0 is the wrapped phase value.
3. Experiments and discussions 3.1. DESPI system for in situ observation
Fig. 2. Illustration of the 1-D unwrapping procedure. The sine and cosine function are given in relative units.
responding to the intensity values of the subtraction Isub and addition Iadd , respectively. The wrapped phase Du calculated from the arctangent of the intensity ratio of the subtraction and addition fringes is shown as the triangle-shaped function, while it takes a saw tooth shape function in the conventional interferometry. Taking into account the assumption we made for the deformation, a solid straight line could be drawn to represents the unwrapped phase value. Although we make an assumption for the object deformation to be monotonically increasing, since the method has a sign ambiguity, there are many practical applications where the assumption is satisfied, such as the tensile experiments described in this study. The first step in the algorithm is to divide the interval into four areas, and assign an index number in each area. The assigned numbers in the areas correspond to the fringe order; therefore the area is indexed as Areas 1–4. The borders of two adjacent areas correspond to the position of the maximum or minimum of the wrapped phase values. These points are determined by scanning the positions of zero gradients from the obtained wrapped phase map expressed by Eq. (4). After correct index numbers have been assigned on the whole area, the phase values at the position ðx; yÞ are unwrapped following the rule:
The method described in Section 2 can be applied to in-plane as well as out-of-plane deformation measurement, because it was considered for the general ESPI. In this paper, we carried out a tensile experiment of metal object, where an inplane sensitive DESPI system was used. The experimental system is illustrated in Fig. 3. The light source was a second harmonic generation of a YAG laser with output power of 200 mW, and wavelength of 532 nm. The light emitted from the laser was expanded by an objective lens L to illuminate the specimen. Here, the upper half part of the expanded beam directly illuminates the object with an angle a, while the lower half part is reflected by the mirror M3 to illuminate the object
Fig. 3. The scheme of the optical set-up for in-plane deformation measurement. M1, M2, and M3 – mirrors, L – objective lens, CCD – camera, DL – data logger, PC1, PC2 – computers, and S – specimen.
40
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
with an angle of a. The propagation vectors of the two branches lie in the plane parallel to the z–y plane and the sensitivity vector lies in the x–y plane. Strictly speaking, the configuration of two divergent beams implies also the presence of sensitivity vector component that lies outside the x–y plane. Taking into account that the magnitude of the sensitivity vector that lies out of x–y plane is very small and that the displacement in the tensile experiment is mainly in-plane, we consider only the y-component of the in-plane displacement. The amount of deformation per unit fringe spacing is given by k=ð2 sin aÞ and was estimated to be 380 nm, with the incident angle a ¼ 45°. The CCD camera was positioned perpendicularly to the surface of the specimen and the interference speckle patterns are taken continuously with a constant acquisition rate. In the experiment, the camera acquisition rate was 2 frames per second (fps). The captured images are subtracted to get correlation fringes with a time interval of 4 s and are displayed on the monitor in real-time. Fifty-six thousand frames were acquired in total during the entire experiment. We performed a tensile experiment of a dumbbell-shaped plate of aluminum alloy A2017 with the effective length of 80 mm, width of 30 mm, and thickness of 5 mm. Before the tensile experiment, the specimens were annealed by exposing in atmosphere of 450 °C for 1 h and slowly cooled to the room temperature. The specimen was stretched on a tensile machine with a constant tensile speed, in a manner that the lower hole of the specimen was caught on a cylinder rod, which was fixed to a cramp grip, and the upper hole was caught on a cylinder, which was fixed on a moving crosshead. Several experiments were performed with different tensile speed of 0.5, 1, and 1.6 lm/s. One of the results from the experiment with the tensile speed of 0.5 lm/s is shown in Fig. 4. Fig. 4 presents the stress curve with the corresponding correlation fringes patterns obtained at different states of the deformation process, which are shown below the stress curve. The stress increases linearly in the elastic deformation state, followed by zigzag variations in the plastic state. The zigzag variation or serration is known as the phenomenon associated with strain aging called Portvin–Le Chartelier
Fig. 4. The stress curve in function of time with the corresponding fringe patterns. They demonstrate the complicated manner in which the deformation develops during the experiment.
effect [16]. The pictures presented demonstrate the development of the plastic deformation in rather complicated manner. Pictures 3 and 4 represent a stable band-like pattern at an angle of approximately 45° to the tensile axis. That pattern suddenly appeared many times during the experiment (when the stress curve changes sharply), and the position where it appeared moved over the specimen upwards or downwards. At the later stage of the plastic deformation state, there were moments at which the angle of the band pattern switches, i.e., the angle changes from approximately +45° to approximately )45° to the tensile axis (pictures 5, 6, and 7). 3.2. Phase analysis by SAM After the frame acquisition, the quantitative analysis based on SAM was performed to obtain the object phase. The procedure included: (1) calculating the bias intensity for every frame and subtracting it from the corresponding frame; (2) calculating the subtraction and addition images and applying filters for removing the speckle noise; (3) calculating the wrapped phase values using addition and subtraction images; and (4) applying the unwrapping procedure. There are two main sources of errors in the SAM with temporal evaluation of bias intensity
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
that reduce the accuracy of the measurement. The first source of errors is the inaccurate determination of bias intensity. Mainly, in real experiments, the accurate determination of the size of the temporal window that should include an integer number of cosine function periods in Eq. (1) is difficult to accomplish automatically, for the whole specimen and for the whole experiment, since it varies in time and space. Considering the fact that the averaging is performed over several intervals of a periodic function, the influence of the fraction of the period over the average value tends to zero with the increasing number of periods inside the window. In our experiment, the temporal averaging window was determined in a manner that it includes at least several periods of the cosine function. Another limitation in the size of the temporal window arises from the fact that the bias intensity is not a constant value. If the temporal window is too wide the bias intensity cannot be regarded as constant inside the window, and the average value received will not be its correct representation. Taking into account the tensile speed and the variation in the bias intensity, it was found that the optimum temporal window is 25 s, which includes on the average more than 6 periods of the cosine function in Eq. (1). The second source of errors comes from the inaccurate determination of the Areas 1–4, etc. Mainly the errors arise when the noise in the images is not filtered out successfully, which will bring inaccurate maximum and minimum positions in the wrapped phase. In order to get a
41
correct result on the phase unwrapping process, application of a filter to minimize the fluctuation due to the speckle noise on the obtained signal becomes inevitable. The noise on the obtained signal is associated with the presence of the high-frequency fluctuation caused by the speckle random phase h, which is uniformly distributed over all values between 0 and 2p. Taking into consideration that h and Du are statistically independent, and the spatial frequency of h is significantly higher than that of Du, a low-pass filter in the spatial domain can be used to minimize noise due to the presence of h. Performing an optimum low-pass filtering, the following approximation of hsinðh þ Du=2Þ sinðDu=2Þi
2 sinðDu=2Þ p
and 2 cosðDu=2Þ p can be achieved with significant accuracy, and consecutively phase value Du can be obtained. The brackets h i denote the spatial filtering operation. In this research, a Gaussian low-pass filter is applied to filter out the high-frequency term due to the random phase h. An important note considered here is that the filter must be adjusted depending on the spatial frequency of the obtained correlation fringes. Fig. 5 gives the results of addition, subtraction, and wrapped phase of a typical band-like pattern that makes angle of approximately 45° to the hcosðh þ Du=2Þ cosðDu=2Þi
Fig. 5. Addition, subtraction, and wrapped phase images with the 3-D plot of the deformation field representing the locally concentrated deformation field. The subtraction interval was 4 s, and temporal averaging window was 25 s.
42
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
Fig. 6. Addition, subtraction, and wrapped phase images with the 3-D plot of the deformation field representing the ‘switching’ of the position of the locally concentrated deformation field from approximately 45° to )45°.
tensile axis, which appeared several times during the experiment. The visibility of the addition correlation fringes is comparable to that of the subtraction correlation images. The 3-D plot of the deformation field after the unwrapping of the phases is also given in this figure. It can clearly be seen that the deformation is concentrated in a very narrow area, which we call locally concentrated strain band. The temporal behavior of this locally concentrated strain band was discussed in more detail in [5,17]. The ‘switching’ of the direction of the band from approximately +45° to )45° to the tensile axis can be seen in Fig. 6, where again addition, subtraction, and wrapped phase images are shown with the corresponding 3-D plot. The amount of deformation after the ‘switching’ is less, but it changes in a pulsating manner as it was discovered and discussed in [5,17]. The interval between the two states is approximately 130 s. The accuracy of the system was evaluated by estimating the standard deviation of the deformation for the area where almost no deformation was detected. It was on the average 50 nm, which corresponds to k=10.
4. Conclusions In this paper, we proposed a dynamic ESPI system for in situ observations and a 2-D subtraction–addition method for evaluating the phase in DESPI correlation images. The bias intensity in the addition speckle correlation images was de-
termined by temporal averaging over a sequence of speckle patterns. The phase map is not obtained with the distinctive 2p phase jumps, but as a modulus of the conventional wrapped phase map and the distinctive jumps disappear. Thus the main limitation of the method is the sign ambiguity and an a priory information for the object under study should be used for the correct unwrapping. However, there are many practical situations, where the proposed method can successfully be applied, as demonstrated in this study. In tensile experiments, the deformation direction is known a priori and is usually monotonically increasing. To perform the phase unwrapping, the folding lines of minimum and maximum of calculated wrapped phase distribution were utilized. Prior to the unwrapping procedure Gaussian filters were applied to reduce the speckle noise. In the present experiment the acquisition rate of the camera is 2 fps, but the system and the method can be applied with highspeed camera, because there is no limitation on the acquisition rate of the camera that is required by the method itself. It was demonstrated that DESPI with the proposed method for obtaining the phase could successfully be applied to dynamic events with an accuracy of approximately k=10. The main objective in developing the DESPI with SAM for obtaining the phase is to expand the possibilities of ESPI in studying dynamic events in rough environmental conditions. The method for receiving the bias intensity allows the application of the system in studying high-speed phenomena. The proposed method is expected to contribute to
V. Madjarova et al. / Optics Communications 212 (2002) 35–43
quantitative analysis of deformation states in different materials, as it is not confined only to tensile experiments.
Acknowledgements This research was partly supported by a Grantin-Aid (13555193) for Science Research of the Ministry of Education, Science, Sport and Culture, Japan.
References [1] P.K. Rastogi, Digital Speckle Pattern Interferometry and Related Techniques, Wiley, Chichester, 2001. [2] D.W. Robinson, C.R. Raed (Eds.), Interferogram Analysis: Digital Fringe Pattern Measurement Techniques, Institute of Physics Publishing, Bristol, 1993.
43
[3] P.K. Rastogi (Ed.), Optical Measurement Techniques and Applications, Artech House, Boston, 1997. [4] P.K. Rastogi (Ed.), Photomechanics, Springer, Berlin, 1999. [5] Suprapedi, S. Toyooka, Opt. Rev. 4 (2) (1997) 284. [6] S. Toyooka, R. Widiastuti, Q. Zhang, H. Kato, Jpn. J. Appl. Phys. 40 (2001) 873. [7] S. Toyooka, X.L. Gong, Jpn. J. Appl. Phys. 34 (1995) L1666. [8] Suprapedi, S. Toyooka, Phys. Mesomech. 1 (1998) 51. [9] X.L. Gong, S. Toyooka, Exp. Mech. 39 (1999) 25. [10] S. Yoshida, Suprapedi, R. Widiastuti, E.T. Triastuti, A. Kusnowo, Opt. Lett. 20 (1995) 755. [11] X.C. de Lega, P. Jacquot, Appl. Opt. 35 (1994) 5115. [12] C. Joenathan, B. Franze, P. Haible, H.J. Tiziani, J. Mod. Opt. 44 (1998) 1975. [13] C. Joenathan, P. Haible, H. Tiziani, Appl. Opt. 38 (1999) 1169. [14] H.O. Saldner, N.E. Molin, K.A. Stetson, Appl. Opt. 35 (1996) 332. [15] J.M. Huntly, G.H. Kaufmann, D. Kerr, Appl. Opt. 38 (1999) 6556. [16] G.E. Dieter, Mechanical Metallurgy, McGraw-Hill, London, 1988. [17] S. Yoshida, S. Toyooka, J. Phys. 13 (2001) 6741.