Digital holography (DH)^[1] provides access to the quantitative measurement of in-depth information of an object. Due to its advantages, DH has been applied widely in various fields, such as microscopic three-dimensional (3D) measurement and imaging, including 3D cell imaging^[2,3], 3D micro-structure measurement^[4], super-resolution measurement^[5], 3D micro-particles tracking^[6,7], etc. Although DH has many advantages compared with conventional holography, its lateral resolution and field of view (FOV) are its major deficiencies due to digital recording devices.

Nowadays, various aperture synthesis methods aiming at increasing DH’s numerical aperture have been proposed to improve DH’s lateral resolution^[8–21]. Most of them sacrifice or merely maintain the FOV instead of enlarging it. Only a few of them in DH achieved simultaneous improvements of both lateral resolution and FOV by adopting the hologram stitching method^[13,16]. Furthermore, the hologram stitching method has proved to be able to improve the axial measurement accuracy of DH as well^[22].

In the hologram stitching method, multiple holograms are recorded by shifting a DH system or a measured object^[13,16]. Then, multiple sub-holograms are stitched together to form an enlarged synthesized hologram. Two key issues in hologram stitching are (1) inconsistent phase errors among the sub-holograms; (2) errors of registration positions of sub-holograms between the real positions and the registered positions. Since the instability of the mechanical movement system is inevitable, the above errors perpetuate, which affect the effectiveness of stitching and could lead to failure of hologram stitching.

To overcome phase inconsistency, various techniques have been introduced^[23–28], including digital phase masks^[23–25], phase error model^[29], and physical phases^[26,28]. However, we would like to point out that, in hologram stitching, the compensation of phase errors is relevant to registration errors in the second issue. Instead, the compensation of phase errors is affected by registration errors. If the registration positions of sub-holograms are incorrect, the phase errors could not be correctly compensated. Hence, both registration errors and phase errors need to be correctly compensated.

To eliminate the registration position errors in the second issue, registration positions among sub-holograms are re-evaluated. Various methods have been proposed to re-evaluate the registration position. Cross-correlation (CC) methods based on similarity tests of overlapping areas in the reconstructed intensity images^[30–33] have been used to re-evaluate the registration position. An image-based method adopting the sharpness as the metric to search for correct registration positions has been proposed in Ref. [34]. Mutual information (MI) has also been proposed as a metric to automatically re-evaluate the registration positions in Ref. [35]. Dai et al. in Ref. [36] mentioned an image profile technique to obtain the location between adjacent sub-holograms, but no experiments were presented, and the images used have not been indicated as phase images or intensity images.

Currently, these above solutions are all based solely on intensity images to estimate registration errors. However, the intensity image is sensitive to the intensity distribution of the illumination light, which affects the correctness of the registration error estimation. Furthermore, the first issue is based on phase images, and the final 3D image of DH is derived solely from phase images. Hence, solving the two issues solely by phase images may have advantages and benefit the final 3D topography image result.

Therefore, in this work, a hologram stitching method solely based on phase images is proposed for aperture synthesis in DH, by which both registration position errors and phase errors among sub-holograms are corrected. Both the simulation and experiment are performed. The results show that, compared with the state-of-the-art intensity-based techniques, the phase-based scheme performs better in compensation of the registration errors among sub-aperture holograms and provides a better synthesized final 3D topography image result.

In this work, a lensless off-axis DH geometry is adopted, shown in Fig. 1. A HeNe laser of 633 nm wavelength and 10 mW power is used as a light source. The light is split into an object illumination beam and a reference beam by a beam splitter. The object beam Es reflected by the object interferes with the reference beam R reflected by a plane mirror to generate a hologram I on the CCD (with a pixel size of 4.65μm×4.65μm and a pixel number of 1280×960) plane. The CCD records the holograms.

After hologram recording, hologram reconstruction is performed. We apply Fourier filtering to the recorded hologram I and extract the object wavefront Es at the hologram plane. The reconstructed image of the object Eo can be calculated from Es by wave propagation from the hologram plane to the image plane by the angular spectrum method^[37] with a reconstruction distance z=119.7mm in this work. Then, the intensity Io and the phase φo of the reconstructed image at the image plane can be derived. Ideally, once the sub-holograms of all sub-apertures are collected by moving the stage, which is a two-axis translation stage (XR50C/M, Thorlabs) in this work, left to right, bottom to top, we can stitch them together by placing them at the designated positions to get a stitched hologram of a full synthetic aperture. After that, we follow the same reconstruction procedures as above.

