Variable shearing holography with applications to phase imaging and metrology

We report variable shear interferometers employing liquid-crystal-based geometric-phase (GP) gratings. Conventional grating-based shear interferometers require two lateral shifts of the gratings to enable phase-shifting capabilities in x- and y- direction and two axial shifts of the gratings to adjust the shear amount in x- and y-direction, i.e., these systems need control of four axes mechanically. Here we show that the GP gratings combined with a pixelated polarization camera give instantaneous-phase shifting so that no mechanical movement is necessary to obtain phase shifts. Furthermore, we show that a single fixed shear with a rotational shear axis provides a more robust selection of shears while requiring the control of only one mechanical axis. We verify this statement using spatial domain and frequency domain criteria. We further show that the resolution of the reconstructed wavefield depends not only on the numerical aperture of the imaging system, the pixel size of the detector, or the spatial coherence of the source but also on the ability to determine the shear amount accurately. To improve this, we report a methodology to accurately detect the shear amounts using the second derivative of the autocorrelation function of the measured holograms, which further relaxes the requirements on mechanical stability.


Introduction:
Single-beam wavefront reconstruction techniques are exciting alternatives to interferometry. These techniques are inherently robust to vibrations and often can be designed in a compact optical setup. Recent research expands those techniques to partial coherent sources, where LEDs are employed in place of lasers. Another interesting field is the area of bio-imaging, where existing microscopes are converted into phase measuring devices using a single-beam phase wavefront reconstruction add-on. A further appealing aspect is that many single-beam systems can be designed to work with no beamsplitter, allowing the sensor to be closer to the sample and enabling simpler optical measurements with higher-numerical aperture. In phase retrieval, there are two representative families of wavefront reconstruction techniques. Firstly, techniques based on deterministic models [1][2][3][4][5], and secondly, techniques that employ an iterative solver that utilizes wave-propagation techniques [6,7] or Fourier relationships [8,9].
With the advent of faster computer systems, new cameras are often developed with a higher pixel count, diminishing many advances of the higher computational power. Optical measurement technologies need to adapt to survive this trend through more innovative algorithms, new measurement concepts, or both.
For instance, in the field of phase retrieval, the convergence may be improved using multi-resolution and relaxation techniques [10]. Another possibility is to alter the illumination [11] or combine deterministic and iterative algorithms [12].
However, the elephant in the room of all methods mentioned above is the inverse problem resulting from the data produced by the unmodulated optical beam. In other words, often, the data is too similar, and significant efforts need to be made to reduce the similarity (e.g., repeating the measurement at a large defocus distance [3]).
Recent advances in computational shear interferometry [13,14] allow actively modulating an incident beam as it interferes with itself. Modulation of the beam increases instantly the diversity in the data. In particular, when employed with phase-shifting techniques [15], it is possible to obtain the complex crossterm of two interfering wavefields directly. Despite this sophistication, the system by Falldorf [13] requires a spatial light modular (SLM), which significantly contributes to the system costs. The wavefront reconstruction technique reported by Falldorf et al. [13] is based on a gradient descent approach for non-symmetric shearing systems. However, compared to the non-symmetric case, deriving the grading of a that non-holomorphic function for the case of a "symmetric" shearing system is not straightforward. Hence other approaches like the alternating projections by Konijnenberg [16] are preferred reconstruction algorithms.
Recent advances in geometric phase (GP) elements lead to the development of GP gratings and GP lenses that can be employed in lateral and radial shearing interferometry [17][18][19], with polarization phase-shifting capabilities [20,21]. The remaining quest is to recover the correct wavefront from the captured phase-shifted interferograms.
This paper proposes a variable shearing interferometer that consists of a series of GP gratings and polarization optics. This setup enables measuring a series of phase-shifted interferograms for a given shear amount, as shown in Figure 1. In a subsequent step, we reconstruct the wavefield using a modified version of the alternating projections algorithm, developed by Konijnenberg [16], and a series of interferograms with different amounts of shear. The discussed systems are the first experimentally demonstrated GP grating-based variable shear interferometers for complex wavefield reconstructions. We discuss various possible configurations of these instruments and demonstrate the success of two distinct shear selection strategies (employing one or two grating pairs). We analyze these findings using the "frequency information density" and (XXXX)X: XXX the "spatial information density" analyses based on transfer function and support function concepts introduced by Servín in [22]. We further show that a GP grating configuration with a fixed shear amount, but an adjustable rotational shear axis (one mechanical rotational axis) is sufficient to obtain robust measurements while maintaining phase-shifting capabilities using a pixelated polarization camera. We also investigate the effect of other error sources (e.g., imperfect shear estimation), which affects the resolution of the reconstructed wavefront. This effect is not commonly known (for coherent systems, the resolution is often a result of numerical aperture, pixel size, or spatial coherence of the source). To mitigate this source of error, we report a shear detection technique using the peaks of the second derivative of the autocorrelation function of the holograms. Measuring the shear amount directly from the hologram is different from relying on a well-calibrated system and increases the system's accuracy.

