
Polarization gratings, also known as geometric phase (GP) gratings, are liquid crystal polymerbased 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 lefthanded 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 firstorder between 3 and 10 degrees. The diffraction angle can be related directly to the wavelength and groove density using the wellknown equation m λ = d sin(θ).
A lateral shear in the spatial domain is achieved when two GP gratings are placed consecutively^{17}, as shown in Fig. 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 Ref. 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 Fig. 2, enabling unique properties that can be used for adjustable shearing interferometry.
Fig. 2 Examples of interferograms for the case of a spherical wave at the distance z = 350 mm with different shear amounts in x and y direction.
As shown in Fig. 2, shearing along the yaxis is possible by flipping the gratings by 90 degrees. An instrument employing a series of four GP gratings enables shearing in both x and ydirection. This configuration shears along the x and yaxis independently. This system uses 4 gratings in total: one pair for xshear and the other pair for yshear.
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 Fig. 3. This configuration requires only one GP gratings pair and is more compact than the other configuration of Fig. 2c.
Fig. 3 Examples of interferograms for the case of a spherical wave at the distance z = 250 mm with constant shear and different grating rotational angles.
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 phaseshifting capability^{23}, i.e., by rotating the axis of the polarizer, it is possible to create a series of phaseshifted interferograms, as shown in Fig. 4.
Fig. 4 a Intensity ditribution of a spherical wavefront with constant amplitude before reaching the linear polarizer. b Intensity after the linear polarizer for four different angles of the polarizer axis to produce a phase shift between interfering waves of 0, π/2, π, and 3π/2.
In the following sections, we propose using GP gratings for multishear variable shearing interferometry with phaseshifting capability. We further evaluate the performance of the instrument with mentioned above configurations using multishear wavefront reconstruction techniques.

