Optimization based evaluation of grating interferometric phase stepping series and analysis of mechanical setup instabilities

The diffraction contrast modalities accessible by X-ray grating interferometers are not imaged directly but have to be inferred from sine like signal variations occurring in a series of images acquired at varying relative positions of the interferometer's gratings. The absolute spatial translations involved in the acquisition of these phase stepping series usually lie in the range of only a few hundred nanometers, wherefore positioning errors as small as 10nm will already translate into signal uncertainties of one to ten percent in the final images if not accounted for. Classically, the relative grating positions in the phase stepping series are considered input parameters to the analysis and are, for the Fast Fourier Transform that is typically employed, required to be equidistantly distributed over multiples of the gratings' period. More general optimization based approaches with relaxed requirements on the sampling positions are often deemed too time consuming in contrast. In the following, a fast converging optimization scheme is presented simultaneously determining the phase stepping curves' parameters as well as the actually performed motions of the stepped grating, including also erroneous rotational motions which are commonly neglected. While the correction of solely the translational errors along the stepping direction is found to be sufficient with regard to the reduction of image artifacts, the possibility to also detect minute rotations about all axes proves to be a valuable tool for system calibration and monitoring. The simplicity of the provided algorithm, in particular when only considering translational errors, makes it well suitable as a standard evaluation procedure also for large image series.


Introduction
X-ray grating interferometry [1,2] facilitates access to new contrast modalities in laboratory X-ray imaging setups and has by now been implemented by many research groups after the seminal publication by Pfeiffer et al. in 2006 [3]. The additional information on X-ray refraction ("differential phase contrast") and ultra small angle scattering ("darkfield contrast") properties of a sample that can be obtained promises both increased sensitivity to subtle material variations as well as insights into the samples' substructure below the spatial resolution of the acquired images.
In contrast to classic X-ray imaging, the absorption, differential phase and darkfield contrasts are not imaged directly but are encoded in sinusoidal intensity variations arising at each detector pixel when shifting the interferometer's gratings relative to each other perpendicular to the beam path and the grating bars. A crucial step in the generation of respective absorption, phase and darkfield images therefore is the analysis of the commonly acquired phase stepping series, which shall be the subject of the present article. Respective examples are shown in Figures 1 and 2.
In principle, the images within such phase stepping series are sampled at about five to ten different relative grating positions equidistantly distributed over multiples of the gratings' period such that the expected sinusoids for each detector pixel may be characterized by standard Fourier decomposition.  Figure 1: Example for a grating interferometric phase stepping series. The intensity variation throughout the series is visually most perceivable in the center. The spatial Moiré fringes are due to imperfectly matched gratings and will translate to reference offset phases of the sinusoid curves found at each detector pixel (cf. Fig. 2). The bottom row shows a corresponding phase stepping curve for the pixel marked orange in the above image series. The sampling positions are subject to an unknown error.
The zeroth order term represents the mean transmitted intensity (as in classic X-ray imaging), while the first order terms encode phase shift and amplitude of the sinusoid. The ratio of amplitude and mean (generally referred to as "visibility") is here related to scattering and provides the darkfield contrast. Higher order terms correspond to deviations from the sinusoid model mainly due to the actual grating profiles and are usually not considered. Given typical grating periods in the range of two to ten micrometers, the actually performed spatial translations lie in the range of 200 to 2000 nanometers. Particular for the smaller gratings, positioning errors as small as 10 to 20 nm imply relative phase errors in the range of five to ten percent, causing uncertainties in the derived quantities in the same order of magnitude. The propagation of noise within the sampling positions onto the extracted signals has e.g. been studied by Revol et al. [4] and first results for the determination of the actual sampling positions from the available image series were shown by Seifert et al. [5] using methods by Vargas et al. [6] from the context of visible light interferometry. Otherwise the problem seems to have not been given much consideration so far.
The present article proposes a simple iterative optimization algorithm both for the fitting of irregularly sampled sinusoids and in particular also for the determination of the actual sampling positions. The use of only basic mathematical operations eases straightforward implementations on arbitrary platforms. Besides uncertainties in the lateral stepping motion, the remaining mechanical degrees of freedom (magnification/expansion and rotations) are also considered. The proposed techniques will be demonstrated on an exemplary data set.