Results:
Shear Selection: Polarization gratings, also known as geometric phase (GP) gratings, are liquid crystal polymer-based optical gratings with a diffraction angle that depends on the polarization and wavelength of light. GP gratings are designed to diffract righthanded circularly polarized light (RHCP) and left-handed circularly polarized light (LHCP) at a negative angle and positive diffraction angle, respectively. The value of this angle depends on the prescription; commercially available gratings with various groove densities result in diffraction angles of the first-order between 3 and 10 degrees. The diffraction angle can be related directly to the wavelength and groove density using the well-known equation m λ = d sin(θ).
A lateral shear in the spatial domain is achieved when two GP gratings are placed consecutively [17], as shown in Figure 1. In that configuration, both identical GP gratings are aligned along the same axis and orientation. In this way, the first GP grating deflects the incoming beams at the angles +/-α, and the second GP grating compensates this deflection by diffraction the beam into -/+ α. The properties of these GP gratings in this configuration are outlined in [17]. Adjusting the distance between the first and the second GP grating allows adjusting the shear amount of the two interfering waves, as shown in Figure 2, enabling unique properties that can be used for adjustable shearing interferometry.
As shown in Figure 2, shearing along the y-axis is possible by flipping the gratings by 90 degrees. An instrument employing a series of four GP gratings enables shearing in both x-and y-direction. This configuration shears along the xand y-axis independently. This system uses 4 gratings in total: one pair for x-shear and the other pair for y-shear.
A different possible configuration is to maintain a fixed shear between two beams and modulate the fringes by rotating the shear axis about the optical axis, as shown in Figure 3. This configuration requires only one GP gratings pair and is more compact than the other configuration of Figure 2c.
Furthermore, as the linearly polarized light at the instrument's input is converted into left and right circular polarized components, there is a need to employ a linear polarizer to make these two waves interfere. Notably, this configuration provides geometric phase-shifting capability [23], i.e., by rotating the axis of the polarizer, it is possible to create a series of phase-shifted interferograms, as shown in Figure 4.
In the following sections, we propose using GP gratings for multi-shear variable shearing interferometry with phaseshifting capability. We further evaluate the performance of the instrument with mentioned above configurations using multishear wavefront reconstruction techniques.
Wavefront reconstruction algorithm for multi-shear variable shear interferometry: The previous section showed how GP gratings could obtain a series of phase-shifted interferograms for various shear amounts. These interferograms allow calculating the complex cross-terms of the two-interfering waves at each shear. The challenge is to reconstruct the complex wavefield from this data. In this section, we present an iterative algorithm based on the approach of Konijnenberg [16] that was inspired by the alternating projection algorithms [9,24].
Consider a general wavefield u with the amplitude A and the phase as u = A • exp(iϕ) (1) which produces the following interferogram , ( ⃗) in shearing, as Where ds ⃗⃗⃗⃗⃗ = [ds x , ds y ] and, ds x and ds y is the shear in xand y-direction, δ is the phase shift induced by the rotating polarizer, and ( ⃗) is the complex cross-term The term ( ⃗) is estimated using phase-shifting techniques, which for the four-bucket algorithm results in The reconstruction algorithm uses the term ( ⃗) that is obtained at different shears ' ⃗⃗⃗⃗ ' to reconstruct the object wavefield through an iterative process. In this work, we employ an alternating projections (AP) approach developed by Konijnenberg [16] for the x-ray community. The algorithm applies two distinct constraint sets, where each constraint allows only a limited family of solutions. The objective of the wavefield reconstruction algorithm is to identify the unique solution that satisfies all constraints (and for all shears). These constraints are satisfied by employing alternating projections between the two solution domains, in which we apply both, socalled "object constraints" and so-called "measurement constraints." Switching between "object constraints" and "measurement constraints" is applied continuously until the algorithm converges to a specific wavefield. The convergence can be conveniently estimated using the following loss function where subscript n denotes that the wavefield uses the shear ⃗⃗⃗⃗⃗⃗⃗ , S is the total number of shears, and for convenience, we define the estimated wavefield after N iterations as ( ⃗) to distinguish it from ideal wavefield ( ⃗). The notation ( ⃗) is used as a general indication to represent estimated wavefields from both constraints, where the notation ( ⃗) denotes the estimate of the wavefield directly after the measurement constraint is applied, and ( ⃗) denotes the estimate of the wavefield directly after the object constraint is applied. For further simplification, the symbols 'n,+' and 'n,-' in the subscript are used to represent positive and negative shears, i.e., ,± ( ⃗) = ( ⃗ ± ⃗⃗⃗⃗⃗⃗⃗ ). Consequently, the estimated complex cross amplitude term can be presented as and can be generated using the estimated wavefield f M n (x ⃗⃗) at any shear ⃗⃗⃗⃗⃗ .
Notably, in the presence of noise, there is no exact solution. However, the wavefield reconstruction algorithm estimates the best fit solution that matches the experiment data. The loss function serves as a metric to be observed for convergence with an increasing number of iterations. Figure 5 shows the steps of the wavefield reconstruction algorithm. After an initial guess, the algorithm continuously applies the so-called object and measurement constraints alternatingly until the loss function reaches the desired level.
The application of those constraints is discussed in more detail below.
Estimation of f M n,± ( x ⃗⃗ ) using the initial guess: This step of the algorithm uses only the initial guess f( x ⃗⃗ ) in the first iteration to estimate f M n,± ( x ⃗⃗ ) . An initial guess is assumed for f( x ⃗⃗ ) (e.g., a spherical wave, a plane wave, or a random wavefield) and is laterally sheared by the same shears used in the experiment. For a wavefield f( x ⃗⃗ ), the shearing is done in the Fourier domain using the following equations.
Assuming that the initial guess wavefield is from the measurement constraint domain, the initial guess and its sheared components are indicated as f M n,± (x ⃗⃗). A wavefield with constant phase and unit amplitude is used as an initial guess for the reconstructions done in this paper.

