Experimental study on subaperture testing with iterative triangulation algorithm

Applying the iterative triangulation stitching algorithm, we provide an experimental demonstration by testing a Φ120mm flat mirror, a Φ1450mm off-axis parabolic mirror and a convex hyperboloid mirror. By comparing the stitching results with the self-examine subaperture, it shows that the reconstruction results are in consistent with that of the subaperture testing. As all the experiments are conducted with a 5-dof adjustment platform with big adjustment errors, it proves that using the above mentioned algorithm, the subaperture stitching can be easily performed without a precise positioning system. In addition, with the algorithm, we accomplish the coordinate unification between the testing and processing that makes it possible to guide the processing by the stitching result. ©2013 Optical Society of America OCIS codes: (120.3180) Interferometry; (120.6650) Surface measurements; (220.4610) Optical fabrication; (220.4840) Testing. References and links 1. J. G. Thunen and O. Y. Kwon, “Full aperture testing with subaperture test optics,” Proc. SPIE 351, 19–27 (1983). 2. W. W. Chow and G. N. Lawrence, “Method for subaperture testing interferogram reduction,” Opt. Lett. 8(9), 468–470 (1983). 3. T. W. Stuhlinger, “Subaperture optical testing: experimental verification,” Proc. SPIE 656, 118–127 (1986). 4. A. Kulawiec, P. Murphy, and M. D. Marco, “Measurement of high-departure aspheres using subaperture stitching with the Variable Optical Null (VON),” Proc. SPIE 7655, 765512-1–765512-4 (2010). 5. C. Supranowitz, C. M. Fee, and P. Murphy, “Asphere metrology using variable optical null technology,” Proc. SPIE 8416, 841604-1– 841604-5 (2012). 6. P. Su, J. H. Burge, and R. E. Parks, “Application of maximum likelihood reconstruction of subaperture data for measurement of large flat mirrors,” Appl. Opt. 49(1), 21–31 (2010). 7. P. Su, Absolute Measurements of Large Mirrors (The University Of Arizona 2008). 8. D. Liu, Y. Yang, C. Tian, Y. Luo, and L. Wang, “Practical methods for retrace error correction in nonnull aspheric testing,” Opt. Express 17(9), 7025–7035 (2009). 9. M. Otsubo, K. Okada, and J. Tsujiuchi, “Measurement of large plane surface shapes by connecting smallaperture interferograms,” Opt. Eng. 33(2), 608–613 (1994). 10. P. F. Zhang, H. Zhao, X. Zhou, and J. J. Li, “Sub-aperture stitching interferometry using stereovision positioning technique,” Opt. Express 18(14), 15216–15222 (2010). 11. M. Novak, C. Zhao, and J. H. Burge, “Distortion mapping correction in aspheric null testing,” Proc. SPIE 7063, 1–8 (2008). 12. C. Zhao, R. A. Sprowl, M. Bray, and J. H. Burge, “Figure measurement of a large optical flat with a Fizeau interferometer and stitching technique,” Proc. SPIE 6293, 1–9 (2006). 13. S. Y. Chen, S. Y. Li, Y. F. Dai, L. Y. Ding, and S. Y. Zeng, “Experimental study on subaperture testing with iterative stitching algorithm,” Opt. Express 16(7), 4760–4765 (2008).


Introduction
Subaperture stitching (SAS) was primarily developed to overcome the aperture size limitations and slope sampling limitations of conventional interferometer, where the primary goal is to obtain a full aperture map from subaperture measurements without having to test the whole part at one time.It's useful especially for large flat mirrors, spherical surfaces with high numerical aperture, large convex surfaces, and aspheric surfaces exceeding the vertical range of the interferometer.
The SAS algorithm has been playing a very important role in the application of the subaperture testing (SAT) method.Many stitching algorithms has been provided and obvious improvements can be observed from the Kwon-Thunen method [1] and the simultaneous fit method [2], to the discrete phase method [3],the QED's variable optical null technique [4,5], and then to Arizona's method with maximum likelihood algorithm considering error of the reference mirror [6,7].Aiming to improving the efficiency of stitching, unifying the processing and testing results well and testing large surfaces with big errors, we proposed a kind of algorithm, iterative triangulation stitching algorithm and a method unifying the processing and testing results.Advantages of the algorithm were claimed and verified through experiments and the performance of the method unifying the processing and testing results was also shown in the experiments.
In this paper, we focus on the iterative triangulation stitching algorithm and experimental demonstration of the high accuracy surface shape measurement and the alignment between the processing and testing maps with the iterative method.The stitching algorithm has been applied for a Φ120mm flat mirror, a Φ1450mm off-axis parabolic mirror and a convex hyperboloid mirror.In the testing, if the full map of the mirror can be tested at one time, the stitching result can be compared with the full aperture testing result to evaluate the stitching accuracy.However in actual measurement, full aperture testing can hardly be achieved especially for large plane mirror or convex mirror.In this situation, we propose a method to evaluate the stitching accuracy by comparing the stitching result with the subaperture testing result which is used to examine the stitching accuracy.This kind of method is also applied in the following experiments.This paper is organized as follows.In Section 2, the basic theory of iterative triangulation stitching algorithm and the method used to unify the coordinates between the processing and testing results are introduced.In Section 3, we apply the above algorithm and method to the actual experiment and the relative testing result is introduced.The conclusion is given in Section 4.