Methods
The task of sinusoid fitting with imprecisely known sampling locations will be partitioned into two separate optimization problems considering either only the sinusoid parameters or only the sampling positions while temporarily fixing the respective other set of parameters. Alternating both optimization tasks will minimize the joint objective function in few iterations. : From left to right: transmission (sinusoid mean), visibility (ratio of sinusoid amplitude and mean) and phase images derived from phase stepping series as shown in Fig. 1. The first rows show acquisitions with and without sample, while the last row shows the sample images (center row) normalized with respect to the empty beam images (first row). Positioning errors in the phase stepping procedure cause the Moiré pattern of the reference phase image to translate into the final results.

Sinusoid Fitting
A simple iterative forward-backprojection algorithm converging to the least squares solution can be derived when representing the sinusoid model as a linear combination of basis functions: with the following identities: The components o, a c , a s may be determined by means of the following iterative scheme (introducing the sample index i and the superscript iteration index k): where the factors of 1/N and 2/N account for the normalization of the respective basis functions (with N being the amount of samples (φ i , y i ) enumerated by i). The scheme reduces to classic Fourier analysis for the case of the abscissas φ i being equidistantly distributed over multiples of 2π and converges within the first iteration in that case. As the update terms to o, a c and a s are proportional to the respective derivatives of the 2 error i (o (k) + a (k) s sin φ i − y i ) 2 , the fixpoint of the iteration will be the least squares fit also in all other cases.
For a stopping criterion, the relative error reduction may be tracked. It is typically found to fall below 0.1% within 10 to 20 iterations given only slightly noisy data (noise sigma three orders of magnitude smaller than sinusoid amplitude) and within less than 10 iterations for most practical cases. For the special case of equidistributed φ i on multiples of 2π, it will immediately drop to 0 after the first iteration. In practice, a fixed amount of iterations in the range of five to fifteen will therefore be adequate as stopping criterion as well.

Phase Step Optimization
An underlying assumption of the previously described least squares fitting procedure is the certainty of the abscissas, i.e. the set of phases φ i at which the ordinates y i have been sampled. As the sampling positions are themselves subject to experimental uncertainties (arising from the mechanical precision 0 /2 3 /2 2 phase noisy data sinusoid fit optimized sampling positions Figure 3: Optimization of the actual (in contrast to the intended) sampling positions by Eq. 13 with respect to a previous sinusoid fit based on the temporary assumption that deviations from the sinusoid model are mainly due to errors on the sampling positions rather than the sampled ordinate values. Averaging of respective phase deviations found for large sets of measurements will finally allow the differentiation of systematic deviations from statistical noise (see also Fig. 6).
of the involved actuators), a further optimization step will be introduced that minimizes the least squares error of the sinusoid fit over the sampling positions φ i . While this procedure obviously results in overfitting when considering only a single phase stepping curve (PSC), it becomes a well-defined error minimization problem when regarding large sets of PSCs sharing the same φ i . In other words, an approach to the minimization of the objective function shall be considered, where j indexes detector pixels.
In order to derive an optimization procedure for the sampling positions, first the fictive case of a perfectly sinusoid PSC with negligible statistical error on the ordinate (the sampled values) shall be considered. Ignoring for now the fact that least squares fits commonly assume only the ordinates to be affected by noise, a least squares fit shall be used to preliminarily determine the parameters of the sinusoid described by the observed data. Assuming then that inconsistencies of the observed data with the model are due to errors on the sampling locations, deviations from their intended positions are given by the data points' lateral distances from the sinusoid curve (cf. Fig. 3). Finally, the actual systematic deviations of the sampling locations can be found by averaging over the respective results for a large set of PSCs sampled simultaneously. This information can be fed back into the original sinusoid fit, which then again allows the refinement of the current estimate of the true sampling positions, finally resulting in an iterative procedure alternatingly optimizing the sinusoid parameters and the actual sampling locations (cf. Algorithm 1).

