## Abstract

A fast, powerful and stable filter based on combined wavelet and Fourier analysis for the elimination of horizontal or vertical stripes in images is presented and compared with other types of destriping filters. Strict separation between artifacts and original features allowing both, suppression of the unwanted structures and high degree of preservation of the original image information is endeavoured. The results are validated by visual assessments, as well as by quantitative estimation of the image energy loss. The capabilities and the performance of the filter are tested on a number of case studies related to applications in tomographic imaging. The case studies include (i) suppression of waterfall artifacts in electron microscopy images based on focussed ion beam nanotomography, (ii) removal of different types of ring artifacts in synchrotron based X-ray microtomography and (iii) suppression of horizontal stripe artifacts from phase projections in grating interferometry.

©2009 Optical Society of America

## 1. Introduction

The reduction of horizontal or vertical stripe artifacts is an important topic in many imaging systems using multiple-detectors, such as in satellite applications [1, 2, 3], spectroradiometry [4, 5], as well as others [6]. Suppression of stripe artifacts is important also in the case of tomographic microscopies, as such artifacts decrease the quality of obtained datasets. Ring artifacts in reconstructed tomographic slices are strictly related to stripe artifacts, since they result from the back-projection of vertical lines or bands observed in the sinogram regime.

Successful removal of stripe artifacts implies, that (a) any horizonal or vertical stripe disappears from an image after filtering, while (b) all structural features, which are different from stripes, and (c) the quantitative values of the image information are optimally preserved. While the human perception is appropriate for visually detecting even minute stripe residues and identifying missing structures, it fails when the quantitative information of an image needs to be assessed. Yet for any automated image evaluation procedure, the quantitative aspect is essential.

Several different methods for stripe artifact reduction have been proposed ranging from moving average filtering [7, 8], histogram matching [9], interpolation approaches [10], frequency filtering with fast Fourier transform (FFT) [3, 11], as well as approaches using wavelets [1, 5, 12]. A large number of algorithms to reduce ring artifacts in tomographic reconstructions, was also suggested [7, 8, 12, 13, 14, 15, 16, 17], in addition to the standard flat field correction [18].

However, none of these techniques seems to completely satisfy the above mentioned requirements (a) - (c). The authors acknowledge that they might show top-quality results in many different cases. However, the general problem which is common to probably any ever existing filter is the fact that, while features that are supposed to be removed are never removed completely, some information to be retained might be seriously damaged, instead. Or in other words, the human capacity to visually discriminate undesired textures usually outperforms the filter selectivity by far. The consequential approach to a solution with reference to the striping filter is providing a methodology to (i) increase tight condensation of striping information into strictly isolated values, and to (ii) provide further grouping of those values into well categorized subclasses permitting even further and more subtle selection of the processed structure.

We developed and present hereinafter a new destriping filter based on wavelet decomposition and Fourier filtering, which tackles both points (i) and (ii) simultaneously, while the other mentioned approaches lack in this aspect. The chosen decomposition approach allows for a more strict separation between artifacts and original features, making the suppression of just the unwanted structure feasible. In this way, stripe and ring artifacts can be removed while preserving the original image information.

Different approaches for such condensing / categorizing of striping information might be possibly based also on other types of transforms. Especially, ridgelet transforms [19, 20, 21] (which are well known for efficiently pinpointing edge features and, in particular, lines of the size of the image), could be envisaged to be used for the condensation of striping artifacts under arbitrary angle. The possible design, potential and drawback of a such approach would have to be evaluated in detail in future studies.

After a short review of the required assumptions in wavelet and Fourier transforms, the new filter is introduced and compared with two approaches from the literature, which are both part of our new concept. The performance of the filter is assessed according to the requirements (a)-(c). In addition, the role of some of the free parameters is examined and the quantitative changes of the image values after filtering are investigated with the objective to keep them as small as possible. Finally, some case studies are presented in the application part (Section 4), with a particular focus on x-ray computed tomography.

#### 1.1. Some considerations about wavelet and Fourier transforms

In both numerical wavelet and Fourier analysis, a discrete signal *f*(*t*) can be approximated by a set of basis functions Γ* _{n}*(

*t*),

*n*∈ {0, ⋯,

*N*}, yielding

where the basis functions are orthogonal to each other, that is, the scalar product of every two Γ_{n}(*t*) meets

Thereby, *f*(*t*) will be unequivocally represented by the coefficients *a _{n}*.

The main motivation for the decomposition of a signal into orthogonal basis functions is the deployment of the original signal information into coefficient classes that specifically group interesting structural patterns. This aspect is particularly attractive for artifacts removal techniques. In addition, orthogonal transforms often yield coefficients, which become partly small or even zero. This characteristic of wavelet and Fourier analysis makes signal decomposition especially attractive for data compression purposes [22]. Moreover, the calculation of the coefficients *a _{n}* is very efficient and achieved in linear time

*o*(

*N*) by simple evaluation of scalar products:

Furthermore, since Fourier and some wavelet transforms are separable, they can be applied successively for each dimension in the case of two or more dimensional signals.

For the discrete Fourier transform (DFT), Γ_{n}(*t*) are sine and cosine functions at discrete frequencies *w _{n}* = 2·

*π*·

*n*, although the complex exponential notation Γ

_{n}(

*t*) =

*e*

^{-i·ωn·t}is often preferred. In the Fourier domain, all structural information of a signal is entirely represented by frequencies, whereas any spatial information is completely distributed over the entire frequency range and thus hidden. In contrast, for the wavelet transform [23, 24, 25], Γ

_{n}(

*t*) are narrow groups of wavelets of different sizes and frequencies and both, spatial and frequency information is retained. For the discrete wavelet transform (DWT), a variety of wavelet types exist, which allow perfect reconstruction of the original signal, in particular the Daubechies (DB) wavelets [26] (Fig.1).

In dyadic, decimated wavelet transform, a single scale wavelet decomposition of a 1D signal *f*(*t*) consists in its fragmentation into a low and a high frequency part. This separation is perfectly reversible. The low frequency part can then be iteratively decomposed in the same way, while leaving the high frequency unchanged. After a decomposition step, the size *N*(*l*) of the resulting low and high frequency parts is halved in order to maintain the total number of coefficients constant. The next decomposition is therefore performed at half the number of coefficients but with the same filter size, i.e. at lower spatial scale, what is equivalent to filtering at double the filter size. A multiscale wavelet transform at a highest decomposition level *L* is therefore the subdivision of *f*(*t*) into one low frequency part, which is represented by one scaling function Φ_{L}(*t*) together with its coefficients, and multiple high frequency parts, which are represented by a set of wavelet functions Ψ_{l}(*t*) and their coefficients at *l* = 1,⋯,*L* different scales, yielding

