1. Introduction

The sub-aperture stitching procedure is very useful to overcome the limited numerical aperture (NA) or to extend vertical resolution of conventional interferometers in optical metrology. In fact, several approaches have been conceived at the aiming to increase the NA of interferometric system and thus in order to enhance spatial resolution or to extend the field of view. Increasing the resolution was achieved by rotating the sample in digital holography in order to record the diffraction wave-front by diverse directions [1]. Diffraction limit was overcome by Kuznetsova et al. [2] by rotating the test object with respect to the optical axis for collecting the rays scattered at wider angles inside the limited aperture of the optical system. Mico et al. demonstrated a scheme for augmenting the resolution of the imaging systems based on employing tilted illumination and common-path interferometry [3,4]. F. Pan et al. achieved the resolution enhanced and speckle suppressed by numerical incoherent superposition in Fresnel off-axis digital holography with different titled illuminations [5]. Instead, some methods were presented which are based on the shift of the CCD to different positions in order to improve the system’s resolution [6]. Martínez-León and Javidi displaced the CCD camera in order to increase both the spatial resolution and the sampling in the recoding procedure [6]. Finally, other resolution enhancement techniques have also been proposed. LeClerc et al. [7] increased synthetically the NA by an on-axis digital heterodyne holography. Besides, spatial resolution was gained by simply using a 1D diffraction grating [8] or a 2D tuneable phase grating [9,10] that allowing to collect light scattered at larger angles. In all the aforementioned approaches usually, it is requested to find a smart procedure to stitching the spectra in the Fourier domain, i.e. the spatial frequency domain, of each recorded hologram in order to get a single synthetic digital hologram with augmented NA. In other cases like in a work by Gyímesi et al. [11] or by Pelagotti et al. [12] the synthetic aperture of the digital hologram is achieved by directly stitching the acquired holograms and matching directly the fringes of each sub-aperture hologram.

On the other side, the so named sub-aperture stitching interferometry (SSI) is proposed for enhancing the spatial resolution or for extending the field of view of conventional interferometers [13]. SSI has been applied in large aperture and large high numerical aperture surface metrology, including planar [14], spherical [15,16], aspheric [17,18] and cylindrical surfaces [19,20]. Usually the full aperture of tested surface is divided into a series of smaller sub-apertures [21]. Then, all of sub-apertures are measured by moving the tested surface. Finally, the full aperture map is obtained by joining these sub-aperture maps. Measurement accuracy is often affected by the alignment errors between the tested object and the rotational axis or between the measuring probe and the rotational axis. In order to minimize alignment errors, a precise platform can be used to deal with these alignment errors. In commercial and very complex interferometers (i.e. QED Technologies) sub-aperture stitching interferometry (SSI) works well by a fully integrated and automated adjustor with multiple degrees of freedom to make each sub-aperture null adjustment thus providing the accurate position and pose of each recorded sub-aperture data [22,23]. A step forward was proposed by Pengfei Zhang et al. by integrating a binocular stereovision system into the SSI, that allowing to obtain accurate position and pose information of the tested surface by measuring the coordinates of intentional marks attached to the fixture [24]. Instead, Xiaokun Wang et al. marked some points on the tested surface to accomplish the alignments between sub-apertures and calibrated the relationship between the coordinates of the tested surface and the camera pixels [25]. On the other hand, it is important to note that the sub-aperture stitching algorithms are also used to compensate the alignment errors. So far, there are many landmark algorithms have been proposed for sub-aperture stitching, including Kwon Thunen method [26] and the simultaneous fit method [27], the discrete phase method [28], the multi-aperture overlap-scanning technique [13], the sub-aperture stitching and localization (SASL) algorithm [29]. These sub-aperture stitching algorithms primarily boil down to solve the least-squares problem, thereby minimizing the relative alignment errors between the adjacent sub-apertures and realizing the stitching treatment. The above algorithms differ from each other by the method for modelling the alignment errors, such as the slope-based [30], Zernike-polynomials-based, analytical model-based [31], ray tracing-based [32], and configuration space-based [33] method. However, the systematic errors introduced by the reference surface error and by the measurement system aberration can affect measurement accuracy. Those errors may accumulate in the sub-aperture measurements and can be amplified during the stitching process. One method to reduce the systematic errors is to provide a direct measurement of the reference surface error by an absolute test and then remove error from each sub-aperture measurement [34]. An alternative method can be removing the reference surface error in the sub-aperture stitching process on the basis of a special algorithm [35]. Lei Huang et al. proposed two-dimensional self-calibration stitching algorithm with liner least squares method, which calibrated the high-order additive system errors by minimizing the de-tilted discrepancies of the overlapped subsets [36,37]. Looking at the existing sub-aperture stitching algorithms, it mainly comes down to solving the least squares problem. The solutions are mainly iterative optimization, reverse reconstruction and so on. And this article is to solve the least squares problem through a new idea. In our previous works, we proposed a sub-aperture stitching interferometry based on the off-axis digital holography [38] and an optimizing stitching strategy, which can simultaneously compensate the alignment errors and the system aberration in the stitching process [39].