Determination of individual phase deviations
Starting with an initial sinusoid fit the error shall be minimized over deviations ∆φ i to φ i while keeping the sinusoid model parameters fixed: Equation 7 is solved when the derivative of the objective function with respect to ∆φ i vanishes: where the last step is a first order Taylor expansion with respect to ∆φ i . This directly leads to the following expression for ∆φ i : where the earlier condition a cos(φ i +∆φ i −φ 0 ) = 0 is approximated to be satisfied when a cos(φ i −φ 0 ) = 0. The restriction to cases with a cos(φ i − φ 0 ) = 0 can be intuitively understood when recalling that cos(φ i − φ 0 ) = 0 implies a maximum or minimum of the sinusoid and a = 0 means that it is constant (φ i independent), in both of which cases there is no sensible choice for ∆φ i = 0. The constraint on the result, ∆φ i π, can simply be taken into account by means of a limiting function such as which provides a linear mapping with slope 1 for ∆φ i m and is bounded at ±m. The choice of m in this case depends on the validity range of the linear approximation of sin(φ i + ∆φ i − φ 0 ) with respect to ∆φ i about φ i − φ 0 , which obviously depends on the magnitude of the curvature of the sinusoid at this point as illustrated in Figure 4. The latter may be accounted for by which reaches its maximum m φi−φ0 = m 0 at φ i − φ 0 = 0 (where sin(φ i − φ 0 ) is actually linear) and smoothly drops to 0 for cos(φ i − φ 0 ) = 0, in which case both the sine and its curvature are maximal and ∆φ i shall remain 0. The actual choice of m 0 can now be based on the validity range of the linear approximation about φ i − φ 0 = 0. The upper bound to m φi−φ0 and thus to m 0 is defined by the range over which sin(x − x 0 ), x ∈ [−m 0 cos 2 (x 0 ), +m 0 cos 2 (x 0 )] is actually invertible (cf. Fig. 4), i.e.: For m 0 = 1.38, the linear approximation used in Eq. 8 deviates by up to 40%. The deviation is reduced to 20% or 5% for m 0 = 1 and m 0 = 1 2 respectively. Combining the above results and introducing for completeness also the detector pixel index j, the following expression for ∆φ ji reducing (and, for sufficiently small ∆φ ji , minimizing) the 2 error (o j + a j sin(φ i + ∆φ ji − φ 0,j ) − y ji ) 2 can be given, choosing m 0 = 1 2 : Figure 3 shows an example of this approximate least squares solution to ∆φ ji . The maximum meaningful range is thus largest for points furthest from these boundaries and reduces to zero exactly at the turning points. Equation 11 approximates this phase dependence of the validity range with a cos 2 function, and Eq. 12 defines the maximum amplitude admissible in order to indeed never exceed turning points.

Noise weighted average of phase deviations
Now that an expression has been derived for the deviations ∆φ ji optimizing the abscissa values for a single PSC given a previous sinusoid fit, the respective results for many PSCs (indexed by j) sharing the same sampling locations may be averaged: using weights w ji factoring in the relative certainty and relevance of the different ∆φ ji . An appropriate choice is where cos 2 (φ i − φ 0,j ) weights the slope dependent error propagation from noisy measurements y ji to ∆φ ji based on the derivative of the sinusoid model at φ i and a 2 j weights the contribution of a particular PSC j to the accumulated 2 error. These considerations lead to where the last step uses the relation α softlimit(x, m) = softlimit(αx, αm) for α ≥ 0. Finally, the above derivations can be combined to an iterative optimization algorithm reducing the accumulated least square error of multiple sinusoid fits (indexed by j) to data points y ji over shared abscissa values φ i as defined by Equation 5. A pseudo code representation is given in Alg. 1, further introducing the relaxation parameter λ k ∈ (0; 1] that may be chosen < 1 in order to damp the adaptions to φ : end for.