The iterative triangulation algorithm
The stitching algorithm is shown in Fig. 1.First we take interferometric measurements and correct distortion to each subaperture, choose a subaperture as a standard and put all the testing data in a global coordinates according to their relative positions.If the SAT is taken in condition of nonnull testing, nonnull testing error should be removed from each SAT result [8].Then we define grid points in the global coordinate and calculate the phase value of the grid points in the overlapping area between subapertures according to the Delaunay triangulation as shown in 2.1.1.After getting the phase value of the grid points in the overlapping area, we calculate the stitching coefficients of each subaperture except the standard one according to our method extended from the one put forward by Masashi Otsubo [9] as shown in 2.1.2.Then we will get the residual map of every two adjacent maps after removing their adjustment terms respectively in their overlapping area.And then we evaluate whether the RMS of the residual meets the requirement.If not, we will calculate the twodimensional cross-correlation between every two adjacent subapertures to get a more accurate positioning to recalculate the stitching coefficients of each subaperture as shown in 2.1.3until the residual meets the requirement.
Finally a full-aperture map is obtained by connecting each subaperture testing result with stitching coefficients removed from its own testing result and the data in the overlapping area is calculated with a stitching factor as shown in 2.14.

Delaunay triangulation interpolation
In actual testing, coordinates of subapertures can be unified in a global coordinate according to the mechanic alignment accuracy or with marked points, targets and so on.Before we calculate the stitching coefficients of each subaperture, the overlapping correspondence in the overlapping area should be found in advance to calculate the stitching coefficients of each subaperture [10].If we assume that all phase data have been transformed to their corresponding global 3D coordinates, the overlapping correspondence problem can be simplified as follows: Before stitching, grid points which represent the same points between subapertures should be defined.As shown in Fig. 2(a), a uniform grid is defined on the X-Y plane.Then, all of the points in the i th and the j th subapertures are projected onto the X-Y plane.For the sake of conciseness, only the grid points falling in the overlapping area are shown as dot points in Fig. 2(a).The coordinates of the dot points in the z direction need to be calculated respectively in each subaperture.Consider the grid points in the overlapping area between the i th and the j th subapertures.The testing result of the i th subaperture is shown in Fig. 2(b) where the triangle is the Delaunay triangulation result to the phase data of the i th subaperture.The triangulation result can be described by Eq. ( 1) where which is used to describe the plane equation of the triangle (the coefficient of z won't be zero).
The dot points in the overlapping area can be found in one triangle and the coordinates of these dot points in the Z direction can be interpolated from the plane function according to the data from the i th subaperture by Eq. ( 2).In the same way, the coordinates of these dot points can also be got in the Z direction according to the testing data from the j th subaperture.
Because the phase distribution in the overlapping area between the adjacent subapertures should be the same, the coefficients of the subapertures can be solved with the method described in Section 2.1.2on the basis of the coordinates of these dot points in each subaperture.

Calculation of stitching coefficients
The content in this part is extended from the algorithm developed by Masashi Otsubo [9].We assume that there are N subaperture measurements and the coordinates of grid points in the z direction have been calculated in each subaperture.For convenience, the N th subaperture is taken as a standard and the map shape of the i th subaperture can be expressed as, where ( , ) i x y Φ is the testing result of the i th subaperture and ( , ) k f x y can be any predefined functions while ik a is the relative coefficient to be fit to each ( , ) k f x y of every subaperture except the standard one.L is the number of terms to be fit.When we perform the stitching to a plane mirror, usually the fitting functions are limited to tip/tilt and piston and the function ( , ) k f x y can be written as: ( , ) ; ( , ) 1.
The algorithm mentioned in the manuscript can be generalized when more terms need to be fit and the fitting functions for each subaperture may not be limited to tip/tilt and piston.
Least squared fitting is used: Equation ( 5) can be transformed to a group of linear equations in form of, P, Q and R are defined as follows, P is a vector in length of ( 1) N L − × row.For convenience, ( 1) ) j L k − ⋅ + th row of the vector.j represents the sequence number of the subaperture while k is the sequence number of predefined functions.That means, and ( 1) ( , )( ( , ) ( , )) Q is a matrix in size of ( 1) N L − × and Q ((j-1)⋅k)((l-1)⋅k') is used to represent the element in the row (( 1) of the matrix.l also represents the sequence number of the subaperture while k′ means the sequence number of predefined functions.Then, and Q ((j-1)⋅k)((l-1)⋅k') can be written as, ( , ) ( , ), ; ( , ) ( , ), .
of the vector and ( 1) where jk a is the coefficient of the predefined fitting function ( , ) k f x y for the j th subaperture.
By solving the linear equations we can obtain the fitting coefficients ik a , and then we can combine the subaperture measurements together to get the full aperture map.In the overlapping areas, a kind of stitching factor is applied to calculate the phase data which is introduced in 2.1.4.

