Optimization Based Evaluation of Grating Interferometric Phase Stepping Series and Analysis of Mechanical Setup Instabilities

Abstract

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 10 nm will already translate into signal uncertainties of 1–10% 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. 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.

1. 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 is the subject of the present article. Respective examples are shown in Figure 1 and Figure 2.

In principle, the images within such phase stepping series are sampled at about 5–10 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 [1,2,3]. 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 2–10 micrometers, the actually performed spatial translations lie in the range of 200–2000 nanometers. Particular for the smaller gratings, positioning errors as small as 10–20 nm imply relative phase errors in the range of 5–10 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 based on a decomposition of the images into three basis images (to be found by means of principal component analysis and subsequent optimization) corresponding to the transmission image and the amplitude modulated sine and cosine of the phase image. Kaeppler et al. [7] optimized the phase stepping positions by means of minimizing an objective function penalizing variations within the visibility and differential phase images that arise from phase stepping errors. For the case of simultaneous grating stepping and object rotation, von Teuffenbach et al. [8] recently reported the inference of stepping errors as part of a maximum likelihood tomographic reconstruction procedure proposed by Ritter et al. [9]. Finally, the authors just became aware of the parallel work by De Marco et al. [10] who investigated the origination of the particular image artifacts caused by phase stepping jitter in order to derive a postprocessing artifact reduction technique.

Especially with regard to the quantitative analysis of darkfield signals [11,12,13], the reduction or elimination of stepping jitter induced errors is highly desirable. This also applies to the directional darkfield [14,15] and anisotropic darkfield tomography techniques [16,17,18,19,20,21], which rely on the comparability of multiple darkfield images to relate signal variations to the anisotropy and orientation of scatterers. As the tomographic methods furthermore require large amounts (${10}^3$–${10}^4$) of projections, phase stepping series processing efficiency is also an issue. With regard to system calibration and monitoring, the inference of the actual stepping motions from the phase stepping series is of interest.

The present article proposes a simple iterative optimization algorithm both for the fitting of irregularly sampled sinusoids and in particular 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 a typical data set.

2. 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. The technique is further extended to an objective function considering spatially varying sampling positions to also model (erroneous) relative grating rotations and magnifications.

2.1. Sinusoid Fitting

The sinusoid is defined as a function $\tilde{y}(\varphi )$ of phase $\varphi$ parameterized by a constant offset o, an amplitude a and a phase offset ${\varphi }_0$:

$$\tilde{y}(\varphi )=o+a\sin (\varphi -{\varphi }_0).$$

(1)

For the purpose of fitting the sinusoid model to given data samples, the equivalent representation as a linear combination of trigonometric basis functions (i.e., as a first order Fourier series) will be more convenient:

$$\begin{array}{cc}\hfill o+a\sin (\varphi -{\varphi }_0)& =o+a\left(\cos {\varphi }_0\sin \varphi -\sin {\varphi }_0\cos \varphi \right)\hfill \\ & =o+a_s\sin \varphi +a_c\cos \varphi \hfill \\ & =(o,a_s,a_c)·{(1,\sin \varphi ,\cos \varphi )}^T,\hfill \end{array}$$

(2)

parameterized by the constant offset o and two auxiliary amplitudes $a_s$ and $a_c$ with the following identities:

$$\begin{array}{cc}\hfill a& =\sqrt{a_s^2+a_c^2}\hfill \\ \hfill \sin {\varphi }_0& =-a_c/a\hfill \\ \hfill \cos {\varphi }_0& =a_s/a\hfill \\ \hfill {\varphi }_0& =\arctan 2(-a_c,a_s).\hfill \end{array}$$

(3)

The parameters o, $a_c$, and $a_s$ of the recast sinusoid model corresponding to the least squares fit to given data samples $({\varphi }_i,y_i)$ enumerated by i may be determined by means of the following iterative scheme (introducing the superscript iteration index k):

$$\begin{array}{cc}\hfill o^{(0)},a_c^{(0)},a_s^{(0)}& =0,0,0\hfill \\ \hfill {\tilde{y}}_i^{(k)}& =o^{(k)}+a_c^{(k)}\cos {\varphi }_i+a_s^{(k)}\sin {\varphi }_i\hfill \\ \hfill o^{(k+1)}& =o^{(k)}+\frac{1}{N}\sum _i(y_i-{\tilde{y}}_i^{(k)})\hfill \\ \hfill a_s^{(k+1)}& =a_s^{(k)}+\frac{2}{N}\sum _i(y_i-{\tilde{y}}_i^{(k)})\sin {\varphi }_i\hfill \\ \hfill a_c^{(k+1)}& =a_c^{(k)}+\frac{2}{N}\sum _i(y_i-{\tilde{y}}_i^{(k)})\cos {\varphi }_i,\hfill \end{array}$$

(4)

