- Email: [email protected]
One-step phase shift 3-D surface profilometry with grating projection
Optic.~ and La~er~ in Engineering 21(1994) 61—75
© 1994 Elsevier Science Limited Printed in Northern Ireland. All rights reserved 0143-8166/94/$7OO
ELSE VIE R
One-Step Phase Shift 3-D Surface Profilometry with Grating Projection Ruowei Gu, Toru Yoshizawa* & Yukitoshi Otani Department of Mechanical Systems Engineering, Tokyo University of Agriculture and Technology, Nakacho 2-24-16, Koganei, Tokyo 184, Japan (Received 30 September 1993; revised version received 14 December 1993; accepted 20 January 1994)
ABSTRACT A 3-D profiling system has been developed using a grating projection apj)r()ach and a one-step phase shift algorithm. In the system, a grating pailein projected on an object surface is acquired by a CCD camera, and the grating’s phase deformation caused by the surface shape is extracted by spatial phase shift processing, which uses only one frame of digital image. One phase value can be calculated from a successive 3-pixel range if the fringe period is set to either 4-pixel width or 8-pixel width. The fixed phase distribution of the system is excluded by a standard plane calibration. With 8-bit input data, the system’s RMS phase accuracy is developed up to 2jr/60 in the experimental examinations.
INTRODUCTION In the field of optical non-contact 3-D measurements, more and more researchers have turned to utilizing phase measuring techniques. The methods presented consist of techniques such as Fourier transform profilometry (FTP),’ phase-shift moire topography (PMT),2 ~sinusoidal spatial phase detection (SPD),5 multi-step phase-shift grating projection (MPGP),ON one-step phase-shift grating projection (OPGP),~etc. The moire method requires precise optical arrangements in order to form flat-contour fringes. It is quite difficult to apply the PMT methods to the crossed-optical-axes systems in non-telecentric cases because the contour fringes are not planar in that case. But by FTP, SPD, MPGP and OPGP methods, 3-D data can be directly extracted from a To whom correspondence should he addressed. 61
K Go, I. ~o.s1,,:a,to, 1.
(hoot
projected grating and the deformed grating phases can be separated into the objective distributions and the fixed optical distributions in the crossed-optical-axes occasions. Iherefore 3-D coordinates can be easily computed after the standard plane calibration of the fixed phase deformations. FIP and SPD are two excellent methods since only one image frame is required, but the phase calculations are complex and need a great deal of gray level data to extract oiie phase value. It IS necessary to use special digital processors to get high processing speed for the two methods. With MPGP approaches, if the mechanical niicro—shitting structure for obtaining several phase-shifted fringe patterns can ensure high locating precision. pure pixel-by-pixel resolution and high accuracy can be obtained. But the system accuracy is usually limited by the mechanical micro-shifter and the measuring speed is difficult to increase. By use of one-step phase shift projection methods OPGP), only one frame of grating image is required. optical crossed— axes deformation can he removed, the precise mechanical shifter is not demanded and all of the phase processes can he in binary-integer-form by using an arc—tangent look—up table. It can also achieve as high a computing accuracy as MPG P but its spatial resolution in the image plai~e of CCD is a little lower than MPGP because it needs at least three adjacent pixels’ gray data to calculate one phase value. It would seem that OPGP is one of the most advantageous approaches for realizing dynamic 3—I) measurements in pr;ictic;tl applications.
ONE-STEP PHASE SHIFT ALGORITHMS In a previous work the authors, referred to the idea ol spatial phase shift methods in inierlerometry ‘~ and established a one—step phase shil graling projection system. Recenik the~’ developed the system for crossed—optical—axes occasions and are now improving it with many new features. The known formula of the one—step phase shift method needs the pitch of input grating patterns to occupy approximately 4 pixels, so I hat the adjacent pixels on the lines vertical to fringe direction cati he regarded as having ~/2 phase shifting values between each other. lheielore the grating pattern can be written as follows if the fringe direction is along vertical rows of ((1) array 1(1/) \\ ide
i
/
tic
Thi
I)
liori,ontal
((/. I) dos
j
(/)(i,
j.:)
(I)
and \ei’tical pixel-counters. respectively.
3-D
surface
profilometrv
63
UI
L
Grating Projector
Reference Plane
w~DiA
CCD Camera Fig.
1.
Diagram
of
/
crossed-optical-axes grating projection measurements.
system
for
3-D
B(i,j) and C(i,j) are background and fringe contrast, respectively, and /(i, j, z) is the deformed grating phase which relates to pixel positions
and object height distributions. Using eqn (1), it can be derived that ~(i,j,
z)
=
[l(i+2,j)-l(i+1,j)] arctan ~ I(i,j) —I(i + 1,j)
~
(2)
where l(i, j), I(i + 1, j) and l(i + 2, j) are 3 adjacent pixels’ gray data, and eqn (2) should satisfy the condition that the background B(i,j) and contrast C(i, j) are not changing within a 3-pixel-range. The crossed-optical-axes grating projection system used is shown in Fig. 1, where the center of a projector’s exit pupil P and the center of a camera lens entrance pupil I determine the orientation of the reference plane’s X-axis, and the origin point 0 of the objective dimensions is at the intersection of the projector’s axis and the camera lens’ axis. The camera lens’ axis should be set vertical to the X-axis so that the CCD’s imaging dimensions can overlap with the reference plane dimensions. The z value is set as positive when the object surface is on the left side of the reference plane. The parameter L is the distance from both P and 1 to the reference plane, and the parameter w is the distance between P and I. Consider one light ray from point P which is projected onto a surface point A with height value —z. The ray PA goes through a point C in the reference plane and the imaging ray Al goes through a point D in the reference plane. The point A’ is the virtual intersection of the ray PA and the grating’s natural image-plane. All the points on the ray PA’ have the same fringe order of the projected
64
I?. Gu, 1. Yoshizawa, Y. Otani
grating. In the telecentric projection case, the ray expressing the same fringe order as the point A’ will he parallel to the axis PG and it will cross the reference plane at the point B. Physically, if two points on the reference plane in Fig. I express the same fringe order of the projected grating but are in different spatial positions, the distance between them is just the factor of the deformed grating’s phase. Therefore, if we regard the fringe periods as 4-pixel width or 8-pixel width, the grating pattern imaged on the reference plane can he written as follows 1(x, v)
=
B (x,
y)
+
((v v)
x cos [2irfx + 2irf0CD(z) + 21Tt)BCC
) + 2r~\Lv~
(3)
where B(v, v) and ((x, v) are background and fringe contrast. respec1 1 or (a is the pixel size equivalent to reference 4a So plane dimensions), CD(z ) is the phase deformation due to object’s tiyely,
f, equals
-
surface height variation and it is easy to find that: CD(:) =
(4)
~—~---
BC(x) is the phase deformation on reference plane due to crossedoptical-axes effects and it is solely dependent Ofl x values. And ~XIis the frequency difference between f~and the averaging practical frequency of acquired fringe patterns. Rewritting eqn (3) 2irfCD(z)
+
2.~f0BC(x)+ 2ir~f v
=
q~(x,V. ~)
(5)
and transforming x, v values to pixel-counting discrete forms, i.e,y ía and v = Jo (I, j are integers), we can change eqn (3) to the following form I(i,j)
When I,,
=
B(i,j) ~ C(i,j)cos [2irfai
+ ~(i,j,
~)J
(6)
(4 pixels per period), eqn (6) is same as eqii (I ). and the
=
4(1
phase /(i. j, :) can he extracted by eqn (2). When f~
,
i.e. the case
of S pixels per period, we can simply transform eqn (2) to the following (/)(1,J,z)=arctan
-
I(i+4,j)~I(i+2,/) I(i,~)~I(t+2,~) --~
—
~ i 4
(7)
This means that the 3 gray data for phase extraction are obtained from
3-D surface profilomerry
65
3 interval positions over a continuous 6-pixel range. If the object’s surface is smoothly changing, eqn (7) will work well, but the spatial resolution will decrease to half of that of the 4 pixel/period case. In order to keep the resolution as high as possible, we have derived a new formula to compute one phase value from a continuous 3-pixel range as follows (/)(i,
I, z) =
[
arctan I (V~ 1). —
L
j)
l(i j)—I(i+2 .
.
.
‘
.
21(z+1,j)—I(i,j)—I(i+2,j)i
1 .
I
Jr — —
4
i
(8)
Equation (8) keeps the same spatial resolution as that of the 4 pixel/period case. If the input picture has a high frequency background pattern or contrast variations, eqn (8) will work better than eqn (7) because eqn (8) uses gray data only in a 3-pixel range. In comparison with multi-step phase shift methods, the way of reducing measurement errors is quite different in one-step phase shift methods. For multi-step methods, the errors caused by mechanical shifting can be averaged by increasing the number of phase-shifted gray data. But in one-step occasions, if we use more gray data in one picture, a longer, continuous pixel-range is used and then the resolution in the imaging plane will be lower and more background noise effects will appear. Therefore the pixel-range for extracting one phase value should be shortest.
PHASE PROCESSING PROCEDURES Phase extraction We realized all the three algorithms shown in eqn (2), eqn (7) and eqn (8) for calculating phase maps from one frame of grating pattern. Depending on the input fringe’s frequency, either 4-pixel-period method or the 8-pixel-period method may be selected. The arc—tangent calculations of the 3 formulas are all executed by using a 16-bit look-up table. Therefore the computing time for extracting a 256 X 256 x 16-bit phase map which is wrapped —Jr to it is controlled to less than 1 second even using a PC/AT computer. It should be illustrated that a shifting constant-phase is attached to each phase value when using eqn (2) or eqn (7) or eqn (8). It is removed by anti-shifting computation while the arc—tangent calculation is performed in the system. This part of the computing program was written in PC Assembler language.
I? Go. i
66
Yoshi~ati’a, Y. ()tani
Phase unwrapping A serious problem for all the phase-shifting techniques is the phase unwrapping after the arc—tangent calculations. If each phase point is unwrapped along some fixed direction, it is always possible to propagate errors from had-quality points to good-quality points. In Ref. II Kwon reported briefly an excellent unwrapping method proposed h~ Shough. It is called ordered phase unwrapping (OPIJ). The idea of this method is that the ‘best’ point among the boundary points of an already-unwrapped area is selected for unwrapping. In our system, we referred to the idea and devoped an OPU procedure using PC Assembler language. For such processing, it needs a quality factor to determine what is the ‘best’ quality point. We made this factor Q, as: =
~
~
I~
--
(~)
Q,
is the quality factor of phase point P, and P,~, is one of 4 adjacent phase points around P, While getting the absolute phase differences, the 2ir jumps should he excluded if they exist. With the quality factor, an optimal order for unwrapping all phase points can he determined automatically and the propagation of errors will he kept as low as possible. The brief flow chart of the ordered phase unwrapping algorithm in our system is indicated in Fig. 2. where
Bad-region recognition In order to get results that only represent good phase points, a had-region recognition procedure has been developed. Ihe had-region recognition can be performed on the it to Jt wrapped phase map because the no-fringe areas and the poor-quality-fringe areas are most --
obvious here during the whole processing. The main feature of the had areas in wrapped phase maps is the phase values jumping randomly in high spatial frequency, so we selected a rectangular area around an inspected phase point and counted the surrounding points whose phase differences to the inspection point satisfy the following condition --
<2Jr
—
~,
(10)
where P is the phase point in the surrounding area, P, is the inspected phase point. and is an allowance of legal amplitude of phase variations. If the number counted is greater than a given threshold, the inspected point will he marked as a had point and will not he processed in phase unwrapping. Equation (10) also considers the case of 2ir jumps during the recognition. ~,
3-D surface profilometry
67
Select a start pixel on central column
1
Find neighbor pixels which are not unwrapped
Get the quality factors of selected pixels and add them in an ordered table if they are not in the table
Table length -l and pick out the best quality’s pixel from the table
If
Yes
_______
table length=-l
End
Unwrap the pixel with one of processed neighbors which has the best quality Fig. 2.
Flow chart of ordered phase unwrapping algorithm.
Correction of fixed optical phase deformations The phase distribution calculated either from eqn (2), eqn (7) or eqn (8) consists of the objective phase, the cross-optical-axes phase and the frequency fraction phase, as shown in eqn (5). Because the last two terms of eqn (5) are not influenced by object shapes, they will remain constant if the optical system is fixed. Therefore we can measure a standard plane as the zero-distribution of the system at first, and store its phase map as the system’s base data 4(x, y). Since the object height distribution is assumed to be zero here, we can get the base distribution from eqn (4) and eqn (5) as 4~(x,y)2Jrf,,BC(x)+2itt~f,,x
(11)
Parameters of eqn (11) are the same as illustrated in eqn (5). We do not need to know the theoretical form of BC(x) because its assignment has already been obtained pixel by pixel in the standard-plane phase map.
K Go, I
68
Yoshizawa.
Y Otani
After phase unwrapping for a measured object, the pure objective phase map can he extracted by subtracting the base phase map from it. Actually, the errors in optical arrangement are also included in the measured base phase map, so the objective phase map is the pure distribution referring to the standard plane. So far, the phase processing system is built and it can ensure the accurate phase measurements by using experimental calibrations of the standard plane. Improvements with grating tilt Originally, grating patterns were set vertical to the CCD’s scanning lines. But after processing, usually some vertical-line-type’s noise pattern with the order of 2it/lO appears on the phase map, and the frequency of it depends on the input fringe’s frequency. One reason for it is attributed to the arc—tangent calculation for phase extraction because it is hard to he accurate near ir/2 or —jr/2. This means that the phases near the two values are extremely sensitive to the A/D sampling errors of sinusoidal fringe shapes. Figure 3(a) shows a phase map of a flat surface using vertical fringe patterns. Since the fringe patterns are almost the same along the vertical rows of CCD array, the sampling errors will also he like that. If the errors on a row are serious, all the row’s phase values will he deformed and even if the phase map is modified by median-value-filter, the noise pattern is not obviously reduced, as shown in Fig. 3(h). If the RMS deviation of a planar phase map is about 2ir/3() originally, it will only he reduced to about 2.ir/50 after 3 X 3 median-value-filter. The noise pattern will confuse the fine distributions of the measured phases. Noting the fact that the fixed optical phases can he excluded by standard plane calibrations, we changed the projected grating tilt slightly on the Y-axis and calculated the phases along CCD’s scanning lines. The phase expressed in eqn (5) will change as q~(x,v, z) where
(Y
=
2iri~C’D(:)+ 2irf,BC(v) + 2irMv + 21rf%’ tan
c~
(L2)
is the angle between the Y-axis and the grating orientation,
and the other parameters are the same as mentioned in eqn (5). For the standard plane measurement, CD(:) = 0, the phase niap becomes (/)/.V,
v, z )
=
2itfBC’(x ) + 2ir~f1v+ 2itf,,v tan
(V
(13)
In eqn (13), rio parameter relates to object height: therefore the effects of grating tilt can also he removed by the standard plane calibration. The distinct feature between vertical grating projection and
69
3-D surface profilo~ne1ry
(a)
(h)
(c)
(d)
Fig. 3. The difference in the computed phase maps of a flat surface measurement comparing the vertical grating projection case with the tilted grating projection case. (a) Raw phase map in vertical grating case: (b) 3 >< 3 filtered result of (a); (c) raw phase map in tilted grating case; (d) 3 x 3 filtered result of (c).
tilted grating projection lies in the results of phase processing. By means of tilted grating projection, the fringe pattern sampled by CCD array is made with more delicate profiles, both in the X-axis and the Y-axis. The calculated phase map does not carry uniform error along vertical rows, and the error pattern is limited to periodic variations both in the X-axis and the Y-axis as shown in Fig. 3(c). Such a noise pattern can be removed by median-value filtering more successfully than in the vertical grating case. Experimental results show that if a planar phase map has the RMS deviation as 2ir/30, but its noise pattern is as in Fig. 3(c), its RMS deviation will be reduced to less than 2ir/100 by only 1 application of a 3 x 3 filter. Figure 3(d) shows the filtering result from
Fig. 3(c) and it is obvious that the noise pattern reduction in tilted grating cases is much quicker than in vertical grating cases.
71)
K
(.i,,
I. Yoshi~atta, 1 Otan,
Here it should he noted that the 3 X 3 median-value-filter processing on phase maps is reasonable because originally each phase value is extracted from 3 adjacent pixels’ gray data, i.e. it is obtained from the objective information in each small area corresponding to 3-pixels’ size. Therefore the system’s X—Y resolution is suitable to he set to the scale of 3 >< 3 pixels.
EXPERIMENTS In our experimental system. the measuring distance shown in Fig. I was set to I. = 110 mm, and the width between the points P and I was set to ii’ = 60 mm. The average grating pitch projected on the standard plane was measured at about 0~4mm and the measured size corresponding to 256 x 256 pixels was 24 mm X 24 mm. Therefore the equivalent pixelsize on the objective reference plane was about 0~094mm. The analog video signal of the deformed grating patterns captured by a C(’D camera was converted to a digital picture with 512 X 512 pixels and 8 hit gray-levels, and transferred to a PC486 microcomputer through an image processing board. All the processing software was written in C and Assembler languages and it was designed for the MS-Windows operating-environment, combining 256 colors or 256 gray-levels bitmap display. As an example, we measured a one hundred Yen Japanese coin which has a pattern varying in a 0’2 mm range. Figure 4 shows the
Fig. 4.
I he di~ii~il gritting pallet ii ol a coin ~viih 1)2 mm heighi \arlailons. It is acquired I morn ihe optical s~siem sho~~ n in Fig. I.
Fig. 5.
Ihe ora’ —coded phase disit thu ilon coiiipiiied tom I ig. 4. inien’.ti~ daLi h~ OI1L’—siej) phase ~.hmii algormihrns. I lie hrighiesi e~ci presenis phase .‘r .imid 11w darkest level pm esenis phase .‘r
3-D surface projilometrv
71
input fringe pattern with the size of 256 X 256 pixels selected from original 512 X 512 pixels digital image. The grating orientation was set with about 10 degrees tilt to the vertical row of the CCD array. Figure 5 shows the gray-coded display of the phase assignment that is wrapped in —Jr to it range corresponding to the range from the darkest gray-level to the brightest gray-level. It was obtained from Fig. 4’s gray data by the one-step phase shift processing, pixel by pixel. Figure 6 shows the bad-region distribution of Fig. S’s phase data and it was recognized by eqn (10) in the case of selecting 5 x 5 pixels as the surrounding area and setting the threshold to 12, i.e. half of the investigated points. The recognition means that the inspected point will be assigned as a bad point if its surrounding illegal-phase-variations are more than the threshold. Figure 7 shows three sequential intermediate situations during the ordered phase unwrapping processing for Fig. S’s phase data. The white areas represent the bad regions which will not be processed, the black areas are the legal areas which will be processed, and the gray areas represent the already unwrapped areas. The gray area grew automatically at the best quality point of the area’s boundary (at any one time), so it ensured the least error propagation in phase unwrapping and all the bad areas were passed regardless of any shape. Figure 8 is the phase map obtained after phase unwrapping, and it contains the object phase, the crossed-optical-axes phase, the frequency fraction phase and the grating tilt phase as expressed in eqn (13). Figure 9(a) shows the result after the standard plane phase calibration of Fig. 8, and the displayed pattern represents the gray-coded phase distribution of
Fig. 6.
Bad-regions of Fig. S’s phase map. They were recognized by expression (10) with 5 X 5 investigating rectangle.
72
1-? (i,,. I. Yo~1,,:oima.
I
~
(l~
( )iiini
}
~ ~
I
~
C
t
‘I..,
~
.‘~:..
~.
.
/
(h)
(a)
Fig. 7.
I lilL IrrlLriiiLdr.lt. ~rl1I,iiI,os at ItL oideicd ph Is~ uo~~rrppiiielI~sr’.’. Ii I ~ pIi,r~~ ~t ii I I IR v~liii. rc,~iori’l~pi~‘LII) the had—mcr1iorl’ \~hich \\crc riot pos.~.‘‘cd lw hi ILk id’. ri. h~ ic,~i1 ilL I’. ~~hrch ~siII tic un~smappLd.ihL ~ ~rtIcrrr i’. IlL .iiic ids tiii’~i rppLd Iic.i ~~hicir gmr’\\s it the hcst—qualit~point oh Ii’ isruod IS each huh rid lu. iiio~sirie‘.hipL. i~shown us (a) —(e) ru ‘~cquciIcc.
the measured coin. Figure 9(b) is the wire—frame plot of Fig. 9( a pl~se distribution. The phase values wrapped it to it are kept in the computer’s memory as a I fi hit binary—integer—form, and after unwrapping es cr\ —
phase point keeps a 24 hit-length in computer memories or disk tiles. We made an accuracy examination for the whole system b~locating a plane surface in the objective space afldl nieasuring its phase distribu— lions at several distances. Through the statistical analyses of all the 251) 25( ) phase points in the center ot the phase maps, s~ e got the ~‘
3-D surface profllometry
Fig. 8.
73
The phase map after phase unwrapping. It contains the object phase and the
fixed optical phase.
~i) (a)
• (b)
Fig. 9.
The object phase map extracted from the standard plane calibration of Fig. 8’s phase data: (a) is the gray-coded display of the phase map and (h) is its 3-D wire-frame
plotting chart. The changing scale of the surface pattern was measured in 02 mm.
distribution of RMS errors along the Z-axis as shown in Fig. 10. The RMS values near the standard plane position were 2,r/160 (in ±FOOmmrange), and up to ±S~00mm from the standard position the RMS deviations were 2ir/60. The reduction of the measurement precision is mainly due to the optical distortion of the projector and the camera lens in defocused imaging occasions. The optical distortion makes the plane surface data increasingly curved when the defocused value increases. Considering that a 2ir value of phase is equivalent to one fringe period p and the P—V deviation is regarded as 3 times the RMS deviation, the surface measuring accuracy can he determined as RMS X 3 X P/2ir. So the system accuracy can he within 20 p~mof the
74
Ii’. (iii, I. Yoshi:aoa. }. Utani 0.040
i~’0.032 2.
•~,r..r,.,,u
-
-.
0.024
-6.00
-4.00
-2.00
0.00
2.00
4.00
600
Plane position along Z axis (Unit: mm) Fig. 10. lhc experimental RMS—error drsirrhutiiin Ii) the whole sssteiii. I) SS:I~ ohiarned by plane surface measurements with scs cral displacements in the stamldIurLt plane ioe~uiiofl.
4
Fig. 11. Ihe 3—I) wire—trame plot of the phase measuremeni at a sample surtace wiiii two known steps. In the result, the high step represented 100 /Lm varration and the ass siep represented So ~ variation.
systeni configurations such that p equals 0’4 mm: the measuring depth is 10 mm and the measuring size is 24 mm X 24 mm. ‘l’o demonstrate the accuracy, we measured a sample surface with known 5)) ~rm and tOO ,um steps on it. Figure 11 shows the wire-frame plotting result of its phase map. The high step represents the 100Mm and the low step
represents the 50
jwi.
3-D surface profllornetry
75
CONCLUSION We have developed the one-step phase shift techniques into 3-D surface measurements, and the one-step phase shifting computation can be executed in 3-pixel-range in the cases of either 4 pixels per fringeperiod or 8 pixels per fringe-period. The areas without fringes or with poor quality phase data can be recognized as bad regions by investigating the irregularity of phase distributions. An optimal method called ordered phase unwrapping was adopted to remove the 2ir jumps in arc—tangent calculated phase maps. The unwrapping procedure always seeks the best-quality point among the boundary phase data of processed areas. It can go around any type of bad region and ensure the least error-propagation in phase unwrapping. The fixed optical phase distributions including crossed-optical-axes phase deformations and system dislocating deformations, can be excluded by a standard plane calibration. This ensures the pure objective phase extraction by experimental techniques, but not by computing calibrations. This is significant in practical applications. In order to restrain the regular noise pattern corresponding to input fringe patterns, the fringe orientation of the projecting grating is set with a small tilt to the CCD array’s vertical rows. In this way, the noise pattern can be controlled both in horizontal and vertical direstions, and can be further reduced by 3 x 3 median-value filtering more rapidly than in the case of vertical gratings. The influence of the grating tilt can also be excluded by standard plane calibrations. Experimental results showed that the RMS accuracy of the phase processing improved up to 2ir/60 in the objective dimension of 24mmX2SmmX 10mm. REFERENCES 1.
2. 3. 4. 5. 6. 7. 8.
9. 10. 11.
Takeda, M. & Mutoh, K., App!. Opt. 22 (1983) 3977. Reid, G. T., Rixon, R. C. & Messer, H. I., Opt. Laser Tech., 16 (1984) 315. Gu, R. & Zhou, R., In Proc. SPIE, Vol. 1230, 1990, 638. Yoshizawa, T. & Tomisawa, T., Opt. Eng. 32 (1993) 1668. Toyooka, S. & Iwaasa, Y., App!. Opt. 25 (1986) 1630. Srinivasan, V., Liu, H. C. & Halioua, M., App!. Opt. 24 (1985) 185. Komatsubara, R. & Yoshizawa, T., J. Japan Soc. Precision Eng. 55 (1989) 85 (in Japanese). Su. X., Bally, G. V. & Vukicevic, D., Opt. Comm. 98 (1993) 141. Gu, R. & Yoshizawa, T., Proc. SPIE, Vol. 1720, 1992, 470. Shough, D. M., Kwon, 0. Y. & Leary, D. F., In Proc. SPIE, Vol. 1221, 1990, 394. Kwon, 0. Y., In Proc. SPIE, Vol. 816, 1987. 196.
© 1994 Elsevier Science Limited Printed in Northern Ireland. All rights reserved 0143-8166/94/$7OO
ELSE VIE R
One-Step Phase Shift 3-D Surface Profilometry with Grating Projection Ruowei Gu, Toru Yoshizawa* & Yukitoshi Otani Department of Mechanical Systems Engineering, Tokyo University of Agriculture and Technology, Nakacho 2-24-16, Koganei, Tokyo 184, Japan (Received 30 September 1993; revised version received 14 December 1993; accepted 20 January 1994)
ABSTRACT A 3-D profiling system has been developed using a grating projection apj)r()ach and a one-step phase shift algorithm. In the system, a grating pailein projected on an object surface is acquired by a CCD camera, and the grating’s phase deformation caused by the surface shape is extracted by spatial phase shift processing, which uses only one frame of digital image. One phase value can be calculated from a successive 3-pixel range if the fringe period is set to either 4-pixel width or 8-pixel width. The fixed phase distribution of the system is excluded by a standard plane calibration. With 8-bit input data, the system’s RMS phase accuracy is developed up to 2jr/60 in the experimental examinations.
INTRODUCTION In the field of optical non-contact 3-D measurements, more and more researchers have turned to utilizing phase measuring techniques. The methods presented consist of techniques such as Fourier transform profilometry (FTP),’ phase-shift moire topography (PMT),2 ~sinusoidal spatial phase detection (SPD),5 multi-step phase-shift grating projection (MPGP),ON one-step phase-shift grating projection (OPGP),~etc. The moire method requires precise optical arrangements in order to form flat-contour fringes. It is quite difficult to apply the PMT methods to the crossed-optical-axes systems in non-telecentric cases because the contour fringes are not planar in that case. But by FTP, SPD, MPGP and OPGP methods, 3-D data can be directly extracted from a To whom correspondence should he addressed. 61
K Go, I. ~o.s1,,:a,to, 1.
(hoot
projected grating and the deformed grating phases can be separated into the objective distributions and the fixed optical distributions in the crossed-optical-axes occasions. Iherefore 3-D coordinates can be easily computed after the standard plane calibration of the fixed phase deformations. FIP and SPD are two excellent methods since only one image frame is required, but the phase calculations are complex and need a great deal of gray level data to extract oiie phase value. It IS necessary to use special digital processors to get high processing speed for the two methods. With MPGP approaches, if the mechanical niicro—shitting structure for obtaining several phase-shifted fringe patterns can ensure high locating precision. pure pixel-by-pixel resolution and high accuracy can be obtained. But the system accuracy is usually limited by the mechanical micro-shifter and the measuring speed is difficult to increase. By use of one-step phase shift projection methods OPGP), only one frame of grating image is required. optical crossed— axes deformation can he removed, the precise mechanical shifter is not demanded and all of the phase processes can he in binary-integer-form by using an arc—tangent look—up table. It can also achieve as high a computing accuracy as MPG P but its spatial resolution in the image plai~e of CCD is a little lower than MPGP because it needs at least three adjacent pixels’ gray data to calculate one phase value. It would seem that OPGP is one of the most advantageous approaches for realizing dynamic 3—I) measurements in pr;ictic;tl applications.
ONE-STEP PHASE SHIFT ALGORITHMS In a previous work the authors, referred to the idea ol spatial phase shift methods in inierlerometry ‘~ and established a one—step phase shil graling projection system. Recenik the~’ developed the system for crossed—optical—axes occasions and are now improving it with many new features. The known formula of the one—step phase shift method needs the pitch of input grating patterns to occupy approximately 4 pixels, so I hat the adjacent pixels on the lines vertical to fringe direction cati he regarded as having ~/2 phase shifting values between each other. lheielore the grating pattern can be written as follows if the fringe direction is along vertical rows of ((1) array 1(1/) \\ ide
i
/
tic
Thi
I)
liori,ontal
((/. I) dos
j
(/)(i,
j.:)
(I)
and \ei’tical pixel-counters. respectively.
3-D
surface
profilometrv
63
UI
L
Grating Projector
Reference Plane
w~DiA
CCD Camera Fig.
1.
Diagram
of
/
crossed-optical-axes grating projection measurements.
system
for
3-D
B(i,j) and C(i,j) are background and fringe contrast, respectively, and /(i, j, z) is the deformed grating phase which relates to pixel positions
and object height distributions. Using eqn (1), it can be derived that ~(i,j,
z)
=
[l(i+2,j)-l(i+1,j)] arctan ~ I(i,j) —I(i + 1,j)
~
(2)
where l(i, j), I(i + 1, j) and l(i + 2, j) are 3 adjacent pixels’ gray data, and eqn (2) should satisfy the condition that the background B(i,j) and contrast C(i, j) are not changing within a 3-pixel-range. The crossed-optical-axes grating projection system used is shown in Fig. 1, where the center of a projector’s exit pupil P and the center of a camera lens entrance pupil I determine the orientation of the reference plane’s X-axis, and the origin point 0 of the objective dimensions is at the intersection of the projector’s axis and the camera lens’ axis. The camera lens’ axis should be set vertical to the X-axis so that the CCD’s imaging dimensions can overlap with the reference plane dimensions. The z value is set as positive when the object surface is on the left side of the reference plane. The parameter L is the distance from both P and 1 to the reference plane, and the parameter w is the distance between P and I. Consider one light ray from point P which is projected onto a surface point A with height value —z. The ray PA goes through a point C in the reference plane and the imaging ray Al goes through a point D in the reference plane. The point A’ is the virtual intersection of the ray PA and the grating’s natural image-plane. All the points on the ray PA’ have the same fringe order of the projected
64
I?. Gu, 1. Yoshizawa, Y. Otani
grating. In the telecentric projection case, the ray expressing the same fringe order as the point A’ will he parallel to the axis PG and it will cross the reference plane at the point B. Physically, if two points on the reference plane in Fig. I express the same fringe order of the projected grating but are in different spatial positions, the distance between them is just the factor of the deformed grating’s phase. Therefore, if we regard the fringe periods as 4-pixel width or 8-pixel width, the grating pattern imaged on the reference plane can he written as follows 1(x, v)
=
B (x,
y)
+
((v v)
x cos [2irfx + 2irf0CD(z) + 21Tt)BCC
) + 2r~\Lv~
(3)
where B(v, v) and ((x, v) are background and fringe contrast. respec1 1 or (a is the pixel size equivalent to reference 4a So plane dimensions), CD(z ) is the phase deformation due to object’s tiyely,
f, equals
-
surface height variation and it is easy to find that: CD(:) =
(4)
~—~---
BC(x) is the phase deformation on reference plane due to crossedoptical-axes effects and it is solely dependent Ofl x values. And ~XIis the frequency difference between f~and the averaging practical frequency of acquired fringe patterns. Rewritting eqn (3) 2irfCD(z)
+
2.~f0BC(x)+ 2ir~f v
=
q~(x,V. ~)
(5)
and transforming x, v values to pixel-counting discrete forms, i.e,y ía and v = Jo (I, j are integers), we can change eqn (3) to the following form I(i,j)
When I,,
=
B(i,j) ~ C(i,j)cos [2irfai
+ ~(i,j,
~)J
(6)
(4 pixels per period), eqn (6) is same as eqii (I ). and the
=
4(1
phase /(i. j, :) can he extracted by eqn (2). When f~
,
i.e. the case
of S pixels per period, we can simply transform eqn (2) to the following (/)(1,J,z)=arctan
-
I(i+4,j)~I(i+2,/) I(i,~)~I(t+2,~) --~
—
~ i 4
(7)
This means that the 3 gray data for phase extraction are obtained from
3-D surface profilomerry
65
3 interval positions over a continuous 6-pixel range. If the object’s surface is smoothly changing, eqn (7) will work well, but the spatial resolution will decrease to half of that of the 4 pixel/period case. In order to keep the resolution as high as possible, we have derived a new formula to compute one phase value from a continuous 3-pixel range as follows (/)(i,
I, z) =
[
arctan I (V~ 1). —
L
j)
l(i j)—I(i+2 .
.
.
‘
.
21(z+1,j)—I(i,j)—I(i+2,j)i
1 .
I
Jr — —
4
i
(8)
Equation (8) keeps the same spatial resolution as that of the 4 pixel/period case. If the input picture has a high frequency background pattern or contrast variations, eqn (8) will work better than eqn (7) because eqn (8) uses gray data only in a 3-pixel range. In comparison with multi-step phase shift methods, the way of reducing measurement errors is quite different in one-step phase shift methods. For multi-step methods, the errors caused by mechanical shifting can be averaged by increasing the number of phase-shifted gray data. But in one-step occasions, if we use more gray data in one picture, a longer, continuous pixel-range is used and then the resolution in the imaging plane will be lower and more background noise effects will appear. Therefore the pixel-range for extracting one phase value should be shortest.
PHASE PROCESSING PROCEDURES Phase extraction We realized all the three algorithms shown in eqn (2), eqn (7) and eqn (8) for calculating phase maps from one frame of grating pattern. Depending on the input fringe’s frequency, either 4-pixel-period method or the 8-pixel-period method may be selected. The arc—tangent calculations of the 3 formulas are all executed by using a 16-bit look-up table. Therefore the computing time for extracting a 256 X 256 x 16-bit phase map which is wrapped —Jr to it is controlled to less than 1 second even using a PC/AT computer. It should be illustrated that a shifting constant-phase is attached to each phase value when using eqn (2) or eqn (7) or eqn (8). It is removed by anti-shifting computation while the arc—tangent calculation is performed in the system. This part of the computing program was written in PC Assembler language.
I? Go. i
66
Yoshi~ati’a, Y. ()tani
Phase unwrapping A serious problem for all the phase-shifting techniques is the phase unwrapping after the arc—tangent calculations. If each phase point is unwrapped along some fixed direction, it is always possible to propagate errors from had-quality points to good-quality points. In Ref. II Kwon reported briefly an excellent unwrapping method proposed h~ Shough. It is called ordered phase unwrapping (OPIJ). The idea of this method is that the ‘best’ point among the boundary points of an already-unwrapped area is selected for unwrapping. In our system, we referred to the idea and devoped an OPU procedure using PC Assembler language. For such processing, it needs a quality factor to determine what is the ‘best’ quality point. We made this factor Q, as: =
~
~
I~
--
(~)
Q,
is the quality factor of phase point P, and P,~, is one of 4 adjacent phase points around P, While getting the absolute phase differences, the 2ir jumps should he excluded if they exist. With the quality factor, an optimal order for unwrapping all phase points can he determined automatically and the propagation of errors will he kept as low as possible. The brief flow chart of the ordered phase unwrapping algorithm in our system is indicated in Fig. 2. where
Bad-region recognition In order to get results that only represent good phase points, a had-region recognition procedure has been developed. Ihe had-region recognition can be performed on the it to Jt wrapped phase map because the no-fringe areas and the poor-quality-fringe areas are most --
obvious here during the whole processing. The main feature of the had areas in wrapped phase maps is the phase values jumping randomly in high spatial frequency, so we selected a rectangular area around an inspected phase point and counted the surrounding points whose phase differences to the inspection point satisfy the following condition --
<2Jr
—
~,
(10)
where P is the phase point in the surrounding area, P, is the inspected phase point. and is an allowance of legal amplitude of phase variations. If the number counted is greater than a given threshold, the inspected point will he marked as a had point and will not he processed in phase unwrapping. Equation (10) also considers the case of 2ir jumps during the recognition. ~,
3-D surface profilometry
67
Select a start pixel on central column
1
Find neighbor pixels which are not unwrapped
Get the quality factors of selected pixels and add them in an ordered table if they are not in the table
Table length -l and pick out the best quality’s pixel from the table
If
Yes
_______
table length=-l
End
Unwrap the pixel with one of processed neighbors which has the best quality Fig. 2.
Flow chart of ordered phase unwrapping algorithm.
Correction of fixed optical phase deformations The phase distribution calculated either from eqn (2), eqn (7) or eqn (8) consists of the objective phase, the cross-optical-axes phase and the frequency fraction phase, as shown in eqn (5). Because the last two terms of eqn (5) are not influenced by object shapes, they will remain constant if the optical system is fixed. Therefore we can measure a standard plane as the zero-distribution of the system at first, and store its phase map as the system’s base data 4(x, y). Since the object height distribution is assumed to be zero here, we can get the base distribution from eqn (4) and eqn (5) as 4~(x,y)2Jrf,,BC(x)+2itt~f,,x
(11)
Parameters of eqn (11) are the same as illustrated in eqn (5). We do not need to know the theoretical form of BC(x) because its assignment has already been obtained pixel by pixel in the standard-plane phase map.
K Go, I
68
Yoshizawa.
Y Otani
After phase unwrapping for a measured object, the pure objective phase map can he extracted by subtracting the base phase map from it. Actually, the errors in optical arrangement are also included in the measured base phase map, so the objective phase map is the pure distribution referring to the standard plane. So far, the phase processing system is built and it can ensure the accurate phase measurements by using experimental calibrations of the standard plane. Improvements with grating tilt Originally, grating patterns were set vertical to the CCD’s scanning lines. But after processing, usually some vertical-line-type’s noise pattern with the order of 2it/lO appears on the phase map, and the frequency of it depends on the input fringe’s frequency. One reason for it is attributed to the arc—tangent calculation for phase extraction because it is hard to he accurate near ir/2 or —jr/2. This means that the phases near the two values are extremely sensitive to the A/D sampling errors of sinusoidal fringe shapes. Figure 3(a) shows a phase map of a flat surface using vertical fringe patterns. Since the fringe patterns are almost the same along the vertical rows of CCD array, the sampling errors will also he like that. If the errors on a row are serious, all the row’s phase values will he deformed and even if the phase map is modified by median-value-filter, the noise pattern is not obviously reduced, as shown in Fig. 3(h). If the RMS deviation of a planar phase map is about 2ir/3() originally, it will only he reduced to about 2.ir/50 after 3 X 3 median-value-filter. The noise pattern will confuse the fine distributions of the measured phases. Noting the fact that the fixed optical phases can he excluded by standard plane calibrations, we changed the projected grating tilt slightly on the Y-axis and calculated the phases along CCD’s scanning lines. The phase expressed in eqn (5) will change as q~(x,v, z) where
(Y
=
2iri~C’D(:)+ 2irf,BC(v) + 2irMv + 21rf%’ tan
c~
(L2)
is the angle between the Y-axis and the grating orientation,
and the other parameters are the same as mentioned in eqn (5). For the standard plane measurement, CD(:) = 0, the phase niap becomes (/)/.V,
v, z )
=
2itfBC’(x ) + 2ir~f1v+ 2itf,,v tan
(V
(13)
In eqn (13), rio parameter relates to object height: therefore the effects of grating tilt can also he removed by the standard plane calibration. The distinct feature between vertical grating projection and
69
3-D surface profilo~ne1ry
(a)
(h)
(c)
(d)
Fig. 3. The difference in the computed phase maps of a flat surface measurement comparing the vertical grating projection case with the tilted grating projection case. (a) Raw phase map in vertical grating case: (b) 3 >< 3 filtered result of (a); (c) raw phase map in tilted grating case; (d) 3 x 3 filtered result of (c).
tilted grating projection lies in the results of phase processing. By means of tilted grating projection, the fringe pattern sampled by CCD array is made with more delicate profiles, both in the X-axis and the Y-axis. The calculated phase map does not carry uniform error along vertical rows, and the error pattern is limited to periodic variations both in the X-axis and the Y-axis as shown in Fig. 3(c). Such a noise pattern can be removed by median-value filtering more successfully than in the vertical grating case. Experimental results show that if a planar phase map has the RMS deviation as 2ir/30, but its noise pattern is as in Fig. 3(c), its RMS deviation will be reduced to less than 2ir/100 by only 1 application of a 3 x 3 filter. Figure 3(d) shows the filtering result from
Fig. 3(c) and it is obvious that the noise pattern reduction in tilted grating cases is much quicker than in vertical grating cases.
71)
K
(.i,,
I. Yoshi~atta, 1 Otan,
Here it should he noted that the 3 X 3 median-value-filter processing on phase maps is reasonable because originally each phase value is extracted from 3 adjacent pixels’ gray data, i.e. it is obtained from the objective information in each small area corresponding to 3-pixels’ size. Therefore the system’s X—Y resolution is suitable to he set to the scale of 3 >< 3 pixels.
EXPERIMENTS In our experimental system. the measuring distance shown in Fig. I was set to I. = 110 mm, and the width between the points P and I was set to ii’ = 60 mm. The average grating pitch projected on the standard plane was measured at about 0~4mm and the measured size corresponding to 256 x 256 pixels was 24 mm X 24 mm. Therefore the equivalent pixelsize on the objective reference plane was about 0~094mm. The analog video signal of the deformed grating patterns captured by a C(’D camera was converted to a digital picture with 512 X 512 pixels and 8 hit gray-levels, and transferred to a PC486 microcomputer through an image processing board. All the processing software was written in C and Assembler languages and it was designed for the MS-Windows operating-environment, combining 256 colors or 256 gray-levels bitmap display. As an example, we measured a one hundred Yen Japanese coin which has a pattern varying in a 0’2 mm range. Figure 4 shows the
Fig. 4.
I he di~ii~il gritting pallet ii ol a coin ~viih 1)2 mm heighi \arlailons. It is acquired I morn ihe optical s~siem sho~~ n in Fig. I.
Fig. 5.
Ihe ora’ —coded phase disit thu ilon coiiipiiied tom I ig. 4. inien’.ti~ daLi h~ OI1L’—siej) phase ~.hmii algormihrns. I lie hrighiesi e~ci presenis phase .‘r .imid 11w darkest level pm esenis phase .‘r
3-D surface projilometrv
71
input fringe pattern with the size of 256 X 256 pixels selected from original 512 X 512 pixels digital image. The grating orientation was set with about 10 degrees tilt to the vertical row of the CCD array. Figure 5 shows the gray-coded display of the phase assignment that is wrapped in —Jr to it range corresponding to the range from the darkest gray-level to the brightest gray-level. It was obtained from Fig. 4’s gray data by the one-step phase shift processing, pixel by pixel. Figure 6 shows the bad-region distribution of Fig. S’s phase data and it was recognized by eqn (10) in the case of selecting 5 x 5 pixels as the surrounding area and setting the threshold to 12, i.e. half of the investigated points. The recognition means that the inspected point will be assigned as a bad point if its surrounding illegal-phase-variations are more than the threshold. Figure 7 shows three sequential intermediate situations during the ordered phase unwrapping processing for Fig. S’s phase data. The white areas represent the bad regions which will not be processed, the black areas are the legal areas which will be processed, and the gray areas represent the already unwrapped areas. The gray area grew automatically at the best quality point of the area’s boundary (at any one time), so it ensured the least error propagation in phase unwrapping and all the bad areas were passed regardless of any shape. Figure 8 is the phase map obtained after phase unwrapping, and it contains the object phase, the crossed-optical-axes phase, the frequency fraction phase and the grating tilt phase as expressed in eqn (13). Figure 9(a) shows the result after the standard plane phase calibration of Fig. 8, and the displayed pattern represents the gray-coded phase distribution of
Fig. 6.
Bad-regions of Fig. S’s phase map. They were recognized by expression (10) with 5 X 5 investigating rectangle.
72
1-? (i,,. I. Yo~1,,:oima.
I
~
(l~
( )iiini
}
~ ~
I
~
C
t
‘I..,
~
.‘~:..
~.
.
/
(h)
(a)
Fig. 7.
I lilL IrrlLriiiLdr.lt. ~rl1I,iiI,os at ItL oideicd ph Is~ uo~~rrppiiielI~sr’.’. Ii I ~ pIi,r~~ ~t ii I I IR v~liii. rc,~iori’l~pi~‘LII) the had—mcr1iorl’ \~hich \\crc riot pos.~.‘‘cd lw hi ILk id’. ri. h~ ic,~i1 ilL I’. ~~hrch ~siII tic un~smappLd.ihL ~ ~rtIcrrr i’. IlL .iiic ids tiii’~i rppLd Iic.i ~~hicir gmr’\\s it the hcst—qualit~point oh Ii’ isruod IS each huh rid lu. iiio~sirie‘.hipL. i~shown us (a) —(e) ru ‘~cquciIcc.
the measured coin. Figure 9(b) is the wire—frame plot of Fig. 9( a pl~se distribution. The phase values wrapped it to it are kept in the computer’s memory as a I fi hit binary—integer—form, and after unwrapping es cr\ —
phase point keeps a 24 hit-length in computer memories or disk tiles. We made an accuracy examination for the whole system b~locating a plane surface in the objective space afldl nieasuring its phase distribu— lions at several distances. Through the statistical analyses of all the 251) 25( ) phase points in the center ot the phase maps, s~ e got the ~‘
3-D surface profllometry
Fig. 8.
73
The phase map after phase unwrapping. It contains the object phase and the
fixed optical phase.
~i) (a)
• (b)
Fig. 9.
The object phase map extracted from the standard plane calibration of Fig. 8’s phase data: (a) is the gray-coded display of the phase map and (h) is its 3-D wire-frame
plotting chart. The changing scale of the surface pattern was measured in 02 mm.
distribution of RMS errors along the Z-axis as shown in Fig. 10. The RMS values near the standard plane position were 2,r/160 (in ±FOOmmrange), and up to ±S~00mm from the standard position the RMS deviations were 2ir/60. The reduction of the measurement precision is mainly due to the optical distortion of the projector and the camera lens in defocused imaging occasions. The optical distortion makes the plane surface data increasingly curved when the defocused value increases. Considering that a 2ir value of phase is equivalent to one fringe period p and the P—V deviation is regarded as 3 times the RMS deviation, the surface measuring accuracy can he determined as RMS X 3 X P/2ir. So the system accuracy can he within 20 p~mof the
74
Ii’. (iii, I. Yoshi:aoa. }. Utani 0.040
i~’0.032 2.
•~,r..r,.,,u
-
-.
0.024
-6.00
-4.00
-2.00
0.00
2.00
4.00
600
Plane position along Z axis (Unit: mm) Fig. 10. lhc experimental RMS—error drsirrhutiiin Ii) the whole sssteiii. I) SS:I~ ohiarned by plane surface measurements with scs cral displacements in the stamldIurLt plane ioe~uiiofl.
4
Fig. 11. Ihe 3—I) wire—trame plot of the phase measuremeni at a sample surtace wiiii two known steps. In the result, the high step represented 100 /Lm varration and the ass siep represented So ~ variation.
systeni configurations such that p equals 0’4 mm: the measuring depth is 10 mm and the measuring size is 24 mm X 24 mm. ‘l’o demonstrate the accuracy, we measured a sample surface with known 5)) ~rm and tOO ,um steps on it. Figure 11 shows the wire-frame plotting result of its phase map. The high step represents the 100Mm and the low step
represents the 50
jwi.
3-D surface profllornetry
75
CONCLUSION We have developed the one-step phase shift techniques into 3-D surface measurements, and the one-step phase shifting computation can be executed in 3-pixel-range in the cases of either 4 pixels per fringeperiod or 8 pixels per fringe-period. The areas without fringes or with poor quality phase data can be recognized as bad regions by investigating the irregularity of phase distributions. An optimal method called ordered phase unwrapping was adopted to remove the 2ir jumps in arc—tangent calculated phase maps. The unwrapping procedure always seeks the best-quality point among the boundary phase data of processed areas. It can go around any type of bad region and ensure the least error-propagation in phase unwrapping. The fixed optical phase distributions including crossed-optical-axes phase deformations and system dislocating deformations, can be excluded by a standard plane calibration. This ensures the pure objective phase extraction by experimental techniques, but not by computing calibrations. This is significant in practical applications. In order to restrain the regular noise pattern corresponding to input fringe patterns, the fringe orientation of the projecting grating is set with a small tilt to the CCD array’s vertical rows. In this way, the noise pattern can be controlled both in horizontal and vertical direstions, and can be further reduced by 3 x 3 median-value filtering more rapidly than in the case of vertical gratings. The influence of the grating tilt can also be excluded by standard plane calibrations. Experimental results showed that the RMS accuracy of the phase processing improved up to 2ir/60 in the objective dimension of 24mmX2SmmX 10mm. REFERENCES 1.
2. 3. 4. 5. 6. 7. 8.
9. 10. 11.
Takeda, M. & Mutoh, K., App!. Opt. 22 (1983) 3977. Reid, G. T., Rixon, R. C. & Messer, H. I., Opt. Laser Tech., 16 (1984) 315. Gu, R. & Zhou, R., In Proc. SPIE, Vol. 1230, 1990, 638. Yoshizawa, T. & Tomisawa, T., Opt. Eng. 32 (1993) 1668. Toyooka, S. & Iwaasa, Y., App!. Opt. 25 (1986) 1630. Srinivasan, V., Liu, H. C. & Halioua, M., App!. Opt. 24 (1985) 185. Komatsubara, R. & Yoshizawa, T., J. Japan Soc. Precision Eng. 55 (1989) 85 (in Japanese). Su. X., Bally, G. V. & Vukicevic, D., Opt. Comm. 98 (1993) 141. Gu, R. & Yoshizawa, T., Proc. SPIE, Vol. 1720, 1992, 470. Shough, D. M., Kwon, 0. Y. & Leary, D. F., In Proc. SPIE, Vol. 1221, 1990, 394. Kwon, 0. Y., In Proc. SPIE, Vol. 816, 1987. 196.