The previous section showed how GP gratings could obtain a series of phaseshifted interferograms for various shear amounts. These interferograms allow calculating the complex crossterms of the twointerfering 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
$$ {\rm{u}}={\rm{A}}\cdot {\rm{exp}}\left({\rm{i}}{\rm{\phi}}\right) $$ (1) which produces the following interferogram
$ {I}_{n,\delta }\left(\overrightarrow{x}\right) $ in shearing, as$$\begin{split} {{\rm{I}}}_{{\rm{n}},{\rm{\delta}}}\left(\overrightarrow{{\rm{x}}}\right)=&{\left{\rm{u}}\left(\overrightarrow{{\rm{x}}}+\overrightarrow{{\rm{d}}{s}_{n}}\right)\right}^{2}+{\left{\rm{u}}\left(\overrightarrow{{\rm{x}}}\overrightarrow{{\rm{d}}{s}_{n}}\right)\right}^{2} \\ & +2{\cal{R}}\left\{{{\rm{C}}}_{{\rm{n}}}\left(\overrightarrow{{\rm{x}}}\right){\rm{exp}}\left({\rm{i}}{\rm{\delta}}\right)\right\}\end{split}$$ (2) Where
$ \overrightarrow{{\rm{d}}{\rm{s}}}=\left[{{\rm{d}}{\rm{s}}}_{{\rm{x}}},{{\rm{d}}{\rm{s}}}_{{\rm{y}}}\right] $ and,$ {{\rm{d}}{\rm{s}}}_{{\rm{x}}} $ and$ {{\rm{d}}{\rm{s}}}_{{\rm{y}}} $ is the shear in x and y direction, δ is the phase shift induced by the rotating polarizer, and$ {C}_{n}\left(\overrightarrow{x}\right) $ is the complex crossterm$$ {C}_{n}\left(\overrightarrow{x}\right)={u}^{*}\left(\overrightarrow{x}\overrightarrow{d{s}_{n}}\right)\cdot u\left(\overrightarrow{x}+\overrightarrow{d{s}_{n}}\right) $$ (3) The term
$ {C}_{n}\left(\overrightarrow{x}\right) $ is estimated using phaseshifting techniques, which for the fourbucket algorithm results in$$ {C}_{n}\left(\overrightarrow{x}\right)=\left({I}_{n,0}\left(\overrightarrow{x}\right){I}_{n,\pi }\left(\overrightarrow{x}\right)\right)+i\left({I}_{n,\frac{\pi }{2}}\left(\overrightarrow{x}\right){I}_{n,\frac{3\pi }{2}}\left(\overrightarrow{x}\right)\right) $$ (4) The reconstruction algorithm uses the term
$ {C}_{n}\left(\overrightarrow{x}\right) $ that is obtained at different shears ‘$ \overrightarrow{{s}_{n}} $ ’ 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 xray 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 socalled “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
$$\begin{split} {\rm{L}} =&\sum\limits_{{\rm{n}} = 1}^{{\rm{S}}}\bigg({\left{{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right){{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)\right}^{2}\\&+{\left{{{\rm{f}}^{{\rm{M}}}_{{\rm{n}},}}}\left(\overrightarrow{{\rm{x}}}\right){{{\rm{f}}^{{\rm{O}}}_{{\rm{n}},}}}\left(\overrightarrow{{\rm{x}}}\right)\right}^{2}\bigg)\end{split} $$ (5) where subscript n denotes that the wavefield uses the shear
$ \overrightarrow{d{s}_{n}} $ , S is the total number of shears, and for convenience, we define the estimated wavefield after N iterations as$ f\left(\overrightarrow{x}\right) $ to distinguish it from ideal wavefield$ u\left(\overrightarrow{x}\right) $ . The notation$ f\left(\overrightarrow{x}\right) $ is used as a general indication to represent estimated wavefields from both constraints, where the notation$ {f}^{M}\left(\overrightarrow{x}\right) $ denotes the estimate of the wavefield directly after the measurement constraint is applied, and$ {f}^{O}\left(\overrightarrow{x}\right) $ 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.,$ {f}_{n,\pm }\left(\overrightarrow{x}\right)=f\left(\overrightarrow{x}\pm \overrightarrow{d{s}_{n}} \right) $ . Consequently, the estimated complex cross amplitude term can be presented as$$ \begin{split} {{\rm{C}}}_{{\rm{e}}{\rm{s}}{\rm{t}},{\rm{n}}}\left(\overrightarrow{{\rm{x}}}\right)=&{{\rm{f}}}^{*}\left(\overrightarrow{{\rm{x}}}\overrightarrow{d{s}_{n}}\right)\cdot f\left(\overrightarrow{{\rm{x}}}+ {\rm{d}}\overrightarrow{{s}_{n}}\right) \\ =&{{{\rm{f}}}_{{\rm{n}},}\left(\overrightarrow{{\rm{x}}}\right)}^{*}\cdot {{\rm{f}}}_{{\rm{n}},+}\left(\overrightarrow{{\rm{x}}}\right) \end{split} $$ (6) and can be generated using the estimated wavefield
$ {{{\rm{f}}}^{{\rm{M}}}}_{{\rm{n}}}\left(\overrightarrow{{\rm{x}}}\right) $ at any shear$ \overrightarrow{{S}_{n}} $ .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.
Fig. 5 shows the steps of the wavefield reconstruction algorithm. After an initial guess, the algorithm continuously applies the socalled object and measurement constraints alternatingly until the loss function reaches the desired level. The application of those constraints is discussed in more detail below.

This step of the algorithm uses only the initial guess
$ {\rm{f}}\left(\overrightarrow{{\rm{x}}}\right) $ in the first iteration to estimate$ {{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},\pm }}\left(\overrightarrow{{\rm{x}}}\right) $ . An initial guess is assumed for$ {\rm{f}}\left(\overrightarrow{{\rm{x}}}\right) $ (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$ {\rm{f}}\left(\overrightarrow{{\rm{x}}}\right) $ , the shearing is done in the Fourier domain using the following equations.$$ {\rm{F}}\left(\overrightarrow{{\rm{k}}}\right)= {\cal{F}}\left\{{\rm{f}}\left(\overrightarrow{{\rm{x}}}\right)\right\} $$ (7) $$ {{{\rm{f}}^{{\rm{M}}}_{{\rm{n}},\pm }}}\left(\overrightarrow{{\rm{x}}}\right)={\cal{F}}^{1}\left\{{\rm{F}}\left(\overrightarrow{{\rm{k}}}\right).{\rm{exp}}\left(\mp 2{\rm{\pi}}\cdot {\rm{i}}\cdot \overrightarrow{{\rm{k}}}\cdot \overrightarrow{{{\rm{d}}{\rm{s}}}_{{\rm{n}}}}\right)\right\} $$ (8) Assuming that the initial guess wavefield is from the measurement constraint domain, the initial guess and its sheared components are indicated as
$ {{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},\pm }}\left(\overrightarrow{{\rm{x}}}\right) $ . A wavefield with constant phase and unit amplitude is used as an initial guess for the reconstructions done in this paper. 
The object constraints are applied by transforming the sheared components of the estimated wavefield to the Fourier domain and compensating the shear for each
$ {{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},\pm }}\left(\overrightarrow{{\rm{x}}}\right) $ to obtain$ {\rm{f}}\left(\overrightarrow{{\rm{x}}}\right) $ using the Fourier Shift theorem. All calculated values for$ {\rm{F}}\left(\overrightarrow{{\rm{k}}}\right) $ are averaged to obtain a more accurate result.$$ {{{\rm{F}}}^{{\rm{M}}}_{{\rm{n}},\pm }}\left(\overrightarrow{{\rm{k}}}\right)= {\cal{F}}\left\{{{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},\pm }}\left(\overrightarrow{{\rm{x}}}\right)\right\} $$ (9) $$\begin{split} {{\rm{F}}}_{{\rm{est}}}\left(\overrightarrow{{\rm{k}}}\right)=& \frac{1}{2{\rm{S}}}\sum\limits_{{\rm{n}}=1}^{{\rm{S}}}\Big({{{\rm{F}}}^{{\rm{M}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{k}}}\right)\cdot {\rm{exp}}\left(2{\rm{\pi}}\cdot {\rm{i}}\cdot \overrightarrow{{\rm{k}}}\cdot \overrightarrow{{{\rm{d}}{\rm{s}}}_{{\rm{n}}}}\right) \\& +{{\rm{F}}^{{\rm{M}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{k}}}\right)\cdot {\rm{exp}}\left(2{\rm{\pi}}\cdot {\rm{i}}\cdot \overrightarrow{{\rm{k}}}\cdot \overrightarrow{{{\rm{d}}{\rm{s}}}_{{\rm{n}}}}\right)\Big) \end{split}$$ (10) Finally,
$ {{\rm{F}}}_{{\rm{e}}{\rm{s}}{\rm{t}}}\left(\overrightarrow{{\rm{k}}}\right) $ 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$ {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},\pm }}\left(\overrightarrow{{\rm{x}}}\right) $ as$$ {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},\pm }}\left(\overrightarrow{{\rm{x}}}\right)={\cal{F}}^{1}\left\{{{\rm{F}}}_{{\rm{e}}{\rm{s}}{\rm{t}}}\left(\overrightarrow{{\rm{k}}}\right)\cdot {\rm{exp}}\left(\mp 2{\rm{\pi}}\cdot {\rm{i}}\cdot \overrightarrow{{\rm{k}}}\cdot \overrightarrow{{{\rm{d}}{\rm{s}}}_{{\rm{n}}}}\right)\right\} $$ (11) 
The measurement constraint equations use the estimated wavefield from the object constraint
$ {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},\pm }}\left(\overrightarrow{{\rm{x}}}\right) $ and the measured complex crossterm$ {{\rm{C}}}_{{\rm{n}}}\left(\overrightarrow{{\rm{x}}}\right) $ (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.$$ {{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)= {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)+\frac{\left({{\rm{C}}}_{{\rm{n}}}\left(\overrightarrow{{\rm{x}}}\right){{{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right)}^{*}\cdot {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)\right)\cdot {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right)}{{\left{{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)\right}^{2}+{\left{{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right)\right}^{2}} $$ (12) $$ {{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right)= {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right)+\frac{{\left({{\rm{C}}}_{{\rm{n}}}\left(\overrightarrow{{\rm{x}}}\right){{{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right)}^{*}\cdot {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)\right)}^{*}\cdot {{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)}{{\left{{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)\right}^{2}+{\left{{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right)\right}^{2}} $$ (13) 
The algorithm uses two sets of constraint equations, object constraint Eq. 10, and measurement constraint Eqs. 12, 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.

The performance of the algorithm was evaluated by simulations using MATLAB. Preliminary simulations were made for a resolution of 256 × 256 pixels. For this simulation, a point source was generated with a wavelength of 600 nm at a zdistance 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 phaseshifted intensities were generated for each shear position. The shears generated for this simulation are shown in Fig. 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 8bit discretization has been applied to the intensities.
Fig. 6 (Simulation) Results of reconstruction with simulation parameters: FoV 256 × 256 pixels, Pixel size 3.45 μm,Wavelength 0.6 μm, Aperture size 90 pixels, Z distance 87.5 mm. a Shear selection. b Ideal phase map. c Reconstructed phase map. d Error map between ideal and reconstructed data.
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 Fig. 6b, c, respectively. The corresponding phase error of the reconstructed wavefront is shown in Fig. 6d. Further investigation of the error map showed the presence of highfrequency artifacts, as highlighted in Fig. 7.
Fig. 7 High frequency artifact from reconstruction algorithm. a Error map zoomed. b Error map profile.
These highfrequency 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 lowpass Fourier filter to mitigate the formation of these artifacts.
Fig. 8 shows error maps generated by mapping the difference between the ideal and reconstructed wavefield before and after applying the filter. A crosssection of the error profile before and after filtering is shown in Fig. 8c, d, respectively. This simple procedure reduces the magnitude of the RMS error by one order of magnitude. Because a Fourier lowpass filter has been employed with a corner frequency of
$ {k}_{c}=\left({4}/{5}\right)\cdot \left({1}/{\rm{pixel\;size}}\right) $ , 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. 
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. Fig. 9a shows the evolution of loss function during reconstructions done in Fig. 8 before the filter was applied.
$$ {\rm{L}} =\sum\limits_{{\rm{n}} = 1}^{{\rm{S}}}\left({\left{{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right){{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},+}}\left(\overrightarrow{{\rm{x}}}\right)\right}^{2}+{\left{{{\rm{f}}}^{{\rm{M}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right){{{\rm{f}}}^{{\rm{O}}}_{{\rm{n}},}}\left(\overrightarrow{{\rm{x}}}\right)\right}^{2}\right) $$ (14) Fig. 9a shows that the loss function converges after 30 iterations of the reconstruction without a Fourier filter. Given the filter’s effectiveness in Fig. 8, we have investigated the use of the filter within the reconstruction algorithm. The sudden surge at the start of the loss function in Fig. 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 inputoutput 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. Fig. 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 highfrequency 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 Ref. 16. For better clarity we summarize the differences to the original algorithm from Ref. 16.
● The data recorded from experiments in Ref. 16 are in the Fourier domain. In this paper, all recorded data from experiments are in the spatial domain.
● Eq. 10 is the modified version of the object constraint equation. The original equation (Eq. 10 in Ref. 16) was defined in the spatial domain. In this paper, this equation is defined in the Fourier domain.
● Eqs. 12,13 are the modified versions of the measurement constraint equations. The original equations (Eq. 15 in Ref. 16) were defined in the Fourier domain. In this paper, these equations are defined in the spatial domain.
● The algorithm in Ref. 16 uses two sets of equations in measurement constraint (Eqs. 15,16 in Ref. 16), roughreconstruction equations for initial convergence (Eq. 16 in Ref. 16) and finereconstruction equations for final refinement (Eq. 15 in Ref. 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 highfrequency artifacts.
● Multiresolution technique^{10} is incorporated into this algorithm to mitigate the effects of initial guess and improve convergence.

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 Fig. 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 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 Fig. 10.
Fig. 10 Comparison of reconstruction results between different configurations. a Shear selection for test 1 in cartesian coordinate configaration. b Shear selection for test 2 in cartesian coordinate configuration. c Shear selection for polar coordinate configuration. d Reconstructed phase map from test 1. e Reconstructed phase map from test 2. f Reconstructed phase map from polar coordinate configuration. g Phase error map from reconstruction using test 1 parameters. h Phase error map from reconstruction using test 2 parameters. i Phase error map from reconstruction using polar coordinate configuration.
Each simulation was performed with 24 shears, where the selections are shown in Figs. 10a−c. The selections analyze two possible xyshearing configurations shown in Fig. 10a, b (see Fig. 2 for configuration) and the fixed shear with variable shear axis shown in Fig. 10c (see Fig. 3 for configuration).
In test 1 (Fig. 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 Fig. 10d, g. When the selection of shears was scattered across one quadrant (Fig. 10b) as in test 2, the reconstruction algorithm had sufficient data to reconstruct most of the pointsource wavefield, as shown in Fig. 10e. However, the error map in Fig. 10h shows a curvature, indicating the possibility of some missing frequencies. This slowly varying error is critical as it produces a socalled 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 Fig. 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 Fig. 10f. The error map shown in Fig. 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.

To better understand the results and characteristics features in the error maps of Fig. 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 nonsymmetric 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
$$ T\left({k}_{x},{k}_{y}\right)=\sum _{n=1}^{S}{{\rm{sin}}}^{2}\left[2\pi \overrightarrow{k}.\overrightarrow{d{s}_{n}}\right] $$ (15) where
$ \overrightarrow{k}=\left[{k}_{x},{k}_{y}\right] $ is the frequency coordinate and$ \overrightarrow{ds}=[d{s}_{x}, $ $ d{s}_{y}] $ 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 Fig. 11.Fig. 11 Transfer function analysis (Frequency information density plots with points below threshold 2): a Test lin cartesian coordinate configuration (34 datapoints below threshold). b Test 2 in cartesian coordinate configuration (21 datapoints below threshold). c Polar coordinate configuration (1 datapoint below threshold). d Test 1 plot in cartesian coordinate configuration zoomed near small frequency regions. e Test 2 plot in cartesian coordinate configuration zoomed near small frequency regions. f Polar coordinate configuration plot zoomed near small frequency regions.
Figs. 11a−c are the frequency domain information density maps for the three shear sets described in Fig. 10. Figs. 11d−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. Fig. 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 45degree angle will have poor reconstruction in the estimated wavefield at longer spatial wavelengths. This conclusion is in agreement with Fig. 10g. A similar behavior is observed in Fig. 11b, e, 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 Fig. 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 Fig. 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.

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 Fig. 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. Fig. 12b shows which areas of the input beam provide phase information to the reconstruction algorithm (regions indicated in yellow).
Fig. 12 Demonstration of data recorded by each shear: a Beam 1 and 2 overlapping. b Summation of overlapped regions indicating regions that have data for reconstruction (yellow regions).
This idea is further extended to multishear interferometry in Fig. 13. Fig. 13a shows the overlap regions from an object beam with a diameter of 600 pixels sheared by 175 pixels along the xdirection. Fig. 13b shows the overlap regions if a second shear is added to the measurement in the orthogonal direction by the same amount. In Fig. 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 Fig. 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). Fig. 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 Fig. 13d has the highest potential to successfully reconstruct the wavefield due to maximum overlap across the whole area.
Fig. 13 Spatial information density distribution for different shear patterns on 1024 × 1024 pixel maps with aperture diameter of 600 pixels. a Single shear of 175 pixels. b Two orthogonal and equal shears of 175 pixels. c 24 shear in polar configuration with a fixed shear of 175 pixels. d 24 shear in polar configuration with a fixed shear of 80.5 pixels.
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 Fig. 10 and the overlap distribution is mapped in Fig. 14. Fig. 14a, b show fewer overlaps along the boundary at 45degree angle. The number of effective measurements drops from 24 to 13. On the contrary, Fig. 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.

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. MonteCarlo 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 Fig. 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 Fig. 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. Figs. 15b−d show the measurement of a USAF target where the Monte Carlo simulation used a shear error of 1, 3, and 5pixels standard deviation, respectively.

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 autocorrelation of the recorded intensity maps. To increase the robustness of this estimation, we have removed the sinusoidal fringe modulation by averaging all four measured phaseshifted intensity patterns. One examplary measurement is shown in Fig. 16a. The secondorder derivative of the autocorrelation map has been computed along one direction to further enhance the peak in the autocorrelation map. The resulting spikes for the case in Fig. 16a are shown in Fig. 16b, and they are prominently identifiable. The secondorder 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.

When working with left and right circular polarized waves, it is convenient to employ pixelated polarization cameras to record four phaseshifted intensities simultaneously. The polarization camera enables singleshot phaseshifting during the measurement process for each shear. The singleshot phaseshifting 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 superpixel 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 leftcircular and rightcircular polarization, it is possible to obtain simultaneously four interferograms at four different phase shifts. The orientation of this pixelated array is shown in Fig. 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 × 1024 pixels.
Fig. 17 Schematic of the polarization camera sensor and intensity maps demonstrating the phase shifts.
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 1/2 superpixels.
Variable shearing holography with applications to phase imaging and metrology
 Light: Advanced Manufacturing , Article number: 16 (2022)
 Received: 16 September 2021
 Revised: 08 February 2022
 Accepted: 15 February 2022 Published online: 25 March 2022
doi: https://doi.org/10.37188/lam.2022.016
Abstract: We report variable shear interferometers employing liquidcrystalbased geometricphase (GP) gratings. Conventional gratingbased shear interferometers require two lateral shifts of the gratings to enable phaseshifting 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 instantaneousphase 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.
Research Summary
Variable Shearing Holography with Geometric Phase Gratings:
Different variable shearing interferometers employing liquidcrystalbased geometricphase (GP) gratings are reported. When combined with pixelated polarization cameras, the GP gratings enable instantaneous phaseshifting capabilities and more robust and compact designs. A comparative analysis of different shearing configurations is reported, which shows that using a single fixed shear but with a variable shear axis provides high accuracy while requiring the control of only one mechanical axis. We also show that, the ability to apply a given shear accurately affects the uncertainty and the resolution of the reconstructed wavefield.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article′s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article′s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.