In this paper, we present a novel sub-apertures stitching approach in DH to extend the field of view of interferometry based on machine learning. According to the computation model of sub-aperture stitching, we construct a network for correcting sub-aperture alignment errors in rotation measurements and systemic aberration of optical path. We use a popular open-source machine learning library, TensorFlow [40], for setting up the network and conducting the network training. By iteration of forward pass and backward pass, the optimal weights of network are generated when loss function minimizes. Then, the aberrations caused by alignment errors and systemic aberration in the sub-aperture maps are abstracted by Zernike polynomial fitting based on these optimal coefficients. Finally, a high-precision full aperture map is recovered through fusing correct sub-aperture maps. Furthermore, we analyze the performance of different learning rates, epochs, and batch sizes in the learning process. The computational model of sub-aperture stitching and corresponding network are introduced in Section 2. Furthermore, a hemisphere optical surface is tested to verify the proposed method, and we analyze the performance of different learning rates, epochs, and batch sizes for getting optimal model, as described in Section 3. Section 4 summarizes the conclusions.

2. Computational model of sub-aperture stitching

In the digital holography sub-aperture stitching interferometry, the measurement process contains two steps. First, the holograms of sub-apertures are sequentially recorded. Then the phase maps of sub-apertures are retrieved by digital holographic reconstruction. Second, an appropriate stitching strategy is employed to synthesize all sub-aperture maps and remove the undesired errors. Thus, the full aperture map is obtained. Figure 1 is the schematic diagram of measuring and stitching of sub-aperture.

 

Fig. 1. Sub-aperture measuring and stitching.

Download Full Size | PDF

In the hologram recording process, the center of curvature of the measured surface is placed at the focus of the illuminating spherical wave, as shown in Fig. 1(a). Thus, all sub-apertures are nulled and tested by rotating the tested surface. After that, the sub-aperture phase maps are retrieved by digital holographic reconstruction. In this process, there are 4 steps: Hologram Apodization, Hologram Filter, Centered Spectrum, and Diffraction Calculation. Specially, Hologram Apodization is to avoid diffraction ripples by multiplying the hologram with a window function. Hologram Filter is to filter out twin image term and zero-order image term in the frequency domain by using a spatial filter. Centered Spectrum is used to compensate the title aberrations due to the off-axis hologram recording. Diffraction Calculation is implemented by using angular spectrum algorithm. Subsequently, based on the position information in the rotation measurement, the sub-aperture phase maps in local coordinate frame (x_l, y_l, z_l) are synthesized to a full aperture map in general coordinate frame (x_g, y_g, z_g), as shown in Fig. 1(b). It is noted that although a high-precision stage is used in the experiment, the sub-aperture maps inevitably suffer from residual position errors coming from uncertainties of manual operations and mechanical structure. As mentioned in reference [38], these misalignments would cause phase aberrations in the sub-aperture maps, including piston aberration, the tilt aberrations, and the defocus aberration. Moreover, the system aberrations are also embedded in the sub-aperture measurements, including the undesired phase curvature and high-order aberrations induced by lens and other components in optical system. In this paper, we employ the machine learning to remove the residual misalignments and the systemic aberration in the stitching process.