with the factors of $1/N$ and $2/N$ accounting for the normalization of the respective basis functions and N being the amount of samples $({\varphi }_i,y_i)$ enumerated by i. The intermediate variables ${\tilde{y}}_i^{(k)}$ denote the respective ordinate values of the sinusoid model for iteration k and abcissas ${\varphi }_i$. The differences $y_i-{\tilde{y}}_i^{(k)}$ thus correspond to the residuals at iteration k. The scheme reduces to classic Fourier analysis for the case of the abscissas ${\varphi }_i$ being equidistantly distributed over multiples of $2\pi$ and converges within the first iteration in that case. As the incremental updates to $o^{(k)}$, $a_s^{(k)}$ and $a_c^{(k)}$ are proportional to the respective derivatives of the ${\ell }_2$ error ${\sum }_i{(o^{(k)}+a_c^{(k)}\cos {\varphi }_i+a_s^{(k)}\sin {\varphi }_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

$${\Delta }_{{\ell }_2}=\frac{\sqrt{{\sum }_i{\left(y_i-{\tilde{y}}_i^{(k-1)}\right)}^2}-\sqrt{{\sum }_i{\left(y_i-{\tilde{y}}_i^{(k)}\right)}^2}}{\sqrt{{\sum }_i{\left(y_i-{\tilde{y}}_i^{(k)}\right)}^2}}$$

(5)

may be tracked. It is typically found to fall below 0.1% within 10–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 ${\varphi }_i$ on multiples of $2\pi$, 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.

2.2. 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 ${\varphi }_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 of the involved actuators), a further optimization step will be introduced that minimizes the least squares error of the sinusoid fit over deviations $\Delta {\varphi }_i$ from the intended sampling positions ${\varphi }_i$. These deviations are expected to be smaller than the typical stepping increment ${\varphi }_{i+1}-{\varphi }_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 actual sampling positions ${\varphi }_i+\Delta {\varphi }_i$. In other words, an approach to the minimization problem

$$o_j,a_j,{\varphi }_{0,j},\Delta {\varphi }_i=\underset{o_j,a_j,{\varphi }_{0,j},\Delta {\varphi }_i}{\mathrm{argmin}}\sum _{i,j}{\left(o_j+a_j\sin ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_{0,j})-y_{ji}\right)}^2$$

(6)

shall be considered, with j indexing phase stepping curves captured by different detector pixels at identical stepping positions ${\varphi }_i+\Delta {\varphi }_i$.

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. Figure 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).

Algorithm 1 Least squares optimization of shared abscissa values ${\varphi }_i$ for simultaneous sinusoid fits to ordinate samples $y_{ji}$ belonging to independent curves j sampled at identical positions ${\varphi }_i$. This represents a special case of Algorithm 2 with spatially invariant sampling phases. The relaxation parameter ${\lambda }_k\in (0;1]$ may be chosen <1 if damping of the updates to ${\varphi }_i^{(k)}$ is desired. For the intermediate argmin operations, see Equations (2)–(5).

1: ${\varphi }_i,y_{ji}$: input data 2: $m_0\leftarrow \frac{1}{2}$                                     ▹ upper limit to $\Delta {\varphi }_{ji},m_0\in (0;1.38]$ 3: ${\varphi }_i^{(0)}\leftarrow {\varphi }_i$ 4: $o_j^{(0)},a_j^{(0)},{\varphi }_{0,j}^{(0)}\leftarrow \underset{o,a,{\varphi }_0}{\mathrm{argmin}}{\sum }_i{\left(o+a\sin ({\varphi }_i^{(0)}-{\varphi }_0)-y_{ji}\right)}^2$                       ▹ initialization 5: for$k=0..k_{max}$do 6:   $\Delta {\varphi }_i^{(k)}\leftarrow \frac{{\sum }_j\mathrm{softlimit}\left(\cos ({\varphi }_i^{(k)}-{\varphi }_{0,j}^{(k)})\left(a_j^{(k)}(y_{ji}-o_j^{(k)})-a_j^2\sin ({\varphi }_i^{(k)}-{\varphi }_{0,j}^{(k)})\right),{(a_j^{(k)})}^2{m}_0{\cos }^4({\varphi }_i^{(k)}-{\varphi }_{0,j}^{(k)})\right)}{{\sum }_j{(a_j^{(k)})}^2{\cos }^2({\varphi }_i^{(k)}-{\varphi }_{0,j}^{(k)})}$ 7:   ${\varphi }_i^{(k+1)}\leftarrow {\varphi }_i^{(k)}+{\lambda }_k\Delta {\varphi }_i^{(k)}$ 8:   $o_j^{(k+1)},a_j^{(k+1)},{\varphi }_{0,j}^{(k+1)}\leftarrow \underset{o_j,a_j,{\varphi }_{0,j}}{\mathrm{argmin}}{\sum }_i{\left(o_j+a_j\sin ({\varphi }_i^{(k+1)}-{\varphi }_{0,j})-y_{ji}\right)}^2$ 9: end for.