For a given set consisting of one scaling function Φ_{L,n} and a few wavelet functions Ψ_{l,n}(*t*) and its translated versions, *f*(*t*) will be represented by the low frequency coefficients *c*
_{L,n} at the scale of level *L*, and by the high frequency coefficients *d*
_{l,n} at the scales of *l* = 1,⋯,*L*.

Due to the reduction of the number of coefficients by a factor of two after each decomposition step, wavelet functions Ψ_{l,n} are attributed with dyadic values. In this way, all wavelets Ψ_{l,n} can be derived from a so called mother wavelet Ψ_{0} with

In the case of a discrete signal *f*(*t _{i}*), dyadic, decimated wavelet analysis can be implemented with the help of quadrature mirror filters, as proposed by Mallat [23], by successively splitting the input signal into a low frequency and a details band by applying a low pass and a high pass filter, respectively. The resulting and downsampled low pass band, which represents the approximation coefficients at a coarser scale, is successively decomposed in the same way by filtering at the next scale. The set of low and high pass filters at different scales is called a filter bank. The data flow of a 1D multiscale wavelet decomposition is sketched in Fig.2.

Looking at a single scale wavelet decomposition step of a 2D signal by using filter banks, the signal *f*(*x*,*y*) is decomposed into a set of four coefficient bands *c _{l}*,

*c*,

_{h}*c*,

_{v}*c*, where

_{d}*c*is the low pass band, and

_{l}*c*,

_{h}*c*,

_{v}*c*the horizontal, vertical and diagonal details bands, respectively. The number of coefficients in each band is 1/4 of the number in the original band. Thereby, 2D multiscale wavelet decomposition of level

_{d}*L*yields

$$+\sum _{l=1}^{L}\sum _{m}\sum _{n}{c}_{v\phantom{\rule[-0ex]{.2em}{0ex}}l,m,n}\xb7{\mathrm{\Psi}}_{v\phantom{\rule[-0ex]{.2em}{0ex}}l,m,n}\left(x,y\right)+\sum _{l=1}^{L}\sum _{m}\sum _{n}{c}_{d\phantom{\rule[-0ex]{.2em}{0ex}}l,m,n}\xb7{\mathrm{\Psi}}_{d\phantom{\rule[-0ex]{.2em}{0ex}}l,m,n}\left(x,y\right)$$

Consequently, the wavelet representation *W* of a 2D signal *f*(*x*,*y*) results in a set of coefficients

which allows recovery of the original information with zero-loss.

In Eq.6, the Φ_{L,m,n}(*x*,*y*) scaling function and the wavelets Ψ_{h l,m,n}(*x*,*y*), Ψ_{v l,m,n}(*x*,*y*), Ψ_{d l,m,n}(*x*,*y*) are 2D functions. However, since the wavelet transform is separable, the wavelet coefficients in Eq.7 can be calculated by successive filtering of the signal in horizontal and vertical directions applying either low (L) or high pass (H) filters (i.e. LL for *c*
_{l L,m,n}, LH for *c*
_{h l,m,n}, HL for *c*
_{v l,m,n}, HH for *c*
_{d l,m,n}).

An example is shown in Fig.3: a wavelet decomposition at level *L* = 4 is applied to the Lena image (Fig.3, left), a well known object in wavelet science [23, 24, 25]. The filter bank corresponding to the Haar wavelet (DB1) (Fig.1, left) was used.

## 2. Concept and implementation of the stripe filter

#### 2.1. Definition of a stripe artifact

An ideal “vertical stripe artifact” on an image *f*(*x*,*y*) is defined as a constant offset *A* of arbitrary width (*x _{e}*−

*x*),

_{a}running over the entire vertical dimension *y*.

Bar-like patterns running only over a limited part of *y* ∈ [*y _{a}*,⋯,

*y*] of the image are not considered as stripe artifacts, but as parts of the “correct” image information and should not be removed. The subsequent considerations are assuming vertical stripes, the case of horizontal ones can be handled analogously.

_{e}Although the proposed filter is specifically designed for the elimination of ideal stripes, its performance will also be discussed for imperfect stripes. It will in particular be shown how the filter parameters need to be tuned to extend the working range of the destripiong procedure to more general stripe artifacts.

#### 2.2. Basic idea

We assume that the original image *f*(*x*,*y*) is impaired by vertical stripe offsets as illustrated in the Lena image to the left of Fig.3.

### a) Wavelet filtering

Eq.6 shows that in 2D multiresolution wavelet decomposition, the vertical details components *c _{v}* are successively detached from all remaining image components. Consequently, the information from vertical stripes is exclusively condensed to

*c*

_{v l,m,n}and to the coefficients of the finally remaining low frequency band

*c*

_{l L,m,n}. Due to the dyadic fractionation of the spatial extensions of the wavelet basis functions (see Eq.5), each successive vertical details band

*c*

_{v l,m,n}

*l*∈ {1,⋯

*L*} is basically comprising a frequency band of dyadically decreasing focal frequency

*F*~ 2

_{f}^{-l}. Accordingly, the amount of detached stripe information enclosed in

*c*

_{v l,m,n}at each decomposition level

*l*depends on the spatial frequency spectrum of the stripes in horizontal direction, which correlates with the stripe width. Hence, the highest decomposition level

*L*required is coupled with the maximum expected stripe width. For a sufficiently large

*L*, the impact of the stripe information to the low pass coefficients

*c*

_{l L,m,n}becomes negligible.

In the example in Fig.3, centre, all information related to the vertical stripes impairing the original Lena image is clearly detached exclusively to the vertical details bands (marked with red frames) of the decomposed image. Blanking of *c*
_{v l,m,n} at some selected *l* and subsequent inverse wavelet transform of the remaining wavelet coefficients will yield a modified version of *f*(*x*,*y*), where vertical stripes are eliminated. This approach, here further refereced as “wavelet filter”, has already been proposed by Torres and Infante [1] and applied to satellite images, which are sometimes heavily damaged by striping noise due to gain and offset differences among the elements of the detector arrays. The authors point out that the truncated wavelet representation *f*′(*x*,*y*) still approximates *f*(*x*,*y*) in a least squares sense (Parseval theorem). However, this kind of stripe elimination is obviously associated with a notable loss of structural information other than the striping noise. The left image of Fig.4 shows the denoised version of the left image of Fig.3 obtained by using this approach. To suppress all stripe residuals, the vertical details *c*
_{v l,m,n} at the first five decomposition levels *l* ∈ {0, ⋯, 4} had to be completely eliminated, yielding a poor quality image. For the satellite images in reference [1], truncation of only the first three levels was sufficient for a successful destriping, leading to a considerably better image quality.