Position calculation of each subaperture with two-dimensional cross-correlation
Before applying the stitching factor described in Section 2.1.4to the subaperture measurements, we need to find out whether the stitching accuracy meets our requirement.In theory, the phase data should be the same in the overlapping area after removing the adjustment terms.So we calculate the residual map of every two adjacent maps to judge whether the RMS of the residual meet the requirement.If it meets the requirement, the full map can be obtained as shown in 2.1.4.If not, the relative position between adjacent subapertures can be described as, ; .
x x dx x y y dy y where 1 1 2 2 ( , ) ( , ) x y x y and are the coordinates of pixels in sub1 and sub2 respectively, ( , ) dx dy the approximate position relationship between them and ( ) Δx,Δy the accurate position relationship to be solved.

( )
Δx,Δy can be calculated by calculating the two-dimensional cross-correlation of the data in sub1 and sub2 as ** F G with the position of the peak ( ) Δx,Δy .

Stitching factor in the overlapping area
Usually average is introduced to all the subapertures which have valid data over this area, however, this value way will leave small steps between adjacent subapertures.
Step is the sudden change of the data in the map which makes the map unsmooth and discontinuous.It usually appears as a line at the edge of the subaperture in the stitching map.To remove the step, low-frequency filter can be used in the full aperture map.However, this method also takes away the high-frequency message of the stitching result, not only in the overlapping area, but also in the non-overlapping area.A kind of stitching factor used in the overlapping area is proposed to overcome the dilemma.
For the sake of conciseness, only two subapertures are considered in the modal as shown in Fig. 3.In Fig. 3, O 1 and O 2 are the centers of each subaperture and P is the point to be calculated in the overlapping area between the subaperture 1 and 2.
In the light of the idea of Lagrange interpolation, we assume that the data approaching the center of the subaperture is more accurate.The phase data of point P can be written as, where 1 z and 2 z are the measurement results of point P in the subaperture 1 and 2 respectively.

( , )
O O x y are the global coordinates of the points P, 0 1 and 0 2 respectively. means the module of the vector inside.
In the general case, considering that the point P is in the overlapping area of N subapertures, in a similar way to Eq. ( 13), the phase date of point P can be written as, where k z is the phase data of the point P in the k th subaperture and k w is the relative weighting factor which can be written as, 1 ( 1 ) For any point in the overlapping area, the sum of the stitching factors should be 1.So Eq. ( 15) should be written as, Based on the Eq. ( 17), the phase data in the overlapping area can be calculated and relative experiment results are described in Part 3.

