- Email: [email protected]
Stitching error calculation method for the slope data on a square subaperture
Optics Communications 380 (2016) 124–133
Contents lists available at ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
Stitching error calculation method for the slope data on a square subaperture Guangrao Guo a, Dahai Li a,n, Kewei E a, Mengyang Li a, Haiping Chen b,n a b
School of Electronics and Information Engineering, Sichuan University, Chengdu 610065, China Research Center of Laser Fusion, China Academy of Engineering Physics, Mianyang 621900, China
art ic l e i nf o
a b s t r a c t
Article history: Received 30 November 2015 Received in revised form 9 May 2016 Accepted 15 May 2016
The data on subaperture in this paper is collected by a Hartmann test device or a Hartmann–Shack wavefront sensor. The slope data on each subaperture of the test optical component is obtained via the linear least squares fitting routine in this paper. Importantly, the error transfer formulas are deduced when the measurement data are polluted by the random noise, which will affect the uncertainty of stitching parameters, such as tilt or defocus. In particular, the accuracy formulas for stitching evaluation, which can be used to calculate the stitching error of each point on each subaperture for both the parallel mode and the serial mode, are derived mathematically. Therefore, this paper provides these formulas to estimate exactly the stitching error in each subaperture if the variance of random noise has been estimated and if the overlapping area is also given in advance. The results of simulation experiments show that the stitching accuracy formulas proposed are verified and they can be used to evaluate stitching accuracy. The reason why the error of the parallel stitching mode is less than that of the serial stitching mode is presented theoretically. & 2016 Elsevier B.V. All rights reserved.
Keywords: Subaperture Slope stitching Parallel stitching Serial stitching Large component Accuracy
1. Introduction With the development of science and technology, large optical components are widely used in space telescopes [1,2], laser fusion systems [3,4] and other aspects. In traditional optical measurements, the large optical surface is tested by a large-scale interferometer. However, the interference testing method requires a flat reference, and most interferometers are sensitive to the measurement environment. As a result, the cost of a large-scale interferometer is greatly increased, and good environment control is required for testing [5,6]. In 1982, the subaperture measuring method was first introduced in interferometric metrology by C.J. Kim at the Optical Sciences Center at the University of Arizona, [7]. Since the method was proposed, it has been widely used in the measurement of large aperture optical elements. More recently, subaperture stitching interferometry is becoming increasingly popular. In 1994, Otsubo proposed a global subaperture stitching method based on stitching parameters; in this method, the misalignment errors is simultaneously fitted and obtained over multiple overlapping subapertures [8]. In 2003, QED Technologies successfully developed an automated subaperture stitching interferometer. This interferometer can successfully test the surface of flat and spherical optical elements with n
Corresponding author. E-mail addresses: [email protected] (D. Li), [email protected] (H. Chen).
http://dx.doi.org/10.1016/j.optcom.2016.05.043 0030-4018/& 2016 Elsevier B.V. All rights reserved.
diameters of less than 200 mm [9]. Moreover, research studies have also focused on the analysis of the accuracy of subaperture stitching. To date, there are three types of methods to evaluate the accuracy of subaperture stitching technique. In the first type, the stitching results are compared with the results of the full aperture interferometry; this method of evaluation is the most simple and straightforward one [10]. In the second type, the effect of random noise and overlapping area on the stitching parameters is qualitatively analyzed. This method can only obtain the error of the stitching parameters between two adjacent subapertures and it cannot obtain the error of each point on each subaperture after the subaperture stitching procedure is completed [8,11]. In the third type, according to the knowledge of statistics, the proportion of the regression sum of squares of the phase difference to the total sum of squares of the phase difference is used to evaluate the accuracy of the stitching [12]. Different from the traditional subaperture stitching method based on wavefront measurements, the subaperture wavefront slope stitching test is based on a Hartmann–Shack wavefront sensor [13]. Hartmann–Shack sensors have the advantages of insensitivity to the air disturbance and vibration, rapid testing and high dynamic range. In addition, subaperture wavefront slope stitching tests can avoid the influence of the piston. Thus, this method has a unique advantage in the testing of a large aperture optical element. Due to these characteristics, an increasing number of scholars are studying the subaperture wavefront slope stitching testing method [14–16]. Currently, research studies on
G. Guo et al. / Optics Communications 380 (2016) 124–133
wavefront slope stitching are mainly focused on the establishment of the mathematical mode and qualitative analysis of the stitching error. However, no mathematical formula is available to evaluate stitching accuracy when random noise exists in measurement process. In this paper, we describe in detail the following: 1) the effect of random noise and overlapping area on the stitching parameters and 2) the error relationship between the stitching parameters and the adjusted slope. Importantly, the error transfer formula for the slope stitching is first determined. On this basis, the respective formulas for solving the accuracy of serial slope stitching mode and parallel slope stitching mode are derived. This paper shows that these formulas are verified by our simulation experiments. The RMS and PV values of the residual error for the reconstruction of the wavefront based on a given slope matrix in our simulation experiments demonstrate that the error of parallel stitching is statistically less than that of serial stitching.
2. Principles As mentioned in the introduction, a large-size component can be measured using the stitching method. In this manner, we can obtain a wider range of measurement area with good resolution using a HS with a limited aperture. The disadvantages of the stitching methods include error propagation, and we must move the sensor or the sample repeatedly to cover the entire surface. Therefore, the slope data on a large aperture consists of multiple slope matrices on different sub-regions. Due to the misalignment, the piston, tilt, and even defocus, will be inevitably introduced into the measurement data when we move the sensor from one place to another. To eliminate the interference of misalignment, the slope data on the overlapping area between two adjacent subapertures are used to obtain the corresponding stitching parameters by a linear least squares fitting routine. As a result, these parameters, such as coefficient of tilt or defocus, are obtained, and then this procedure is repeated to the end of the stitching task until all these parameters are obtained. Next, these parameters are treated to construct the stitched slope data on the whole component. For simplicity of discussion, as shown in Fig. 1, we assume that the two adjacent subapertures W1 and W2 have a common area due to their overlapping area, and the measured phase distributions over the two areas are denoted as Φ1 and Φ 2, respectively. The connection of the two adjacent areas is considered first. Let W1 be the reference plane, and let the relationship between the phase distributions Φ1 and Φ 2 in the common area fulfill the following equation:
Φ2 = Φ1 + P + Tx x + Ty y,
(1)
where P is the piston, and Tx and Ty are the coefficients of the tilt in the x and y directions, respectively. Considering the derivation of x -tilt and y -tilt in Eq. (1), we have
G2x − G1x = Tx G2y − G1y = Ty,
125
(2)
where G1x = ∂Φ1/∂x and G1y = ∂Φ1/∂y are the x -tilt and y -tilt measured slopes, respectively, over the area W1, and G2x = ∂Φ 2/∂x and G2y = ∂Φ 2/∂y are the x -tilt and y -tilt measured slopes over the area W2, respectively. The parameter P disappears in Eq. (2) after taking the derivative. In other words, the piston does not influence the measured slope on subapertures; thus, we only need to find the coefficient Tx of x -tilt and coefficient Ty of y -tilt. In general, a minimum number of sampling points in the common area are chosen to accomplish the calculation of the values of Tx and Ty within an acceptable stitching error because the existence of random noise during measurement will affect the calculation accuracy of parameters. Suppose there are n sampling points in their common area; in this case, we can obtain the corresponding equations from Eq. (2) as follows:
⎡ G2x (x1,y ) − G1x (x1, y ) ⎤ ⎡ ⎤ 1 1 ⎥ ⎢ 1⎥ ⎢ ⎢ G2x (x2,y2 ) − G1x (x2 , y2 ) ⎥ = ⎢ 1⎥ ⎡⎣ T ⎤⎦, ⎥ ⎢ ⋮⎥ x ⎢ ⋯ ⎥ ⎢ ⎥ ⎢ ⎣ G2x (xn , yn ) − G1x (xn , yn )⎦ ⎣ 1⎦
(3a)
⎡ G2y (x1,y ) − G1y (x1, y ) ⎤ ⎡ ⎤ 1 1 ⎥ ⎢ 1⎥ ⎢ ⎢ G2y (x2,y2 ) − G1y (x2 , y2 ) ⎥ = ⎢ 1⎥ ⎡ T ⎤. ⎥ ⎢ ⋮⎥ ⎣ y ⎦ ⎢ ⋯ ⎥ ⎢ ⎥ ⎢ ⎣ G2y (xn , yn ) − G1y (xn , yn )⎦ ⎣ 1⎦
(3b)
Accordingly, Tx and Ty are the averages of the differences in the left-hand side of Eqs. (3a) and (3b), Next, the adjusted results of G2x and G2y can be solved by using Eq. (2). Subsequently, the adjusted results of G2x and G2y on the W2 can be connected with the results of G1x and G1y on the W1 successfully; as a result, we obtain the stitched slope on the two adjacent areas. Next, consider another case. Both tilt and defocus are introduced in the measured slope data, and the relationship between the phase distributions Φ1 and Φ 2 in their common area fulfills the following equation:
Φ2 − Φ1 = P + Tx x + Ty y + D (x2 + y2 ),
(4)
where D is the coefficient of defocus. Likewise, considering the derivation of x -tilt and y -tilt in Eq. (4), we have
G2x − G1x = Tx + 2Dx x G2y − G1y = Ty + 2Dy y .
(5)
Similarly, suppose there are n sampling points in the common area. From Eq. (5), we obtain
⎡ G2x (x1,y ) − G1x (x1, y ) ⎤ ⎡ 1, x ⎤ 1 1 1 ⎥ ⎢ ⎢ ⎥ ⎢ G2x (x2,y2 ) − G1x (x2 , y2 ) ⎥ = ⎢ 1, x2 ⎥ ⎡⎢ Tx ⎤⎥, ⎥ ⎢ ⋮ ⎥ ⎣ 2D x ⎦ ⎢ ⋯ ⎥ ⎢ ⎢ ⎥ ⎣ G2x (xn , yn ) − G1x (xn , yn )⎦ ⎣ 1, xn ⎦
(6a)
⎡ G2y (x1,y ) − G1y (x1, y ) ⎤ ⎡ 1, y1 ⎤ 1 1 ⎥⎡ ⎥ ⎢ ⎢ ⎤ ⎢ G2y (x2,y2 ) − G1y (x2 , y2 ) ⎥ = ⎢ 1, y2 ⎥ ⎢ Ty ⎥. ⎥ ⎢ ⋮ ⎥ ⎢⎣ 2Dy ⎥⎦ ⎢ ⋯ ⎥ ⎥ ⎢ ⎢ ⎣ G2y (xn , yn ) − G1y (xn , yn )⎦ ⎣ 1, yn ⎦
(6b)
According to the least-squares method, we obtain the coefficients of Tx , Ty and Dx , Dy in Eq. (6a) and (6b). Next, the adjusted results of G2x and G2y can be obtained by using Eq. (5). As a result, the stitched slope on the two adjacent areas is obtained. 3. Error analysis Fig. 1. Schematic of subaperture stitching.
For an actual measurement, the slope data measured by a
126
G. Guo et al. / Optics Communications 380 (2016) 124–133
Hartmann–Shack sensor could be polluted by random noise to a certain extent. The accuracy of these coefficients related to the stitching parameters obtained by using Eqs. (3) and (6) are influenced by noise. In addition, the accuracy of the coefficients is also affected by the number of sample points in their common area. If these parameters have some errors, then the decrease of the accuracy of the adjusted slope is inevitable. In this paper, the influence of the two factors on the accuracy of stitched slope data is discussed mathematically in detail. The error analysis is conducted in the two respective cases of with or without defocus by using (Eqs. (3) and 6).
Therefore, each point in G″2x (xi , yi ) will have Gaussian noise with the variance described in Eq. (10). This result shows that the variance of T1x on the overlapping area between A1 and A2 is not equal to that of the adjusted slope on the other area of subaperture A2. In addition, the uncertainty on other area of subaperture A2 is increased by the variance of the T1x on the shaded area. Using the approach of data stitching between the subaperture A1 and A2 described above, the tilt coefficient T2x on the subaperture A3 can be written as
3.1. Slope data including the tilt coefficient
T2x =
We consider the slope data in the x direction as an example and suppose that the subapertures A1, A2 , ⋯ , AN in the measured surface are square shaped, as shown in Fig. 2; A1 is defined as the standard reference subaperture and the arrow indicates the path of stitching. The size of the slope matrix in each subaperture is W × W , and the shaded area with the size of the slope matrix W × M is the common area of the two adjacent subapertures. Let the two slope matrices in the x direction on the two adjacent subapertures A1 and A2 be G1x and G2x , respectively. We assume the RMS (Root Mean Square) value of the random noise in the measured slope of each subaperture is σ . In fact, in a stable test environment, the random noise can be statistically estimated through multiple measurements. For the stitching between the subaperture A1 and A2 , the unknown parameter T1x (the tilt coefficient in the x direction on A2) can be easily written according to the Eq. (3a) as follows:
Because the matrices G3x (xi , yi ), G′2x (xi , yi ) and T1x are mutually independent, the variance of T2x is
T1x =
1 WM
WM
∑ ⎡⎣ G2x (xi , yi ) − G1x (xi , yi ) ⎤⎦. i=1
(7)
Note that the slope matrix G1x (xi , yi ) and G2x (xi , yi ) are mutually independent, so the variance of T1x on the shaded area can be expressed as follows:
σT21x =
2σ 2 . WM
(8)
In this paper, to distinguish the variance of stitching parameters on the overlapping area (the shaded) between A1 and A2 from that of adjusted slope on the other area of subaperture A2, in which the shaded area is not included, we define two new slope matrices on the remaining area of subaperture A2 as G′2x (xi , yi ) and G″2x (xi , yi ) , which correspond to the portion of matrix G2x (xi , yi ) and the adjusted results, respectively. After removing the tilt coefficient T1x from the G′2x (xi , yi ), we obtain the slope matrix G″2x (xi , yi ), which can be written as
G″2x (xi , yi ) = G′2x (xi , yi ) − T1x
(9)
Note that the matrices G′2x (xi , yi ) and T1x are mutually independent as well; thus, the variance of the adjusted slope matrix G″2x (xi , yi ) can be obtained as follows:
σ 2G 2″x (xi, yi ) = σ 2 +
2 2 σ WM
(10)
WM ∑i = 1 ⎡⎣ G3x (xi , yi ) − G′2x (xi , yi ) ⎤⎦
σT22x =
WM
+ T1x
4σ 2 WM
(11)
(12)
Likewise, after removing the tilt coefficient T2x from G′3x (xi , yi ), we obtain the adjusted slope matrix G″3x (xi , yi ) as follows:
G″3x (xi , yi ) = G′3x (xi , yi ) − T2x
(13)
Note that the two matrices G′3x and T2x are mutually independent; thus, the variance of G″3x (xi , yi ) can be obtained as follows,
σ 2G3″x (xi, yi ) = σ 2 +
4σ 2 WM
(14)
Accordingly, based on the two Eqs. (10) and (14), we can see that the uncertainty of tilt coefficients due to the random noise will transfer from one to the next subaperture and is increased gradually for the subsequent one to the last one. As a result, the stitching accuracy decreases. Therefore, for the stitching of N subapertures, which are stitched sequentially from one to the next subaperture, we can conclude the general equations for the calculation of the variance of T(N − 1) x and the G″Nx (xi , yi ) separately are as follows,
σT2(N − 1) x = σ G2″
Nx (xi, yi )
2 (N − 1) σ 2 , WM
= σ2 +
(15)
2 (N − 1) σ 2 . WM
(16)
Here, Eqs. (15) and (16) show that the number of subapertures and the overlapping area should be properly adopted according to the requirement of accuracy for the slope measurement when we conduct the stitching by use of slope data provided by Hartmann– Shack sensor. Similarly, for the discussion regarding the measured slope data in the y direction, we can also obtain similar equations as Eqs. (15) and (16). Eq. (16) is the basic error transfer formula proposed in this paper for the slope stitching. By using Eq. (16), we can calculate the stitching error on each subaperture if the random noise is measured in advance and the overlapping area is given. 3.2. Slope data including tilt and defocus coefficients Sometimes, the defocus may disturb the measurement besides the tilt. For this situation, from the statistical theory, the Eq. (6a) is suitable for linear regression model [17], and the slope data in the x direction can be written as
li = 2 × βxi + β0 + εi , Fig. 2. Schematic of multiple subapertures stitching.
i = 1, 2, ⋯, n.
(17)
Here, li is defined as li = G2x (xi , yi ) − G1x (xi , yi ), and n is the
G. Guo et al. / Optics Communications 380 (2016) 124–133
number of sampling points in the common area. β0 and β are the regression constant and the regression coefficient, respectively. εi represents the influence of all the other random variables on the slope data li . Suppose Dx and Tx are the least squares estimation of the parameter β and β0 , respectively; in this case, the linear regression equation can be written as follows:
^ li = 2 × Dx xi + Tx
(18)
from the G′2x (xi , yi ), we can write G″2x (xi , yi ) on the other area of subaperture A2 as follows:
^ G″2x (xi , yi ) = G′2x (xi , yi ) − li .
σ G2″
2x (xi, yi )
^ L − Xb = V .
(19)
(
σT2x
(21)
^ According to statistical theory, the variance of li in Eq. (18) can be expressed as
σ ^2 = 4σ D2x x2 + σT2x + 4xCov (Dx , Tx ),
(22)
li
where
Dx
σ2
and
are the variances of the Dx and Tx , re-
Tx
)
6σ 2 , t 2W (M 3 − M ) ⎞ 2⎛ 1 12x¯ 2 = + 2 3 ⎜ ⎟ σ 2, W⎝M t (M − M ) ⎠
(20)
)−1X T L.
σ2
(
⎤ ⎥ 2 σ + σ 2. 2⎥ ⎦
σ D2x =
x1 ⎤ ⎡ v1⎤ ⎥ ⎢v ⎥ x2 ⎥ ^ ⎡ Tx ⎤ 2 ,b=⎢ ⎥, V = ⎢ ⎥, ⋮⎥ ⎣ 2D x ⎦ ⎢ ⋮⎥ ⎥ ⎣ vn ⎦ xn ⎦
and vi , which is the difference between the test value li and the ^ regression value li , represents the residual. Suppose the RMS value of the random noise in the measured slope distribution is σ , i.e., each point in G1x and G2x has measurement noise with the RMS value of σ , and the RMS value of li with Gaussian distribution is ^ 2 σ .The unknown b can be obtained by the following equation:
^ b = XTX
⎡ (xi − x¯ )2 ⎢1 = 2⎢ + 1 n n n ∑i = 1 xi2 − n ∑i = 1 xi ⎣
(30)
Further, assume that t is the physical interval between two sampling points in the x direction; thus, we can separately rewrite Eqs. (25), (26) and (30) as follows:
where
⎡ l1 ⎤ ⎡1 ⎢ ⎥ ⎢ 1 l2 ⎥ ⎢ L= , X=⎢ ⎢ ⋮⎥ ⎢⋮ ⎢ ⎥ ⎢ ⎣1 ⎣ ln ⎦
(29)
^ Note that G′2x (xi , yi ) and li are mutually independent, the variance of G″2x (xi , yi ) can be obtained as
The matrix form of the corresponding error equations is given by
127
σG2 2x (xi, yi ) ″
2⎡ 1 12 (x − x¯ )2 ⎤ 2 = + 2 3i ⎢ ⎥ σ + σ 2. W⎣M t (M − M ) ⎦
(31) (32) (33)
Eq. (33) demonstrates that the error of slope stitching is related to the random noise, the number of sampling points in common area, and the location of the sampling points. Likewise, the uncertainties of the tilt and defocus coefficients will be both transferred and are increased gradually from one to the next subaperture. As a result, the stitching accuracy will suffer more than for the case of only the tilt/tip existing. In addition, it is also more complicated to calculate the stitching accuracy on each subaperture. Similarly, the discussion of the slope stitching accuracy in the y direction will have the same result as Eq. (33). When we use the Hartmann–Shack wavefront sensor to measure the slope data on subapertures, generally, only the tilt coefficient is introduced into the measurement. Therefore, in the following analysis, the defocus is neglected.
spectively. Cov (Dx , Tx ) is the covariance between Dx and Tx . Solving Eq. (21), we obtain
(
Tx =
n
n
) (
n ∑i = 1 xi li − ∑i = 1 xi
Dx =
n
n
n ∑i = 1 li
n ∑i = 1 xi
2 [n ∑i = 1 xi2 − ( ∑i = 1 xi
(
n ∑i = 1 xi2
)(
n l i=1 i 2 ) ]
)( ∑
)−( − (∑
)( )
n ∑i = 1 xi li
2 n x i=1 i
n n ∑i = 1 xi2
4. Stitching mode
) (23)
) (24)
Therefore, the respective variances of Dx and Tx and the covariances between Dx and Tx can be obtained as follows,
σ D2x =
σT2x
1 n 2 ∑ n x i2 − ∑ n x i i=1 i=1
(
2
)
σ2 (25)
⎛ x¯ 2 ⎜1 = 2⎜ + n 1 n 2 n ∑i = 1 xi − n ∑i = 1 xi ⎝
(
Cov (Dx , Tx ) =
⎞ ⎟ 2 σ 2⎟ ⎠
)
n − ∑i = 1 xi n
(
2
n
)
n ∑i = 1 xi2 − ∑i = 1 xi
(26)
σ2 (27)
Substituting Eqs. (25), (26) and(27) into Eq. (22), the variance of
^ li is expressed as
⎡ (x − x¯ )2 ⎢1 σ ^2 = 2 ⎢ + 1 n n li n ∑i = 1 xi2 − n ∑i = 1 xi ⎣
(
n
⎤ ⎥ 2 σ . 2⎥ ⎦
)
(28)
Where x¯ = ∑i = 1 xi . Based on the above-described treatment of the stitching with tilt coefficient, after removing the tilt and defocus coefficients
Generally, at present, the stitching mode is serial stitching, as shown in Fig. 3(a), which has the advantage of a simple mechanical structure and easier operation. However, serial stitching requires more time to finish the stitching task, and, as we know, the error will accumulate gradually and cause a large stitching error to occur in the measurement. In this paper, the parallel stitching mode (as shown in Fig. 3(b)) is proposed and discussed. Unlike the serial stitching, it can stitch multi-subapertures simultaneously, and the parallel stitching method requires less time to finish a stitching task with reduced loss of accuracy. Thus, it is necessary to apply an appropriate stitching mode to achieve our goal with a smaller error, according to the path and sequence of stitching. By considering an example, the error of the serial stitching mode and parallel stitching mode will be analyzed in detail. For simplicity, a component model with 3 × 3 subapertures is established, as shown in Fig. 3. The size of the slope matrix in each subaperture is 10 × 10. In addition, the size of the slope matrix in the overlapping area is 2 × 10 or 10 × 2. Each element (point) of slope matrix in each subaperture has random noise of mean equal to 0 and RMS value of σ . 4.1. Serial stitching As shown in Fig. 3(a), suppose that the subaperture A1 is the reference aperture. The stitching treatment between the adjacent A1, A2 and A3 is the same as shown in Fig. 2 from the left to the
128
G. Guo et al. / Optics Communications 380 (2016) 124–133
Fig. 3. The two modes of 3 × 3 subapertures stitching. (a) Serial stitching (b) Parallel stitching.
right. Moreover, the stitching error can be calculated by using Eq. (16). Thus, we only need to discuss the stitching between the subaperture A4 and A5. For simplicity of calculation, we ignore the common area (the shaded area in Fig. 3(a)) between subapertures A2, A3, and A4 when the stitching between A4 and A3 is implemented; thus, the overlapping area becomes the matrix size of 2 × 8. According to Eq. (11), the tilt coefficient on A4 can be expressed as
T3x =
1 16
16
∑ ( G4x (xi , yi ) − G′3x (xi , yi ) ) + T2x. i=1
(34)
As a result, the variance of T3x can be determined to be equal to 0.325σ 2 by use of Eq. (15). Further, after canceling the tilt coefficient T3x from the G′4x (xi , yi ), we obtain the adjusted G″4x (xi , yi ), namely, G″4x (xi , yi ) = G′4x (xi , yi ) − T3x . Likewise, the variance of G″4x (xi , yi ) can be determined to be equal to 1.325σ 2 by use of Eq. (16). In the same way, when the stitching of slope data on A5, A2 and A4 is implemented, we ignore the data in the common area ( 2 × 2 sample points, i.e., the shaded area in Fig. 3(a)) between A2 , A4 and A5, i.e., the tilt coefficient on A5 can be expressed as T4x =
16 ∑16 i = 1 (G5x (xi, yi ) − G′4x (xi, yi ) + T3x ) + ∑k = 1 (G5x (xk , yk ) − G′2x (xk , yk ) + T1x ) 32
(35)
where subscript i denotes the overlap between A5 and A4 , and subscript k denotes the overlap between A5 and A2. For approximate calculation, we suppose the variance of each point on the entire subaperture A2 is σ 2G 2x (xi, yi ) ; in addition, because T1x and T3x ″ are not statistically independent, the variance of T4x is given by
σT24x =
σ2 1 1 1 + σT23x + σT21x + Cov (T1x, T3x ). 16 4 4 2
(36)
where Cov (T1x, T3x ), which is the covariance between the T1x and T3x , is expressed as Cov (T1x, T3x ) = ⎛ ⎞ ∑16 (G4x (xi, yi ) − G′3x (xi, yi )) + ∑16 k = 1 [(G3x (xk , yk ) − G′2x (xk , yk )) + T1x ] ⎟ Cov ⎜ T1x, i = 1 ⎜ ⎟ 16 ⎝ ⎠
(37)
Because G4x (xi , yi ) , G′3x (xi , yi ) , G′2x (xk , yk ) and T1x are statistically independent, the Cov (T1x, T3x ) is equal to σT21x . Substituting the quantity into Eq. (36), we obtain the variance of T4x , which is equal to 0.219 σ 2. Likewise, after removing the tilt coefficient T4x from G′5x (xi , yi ), we obtain the adjusted G″5x (xi , yi ), namely, G″5x (xi , yi ) =
G′5x (xi , yi ) − T4x ; thus, the variance of G″5x (xi , yi ) can be solved and is found to be equal to 1.219σ 2. As mentioned above, the stitching for the data on A7 is found to be similar to that on A4 , and the stitching for the data on A6 , A8 and A9 is similar to that on A5. Thus, we can approximately determine the stitching error of each point of the stitched slope on the test component if the variance of random noise has been estimated and the overlapping area is obtained in advance. Without exception, the method and formulas mentioned above are also applicable to the serial stitching with more subapertures. 4.2. Parallel stitching For parallel stitching, as shown in Fig. 3(b), the center subaperture A1 is considered as the reference aperture. If these subapertures have the same distance from the reference aperture, then their adjusted slope will have an identical stitching error. Thus, we only must discuss the stitching for the data between the four subapertures on the top right ( A1, A2a , A2b , A3a ). The tilt coefficient on A2a and A2b are defined as T1ax and T1bx , respectively, according to Eq. (7): T1a x = ∑20 i = 1 (G2a x (xi , yi ) − G1x (xi , yi )) /20, 20
= ∑i = 1 (G2b x (xi , yi ) − G1 x (xi , yi )) /20 . Due to their dependence, the covariance between T1a x and T1b x is expressed as Cov (T1a x, T1b x ). Because G2a x (xi , yi ), G2b x (xi , yi ) and G1 x (xi , yi ) are statistically independent on the other area of themselves and there are 2 × 2 points in the common area (the shaded area in Fig. 3(b)), accordingly, Cov (T1a x, T1b x ) can be solved as follows: T1b
x
⎡ ∑ 4 G (x , y ) ∑ 4 G (x , y ) ⎤ 1x i i 1x i i ⎥ σ2 = Cov (T1a x, T1b x ) = Cov ⎢ i = 1 , i =1 . ⎢ ⎥ 100 20 20 ⎣ ⎦
(38)
For simplicity of calculation, when the slope data on A3a is separately stitched with that on A2a and A2b , we already ignore the data on the common area ( 2 × 2 sample points, i.e., the shaded area in Fig. 3(b)) between A2a , A2b , and A3a . Thus, the tilt coefficient on A3a is given by T2a x =
16 16 ∑i = 1 (G3a x (x i , yi ) − G′2a x (x i , yi )) + ∑i = 1 (G3b x (x i , yi ) − G′2b x (x i , yi ))
32 +
T1a x T + 1b x 2 2
Because only T1a x and T1b parameter T2a x is given by
(39) x
are dependent, the variance of
G. Guo et al. / Optics Communications 380 (2016) 124–133
Fig. 4. The statistically simulated slope error of serial stitching.
σT22a x =
σT2 σT2 1 σ2 + 1a x + 1b x + Cov (T1a x, T1b x ) 16 4 4 2
(40)
Substituting Eq. (38) into Eq. (40), we obtain the variance of the parameter T2a x to be equal to 0.118 σ 2. Likewise, G′3a x (xi , yi ) is the slope data on the other area of A3a . After canceling the tilt coefficient T2a x from G′3a x (xi , yi ). we obtain Considering G″3a x (xi , yi ) = G′3a x (xi , yi ) − T2a x . G″3a x (xi , yi ): G′3a x (xi , yi ) and T2a x are statistically independent, the variance of G″3a x (xi , yi ) is solved to be equal to 1.118σ 2. Similarly, the stitching error of each point in the stitched slope can be calculated or predicted by this method if the variance of random noise has been estimated and the overlapping area is obtained. In addition, as in the case of serial stitching, the parallel stitching formulas above are also applicable to the case with more subapertures. 4.3. Simulation To verify the stitching error formulas for the two stitching modes, a simulation experiment of subapertures stitching is implemented. As mentioned above, a component with 3 × 3 subapertures is discussed, the size of the slope matrix is 10 × 10 for each subaperture, and the size of the slope matrix in the overlapping area is 2 × 10 or 10 × 2. Next, the error of each point of the stitched slope on a component is calculated using the error formulas proposed in this paper. The serial stitching simulation experiment procedure is as follows: 1. Generate a matrix of slope data with the order of magnitude of 1 × 10−6rad for a test optical component first. 2. A Gaussian random noise matrix with mean of 0 and RMS value of σ = 2.423 × 10−7rad is added to the above slope data of each subaperture. 3. An amount of tilt, which is used to simulate the misalignment of a slope sensor, is added to the slope data of each subaperture (except the subaperture A1). 4. Stitch the slope data on the all subapertures in serial stitching mode. Repeating the steps 2–4, the stitching for the 9 regions that constitute the entire stitched slope is implemented by using serial stitching mode. Next, the RMS of noise at each slope data in each region is statistically calculated after repeating the trials 100,000 times, as shown in Fig. 4, in which the color of each point reflects the value of noise, but they are nearly the same in one region. In addition, the average of RMS of noise at all points in each region is
129
compared with the theoretical value from our formulas (Eqs. (34)– (37) and (16)) proposed in the paper, as shown in Table 1. It can be observed that the statistical value on the different regions is almost completely equal to the predicted value on the corresponding regions by our formulas. In addition, the error of each point in the fourth region but not in the ninth region has the largest stitching error among the nine regions. Therefore, the result of this paper provides us with formulas to exactly determine the stitching error on each subaperture if the variance of random noise has been estimated and the overlapping area is given. Instead of serial stitching, the simulation procedure mentioned above is implemented for parallel stitching. Likewise, after repeating steps 2–4, the RMS of noise at each slope data in each region is statistically calculated after repeating the trials 100,000 times, as shown in Fig. 5. In addition, the average of RMS of noise at all points in each region is compared with the theoretical value from our formulas (Eqs. (38)–,(40) and (16)) proposed in the paper, as shown in Table 2. In addition, the error of each point in the third region has the largest error value among the nine regions. As a result, the formulas provided in this paper can be used to estimate the error in each subaperture if the variance of random noise has been estimated and the overlapping area is given in advance. Moreover, we can see that the error obtained using the parallel stitching mode is smaller than that obtained using the serial stitching mode. For further comparison, using the weighted average error of each point in the stitched slope as the reference to evaluate the stitching accuracy for the two modes, we define the error as N
error =
∑ σ Gi″ × pi .
(41)
i=1
where N is the number of the sub-regions in the component, σ Gi″ is the theoretical variance value of the error on the ith region, and pi is the ratio of the area of the ith region to that of the entire component. So thus, the error of the serial stitching and parallel stitching are 2.63 × 10−7rad and 2.53 × 10−7rad from Eq. (41), respectively. This result shows that the parallel stitching mode has a higher accuracy than the serial stitching mode. For an element with 3 × 3 subapertures and two rows or two columns overlapping area, with the increase of the random noise from 0 to 5 × 10−7rad , the weighted average errors of the two stitching modes at different noise are obtained using Eq. (41) and are illustrated in Fig. 6. We can see that the errors of both the serial stitching and parallel stitching linearly increase with increasing random noise, and the error of the serial stitching is greater than that of the parallel stitching. Similarly, given the random noise of 2.423 × 10−7rad , when the overlapping area between adjacent subapertures increases from 2 to 9 rows or columns, the errors of the two stitching modes with the different overlapping area are determined by Eq. (41) and are Table 1 Comparison of the mean of statistical errors and theoretical values in 9 regions by serial stitching. Region
The statistical value (100,000 times)
The theoretical value
1
2.42 × 10−7
2.42 × 10−7
2
2.54 × 10−7
2.54 × 10−7
3
2.66 × 10−7
2.66 × 10−7
4
2.79 × 10−7
2.79 × 10−7
5
2.66 × 10−7
2.68 × 10−7
6
2.56 × 10−7
2.56 × 10−7
7
2.70 × 10−7
2.70 × 10−7
8
2.69 × 10−7
2.69 × 10−7
9
2.75 × 10−7
2.74 × 10−7
130
G. Guo et al. / Optics Communications 380 (2016) 124–133
Fig. 5. The simulated slope error by using the parallel stitching mode.
Table 2 Comparison of the mean of statistical errors and theoretical values in 9 regions by parallel stitching. Region
The statistical value (100,000 times)
The theoretical value
1
2.42 × 10−7
2.42 × 10−7
2
2.54 × 10−7
2.54 × 10−7
3
2.56 × 10−7
2.56 × 10−7
Fig. 6. Comparison of the errors for two stitching modes when the random noise is increasing.
shown in Fig. 7. The errors of both the serial stitching and parallel stitching decrease with the increase of the overlapping area, and the error of serial stitching is greater than that of parallel stitching.
5. Wavefront reconstruction of A component In the following simulation, assume that there is a component with the sizes of 420 mm 420 mm and the slope matrices of the whole component in x - and y -directions have dimensions of 42 × 42. Based on the slope matrices without noise, the wavefront of this component is reconstructed by a zonal reconstruction method using the Southwell's slope geometry [18] and is shown in Fig. 8. For testing our method, the component is divided into 5 × 5 subapertures. The size of the slope matrix of each subaperture is 10 × 10, and the size of the overlapping area is 2 × 10 or 10 × 2. Each point of the slope matrix of each subaperture has random noise whose mean is 0 and RMS value s is 2.423 × 10−7rad . Next, the slope is stitched by the serial stitching and the parallel stitching, respectively. Subsequently, the wavefront is reconstructed from the stitched slope data. The simulation results are shown in Figs. 9 and 10 as follows, Wavefront residual refers to the reconstruction from the stitched slope minus that from the original slope. From Figs. 9 and 10, we can see the reconstructed wavefront from the stitched slope by both the stitching modes are nearly similar to that from the original slope. In addition, the residual data show the accuracy of parallel stitching is better than that of serial stitching. To further verify the conclusion, 1000 simulation trials for wavefront reconstruction were implemented by using the two stitching modes. In each trial, the PV and RMS of the wavefront residual from the parallel stitching mode are subtracted from that of the serial stitching mode, and the distribution of their differences are shown in Fig. 11(a) and (b), respectively. The number of “*” points above the red straight line is observed to be much higher than that under this line in the two figures. From the statistical results, we can conclude that the error of the serial stitching is higher than that of the parallel stitching.
Fig. 7. The comparison of the errors for the two stitching modes when the overlapping area is increasing.
6. Conclusions In this paper, the mathematical model of the slope stitching based on the least squares fitting routine is presented. Using the model the stitching parameters, such as tilt and defocus, can be determined. After the parameters are removed, the surface figure of
Fig. 8. Wavefront construction from the original slope data.
G. Guo et al. / Optics Communications 380 (2016) 124–133
131
Fig. 9. The simulated results from the parallel stitching mode: (a) wavefront estimation from stitching slope and (b) wavefront residual.
Fig. 10. The simulated results from the serial stitching mode: (a) wavefront estimation from stitching slope and (b) wavefront residual.
Fig. 11. Statistical comparison of the two modes: (a) the PV difference of wavefront residual and (b) the RMS difference of wavefront residual.
the test component can be retrieved using a zonal algorithm or a model approach. More importantly, the error transfer formulas are determined when the variance of random noise has been estimated
and the overlapping area is given in advance. Our discussion indicated that the uncertainty of tilt or defocus due to the random noise will transfer from one subaperture to the next subaperture and
132
G. Guo et al. / Optics Communications 380 (2016) 124–133
is increased gradually; such uncertainty will cause the decrease of stitching accuracy. Mathematically, the formulas for stitching accuracy evaluation can be used to calculate the stitching error of each point on each subaperture for the parallel mode and the serial mode. The results of the simulation experiments demonstrated that all the stitching accuracy formulas proposed are verified and that these formulas have the ability to evaluate the stitching accuracy. In practice, the number of subapertures and the overlapping areas should be properly adjusted according to the requirement of accuracy for the slope measurement when we conduct the stitching by use of slope data provided by a Hartmann–Shack sensor. Moreover, the accuracy of these two stitching modes is compared, and the simulation experiment and theoretical analysis show that the error of parallel stitching is less than that of serial stitching.
n
n x2 i=1 i
(∑
Tx =
n
) (
n l i=1 i 2 ) ]
)( ∑
n
2 [n ∑i = 1 xi2 − ( ∑i = 1 xi n l i=1 i n ∑ n i = 1 xi2
)( ∑
) (A.4)
n x i=1 i n ∑ ( i = 1 xi )2
) − (∑ −
n xl i=1 i i
)( ∑
) (A.5)
From the simplification of the formula of Eq. (A.4), we obtain n ⎡ n ⎤ 1 ∑i = 1 ⎣ (nxi − ∑i = 1 xi ) li ⎦ n n 2 2 2 n ∑i = 1 xi − ( ∑i = 1 xi )
Dx =
(A.6)
The process of solving the variance of Dx is as follows:
σ D2x =
1 4
Funding information
= This work was supported by the National Defense Basic Research Program of China (JSJC2013212C002) and the National High Technology Research and Development Program of China (2015AA015902).
=
n ⎡ n ∑i = 1 ⎢⎣ nxi − ∑i = 1 xi
⎤ × 2σ 2⎥⎦ 2 ⎤2 ⎡ n n 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
2
)
(
)
2 2⎤ ⎡ 2 n 2 n n n ∑i = 1 xi − 2n ∑i = 1 xi + n ∑i = 1 xi ⎥⎦ × 2σ 2 1 ⎢⎣ 2 ⎤2 ⎡ 4 n n 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
n 1 2 n ∑n x 2 − ∑n xi i=1 i=1 i
(
) (
2
)
(
)
)
σ2 (A.7)
Simplifying the formula of Eq. (A.5), we obtain
Appendix A From Eq. (20), we obtain
(
n
(
n ∑i = 1 xi li − ∑i = 1 xi
Dx =
n n n ∑i = 1 ⎡⎣ ( ∑i = 1 xi2 − xi ∑i = 1 xi ) li ⎤⎦
Tx =
n ⎤ ⎡ n ⎢ ∑ xi2 − ∑ xi ⎥ ⎢ 1 −1 i=1 i=1 ⎥ XTX = ⎥ ⎢ n n n 2 2 n ∑i = 1 xi − ( ∑i = 1 xi ) ⎢ ⎥ x n ⎥ ⎢− ∑ i ⎦ ⎣ i=1
n
n
n ∑i = 1 xi2 − ( ∑i = 1 xi )2
(A.8)
The process of solving for the variance of
)
σT2x
(A.1)
⎡ n ⎤ ⎢ ∑ li ⎥ ⎢ i=1 ⎥ XTL = ⎢ n ⎥ ⎥ ⎢ ⎢ ∑ xi ⎥ ⎣ i=1 ⎦
=
(
)
2 ⎤2 ⎡ n n 2 ⎣⎢ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
⎡ n 2σ 2 × ⎣⎢ n ∑i = 1 xi2
(
=
)
2
)
−
n ∑i = 1 xi2 2 ⎤2 ⎡ n n 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
⎡⎛ n ⎞⎛ n ⎞ ⎛ n ⎞⎛ n ⎞⎤ ⎢ ⎜ ∑ x 2 ⎟ ⎜ ∑ li ⎟ − ⎜ ∑ x i ⎟ ⎜ ∑ x i li ⎟ ⎥ i ⎟⎜ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎥ ⎢⎝ ⎠⎝ i=1 ⎠ ⎝ i=1 ⎠⎝ i=1 ⎠ 1 i=1 ^ ⎥ ⎢ b= n n 2 2 n n n ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎥ n ∑i = 1 xi − ( ∑i = 1 xi ) ⎢ ⎥ ⎢ n ⎜⎜ ∑ xi li ⎟⎟ − ⎜⎜ ∑ xi ⎟⎟ ⎜⎜ ∑ li ⎟⎟ ⎥⎦ (A.3) ⎢⎣ ⎝ i=1 ⎠ ⎝ i=1 ⎠⎝ i=1 ⎠
n 2 ∑i = 1 xi2
⎡ n 2 ⎣⎢ n ∑i = 1 xi
(A.2)
Solving Eq. (21), we obtain
is as follows:
n n n ∑i = 1 ⎡⎣ ∑i = 1 xi2 − xi ∑i = 1 xi2 × 2σ 2⎤⎦
(
=
Tx
n x2 i=1 i
(∑ − (∑
n
+ ∑i = 1 xi2
n x2 i=1 i
)( ∑
) ⎤⎦⎥
2 ⎤2 n x ⎥ i=1 i ⎦
)
× 2σ 2
)
⎛ x¯ 2 ⎜1 = 2⎜ + 1 n n 2 n ∑i = 1 xi − n ∑i = 1 xi ⎝
(
⎞ ⎟ 2 σ 2⎟ ⎠
)
(A.9)
The process of solving for the covariance of Dx and Tx is as follows:
Dx and Tx can be obtained as
n ⎡ n n ⎛ ⎤⎞ n n 2 ⎜ 1 ∑i = 1 nxi − ∑i = 1 xi li ∑i = 1 ⎣ ∑i = 1 xi − xi ∑i = 1 xi li ⎦ ⎟ cov (Dx , Tx ) = cov ⎜ , ⎟⎟ 2 n n ⎜ 2 n ∑n x 2 − ∑n x 2 n ∑i = 1 xi2 − ∑i = 1 xi ⎝ ⎠ i=1 i i=1 i n ⎡ n n n n ⎡ n n n 2 2⎤ 2 2⎤ 1 ∑i = 1 ⎣ nxi − ∑i = 1 xi ∑i = 1 xi − xi ∑i = 1 xi σ li ⎦ 1 ∑i = 1 ⎣ nxi − ∑i = 1 xi ∑i = 1 xi − xi ∑i = 1 xi σ li ⎦ = = 2 ⎤2 2 ⎤2 ⎡ ⎡ 2 2 n n n n 2 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦ ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
(
(
(
) )
(
)(
1 2
)
)
(
=
)
(
)
)
(
)
(
⎡ n n n n n 2σ 2 × ⎢⎣ n ∑i = 1 xi )( ∑i = 1 xi2 − n ∑i = 1 xi )( ∑i = 1 xi2 − n ∑i = 1 xi 2 ⎤2 ⎡ n n 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
)(
) (
(
)
n x2 i=1 i
)( ∑
n x i=1 i
) + (∑
)
2⎤ n x ⎦ i=1 i ⎥
)( ∑
)
n
=
− ∑i = 1 xi n n ∑i = 1 xi2
(
n
− ∑i = 1 xi )
2
)
σ2 (A.10)
G. Guo et al. / Optics Communications 380 (2016) 124–133
References [1] B.E. Catanzaro, J.A. Thomas, E.J. Cohen, Comparison of full-aperture interferometry to subaperture stitched interferometry for a large-diameter fast mirror, in: Presented at the International Symposium on Optical Science and Technology, 2001. [2] Chunyu Zhao, J.H. Burge. Stitching of off-axis sub-aperture null measurements of an aspheric surface, in: Proceedings of SPIE – The International Society for Optical Engineering, 2008, pp. 706316–706316-7. [3] Michael Bray, Stitching interferometer for large optics: recent developments of a system, in: Presented at the Third International Conference on Solid State Lasers for Application to Inertial Confinement Fusion, 1999. [4] Yingjie Yu, Mingyi Chen, Correlative stitching interferometer and its key techniques, in: Proceedings of International Symposium on Optical Science and Technology. International Society for Optics and Photonics, 2002. [5] B. Koehler, M. Kraus, J. Moresmau, K. Wirenstrand, M. Duchateau, P. Duhoux, R. Karban, C. Flebus, E. Gabriel, and O. Pirnay, The VLTI Auxiliary Telescopes: commissioning of AT1, in: Proceedings of SPIE Astronomical Telescopes þ Instrumentation International Society for Optics and Photonics, 2004, p. 600. [6] Franziska Harms, Jürgen Wolf, Patrick Waddell, Edward Dunham, Andreas Reinacher, Ulrich Lampater, Holger Jakob, Lisa Bjarke, Sybil Adams, Randy Grashuis, Allan Meyer, Kenneth Bower, Keith Schweikhard, Thomas Keilig, On sky testing of the SOFIA telescope in preparation for the first science observations, in: Proceedings of SPIE Optical Engineering þ Applications. International Society for Optics and Photonics, 2009, pp. 745303–745303-15. [7] C.J. Kin, J.C. Wyant, Subaperture test of a large flat or a fast aspheric surface, J. Opt. Soc. Am. 71 (1981) 1587.
133
[8] M. Otsubo, K. Okada, J. Tsujiuchi, Measurement of large plane surface shapes by connecting small-aperture interferograms, Opt. Eng. 33 (2) (1994) 608. [9] P. Murphy, J. Fleig, G. Forbes, D. Miladinovic, G. Devries, S. O'Donohue, Subaperture stitching interferometry for testing mild aspheres, in: Presented at the SPIE Optics þ Photonics, 2006. [10] ChaoWen Liang, HungSheng Chang, PoChih Lin, ChengChung Lee, YiChun Chen, Vibration modulated subaperture stitching interferometry, Opt. Express 21 (15) (2013) 18255. [11] Rongzhu Zhang, C. Yang, Q. Xu, Precision analyses of a stitching interferometer testing system, Appl. Opt. 45 (11) (2006) 2399–2403. [12] Guangda Zhan, Xudong Xu, Zhengxiang Shen, Xiaoqiang Wang, ZhanShan Wang, Huasong Liu, Yiqin Ji, A global subaperture stitching algorithm and its precision evaluation, 2012-01-01, 2012. [13] D.R. Neal, J. Copland, D.A. Neal, Shack–Hartmann wavefront sensor precision and accuracy, 2002-01-01, 2002. [14] D.R. Neal, Apparatus and method for evaluating a target larger than a measuring aperture of a sensor. [15] T.D. Raymond, D.R. Neal, D.M. Topa, T.L. Schmitz, High-speed noninterferometric nanotopographic characterization of Si wafer surfaces, Nanoscale Opt. Appl. 4809 (2002) 208–216. [16] A.S.Y. He, G. Tang, Subaperture stitching based on Hartmann wavefront sensor, in: Proceedings of SPIE - The International Society for Optical Engineering, 2010. [17] J. Rice, Mathematical statistics and data analysis Cengage Learning, 2006. [18] W. Southwell, Wave-front estimation from wave-front slope measurements, J. Opt. Soc. Am. 70 (8) (1980) 998–1006.
Contents lists available at ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
Stitching error calculation method for the slope data on a square subaperture Guangrao Guo a, Dahai Li a,n, Kewei E a, Mengyang Li a, Haiping Chen b,n a b
School of Electronics and Information Engineering, Sichuan University, Chengdu 610065, China Research Center of Laser Fusion, China Academy of Engineering Physics, Mianyang 621900, China
art ic l e i nf o
a b s t r a c t
Article history: Received 30 November 2015 Received in revised form 9 May 2016 Accepted 15 May 2016
The data on subaperture in this paper is collected by a Hartmann test device or a Hartmann–Shack wavefront sensor. The slope data on each subaperture of the test optical component is obtained via the linear least squares fitting routine in this paper. Importantly, the error transfer formulas are deduced when the measurement data are polluted by the random noise, which will affect the uncertainty of stitching parameters, such as tilt or defocus. In particular, the accuracy formulas for stitching evaluation, which can be used to calculate the stitching error of each point on each subaperture for both the parallel mode and the serial mode, are derived mathematically. Therefore, this paper provides these formulas to estimate exactly the stitching error in each subaperture if the variance of random noise has been estimated and if the overlapping area is also given in advance. The results of simulation experiments show that the stitching accuracy formulas proposed are verified and they can be used to evaluate stitching accuracy. The reason why the error of the parallel stitching mode is less than that of the serial stitching mode is presented theoretically. & 2016 Elsevier B.V. All rights reserved.
Keywords: Subaperture Slope stitching Parallel stitching Serial stitching Large component Accuracy
1. Introduction With the development of science and technology, large optical components are widely used in space telescopes [1,2], laser fusion systems [3,4] and other aspects. In traditional optical measurements, the large optical surface is tested by a large-scale interferometer. However, the interference testing method requires a flat reference, and most interferometers are sensitive to the measurement environment. As a result, the cost of a large-scale interferometer is greatly increased, and good environment control is required for testing [5,6]. In 1982, the subaperture measuring method was first introduced in interferometric metrology by C.J. Kim at the Optical Sciences Center at the University of Arizona, [7]. Since the method was proposed, it has been widely used in the measurement of large aperture optical elements. More recently, subaperture stitching interferometry is becoming increasingly popular. In 1994, Otsubo proposed a global subaperture stitching method based on stitching parameters; in this method, the misalignment errors is simultaneously fitted and obtained over multiple overlapping subapertures [8]. In 2003, QED Technologies successfully developed an automated subaperture stitching interferometer. This interferometer can successfully test the surface of flat and spherical optical elements with n
Corresponding author. E-mail addresses: [email protected] (D. Li), [email protected] (H. Chen).
http://dx.doi.org/10.1016/j.optcom.2016.05.043 0030-4018/& 2016 Elsevier B.V. All rights reserved.
diameters of less than 200 mm [9]. Moreover, research studies have also focused on the analysis of the accuracy of subaperture stitching. To date, there are three types of methods to evaluate the accuracy of subaperture stitching technique. In the first type, the stitching results are compared with the results of the full aperture interferometry; this method of evaluation is the most simple and straightforward one [10]. In the second type, the effect of random noise and overlapping area on the stitching parameters is qualitatively analyzed. This method can only obtain the error of the stitching parameters between two adjacent subapertures and it cannot obtain the error of each point on each subaperture after the subaperture stitching procedure is completed [8,11]. In the third type, according to the knowledge of statistics, the proportion of the regression sum of squares of the phase difference to the total sum of squares of the phase difference is used to evaluate the accuracy of the stitching [12]. Different from the traditional subaperture stitching method based on wavefront measurements, the subaperture wavefront slope stitching test is based on a Hartmann–Shack wavefront sensor [13]. Hartmann–Shack sensors have the advantages of insensitivity to the air disturbance and vibration, rapid testing and high dynamic range. In addition, subaperture wavefront slope stitching tests can avoid the influence of the piston. Thus, this method has a unique advantage in the testing of a large aperture optical element. Due to these characteristics, an increasing number of scholars are studying the subaperture wavefront slope stitching testing method [14–16]. Currently, research studies on
G. Guo et al. / Optics Communications 380 (2016) 124–133
wavefront slope stitching are mainly focused on the establishment of the mathematical mode and qualitative analysis of the stitching error. However, no mathematical formula is available to evaluate stitching accuracy when random noise exists in measurement process. In this paper, we describe in detail the following: 1) the effect of random noise and overlapping area on the stitching parameters and 2) the error relationship between the stitching parameters and the adjusted slope. Importantly, the error transfer formula for the slope stitching is first determined. On this basis, the respective formulas for solving the accuracy of serial slope stitching mode and parallel slope stitching mode are derived. This paper shows that these formulas are verified by our simulation experiments. The RMS and PV values of the residual error for the reconstruction of the wavefront based on a given slope matrix in our simulation experiments demonstrate that the error of parallel stitching is statistically less than that of serial stitching.
2. Principles As mentioned in the introduction, a large-size component can be measured using the stitching method. In this manner, we can obtain a wider range of measurement area with good resolution using a HS with a limited aperture. The disadvantages of the stitching methods include error propagation, and we must move the sensor or the sample repeatedly to cover the entire surface. Therefore, the slope data on a large aperture consists of multiple slope matrices on different sub-regions. Due to the misalignment, the piston, tilt, and even defocus, will be inevitably introduced into the measurement data when we move the sensor from one place to another. To eliminate the interference of misalignment, the slope data on the overlapping area between two adjacent subapertures are used to obtain the corresponding stitching parameters by a linear least squares fitting routine. As a result, these parameters, such as coefficient of tilt or defocus, are obtained, and then this procedure is repeated to the end of the stitching task until all these parameters are obtained. Next, these parameters are treated to construct the stitched slope data on the whole component. For simplicity of discussion, as shown in Fig. 1, we assume that the two adjacent subapertures W1 and W2 have a common area due to their overlapping area, and the measured phase distributions over the two areas are denoted as Φ1 and Φ 2, respectively. The connection of the two adjacent areas is considered first. Let W1 be the reference plane, and let the relationship between the phase distributions Φ1 and Φ 2 in the common area fulfill the following equation:
Φ2 = Φ1 + P + Tx x + Ty y,
(1)
where P is the piston, and Tx and Ty are the coefficients of the tilt in the x and y directions, respectively. Considering the derivation of x -tilt and y -tilt in Eq. (1), we have
G2x − G1x = Tx G2y − G1y = Ty,
125
(2)
where G1x = ∂Φ1/∂x and G1y = ∂Φ1/∂y are the x -tilt and y -tilt measured slopes, respectively, over the area W1, and G2x = ∂Φ 2/∂x and G2y = ∂Φ 2/∂y are the x -tilt and y -tilt measured slopes over the area W2, respectively. The parameter P disappears in Eq. (2) after taking the derivative. In other words, the piston does not influence the measured slope on subapertures; thus, we only need to find the coefficient Tx of x -tilt and coefficient Ty of y -tilt. In general, a minimum number of sampling points in the common area are chosen to accomplish the calculation of the values of Tx and Ty within an acceptable stitching error because the existence of random noise during measurement will affect the calculation accuracy of parameters. Suppose there are n sampling points in their common area; in this case, we can obtain the corresponding equations from Eq. (2) as follows:
⎡ G2x (x1,y ) − G1x (x1, y ) ⎤ ⎡ ⎤ 1 1 ⎥ ⎢ 1⎥ ⎢ ⎢ G2x (x2,y2 ) − G1x (x2 , y2 ) ⎥ = ⎢ 1⎥ ⎡⎣ T ⎤⎦, ⎥ ⎢ ⋮⎥ x ⎢ ⋯ ⎥ ⎢ ⎥ ⎢ ⎣ G2x (xn , yn ) − G1x (xn , yn )⎦ ⎣ 1⎦
(3a)
⎡ G2y (x1,y ) − G1y (x1, y ) ⎤ ⎡ ⎤ 1 1 ⎥ ⎢ 1⎥ ⎢ ⎢ G2y (x2,y2 ) − G1y (x2 , y2 ) ⎥ = ⎢ 1⎥ ⎡ T ⎤. ⎥ ⎢ ⋮⎥ ⎣ y ⎦ ⎢ ⋯ ⎥ ⎢ ⎥ ⎢ ⎣ G2y (xn , yn ) − G1y (xn , yn )⎦ ⎣ 1⎦
(3b)
Accordingly, Tx and Ty are the averages of the differences in the left-hand side of Eqs. (3a) and (3b), Next, the adjusted results of G2x and G2y can be solved by using Eq. (2). Subsequently, the adjusted results of G2x and G2y on the W2 can be connected with the results of G1x and G1y on the W1 successfully; as a result, we obtain the stitched slope on the two adjacent areas. Next, consider another case. Both tilt and defocus are introduced in the measured slope data, and the relationship between the phase distributions Φ1 and Φ 2 in their common area fulfills the following equation:
Φ2 − Φ1 = P + Tx x + Ty y + D (x2 + y2 ),
(4)
where D is the coefficient of defocus. Likewise, considering the derivation of x -tilt and y -tilt in Eq. (4), we have
G2x − G1x = Tx + 2Dx x G2y − G1y = Ty + 2Dy y .
(5)
Similarly, suppose there are n sampling points in the common area. From Eq. (5), we obtain
⎡ G2x (x1,y ) − G1x (x1, y ) ⎤ ⎡ 1, x ⎤ 1 1 1 ⎥ ⎢ ⎢ ⎥ ⎢ G2x (x2,y2 ) − G1x (x2 , y2 ) ⎥ = ⎢ 1, x2 ⎥ ⎡⎢ Tx ⎤⎥, ⎥ ⎢ ⋮ ⎥ ⎣ 2D x ⎦ ⎢ ⋯ ⎥ ⎢ ⎢ ⎥ ⎣ G2x (xn , yn ) − G1x (xn , yn )⎦ ⎣ 1, xn ⎦
(6a)
⎡ G2y (x1,y ) − G1y (x1, y ) ⎤ ⎡ 1, y1 ⎤ 1 1 ⎥⎡ ⎥ ⎢ ⎢ ⎤ ⎢ G2y (x2,y2 ) − G1y (x2 , y2 ) ⎥ = ⎢ 1, y2 ⎥ ⎢ Ty ⎥. ⎥ ⎢ ⋮ ⎥ ⎢⎣ 2Dy ⎥⎦ ⎢ ⋯ ⎥ ⎥ ⎢ ⎢ ⎣ G2y (xn , yn ) − G1y (xn , yn )⎦ ⎣ 1, yn ⎦
(6b)
According to the least-squares method, we obtain the coefficients of Tx , Ty and Dx , Dy in Eq. (6a) and (6b). Next, the adjusted results of G2x and G2y can be obtained by using Eq. (5). As a result, the stitched slope on the two adjacent areas is obtained. 3. Error analysis Fig. 1. Schematic of subaperture stitching.
For an actual measurement, the slope data measured by a
126
G. Guo et al. / Optics Communications 380 (2016) 124–133
Hartmann–Shack sensor could be polluted by random noise to a certain extent. The accuracy of these coefficients related to the stitching parameters obtained by using Eqs. (3) and (6) are influenced by noise. In addition, the accuracy of the coefficients is also affected by the number of sample points in their common area. If these parameters have some errors, then the decrease of the accuracy of the adjusted slope is inevitable. In this paper, the influence of the two factors on the accuracy of stitched slope data is discussed mathematically in detail. The error analysis is conducted in the two respective cases of with or without defocus by using (Eqs. (3) and 6).
Therefore, each point in G″2x (xi , yi ) will have Gaussian noise with the variance described in Eq. (10). This result shows that the variance of T1x on the overlapping area between A1 and A2 is not equal to that of the adjusted slope on the other area of subaperture A2. In addition, the uncertainty on other area of subaperture A2 is increased by the variance of the T1x on the shaded area. Using the approach of data stitching between the subaperture A1 and A2 described above, the tilt coefficient T2x on the subaperture A3 can be written as
3.1. Slope data including the tilt coefficient
T2x =
We consider the slope data in the x direction as an example and suppose that the subapertures A1, A2 , ⋯ , AN in the measured surface are square shaped, as shown in Fig. 2; A1 is defined as the standard reference subaperture and the arrow indicates the path of stitching. The size of the slope matrix in each subaperture is W × W , and the shaded area with the size of the slope matrix W × M is the common area of the two adjacent subapertures. Let the two slope matrices in the x direction on the two adjacent subapertures A1 and A2 be G1x and G2x , respectively. We assume the RMS (Root Mean Square) value of the random noise in the measured slope of each subaperture is σ . In fact, in a stable test environment, the random noise can be statistically estimated through multiple measurements. For the stitching between the subaperture A1 and A2 , the unknown parameter T1x (the tilt coefficient in the x direction on A2) can be easily written according to the Eq. (3a) as follows:
Because the matrices G3x (xi , yi ), G′2x (xi , yi ) and T1x are mutually independent, the variance of T2x is
T1x =
1 WM
WM
∑ ⎡⎣ G2x (xi , yi ) − G1x (xi , yi ) ⎤⎦. i=1
(7)
Note that the slope matrix G1x (xi , yi ) and G2x (xi , yi ) are mutually independent, so the variance of T1x on the shaded area can be expressed as follows:
σT21x =
2σ 2 . WM
(8)
In this paper, to distinguish the variance of stitching parameters on the overlapping area (the shaded) between A1 and A2 from that of adjusted slope on the other area of subaperture A2, in which the shaded area is not included, we define two new slope matrices on the remaining area of subaperture A2 as G′2x (xi , yi ) and G″2x (xi , yi ) , which correspond to the portion of matrix G2x (xi , yi ) and the adjusted results, respectively. After removing the tilt coefficient T1x from the G′2x (xi , yi ), we obtain the slope matrix G″2x (xi , yi ), which can be written as
G″2x (xi , yi ) = G′2x (xi , yi ) − T1x
(9)
Note that the matrices G′2x (xi , yi ) and T1x are mutually independent as well; thus, the variance of the adjusted slope matrix G″2x (xi , yi ) can be obtained as follows:
σ 2G 2″x (xi, yi ) = σ 2 +
2 2 σ WM
(10)
WM ∑i = 1 ⎡⎣ G3x (xi , yi ) − G′2x (xi , yi ) ⎤⎦
σT22x =
WM
+ T1x
4σ 2 WM
(11)
(12)
Likewise, after removing the tilt coefficient T2x from G′3x (xi , yi ), we obtain the adjusted slope matrix G″3x (xi , yi ) as follows:
G″3x (xi , yi ) = G′3x (xi , yi ) − T2x
(13)
Note that the two matrices G′3x and T2x are mutually independent; thus, the variance of G″3x (xi , yi ) can be obtained as follows,
σ 2G3″x (xi, yi ) = σ 2 +
4σ 2 WM
(14)
Accordingly, based on the two Eqs. (10) and (14), we can see that the uncertainty of tilt coefficients due to the random noise will transfer from one to the next subaperture and is increased gradually for the subsequent one to the last one. As a result, the stitching accuracy decreases. Therefore, for the stitching of N subapertures, which are stitched sequentially from one to the next subaperture, we can conclude the general equations for the calculation of the variance of T(N − 1) x and the G″Nx (xi , yi ) separately are as follows,
σT2(N − 1) x = σ G2″
Nx (xi, yi )
2 (N − 1) σ 2 , WM
= σ2 +
(15)
2 (N − 1) σ 2 . WM
(16)
Here, Eqs. (15) and (16) show that the number of subapertures and the overlapping area should be properly adopted according to the requirement of accuracy for the slope measurement when we conduct the stitching by use of slope data provided by Hartmann– Shack sensor. Similarly, for the discussion regarding the measured slope data in the y direction, we can also obtain similar equations as Eqs. (15) and (16). Eq. (16) is the basic error transfer formula proposed in this paper for the slope stitching. By using Eq. (16), we can calculate the stitching error on each subaperture if the random noise is measured in advance and the overlapping area is given. 3.2. Slope data including tilt and defocus coefficients Sometimes, the defocus may disturb the measurement besides the tilt. For this situation, from the statistical theory, the Eq. (6a) is suitable for linear regression model [17], and the slope data in the x direction can be written as
li = 2 × βxi + β0 + εi , Fig. 2. Schematic of multiple subapertures stitching.
i = 1, 2, ⋯, n.
(17)
Here, li is defined as li = G2x (xi , yi ) − G1x (xi , yi ), and n is the
G. Guo et al. / Optics Communications 380 (2016) 124–133
number of sampling points in the common area. β0 and β are the regression constant and the regression coefficient, respectively. εi represents the influence of all the other random variables on the slope data li . Suppose Dx and Tx are the least squares estimation of the parameter β and β0 , respectively; in this case, the linear regression equation can be written as follows:
^ li = 2 × Dx xi + Tx
(18)
from the G′2x (xi , yi ), we can write G″2x (xi , yi ) on the other area of subaperture A2 as follows:
^ G″2x (xi , yi ) = G′2x (xi , yi ) − li .
σ G2″
2x (xi, yi )
^ L − Xb = V .
(19)
(
σT2x
(21)
^ According to statistical theory, the variance of li in Eq. (18) can be expressed as
σ ^2 = 4σ D2x x2 + σT2x + 4xCov (Dx , Tx ),
(22)
li
where
Dx
σ2
and
are the variances of the Dx and Tx , re-
Tx
)
6σ 2 , t 2W (M 3 − M ) ⎞ 2⎛ 1 12x¯ 2 = + 2 3 ⎜ ⎟ σ 2, W⎝M t (M − M ) ⎠
(20)
)−1X T L.
σ2
(
⎤ ⎥ 2 σ + σ 2. 2⎥ ⎦
σ D2x =
x1 ⎤ ⎡ v1⎤ ⎥ ⎢v ⎥ x2 ⎥ ^ ⎡ Tx ⎤ 2 ,b=⎢ ⎥, V = ⎢ ⎥, ⋮⎥ ⎣ 2D x ⎦ ⎢ ⋮⎥ ⎥ ⎣ vn ⎦ xn ⎦
and vi , which is the difference between the test value li and the ^ regression value li , represents the residual. Suppose the RMS value of the random noise in the measured slope distribution is σ , i.e., each point in G1x and G2x has measurement noise with the RMS value of σ , and the RMS value of li with Gaussian distribution is ^ 2 σ .The unknown b can be obtained by the following equation:
^ b = XTX
⎡ (xi − x¯ )2 ⎢1 = 2⎢ + 1 n n n ∑i = 1 xi2 − n ∑i = 1 xi ⎣
(30)
Further, assume that t is the physical interval between two sampling points in the x direction; thus, we can separately rewrite Eqs. (25), (26) and (30) as follows:
where
⎡ l1 ⎤ ⎡1 ⎢ ⎥ ⎢ 1 l2 ⎥ ⎢ L= , X=⎢ ⎢ ⋮⎥ ⎢⋮ ⎢ ⎥ ⎢ ⎣1 ⎣ ln ⎦
(29)
^ Note that G′2x (xi , yi ) and li are mutually independent, the variance of G″2x (xi , yi ) can be obtained as
The matrix form of the corresponding error equations is given by
127
σG2 2x (xi, yi ) ″
2⎡ 1 12 (x − x¯ )2 ⎤ 2 = + 2 3i ⎢ ⎥ σ + σ 2. W⎣M t (M − M ) ⎦
(31) (32) (33)
Eq. (33) demonstrates that the error of slope stitching is related to the random noise, the number of sampling points in common area, and the location of the sampling points. Likewise, the uncertainties of the tilt and defocus coefficients will be both transferred and are increased gradually from one to the next subaperture. As a result, the stitching accuracy will suffer more than for the case of only the tilt/tip existing. In addition, it is also more complicated to calculate the stitching accuracy on each subaperture. Similarly, the discussion of the slope stitching accuracy in the y direction will have the same result as Eq. (33). When we use the Hartmann–Shack wavefront sensor to measure the slope data on subapertures, generally, only the tilt coefficient is introduced into the measurement. Therefore, in the following analysis, the defocus is neglected.
spectively. Cov (Dx , Tx ) is the covariance between Dx and Tx . Solving Eq. (21), we obtain
(
Tx =
n
n
) (
n ∑i = 1 xi li − ∑i = 1 xi
Dx =
n
n
n ∑i = 1 li
n ∑i = 1 xi
2 [n ∑i = 1 xi2 − ( ∑i = 1 xi
(
n ∑i = 1 xi2
)(
n l i=1 i 2 ) ]
)( ∑
)−( − (∑
)( )
n ∑i = 1 xi li
2 n x i=1 i
n n ∑i = 1 xi2
4. Stitching mode
) (23)
) (24)
Therefore, the respective variances of Dx and Tx and the covariances between Dx and Tx can be obtained as follows,
σ D2x =
σT2x
1 n 2 ∑ n x i2 − ∑ n x i i=1 i=1
(
2
)
σ2 (25)
⎛ x¯ 2 ⎜1 = 2⎜ + n 1 n 2 n ∑i = 1 xi − n ∑i = 1 xi ⎝
(
Cov (Dx , Tx ) =
⎞ ⎟ 2 σ 2⎟ ⎠
)
n − ∑i = 1 xi n
(
2
n
)
n ∑i = 1 xi2 − ∑i = 1 xi
(26)
σ2 (27)
Substituting Eqs. (25), (26) and(27) into Eq. (22), the variance of
^ li is expressed as
⎡ (x − x¯ )2 ⎢1 σ ^2 = 2 ⎢ + 1 n n li n ∑i = 1 xi2 − n ∑i = 1 xi ⎣
(
n
⎤ ⎥ 2 σ . 2⎥ ⎦
)
(28)
Where x¯ = ∑i = 1 xi . Based on the above-described treatment of the stitching with tilt coefficient, after removing the tilt and defocus coefficients
Generally, at present, the stitching mode is serial stitching, as shown in Fig. 3(a), which has the advantage of a simple mechanical structure and easier operation. However, serial stitching requires more time to finish the stitching task, and, as we know, the error will accumulate gradually and cause a large stitching error to occur in the measurement. In this paper, the parallel stitching mode (as shown in Fig. 3(b)) is proposed and discussed. Unlike the serial stitching, it can stitch multi-subapertures simultaneously, and the parallel stitching method requires less time to finish a stitching task with reduced loss of accuracy. Thus, it is necessary to apply an appropriate stitching mode to achieve our goal with a smaller error, according to the path and sequence of stitching. By considering an example, the error of the serial stitching mode and parallel stitching mode will be analyzed in detail. For simplicity, a component model with 3 × 3 subapertures is established, as shown in Fig. 3. The size of the slope matrix in each subaperture is 10 × 10. In addition, the size of the slope matrix in the overlapping area is 2 × 10 or 10 × 2. Each element (point) of slope matrix in each subaperture has random noise of mean equal to 0 and RMS value of σ . 4.1. Serial stitching As shown in Fig. 3(a), suppose that the subaperture A1 is the reference aperture. The stitching treatment between the adjacent A1, A2 and A3 is the same as shown in Fig. 2 from the left to the
128
G. Guo et al. / Optics Communications 380 (2016) 124–133
Fig. 3. The two modes of 3 × 3 subapertures stitching. (a) Serial stitching (b) Parallel stitching.
right. Moreover, the stitching error can be calculated by using Eq. (16). Thus, we only need to discuss the stitching between the subaperture A4 and A5. For simplicity of calculation, we ignore the common area (the shaded area in Fig. 3(a)) between subapertures A2, A3, and A4 when the stitching between A4 and A3 is implemented; thus, the overlapping area becomes the matrix size of 2 × 8. According to Eq. (11), the tilt coefficient on A4 can be expressed as
T3x =
1 16
16
∑ ( G4x (xi , yi ) − G′3x (xi , yi ) ) + T2x. i=1
(34)
As a result, the variance of T3x can be determined to be equal to 0.325σ 2 by use of Eq. (15). Further, after canceling the tilt coefficient T3x from the G′4x (xi , yi ), we obtain the adjusted G″4x (xi , yi ), namely, G″4x (xi , yi ) = G′4x (xi , yi ) − T3x . Likewise, the variance of G″4x (xi , yi ) can be determined to be equal to 1.325σ 2 by use of Eq. (16). In the same way, when the stitching of slope data on A5, A2 and A4 is implemented, we ignore the data in the common area ( 2 × 2 sample points, i.e., the shaded area in Fig. 3(a)) between A2 , A4 and A5, i.e., the tilt coefficient on A5 can be expressed as T4x =
16 ∑16 i = 1 (G5x (xi, yi ) − G′4x (xi, yi ) + T3x ) + ∑k = 1 (G5x (xk , yk ) − G′2x (xk , yk ) + T1x ) 32
(35)
where subscript i denotes the overlap between A5 and A4 , and subscript k denotes the overlap between A5 and A2. For approximate calculation, we suppose the variance of each point on the entire subaperture A2 is σ 2G 2x (xi, yi ) ; in addition, because T1x and T3x ″ are not statistically independent, the variance of T4x is given by
σT24x =
σ2 1 1 1 + σT23x + σT21x + Cov (T1x, T3x ). 16 4 4 2
(36)
where Cov (T1x, T3x ), which is the covariance between the T1x and T3x , is expressed as Cov (T1x, T3x ) = ⎛ ⎞ ∑16 (G4x (xi, yi ) − G′3x (xi, yi )) + ∑16 k = 1 [(G3x (xk , yk ) − G′2x (xk , yk )) + T1x ] ⎟ Cov ⎜ T1x, i = 1 ⎜ ⎟ 16 ⎝ ⎠
(37)
Because G4x (xi , yi ) , G′3x (xi , yi ) , G′2x (xk , yk ) and T1x are statistically independent, the Cov (T1x, T3x ) is equal to σT21x . Substituting the quantity into Eq. (36), we obtain the variance of T4x , which is equal to 0.219 σ 2. Likewise, after removing the tilt coefficient T4x from G′5x (xi , yi ), we obtain the adjusted G″5x (xi , yi ), namely, G″5x (xi , yi ) =
G′5x (xi , yi ) − T4x ; thus, the variance of G″5x (xi , yi ) can be solved and is found to be equal to 1.219σ 2. As mentioned above, the stitching for the data on A7 is found to be similar to that on A4 , and the stitching for the data on A6 , A8 and A9 is similar to that on A5. Thus, we can approximately determine the stitching error of each point of the stitched slope on the test component if the variance of random noise has been estimated and the overlapping area is obtained in advance. Without exception, the method and formulas mentioned above are also applicable to the serial stitching with more subapertures. 4.2. Parallel stitching For parallel stitching, as shown in Fig. 3(b), the center subaperture A1 is considered as the reference aperture. If these subapertures have the same distance from the reference aperture, then their adjusted slope will have an identical stitching error. Thus, we only must discuss the stitching for the data between the four subapertures on the top right ( A1, A2a , A2b , A3a ). The tilt coefficient on A2a and A2b are defined as T1ax and T1bx , respectively, according to Eq. (7): T1a x = ∑20 i = 1 (G2a x (xi , yi ) − G1x (xi , yi )) /20, 20
= ∑i = 1 (G2b x (xi , yi ) − G1 x (xi , yi )) /20 . Due to their dependence, the covariance between T1a x and T1b x is expressed as Cov (T1a x, T1b x ). Because G2a x (xi , yi ), G2b x (xi , yi ) and G1 x (xi , yi ) are statistically independent on the other area of themselves and there are 2 × 2 points in the common area (the shaded area in Fig. 3(b)), accordingly, Cov (T1a x, T1b x ) can be solved as follows: T1b
x
⎡ ∑ 4 G (x , y ) ∑ 4 G (x , y ) ⎤ 1x i i 1x i i ⎥ σ2 = Cov (T1a x, T1b x ) = Cov ⎢ i = 1 , i =1 . ⎢ ⎥ 100 20 20 ⎣ ⎦
(38)
For simplicity of calculation, when the slope data on A3a is separately stitched with that on A2a and A2b , we already ignore the data on the common area ( 2 × 2 sample points, i.e., the shaded area in Fig. 3(b)) between A2a , A2b , and A3a . Thus, the tilt coefficient on A3a is given by T2a x =
16 16 ∑i = 1 (G3a x (x i , yi ) − G′2a x (x i , yi )) + ∑i = 1 (G3b x (x i , yi ) − G′2b x (x i , yi ))
32 +
T1a x T + 1b x 2 2
Because only T1a x and T1b parameter T2a x is given by
(39) x
are dependent, the variance of
G. Guo et al. / Optics Communications 380 (2016) 124–133
Fig. 4. The statistically simulated slope error of serial stitching.
σT22a x =
σT2 σT2 1 σ2 + 1a x + 1b x + Cov (T1a x, T1b x ) 16 4 4 2
(40)
Substituting Eq. (38) into Eq. (40), we obtain the variance of the parameter T2a x to be equal to 0.118 σ 2. Likewise, G′3a x (xi , yi ) is the slope data on the other area of A3a . After canceling the tilt coefficient T2a x from G′3a x (xi , yi ). we obtain Considering G″3a x (xi , yi ) = G′3a x (xi , yi ) − T2a x . G″3a x (xi , yi ): G′3a x (xi , yi ) and T2a x are statistically independent, the variance of G″3a x (xi , yi ) is solved to be equal to 1.118σ 2. Similarly, the stitching error of each point in the stitched slope can be calculated or predicted by this method if the variance of random noise has been estimated and the overlapping area is obtained. In addition, as in the case of serial stitching, the parallel stitching formulas above are also applicable to the case with more subapertures. 4.3. Simulation To verify the stitching error formulas for the two stitching modes, a simulation experiment of subapertures stitching is implemented. As mentioned above, a component with 3 × 3 subapertures is discussed, the size of the slope matrix is 10 × 10 for each subaperture, and the size of the slope matrix in the overlapping area is 2 × 10 or 10 × 2. Next, the error of each point of the stitched slope on a component is calculated using the error formulas proposed in this paper. The serial stitching simulation experiment procedure is as follows: 1. Generate a matrix of slope data with the order of magnitude of 1 × 10−6rad for a test optical component first. 2. A Gaussian random noise matrix with mean of 0 and RMS value of σ = 2.423 × 10−7rad is added to the above slope data of each subaperture. 3. An amount of tilt, which is used to simulate the misalignment of a slope sensor, is added to the slope data of each subaperture (except the subaperture A1). 4. Stitch the slope data on the all subapertures in serial stitching mode. Repeating the steps 2–4, the stitching for the 9 regions that constitute the entire stitched slope is implemented by using serial stitching mode. Next, the RMS of noise at each slope data in each region is statistically calculated after repeating the trials 100,000 times, as shown in Fig. 4, in which the color of each point reflects the value of noise, but they are nearly the same in one region. In addition, the average of RMS of noise at all points in each region is
129
compared with the theoretical value from our formulas (Eqs. (34)– (37) and (16)) proposed in the paper, as shown in Table 1. It can be observed that the statistical value on the different regions is almost completely equal to the predicted value on the corresponding regions by our formulas. In addition, the error of each point in the fourth region but not in the ninth region has the largest stitching error among the nine regions. Therefore, the result of this paper provides us with formulas to exactly determine the stitching error on each subaperture if the variance of random noise has been estimated and the overlapping area is given. Instead of serial stitching, the simulation procedure mentioned above is implemented for parallel stitching. Likewise, after repeating steps 2–4, the RMS of noise at each slope data in each region is statistically calculated after repeating the trials 100,000 times, as shown in Fig. 5. In addition, the average of RMS of noise at all points in each region is compared with the theoretical value from our formulas (Eqs. (38)–,(40) and (16)) proposed in the paper, as shown in Table 2. In addition, the error of each point in the third region has the largest error value among the nine regions. As a result, the formulas provided in this paper can be used to estimate the error in each subaperture if the variance of random noise has been estimated and the overlapping area is given in advance. Moreover, we can see that the error obtained using the parallel stitching mode is smaller than that obtained using the serial stitching mode. For further comparison, using the weighted average error of each point in the stitched slope as the reference to evaluate the stitching accuracy for the two modes, we define the error as N
error =
∑ σ Gi″ × pi .
(41)
i=1
where N is the number of the sub-regions in the component, σ Gi″ is the theoretical variance value of the error on the ith region, and pi is the ratio of the area of the ith region to that of the entire component. So thus, the error of the serial stitching and parallel stitching are 2.63 × 10−7rad and 2.53 × 10−7rad from Eq. (41), respectively. This result shows that the parallel stitching mode has a higher accuracy than the serial stitching mode. For an element with 3 × 3 subapertures and two rows or two columns overlapping area, with the increase of the random noise from 0 to 5 × 10−7rad , the weighted average errors of the two stitching modes at different noise are obtained using Eq. (41) and are illustrated in Fig. 6. We can see that the errors of both the serial stitching and parallel stitching linearly increase with increasing random noise, and the error of the serial stitching is greater than that of the parallel stitching. Similarly, given the random noise of 2.423 × 10−7rad , when the overlapping area between adjacent subapertures increases from 2 to 9 rows or columns, the errors of the two stitching modes with the different overlapping area are determined by Eq. (41) and are Table 1 Comparison of the mean of statistical errors and theoretical values in 9 regions by serial stitching. Region
The statistical value (100,000 times)
The theoretical value
1
2.42 × 10−7
2.42 × 10−7
2
2.54 × 10−7
2.54 × 10−7
3
2.66 × 10−7
2.66 × 10−7
4
2.79 × 10−7
2.79 × 10−7
5
2.66 × 10−7
2.68 × 10−7
6
2.56 × 10−7
2.56 × 10−7
7
2.70 × 10−7
2.70 × 10−7
8
2.69 × 10−7
2.69 × 10−7
9
2.75 × 10−7
2.74 × 10−7
130
G. Guo et al. / Optics Communications 380 (2016) 124–133
Fig. 5. The simulated slope error by using the parallel stitching mode.
Table 2 Comparison of the mean of statistical errors and theoretical values in 9 regions by parallel stitching. Region
The statistical value (100,000 times)
The theoretical value
1
2.42 × 10−7
2.42 × 10−7
2
2.54 × 10−7
2.54 × 10−7
3
2.56 × 10−7
2.56 × 10−7
Fig. 6. Comparison of the errors for two stitching modes when the random noise is increasing.
shown in Fig. 7. The errors of both the serial stitching and parallel stitching decrease with the increase of the overlapping area, and the error of serial stitching is greater than that of parallel stitching.
5. Wavefront reconstruction of A component In the following simulation, assume that there is a component with the sizes of 420 mm 420 mm and the slope matrices of the whole component in x - and y -directions have dimensions of 42 × 42. Based on the slope matrices without noise, the wavefront of this component is reconstructed by a zonal reconstruction method using the Southwell's slope geometry [18] and is shown in Fig. 8. For testing our method, the component is divided into 5 × 5 subapertures. The size of the slope matrix of each subaperture is 10 × 10, and the size of the overlapping area is 2 × 10 or 10 × 2. Each point of the slope matrix of each subaperture has random noise whose mean is 0 and RMS value s is 2.423 × 10−7rad . Next, the slope is stitched by the serial stitching and the parallel stitching, respectively. Subsequently, the wavefront is reconstructed from the stitched slope data. The simulation results are shown in Figs. 9 and 10 as follows, Wavefront residual refers to the reconstruction from the stitched slope minus that from the original slope. From Figs. 9 and 10, we can see the reconstructed wavefront from the stitched slope by both the stitching modes are nearly similar to that from the original slope. In addition, the residual data show the accuracy of parallel stitching is better than that of serial stitching. To further verify the conclusion, 1000 simulation trials for wavefront reconstruction were implemented by using the two stitching modes. In each trial, the PV and RMS of the wavefront residual from the parallel stitching mode are subtracted from that of the serial stitching mode, and the distribution of their differences are shown in Fig. 11(a) and (b), respectively. The number of “*” points above the red straight line is observed to be much higher than that under this line in the two figures. From the statistical results, we can conclude that the error of the serial stitching is higher than that of the parallel stitching.
Fig. 7. The comparison of the errors for the two stitching modes when the overlapping area is increasing.
6. Conclusions In this paper, the mathematical model of the slope stitching based on the least squares fitting routine is presented. Using the model the stitching parameters, such as tilt and defocus, can be determined. After the parameters are removed, the surface figure of
Fig. 8. Wavefront construction from the original slope data.
G. Guo et al. / Optics Communications 380 (2016) 124–133
131
Fig. 9. The simulated results from the parallel stitching mode: (a) wavefront estimation from stitching slope and (b) wavefront residual.
Fig. 10. The simulated results from the serial stitching mode: (a) wavefront estimation from stitching slope and (b) wavefront residual.
Fig. 11. Statistical comparison of the two modes: (a) the PV difference of wavefront residual and (b) the RMS difference of wavefront residual.
the test component can be retrieved using a zonal algorithm or a model approach. More importantly, the error transfer formulas are determined when the variance of random noise has been estimated
and the overlapping area is given in advance. Our discussion indicated that the uncertainty of tilt or defocus due to the random noise will transfer from one subaperture to the next subaperture and
132
G. Guo et al. / Optics Communications 380 (2016) 124–133
is increased gradually; such uncertainty will cause the decrease of stitching accuracy. Mathematically, the formulas for stitching accuracy evaluation can be used to calculate the stitching error of each point on each subaperture for the parallel mode and the serial mode. The results of the simulation experiments demonstrated that all the stitching accuracy formulas proposed are verified and that these formulas have the ability to evaluate the stitching accuracy. In practice, the number of subapertures and the overlapping areas should be properly adjusted according to the requirement of accuracy for the slope measurement when we conduct the stitching by use of slope data provided by a Hartmann–Shack sensor. Moreover, the accuracy of these two stitching modes is compared, and the simulation experiment and theoretical analysis show that the error of parallel stitching is less than that of serial stitching.
n
n x2 i=1 i
(∑
Tx =
n
) (
n l i=1 i 2 ) ]
)( ∑
n
2 [n ∑i = 1 xi2 − ( ∑i = 1 xi n l i=1 i n ∑ n i = 1 xi2
)( ∑
) (A.4)
n x i=1 i n ∑ ( i = 1 xi )2
) − (∑ −
n xl i=1 i i
)( ∑
) (A.5)
From the simplification of the formula of Eq. (A.4), we obtain n ⎡ n ⎤ 1 ∑i = 1 ⎣ (nxi − ∑i = 1 xi ) li ⎦ n n 2 2 2 n ∑i = 1 xi − ( ∑i = 1 xi )
Dx =
(A.6)
The process of solving the variance of Dx is as follows:
σ D2x =
1 4
Funding information
= This work was supported by the National Defense Basic Research Program of China (JSJC2013212C002) and the National High Technology Research and Development Program of China (2015AA015902).
=
n ⎡ n ∑i = 1 ⎢⎣ nxi − ∑i = 1 xi
⎤ × 2σ 2⎥⎦ 2 ⎤2 ⎡ n n 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
2
)
(
)
2 2⎤ ⎡ 2 n 2 n n n ∑i = 1 xi − 2n ∑i = 1 xi + n ∑i = 1 xi ⎥⎦ × 2σ 2 1 ⎢⎣ 2 ⎤2 ⎡ 4 n n 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
n 1 2 n ∑n x 2 − ∑n xi i=1 i=1 i
(
) (
2
)
(
)
)
σ2 (A.7)
Simplifying the formula of Eq. (A.5), we obtain
Appendix A From Eq. (20), we obtain
(
n
(
n ∑i = 1 xi li − ∑i = 1 xi
Dx =
n n n ∑i = 1 ⎡⎣ ( ∑i = 1 xi2 − xi ∑i = 1 xi ) li ⎤⎦
Tx =
n ⎤ ⎡ n ⎢ ∑ xi2 − ∑ xi ⎥ ⎢ 1 −1 i=1 i=1 ⎥ XTX = ⎥ ⎢ n n n 2 2 n ∑i = 1 xi − ( ∑i = 1 xi ) ⎢ ⎥ x n ⎥ ⎢− ∑ i ⎦ ⎣ i=1
n
n
n ∑i = 1 xi2 − ( ∑i = 1 xi )2
(A.8)
The process of solving for the variance of
)
σT2x
(A.1)
⎡ n ⎤ ⎢ ∑ li ⎥ ⎢ i=1 ⎥ XTL = ⎢ n ⎥ ⎥ ⎢ ⎢ ∑ xi ⎥ ⎣ i=1 ⎦
=
(
)
2 ⎤2 ⎡ n n 2 ⎣⎢ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
⎡ n 2σ 2 × ⎣⎢ n ∑i = 1 xi2
(
=
)
2
)
−
n ∑i = 1 xi2 2 ⎤2 ⎡ n n 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
⎡⎛ n ⎞⎛ n ⎞ ⎛ n ⎞⎛ n ⎞⎤ ⎢ ⎜ ∑ x 2 ⎟ ⎜ ∑ li ⎟ − ⎜ ∑ x i ⎟ ⎜ ∑ x i li ⎟ ⎥ i ⎟⎜ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎥ ⎢⎝ ⎠⎝ i=1 ⎠ ⎝ i=1 ⎠⎝ i=1 ⎠ 1 i=1 ^ ⎥ ⎢ b= n n 2 2 n n n ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎥ n ∑i = 1 xi − ( ∑i = 1 xi ) ⎢ ⎥ ⎢ n ⎜⎜ ∑ xi li ⎟⎟ − ⎜⎜ ∑ xi ⎟⎟ ⎜⎜ ∑ li ⎟⎟ ⎥⎦ (A.3) ⎢⎣ ⎝ i=1 ⎠ ⎝ i=1 ⎠⎝ i=1 ⎠
n 2 ∑i = 1 xi2
⎡ n 2 ⎣⎢ n ∑i = 1 xi
(A.2)
Solving Eq. (21), we obtain
is as follows:
n n n ∑i = 1 ⎡⎣ ∑i = 1 xi2 − xi ∑i = 1 xi2 × 2σ 2⎤⎦
(
=
Tx
n x2 i=1 i
(∑ − (∑
n
+ ∑i = 1 xi2
n x2 i=1 i
)( ∑
) ⎤⎦⎥
2 ⎤2 n x ⎥ i=1 i ⎦
)
× 2σ 2
)
⎛ x¯ 2 ⎜1 = 2⎜ + 1 n n 2 n ∑i = 1 xi − n ∑i = 1 xi ⎝
(
⎞ ⎟ 2 σ 2⎟ ⎠
)
(A.9)
The process of solving for the covariance of Dx and Tx is as follows:
Dx and Tx can be obtained as
n ⎡ n n ⎛ ⎤⎞ n n 2 ⎜ 1 ∑i = 1 nxi − ∑i = 1 xi li ∑i = 1 ⎣ ∑i = 1 xi − xi ∑i = 1 xi li ⎦ ⎟ cov (Dx , Tx ) = cov ⎜ , ⎟⎟ 2 n n ⎜ 2 n ∑n x 2 − ∑n x 2 n ∑i = 1 xi2 − ∑i = 1 xi ⎝ ⎠ i=1 i i=1 i n ⎡ n n n n ⎡ n n n 2 2⎤ 2 2⎤ 1 ∑i = 1 ⎣ nxi − ∑i = 1 xi ∑i = 1 xi − xi ∑i = 1 xi σ li ⎦ 1 ∑i = 1 ⎣ nxi − ∑i = 1 xi ∑i = 1 xi − xi ∑i = 1 xi σ li ⎦ = = 2 ⎤2 2 ⎤2 ⎡ ⎡ 2 2 n n n n 2 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦ ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
(
(
(
) )
(
)(
1 2
)
)
(
=
)
(
)
)
(
)
(
⎡ n n n n n 2σ 2 × ⎢⎣ n ∑i = 1 xi )( ∑i = 1 xi2 − n ∑i = 1 xi )( ∑i = 1 xi2 − n ∑i = 1 xi 2 ⎤2 ⎡ n n 2 ⎢⎣ n ∑i = 1 xi − ∑i = 1 xi ⎥⎦
(
)(
) (
(
)
n x2 i=1 i
)( ∑
n x i=1 i
) + (∑
)
2⎤ n x ⎦ i=1 i ⎥
)( ∑
)
n
=
− ∑i = 1 xi n n ∑i = 1 xi2
(
n
− ∑i = 1 xi )
2
)
σ2 (A.10)
G. Guo et al. / Optics Communications 380 (2016) 124–133
References [1] B.E. Catanzaro, J.A. Thomas, E.J. Cohen, Comparison of full-aperture interferometry to subaperture stitched interferometry for a large-diameter fast mirror, in: Presented at the International Symposium on Optical Science and Technology, 2001. [2] Chunyu Zhao, J.H. Burge. Stitching of off-axis sub-aperture null measurements of an aspheric surface, in: Proceedings of SPIE – The International Society for Optical Engineering, 2008, pp. 706316–706316-7. [3] Michael Bray, Stitching interferometer for large optics: recent developments of a system, in: Presented at the Third International Conference on Solid State Lasers for Application to Inertial Confinement Fusion, 1999. [4] Yingjie Yu, Mingyi Chen, Correlative stitching interferometer and its key techniques, in: Proceedings of International Symposium on Optical Science and Technology. International Society for Optics and Photonics, 2002. [5] B. Koehler, M. Kraus, J. Moresmau, K. Wirenstrand, M. Duchateau, P. Duhoux, R. Karban, C. Flebus, E. Gabriel, and O. Pirnay, The VLTI Auxiliary Telescopes: commissioning of AT1, in: Proceedings of SPIE Astronomical Telescopes þ Instrumentation International Society for Optics and Photonics, 2004, p. 600. [6] Franziska Harms, Jürgen Wolf, Patrick Waddell, Edward Dunham, Andreas Reinacher, Ulrich Lampater, Holger Jakob, Lisa Bjarke, Sybil Adams, Randy Grashuis, Allan Meyer, Kenneth Bower, Keith Schweikhard, Thomas Keilig, On sky testing of the SOFIA telescope in preparation for the first science observations, in: Proceedings of SPIE Optical Engineering þ Applications. International Society for Optics and Photonics, 2009, pp. 745303–745303-15. [7] C.J. Kin, J.C. Wyant, Subaperture test of a large flat or a fast aspheric surface, J. Opt. Soc. Am. 71 (1981) 1587.
133
[8] M. Otsubo, K. Okada, J. Tsujiuchi, Measurement of large plane surface shapes by connecting small-aperture interferograms, Opt. Eng. 33 (2) (1994) 608. [9] P. Murphy, J. Fleig, G. Forbes, D. Miladinovic, G. Devries, S. O'Donohue, Subaperture stitching interferometry for testing mild aspheres, in: Presented at the SPIE Optics þ Photonics, 2006. [10] ChaoWen Liang, HungSheng Chang, PoChih Lin, ChengChung Lee, YiChun Chen, Vibration modulated subaperture stitching interferometry, Opt. Express 21 (15) (2013) 18255. [11] Rongzhu Zhang, C. Yang, Q. Xu, Precision analyses of a stitching interferometer testing system, Appl. Opt. 45 (11) (2006) 2399–2403. [12] Guangda Zhan, Xudong Xu, Zhengxiang Shen, Xiaoqiang Wang, ZhanShan Wang, Huasong Liu, Yiqin Ji, A global subaperture stitching algorithm and its precision evaluation, 2012-01-01, 2012. [13] D.R. Neal, J. Copland, D.A. Neal, Shack–Hartmann wavefront sensor precision and accuracy, 2002-01-01, 2002. [14] D.R. Neal, Apparatus and method for evaluating a target larger than a measuring aperture of a sensor. [15] T.D. Raymond, D.R. Neal, D.M. Topa, T.L. Schmitz, High-speed noninterferometric nanotopographic characterization of Si wafer surfaces, Nanoscale Opt. Appl. 4809 (2002) 208–216. [16] A.S.Y. He, G. Tang, Subaperture stitching based on Hartmann wavefront sensor, in: Proceedings of SPIE - The International Society for Optical Engineering, 2010. [17] J. Rice, Mathematical statistics and data analysis Cengage Learning, 2006. [18] W. Southwell, Wave-front estimation from wave-front slope measurements, J. Opt. Soc. Am. 70 (8) (1980) 998–1006.