### b) FFT filtering

In the spatial frequency domain *F*(*x*̂,*ŷ*) of *f*(*x*,*y*), ideal vertical stripes include high frequency parts in the horizontal direction *x*̂, while in the vertical direction *ŷ* after the 2D FFT, the stripe offset yields Dirac delta functions *δ*(*ŷ*) (i.e. unit impulse for discrete functions) at all *x*̂. In other words, there are no frequency components stemming from vertical stripes in *ŷ* ≠ 0. Consequently by eliminating the Fourier coefficients *F*(*x*̂,*ŷ*) of *f*(*x*,*y*) at all *x*̂ for *ŷ* = 0, the entire information arising from ideal vertical stripes will be erased. The concept of damping the values of these coefficients is in fact the principal concept [4, 27] or an important part [3, 28] at the base of several existing stripe removal procedures.

For this purpose, a simple approach in the Fourier space is the application of a bandpass filter around *ŷ* ≈ 0. For instance, a selective damping of *F*(*x*̂,*ŷ*) on the *x*̂-axis can be obtained by multiplication of the FFT coefficients with a Gaussian function *g*(*x*̂,*ŷ*),

The value of *σ* determines the width of the filter in *ŷ*-direction and is selected according to the expected deviation from the vertical of the stripes in *x*-direction in the spatial domain. Note that for this specific FFT filter, the Fourier coefficients at *F*(*x*̂ = 0,*ŷ*) represent average offsets over the entire image width *x* and therefore should not be removed. For this reason, damping is applied for ∀*x*̂ ≠ 0 only.

This type of FFT stripe filter, here further referenced as “FFT filter”, has been applied to the Lena image in Fig.3 (left). Since the two colored stripe artifacts are exactly vertical, in this case a small *σ* → 0 has been selected so that essentially, *F*(*x*̂,*ŷ*) was set to zero at *ŷ* = 0 only. The resulting destriped image is displayed in the centre of Fig. 4.

### c) Combined wavelet-FFT filtering

The new filtering technique proposed in this paper achieves, in the first step, a tight condensation of the stripe information by combining the wavelet and FFT transforms. It then performs a damping of the relevant coefficients by using a Gaussian function.

In the first step, the original image *f*(*x*,*y*) (e.g. the left image in Fig.3) is wavelet decomposed into *W* = {*c*
_{l L,m,n}, *c*
_{h l,m,n}, *c*
_{v l,m,n}, *c*
_{d l,m,n}},*l* ∈ {1,⋯,*L*} in order to separate the structural information into horizontal, vertical and diagonal details bands at different resolution scales.

Subsequently, the bands containing the stripe information (e.g. *c*
_{v l,m,n} for vertical stripes, see portions enframed in red in the centre image of Fig.3) are FFT transformed to further tighten the stripe information into narrow bands (e.g. to the *x*̂-axis in the FFT domain in case of vertical stripes, see enframed horizontal regions in the right image of Fig.3). Vertical 1D FFT transforms of all columns in the subregions *c*
_{v l,m,n} are sufficient to condensate the stripe information to the *x*̂-axis. A further transform in the horizontal direction does not result in an additional compaction. For the FFT filter in section 2.2, the additional transform in *x*̂ direction is necessary to be able to condense the global offset values around *F*(0,0). This second transform is instead not necessary for the combined wavelet-FFT filter, because for *c*
_{v l,m,n}, this global offset value is close to zero since these coefficients contain high frequency information only. For the proposed filter, the global offset values are practically entirely detached to the coefficients *c*
_{l L,m,n}.

The stripe information condensed in such a way is then eliminated by multiplication with a Gaussian damping function *g*(*x*̂,*ŷ*) according to Eq.9. Finally, the destriped image *f*′(*x*,*y*) is reconstructed from the filtered coefficients. The result is displayed in Fig.4 on the right image.

The proposed new filter is further referenced as “wavelet-FFT filter”.

#### 2.3. Algorithm

The core algorithm performing the new “wavelet-FFT filter” is presented in Fig.5. It is written in Matlab command language and envisaged for the case of vertical stripe artifacts removal.

This algorithm, in addition to the input image *’ima’*, employs three parameters for the filtering process: the highest decomposition level *L* (*’decNum’*), the wavelet type (*’wname’*) and the damping factor *σ* (*’sigma’*) from Eq.9. In our approach, *σ* is kept uniform for each decomposition level *l* ∈ {0,⋯,*L*}. Alternatively, *σ* could also be adjusted according to the resolution at each *l*.

The algorithm consists of three distinct parts. At lines 4-6, the wavelet decomposition is calculated by recursive splitting of the original image and of the low resolution coefficients *c*
_{l l,m,n} from the former decomposition level *l* – 1 into {*c*
_{l l,m,n}, *c*
_{h l,m,n}, *c*
_{v l,m,n}, *c*
_{d l,m,n}}. The boundary conditions are met by symmetric mirroring. For all decomposition levels *l* ∈ {0, ⋯, *L*}, the vertical details coefficients *c*
_{v l,m,n} are then Fourier transformed to the coefficients *ĉ*
_{v l,m,n} by a columnwise 1D FFT (lines 11-12). Subsequently at lines 15-16, these coefficients are multiplied with a Gaussian function, to eliminate those close to the *x*̂-axis in the Fourier domain. After such damping, the coefficients *ĉ*′_{v l,m,n} are transformed back to the wavelet space (line 19). Finally, the destriped image is reconstructed at lines 23-27 from the refined wavelet coefficients *c*′_{v l,m,n}. The resizing command at line 25 is necessary for the Matlab function ’*idwt2*’ to enforce the correct dimensions for the gradually reconstructed *c*′_{l l,m,n}. The resulting image is the destriped version of *’ima’*.

## 3. Discussion of wavelet-Fourier filter characteristics