Algorithm 2 Simultaneous least squares optimization of abscissa values ${\varphi }_{ji}$ and sinusoid fits to ordinate samples $y_{ji}$ belonging to independent curves j sampled at positions ${\varphi }_{ji}={\varphi }_i(j)$ with ${\varphi }_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 ${\varphi }_i(j)$.

1: ${\varphi }_i,y_{ji}$: input data 2: $m_0\leftarrow \frac{1}{2}$                                       ▹ upper limit to $\Delta {\varphi }_{ji},m_0\in (0;1.38]$ 3: ${\varphi }_{ji}^{(0)}\leftarrow {\varphi }_i\forall j$                             ▹ initialization of sampling phases with intended values 4: $o_j^{(0)},a_j^{(0)},{\varphi }_{0,j}^{(0)}\leftarrow \underset{o,a,{\varphi }_0}{\mathrm{argmin}}{\sum }_i{\left(o+a\sin ({\varphi }_i^{(0)}-{\varphi }_0)-y_{ji}\right)}^2$                      ▹ initial sinusoid fits 5: for$k=0..k_{max}$do 6:   $\Delta {\varphi }_{ji}^{(k)}\leftarrow \left\{\begin{array}{cc}0\hfill & a_j^{(k)}\cos ({\varphi }_{ji}^{(k)}-{\varphi }_{0,j}^{(k)})=0\hfill \\ \mathrm{softlimit}\left(\frac{(y_{ji}-o_j^{(k)})/a_j^{(k)}-\sin ({\varphi }_{ji}^{(k)}-{\varphi }_{0,j}^{(k)})}{\cos ({\varphi }_{ji}^{(k)}-{\varphi }_{0,j}^{(k)})},m_0{\cos }^2({\varphi }_{ji}^{(k)}-{\varphi }_{0,j}^{(k)})\right)\hfill & \mathrm{else}\hfill \end{array}\right.$ 7:   $w_{ji}^{(k)}\leftarrow a_j^{(k)2}{\cos }^2({\varphi }_{ji}^{(k)}-{\varphi }_{0,j}^{(k)})$ 8:   $\Delta {\varphi }_i^{(k)},{\nabla }_h{\varphi }_i^{(k)},{\nabla }_v{\varphi }_i^{(k)},{\nabla }_{hv}{\varphi }_i^{(k)},{\nabla }_{h^2}{\varphi }_i^{(k)}\leftarrow \underset{\Delta {\varphi }_i,{\nabla }_h{\varphi }_i,{\nabla }_v{\varphi }_i,{\nabla }_{hv}{\varphi }_i,{\nabla }_{h^2}{\varphi }_i}{\mathrm{argmin}}{\sum }_j{w}_{ji}{\left(\Delta {\varphi }_i(j)-\Delta {\varphi }_{ji}^{(k)}\right)}^2$                                               ▹ for $\Delta {\varphi }_i(j)$, cf. Equation (19) 9:   ${\varphi }_{ji}^{(k+1)}\leftarrow {\varphi }_{ji}^{(k)}+\Delta {\varphi }_i^{(k)}+{\nabla }_h{\varphi }_i^{(k)}(h-h_0)+{\nabla }_v{\varphi }_i^{(k)}(v-v_0)+{\nabla }_{hv}{\varphi }_i^{(k)}(h-h_0)(v-v_0)+{\nabla }_{h^2}{\varphi }_i^{(k)}{(h-v_0)}^2$ 10:   $o_j^{(k+1)},a_j^{(k+1)},{\varphi }_{0,j}^{(k+1)}\leftarrow \underset{o_j,a_j,{\varphi }_{0,j}}{\mathrm{argmin}}{\sum }_i{\left(o_j+a_j\sin ({\varphi }_{ji}^{(k+1)}-{\varphi }_{0,j})-y_{ji}\right)}^2$ 11: end for.

2.2.1. Determination of Individual Phase Deviations

Starting with an initial sinusoid fit

$$o,a,{\varphi }_0=\underset{o,a,{\varphi }_0}{\mathrm{argmin}}\sum _i{\left(o+a\sin ({\varphi }_i-{\varphi }_0)-y_i\right)}^2,$$

(7)

to given data samples $({\varphi }_i,y_i)$, the residual sum of squares shall be minimized over deviations $\Delta {\varphi }_i$ to the abscissas ${\varphi }_i$ while keeping the sinusoid model parameters o, a and ${\varphi }_0$ fixed:

$$\Delta {\varphi }_i=\underset{\Delta {\varphi }_i}{\mathrm{argmin}}\sum _i{\left(o+a\sin ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_0)-y_i\right)}^2.$$