Alignment between processing and testing coordinates
In the actual processing, testing results are usually used to guide processing.However, the testing results are in the testing coordinate which is in unit of pixel while the processing coordinate is in unit of millimeter.The alignment needs to be accomplished between these two kinds of coordinates and then testing and processing can be unified.Marked points are used to accomplish the alignment as shown in Fig. 4. In actual SAT, distortion should be removed from each phase map before stitching [11,12].After distortion correction, the coordinates of targets in the testing results can be described as ( , ), 1,2,3 , i i x y i =  while the relative coordinates in the processing coordinate system can be written as ( , ), 1,2,3 where ( , ) x y t t is the relative translation between the coordinates of testing and processing and θ is the relative rotation degree between them.s is the magnification to each pixel in the testing results.After correction of distortion to each phase map, s should be the same to every point.
Then the relationship between the coordinates before and after the transformation should be: where ' '

( , )
i i X Y is the coordinate of pixel ( , )   i i x y after alignment between testing and processing.By minimizing the sum, ' where N is the number of marked points, matrix T can be solved and the relationship between the processing and testing coordinates can be obtained.

Φ120mm flat mirror
In this experiment, the clear aperture of the flat mirror is about 120mm.The SAT is carried out with a Φ150mm interferometer and 5-dof adjustment platform [Fig.5].Both translations (X, Y and Z) and rotational tables (yaw and pitch) are adjusted manually.As a verification experiment, four subapertures are measured with the interferometer (up, down, left and right subapertures).The SAT results are shown in Fig. 6.At the same time, the full aperture testing to the flat mirror can be accomplished with the same interferometer.Figure 7 gives the full aperture testing result.In the experiment, the testing accuracy is within 2/1000λ.The alignments between subapertures in both X and Y directions are accomplished with the marked point.Figures 8 and 9 give the stitching results with traditional stitching algorithm [9] and the iterative triangulation algorithm respectively.In the experiment, the stopping criterion of the stitching is that the RMS of the residual map between every two adjacent subapertures is less than 1.5nm.After 3 circles of iteration, the relative stitching result is calculated.To compare the stitching results between two stitching algorithms, we calculate the residual maps between the full aperture testing result and the stitching results with relative algorithms respectively.By subtracting the data between stitching result and full aperture testing result point by point, the residual map can be got.Figures 10 and 11 give the residual maps with these two kinds of algorithm.From Fig. 10, obvious steps in the residual map calculated with traditional stitching algorithm can be observed at the edge of the subaperture.However, the residual errors between the full aperture testing result and the iterative triangulation algorithm stitching result (Fig. 11) seem to be uniformly distributed, and there is no step in the residual map which means that the stitching result is smooth and continuous.
The PV and RMS errors of the residual map with tradition algorithm are 0.038λ and 0.003λrespectively while the PV and RMS errors of the residual map with the iterative triangulation algorithm are 0.020λ and 0.002λ respectively which are better than the traditional method.At the same time, the PV and RMS errors of the residual map make the stitching algorithm and experiment result credible.

Φ1450mm off-axis parabolic mirror
In this experiment, the clear aperture of the off-axis parabolic mirror is about 1450 mm and compensator is used to accomplish the test.As the mirror is in the roughly polishing period, there are relatively big figure errors in the mirror surface measurement, and local interference fringes are so densely packed in the full aperture test that the map of the mirror can't be figured out analytically very well.However, interference fringes can be observed more clearly with local amplification.So, the goal of full aperture map reconstruction can be accomplished by local amplification stitching.A simple 5-dof platform is built for the subaperture testing of the mirror as shown in Fig. 12.Both translations (X, Y and Z) and rotational tables (yaw and pitch) are adjusted manually.The alignments between subapertures in both X and Y directions are accomplished with targets.In the experiment, the testing accuracy is within 3/1000λ.In total, 3 subapertures (up, middle and down) are tested and the measured results are shown in Fig. 13.After correcting distortion to each subaperture, according to the position of targets and parameters of optical test, based on the above algorithm, the stitching result is obtained as shown in Fig. 14 with PV and RMS values respectively.The stopping criterion of the stitching is that the RMS of the residual map between every two adjacent subapertures is less than 2nm.After 9 circles of iteration, the relative stitching result is calculated as in Fig. 14.There is an impressive advantage which can be got from the stitching result.In Fig. 14, eight regular circles are distributed on the stitching map.Since the targets are not reflected, the testing result from interferometer is shown in Fig. 15.The data on edge of the target is a little lower than the data outside it, which leads to the target shown in the stitching map as a regular circle with the stitching algorithm mentioned above.Figure 16 gives out the stitching map after replacing the data in the target with interpolation.Table 1.Coordinates of targets in the testing and processing coordinate system respectively.
Fig. 18.Surface map after different processing cycles.
According to the alignment result with above method, after different processing cycles, the surface map of the mirror are shown as in Fig. 18.It can be got from Fig. 18 that both the PV and RMS become better after a new processing cycle which proves that the processing is convergent and the alignment result is convincing.Usually the stitching accuracy can be evaluated by comparing the stitching result with the full aperture testing result as shown in Section 3.1 [13].However sometimes it is difficult to get the full map at one time especially for large flat mirror and convex mirror.In this case, to examine the stitching accuracy, we choose a subaperture different from the subapertures used for stitching as a self-examine subaperture [Fig.19(a)].After distortion correction to the subaperture, the residual map between stitching and self-examine subaperture maps is shown in [Fig.19(b)].The RMS error of the residual map is 0.003λ.