#### 3.1. Conceptual objection: Fourier transform of the wavelet domain?

Considering the presented wavelet-FFT filtering concept, the question about the sense of a Fourier transform within the wavelet domain arises. From a philosophical standpoint, a Fourier transform of wavelet coefficients has a rather unclear meaning. In the context of the stripe filter, the wavelet decomposition should, however, rather be interpreted from the signal analysis point of view, which is the filtering of a signal with low and highpass filters. This filtering approach allows to specifically select for further “curing” only those bands affected by stripe artifacts, while leaving the remaining bands unchanged. The special property of this filtering procedure, common to orthogonal transforms, is its reversibility. Using the filtering interpretation, the detail bands can be thought as appearing in the spatial domain without need for the misleading term “wavelet domain”. Accordingly, the subsequent Fourier transform converts them from the spatial into the frequency space only.

#### 3.2. Tight condensation of stripe information

The first innovative feature characterizing the high performance of the proposed wavelet-FFT filter is the tight condensation of vertical (or horizontal) stripes into strictly isolated *ĉ*
_{v l,m,n} subsets. This characteristic is required to achieve the requirements (a) - (c) in Section 1, and is visualized in Fig. 6. The image to the left shows a pattern consisting exclusively of perfect vertical colored stripes. The centre image displays the resulting wavelet coefficients after a multiscale decomposition up to level *L* = 4 obtained using the Haar wavelet (DB1). The entire image information is now condensed in the vertical details band *c*
_{v l,m,n} and the remaining low pass band *c*
_{l L,m,n}. The latter can be further decomposed until the stripe characteristics disappear, or due to the successive downsampling, until some few coefficients which are containing the overall background offset is remaining only. All other coefficients *c*
_{h l,m,n} and *c*
_{d l,m,n} are zero except for some values at the boundaries which occur due to border effects. After FFT transform of *c*
_{v l,m,n} for all *l* yielding *ĉ*
_{v l,m,n} (right image of Fig.6), the only non-zero coefficients are localized on the abscissa and on the unchanged low pass band *c*
_{l L,m,n}.

Vice versa, these non-zero coefficients *ĉ*
_{v l,m,n=0} contain pure stripe information only. This fact is demonstrated in Fig.7, where four artificial sets of coefficients in the FFT domain have been created by assuming some random complex values for *ĉ*
_{v l,m,n} =𝓝 at some horizontal lines 𝓝, while assuming zero values for all remaining coefficients. Such randomly generated coefficient sets are then reconstructed to the spatial domain. Fig.7 displays the arising artificial images for 𝓝 = 0 (leftmost image), 𝓝 = ±1, (2nd leftmost), 𝓝 = ±2, (2nd rightmost), and at 𝓝 = ±3 (rightmost). Obviously, only the image to the left, which was built from *ĉ*
_{v l,m,n=0} is exclusively consisting of perfect stripes as defined in Section 2.1, while for increasing |𝓝| > 0, the pattern progressively deviates from the strict stripe definition.

With respect to the wavelet-FFT filter, the range {0, ⋯, 𝓝} of damped coefficients directly depends on *σ* (Eq.9): while for small *σ* → 0, only stripes in the strict sense (i.e. at 𝓝 = 0) will be erased, partial striping will be damped as well for higher *σ* (i.e. particularly for *σ* > 1).

For instance, the Lena image of Fig. 3 (left) and of Fig. 9 (top left) are both ravaged by perfect stripes in the sense of Section 2.1, wherefore damping was performed at *σ* → 0. Larger *σ* would not contribute to an improvement in destriping, but rather cause original image information to be further eliminated. In contrast, imperfect stripes (e.g. stripes that are not perfectly straight, or with an offset value *A*, which is not constant over the whole length) must be eliminated at increased *σ*. For example, the stripes in Fig. 11 (top left — further explanation is provided in Section 4.1) are slightly crooked, and furthermore contain varying offset values. Thus, filtering has been performed at *σ* = 7, whereby the overall energy changes (see Section 3.8) have got considerably larger and cannot be exclusively ascribed to the removed stripes.

#### 3.3. Separation of stripe information

The destriping procedure, i.e. Fourier filtering and coefficient damping, is applied only to some selected detail bands of the wavelet coefficients. This is the key reason for the high performance of this wavelet-FFT filter in removing stripe artifact. In fact, in this way, each processed detail band holds structural information within a bounded scale only therefore enabling stripe filtering to be optimized to the specific scale. In addition, a large number of detail bands remains completely unchanged guaranteeing the preservation of all structural features outside the detail bands affected by stripe artifacts.

After each single wavelet decomposition step, the remaining low pass coefficients *c*
_{l l,m,n} still contain the self-similar full image information, yet at a lower resolution. In contrast, all other bands only contain details, i.e. high pass information for getting increased resolution. Hence, the mean image values at lower resolution are still completely preserved when coefficients from the details bands are removed. This is not the case if coefficient damping is performed in the Fourier space only (see Section 2.2), since these coefficients also contain quantitative information about regional offset values, which will be affected during the filtering procedure.

The power of the strict discrimination between stripe information and other structures is demonstrated in the extreme example in Fig.8. The Lena image (1) was first scaled by a value of ≈ 1/50 in order to further fade out the original structures relative to the stripes. An image consisting only of colored horizontal and vertical stripes (2) was then added to it. In the resulting image (3), the stripe noise completely covers all other information. The image corruption is so strong that the Lena or any other structure in the original image cannot be visually recognized.

Successive destriping in the horizontal and vertical directions using the proposed wavelet-FFT filter at *σ* → 0, *L* = 8 and DB42 yields image (4), which shows a high preservation of the structural information of the original Lena. Major discrepancies are observed in the color structure. Shifts in the lower spatial frequencies are responsible for these differences. Such shifts arise from the mean values of the heavy stripe pattern, which are not easily removed and overcast the mean values of the Lena image (especially due to its scaling and therefore minimal contribution to the overall image). The difference between the original Lena (1) and the destriped version of the Lena ravaged by stripes (4) is displayed in image (5). The residual image shows no evidence of any structural information related to the original Lena motive, except for some horizontal or vertical stripe-like features (e.g. the doorframe to the left). It therefore confirms the high structural preservation capabilities of the new proposed wavelet-FFT filter, even when stripe contamination is so strong that the original image is not visually recognizable.

#### 3.4. Computation time