(8)

For reasons of better readability, the detector pixel index (or equivalently PSC index) j of the sinusoid parameters o, a, ${\varphi }_0$ and the data samples $y_i$ has been omitted here and is explicitly added again to the final results at the end of this subsection.

Equation (8) is solved when the derivative of the objective function with respect to $\Delta {\varphi }_i$ vanishes:

$$\begin{array}{c}0=\frac{\mathrm{d}}{\mathrm{d}\Delta {\varphi }_i}\sum _i{\left(o+a\sin ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_0)-y_i\right)}^2\hfill \\ 0=2\left(o-y_i+a\sin ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_0)\right)a\cos ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_0)\hfill \\ \mathrm{for}a\cos ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_0)\ne 0:\hfill \\ 0=\left(o-y_i+a\sin ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_0)\right)\hfill \\ \mathrm{for}\Delta {\varphi }_i\ll \pi :\hfill \\ 0\approx \left(o-y_i+a\sin ({\varphi }_i-{\varphi }_0)+\Delta {\varphi }_{i}a\cos ({\varphi }_i-{\varphi }_0)\right)\hfill \end{array}$$

(9)

where the last step is a first order Taylor expansion with respect to $\Delta {\varphi }_i$. This directly leads to the following expression for $\Delta {\varphi }_i$:

$$\Delta {\varphi }_i\approx \frac{1}{\cos ({\varphi }_i-{\varphi }_0)}\left(\frac{y_i-o}{a}-\sin ({\varphi }_i-{\varphi }_0)\right)\mathrm{for}\Delta {\varphi }_i\ll \pi \mathrm{and}a\cos ({\varphi }_i-{\varphi }_0)\ne 0,$$

(10)

where the earlier condition $a\cos ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_0)\ne 0$ is approximated to be satisfied when $a\cos ({\varphi }_i-{\varphi }_0)\ne 0$. The restriction to cases with $a\cos ({\varphi }_i-{\varphi }_0)\ne 0$ can be intuitively understood when recalling that $\cos ({\varphi }_i-{\varphi }_0)=0$ implies a maximum or minimum of the sinusoid and $a=0$ means that it is constant (${\varphi }_i$ independent), in both of which cases there is no sensible choice for $\Delta {\varphi }_i\ne 0$. The constraint on the result, $\Delta {\varphi }_i\ll \pi$, can simply be taken into account by means of a limiting function parameterized by a maximal absolute value $m\ge 0$ such as

