Two-dimensional stitching interferometry for self-calibration of high-order additive systematic errors

Stitching interferometry is performed by collecting interferometric data from overlapped sub-apertures and stitching these data together to provide a full surface map. The propagation of the systematic error in the measured subset data is one of the main error sources in stitching interferometry for accurate reconstruction of the surface topography. In this work, we propose, using the redundancy of the captured subset data, two types of two-dimensional (2D) self-calibration stitching algorithms to overcome this issue by in situ estimating the repeatable high-order additive systematic errors, especially for the application of measuring X-ray mirrors. The first algorithm, called CS short for “Calibrate, and then Stitch”, calibrates the high-order terms of the reference by minimizing the de-tilted discrepancies of the overlapped subsets and then stitches the reference-subtracted subsets. The second algorithm, called SC short for “Stitch, and then Calibrate”, stitches a temporarily result and then calibrates the reference from the de-tilted discrepancies of the measured subsets and the temporarily stitched result. In the implementation of 2D scans in xand y-directions, step randomization is introduced to generate nonuniformly spaced subsets which can diminish the periodic stitching errors commonly observed in evenly spaced subsets. The regularization on low-order terms enables a highly flexible option to add the curvature and twist acquired by another system. Both numerical simulations and experiments are carried out to verify the proposed method. All the results indicate that 2D high-order repeatable additive systematic errors can be retrieved from the 2D redundant overlapped data in stitching interferometry. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
The surface topographic measurement based on interferometry is a major tool for optical surfaces inspection [1,2]. As the accuracy and precision requirement is getting higher while the measuring range and size being larger, the "standard" interferometer Field Of View (FOV) becomes insufficient for some industrial applications, e.g. to measure large or long and curved mirrors. To access a bigger measuring size with a limited interferometer FOV, the Sub-Aperture (SA) stitching interferometry [3][4][5] was developed as a low cost solution to flexibly extend the measuring range and size, while preserving the native lateral resolution. Stitching interferometry is one of the high-precision optical surface measurement techniques and it is applied in many fields [3,4,[6][7][8][9][10].
Stitching interferometry acquires many highly-overlapped SAs. Using the high redundancy of the overlapped data in the acquisition process, it is possible to retrieve additional information [11][12][13]. In some applications, the Surface Under Test (SUT) has the same quality as the reference surface of the interferometer [8,9]. Therefore, it is a practical and meaningful task to carefully design the stitching strategy from data acquisition to stitching algorithm. By doing this, we can fully use the redundant measurement data to get an accurate estimation of the systematic error. Polack et. al. use Legendre polynomial interpolation of the reference to solve the 1D stitching problem in the case where the stitching steps are not an integer number of pixels [14]. Nicolas et. al. provided a solution to reconstruct a one-dimensional (1D) profile of the surface and extract the 1D additive systematic error in the stitching process [15]. These 1D solutions [14,15] based on Legendre polynomials or cubic B-splines offer the capability of using the redundant stitching data in a practical manner and paths a way for two-dimensional (2D) cases. The 2D self-calibration stitching method for the general optics has been studied [16][17][18][19][20], which show the effectiveness in reducing the reference errors.
In this work, we present two types of 2D stitching algorithms capable of estimating high-order additive systematic errors of the interferometer with linear least squares method, especially for measuring X-ray mirrors. The scanning steps during the data acquisition are specially designed for reference retrieval to avoid unwanted periodic errors. We validate the proposed method with simulation using a known reference as a benchmark. In parallel, experiments are also carried out to demonstrate the feasibility of the proposed method to retrieve reference during the stitching measurement process. The retrieved 2D reference map is compared with the calibrated reference data. The capability and limitations of the proposed method are discussed for possible future improvements.

Mathematical model of 2D data stitching
In order to clearly describe the stitching problem, a mathematical model of the 2D stitching is established. As shown in Fig. 1, the measurement of nth subset among the total N subsets m n (x, y) is taken by an interferometer at location (x n , y n ) during the stitching acquisition process.
Stitching interferometry acquires many highly-overlapped SAs. Using the high redundancy of the overlapped data in the acquisition process, it is possible to retrieve additional information [11][12][13]. In some applications, the Surface Under Test (SUT) has the same quality as the reference surface of the interferometer [8,9]. Therefore, it is a practical and meaningful task to carefully design the stitching strategy from data acquisition to stitching algorithm. By doing this, we can fully use the redundant measurement data to get an accurate estimation of the systematic error. Polack et. al. use Legendre polynomial interpolation of the reference to solve the 1D stitching problem in the case where the stitching steps are not an integer number of pixels [14]. Nicolas et. al. provided a solution to reconstruct a one-dimensional (1D) profile of the surface and extract the 1D additive systematic error in the stitching process [15]. These 1D solutions [14,15] based on Legendre polynomials or cubic B-splines offer the capability of using the redundant stitching data in a practical manner and paths a way for two-dimensional (2D) cases. The 2D self-calibration stitching method for the general optics has been studied [16][17][18][19][20], which show the effectiveness in reducing the reference errors.
In this work, we present two types of 2D stitching algorithms capable of estimating high-order additive systematic errors of the interferometer with linear least squares method, especially for measuring X-ray mirrors. The scanning steps during the data acquisition are specially designed for reference retrieval to avoid unwanted periodic errors. We validate the proposed method with simulation using a known reference as a benchmark. In parallel, experiments are also carried out to demonstrate the feasibility of the proposed method to retrieve reference during the stitching measurement process. The retrieved 2D reference map is compared with the calibrated reference data. The capability and limitations of the proposed method are discussed for possible future improvements.

Mathematical model of 2D data stitching
In order to clearly describe the stitching problem, a mathematical model of the 2D stitching is established. As shown in Fig. 1, the measurement of nth subset among the total N subsets m n (x, y) is taken by an interferometer at location (x n , y n ) during the stitching acquisition process. The nth subset can be described by the expression where x n and y n are the in-plane translation amount in x and y directions to move the nth subset z(x + x n , y + y n ) inside a single SA. In general, r(x, y) stands for the repeatable additive systematic error of the interferometer, and here it is mainly represented by the interferometer reference error inside the chosen SA. The geometric parameters t n = [a n , b n , c n ] denote the x-tilt, y-tilt and piston of the nth subset with t n · [x, y, 1] = a n x + b n y + c n from the motion errors of the translations stage and the angular adjustment to null fringes. The last term n a (x, y) stands for the additive random noise in the measurement. The nth subset can be described by the expression where x n and y n are the in-plane translation amount in x and y directions to move the nth subset z(x + x n , y + y n ) inside a single SA. In general, r(x, y) stands for the repeatable additive systematic error of the interferometer, and here it is mainly represented by the interferometer reference error inside the chosen SA. The geometric parameters t n = [a n , b n , c n ] denote the x-tilt, y-tilt and piston of the nth subset with t n · [x, y, 1] = a n x + b n y + c n from the motion errors of the translations stage and the angular adjustment to null fringes. The last term n a (x, y) stands for the additive random noise in the measurement.

Two self-calibration stitching algorithms
Two types of self-calibration stitching algorithms are presented. In the first self-calibration stitching algorithm, the reference r is first calibrated and then stitch the reference-removed subsets m − r into a whole piece of z. We name this algorithm CS which is short for "Calibrate, and then Stitch." The second self-calibration stitching algorithm, called SC short for "Stitch, and then Calibrate", stitches the ith updated reference-removed subsets m − r i to get z i first, and then estimate the reference r i+1 for the next iteration. It is interesting to notice that the iterative operation has also been studied and applied for the reference estimation with three-flat test [21,22]. Both of CS and SC algorithms can be applied to any kind of underlying stitching strategy. Since we solve the stitching problem based on matrix operation in this work, the proposed self-calibration stitching algorithms (CS and SC) will be mainly addressed to cooperate with the pixel-relation-based stitching strategy and the subset-relation-based stitching strategy [23]. In pixel-relation-based stitching strategy, all relations between the measured data m, t, z, and r (if involved) are established at each pixel, while the subset-relation-based stitching strategy establishes those dependant equations with each subset.
When the CS algorithm cooperating with the pixel-relation-based stitching strategy named CS-P as illustrated in Fig. 2(a), we can consider that the SUT shape is cancelled out in the overlapped height difference d. Similar to the 1D case [13], its relations with t and r can be established by matching the corresponding pixels as where G is the geometric relation matrix. In fact, we are not interested in the tip, tilt and piston of the stitched z, so the matrix G can only have 3(N − 1) columns by simply regularizing t 1 = [0, 0, 0]. D is the pixel relation matrix carefully composed to match pixels in the sheared reference maps. D has M columns, when there are M pixels in an SA. In this situation, t and r can be estimated together directly with no iteration as After that, the stitching result z can be obtained by simply merging the reference-subtracted data m − r with the known t.
If the subset-relation-based stitching strategy is used with the CS algorithm, named CS-S, the initial reference r 0 can be set as zeros. From the de-tilted height differenced i in the overlaps with m − r i in the ith iteration, only the reference is estimated through a compensatory update r i+1 = r i +r i , where the reference compensation amountr i is estimated from the de-tilted discrepanciesd i by D ·r i −d i → min.
(4) By solving this linear least squares problem,r i is calculated usinĝ height differenced i again in the next iteration unit it meets the termination criterion. The resultant height z is stitched from the reference-removed data m − r i+1 with the subset-relation-based stitching strategy. The self-calibration stitching problem is divided into two sub-problems in the CS-S algorithm: the reference estimation and the stitching with a known reference. The essential idea of the SC algorithm in Fig. 2(c) is based on the following consideration. If the estimated reference r is correct, the discrepancy d s between the reference-removed subsets m − r and the stitched height z in corresponding SA should only be a tilted plane with random noise. The superscribe "s" denotes the operation in SA. The discrepancy d s can be de-tilted and symbolized asd s . We can make iterative estimations of the reference r by minimizing this de-tilted discrepanciesd s with a least squares method. As illustrated in Fig. 2(c), the reference r 0 is initialized with zeros. For the ith iteration, z i is stitched from the reference-subtracted subsets m − r i . The reference compensation amountr i is then estimated from the de-tilted discrepancieŝ d s i between each m n , n ∈ [1, N] and the stitched height z i with the corresponding pixels inside the nth SA as where D s is the relation matrix. In this matrix based algorithm, some matrix operations are carefully taken when composing D s to match pixels in all subsets with the ones of the stitched map in the corresponding SAs. The size of D s can be fixed as M × M by collapsing matrices describing different SAs into one. As the size of D s becomes independent of N, this operation will be useful when stitching a large number of SAs. By solving this linear least squares problem, r i is calculated usingr Similarly, once i is larger than i max , or the standard deviation ofr i is smaller than thr, it ends up with the resultant z i and r i . Otherwise, the loop continues with the newly updated the estimated reference r i+1 = r i +r i . Otherwise, the loop continues with the newly updated reference r i+1 which will be removed from the raw measurement data m to re-calculate de-tilted height differenced i again in the next iteration unit it meets the termination criterion. The resultant height z is stitched from the reference-removed data m − r i+1 with the subset-relation-based stitching strategy. The self-calibration stitching problem is divided into two sub-problems in the CS-S algorithm: the reference estimation and the stitching with a known reference.
The essential idea of the SC algorithm in Fig. 2(c) is based on the following consideration. If the estimated reference r is correct, the discrepancy d s between the reference-removed subsets m − r and the stitched height z in corresponding SA should only be a tilted plane with random noise. The superscribe "s" denotes the operation in SA. The discrepancy d s can be de-tilted and symbolized asd s . We can make iterative estimations of the reference r by minimizing this de-tilted discrepanciesd s with a least squares method. As illustrated in Fig. 2(c), the reference r 0 is initialized with zeros. For the ith iteration, z i is stitched from the reference-subtracted subsets m − r i . The reference compensation amountr i is then estimated from the de-tilted discrepancieŝ d s i between each m n , n ∈ [1, N] and the stitched height z i with the corresponding pixels inside the nth SA as where D s is the relation matrix. In this matrix based algorithm, some matrix operations are carefully taken when composing D s to match pixels in all subsets with the ones of the stitched map in the corresponding SAs. The size of D s can be fixed as M × M by collapsing matrices describing different SAs into one. As the size of D s becomes independent of N, this operation will be useful when stitching a large number of SAs. By solving this linear least squares problem, r i is calculated usingr Similarly, once i is larger than i max , or the standard deviation ofr i is smaller than thr, it ends up with the resultant z i and r i . Otherwise, the loop continues with the newly updated r i+1 = r i +r i in the next iteration. As a result, the interferometric measurements can be best explained by the stitched result z and the estimated reference r using these two types of self-calibration stitching algorithms. In the next section, we will address some possible ambiguity issues in this stitching model and self-calibration estimation followed by our corresponding strategies to resolve them.

Ambiguities and regularization in data acquisition and processing
As shown in Eq. (1), the measured m n (x, y) is a combination of three unknowns (SUT shape z, reference r, and geometric parameters t). Ambiguities may happen under certain conditions which makes the stitching problem become ill-posed. To resolve the issues, some regularization needs to be applied in the stitching strategy.

Regularization for periodic errors in two directions
In many stitching applications, the xand y-translation step sizes ∆x and ∆y are selected as a uniform value in one direction during the stitching process. However, using these uniform stepped stitching, the stitched height map will suffer with periodic errors in both SUT stitching and reference reconstruction, if we want to reconstruct the reference from the redundant data.
Let's consider there is situation in which z(x, y) and r(x, y) are replaced by another pair of SUT shape and reference shape z * (x, y) and r * (x, y) given by the following expression where p(x) and q(y) are periodic functions with periods T x and T y in x and y as If the step sizes are uniform, x n = n · ∆x and y n = n · ∆y. Then choosing T x = ∆x/n x with n x ∈ N and T y = ∆y/n y with n y ∈ N, we have x n = nn x · T x and y n = nn y · T y . Since nn x ∈ N and nn y ∈ N, , and q(y − y n ) = q(y − nn y · T y ) = q(y). According to the equations above, we have While n is an arbitrary index in [1, N], different pairs of SUT and reference (z(x, y) with r(x, y)) and (z * (x, y) with r * (x, y)) yield to the same measurement data m n (x, y). It is important to note that it is impossible for the stitching algorithm to distinguish one pair of solutions from another, since all of them can perfectly explain the measurement data everywhere at N subset locations. This periodic error issue is due to the uniform stitching steps, in order to overcome this problem, we suggest to make randomized nonuniform steps for both xand y-directions during the data acquisition process to make this ill-posed problem solvable. As illustrated in Fig. 3, the periodic repeated pattern in Fig. 3(a) is "broken" by the nonuniform steps by randomization in Fig. 3(b). Fig. 3. In 2D stitching interferometery, uniform steps (a) and randomized nonuniform steps (b) are two kinds of possible SA stepping strategies. We suggest using nonuniform steps (b) to avoid the periodic errors when estimating additive systematic errors.

(a) Sub-apertures with uniform steps (b) Sub-apertures with nonuniform steps
in x-and y-directions to cover the entire SUT. We can make sure at least one randomization solution is available by keeping the step number larger than the absolute value of the residual pixel number. Row-by-row or column-by-column vectors of potential step size are added with random variations. Because the random variations are generated in a certain range determined from the previous steps, the step sizes in each vector are not absolutely random. Consequently, we randomly sample the step sizes from one calculated step vector to make the resultant steps completely random.
With this regularization implemented in data acquisition, the ambiguous periodic errors created on SUT and reference by the self-calibration stitching algorithms can be significantly reduced.

Regularization for curvature and twist ambiguities
Despite the periodic errors when using uniform stitching steps, the curvature and twist are also ambiguous in the mathematical model described in Eq. (1). We consider three second order terms (x-curvature, y-curvature and twist) adding to z(x, y) and r(x, y) to compose the new shapes of SUT z * (x, y) and reference r * (x, y) as where the coefficients C 11 , C 12 , and C 22 ∈ R can be arbitrary values. The linear terms a n x + b n y + c n in Eq. (1) are replaced by a * n x + b * n y + c * n whose coefficients are a * n = a n + 2C 11 x n + C 12 y n , The measured data acquired by the interferometer will be given by Fig. 3. In 2D stitching interferometery, uniform steps (a) and randomized nonuniform steps (b) are two kinds of possible SA stepping strategies. We suggest using nonuniform steps (b) to avoid the periodic errors when estimating additive systematic errors.
In our implementation, the stitching step sizes are always integer pixels, and we first calculate the preferred step sizes by the user-preferred overlapping ratios and the desired stitching length in xand y-directions to cover the entire SUT. We can make sure at least one randomization solution is available by keeping the step number larger than the absolute value of the residual pixel number. Row-by-row or column-by-column vectors of potential step size are added with random variations. Because the random variations are generated in a certain range determined from the previous steps, the step sizes in each vector are not absolutely random. Consequently, we randomly sample the step sizes from one calculated step vector to make the resultant steps completely random.
With this regularization implemented in data acquisition, the ambiguous periodic errors created on SUT and reference by the self-calibration stitching algorithms can be significantly reduced.

Regularization for curvature and twist ambiguities
Despite the periodic errors when using uniform stitching steps, the curvature and twist are also ambiguous in the mathematical model described in Eq. (1). We consider three second order terms (x-curvature, y-curvature and twist) adding to z(x, y) and r(x, y) to compose the new shapes of SUT z * (x, y) and reference r * (x, y) as where the coefficients C 11 , C 12 , and C 22 ∈ R can be arbitrary values. The linear terms a n x+b n y+c n in Eq. (1) are replaced by a * n x + b * n y + c * n whose coefficients are a * n = a n + 2C 11 x n + C 12 y n , The measured data acquired by the interferometer will be given by Again, different solutions with arbitrary coefficients C 11 , C 12 , and C 22 yield to identical measurement m n (x, y) for every subset, as n is arbitrary. Infinite combinations of SUT, reference, and motions can perfectly explains the measured data. As a result, with the linear translations in the stitching motion, there is no self-calibration stitching algorithm able to separate the 2nd order terms on the SUT and the reference surfaces. Other measurement setups such as skip-flat test could become possible ways to overcome this issue, but it is out of the scope of this work. Unlike the periodic errors in Section 1 coupling with SUT and reference which occurs under uniform step condition only, the curvature and twist ambiguity always exists in the stitching model and couples not only with the curvature and twist on SUT and reference but also with the mechanical motions. Since this ambiguity exists without any extra hypothesis, it implies that the information from the stitching data acquisition is not enough to solve this ambiguity problem. A possible regularization can be a better knowledge of the reference curvature and twist or the mechanical motion. For example, the real curvature and twist of the reference can be utilized in the stitching algorithm, or the motion angles and straightness can be traced to reduce the uncertainties of the curvature and twist. In this work, our stitching algorithm is regularized to not estimate the piston, x-tilt, y-tilt, x-curvature, twist, and y-curvature for the reference shape.
The estimation of r in Eq. (3) is implemented by using the extended matrices G ex and D ex , as well as the extended vector d ex as where the extensions are d exi = Similarly, the estimation ofr i in Eq. (5) becomeŝ Considering to the large size and the sparsity of the matrices, e.g. G ex and D ex , sparse matrices should be used in the calculation to reduce the memory cost. The function LSMR() [24,25] is recommended to solve large sparse least-squares problems. The estimation ofr i in Eq. (7) can be regularized aŝ where D s ex = The regularization in Eqs. (19)-(21) enables a highly flexible option to add the curvature and twist acquired by another system. The high-order additive systematic errors are estimated with the redundant data.

Simulation-based numerical case study
A stitching measurement of a 40mm × 10mm sample is simulated to verify the proposed iterative self-calibration algorithms. The height distribution of the SUT in Fig. 4(a) is generated as xy p x p y + cos 2π( The SA of the interferometer is 64 × 64 pixels with 0.12 mm lateral resolution. The user-defined overlapping ratio is set as 80% for both directions. As discussed in Section 3.1, randomized step sizes are performed to resolve the periodic errors. The range for the randomization is ±3 pixels. In total, 3 × 23 subsets of interferometric height maps are simulated with the reference r shown in Fig. 4(b), 1 µrad RMS angular adjustment in geometric parameters t, and 0.5 nm RMS additive normally distributed random noises n a . The CS-P, CS-S, and SC algorithms are performed to demonstrate the performance of the proposed 2D self-calibration stitching.
The non-iterative CS-P result is shown in Fig. 5. The stitching result in Fig. 5(a) and estimated reference in Fig. 5(b) are very close to the simulated SUT and reference in Fig. 4. The up to the 2nd order removed stitching error in Fig. 5(c) is about 0.26 nm RMS, and corresponding reference estimation error is 0.07nm RMS as shown in Fig. 5(d). The results using CS-S algorithm are shown in Fig. 6. The first estimated reference r 1 in Fig. 6(b) is similar to the simulated truth with only 0.11 nm RMS in Fig. 6(c). However, low-frequency variations are shown on the stitching error in Fig. 6(d). With the compensation r 2 = r 1 +r 1 in the first iteration, the reference estimation error is 0.07 nm RMS in Fig. 6(g) and the stitching error is only 0.21 nm RMS in Fig. 6(h). Since the RMS value of the next compensation amountr 2 in Fig. 6(i) is less than the preset threshold thr = 0.01 nm, the iteration terminates with the final stitching error shown in Fig. 6(l). In fact, the stitching calculations in iterations shown in Figs. 6(d) and 6(h) are not necessary in CS-S algorithm, and only the reference is estimated through iterative compensation.
By applying the SC algorithm, the reference is estimated during the iterative stitching process in Fig. 7. As shown in Figs. 7(a), 7(e), and 7(i), the RMS value of the up to the 2nd order The SA of the interferometer is 64 × 64 pixels with 0.12 mm lateral resolution. The user-defined overlapping ratio is set as 80% for both directions. As discussed in Section 1, randomized step sizes are performed to resolve the periodic errors. The range for the randomization is ±3 pixels. In total, 3 × 23 subsets of interferometric height maps are simulated with the reference r shown in Fig. 4(b), 1 µrad RMS angular adjustment in geometric parameters t, and 0.5 nm RMS additive normally distributed random noises n a . The CS-P, CS-S, and SC algorithms are performed to demonstrate the performance of the proposed 2D self-calibration stitching.
The non-iterative CS-P result is shown in Fig. 5. The stitching result in Fig. 5(a) and estimated reference in Fig. 5(b) are very close to the simulated SUT and reference in Fig. 4. The up to the 2nd order removed stitching error in Fig. 5(c) is about 0.26 nm RMS, and corresponding reference estimation error is 0.07nm RMS as shown in Fig. 5(d).
where the amplitude A = 200 × 10 −9 m, and the periods p x = p y = 16 × 10 −3 m. The SA of the interferometer is 64 × 64 pixels with 0.12 mm lateral resolution. The user-defined overlapping ratio is set as 80% for both directions. As discussed in Section 3.1, randomized step sizes are performed to resolve the periodic errors. The range for the randomization is ±3 pixels. In total, 3 × 23 subsets of interferometric height maps are simulated with the reference r shown in Fig. 4(b), 1 µrad RMS angular adjustment in geometric parameters t, and 0.5 nm RMS additive normally distributed random noises n a . The CS-P, CS-S, and SC algorithms are performed to demonstrate the performance of the proposed 2D self-calibration stitching.
The non-iterative CS-P result is shown in Fig. 5. The stitching result in Fig. 5(a) and estimated reference in Fig. 5(b) are very close to the simulated SUT and reference in Fig. 4. The up to the 2nd order removed stitching error in Fig. 5(c) is about 0.26 nm RMS, and corresponding reference estimation error is 0.07nm RMS as shown in Fig. 5(d). The results using CS-S algorithm are shown in Fig. 6. The first estimated reference r 1 in Fig. 6(b) is similar to the simulated truth with only 0.11 nm RMS in Fig. 6(c). However, low-frequency variations are shown on the stitching error in Fig. 6(d). With the compensation r 2 = r 1 +r 1 in the first iteration, the reference estimation error is 0.07 nm RMS in Fig. 6(g) and the stitching error is only 0.21 nm RMS in Fig. 6(h). Since the RMS value of the next compensation amountr 2 in Fig. 6(i) is less than the preset threshold thr = 0.01 nm, the iteration terminates with the final stitching error shown in Fig. 6(l). In fact, the stitching calculations in iterations shown in Figs. 6(d) and 6(h) are not necessary in CS-S algorithm, and only the reference is estimated through iterative compensation.
By applying the SC algorithm, the reference is estimated during the iterative stitching process in Fig. 7. As shown in Figs. 7(a), 7(e), and 7(i), the RMS value of the up to the 2nd order The results using CS-S algorithm are shown in Fig. 6. The first estimated reference r 1 in Fig. 6(b) is similar to the simulated truth with only 0.11 nm RMS in Fig. 6(c). However, low-frequency variations are shown on the stitching error in Fig. 6(d). With the compensation r 2 = r 1 +r 1 in the first iteration, the reference estimation error is 0.07 nm RMS in Fig. 6(g) and the stitching error is only 0.21 nm RMS in Fig. 6(h). Since the RMS value of the next compensation amountr 2 in Fig. 6(i) is less than the preset threshold thr = 0.01 nm, the iteration terminates with the final stitching error shown in Fig. 6(l). In fact, the stitching calculations in iterations shown in Figs. 6(d) and 6(h) are not necessary in CS-S algorithm, and only the reference is estimated through iterative compensation.
By applying the SC algorithm, the reference is estimated during the iterative stitching process in  removed stitching error is getting lower, as the reference r is updated during the stitching process. The stitching error with repeating patterns shown in Fig. 7(a) is very typical if the reference is not well calibrated. The reference compensation amountr 0 in Fig. 7(b) is estimated from the stitched result z 0 in this step and the raw SA measures. The updated reference r 1 = r 0 +r 0 in Fig. 7(c) has a reference estimation errors (terms up to the 2nd order removed) only at the sub-nanometer level with a period correlated to the average step size. In the 1st iteration, the stitching errors are reduced from 3.14 nm RMS down to 0.29 nm RMS, which also indicates that the previous reference estimation are effective. The remaining errors are small but still evident with a periodic shape shown in Fig. 7(e). The estimatedr 1 in Fig. 7(f) is very similar to the previous reference estimation error in Fig. 7(d). After updating the reference r 2 = r 1 +r 1 , the high-order discrepancy from the true value is almost random noise as shown in Fig. 7(h). A small random noise like stitching error(terms up to the 2nd order removed) can be expected in the next iteration as demonstrated in Fig. 7(i) with 0.2 nm RMS. The stitching algorithm converges after the 2nd iteration, as the reference compensation amountr 2 in Fig. 7(j) is less than the preset RMS threshold thr = 0.01 nm. The simulation demonstrates that the proposed method is effective in estimating the reference r and stitching the SUT shape z. Next, we will show a real stitching interferometry measurement to verify the feasibility of the proposed method.

Experimental verification
Some experiments are implemented on the interferometric stitching platform developed at NSLS-II [23] to test and verify the proposed method. The setup of the stitching interferometer (see Fig. 8) is composed of a Fizeau interferometer and several translation and rotation stages to adjust the relative position and angle between the SUT and the interferometer reference surface. As marked in Fig. 8, we only need the xand y-motions, θ x and θ y rotations to perform 2D scans and to null the fringes at each subset location in this experiment. The SUT is a 190-mm-long flat silicon mirror. The central 180 mm × 10 mm area is the region of interest. We use a 256 × 64 pixels rectangular mask window as the interferometer SA, and the pixel lateral resolution is 0.12 mm. The user-defined overlap ratios are set as 80% for both xand y-directions. The software algorithm automatically calculates the random step variation to cover the whole region of interest. It ensures that the random step size variations are within the preset limits ([−2, +2] pixels) and the sizes are in integer pixel to avoid the sub-pixel matching issue. As a result, the actual overlap ratio is about 81% in x-direction and 85% in y-direction. In total, 3 × 26 subsets are captured to fully cover the region of interest. We take 64-averaged measurement at each SA. The repeatability of a 64-averaged measurement in our environment is about 0.2 nm RMS. The tip-tilt threshold is set as 1 µrad to null fringes before capturing each subset. The RMS threshold for reference compensation thr = 0.01 nm.
In order to make a quantitative evaluation of the proposed self-calibration, a well-polished 240mm×40mm flat X-ray mirror with about 0.6 nm RMS is used to calibrate the reference by averaging several uncorrelated regions. This calibrated reference ( Fig. 9(a)) is used as the benchmark to judge the estimated reference during iterations. The stitched SUT shape in Fig. 9(b) is used as the benchmark to evaluate the stitching errors.

Experimental verification
Some experiments are implemented on the interferometric stitching platform developed at NSLS-II [23] to test and verify the proposed method. The setup of the stitching interferometer (see Fig. 8) is composed of a Fizeau interferometer and several translation and rotation stages to adjust the relative position and angle between the SUT and the interferometer reference surface. As marked in Fig. 8, we only need the x-and y-motions, θ x and θ y rotations to perform 2D scans and to null the fringes at each subset location in this experiment. The SUT is a 190-mm-long flat silicon mirror. The central 180 mm × 10 mm area is the region of interest. We use a 256 × 64 pixels rectangular mask window as the interferometer SA, and the pixel lateral resolution is 0.12 mm. The user-defined overlap ratios are set as 80% for both xand y-directions. The software algorithm automatically calculates the random step variation to cover the whole region of interest. It ensures that the random step size variations are within the preset limits ([−2, +2] pixels) and the sizes are in integer pixel to avoid the sub-pixel matching issue. As a result, the actual overlap ratio is about 81% in x-direction and 85% in y-direction. In total, 3 × 26 subsets are captured to fully cover the region of interest. We take 64-averaged measurement at each SA. The repeatability of a 64-averaged measurement in our environment is about 0.2 nm RMS. The tip-tilt threshold is set as 1 µrad to null fringes before capturing each subset. The RMS threshold for reference compensation thr = 0.01 nm.
(b) (a) Fig. 9. By subtracting a well-calibrated reference (a) from the captured subsets, the stitching result of the SUT (b) will be used as the benchmark for the stitching error evaluation.
In order to make a quantitative evaluation of the proposed self-calibration, a well-polished 240mm×40mm flat X-ray mirror with about 0.6 nm RMS is used to calibrate the reference by averaging several uncorrelated regions. This calibrated reference ( Fig. 9(a)) is used as the benchmark to judge the estimated reference during iterations. The stitched SUT shape in Fig.  9(b) is used as the benchmark to evaluate the stitching errors.

Self-calibration stitching results
We present self-calibration stitching results with CS-P and SC algorithms. The CS-P estimated reference is shown in Fig. 10(a) with a 0.11 nm RMS high-order reference estimation error in Fig. 10(b). With the reference r, the stitched height z is shown in Fig. 10 Fig. 9. By subtracting a well-calibrated reference (a) from the captured subsets, the stitching result of the SUT (b) will be used as the benchmark for the stitching error evaluation.

Self-calibration stitching results
We present self-calibration stitching results with CS-P and SC algorithms. The CS-P estimated reference is shown in Fig. 10(a) with a 0.11 nm RMS high-order reference estimation error in Fig. 10(b). With the reference r, the stitched height z is shown in Fig. 10(c). The iterative stitching and reference estimation results with the SC algorithm are shown in Fig. 11. The calculation starts with r 0 = 0. The stitched z 0 shows repeating error patterns due to the lack of knowledge on the reference. It is obvious that their period is about the averaged step size.
The reference compensation amountr 0 is estimated from the difference between m and z 0 . The main topography of the reference r 1 = r 0 +r 0 is updated via this estimation. However, comparing to the calibrated reference, there are periodic errors shown in Fig. 11(d), which results in obvious waviness errors on the stitched height map in the next iteration shown in Fig. 11(e). In SC iteration 1, The reference compensation amountr 1 is estimated as shown in Fig. 11(f). After this compensation, the discrepancy between the updated reference r 2 = r 1 +r 1 and the calibrated reference benchmark is very small (only a 0.11 nm RMS) with no obvious periodic error shown in Fig. 11(h). The periodic patterns on the stitched z 2 is reduced (see Fig. 11(i)) in benchmark to judge the estimated reference during iterations. The stitched SUT shape in Fig.  9(b) is used as the benchmark to evaluate the stitching errors.

Self-calibration stitching results
We present self-calibration stitching results with CS-P and SC algorithms. The CS-P estimated reference is shown in Fig. 10(a) with a 0.11 nm RMS high-order reference estimation error in Fig. 10(b). With the reference r, the stitched height z is shown in Fig. 10(c). The iterative stitching and reference estimation results with the SC algorithm are shown in Fig.  11. The calculation starts with r 0 = 0. The stitched z 0 shows repeating error patterns due to the lack of knowledge on the reference. It is obvious that their period is about the averaged step size. the SC iteration 2. Finally, the calculatedr 2 in Fig. 11(j) is less than the preset RMS threshold thr = 0.01 nm. The iteration ends up with the final stitched SUT height z 2 and the estimated reference r 2 . The final results of the SC algorithm in Fig. 11 is very similar to the result of the CS-P algorithm in Fig. 10. The following studies are based on results of the SC algorithm.

Repeatability study
To study the repeatability of the proposed self-calibration stitching method, 10 repeating scans are performed. These 10 SC algorithm stitched SUT maps are shown in Fig. 12(a). To highlight the details of the height variation on the SUT, the color range is only within [−3σ, +3σ] where σ is the RMS value of the stitched height values.
Illustrated with the same color map, Fig. 12(b) shows their discrepancies from the average map. All the discrepancies give an overall repeatability value less than 0.1 nm RMS.
The 10 estimated reference maps are displayed in Fig. 13(a) with only 0.015 nm RMS overall repeatability in their discrepancies shown in Fig. 13(b), which illustrates excellent repeatability of the reference estimation.

Comparison with the stitching result using the calibrated reference
The stitched SUT height with the proposed SC algorithm (Fig. 11(i)) is evaluated by comparing with the stitching result using the calibrated reference as benchmark (Fig. 9). The stitching result without subtracting the systematic error ( Fig. 11(a)) is also evaluated as a contrast. The comparison results are shown in Fig. 14.
The obvious periodic error in Fig. 14(a) is one of the drawbacks to use a stitching mechanism with systematic errors. Our proposed self-calibration method provide a way to in situ estimate the repeatable high-order additive systematic error during the analysis of the measurement data. As shown in Fig. 14(b), the artificial waviness diminishes in our stitched results without any pre-knowledge of the systematic error.

Self-consistency study using different SA windows
We perform another measurement of the same SUT using a new SA with the same size (256 × 64 pixels) but in a different location on our 4-inch pupil interferometer. The previous SA is labelled as SA 1 and the new one as SA 2 .
Comparing to SA 1 , the SA 2 window is shifted 128 pixels to the right as shown in Fig. 15. This shift will allow us to check if the self-calibration stitching method will give self-consistent reference results in their common area SA 1 ∩ SA 2 .  Fig. 11. With iterative compensations to the reference r in the SC algorithm, the SUT height z is stitched with better visualization (less periodic errors). Fig. 11. With iterative compensations to the reference r in the SC algorithm, the SUT height z is stitched with better visualization (less periodic errors).
reference r 2 . The final results of the SC algorithm in Fig. 11 is very similar to the result of the CS-P algorithm in Fig. 10. The following studies are based on results of the SC algorithm.

Repeatability study
To study the repeatability of the proposed self-calibration stitching method, 10 repeating scans are performed. These 10 SC algorithm stitched SUT maps are shown in Fig. 12(a). To highlight the details of the height variation on the SUT, the color range is only within [−3σ, +3σ] where σ is the RMS value of the stitched height values. Illustrated with the same color map, Fig. 12(b) shows their discrepancies from the average map. All the discrepancies give an overall repeatability value less than 0.1 nm RMS.

Repeatability study
To study the repeatability of the proposed self-calibration stitching method, 10 repeating scans are performed. These 10 SC algorithm stitched SUT maps are shown in Fig. 12(a). To highlight the details of the height variation on the SUT, the color range is only within [−3σ, +3σ] where σ is the RMS value of the stitched height values. Illustrated with the same color map, Fig. 12(b) shows their discrepancies from the average map. All the discrepancies give an overall repeatability value less than 0.1 nm RMS.   Fig. 13. The estimated reference (a) and the discrepancies (b) from their average of these 10 repeating scans indicate the reference estimation is very repeatable.
The 10 estimated reference maps are displayed in Fig. 13(a) with only 0.015 nm RMS overall repeatability in their discrepancies shown in Fig. 13(b), which illustrates excellent repeatability of the reference estimation.

Comparison with the stitching result using the calibrated reference
The stitched SUT height with the proposed SC algorithm (Fig. 11(i)) is evaluated by comparing with the stitching result using the calibrated reference as benchmark (Fig. 9). The stitching result without subtracting the systematic error ( Fig. 11(a)) is also evaluated as a contrast. The comparison results are shown in Fig. 14. with the proposed self-calibration stitching method indicates its effectiveness. These two stitching error maps are plotted with surfaces in the same coordinate system with a relative shift in y-direction for better comparison.
The obvious periodic error in Fig. 14(a) is one of the drawbacks to use a stitching mechanism with systematic errors. Our proposed self-calibration method provide a way to in situ estimate the repeatable high-order additive systematic error during the analysis of the measurement data. As shown in Fig. 14(b), the artificial waviness diminishes in our stitched results without any pre-knowledge of the systematic error.
5.4. Self-consistency study using different SA windows with the proposed self-calibration stitching method indicates its effectiveness. These two stitching error maps are plotted with surfaces in the same coordinate system with a relative shift in y-direction for better comparison.
As shown in Fig. 14(b), the artificial waviness diminishes in our stitched results without any pre-knowledge of the systematic error.

Self-consistency study using different SA windows
We perform another measurement of the same SUT using a new SA with the same size (256 × 64 pixels) but in a different location on our 4-inch pupil interferometer. The previous SA is labelled as SA 1 Fig. 15. Experiments with a 50% overlapping area in common between SA 1 and SA 2 are designed to study the self-consistency of the 2D self-calibration stitching method. Fig. 15. Experiments with a 50% overlapping area in common between SA 1 and SA 2 are designed to study the self-consistency of the 2D self-calibration stitching method.
From Fig. 16, it is not difficult to notice that the stitched SUT shapes using different SAs are very close to each other with only 0.12 nm RMS height difference on the 180 mm × 10 mm surface area. On the other hand, the estimated reference maps in Fig. 17 also indicate the self-consistency of the stitching method as the reference maps in the common area of the two SAs are similar to each other with only 0.10 nm RMS difference, and this difference is dominated by the second order terms.
Comparing to SA 1 , the SA 2 window is shifted 128 pixels to the right as shown in Fig. 15. This shift will allow us to check if the self-calibration stitching method will give self-consistent reference results in their common area SA 1 ∩ SA 2 . From Fig. 16, it is not difficult to notice that the stitched SUT shapes using different SAs are very close to each other with only 0.12 nm RMS height difference on the 180 mm × 10 mm surface area. On the other hand, the estimated reference maps in Fig. 17 also indicate the self-consistency of the stitching method as the reference maps in the common area of the two SAs are similar to each other with only 0.10 nm RMS difference, and this difference is dominated by the second order terms.  By removing the 2nd order terms, the difference inside the common area drops to 0.03 nm RMS, which shows excellent self-consistency of the proposed method. Since the terms up to the 2nd order are not involved in the estimation as mentioned in Eq. (15), the small curvature difference in the common area between the two estimated reference maps is due to the different height variations in the non-overlapping regions of these two SAs which result in different Comparing to SA 1 , the SA 2 window is shifted 128 pixels to the right as shown in Fig. 15. This shift will allow us to check if the self-calibration stitching method will give self-consistent reference results in their common area SA 1 ∩ SA 2 . From Fig. 16, it is not difficult to notice that the stitched SUT shapes using different SAs are very close to each other with only 0.12 nm RMS height difference on the 180 mm × 10 mm surface area. On the other hand, the estimated reference maps in Fig. 17 also indicate the self-consistency of the stitching method as the reference maps in the common area of the two SAs are similar to each other with only 0.10 nm RMS difference, and this difference is dominated by the second order terms.  Fig. 17. The reference estimations in these two 50% overlapped SAs (a) give self-consistent results in their common area (b) and (c). There is only a tiny height difference in their common area with the tilt terms removed (b) or up to the 2nd order terms removed (c). difference in the common area between the two estimated reference maps is due to the different height variations in the non-overlapping regions of these two SAs which result in different curvature values in each SA.

Discussion
We have proposed two types of 2D self-calibration stitching algorithms (the CS and SC algorithms) for stitching interferometry. In this work, the calculation is fully based on pixel-to-pixel relations, since the stitching steps are exact integer pixel numbers. However, the proposed algorithms can be easily extended and applied to an interpolation representation of the reference surface with a linear combination of basis functions when the stitching steps are not integer pixels and the interpolation is unavoidable. In that case, the relation matrices in the CS-P and CS-S algorithms are no longer that sparse which could yield heavier computation issues.
The merits and limitations of the proposed CS-P, CS-S, and SC algorithms are discussed below.
In pixel-relation-based stitching strategy, the CS-P algorithm directly estimate the reference and stitched height at the same time. Although it can provide a direct least squares solution with LSMR() function for data with a reasonable size shown in the experiment, it needs to pay careful attentions to the tolerances in LSMR() and the weights of the extended equations to get accurate stitching and estimation results.
Different from the CS-P algorithm, the CS-S and SC algorithms separate the reference estimation and height stitching, which reduces the computational dimension of the problem. The CS-S algorithm estimates the reference first and then uses it for shape stitching. In CS-S algorithm flow, the reference estimation has less influence from the stitched SUT result, and it has a better efficiency in utilizing the estimated reference comparing to the SC algorithm. The merit of SC algorithm is its capability of keeping the size of the relation matrix D s in Eq. (6) fixed as M × M, which does not depend on the subset number N. Therefore, it has a unique advantage in dealing with stitching with large number of subsets. It is worth to note that the computation in assembling the matrix D s can be time consuming for a limited computing resources, if M becomes a huge number.

Conclusion
In our work, the repeatable high-order additive systematic errors can be estimated from the redundant data acquired in xand y-directions stepping with big overlaps of neighboring subsets. The periodic errors due to the uniform steps in data acquisition are carefully studied. Randomized nonuniform step sizes are suggested to "break" the artificial periodicity during the data acquisition. Regularization on curvature and twist are introduced in the algorithms to potentially preserve the measured reference curvature and twist. The proposed algorithms have been cross-checked and verified using both simulation and experimental data. Further study on the memory cost and speeding up the computation will be our future work to resolve this practical issue.