The schematic diagram of the computation model is shown in Fig. 2. The measured phase map of a sub-aperture can be divided into two parts and described as

 

Fig. 2. The computation model of sub-aperture stitching.

Download Full Size | PDF

Further, according to the sub-aperture measurement process, the phase map of P_{sub_fitting} can also be described as

Finally, when these phase maps in Eq. (2) are fitted by using the Zernike polynomials, the Eq. (2) is rewritten as

If a set of appropriate coefficients in Eq. (3) are obtained from the sub-aperture measurements, the full aperture map can be retrieved by Zernike polynomial fitting. Here, we perform this optimal process by using machine learning. First, a network is built according to the above computational model, as shown in Fig. 3. The network consists of an input layer, a multiple layer, an addition layer, and an output layer. The input layer contains 4 parts, including the phase value of each point in sub-apertures, P_{sub_measured}(x_g,y_g), 36-term Zernike polynomials with corresponding global coordinates of each point, Zernk_i(x_g,y_g), i=1,2,…,36, 4-term Zernike polynomials with global coordinates of each point, Zernk_i(x_g,y_g), i=1,2,…,4, and 36-term Zernike polynomials with local coordinates of each point, Zernk_i(x_l,y_l), i=1,2,…,36. The multiplication layer performs the multiplications of these 3 sets of Zernike and 3 set of weights of network, such as Weight_i^ts, Weight_i^ae, and Weight_i^se. There are 3 channels to fit the phase maps of the tested surface, alignment errors, and system errors, respectively. An addition layer is used to add the 3 terms of multiple layer. The loss layer is to calculate the difference between the current predictions of P_{sub_fitting}, and the actual measurements of P_{sub_measured}. The output layer outputs the minimum value of loss layer. It is noted that weights of network are actually the unknown polynomial coefficients and can be learned by network training.

 

Fig. 3. The network for learning the optimal coefficients.

Download Full Size | PDF

In the training process, the initial weights of network are set to a random value between 0 and 1. A forward pass is implemented, in which the input data travels through all layers and generates output values of loss function. Then, a backward pass is performed to update the learnable weights according to the difference between the predictions and measurements. So, the network can be trained using steepest descent exactly as the backpropagation algorithm in neural networks [43]. Specifically, the learning algorithm carries out the minimization of loss function:

The optimization of weights is carried out iteratively by calculating the gradient, which is part of the backward pass of learning approach, with respect to loss function. Especially, Adam optimizer is used to perform this recursive computation of the gradient based on the value of the loss function considering the efficiency and adaptability of the algorithm [43]. Then, the weights of the network are updated according to the first and second moment estimates of gradient. Taking into account of the speed and convergence of iterative optimization, we assess two ways of calculating the loss function and updating the weights of the network. One is to calculate the loss function with all sub-aperture measurements and network predictions and update the network weights based on the average loss value. Another is to calculate the loss function and update network weights one by one, until all the sub-aperture measurements and network predictions have been involved. Furthermore, we investigate an appropriate batch size to achieve high-efficiency of learning process. In principle, the learning rate and the epoch number determine the accuracy and efficiency of machine learning. Therefore, we investigate the significance of the two parameters in the network training.

After network training, 3 sets of optimal weights are obtained. Based on these coefficients, the real profile of whole tested surface can be retrieved. The synthesis procedure is shown in the Fig. 4.

 

Fig. 4. The complete full aperture phase map.

Download Full Size | PDF