An outstanding advantage of the proposed filter is its minimal computational effort. The proposed algorithm consists of a discrete wavelet transform (DWT), a discrete Fourier transform, generally implemented by a fast Fourier transform (FFT), and damping of some coefficients, which is basically a scalar multiplication with a Gaussian function (Section 2.2). All three tasks can be performed with *O*(*N*) operations (*N*: number of image pixels), therefore the entire filtering process runs in linear time. Since the wavelet decomposition is performed by means of 1D wavelet filtering, the computation time depends on the size *l _{w}* of the discrete wavelets with

*O*(

*l*).

_{w}#### 3.5. Translation invariance

In the proposed filtering process, a dyadic, decimated wavelet analysis is performed, followed by an FFT of the coefficients corresponding to vertical details. This process is not translation invariant, that is, the result of a translated signal is not the translated result. However, since the stripe removal process involves a damping of the coefficients *ĉ*
_{v l,m,n}, which is uniform along the entire abscissa (Eq.9), the striping information is completely removed for all locations in the case of ideal stripes (where *g*(*x*̂,0) = 0), or evenly damped for each class of non-ideal stripes (Fig.7).

#### 3.6. Dyadic detail levels

The dyadic scales of the detail levels offer the possibility to selectively apply the destriping procedure to artifacts of specific widths. In fact, the width of removed stripes is associated with the levels of decomposition subjected to filtering (Section 2.2). If the width range of possible stripe artifacts in an image is narrow, the decomposition levels submitted to the damping procedure can be reduced and unnecessary loss of valid information avoided.

In Fig.9, the original Lena image (bottom rightmost image) was ravaged by seven stripes of different colors and sizes (24, 16, 12, 8, 4, 2, 1 pixel, top leftmost image). Subsequently, it was filtered using the DB15 wavelet and *σ* → 0, by assuming different highest decomposition levels *L* = 0,⋯,7. The remaining *c*
_{l L,m,n} for *l* > *L* were left untouched.

While the smallest stripe to the right of the image (blue) is not visible any more after *L* = 3 (top rightmost), green shadows of the largest stripe to the left can still be well perceived at *L* = 5.

Fig. 9 also reveals that in some cases, specifically for images with blurred stripes, it might even be favourable to perform the filtering on a subset of the decomposition numbers *l* = {*l _{a}*, ⋯,

*l*} only. In particular, an image with the same type of stripes as the top rightmost image would have to be filtered at the decompositions

_{e}*l*= {4, ⋯, 7} only.

#### 3.7. Assessment of filter quality

For reliable validation of the filter quality, both, qualitative and quantitative aspects are decisive. Qualitatively, the destriped images need to be free from stripes, while all other image details have to be preserved. Quantitatively, local mean values of the filtered image away from stripes must be maintained — an imporant requirement for quantitative image analysis.

We propose two different methods for the qualitative and quantitative assessment of the performance of our new filter. Visual inspection is probably the most adequate and reliable approach to evaluate the destriping success and to validate the detail preservations. Though, the human eye highly adapts to ambient changes and thus is not reliable to assess the quantitative stability. An objective well known measure in signal processing for a quantitative evaluation of the results is the total energy of a signal *s _{x}*,

*x*∈ {0, ⋯,

*N*} (from Parseval’s theorem), which is defined as the sum of squared moduli

In the ideal case, energy changes should arise solely from stripe suppression. In practice, this is not the case, whereby the combination of visual rating with a minimal relative change of *ε _{s}* between the original and the processed image is supposed to be a robust technique for the assessment of the filter performance and the result quality.

The loss of the energy can be expressed by the energy of the difference of the original (*s _{o}*) and the filtered (

*s*) images, relative to the original one, resulting in the relative mean square error (MSE):

_{f}This relative energy metric (Eq. 11) was applied to the images of Fig.4 and the relative energy difference before and after filtering calculated. The resulting energy changes are:

wavelet: 2.275 [%] | FFT: 2.854 [%] | wavelet-FFT: 0.597 [%] |

The Lena image processed using the proposed wavelet-FFT stripe filter shows the better visual quality while achieving four times less energy change.

#### 3.8. Choice of wavelets

For the same example as in the previous section, different types and sizes of wavelets have been applied. The use of Daubechies (DB) wavelets [26] shows visually favourable results even at low energy changes *ε _{s}* (Section 3.7). The reason for this has not been explored and may presumably be accompanied with the wavelet smoothness. In Fig. 10, the energy changes in [%] after filtering with Daubechies wavelets of different size are displayed. Three images have been treated: the Lena (see Section 1), cement particles (see Section 4.1), and a sinogram (see Section 4.2). The values for the damping coefficients

*σ*and for the highest decomposition level

*L*have been selected as small as possible, but large enough to make visible stripe remnants disappear. This goal was achieved at damping coefficients of

*σ*= 1.5 (Lena),

*σ*= 7.0 (cement),

*σ*=1.5 (sinogram), and for highest decomposition levels of

*L*= 4,

*L*= 8,

*L*= 4, respectively. While these values have been kept constant, the wavelet filters have been varied from DB1 to DB43.

For all three images, the energy changes *ε _{s}* are tendentially decreasing with increasing wavelet size. The use of the larger wavelet sizes however results in increasing calculating effort. By using a wavelet size between 10 and 15, a reasonable compromise between minimal

*ε*and low processing time was found. Visual assessments agree with these findings. Since the visual differences are subtle, no further illustrations are provided.

_{s}The largest *ε _{s}* occur for the cement image. This is due to (1) the large damping coefficient and (2) highest decomposition level, and (3) to the heavy stripe artifacts.

## 4. Applications

#### 4.1. Waterfall effect in focused ion beam nanotomography (FIB-nt)

FIB-nt is a technique that utilises dual (electron and ion) beam microscopy for 3D imaging with nanometric resolution. The 3D datasets are produced by a repetitive ion milling and electron imaging process. Further details can be read from the literature [29].

A severe problem in FIB-nt is the possible contamination of electron images by vertical stripes, the so called waterfall effect. This artifact is caused by the occurrence of material phases in the sample microstructure that exhibit a variable resistance to ion milling. It is especially pronounced in cryo FIB-nt [30], where inhomogeneous regions occur in the metalorganic protective layer used in FIB-nt to avoid charging effects. The example in Fig. 11 (top left) shows a single slice from a 3D volume of unhydrated particles from cement paste with a total size of 17.1 × 14.6*μm*, recorded by cryo FIB-nt from Holzer et. al. [30, 31].