Φ130mm convex hyperboloid mirror
In this experiment, the clear aperture of the convex hyperboloid mirror is about 130mm.Conic constant κ is −1.8128 and the vertex radius of curvature R is 1227.65mm.A simple 5dof platform is built for the SAT of the mirror as shown in Fig. 20 and the alignments between subapertures in both X and Y directions are accomplished with targets.In the experiment, the testing accuracy is within 3/1000λ.A Φ150mm aperture interferometer and a standard transmission sphere with F = 11 are chosen to accomplish the SAT.The radius of each subaperture r is about: Five subapertures are designed to finish the full aperture testing as shown in Fig. 21 and the measured results are shown in Fig. 22.In the measurement to the convex hyperboloid mirror with a standard transmission sphere, extra aberrations (nonnull testing error) will be added to the testing result [8].For the central subaperture, the extra aberrations behave as a combination of power and spherical aberrations (Fig. 23).For the off-axis subaperture, take subaperture 2 as an example.The behavior of the extra aberrations is shown in Fig. 24.In the stitching, the stopping criterion is that the RMS of the residual map between every two adjacent subapertures is less than 4nm.After 17 times of iteration, the stitching result is shown in Fig. 25. Figure 26 gives out the stitching map after replacing the data in the targets with interpolation.compensation using iterative triangulation algorithm and a method accomplishing the alignment between processing coordinate system and testing coordinate system.To evaluate the stitching accuracy for mirrors that can't be tested in full aperture one time, self-examine subaperture can be measured to get the residual map as shown above to fulfill the stitching accuracy evaluation.It can be proved from the above experimental results that the reported method can obtain the reconstructed full-aperture surface map with satisfactory accuracy without precise mechanical adjustment device.As the experimental study is only for mild aspheres now, further theoretical and experimental developments of the reported method on strong aspheres, especially strong convex aspheres and freeform surfaces will be taken to accomplish the SAT to such surfaces with or without any compensation.

Fig. 2 .
Fig. 2. Calculation in the overlapping area.(a) sketch of the projection on the X-Y plane; (b) sketch of Delaunay triangulation.
y , z ), (x , y , z ) and 3 3 3 (x , y , z ) are the coordinates of the three points making up of the fundamental triangle in the triangulation.a, b and d are coefficients of the equation, z ax by d = − − −

Fig. 3 .
Fig. 3. Schematic diagram of calculation with the stitching factor.

Fig. 4 .
Fig. 4. The configuration for alignment and correction of mapping distortion.

Fig. 10 .
Fig. 10.Residual map between traditional algorithm and full aperture testing result.

Fig. 11 .
Fig. 11.Residual map between iterative algorithm and full aperture testing result.

Fig. 16 .
Fig. 16.Stitching result after replacing the data in the targets.
(a) self-examine subaperture (b) residual map (a) surface map after one processing cycle (b) surface map after two processing cycles.

Fig. 26 .
Fig. 26.Stitching result after replacing the data in the targets.As the same with Section 3.2, the processing coordinates of the targets are got from Three Coordinate Measuring Machine as shown in Fig. 27 and the relative testing coordinates of the targets can be read from the stitching result [Fig.25].Coordinates of targets in the processing and testing coordinate systems are shown in Table 2.The alignment coefficients can be found by Eqs.(20-21) and the results are as follows: 0.1258, 64.8861mm, 64.9580mm, 0.0745 .