Inhomogeneous sampling phase deviations
Up to now, it has been assumed that deviations from the intended phase stepping positions are due to purely translational uncertainties in the relative motion of the involved gratings, resulting in offsets ∆φ i of the actual from the intended sampling phases that are homogeneous throughout the whole detection area. When also considering relative grating period changes (e.g. due to either thermal expansion or motion induced changes in magnification) and rotary motions of the interferometer's gratings relative to each other (e.g. due to backlashes within the mechanical actuators), the effective sampling phases at each phase step may exhibit gradients over the detection area. Given the small grating periods (micrometer scale) compared to the total extents of the gratings (centimeter scale), both tilts in the sub-microrad range and relative period changes in the range of 10 −7 will already manifest themselves in observable gradients. The corresponding optimization problem regarding also gradients is, analog to Equation 5, given by: where the spatial dependence of the phase deviations ∆φ i has been accounted for by a functional dependence on the detector pixel index j. Given the expected gradients (as illustrated in Figure 5), ∆φ i (j) has the following form: h φ i quantifying the respective gradients in horizontal and vertical direction as well as the mixed term and the curvature in horizontal direction, and h and v being spatial detector pixel indices related to the linear pixel index j through the amount N h of pixels within one detector row: The constant offsets h 0 and v 0 characterize the detector center. The optimization of the extended objective function in Equation 17 can be performed analog to that of Eq. 5 with the only difference lying in the evaluation of the spatial phase difference maps ∆φ ji translation const. offset defined by Eq. 13 (an example of which is shown in Figure 6). The weighted average derived in the previous section in order to determine the homogeneous offset ∆φ i can be extended to a generalized linear least squares fit of the model ∆φ i (j) = ∆φ i (h(j), v(j)) defined by Eq. 18 to the local estimates ∆φ ij (Eq. 13, Fig. 6), also taking the weights defined by Eq. 15 into account. Said procedure is stated more formally in Algorithm 2.
Basic geometric considerations neglecting higher order interrelations of the considered effects (e.g. rotation and effective period change) result in the following relations between the observable parameters h φ i and relative translatory and rotatory motions of the interferometer's gratings. In order to relate various magnification changes to spatial motions based on the intercept theorem, an assumption has to be made as to which of the gratings actually moved. Here, the grating that is mounted on the linear phase stepping actuator is assumed to be the cause of all relative motions of both gratings also including tilts and rotations. The "source-grating distance" in the following equations will thus refer to the stepped grating. ∆φ i quantifies the translational error analog to Section 2.2.2. In contrast to the previous section, the present model distinguishes between homogeneous phase deviations induced by translation and the mean component induced by the ∇ 2 h φ i (h − h 0 ) 2 term in the case of non-vanishing curvature of the spatial phase deviation.
The vertical gradient parameter ∇ v φ i is related to a relative rotation η of both gratings about the optical axis: where the "effective grating period" refers to the projected period length at the location of the detector, which should be identical for both interferometer gratings (not considering the optional additional coherence grating close to the X-ray source). The horizontal gradient parameter ∇ h φ i is related to a relative mismatch in effective grating periods of the gratings either due to relative distance changes along the optical axis or due to actual expansions (e.g. thermally induced): relative period mismatch = effective period difference effective grating period = ∇ h φ i 2π effective grating period detector pixel pitch .
Algorithm 2 Simultaneous least squares optimization of abscissa values φ ji and sinusoid fits to ordinate samples y ji belonging to independent curves j sampled at positions φ ji = φ i (j) with φ i (j) being a slowly varying polynomial with respect to the spatial coordinates h(j) and v(j) accounting for the expected effects due to translations, magnification and rotations of an interferometer's gratings. The procedure reduces to Algorithm 1 when considering only the zeroth order term of φ i (j).
1: φ i , y ji : input data 2: m 0 ← 1 2 upper limit to ∆φ ji , m 0 ∈ (0; 1.38] When assuming relative grating period mismatches to be caused by changes in magnification due to translations of one of the gratings along the optical axis, the following relation applies to first order: translation distance = (relative period mismatch) (source-grating distance) 2 source-detector distance = ∇ h φ i 2π effective grating period detector pixel pitch (source-grating distance) 2 source-detector distance .
The change ∇ hv φ i of the horizontal gradient throughout the vertical direction corresponds to a relative change in magnification from top to bottom, e.g. due to a tilt θ of one of the gratings about the horizontal axis. Using the above relation between magnification changes and spatial displacements, the tilt θ about the horizontal axis is related to ∇ hv φ i approximately by tan θ = ∇ hv φ i 2π effective grating period detector pixel pitch (detector pixel pitch) source-grating distance source-detector distance (effective grating period)(source-grating distance) (detector pixel pitch) 2 .