$$\begin{array}{cc}\hfill \mathrm{softlimit}(\Delta {\varphi }_i,m\ge 0)& =\left\{\begin{array}{cc}0\hfill & m=0\hfill \\ m\tanh \left(\frac{\Delta {\varphi }_i}{m}\right)\hfill & \mathrm{else}\hfill \end{array}\right.,\hfill \end{array}$$

(11)

which converges to the identity function for $\Delta {\varphi }_i\ll m$ and is bounded at $\pm m$. The choice of m in this case depends on the validity range of the linear approximation of $\sin ({\varphi }_i+\Delta {\varphi }_i-{\varphi }_0)$ with respect to $\Delta {\varphi }_i$ about ${\varphi }_i-{\varphi }_0$, which obviously depends on the magnitude of the curvature of the sinusoid at this point, as illustrated in Figure 4. This may be accounted for by introducing a suitable ${\varphi }_i-{\varphi }_0$ dependence to m denoted by the respective argument:

$$m({\varphi }_i-{\varphi }_0)=m_0{\cos }^2({\varphi }_i-{\varphi }_0).$$

(12)

$m({\varphi }_i-{\varphi }_0)$ reaches its largest value $m_0$ at ${\varphi }_i-{\varphi }_0=0$ (where $\sin ({\varphi }_i-{\varphi }_0)$ is actually linear) and smoothly drops to 0 for $\cos ({\varphi }_i-{\varphi }_0)=0$, in which case both the sine and its curvature are extremal and $\Delta {\varphi }_i$ shall and will be limited to 0. The upper bound for $m({\varphi }_i-{\varphi }_0)$ and thus for the choice of $m_0$ is defined by the range of values $\Delta {\varphi }_i\in [-m_0{\cos }^2({\varphi }_i-{\varphi }_0),+m_0{\cos }^2({\varphi }_i-{\varphi }_0)]$ over which $\sin (\Delta {\varphi }_i+{\varphi }_i-{\varphi }_0)$ is actually invertible (cf. Figure 4). The sine function $\sin (\vartheta )$ is locally invertible over intervals of $\vartheta \in [n\pi -\frac{\pi }{2},n\pi +\frac{\pi }{2}]$ with n being an integer number and $\vartheta$ denoting the argument ($\vartheta =\Delta {\varphi }_i+{\varphi }_i-{\varphi }_0$). For the determination of $m_0$ it is sufficient to consider the case $n=0$:

$$\begin{array}{cc}\hfill m_0{\cos }^2(\vartheta )& \le \frac{\pi }{2}-\left|\vartheta \right|\forall \vartheta \in [-\frac{\pi }{2},+\frac{\pi }{2}]\hfill \\ \hfill m_0& \lesssim 1.38.\hfill \end{array}$$

(13)

For $m_0=1.38$, the linear approximation used in Equation (9) deviates by up to $40\%$. The deviation is limited to $20\%$ or $5\%$ for $m_0=1$ and $m_0=\frac{1}{2}$, respectively.

Combining the above results and reintroducing the detector pixel index j, the following expression for pixel (j) and phase step (i) wise phase stepping deviations $\Delta {\varphi }_{ji}$ reducing (and, for sufficiently small $\Delta {\varphi }_{ji}$, minimizing) the ${\ell }_2$ error ${\left(o_j+a_j\sin ({\varphi }_i+\Delta {\varphi }_{ji}-{\varphi }_{0,j})-y_{ji}\right)}^2$ can be given, choosing $m_0=\frac{1}{2}$:

$$\Delta {\varphi }_{ji}\approx \left\{\begin{array}{cc}0\hfill & a_j\cos ({\varphi }_i-{\varphi }_{0,j})=0\hfill \\ \mathrm{softlimit}\left(\frac{\left(\frac{y_{ji}-o_j}{a_j}-\sin ({\varphi }_i-{\varphi }_{0,j})\right)}{\cos ({\varphi }_i-{\varphi }_{0,j})},\frac{1}{2}{\cos }^2({\varphi }_i-{\varphi }_{0,j})\right)\hfill & \mathrm{else}\hfill \end{array}\right..$$

(14)

Figure 3 shows an example of this approximate least squares solution to $\Delta {\varphi }_{ji}$.

2.2.2. Noise Weighted Average of Phase Deviations

Now that an expression has been derived for the deviations $\Delta {\varphi }_{ji}$ optimizing the abscissa values for individual PSCs indexed by j given previous sinusoid fits, the respective results for all PSCs sharing the same sampling locations may be averaged, such that the mean deviations $\Delta {\varphi }_i$ essentially correspond to the actual systematic phase stepping errors (whereas statistical noise will mostly cancel out):

$$\Delta {\varphi }_i=\frac{{\sum }_j{w}_{ji}\Delta {\varphi }_{ji}}{{\sum }_j{w}_{ji}},$$

(15)

using weights $w_{ji}$ factoring in the relative certainty and relevance of the individual $\Delta {\varphi }_{ji}$. An appropriate choiceis

$$w_{ji}=a_j^2{\cos }^2({\varphi }_i-{\varphi }_{0,j}),$$

(16)

where ${\cos }^2({\varphi }_i-{\varphi }_{0,j})$ weights the slope dependent error propagation from noisy measurements $y_{ji}$ to $\Delta {\varphi }_{ji}$ based on the derivative of the sinusoid model at phase step i and $a_j^2$ weights the contribution of a particular PSC j to the accumulated ${\ell }_2$ error. These considerations leadto

$$\begin{array}{cc}\hfill \Delta {\varphi }_i& =\frac{{\sum }_j{a}_j^2{\cos }^2({\varphi }_i-{\varphi }_{0,j})\Delta {\varphi }_{ji}}{{\sum }_j{a}_j^2{\cos }^2({\varphi }_i-{\varphi }_{0,j})}\hfill \\ & =\frac{{\sum }_j{a}_j^2{\cos }^2({\varphi }_i-{\varphi }_{0,j})\mathrm{softlimit}\left(\frac{\left(\frac{y_{ji}-o_j}{a_j}-\sin ({\varphi }_i-{\varphi }_{0,j})\right)}{\cos ({\varphi }_i-{\varphi }_{0,j})},\frac{1}{2}{\cos }^2({\varphi }_i-{\varphi }_{0,j})\right)}{{\sum }_j{a}_j^2{\cos }^2({\varphi }_i-{\varphi }_{0,j})}\hfill \\ & =\frac{{\sum }_j\mathrm{softlimit}\left(\cos ({\varphi }_i-{\varphi }_{0,j})\left(a_j(y_{ji}-o_j)-a_j^2\sin ({\varphi }_i-{\varphi }_{0,j})\right),\frac{1}{2}{a}_j^2{\cos }^4({\varphi }_i-{\varphi }_{0,j})\right)}{{\sum }_j{a}_j^2{\cos }^2({\varphi }_i-{\varphi }_{0,j})},\hfill \end{array}$$

(17)

where the last step uses the relation $\alpha \mathrm{softlimit}(x,m)=\mathrm{softlimit}(\alpha x,\alpha m)$ for arbitrary real-valued arguments x and m and factors $\alpha \ge 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 ${\varphi }_i$ as defined by Equation (6). A pseudo code representation is given in Algorithm 1, further introducing the relaxation parameter ${\lambda }_k\in (0;1]$ that may be chosen <1 in order to damp the adaptions to ${\varphi }_i^{(k)}$ if desired. The intermediate sinusoid fits may be accomplished using the iterative algorithm described in the previous section.

2.2.3. 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 $\Delta {\varphi }_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 gradients is, analog to Equation (6), given by:

$$o_j,a_j,{\varphi }_{0,j},\Delta {\varphi }_i(j)=\underset{o_j,a_j,{\varphi }_{0,j},\Delta {\varphi }_i(j)}{\mathrm{argmin}}\sum _{i,j}{\left(o_j+a_j\sin ({\varphi }_i+\Delta {\varphi }_i(j)-{\varphi }_{0_j})-y_{ji}\right)}^2$$

(18)

where the spatial dependence of the phase deviations $\Delta {\varphi }_i$ has been accounted for by a functional dependence on the detector pixel index j. Given the expected gradients (as illustrated in Figure 5), the spatially varying phase deviations $\Delta {\varphi }_i(j)$ have the following form:

$$\Delta {\varphi }_i(j)=\Delta {\varphi }_i+{\nabla }_h{\varphi }_i(h-h_0)+{\nabla }_v{\varphi }_i(v-v_0)+{\nabla }_{hv}{\varphi }_i(h-h_0)(v-v_0)+{\nabla }_{h^2}{\varphi }_i{(h-h_0)}^2,$$

(19)

with the coefficients ${\nabla }_h{\varphi }_i$, ${\nabla }_v{\varphi }_i$, ${\nabla }_{hv}{\varphi }_i$ and ${\nabla }_{h^2}{\varphi }_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:

$$\begin{array}{cc}\hfill j& =v{N}_h+h\hfill \\ \hfill h& =j\mathrm{mod}{N}_h\hfill \\ \hfill v& =(j-h)/N_h.\hfill \end{array}$$

(20)

The constant offsets $h_0$ and $v_0$ characterize the detector center.

The optimization of the extended objective function in Equation (18) can be performed analog to that of Equation (6) with the only difference lying in the evaluation of the spatial phase difference maps $\Delta {\varphi }_{ji}$ defined by Equation (14) (an example of which is shown in Figure 6). The weighted average derived in the previous section to determine the homogeneous offset $\Delta {\varphi }_i$ can be extended to a generalized linear least squares fit of the model $\Delta {\varphi }_i(j)=\Delta {\varphi }_i(h(j),v(j))$ defined by Equation (19) to the local estimates $\Delta {\varphi }_{ji}$ (Equation (14), Figure 6), also taking the weights defined by Equation (16) 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 ${\nabla }_h{\varphi }_i$, ${\nabla }_v{\varphi }_i$, ${\nabla }_{hv}{\varphi }_i$, ${\nabla }_{h^2}{\varphi }_i$ and relative translatory and rotatory motions of the interferometer’s gratings. 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.

$\Delta {\varphi }_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 ${\nabla }_{h^2}{\varphi }_i{(h-h_0)}^2$ term in the case of non-vanishing curvature of the spatial phase deviation.

The vertical gradient parameter ${\nabla }_v{\varphi }_i$ is related to a relative rotation $\eta$ of both gratings about the optical axis:

$$\begin{array}{cc}\hfill \tan \eta & =\frac{{\nabla }_v{\varphi }_i}{2\pi }\frac{\mathrm{effective}\mathrm{grating}\mathrm{period}}{\mathrm{detector}\mathrm{pixel}\mathrm{pitch}},\hfill \end{array}$$

(21)

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 ${\nabla }_h{\varphi }_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):

$$\mathrm{relative}\mathrm{period}\mathrm{mismatch}=\frac{\mathrm{effective}\mathrm{period}\mathrm{difference}}{\mathrm{effective}\mathrm{grating}\mathrm{period}}=\frac{{\nabla }_h{\varphi }_i}{2\pi }\frac{\mathrm{effective}\mathrm{grating}\mathrm{period}}{\mathrm{detector}\mathrm{pixel}\mathrm{pitch}}.$$

(22)

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:

$$\begin{array}{cc}\hfill \mathrm{translation}\mathrm{distance}& =(\mathrm{relative}\mathrm{period}\mathrm{mismatch})\frac{{(\mathrm{source}-\mathrm{grating}\mathrm{distance})}^2}{\mathrm{source}-\mathrm{detector}\mathrm{distance}}\hfill \\ & =\frac{{\nabla }_h{\varphi }_i}{2\pi }\frac{\mathrm{effective}\mathrm{grating}\mathrm{period}}{\mathrm{detector}\mathrm{pixel}\mathrm{pitch}}\frac{{(\mathrm{source}-\mathrm{grating}\mathrm{distance})}^2}{\mathrm{source}-\mathrm{detector}\mathrm{distance}}.\hfill \end{array}$$

(23)

The change ${\nabla }_{hv}{\varphi }_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 $\theta$ of one of the gratings about the horizontal axis. Using the above relation between magnification changes and spatial displacements, the tilt $\theta$ about the horizontal axis is related to ${\nabla }_{hv}{\varphi }_i$ approximately by

$$\begin{array}{cc}\hfill \tan \theta & =\frac{{\nabla }_{hv}{\varphi }_i}{2\pi }\frac{\mathrm{effective}\mathrm{grating}\mathrm{period}}{\mathrm{detector}\mathrm{pixel}\mathrm{pitch}}{\left((\mathrm{detector}\mathrm{pixel}\mathrm{pitch})\frac{\mathrm{source}--\mathrm{grating}\mathrm{distance}}{\mathrm{source}--\mathrm{detector}\mathrm{distance}}\right)}^{-1}\hfill \\ & \hfill \times \frac{{(\mathrm{source}--\mathrm{grating}\mathrm{distance})}^2}{\mathrm{source}-\mathrm{detector}\mathrm{distance}}\\ & =\frac{{\nabla }_{hv}{\varphi }_i}{2\pi }\frac{(\mathrm{effective}\mathrm{grating}\mathrm{period})(\mathrm{source}-\mathrm{grating}\mathrm{distance})}{{(\mathrm{detector}\mathrm{pixel}\mathrm{pitch})}^2}.\hfill \end{array}$$

(24)

A non-vanishing curvature ${\nabla }_{h^2}{\varphi }_i$ arises in the case of a rotary motion about the vertical axis (slant) and is analogously related to the slant angle $\phi$ to first order by

$$\tan \phi =\frac{{\nabla }_{h^2}{\varphi }_i}{2\pi }\frac{(\mathrm{effective}\mathrm{grating}\mathrm{period})(\mathrm{source}--\mathrm{grating}\mathrm{distance})}{{(\mathrm{detector}\mathrm{pixel}\mathrm{pitch})}^2}.$$

(25)

2.3. Experimental Setup

The interferometer used to demonstrate the described analysis techniques consists of a set of three gratings (coherence grating G0, phase grating G1 and absorption grating G2) manufactured by microworks GmbH (Karlsruhe, Germany) for a design energy of 54keV and mounted in an in house cone beam micro-CT setup comprising a commercial microfocus X-ray source operated at 80 kV acceleration voltage as well as a commercial flat panel detector with a pixel pitch of 74.8 µm. The grating periods are 4.8 µm (G0), 2.4 µm (G1) and 4.8 µm (G2), and the gratings are placed at 100 cm, 125 cm and 150 cm distance from the X-ray source, respectively. The field of view is limited by the diameter of G2 of approximately 10 cm. The G1 grating is mounted on a piezo driven linear actuator responsible for the phase stepping. The detector is placed right behind the absorption grating in about 155 cm distance from the source. Images are averaged over 10 exposures of 0.5 s each and are cropped to a region of 850 × 850 px${}^2$ corresponding to an area of 6.4 × 6.4 cm${}^2$ on the detector. A piece of plastic hose of 2 cm diameter placed about 143 cm distance from the source serves as the sample object.

3. 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 ${\varphi }_{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 Equation (17) for all 15 frames of the phase stepping series. Figure 6 provides an example for the first frame of the series. By iterating the sinusoid fits and the corrections to the sampling positions to reduce the overall least squares error (cf. Equation (6)) 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 Figure 9. Correspondingly, the deviations from the intended phase stepping positions are almost completely deduced within the first iteration, as shown in Figure 7 (bottom). Nevertheless, complete suppression of the Moiré artifacts within the final images requires some further iterations, as illustrated in Figure 10.

In addition to the mean deviations of the phase steps from the intended positions, spatial gradients throughout the detection area have also been considered (cf. Algorithm 2). Figure 11 shows the respective differential deviations from the intended phase steps between the first nine frames of the phase stepping series, normalized to the nominal homogeneous phase stepping increment of $2\pi /5$. The mean contributions of each component are listed in

Table 1. Root mean square contributions of the mean, gradient and curvature components of $\Delta {\varphi }_i(j)$ to the sampling phase deviations found for the present phase stepping series in units of radians. The homogeneous error $\Delta {\varphi }_i$ is by far the dominating effect. The contributions of ${\nabla }_{hv}{\varphi }_i$ and ${\nabla }_{h^2}{\varphi }_i$ range in the order of magnitude of the expected noise level of ${10}^{-4}\mathrm{rad}$ (cf. Equations (26) and (27)).

$\sqrt{\overline{\Delta {\mathit{\varphi }}_{\mathit{i}}^2}}$	$\sqrt{\overline{{\left({\mathit{\nabla }}_{\mathit{h}}{\mathit{\varphi }}_{\mathit{i}}(\mathit{h}-{\mathit{h}}_0)\right)}^2}}$	$\sqrt{\overline{{\left({\mathit{\nabla }}_{\mathit{v}}{\mathit{\varphi }}_{\mathit{i}}(\mathit{v}-{\mathit{v}}_0)\right)}^2}}$	$\sqrt{\overline{{\left({\mathit{\nabla }}_{\mathbf{hv}}{\mathit{\varphi }}_{\mathit{i}}(\mathit{h}-{\mathit{h}}_0)(\mathit{v}-{\mathit{v}}_0)\right)}^2}}$	$\sqrt{\overline{{\left({\mathit{\nabla }}_{{\mathit{h}}^2}{\mathit{\varphi }}_{\mathit{i}}{(\mathit{h}-{\mathit{h}}_0)}^2\right)}^2}}$
$1.2\times {10}^{-1}$	$2.4\times {10}^{-3}$	$3.6\times {10}^{-3}$	$3.6\times {10}^{-4}$	$7.5\times {10}^{-4}$

. While the homogeneous error of $0.1\mathrm{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. Figure 8).

Figure 12 shows variations in the relative alignment of the gratings derived from the inhomogeneous phase stepping analysis by means of Equations (21)–(25). 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 12 indicate rotations about an off-center pivot point about ${10}^{-1}\mathrm{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 (Figure 12, 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 well explicable by temperature variations in the order of magnitude of ${10}^{-1}\mathrm{K}$ given a thermal expansion coefficient in the order of magnitude of ${10}^{-6}{\mathrm{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 Figure 12).

As a crude error assessment, the standard error of the mean phase deviation can be estimated from the sinusoid fits’ root mean square error:

$$\begin{array}{cc}\hfill 31{\sigma }_{\mathrm{mean}}& \approx \frac{1}{\sqrt{\mathrm{contributing}\mathrm{detector}\mathrm{pixels}}}\frac{\sin \mathrm{usoid}\mathrm{fit}\mathrm{RMSE}}{\mathrm{mean}\sin \mathrm{usoid}\mathrm{amplitude}}\hfill \\ & \approx \sqrt{\frac{2}{\mathrm{detector}\mathrm{pixels}}}\frac{\sin \mathrm{usoid}\mathrm{fit}\mathrm{RMSE}}{\mathrm{mean}\sin \mathrm{usoid}\mathrm{amplitude}}.\hfill \end{array}$$

(26)

The latter corresponds for the present data set to about $6\%$ of the mean observed sinusoid amplitude, which directly translates to $6\times {10}^{-2}\mathrm{rad}$ on the abscissa. Given the amount of detector pixels contributing to the least squares fits of $\Delta {\varphi }_i(j)$ within each frame of the phase stepping series, a standard error in the order of magnitude of

$${\sigma }_{\mathrm{mean}}\approx {10}^{-4}\mathrm{rad}$$

(27)

results. This implies that, according to the results given in Table 1, the tilt and slant contributions ${\nabla }_{hv}\varphi$ and ${\nabla }_{h^2}\varphi$ are close to the expected noise level for the present case.

4. Discussion and Conclusions

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. The additional effort (with respect to [5,6]) of explicitly optimizing phases rather than Fourier coefficients allows for a straight forward extension of the optimization algorithm also in the case of spatially varying phase stepping increments due to further mechanical degrees of freedom besides the intended translatory stepping motion. By division of the full optimization task into three easily tractable subproblems (phase stepping curve fitting, identification of sampling position deviations and fitting of the latter to an expected spatial model), the use of generic nonlinear optimization algorithms as, e.g., used in [7] is avoided. Of these subtasks, only the second one is actually nonlinear, and a significant portion of the present article has been devoted to its approximate linearization (also taking noise propagation into consideration). The problem of simultaneous phase stepping curve evaluation and (spatially varying) sampling position determination is thus reduced to the iterative alternation of two generalized linear least squares problems and one linear approximation to a nonlinear problem. Due to this almost linear nature, sufficient convergence is achieved within less than five iterations for the presented example. Although no rigorous convergence analysis has been performed, the outlined structure of the optimization problem does not raise severe concerns regarding its stability. In cases of doubt, convergence can simply be slowed down by means of the relaxation parameter ${\lambda }_k$.

For the presented example data set, mean phase stepping errors of up to 10% of the nominal step width have been found, and their correction results in both a considerable reduction of the overall root mean square error by almost a factor of two and 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\mathrm{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. This crosstalk between sub-microrad rotations and effective translations further 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\mathrm{s}$ per phase stepping series (using graphics processors), it is well suitable as a standard processing method also for large image series.