We have destriped the FIB-nt image using three types of filters (Section 2.2): wavelets only (top right image), FFT only (bottom left), and wavelet-FFT (bottom right). Visual inspection clearly shows the superior performace of the combined wavelet-FFT filter. The respective relative energy changes are 0.57 % for the wavelet-FFT filter, 1.40 % for the wavelet filter, and 1.52 % for the FFT filter. The optimum results with the combined filter were achieved when the DB42 wavelet was used up to the highest decomposition level *L* = 8, and with a damping coefficient of *σ* = 7. To minimize the overall structural loss, the stripe removal with the wavelet-only method was performed at *L* = 3 only.

Waterfall artifacts in FIB-nt represent an especially challenging application, as observed stripes deviate from the ideal description of Eq.2.1, that is, they are (1) not perfectly vertical, and (2) their offset values *A* are fluctuating. Hence, the original image in Fig.11 (top left) is a good example for showing the effect of the damping coefficient *σ*, and the highest decomposition level *L*. In Fig.12, different *σ*= 1,2,4,8 (rows) and *L* = 2,4,6,8 (columns) have been applied.

The results clearly show that unlike for perfect stripes, where *σ* → 0 is sufficient (see Fig.9) and the filtering must be adjusted according to *L* only, in this more general case the values for *σ* need to be larger for the complete elimination of imperfect stripes. Due to some broad artifacts, stripes don’t completely disappear below *L* = 8.

#### 4.2. Ring artifacts in synchrotron based X-ray Tomographic Microscopy

### a) Origin

In synchrotron based X-ray Tomographic Microscopy [32], a sample is rotated in a stationary parallel monochromatic beam and 2D radiographs for many angular positions over 180° are recorded with a high-resolution area detector (e.g. CCD camera). The information of the 2D projections is usually rearranged in sinograms before reconstruction. Each sinogram contains the information of a particular horizontal detector line for all angular positions. A single detector element records, for every angular position, the intensity of one particular ray at one particular position in the beam. The information of this detector pixel back-projects, in the reconstructed slice, onto a half-circle centered on the rotation axis. Defective detector elements (e.g. dead pixels in a CCD) with non linear responses to incoming intensity will therefore appear in the reconstructions as sharp rings with a width of one pixel (e.g. Fig.13(c), see [33]). Similar artifacts also arise from dusty or damaged scintillator screens. Miscalibrated detector pixels, e.g. due to beam instabilities not completely taken into account by a normalization correction, give instead rise to wider and less marked rings (e.g. Fig.20(a)).

Minimization of ring artifacts by using adequate scanning protocols, high quality scintillator screens and CCD cameras is possible. It is, however, difficult to completely avoid such artifacts and therefore achieve highest quality reconstruction solely by experimental measures. Software routines are required to remove the remaining distortions and thus allow full qualitative and quantitative exploitation of the acquired information.

### b) Removal of sharp and marked ring artifacts

Dead pixels in a CCD chip and damaged scintillator screens are responsible for sharp and marked vertical stripes in sinograms, which back-project to half circles (ring artifacts) in tomographic reconstructions. These artifacts are clearly visible in the example in Fig.13 and Fig.14. We used the wavelet-FFT filter to attenuate these features and improve the reconstruction quality. We tested different parameter combinations (wavelet type, damping coefficient *σ* and decomposition level *L*) aiming at the best image quality, assessed visually and attested with line profiles, and a minimal relative energy change. We experimented in particular with different wavelet types and found that with DB25 the best compromise between image quality and energy loss is achieved. For Daubechies wavelets of a smaller size (down to DB1), larger discrepancies in the line profiles are observed; for larger wavelet sizes (up to DB40), artifacts in the reconstructions appear. A decomposition level *L* larger than 6 needs to be avoided for this example, because details of the image were also starting to be filtered out. Decomposition level smaller than 5 are instead not enough to removed all ring artifacts. We therefore chose *L* = 5. Based on visual inspection, the decomposition level l=0 was left unchanged. In fact, FFT filtering all decomposition levels leads to a ”windmill” artifact in the center of the image, judged to have a larger disturbing impact than the remaining faint tiny rings when l=0 is untouched. Finally for *σ* smaller than 2.4, not all ring artifacts were properly removed, larger *σ* resulted in a significant increase in relative energy loss. The images shown in Fig.13(b), (c) and Fig.14(b) have therefore been obtained with a DB25 wavelet and filtering of the decomposition levels 1 to 5 with *σ* = 2.4. The relative energy change is 0.018% only. Visually, the quality improvement is clear. Vertical stripes in the sinograms (Fig.13(b)) have disappeared as well as the marked rings in the reconstructed slice (Fig.13(d)). Removal of these artifacts from the reconstructions reveals new structural details (e.g. arrows in Fig.14(b)), which were only hardly discernible in the corrupted image. Horizontal and vertical line profiles through the center of the image (Fig.15) confirm the minimal changes to the reconstructed slices for regions not affected by the ring artifacts.

Similar results are obtained in a second example, where the artifacts consist of a wider strong ring (Fig.16(a)). In this case, we obtained the best results with the DB30 wavelet and FFT filtering of the wavelet components 1 to 5 (Fig.16(b)) or 0 to 6 (Fig.16(c)) with *σ* = 2.4. A wavelet decomposition up to level *L* = 6 allows to completely remove the ring artifact (Fig.17(c)) unraveling the hidden structure, though at the expenses of little (windmill) deformation and grey level drop at the center of the image (Fig.18(c), 19(a) and 19(b)). Choosing a *L* lower than 6 has less influence on the grey level at the image center (Fig.18(b), 19(a) and 19(b)), leaves though a faint dark shadow where the ring originally was (Fig.17(b)), nonetheless revealing the lost features. The related relative energy change is small (0.17% and 0.14%, respectively): as observed in the line profiles (Fig.19) changes are minimal away from the image center and the artifact.

### c) Removal of wide and faint ring artifacts