(23)
A non-vanishing curvature ∇ 2 h φ i arises in case of a rotary motion about the vertical axis (slant) and is analogously related to the slant angle ϕ to first order by (effective grating period)(source-grating distance) (detector pixel pitch) 2 .  Figure 6: Example of phase differences (left) obtained from Equation 13 (cf. also Fig. 3) for the first frame of the phase stepping series shown in Fig. 1 with corresponding importance weights (right) as defined by Eq. 15. White regions (0 on the colorbar) in the left hand side ∆φ map correspond to samples close to or exactly on turning points of the fitted sinusoids, where phase differences cannot be effectively determined. These areas get little weighting in the determination of the average phase deviation as can be seen by the corresponding dark fringes in side weighting map on the right.

Experiment and Results
Phase stepping series of 15 images sampled at varying relative grating shifts uniformly distributed over three grating periods have been acquired both with and without sample in the beam path. Figure 1 shows the first five frames of the empty beam series. The resulting phase stepping curves at each pixel (indexed by j) of the detector have been evaluated using a least squares fit to a sinusoid model parameterized by mean o j , amplitude a j and phase offset φ 0,j under the initial assumption of perfectly stepped gratings. These preliminary results are shown in Figure 2 and correspond to those obtained by classic Fourier analysis of the phase stepping curves. Deviations of the sampling positions from the intended ones are then determined based on systematic deviations of the sampled data from the fitted sinusoids by means of Eq. 16 for all 15 frames of the phase stepping series. Figure 6 shows an exemplary result for the first frame of the series. Iterating the sinusoid fits and the corrections to the sampling positions in order to reduce the overall least squares error (cf. Equation 5) by means of Algorithm 1, the sampling positions' deviations are found as shown in Figure 7. Figure 8 shows the reduction of Moiré modulated systematic errors in the final results, i.e. the transmission, visibility and differential phase images. The root mean square error is reduced by almost a factor of two in the present example and is already close to convergence after the first iteration as can be seen in Fig. 9.
In addition to the mean deviations of the phase steps from the intended positions, spatial gradients throughout the detection area have been also considered (cf. Alg. 2). Figure 10 shows the respective differential deviations from the intended phase steps between the first 9 frames of the phase stepping series, normalized to the nominal homogeneous phase stepping increment of 2π/5. The mean contributions of each component are listed in Table 1. While the homogeneous error of 0.1 rad ranges within 10% of the nominal step size (or 2% of the grating period), the remaining effects are two to three orders of magnitude smaller. The root mean square error of the sinusoid fits for the whole phase stepping series  is reduced by 0.1% relative to the optimization considering only homogeneous phase step deviations as shown in Figure 9. Consequently, the derived images (not shown) are visually equivalent to those obtained previously (cf. Fig. 8). Figure 11 shows variations in the relative alignment of the gratings derived from the inhomogeneous phase stepping analysis by means of Equations 20-24. Besides deviations from the nominal linear motion of the gratings, minute rotations as well as subtle changes in relative magnification can also be detected. The correlation between grating rotations and translational errors visible in the left hand side graph in Figure 11 indicate rotations about an off-center pivot point about 10 −1 m below the grating center, which is consistent with the actual placement of the phase stepping actuator in the experimental setup. Although the analog correlation between tilts about the horizontal axis and translation (along the optical axis) induced variations in magnification is much less pronounced (Fig. 11, right hand side), the mean trend and magnitude are also consistent with the assumption of a pivot point below the field of view. However, the observed magnitude (10 −7 ) of the relative mismatch of the effective grating periods is as well explicable by temperature variations in the order of magnitude of 10 −1 K given a thermal expansion coefficient in the order of magnitude of 10 −6 K −1 for the typical wafer materials silicon and graphite. Finally, the phase stepping inhomogeneities further suggest rotational motions about the vertical axis on the microrad scale (also shown in Fig. 11).
As a crude error assessment, the standard error of the mean phase deviation can be estimated from the sinusoid fits' root mean square error: contributing detector pixels sinusoid fit RMSE mean sinusoid amplitude ≈ 2 detector pixels sinusoid fit RMSE mean sinusoid amplitude .
The latter corresponds for the present data set to about 6% of the mean observed sinusoid amplitude, which directly translates to 6×10 −2 rad on the abscissa. Given the amount of detector pixels contributing to the least squares fits of ∆φ i (j) within each frame of the phase stepping series, a standard error in the order of magnitude of σ mean ≈ 10 −4 rad results. This implies that, according to the results given in Table 1 Figure 10: Relative deviations of the actual phase steps from the intended step width of 2π 5 between the first nine frames of the phase stepping series ( 5 2π (∆φ i+1 (j) − ∆φ i (j) − 2π 5 )). The variations ∆φ i (j) have been determined by optimization of Eq. 17 assuming the spatial dependence defined by Eq. 18 (see also Algorithm 2)  Figure 11: Quantitative results derivable from the inhomogeneities in the phase stepping deviations ∆φ i (j) (Eq. 18). On the left, the change in tilt angle about the optical axis from frame to frame within the phase stepping series is shown along with the accompanying linear motion error. Rotation correlated translations indicate an off-center pivot point. On the right hand side, the relative grating scaling error is shown along with the found tilt and slant about the horizontal and vertical axis. These quantities represent deviations from the mean grating alignment throughout the phase stepping series. The tilt and slant angles range close to the expected noise level (cf. Table 1 and Eq. 26).