Application of the Object constraint:
The object constraints are applied by transforming the sheared components of the estimated wavefield to the Fourier domain and compensating the shear for each f M n,± (x ⃗⃗) to obtain f( x ⃗⃗ ) using the Fourier Shift theorem. All calculated values for F( k ⃗⃗ )are averaged to obtain a more accurate result.

ACCEPTED ARTICLE PREVIEW
(XXXX)X: XXX Finally, F est (k ⃗⃗ ) is transformed back into the spatial domain using the inverse Fourier transform. However, because the measurement constraints (applied in the second part of the iteration) require the wavefield for different shears, we apply again the Fourier shift theorem to obtain f O n,± (x ⃗⃗) as Application of the Measurement constraint: The measurement constraint equations use the estimated wavefield from the object constraint f O n,± (x ⃗⃗) and the measured complex cross-term C n (x ⃗⃗) (obtained from phase shifting) to estimate the wavefield solution in the measurement constraint domain. The following equations are a modified version of the equations used by Konijnenberg in [16] to estimate the wavefield.
Behavior of the algorithm: The algorithm uses two sets of constraint equations, object constraint (Eq. 10), and measurement constraint (Eqs. 12 and 13). Each set of constraint equations has a solution set. The algorithm starts with an initial guess projected onto the object constraint solution set by using the respective equations. This produces an estimated wavefield which is a projection on the object constraint solution set. This wavefield is then projected onto the measurement constraint solution set using the respective equations, and the process repeats until convergence. Ideally, the solution we are looking for is the one that is at the intersection of these two solution sets. However, in practice, due to noise and other external factors, such intersection does not exist, and the algorithm converges when the projections oscillate between the two solution sets with solutions at the closest proximity to each other.
Simulation using synthetic data: The performance of the algorithm was evaluated by simulations using MATLAB. Preliminary simulations were made for a resolution of 256 x 256 pixels. For this simulation, a point source was generated with a wavelength of 600 nm at a z-distance of 87.5 mm and was sampled with a pixel size of 3.45 μm. The reconstruction algorithm uses FFTs to shear the estimated wavefields at different steps in the reconstruction. The use of FFTs, in turn, implicitly imposes a periodic boundary condition. Therefore, we use an aperture for both simulations and experiments to mitigate problems during reconstruction that might arise from periodic properties of FFTs. For this reason, a circular aperture mask with a diameter of 180 pixels was applied. Future reconstruction algorithms may exclude samples that exceed the computational window, especially at larger shear. Those samples need to be reconstructed with measurements taken at smaller shear amounts.
A series of four phase-shifted intensities were generated for each shear position. The shears generated for this simulation are shown in Figure 6a. Additive white Gaussian noise (AWGN) with a standard deviation of 20% of the peak intensity has been added to the generated intensities to model the effect of noise in these simulations. Afterward, to realistically model the effect of cameras, an 8-bit discretization has been applied to the intensities.
For this simulation, a total of 48 interferograms have been generated corresponding to 12 different shears with 4 phaseshifted interferograms. The ideal wavefront and reconstructed wavefront are shown in Figures 6b and 6c, respectively. The corresponding phase error of the reconstructed wavefront is shown in Figure 6d. Further investigation of the error map showed the presence of high-frequency artifacts, as highlighted in Figure 7.
These high-frequency artifacts occur at Nyquist frequency and appear to be a product of the reconstruction algorithm. Further analysis is required to investigate the source of these artifacts. For the reconstructions discussed in this paper, we use a low-pass Fourier filter to mitigate the formation of these artifacts. Figure 8 shows error maps generated by mapping the difference between the ideal and reconstructed wavefield before and after applying the filter. A cross-section of the error profile before and after filtering is shown in Figures 8c and 8d, respectively. This simple procedure reduces the magnitude of the RMS error by one order of magnitude. Because a Fourier low-pass filter has been employed with a corner frequency of = ( , an artifact is produced along the boundary of the aperture. However, other filtering techniques may be applicable depending on the application, and even Zernike polynomials may be fitted to the data. The use of different filtering techniques is beyond the scope of this work.

Evolution of loss function:
As mentioned earlier, the evolution of the loss function indicates how close the algorithm is to achieving the best solution or achieving convergence. The loss function is evaluated in these experiments by calculating the change between the solutions generated by the two constraints on positive and negative sheared estimated wavefields. The loss function was already defined in Eq. 5 and is reiterated here in Eq. 14 for convenience. Figure 9a shows the evolution of loss function during reconstructions done in Figure 8 before the filter was applied.
Figure 9a shows that the loss function converges after 30 iterations of the reconstruction without a Fourier filter. Given the filter's effectiveness in Figure 8, we have investigated the use of the filter within the reconstruction algorithm. The sudden surge at the start of the loss function in Figure 9 is a characteristic feature of phase retrieval and reconstruction algorithms. The temporary increase of the value of the Loss function is expected because the solver climbs out of the local minimum. Other phase retrieval algorithms such as the "hybrid input-output algorithm" [8] show a similar behavior.
Two cases were considered for this study: (i) using the filter after every 5 iterations and (ii) applying the filter in every iteration. Figure 9b shows that applying a filter enables a drop in the convergence values in both cases, thus indicating an improvement in the reconstruction process. Interestingly, the reconstruction algorithm for case (i) consistently converges back to the solution containing the high-frequency artifact. Therefore, we recommend employing a Fourier filter within the wavefront reconstruction algorithm to suppress all undesired artifacts continuously.
Please note, the algorithm described is a modified version of the alternating projections algorithm developed by Konijnenberg et al. in [16]. For better clarity we summarize the differences to the original algorithm from [16].
• The data recorded from experiments in [16] are in the Fourier domain. In this paper, all recorded data from experiments are in the spatial domain. were defined in the Fourier domain. In this paper, these equations are defined in the spatial domain. • The algorithm in [16] uses two sets of equations in measurement constraint (Eqs. 15 and 16 in [16]), roughreconstruction equations for initial convergence (Eqs. 16 in [16]) and fine-reconstruction equations for final refinement (Eqs. 15 in [16]). This paper only uses the modified versions of the fine reconstruction equations. • After initial convergence, we include a Fourier filter in the measurement constraint equations to mitigate the effects of high-frequency artifacts. • Multi-resolution technique [10] is incorporated into this algorithm to mitigate the effects of initial guess and improve convergence.

Comparison of the two configurations: (Simulation)
In the cartesian coordinate shear system, there is the freedom to choose the shear in the X and Y direction independently. However, the optical configurations shown in Figure 2 are restricted to one quadrant due to geometrical constrains. In principle, this limitation may be overcome using a 4f-optical setup (obtain positive and negative shears for a grating pair). However, this configuration increases the size of the system. There is freedom in the polar coordinate shear system to choose shears across all four quadrants, but these shears lie on a circle with a constant radius. When these shear systems are implemented experimentally, the shear amount cannot be controlled directly and estimated.
The choice of shear selection has an impact on the reconstruction of the wavefield. In general, increasing the number of shears provides more information on the object (XXXX)X: XXX wavefield for the algorithm to reconstruct the wavefield. However, improper selection of shears could also lead to the algorithm not having sufficient frequencies to reconstruct the object wavefield, thus resulting in artifacts. We want to highlight this typical behavior of this system via simulation using a few exemplary sets of shears, as shown in Figure 10.
Each simulation was performed with 24 shears, where the selections are shown in Figures 10a, 10b, and 10c. The selections analyze two possible xy-shearing configurations shown in Figure 10a and 10b (see Figure 2 for configuration) and the fixed shear with variable shear axis shown in Figure  10c (see Figure 3 for configuration).
In test 1 (Figure 10a), the shears have been selected along an approximately linear trend. The reconstructed phase map from this shear selection showed that some of the frequencies were not being reconstructed correctly, which is a consequence of the ambiguity that arises from the measured data. The reconstruction results are shown in Figures 10d and  10g. When the selection of shears was scattered across one quadrant (Figure 10b) as in test 2, the reconstruction algorithm had sufficient data to reconstruct most of the point-source wavefield, as shown in Figure 10e. However, the error map in Figure 10h shows a curvature, indicating the possibility of some missing frequencies. This slowly varying error is critical as it produces a so-called form error that is difficult to filter out when observing aperture size measurement objects (e.g., mirrors or lenses under test). This error may be less critical when the sample is small (e.g., in the case of a cell under a microscope). The shear selection shown in Figure 10c is based on a variable shearing axis but with a fixed shear amount. The latter system allows accessing all four quadrants, and the shears can be selected at all possible rotation angles of the shear axis. The reconstruction is shown in Figure 10f. The error map shown in Figure 10i shows that the error map is dominated only by highfrequency noise patterns, which can be filtered more easily. These results demonstrate that selecting shears across all quadrants results in better reconstruction because it provides more information for a successful reconstruction.
(XXXX)X: XXX Transfer function analysis (see also [22]), to estimate frequency domain information densities: To better understand the results and characteristics features in the error maps of Figure 10, further study was done using transfer function analysis to investigate the effectiveness of the selected shear sets to transmit information for wavefield reconstruction. Extensive studies have been done by Falldorf [25] for non-symmetric shears and Servin [22] for symmetric shears in understanding the frequency response of spatially shearing systems. Inspired by Servin's description [22] for the transfer function, we define the frequency information density as where ⃗⃗ = [ , ] is the frequency coordinate and ⃗⃗⃗⃗⃗ = [ , ] is the shear in x-and y-directions. This transfer function analysis estimates the "frequency information density" at each spatial frequency coordinate for a given set of shears. These functions were evaluated and mapped for all three configurations as shown in Figure 11. Figures 11(a), (b) and (c) are the frequency domain information density maps for the three shear sets described in Figure 10. Figures 11 (d), (e) and (f) are the same maps but zoomed in near low spatial frequencies to visualize better the behavior at these regions, where regions below the threshold of 2 are highlighted in red. Figure 11d shows a significantly high amount of frequency sets that have transfer function values below the threshold. These cluster appear to dominate in regions where kx ≈ ky. This indicates that a large number of frequencies along the 45-degree angle will have poor reconstruction in the estimated wavefield at longer spatial wavelengths. This conclusion is in agreement with Figure 10g. A similar behavior is observed in Figures 11b and 11e, however, the number of frequency sets with transfer function values below the threshold of 2 are located at both low and high frequencies, with a minimal presence near low spatial frequencies. The locations of these cluster indicate that very low spatial frequencies will appear in the error map but will have better reconstruction than the previous shear set. This is observed in Figure 10h, but there are also high frequency errors present due to low values of the transfer function at higher spatial frequencies. The polar configuration in Figure 11f has no transfer function values below the threshold of 2, except for the (trivial) frequency at (0,0). This indicates good reconstruction of the estimated wavefield at all spatial frequencies.
(XXXX)X: XXX Spatial support analysis (see also [22]), to estimate spatial domain information densities: The following section demonstrates how spatial overlaps between different shears impact the reconstruction of the final wavefield. When an object wavefield passes through the gratings, it splits in two, represented as beam 1 and 2 in Figure  12a. Depending on the shear amount, there is an overlap region between the two beams (overlap region is indicated in yellow). Since the overlap region produces an interference pattern, the phase information from the overlap region feeds into the algorithm. Figure 12b shows which areas of the input beam provide phase information to the reconstruction algorithm (regions indicated in yellow). This idea is further extended to multi-shear interferometry in Figure 13. Figure 13a shows the overlap regions from an object beam with a diameter of 600 pixels sheared by 175 pixels along the x-direction. Figure 13b shows the overlap regions if a second shear is added to the measurement in the orthogonal direction by the same amount. In Figure 13b the yellow areas indicate regions with two overlaps, green areas indicate regions with one overlap, and the dark blue areas indicate regions with no overlap. This idea is extended over 24 shears in Figure 13c, showing a significantly higher count of overlaps (or number of effective measurements). This case also shows an example of a situation where the algorithm will not completely reconstruct all spatial frequencies near the wavefront's center because the chosen fixed shear is too large (175 pixels). Figure 13d shows an alternative set of 24 shears (i.e., with a small shear amount compared to (c) with 80.5pixels). There, all of the spatial locations have been measured at least 20 times. The case presented in Figure 13d has the highest potential to successfully reconstruct the wavefield due to maximum overlap across the whole area. It should be noted that (as discussed by Servin [22]) the wavefront may be reconstructed even at regions where there is no overlap. This phenomenon is especially true for lower spatial frequencies because they do not need to be sampled everywhere. However, the reconstruction is increasingly difficult at higher spatial frequencies at a given location if no interference exists. The same concept is applied to the three configurations / shear sets discussed in Figure 10 and the overlap distribution is mapped in Figure 14. Figures 14 (a) and (b) show fewer overlaps along the boundary at 45-degree angle. The number of effective measurements drops from 24 to 13. On the contrary, Figure 14c shows uniform distribution of overlaps covering the majority area of the object beam. At the border, there are still 20 effective measurements recorded, which indicates better reconstruction across the whole area of the object.

Effect of pixel shear errors on the resolution of the reconstruction (Monte Carlo simulation):
Any error present in the estimated shear value propagate through the reconstruction process and result in an error in the final reconstructed wavefield. This mismatch between actual and assumed shear amounts can result in error artifacts in the reconstructed wavefield or missing frequencies from the original wavefield. Monte-Carlo simulations were carried out to evaluate the effect of errors in the estimated shear pixels in reconstructing the wavefield by adding normal distributed random errors to the shear values in X and Y direction for the shear configuration in Figure 10c (100 simulations for each case). Furthermore, a total of 20% intensity noise (of the peak value) has been added to all interferograms to include the effects of noise in realistic measurement conditions. This series of simulations provides an overview of how the phase errors in reconstructed wavefields are affected by the imperfect shears. The results of the series of simulations have been summarized in a box plot, as shown in Figure 15. As expected, when observing the RMS error in the reconstructed phase, there is a clear correlation with the increasing amount of shear error. However, it is interesting to notice that the different spatial frequencies are affected in different ways; as the standard deviation of the shear error increases, the ability of the algorithm to reconstruct higher spatial frequencies decreases. Figures 15 (b), (c), and (d) show the measurement of a USAF target where the Monte Carlo simulation used a shear error of 1-, 3-, and 5-pixels standard deviation, respectively.
Additional processing of measurement data needed for experimental data: Shear detection: The reconstruction algorithm requires the shear amount to be known. As concluded from the Monte Carlo simulations, the better the estimation of the shear amount, the lower the reconstruction error and the higher the spatial resolution. In this work, we have estimated the shear amount by calculating the auto-correlation of the recorded intensity maps. To increase the robustness of this estimation, we have removed the sinusoidal fringe modulation by averaging all four measured phase-shifted intensity patterns. One examplary measurement is shown in Figure 16a. The second-order derivative of the auto-correlation map has been computed along one direction to further enhance the peak in the autocorrelation map. The resulting spikes for the case in Figure  16a are shown in Figure 16b, and they are prominently identifiable. The second-order derivative of the autocorrelation map comprises three peaks: a center peak and a peak on either side of the center peak, which is located at twice the value of the actual positive and negative shear amounts. The shear amount is estimated by tracking one of the side peak positions to the center peak. The difference in X and Y coordinates of the side and center peak gives twice the shear value for each averaged intensity map. The actual shear can be detected down to an accuracy of 0.5 pixels.

Additional compensation techniques for the polarization camera:
When working with left and right circular polarized waves, it is convenient to employ pixelated polarization cameras to record four phase-shifted intensities simultaneously. The polarization camera enables single-shot phase-shifting during the measurement process for each shear. The single-shot phase-shifting capability offers an advantage in reducing the number of translation axes, enhanced stability, and mitigation of alignment errors. The reduction of the number of translation axes further ensures repeatability.
Commonly, each super-pixel on the detector of the polarization camera has four pixels that record the image at four different polarizer angles. As the two incident waves have a left-circular and right-circular polarization, it is possible to obtain simultaneously four interferograms at four different phase shifts. The orientation of this pixelated array is shown in Figure 17, which also shows an example of a defocused point source measured at a single fixed shear with 0°, 90°, 180°, and 270° phase shifts. The recorded intensities have a resolution of 1024 x 1024 pixels.
One critical correction is needed before the phase can be calculated. The individual frames are often clustered in a superpixel and can be separated without difficulty. However, all four frames are not located at the same spatial position. An additional interpolation is needed (e.g., including the Fourier shift theorem) to obtain all four interferograms at the same spatial location. In this work, we have compensated this spatial mismatch for the interferograms I2, I3, and I4 (but not I1), to obtain an interferogram that is shifted by ½ super-pixels.

Experimental setup:
Cartesian Coordinate Shear configuration: Figure 18 shows the experimental setup for the Cartesian coordinate shear setup, along with a schematic showing the arrangement of optical components for the setup. The light source illuminating the object is a class 3B, 520 nm wavelength, 30 mW green laser. The object used for the experiments is a USAF target. Two sets of lenses with 4f configuration were used to keep the aperture 1 and the object in focus. Aperture 2 was used to limit the effects of stray light, back reflections, and the higher-order diffraction orders. A polarization camera was used to capture the four phase-shifted interferograms simultaneously (FLIR Blackfly-S model BFS-U3-51S5P-C). The grating arrangement is shown in detail in Figure 19.
The cartesian coordinate shear system in Figure 19a has four polarization gratings, two on the left-hand side (blue) to shear the wavefield in the vertical direction (y-direction), and the two on the right-hand side (green) to shear the wavefield in the (XXXX)X: XXX horizontal direction (x-direction). Two leadscrew stages are connected to the two outermost X and Y polarization gratings which translate along the Z-axis to create the shears. The leadscrew stages utilize stepper motors and encoders to accurately determine the stages' position and limit switches for homing the system. Additive manufacturing (3D printing) was used to construct the stage, and accompanying electrical circuitry and software were developed for the closed-loop control system. The two inner gratings are fixed. Since only one grating can move in each grating pair, the shears are limited to only one quadrant (see Figure 10a and 10b).
An experiment was carried out to validate the previous findings regarding the shear configuration in Figure 10a. The results for the first configuration are shown in Figure 20, where a total of 24 different shear amounts have been applied (see Figure 20a) using an instrument that can independently adjust the x-and y-shears. The object used in this experiment is a simple section of an amplitude USAF target, where the (XXXX)X: XXX remaining sections are blacked out due to aperture 1 (see Figure 18a).
The shears were estimated from experimental data using the auto-correlation method, as discussed in the previous section. Once the shears were estimated, they were given as input to the algorithm to reconstruct the wavefield. A multi-resolution technique [10] was implemented within the reconstruction algorithm to mitigate the effects of initial guesses and improve the convergence. For our dataset with 1024 x 1024 pixels, a 5level multi-resolution reconstruction starting from the coarsest grid with 64 x 64 pixels was implemented. The estimated wavefield is then up-sampled by a factor of 2, and the reconstruction process happens at a higher resolution (i.e., 128 x 128 pixels) using the newly up-sampled wavefield as the initial guess. This process repeats until the estimated wavefield reaches the original resolution, i.e., 1024 x 1024 pixels. The results for this 5-layer multi-resolution approach are shown in Figure 20.
The reconstructed amplitude map (Figure 20b) shows the object features that would be expected. However, the phase map in Figure 20c shows the presence of parasitic signals. This reconstruction is expected given the outcome of previous simulations, see Figure 10. This is a direct consequence of the shear selection: most of the shears selected for this experiment are dominant along the -45° angle. This selection has generated parasitic signals along the 45° angle, as seen in Figure 20c.

Constant shearing Polar coordinate shear configuration:
The other shear configuration that produced consistently better results was the configuration of a fixed shear with a variable shear axis, see Figure 10c. An experimental implementation is shown in Figure 21, where the xy-shear system in Figure 18 is replaced with a simple two grating set up on a rotational mount. Figure 22 shows a more detailed view of the shearing mechanism. Notably, the 2-grating setup can be realized in a more compact optical setup; however, we chose to keep the 4f configuration to compare the two setups better. The results for the configuration in Figure 21, using 24 different shear amounts, are shown in Figure 23. Similar to the case of Figure 20, a portion of the USAF target has been selected for the measurement, and the field of view is limited by aperture 1 (see Figure 21). Figures 23a and 23b show that the solver successfully reconstructs both the amplitude and the phase map, free of parasitic fringes. For better visualization, a mask has been generated (see Figure 23c) using the amplitude information and applied to the phase data (see Figure 23d).
A further experiment has been conducted to verify that the instrument in Figure 21 can capture and reconstruct the wavefield at a defocused distance. For this purpose, we shift the camera's position by a distance of ~170mm, so that both the object and the aperture get defocused by the same amount.
The reconstructed wavefield is shown in Figures 24a and 24c. After reconstruction, the wavefield was numerically refocused to the visually best reconstruction distance (z=162.94 mm) using the Angular Spectrum method [26,27]. The results of the refocused wavefield are shown in Figures 24b and 24d. For better visualization, various regions of the reconstructions have been zoomed in and are presented in Figure 25.

Materials and methods:
The light source used for the experiment is an FC-520-040-PM laser system (520 nm, 30 mW) from CivilLaser. The object sample is a combined USAF 1951 and Dot Grid Target (Edmund Optics #62-465). The GP gratings have been purchased from Edmund Optics (#12-677, 96% diffraction efficiency at 550nm, 159grooves/mm, 25 x 25 x 0.45 mm). The light illuminating the object is linearly polarized to produce 50% RHCP and 50% LHCP at the GP grating, maximizing the fringe visibility. The camera used for imaging is a FLIR Blackfly-S model BFS-U3-51S5P-C polarization camera. The lenses used for the 4f-configuration have an effective focal length of 50mm.

Conclusion:
In this work, a GP grating-based shearing holography system has been presented. The GP gratings used in this work have the unique ability to have a polarization-sensitive diffraction angle, where the RHCP and LHCP light beams are diffracted in two different directions. Figures 1, 2, and 3 show how these properties can be exploited with a linear polarizer to create an adjustable shearing interferometer. This paper presented two distinct configurations, (i) one that allows for adjustable xy-shearing, and (ii) one that keeps the shear amount constant but rotates the shear axis. Another aspect of the GP grating-based shearing interferometer is that a series of phaseshifted interferograms can be obtained by either rotating the polarizer (see Figure 4) or using a pixelated polarization camera (see Figure 17).
The ability to produce a series of phase-shifted interferograms obtained under various shear amounts opens the possibility to obtain high fidelity data that can be processed in a wavefront reconstruction algorithm. In this work, we have implemented the approach of Konijnenberg [16] outlined in Figure 5 but added a filtering procedure to remove algorithmspecific artifacts (see Figures 6d, 7, 8, and 9). We further investigated the effect of the shear selection. Figure 10 showed that choosing the shears to be on a circle in the sx-sy plane provides consistently good results, where the error is of highfrequency nature. These findings are supported using Servin's [22] transfer function analysis and spatial domain support that were referred to as "frequency domain information density" (see Figure 11) and "spatial domain information density (see Figure 14). One critical step is the estimation of the exact shear amount. Any discrepancy between the actual shear amount and the values used within the wavefront reconstruction algorithms will increase, to some degree, the RMS in the reconstructed phase, see Figure 15. However, the results in Figure 15 also show that an imperfect shear estimate value has a strong impact on the resolution of the reconstructed wavefield (even if the RMS increases only marginally by 10%, for the case of typical levels of intensity noise). To ensure accurate results, we have reported a shear estimation technique using the second derivative of the autocorrelation function, where the spikes of that function estimate shear amount down to an accuracy of 0.5 pixels, as shown in Figure 16. The reconstruction algorithm presented in this paper does not limit the possible applied shears. The algorithm can process sub-pixel shearing as the algorithm primarily operates on FFTs. However, the shear detection technique presented in this paper is limited in identifying shears. The procedure detects twice the values of shears for a sheared interferogram as integer, and hence the shears can be estimated down to integer multiples of 0.5 pixels. More sophisticated shear detection algorithms (using weighted mass algorithms or techniques based on ideas from Guizar-Sicairos M. et al. [28]) might enable finer shearing detection down to 0.01 pixels and are subject to future research.
The previous aspects have been verified experimentally for both shearing configurations, as shown in Figures 18 and 19 for the xy-shearing component and Figures 21 and 22 for the configuration that rotates the shear axis. As expected, the system configuration that keeps both gratings fixed while the grating pair rotates about the optical axis provides the visually best reconstruction result with no parasitic fringes, as shown in Figure 20 and 23. A further advantage of the setup with the variable shear axis is that it requires only two GP gratings and can be, in principle, incorporated into a more compact setup. There are also several motorized rotating mounts available that could automate this process.
Finally, we demonstrate the feasibility of this approach by reconstructing the wavefield of a defocused object at the distance z=162.94 mm. Refocusing using the Angular Spectrum method [26,27] produces Figures 24 and 25.
These findings show the advantage of GP grating-based shear interferometers compared to conventional dual grating (or dual double-grating) configurations. A conventional dualaxis grating-based variable shearing interferometer with phase shift capabilities requires two independent axial shifts of the x-and y-gratings (to adjust the shear) and lateral shifts of the gratings in x-and y-direction (for phase-shifting). Having to control four mechanical axes independently is a significant system complexity, where the freedom to adjust all parameters appropriately is often compromised due to space constraints. In contrast, the results of this work show that a GP grating configuration with a fixed shear amount, but an adjustable rotational shear axis can provide the shear selections for robust measurements while maintaining phase-shifting capabilities using a pixelated polarization camera. In other words, robust variable shearing interferometry can be made possible by controlling only one mechanical axis (instead of four). This finding has implications for future compact interferometer designs.