First, the full aperture map of P_{sub_fitting} can be recovered by Zernike polynomial fitting. The full aperture profile of P_{full_profile} is also retrieved, which describes the real shape of the whole tested surface with the low-frequency information. The sup-aperture residuals of P_{sub_residuals} can be obtained by subtracting P_{sub_fitting} from P_{sub_measured} for each sub-aperture. Further, by synthesizing all P_{sub_residuals}, the full aperture residuals of P_{full_residuals} are retrieved, which presents the high-frequency information of the whole tested surface, as shown in Fig. 4. Finally, the complete full-aperture map of P_{full_result} is expressed as

3. Experimental verification and result analysis

In order to verify the proposed approach, a convex lens/hemisphere with a diameter of 25.4 mm and a radius curvature of 13.1 mm is measured by digital holographic sub-aperture stitching interferometry. The schematic diagram of the optical system is shown in Fig. 5. A linear He–Ne laser (632.8nm, 5mW) is used as the light source. The combination of the polarizing beam splitter (PBS) and the half-wave plate (λ/2) is used for the adjustment of the intensity ratio of the beams in the reference arm and the object arm. In the reference arm a half-wave plate (λ/2) is introduced to obtain parallel polarization states in both arms at the exit of the interferometer. The beam expanders (BE1, BE2), including pinholes for spatial filtering, are introduced to produce plane waves. For reflection testing, a custom-made complex lens (CL) with high NA is inserted between the BS and the tested sphere. The reflected wave-front, called the object wave, provides information about the figure and irregularity of the surface under test. The object wave O and the reference wave R are recombined by the BS and interfere at the exit of the interferometer, where a black and white CMOS camera (the resolution of 2048×2048 pixels, the pixel size of 5.5×5.5 µm) records the hologram intensity. In order to create off-axis holograms, the orientation of the mirror (M4) which reflects the reference wave is adjusted so that the reference wave reaches the CMOS with a tilt angle, while the object wave propagates perpendicularly to the hologram plane. In this experiment, we arrange 24 off-axis sub-apertures on 3 rings around the central one for covering the full aperture of the tested surface and providing sufficient overlaps, as illustrated in Fig. 6(a-1). In the training process, the Zernike fitting is recursively implemented based on the sub-aperture measurements, especially in the global coordinates. Zernike polynomials are a complete set of polynomials that are orthogonal or orthonormal in a continuous fashion over the interior of a unit circle. To ensure orthogonality and accuracy, it is necessary that enough sampling points are involved in fitting discrete data. Considering the measurement efficiency and calculation time, we choose about 40% overlap between the adjacent sub-apertures. The sub-apertures scheme along the meridians and parallels distributed uniformly on the spherical surface. The off-axis angles are 15°, 30°, and 55° respectively. Because of the big off-axis angle at the edge, the projection of the outline of the outmost sub-aperture on OXY plane is no longer a circle. For the hemisphere, the projections of the sub-apertures that are not adjacent in the OXY plane will overlap. So, we expand the spherical sub-apertures along the meridian. First, the spherical peak is made as the origin, and then the other points are distributed along the meridian direction; the expanded sub-aperture distribution is shown in Fig. 6(a-2). In the measuring process, the digital holograms of different sub-apertures are recorded individually by rotating a stage with multiple degrees of freedom against the symmetrical axis of the lens. Then, the measured map of each sub-aperture was numerically retrieved by using digital holographic reconstruction algorithm [44]. After that, we applied the proposed approach to recover the full aperture map of the lens. For verification of the proposed approach, the lens is measured also by a commercial interferometer (i.e. QED sub-aperture stitching interferometer).

 

Fig. 5. The schematic diagram of the optical system.

Download Full Size | PDF

 

Fig. 6. The training results with different learning rate for 24 sub-apertures. (a-1) is the spatial distribution of 24 sub-apertures and (a-2) is an expanded sub-aperture distribution on the plane, (c) and (d) are partial enlargements of the (b), respectively.

Download Full Size | PDF

For the training model, there are several crucial parameters related to the convergence of training, which are learning rate, batch size, and epoch number. Here, batch size corresponds to the number of sub-aperture maps in one forward/backward pass or one iteration. Therefore, the number of iterations in one epoch is equal to N-batchSize+1, where N is a total number of sub-aperture and it is equal to 24 in this experiment.