However, in practice, the designated positions have errors, which make it difficult to generate a correct synthetic aperture hologram by directly synthesizing them. Furthermore, the temporal and spatial variations in the acquisition process also cause errors^[38]. Three major errors focused in this work are: •Phase error in each sub-aperture hologram: relative changes of the reference field to the object field caused by vibration and temperature fluctuations.•Registration error: the difference between real positions of sub-holograms and their designated positions in in-plane directions (x and y directions) caused by the inaccuracy of the mechanical movement of the translation stage.•Phase error between sub-holograms: phase error between different sub-holograms caused by the instability of mechanical movement of the translation stage in the out-of-plane direction (z-direction).

The procedures of the proposed aperture synthesis are presented as follows. Firstly, multiple sub-aperture holograms with overlapping regions are recorded, and the optical field Eo of each sub-aperture is acquired by numerical reconstruction. Secondly, the first major error discussed above, the phase errors φap in each sub-hologram, is estimated and compensated by the method described in Refs. [11–15,28] in this work. After that, the second major error, the registration error (e^up,e^vp), is estimated to get the correct registration position of adjacent sub-apertures [note that the spatial coordinates of the object plane and detector plane are defined as (x, y) and (u, v), respectively]. Then, the third major error, the phase error φepq between adjacent pth and qth sub-apertures, is estimated and compensated. We want to emphasize that the compensation accuracy of phase error φepq depends on the estimation accuracy of the registration error (e^up,e^vp). After the compensation of all three types of errors above, a large hologram is synthesized. Then, we reconstruct the synthesized optical field Eo, the intensity image Io, and phase image φo.

The procedures of the proposed method aiming to estimate the registration errors (e^up,e^vp) include four steps: Step 1. Set up a set of all possible registration errors (eupk,evpl) for the pth sub-aperture, where k and l are integer indices, satisfying −K≤k≤K and −L≤l≤L. K and L are the maximum errors in the u direction and v direction, respectively, which are determined by the movement accuracy of the translation stage.Step 2. The phase images φop(ui,vj) and φoq(ui′,vj′) of two neighboring sub-apertures (note that in the overlapping area with I×J pixels) are used to calculate the constant phase error φIJpq between the pth and qth sub-apertures by Eq. (1): φIJpq=1I×J[∑i=1I∑j=1JφIJp(ui,vj)−φIJq(ui′,vj′)].(1)The phase of the qth sub-aperture is compensated as φoq′(ui′,vj′)=φoq(ui′,vj′)+φIJpq.(2)Step 3. Assuming the qth sub-aperture moves with the designated distances du and dv with registration errors (eupk, evpl) relative to the pth sub-aperture in the u direction and v direction, respectively, we have ui=ui′+du+eupk,(3)vj=vj′+dv+evpl,(4)where ui′ and vj′ are positions of the pth sub-aperture in the u direction and v direction. The similarities between phase images φIJp(ui,vj) and φoq′(ui′,vj′) are calculated by various methods listed in Refs. [35,39,40], such as CC criterion, normalized CC (NCC) criterion, zero-NCC (ZNCC) criterion, MI criterion, sum of squared differences (SSD) criterion, normalized SSD (NSSD) criterion, and zero-NSSD (ZNSSD) criterion. We want to emphasize that the correlation coefficient CNSSD calculated by NSSD criterion can be derived from CNSSD=2(1−CNCC), and the correlation coefficient CZNSSD calculated by ZNSSD criterion can be derived from CZNSSD=2(1−CZNCC). By traversing all possible registration errors (eupk,evpl), a set of the corresponding similarity indicators Ckl is obtained.Step 4. The optimal estimation of the registration error (e^up,e^vp) is found by minimum or maximum values of the set Ckl obtained in step 3. For CC criterion, NCC criterion, ZNCC criterion, and MI criterion, the indicator of the similarity Ckl is proportional to the similarity extent. For SSD criterion, NSSD criterion, and ZNSSD criterion, the indicator is inversely proportional to the similarity extent.

The above algorithm with steps 1 to 4 is extended to the next neighboring sub-aperture one by one for aperture synthesis of a final whole aperture.

The process of state-of-the-art methods is based on the similarity of intensity images in the overlapping regions of nearby sub-apertures. It normally includes all steps except step 3 in the proposed method. Nevertheless, there are limitations in the similarity estimation by intensity images. The most important one is the uneven intensity distribution caused by the uneven distribution of laser light, non-parallel between the coordinate system of the CCD and the measured object, and the off-axis angle between the reference and object lights. These factors affect the estimation accuracy.