Miscalibrated detector pixels, e.g. due to beam instabilities not completely taken into account by a normalization correction, give rise to wide and faint rings. An extreme example is shown in Figs.20 and 21: the original image (Fig.20(a)) is so strongly contaminated by ring artifacts, that the underlying structure is hardly discernible and a quantitative analysis completely impossible. Sample details can however be unraveled by the wavelet-FFT filter. In particular, using a DB15 wavelet to filter the 0 to 4 components with a *σ* = 2.4, encouraging results have been obtained (Fig.20(b)). The choice of these parameters permits the removal of most rings from the reconstructed slice, with minimal distortion at the image center. Higher order wavelets have a stronger effect on the center of the image, while at least the first 5 components need to be filtered to be able to remove the observed rings. Taking more components into account results into a larger energy loss, without an appreciable improvement in the image quality. Due to the strong image corruption, the observed relative energy change of 18% for the chosen optimal parameters is significantly higher than for the previous examples.

#### 4.3. Horizontal stripe artifacts in phase projections obtained with grating based X-ray phase contrast imaging

Phase sensitive X-ray imaging methods provide increased contrast over conventional absorption based imaging, and therefore new and otherwise inaccessible information. The Differential Phase Contrast (DPC) imaging technique [36] uses a grating interferometer and a phase-stepping approach to separate the phase information from other contributions. Projections of the phase information obtained through integration of DPC radiographs show characteristic horizontal stripes, arising from beam instabilities and sub-optimal background normalization. An example is shown in Fig.22(a).

Such artifacts can successfully be minimized with the new filter. Best results have been obtained with the DB30 wavelet *L* = 3 and *σ* = 3 (Fig.22(b)). The use of other wavelet lengths did not result in an appreciable difference in the images. The chosen wavelet type produced the minimal change in the relative energy confirmed by the smallest discrepancies between the original and filtered projection in line profiles. Damping components of higher order than *L* = 3 caused a large change in the relative energy not justified by a significant improvement in the image quality. Finally, an increase of *σ* was connected to a smoothing of the line profile and fading of small stripes. The chosen parameters represent a good compromise between optimal artifact suppression and high preservation of the original image information. The relative energy change was 0.28%.

#### 4.4. Conclusions

The presented combined wavelet-FFT filter represents a powerful approach for a wide range of stripe artifact removal problems. It is designed for ideal stripe artifacts consisting of straight horizontal or vertical constant offsets traversing the entire image. Tight contraction and classification of such stripe information is achieved by successive wavelet and Fourier transforms guaranteeing high preservation of all structural information different from ideal stripes during the filtering procedure.

Three free parameters allow to control the filter behavior and extend its applicability to artifacts deviating from ideal stripes: (1) the choice of the wavelet, (2) the highest decomposition level *L*, and (3) the damping factor *σ*. Optimal results have been achieved for large Daubechies wavelets (i.e. large number of vanishing moments, at least larger than DB3). The optimal choice of the highest decomposition level *L* depends on the largest stripe width, as well as on the sheerness of the slopes at the stripe borders. As long as satisfactory stripe removal can be achieved, small *L* are preferred in order to preserve the quantitative overall offset values as best as possible. The optimal damping factor *σ* depends on how ideal the stripes are in terms of straightness and offset constancy. While *σ* must be increased for slightly distorted stripes and for stripes with varying offset, it can instead be lowered to *σ* = 0 for ideal artifacts.

The outstanding benefits of this new filter are 1) the possibility to adjust the parameters to selectively choose the bandwidth of the stripe information, which is to be removed. Moreover 2), the filter is softly responding to any parameter changes, therefore exhibiting a good stability. Finally 3), the filter performs at linear calculation time.

In addition to qualitative visual assessment, the filter performance was validated quantitatively by evaluating the changes of the signal energy before and after filtering. If only few stripe artifacts are present, a relatively small loss of signal energy is usually achieved, indicating that most of the remaining image information is being preserved in a quantitative sense.

Several examples of color and gray level images containing both pronounced and subtle horizontal and / or vertical stripes have been presented in the field of tomographic microscopies, in particular in cryo FIB-nt and synchrotron microtomography. These various examples demonstrate a wide range of applicability and a good adaptability of the filter to different situations and requirements. Compared to previous approaches, the new technique shows a very high selectivity.

A java application of the filter can be freely downloaded from ftp://ftp.empa.ch/pub/empa/outgoing/BeatsRamsch/stripeFilter/xStripes.jar

## Acknowledgments

The authors would like to acknowledge the comments and discussions with Pietro Lura (Lab of Concrete and Construction Chemistry, Empa).

## References and links

**1. **J. Torres and S.O. Infante, “Wavelet Analysis for the Elimination of Striping Noise in Satellite Images,” Opt. Eng. **40**, 1309–1314 (2001). [CrossRef]

**2. **M.J. Oimoen, “An Effective Filter for Removal of Production Artifacts in U.S. Geological Survey 7.5-Minute Digital Elevation Models,” Proc. 14th Intl Conf. Appl. Geol. Remote Sens., November, Las Vegas, N USAV, 311–319 (2000).

**3. **J.G. Liu and G.L.K. Morgan, “FFT Selective and Adaptive Filtering for Removal of Systematic Noise in ETM+ Imageodesy Images,” IEEE Trans. Geosci. Remote Sens. **44**, 3716–3724 (2006). [CrossRef]

**4. **J. Chen, Y. Shao, H. Guo, W. Wang, and B. Zhu, “Destriping CMODIS Data by Power Filtering,” IEEE Trans. Geosci. Remote Sens. **41**, 2119–2124 (2003). [CrossRef]

**5. **J. Chen, H. Lin, Y. Shao, and L. Yang, “Oblique Striping Removal in Remote Sensing Imagery Based on Wavelet Transform,” Int. J. Remote Sens. **27**, 1717–1723 (2006). [CrossRef]

**6. **M. Albani and B. Klinkenberg, “A Spatial Filter for the Removal of Striping Artifacts in Digital Elevation Models,” Photogramm. Engineering & Remote Sens. **69**, 755–765 (2003).

**7. **M. Boin and A. Haibel, “Compensation of Ring Artefacts in Synchrotron Tomographic Images,” Opt. Express **14**, 12071–12075 (2006). [CrossRef] [PubMed]

**8. **Z. Cai, “Ringing Artefact Reduction Using Adaptive Averaging Filtering,” Proc. ISCE, UK, June, 156–159 (2004).

**9. **P. Rakwatin, W. Takeuchi, and Y. Yasuoka, “Stripe Noise Reduction in MODIS Data by Combining Histogram Matching with Facet Filter,” IEEE Trans. Geosci. Remote Sens. **45**, 1844–1856 (2007). [CrossRef]