We first discuss the training results with different learning rates. In Fig. 6(b), the loss function curves with 500 epochs and 24 batch sizes are shown with different learning rates. It can be seen that each curve converges under a finite iteration, and the lower learning rate decays the loss slower and it can’t converge to the minimum. On the other hand, although the higher learning rate is able to decay the loss faster, the value of the loss is easy to get stuck and emerge shock, as shown in Figs. 6(c) and 6(d), which are partial enlargements of Fig. 6(b). Based on the above results, we choose the learning rate of 0.1 and start training to get the desired result.

To find a suitable batch size, we investigate the results with different ‘batch Size’ in the training process. Here, the batch size is selected as 1, 6, 12, and 24 respectively. If the batch size equals to 1, it is meant that the loss function is calculated by putting only one sub-aperture into network. If the batch size equals to 24, it is meant that the loss function is calculated with simultaneously putting all sub-apertures into network. It is noted that each batch moves only one sub-aperture for each iteration calculation to make the overlap between adjacent batches as large as possible. Table. 1 lists the results with the four batch sizes and the results of QED. Here, the number of iterations of the four groups are the same. It can be seen that the results obtained from the lower batch size differ greatly from the QED results, and then the peaks and valleys (PV) and RMS values are closest to QED results when the batch size is equal to 12.

Table 1. Measurement results of the four batches and QED.

View Table | View all tables in this article

Figure 7 shows the phase maps and Zernike coefficient bar graphs of different batch sizes, and it can be seen that the results of all 4 different batch sizes are consistent in QED result. In general, the larger the batch size is, the smaller the error is, and the closer training result is the QED result. It also proves the above view that the method with the smaller batch size is faster, but the convergence is not good. Also, the best training efficiency is acquired when the batch size is equal to 12.

 

Fig. 7. The results with 4 different batch sizes. (a-1) the full-aperture map of QED interferometer, (b-1)-(e-1) the full-aperture maps with the batch of 1, 6, 12, and 24, (a-2) Zernike coefficients based on QED results, (b-2)-(e-2) Zernike coefficient based on the network training results with 4 different batch sizes.

Download Full Size | PDF

 

Fig. 8. The curves of the loss function with 4 different batches. (a) - (d) the converge curves with the batch size of 24, 12, 6 and 1, respectively. (a-1)–(d-1) are the partially enlarged view of the corresponding position on the left.

Download Full Size | PDF

Figure 8 shows the loss curves for 4 different batch sizes, where the horizontal axis represents the number of iterations and the vertical axis represents the loss value at a log scale. Figures 8(a) –8(d) correspond to the loss convergence curves of the batch size equal to 24, 12, 6 and 1. Figures 8(a-1), 8(b-1), 8(c-1) and 8(d-1) correspond to the enlarged view in the left dotted line frame, and it can be seen that the curve of graph Fig. 8(a) is relatively smooth, while the graph Figs. 8(b)–8(d) shows oscillation during a different period, which are 26, 38, 48 iterations respectively. The larger the lot size is, the smoother the convergence curve is and the faster the convergence is. However, when the batch size is equal to 24, it is very time consuming and high-speed computer hardware is needed. Therefore, considering the effects of training efficiency and accuracy, in this experiment the parameter selected for training the model in the batch size is 12.

Thirdly, we investigate the results with different epochs. The training results are shown in Fig. 9 at the different epochs from 1 to 500, the graph Figs. 9(a-2)–9(f-2) are section view plot of Figs. 9(a-1)–9(f-1), respectively. It can be seen that the results are indeed approaching towards QED result with the first 300 epochs. After training with more epochs, the results tend to be stable and close to the QED result.

 

Fig. 9. The results in the training process. (a) - (f) the training results with epoch number of 1, 100, 200, 300, 400, and 500.

Download Full Size | PDF