In holographic synthesis, if the registration position is correct, the difference between the synthesized and the pre-synthesized images should be minimal. To quantitatively evaluate the registration error estimated by different methods, an indicator, degree of image distortion (DID), is introduced herein, which is defined as DID=10log10∑(X−X¯)2∑(X−Y)2,(5)where X and Y are the pre-synthesized images and the synthesized images, respectively. X¯ is a mean of X. This indicator evaluates the difference between these two images. A larger DID value indicates a smaller difference between X and Y: when X=Y, DID→∞.

The procedure to verify the effectiveness of DID is shown in Fig. 2. A hologram measured in a real experiment is utilized to perform the simulation. Two sub-holograms with overlapping areas are truncated from it. We assign different registration positions and stitch them together. Then, we reconstruct the object field extracted from the synthesized hologram. Finally, the DID of the phase images before and after stitching is calculated. The DID of different assumed registration errors is shown in Fig. 3. The DID indicator reaches its maximum when the registration error is (0,0). It decreases with the increase of the registration error. Thus, DID proves to be an effective indicator to evaluate registration errors.

The experiment setup is shown in Fig. 1. The United States Air Force (USAF) target is used as a specimen to examine the proposed method. Nine sub-holograms (3×3) covering different parts of the USAF target are recorded by moving the translation stage and transferred into the computer. The registration errors are estimated by the proposed method discussed in subsection 3.1 as well as the state-of-the-art technique^{[30–32,34,35]}. Sub-holograms are synthesized to form a large hologram with error compensated registration positions. Then, they are reconstructed by the angular spectrum method^[37] with a reconstruction distance of z=119.7mm. Then, the intensity Io and phase φo of the reconstructed image at the image plane can be derived. The DID values of the synthesized phase images obtained by the state-of-the-art techniques and the proposed method are compared in Table 1 with different correlation criteria^{[30–32,34,35]}.

The synthesized phase images (the center part with 1280×1280 pixels) by different methods are presented in Fig. 4, with their corresponding DID values shown as well. From Figs. 4(a)–4(f), the DID increases from 44.20, 46.18, 47.61, 49.14, and 49.34 to the maximum value of 49.49 (found by traversing all possible positions for all holograms). The quality of the synthesized image increases as well. Note that the phase jumps (as the upper right corner image of each image shows) become smaller, and the high-resolution bars located in the center of the images (as the upper left corner image of each image shows) are better resolved. It is consistent with the results of the simulation performed above. Therefore, the rule that the higher the DID, the better the synthesized quality, is still applicable.

The synthesized images of the proposed methods shown in Figs. 4(d) and 4(e) are better than the synthesized images of the state-of-the-art methods shown in Figs. 4(a)–4(c) in phase jumps and in lateral resolutions in the center of images, where high-resolution bars are located. It is hard to find any difference between Figs. 4(e) and 4(f) through the method of subjective observation. Thus, the registration positions estimated by the proposed method could be considered as the optimal registration positions (the maximum DID value). To highlight the above difference between the proposed and state-of-the-art methods, their best results are compared in Fig. 5. Figure 5(a) is the highlight of Fig. 4(c), while Fig. 5(b) is the highlight of Fig. 4(e). In the right corner of Figs. 5(a) and 5(b), the high-resolution bars are better resolved in Fig. 5(b) than in Fig. 5(a). Let us further focus on their profiles in detail in Figs. 5(a) and 5(b), where only part of the profile was selected (marked with a dotted line). As shown in Fig. 5(c), the red dotted profiles of the state-of-the-art method present more severe phase jumps than the blue solid profiles of the proposed method. Hence, through subjective observations, the synthesized images of the proposed methods are better than the synthesized images of the state-of-the-art methods in both phase jumps and resolution.

The state-of-the-art method with SSD yields the worst result with the lowest DID value. This may be due to the inconsistency of the light intensity in the overlapping regions, which is caused by the uneven intensity distribution of the illumination light. By using the state-of-the-art method with NCC and ZNCC, the results improved as shown in Table 1. This is because, in the latter two criteria, where the intensities are normalized, they are less sensitive to the scale changes of intensity in the overlapping area. On the contrary, the phase images do not suffer such a problem, because the phase image is the reflection of the object surface shape unrelated to the distribution of the light illumination. Thus, the results obtained via the proposed method with all criteria are extremely similar and better than that obtained by the state-of-the-art method.

Therefore, based on the above discussions, we conclude that both subjective observations and objective quantitative DID evaluations verify that our proposed methods based solely on phase images are better than the state-of-the-art methods. Also, we note that though the off-axis DH was applied in this work, any geometry of DH, which is capable of obtaining an accurate phase map from a hologram, is suitable for the proposed method.