**10. **F. Tsai, S. Lin, J. Rau, L. Chen, and G. Liu, “Destriping Hyperion Imagery Using Spline Interpolation,” Proc. 26th Asian Conf. Remote Sensing, November, Hanoi, Vietnam (2005).

**11. **K. Arrell, S. Wise, J. Wood, and D. Donoghue, “Spectral Filtering as a Method of Visualising and Removing Striped Artefacts in Digital Elevation Data,” Earth Surface Processes and Landforms **33**, 943–961 (2007). [CrossRef]

**12. **X. Tang, R. Ning, R. Yu, and D. Conover, “Cone Beam Volume CT Image Artifacts Caused by Defective Cells in X-Ray Flat Panel Imagers and the Artifact Removal Using a Wavelet-Analysis-Based Algorithm,” Med. Phys. **28**, 812–825 (2001). [CrossRef] [PubMed]

**13. **G. Kowalski, “Suppression of Ring Artefacts in CT Fan-Beam Scanners,” IEEE Trans. Nucl. Sci. **NS-25**, 1111–1116 (1978). [CrossRef]

**14. **A. Lyckegaard, G. Johnson, and P. Tafforeau, “Correction of Ring Artefacts in X-ray Tomographic Images,” 3D IMS, 1st Conf. 3D-Imaging on Materials and Systems, sept.8–12, Carans-Maubuisson, 101 (2008).

**15. **M. Axelsson, S. Svensson, and G. Borgefors, “Reduction of Ring Artifacts in High Resolution X-Ray Microtomography Images,” Pattern. recog. **4174**, 61–70 (2006). [CrossRef]

**16. **J. Sijbers and A. Postnov, “Reduction of Ring Artefacts in High Resolution Micro-CT Reconstructions,” Phys. Med. Biol. **49**, N247–N253 (2004). [CrossRef] [PubMed]

**17. **C. Raven, “Numerical Removal of Ring Artifacts in Microtomography,” Rev. Sci. Instrum. **69**, 2978–2980 (1998). [CrossRef]

**18. **G.T. Herman, “Image Reconstruction from Projections-the Fundamentals of Computed Tomography,” Academic Press New York (1980).

**19. **D.L. Donoho and G. Flesia, “Digital ridgelet transform based on true ridge functions,”
Beyond Wavelets, J. Stoeck-ler, and G. Welland (Eds), Academic Press, 1–33 (2001).

**20. **E. J. Candés and D. L. Donoho, “Curvelets - a surprisingly effective nonadaptive representation for objects with edges,” Curves and Surfaces,
L. L. Schumaker et al. (eds), Vanderbilt University Press, Nashville, TN, 1–10 (1999).

**21. **S. Foucher, V. Gouaillier, and L. Gagnon, “Global semantic classification of scenes using ridgelet transform,” IS&T/SPIE Symposium on Electronic Imaging: SPIE Human Vision **5292**, San Jose, 402–413 (2004).

**22. **C. Christopoulos and T. Ebrahimi, “The JPEG2000 Still Image Coding System: An Overview,” IEEE Trans. Consumer Electronics **47**, 1103–1127 (2000). [CrossRef]

**23. **S. Mallat, “A Wavelet Tour of Signal Processing,” Academic Press, Second Edition, ISBN: 978-0-12-466606-1, 1-end (1999).

**24. **S.G. Mallat, “A Theory for Multiresolution Signal Decomposition: The Wavelet Representation,” IEEE Trans. Pattern Analysis and Machine Intelligence **11**, 674–693 (1989). [CrossRef]

**25. **W. Bäni, “Wavelets. Eine Einführung für Ingenieure,” Oldenbourg, ISBN: 978-3-486-57706-8, 1-end (2005).

**26. **I. Daubechies, “Ten lectures on Wavelets,” SIAM, Philadelphia, PA (1992).

**27. **Z. Zhang, Z. Shi, W. Guo, and S. Huang, “Adaptively Image De-striping through Frequency Filtering,” ICO20: Opt. Inf. Proc., Proc. SPIE **6027**, 60273V (2006). [CrossRef]

**28. **Z. Wang and Y. Fu, “Frequency-domain Regularized Deconvolution for Images with Stripe Noise,” ICIG Proc. 4th Intl Conf. Image and Graphics, 110–115 (2007).

**29. **L. Holzer, F. Indutnyi, P. Gasser, B. Münch, and M. Wegmann, “Three-Dimensional Analysis of Porous BaTiO_{3} Ceramics Using FIB Nanotomography,” J. Microsc. **216** (1), 84–95 (2004). [CrossRef] [PubMed]

**30. **L. Holzer, Ph. Gasser, A. Kaech, M. Wegmann, A. Zingg, R. Wepf, and B. Münch, “Cryo-FIB-nanotomography for Quantitative Analysis of Particle Structures in Cement Suspensions,” J. Microsc. **227**, 216–228 (2007). [CrossRef] [PubMed]

**31. **A. Zingg, L. Holzer, A. Kaech, and F. Winnefeld, “The Microstructure of Dispersed and Non-dispersed Fresh Cement Pastes-New In-sight by Cryo-Microscopy,” Cement and Concrete Research **38**, 522–529 (2008). [CrossRef]

**32. **A.C. Kak and M. Slaney, “Principles of computerized tomographic imaging,” Society for Industrial and Applied Mathematics, Philadelphia, PA (2001).

**33. **L.A. Shepp and J.A. Stein, “Simulated reconstruction artifacts in computerized X-ray tomography,” Reconstruction Tomography in Diagnostic Radiology and Nuclear Medicine, University Park Press, Baltimore, 33–48 (1977).

**34. **M. Stampanoni, A. Groso, A. Isenegger, G. Mikuljan, Q. Chen, D. Meister, M. Lange, R. Betemps, S. Henein, and R. Abela, “TOMCAT: A beamline for TOmographic Microscopy and Coherent rAdiology experimenTs.,” Synchrotron Radiation Instrumentation **879**, 848–851 (2007).

**35. **D. Prêt, F. Villiéras, I. Bihannic, S. Sammartino, L. Michot, M. Pelletier, and M. Stampanoni “Influence of hydration on montmorillonite organization : a multiscale study based on synchrotron X-ray tomography, diffraction and small angle scattering of neutrons,” in prep. for Langmuir.

**36. **T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “Quantitative X-ray phase imaging with a grating interferometer,” Opt. Express **13**, 6296–6304 (2005). [CrossRef] [PubMed]