Table 2 shows the PV and RMS of the training results with different epochs, it can be seen that the value of first three groups of data changes greatly, and the values of last three sets tend to be stable, indicating that the results are convergent and accurate.

Table 2. The value of PV and RMS in different epochs.

View Table | View all tables in this article

It can be seen from Fig. 9 and Table 2 that when the batch is larger than 300, the results begin to converge and tend to be close the QED result. When the epoch is 500, the loss value is minimum and the result is similar to the QED result. So we chose 500 batches as the research object. Finally, the result of the proposed method with the epoch number of 500, the batch size of 12, and the learning rate of 0.1 is shown in Fig. 10(a), with PV of 0.438λ and RMS of 0.087λ. The result of QED interferometer is shown in Fig. 10(c), with PV of 0.518λ and RMS of 0.103λ. The cross-sectional curves are shown in Figs. 10(b) and 10(d), respectively. It is found that the result of the proposed approach has a basic consistency with the result of QED.

 

Fig. 10. (a) full-aperture map by proposed method, (c) full-aperture map by QED interferometer, (b) and (d) are section curves of (a) and (c), respectively.

Download Full Size | PDF

Furthermore, we extracted the Zernike polynomial coefficients from the two results and show them in Fig. 11(c) with the side views of the surface errors shown in Fig. 11(a) and Fig. 11(b), respectively. It can be seen that these Zernike polynomial coefficients appear nearly identical. Although there are some subtle differences, it would be mainly attributed to the aberration compensation. In addition, it can be seen from Fig. 11(c) that the coefficients of the 9th, 16th, and 25th terms are larger than other coefficients. These three terms correspond to the third-order spherical aberration, the fifth-order spherical aberration, and the seventh-order spherical aberration, indicating that the measured aberration is mainly characterized by spherical aberration, and the third-order spherical aberration has the largest contribution. The fifth-order and seventh-order spherical aberration have a slightly smaller impact.

 

Fig. 11. (a) and (b) are the oblique plots of the proposed method and QED interferometer, (c) is the Zernike coefficient bar graph of the two methods.

Download Full Size | PDF

4. Conclusion

We model the sub-aperture stitching process by using a network and the full aperture phase map is spliced by training the model, while the sub-aperture alignment errors and the systemic aberration caused by the optical systems are totally compensated. In the network, the input are the sub-aperture phase maps. The multiplication layer is built with Zernike polynomials and coefficients, having three channels for fitting the surface profile of sub-aperture and the aberration related to the alignment errors and the measuring system. These Zernike coefficients are considered as learnable weights in the training process. The output of the network is modeled as the loss function we aim to minimize and the output of output layer are the optimal weights. The full aperture map without the alignment errors and the systemic aberration are obtained by Zernike polynomial fitting based on these optimal coefficients. In order to optimize the training model, parameters that are crucial to the convergence of model training, namely learning rate, batch size, and epoch number, are analyzed. The experimental verification is carried out, in which a convex lens is measured by digital holographic sub-aperture stitching interferometry and a commercial interferometer (i.e. QED interferometer), respectively. In conclusion, we demonstrate that sub-aperture stitching interferometry is accurately achieved by modelling and eliminating the systematic error and the alignment errors in network training. Moreover, it can be seen that the advantage of our method is that it can handle the stitching of a large number of sub-apertures by choosing different batch sizes for experiments to avoid being limited by computer hardware. We compared results obtained by our holographic system with the measurement achieved by an industrial but more complex interferometer. Results from our experimental measurements show that the relative deviations of PV and RMS are 0.08λ and 0.016λ respectively thus demonstrating the feasibility and accuracy of the proposed new approach by a simpler and more flexible measurement system based on digital holography. Moreover, we demonstrate for the first time, to the best of our knowledge the application of machine learning for assisting the sub-apertures stitching processes in holographic and interferometric systems. The proposed method has the potential for addressing rapid measurement of surface quality in realistic workshop conditions with high-precision and low cost when compared to the state of the art instruments.