Discussion and Conclusion
A fast converging iterative algorithm for the joint optimization of both the sinusoid model parameters and the actual sampling locations for the evaluation of grating interferometric phase stepping series has been proposed. Within the optimization procedure, spatially varying phase stepping increments due to further mechanical degrees of freedom besides the intended translatory stepping motion can also be accounted for. Mean phase stepping errors of up to 10% of the nominal step width have been found for a typical data set, and their correction results both in a considerable reduction of the overall root mean square error by almost a factor of two as well as a significant visual improvement of the final images. Although higher order effects are observable, their contribution was found to be two to three orders of magnitude smaller than that of the mean stepping error, and their correction thus did not contribute to further improvements in visual image quality in the present case. However, the higher order deviations allow the detection of minute motions of the gratings and thus provide a valuable tool for the monitoring and debugging of experimental setups. First order approximations for the relations between spatial phase variations and mechanical degrees of freedom of the moved grating have been given. For the present data set, linear motion errors up to 25 nm as well as rotational motions on the microrad scale have been inferred from the phase stepping series. While the tilt and slant angles about the horizontal and vertical axes respectively have been found to range close to the expected noise level and should rather be interpreted as upper limits to actual motions, magnification changes in the range of 10 −7 and sub-microrad rotations about the optical axis were well detectable. The expected correlations between rotation and translation due to an off-center pivot point further support the plausibility of the results. Especially the crosstalk between sub-microrad rotations and effective translations indicates that noticeable phase stepping errors will be almost inevitable even for very carefully designed experiments, wherefore an optimization based evaluation of the phase stepping series as proposed in Algorithm 1 is generally advisable. With processing speeds in the range of 0.1 s per phase stepping series, it is well suitable as a standard processing method also for large image series.