## Abstract

We introduce a new X-ray speckle-vector tracking method for phase imaging, which is based on the wavelet transform. Theoretical and experimental results show that this method, which is called wavelet-transform-based speckle-vector tracking (WSVT), has stronger noise robustness and higher efficiency compared with the cross-correlation-based method. In addition, the WSVT method has the controllable noise reduction and can be applied with fewer scan steps. These unique features make the WSVT method suitable for measurements of large image sizes and phase shifts, possibly under low-flux conditions, and has the potential to broaden the applications of speckle tracking to new areas requiring faster phase imaging and real-time wavefront sensing, diagnostics, and characterization.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Modern synchrotron radiation facilities and X-ray free electron lasers (XFEL) provide high-brightness X-ray beams that enable advanced techniques in a variety of applications in medical, material, biological, and physical sciences [1,2]. Among all applications, X-ray phase sensing techniques have been attracting increasing attention in phase-contrast imaging as well as in at-wavelength metrology for X-ray optics, because of their high sensitivity, high spatial resolution, and high penetration power through samples [3–5]. During past decades, many phase imaging techniques were developed based on either propagation method, such as X-ray coherent diffraction imaging (XCDI) [6] , or phase gradient methods, such as X-ray grating-based imaging (GBI) [7–11] and speckle-based imaging (SBI) method [12–14]. Compared with XCDI methods using far-field diffraction, the GBI and SBI methods have a lower spatial resolution but a much larger field of view, simpler setup, and lower requirements on X-ray flux and coherence.

The GBI method has been successfully applied to wavefront sensing, phase imaging, and X-ray at-wavelength metrology [7,15–18]. High phase sensitivity becomes feasible with proper calibration and optimized data processing [19]. However, the fabrication of high-quality high-aspect-ratio gratings is challenging, and relatively high beam coherence is required to achieve high sensitivity. To overcome these bottlenecks, the near-field SBI method was developed and successfully demonstrated on sources ranging from synchrotrons to table-top X-ray sources [20–23].

Several variations of the SBI method have been developed using different data acquisition procedures and reconstruction algorithms [12,14,22,24–29]. Experimentally, SBI measurements can be categorized into two types, the so-called single-image and speckle scanning methods. The former represents a data set with a single image with the sample plus a single reference image without the sample - making it rather a two-shot acquisition. In both the sample and reference images, a speckle generator, usually a piece of sandpaper or membrane filter, is always present in the beam to generate the speckle pattern. The speckle scanning method requires the movement of the speckle generator or sample in 1-D or 2-D grids, again one scan with the sample in place and one reference scan. In general, the speckle scanning method provides higher spatial resolution and phase sensitivity comparing with the single-image setup. Another application of the SBI methods is to measure the wavefront of the beam in the so-called absolute wavefront sensing mode, which still requires two data sets acquired at two speckle-to-detector positions along the optical axis [5,13,29].

The central data analysis process of SBI methods is to find local spatial shifts of the speckle pattern, and thus the phase angles caused by the sample. Most of existing methods are based on cross-correlation algorithms which are accurate but computationally expensive. For a typical $1024\times 1024$ pixels image with a calculation window of $10\times 10$ pixels, over 100 million cross-correlation calculations are necessary to give accurate results.

Recently, methods based on unified modulated pattern analysis (UMPA) and geometric flow were proposed for speckle-based imaging [25,27,28]. The UMPA method can work with less scanning positions and achieve high sensitivity, but the spatial resolution is limited due to the correlation operation using the template window. Instead of using correlation, the geometric flow method, which is based on the transport of intensity equation (TIE), has advantages in calculation speed and single-shot measurement ability. However, it assumes a weak absorption and scattering sample and the results can be highly affected by the intensity fluctuation. In this work, we focus the discussion and comparison with the correlation-based methods, which has a wider application range.

To improve the efficiency of the correlation-based SBI methods, similarity-comparison based methods can be used for tracking speckle shifts instead of cross-correlation algorithms. Methods based on singular value decomposition (SVD), discrete Fourier transform (DFT), and discrete wavelet transform (DWT) have been already applied in image pattern recognition and time-series similarity comparison [30–34]. The latter two algorithms are both high-speed transforms. The DWT method is the most efficient with a time scale proportional to $n\log (n)$ comparing to that of DFT ($n^2$) for a data length of $n$. In addition, the wavelet transform can provide both spatial frequency and location information, while the Fourier transform gives frequency information only. This feature enables the DWT to provide a trade-off between local information and spatial resolution, which is preferred for speckle tracking in terms of noise robustness.

In this work, a wavelet-transform-based speckle vector tracking (WSVT) method is proposed to improve the SBI efficiency and broaden its application in phase sensing measurements. In Sec. 2. we will present mathematical details of the correlation-based speckle vector tracking (CSVT) and the WSVT algorithms. Performance evaluations and comparisons of the two methods are shown with experimental results in Sec. 3.

## 2. X-ray speckle-based imaging method

Among the many variations of SBI techniques [29], we take the X-ray speckle vector tracking (XSVT) method as an example but note that the discussion can also be extended to single-shot or other scanning methods. As shown in Fig. 1, two data sets, $S^s (x, y, t)$ and $S^r (x, y, t)$, of images are recorded at $n$ different transverse positions of the speckle generator (indexed as $t \in [1, n]$). Here ($x$, $y$) represents the detector coordinates. The sample data set, $S^s (x, y, t)$, is taken with the sample in the beam at a certain sample-to-detector distance, while $S^r (x, y, t)$ is the reference data set taken without the sample at the same detector position. Alternatively, in the absolute wavefront sensing mode, $S^r (x, y, t)$ and $S^s (x, y, t)$ are two data sets taken at different positions along the optical axis [5,13,29].

The scan positions of the speckle generator are the same for both sample and reference data sets for a given $t$. In order to obtain the speckle displacements, 2D scan, either with randomly chosen positions or following a linear or exponential 2D position grid, is used to obtain better accuracy and improve noise robustness compared with 1D scan. The maximum lateral displacement of the 2D scan pattern has to be larger than the largest possible shifts of the speckle pattern in each transverse direction.

The core procedure of the XSVT data analysis is to find the relative displacement (or shift) of each pixel (or a group of pixels within a local window) between the two data sets. This is performed by comparing the 1D intensity scan profiles $I^r_{x,y}(t)$ and $I^s_{x,y}(t)$ of each pixel at $(x,y)$ of the detector extracted from both the reference and sample data sets, respectively (see Fig. 1). The algorithm to carry out this step is so essential to the XSVT method that it not only changes the data analysis efficiency (see Sec. 3.2) but also has impact on the data acquisition of the measurements (see Sec. 3.4).

The sample speckle image can be viewed as a distorted pattern from the reference image. A pixel on the sample image at $(x, y)$ is moved from a position $(x',y')$ on the reference image, given by

where $\Delta x_0$ and $\Delta y_0$ are the displacement of the pixel. $\Delta x_0$ and $\Delta y_0$ are related to the gradient of the wavefront phase $\Phi$ in the horizontal $x$ and vertical $y$ directions, respectively, or with the X-ray wavelength $\lambda$, and $D$ is the propagation distance depending on the geometry. The parameter $a$ is a geometric scaling factor [29]. When the sample is placed downstream of the speckle generator, $D$ is the sample-to-detector distance and $a$ = 1. When the sample is upstream of the speckle generator, $D$ is the speckle-generator-to-detector distance and $a = L/(L + z)$, where $z$ is the sample-to-speckle-generator distance and $L$ is the source-to-sample distance.Once the differential phases $\Phi '_x$ and $\Phi '_y$ are extracted from Eq. (2), the 2D phase profile $\Phi (x, y)$ can be obtained through integration using the Frankot-Chellappa method [35] or any other 2D numerical integration routines [36,37].

#### 2.1 Correlation-based speckle vector tracking (CSVT)

The CSVT method uses the cross-correlation algorithm to determine the similarities of single-pixel signals $I_{x,y}(t)$ from the sample and reference data sets. Displacements $\Delta x_0$ and $\Delta y_0$ of a single pixel at $(x, y)$ on the detector are obtained by maximizing the cross-correlation coefficient between the intensity profiles $I^s_{x,y}(t)$ from the sample data and $I^r_{x', y'}(t)$ from the reference data, or

Traditionally, the process described in Eq. (3) is carried out by comparing each single-pixel profile from the sample data with profiles within a searching window around the same pixel from the reference data. The searching window size, given by $N^w_x \times N^w_y$, needs to be large enough to cover the largest possible pixel displacement, which is sample dependent. For an image area of $N_x \times N_y$ pixels, a total number of $N^w_x \times N^w_y \times N_x \times N_y$ 1D cross-correlation operations are required, which is extremely computational expensive.

It is possible to improve the data-analysis efficiency and resolution using a modified CSVT algorithm. Instead of performing 1D cross-correlations of individual pixel profiles in the searching window, one can use matrix operation to obtain the cross-correlation coefficients for the whole window in one step as

The pixel displacement can be obtained by finding the maximum of the coefficient matrix $C(p,q)$. All results with CSVT shown in this work used the vectorized method. Moreover, sub-pixel resolution is achievable with different minimization and interpolation routines. One of such methods will be given in Sec. 2.3. One should note that the pixel or sub-pixel resolution here refers to the resolution of finding speckle displacements, which determines the angular sensitivity of the method. On the other hand, the spatial resolution of the reconstructed phase image is limited by the effective pixel resolution of the detector system.

#### 2.2 Wavelet-transform-based speckle vector tracking (WSVT)

We introduce here the WSVT method not only to improve the efficiency of data analysis and measurements but also to enhance the noise robustness. A brief introduction of wavelet transform and its comparison with Fourier transform is given first to highlight its fundamental advantage over the XSVT method.

The Fourier transform is the most powerful tool to decompose a signal (spatial signal in this work) into frequency components. Each frequency component covers the entire spatial range of the signal. However, it does not tell where a particular frequency component appears. For XSVT methods, the local information (either along the transverse axes $x$ and $y$ or the scanning axis $t$) of frequency components can also be helpful in the speckle tracking process, which is what the wavelet transform is designed for.

The wavelet transform [38,39] is based on projecting a signal to a set of orthonormal basis functions, named wavelets. By applying the discrete wavelet transform (DWT), the 1D single-pixel signal $I(t)$ can be decomposed into

There are several choices of orthonormal basis functions, such as Haar, Daubechies, Symlets, and so on, available for different applications [39]. The mathematical differences of these bases are beyond the scope of this work. Instead, we only focus on the physical and computational aspects of the DWT. Let’s take an example where $I(t)$ has $n = 2^7 = 128$ points [solid curve in Fig. 2(a)], and the Haar basis function is used for the illustration purpose only. Note that the number of points does not need to be a power of 2. The DWT decomposition is an iterative process. The first level DWT decomposition generates $W_{\phi }(j_1)$ and $W_{\psi }(j_1)$ which each contains $n/2 = 64$ points. The approximation coefficient $W_{\phi }(j_1)$ contains the low-frequency information while the coefficient $W_{\psi }(j_1)$ describes details of the signal. At the same time, both coefficients still retain the local information in contrary to the Fourier transform. If one takes only the $W_{\phi }(j_1)$ part and inverse transform to the $t$ space, the low-pass filtered signal [dotted curve in Fig. 2(a)] can be reconstructed. The first level coefficient $W_{\phi }(j_1)$ can be further decomposed to $W_{\phi }(j_2)$ and $W_{\psi }(j_2)$ with $n/4 = 32$ points. The process is carried out iteratively until the highest order $m$ ($m = 6$ in this case) when $W_{\phi }(j_m)$ has only two points. The inverse wavelet transform is thus a reversed iterative process. Computationally, wavelet coefficients $W_{\phi }$ and $W_{\psi }$ are stored together in an array as shown in Fig. 2(b) indicating the number of points and position of each component. Note that only $W_{\phi }(j_m)$ is stored, but $W_{\psi }(j)$ from all levels are recorded.

The levels of DWT decomposition and number of points in each coefficient also depend on the choice of basis functions and the padding of the data. In this work, we use the python function *wavedec* from the *pywavelets* package with the Daubechies basis (db2) and zero padding mode. For a 1D signal with 100 points and a fifth level decomposition, the numbers of points are 6, 6, 9, 15, 27, and 51 for coefficients $W_{\phi }(j_5)$, and $W_{\psi }(j_5)$ to $W_{\psi }(j_1)$, respectively.

Similar to Fourier transform, one can cut off the low-order detail coefficients from the end of the coefficient array [see Fig. 2(b)] and reconstruct the signal with different levels of details as shown in Fig. 2(a). From this point of view, the DWT process is only slightly different from the Fourier transform. However, since DWT is an orthogonal transform, one can use the Euclidean distance, which is preserved before and after the wavelet transform [33], to compute the similarities between signals. Therefore, the displacement map can be solved by minimizing the square of Euclidean distance for each pixel as

The WSVT method without cutting off any detail coefficients will give identical results as CSVT. The main advantage of WSVT is that by applying detail coefficient cutoff, the noise reduction and computation efficiency can be both improved (see Sec. 3 for details). In addition, the WSVT method can also be applied to single-shot X-ray speckle tracking where a template window was used to find the correlation map, which will also significantly accelerate the reconstruction process.

#### 2.3 Sub-pixel displacement searching

As shown in Eqs. (3) and (6) for CSVT and WSVT, respectively, the displacements $\Delta x_0$ and $\Delta y_0$ are discrete in pixel numbers which limits the displacement resolution to be the single pixel size. In order to achieve sub-pixel resolution, interpolation or gradient estimation method can be used. Here the gradient estimation method is used as an example.

Assuming $R(u,v)$ is the evaluation function (cross-correlation coefficients for CSVT or square of Euclidean distances for WSVT) obtained for a window of pixels and $(u_0, v_0)$ is the pixel position of the extremum (maximum for CSVT and minimum for WSVT). The sub-pixel extremum value of $R(u,v)$ can be expanded in Taylor series as

Using the condition that the derivatives of $R(u+\delta u,v+\delta v)$ of variables $\delta u$ and $\delta v$ are both zero at the extremum position, the sub-pixel displacements can be obtained as

## 3. Experiments and results

XSVT measurements were carried out at the X-ray optics testing beamline 1-BM [40] at the Advanced Photon Source (APS) using a setup shown in Fig. 3. The bending magnet radiation was monochromatized by a Si(111) double-crystal monochromator to X-ray energy of 14 keV. Beryllium (Be) refractive lenses were used as test samples located at 34 m from the source. A membrane filter (speckle generator) with an average pore size of 5 µm was placed at 272 mm downstream of the sample. Thus the $a$ factor in Eq. (2) is 0.992. The distance from the speckle generator to the detector system was $D = 0.5$ m. The detector system (see Fig. 3) contains a 100 µm thick LuAG:Ce scintillator (SC), a 45 degree reflection mirror (M), a $10\times$ objective lens (OL), a tube lens (TL), and an Andor Neo sCMOS camera with $2560 \times 2180$ pixels of 6.5 µm pixel size. The effective pixel size of all recorded images is 0.65 µm.

First, a 2D paraboloid Be lens with 200 µm apex radius of curvature was used as sample, which has a theoretical focal length of $57.5$ m at 14 keV. The sample and reference data sets were both taken with the speckle generator scanned in a $10 \times 10$ mesh ($n = 100$) with a step size of 3 µm in both horizontal and vertical directions. To avoid information missing at a particular frequency by using a regular grid, a random shift up to 0.3 µm was added to each scan step, identical in sample and reference scan.

All data analysis was performed on a laptop with eight cores of i7 CPU and 16 GB RAM. The code was written in python but with the matrix operation optimized in C++. This C++ optimized CSVT code already has improved efficiency comparing with publicly available software [41]. To provide a fair comparison, the CSVT and WSVT codes are only different in the central algorithm part. The code used in this work is available in the GitHub repository [42].

#### 3.1 WSVT benchmarking

Figure 4 shows the reconstructed phase results from both CSVT and WSVT methods following procedures described in section 2. A searching window size of $N^w_x = N^w_y = 40$ was used for each pixel in both analysis. The WSVT method was carried out with the Daubechies (db2) wavelet transform [39] and detail coefficients cut off at the first level (C2), where $W_{\psi }(j_1)$ and $W_{\psi }(j_2)$ were removed.

The reconstructed phase profiles of the single lens are shown in Figs. 4(a) and 4(b) using CSVT and WSVT with the residual phase error profiles (removing the ideal parabolic phase) given in Figs. 4(d) and 4(e), respectively. The agreement between results from the two methods is demonstrated in Fig. 4(c), where the relative phase difference, $(\Phi _{\textrm{WSVT}} - \Phi _{\textrm{CSVT}})/|\Phi _{\textrm{CSVT}}|_{max}$, between the two methods is below $0.3\%$. To clearly show the phase difference, horizontal and vertical line cuts of Figs. 4(d) and 4(e) are compared in Fig. 4(f). The RMS phase difference between CSVT and WSVT is $0.005 \lambda$ over the entire area in Figs. 4(d) and 4(e). The phase sensitivity of the speckle tracking setup in this work can be experimentally determined by using two sets of reference images taken one after the other, which will contain systematic errors from both beam stability and motor repeatability. The RMS phase errors obtained with WSVT and CSVT methods are both 0.007$\lambda$. Therefore, the agreement between WSVT and CSVT algorithms is beyond the measurement systematic error.

The spatial resolution of the WSVT and CSVT methods is expected to be the same, since both have the same experimental conditions including the speckle pattern size, scan steps, and the resolution of the detector system. In this work, the spatial resolution was dominated by the detector system, which is estimated to be 2.2 µm considering the thickness effect of the scintillator [43].

#### 3.2 Efficiency comparison

The one major advantage of the WSVT method is the improved calculation speed. Here we compare the computation time of CSVT and WSVT with different cutoff levels of wavelet detail coefficients [C1: removed $W_{\psi }(j_1)$; C2: removed $W_{\psi }(j_1)$ and $W_{\psi }(j_2)$; C3: removed $W_{\psi }(j_1)$, $W_{\psi }(j_2)$ and $W_{\psi }(j_3)$]. Note that the WSVT without cutoff is identical to CSVT but with slightly increase in calculation time due to the additional wavelet transform step, and thus should be avoided.

Table 1 summarizes results for different image area sizes $N_x \times N_y$ (with a fixed searching window size of $40 \times 40$ pixels) and different searching window sizes $N^w_x \times N^w_y$ (with a fixed image size of $512 \times 512$ pixels). The WSVT method is certainly faster than the CSVT analysis because of the reduction in data volume through the wavelet transform and cutoff.

The time cost increases linearly with the detector image size (sample size). The WSVT analysis with C2 and C3 cutoff can be more than a factor of two faster than CSVT at large image sizes (e.g., $N_x=N_y \ge 1024$). The speed advantage of WSVT is more pronounced in cases of large searching window sizes, which are needed when the sample-induced phase shift is large (e.g., strongly focusing optics as shown in Sec. 3.4). The speed of WSVT-C2 and WSVT-C3 are nearly three and four times faster than CSVT, respectively, for $N^w_x=N^w_y = 160$. In general, the WSVT method shows much higher efficiency for conditions with large image sizes and searching window sizes.

WSVT analysis with higher-order cutoffs is sure faster because of the smaller data volume. However, the cutoff level needs to be low enough to ensure accuracy. Figure 5 shows the residual phase results with different levels of wavelet coefficient cutoff. Because of padding and the db2 basis function, WSVT with all coefficients (114 points) has slightly higher data volume than CSVT (100 points), while WSVT-C1, C2, and C3 reduce the data points to 63, 36, and 21, respectively. The phase results in Figs. 5(b) and 5(c) are almost identical to that in (a). Truncating detail coefficients up to the second level does not imperil the result accuracy in this case. However, Fig. 5(d) shows that the C3 cutoff will give too large phase errors. The RMS phase errors between Fig. 5 and the result from CSVT (194 s calculation time) are (a) $10^{-14} \lambda$, (b) 0.007 $\lambda$, (c) 0.01 $\lambda$, and (d) 0.17 $\lambda$, with the calculation time of (a) 216 s, (b) 145 s, (c) 109 s, and (d) 89 s, respectively. Of course, this trade-off between accuracy and speed depends on the experimental condition and data quality. We show next that the WSVT is preferred in dealing with high noise data, where WSVT diplays higher accuracy at faster speed as compared to CSVT.

#### 3.3 Noise robustness comparison

Thanks to the special features of the wavelet transform, the displacement searching process has an intrinsic noise filtering effect, which is another major advantage of the WSVT method. To study the noise robustness of both CSVT and WSVT methods, different levels of Gaussian noise were added to the signal of each pixel mimicking the scattering background noise and the dark detector noise. The noise level was controlled by the Gaussian $\sigma$ value in terms of percentage of the signal.

The phase reconstruction results are compared in Fig. 6. The result of WSVT analysis with no coefficient cutoff and without added noise was taken as the reference, which is also equivalent to the CSVT result without noise. The phase reconstruction error using WSVT-C2 is much less than that of CSVT, especially at high noise levels ($\sigma \ge 0.1$). The RMS reconstruction errors of WSVT-C2 are $0.01 \lambda$ and $0.037 \lambda$ for noise levels of $\sigma = 0.1$ and 0.15, respectively, while the RMS errors of the CSVT method are $0.025 \lambda$ and $0.058 \lambda$. The strong noise robustness of WSVT can expand the application of the XSVT technique into low signal conditions (e.g., using a lab source) as well as the potential use of fast detectors with low dynamic range.

#### 3.4 Scan step reduction

The high-resolution and noise robustness of the speckle vector tracking technique, in general, relies on the redundancy of image collection. The XSVT data is usually acquired with a 2D scan, which can be totally random 2D scan, a uniform linear 2D scan with (as in this work) or without random step fluctuations or a non-uniform 2D scan with logarithmic step sizes. Nevertheless, the total number of scan steps ($n = n_x \times n_y$) must be large enough so that each pixel vector can be distinguishable from its neighboring pixels. Increasing $n$ will improve accuracy but at the cost of increasing measurement and data analysis time. Slower measurements will also be more vulnerable to the experimental instabilities such as the vibration and drift of optics. In this section, we study the possibility of reducing the number of scan steps to speed up experiments without compromising data quality and data processing speed using the WSVT and CSVT methods.

With the original $n = 100$ ($10 \times 10$), Fig. 5 in Sec. 3.2 shows that the WSVT-C2 analysis can still provide accurate results, which means the number of data points can be reduced if no noise is present in the data. To valid this prediction, WSVT-C1 analysis was applied to different $n$ steps ranging from $3 \times 3$ to $10\times {10}$ and compared to CSVT analysis. The relative reconstruction errors, taking the $10\times {10}$ case as the reference, are plotted as a function of $n$, as shown in Fig. 7. Both WSVT-C1 and CSVT methods show similar errors with the same number of scan steps. A scan number of $n = 5^2 = 25$ is adequate to provide accurate results in this particular case, which is consistent with our prediction above. Furthermore, one should note that even with such a small number of scan steps, WSVT-C1 managed to reduce the data by an additional factor of two through cutting off $W_{\psi }(j_1)$, thanks to the robust noise reduction feature of the WSVT method.

#### 3.5 Strongly focusing optics

Strongly focusing optics are the most challenging samples to measure with any at-wavelength metrology tools. The large phase angle created by the optics generates substantial speckle pattern displacement and distortion on the detector [5]. As a result, the WSVT analysis requires large searching window sizes, which means a dramatic increase in calculation time. On the other hand, to deliver micro- to nano-focusing beam with high-quality wavefront, phase angle errors caused by these optics need to be in the microradian or even nanoradian level, which is the desired angular sensitivity for the at-wavelength metrology techniques. In this section, we show that the WSVT method can handle these strongly focusing optics with high efficiency and accuracy.

Because of the small decrement of X-ray refraction index, a combination of multiple lenses, namely the compound refractive lens (CRL), is necessary to provide strong focusing. In this work, a CRL containing 18 Be lenses with 200 µm apex radius of curvature was used to give a theoretical focal length of 3.2 m at 14 keV. The speckle-generator-to-detector distance was reduced to $D =$ 310 mm in order to limit the searching window size at the expense of a lower angular sensitivity.

The WSVT analysis was carried out with $1700 \times 1700$ pixels corresponding to an image size of $1.1 \times 1.1~\rm {mm}^2$ with a searching window size of $200 \times 200$ pixels. The coefficient cutoff was selected to the C2 level to provide high speed and adequate precision. The phase reconstruction results are shown in Fig. 8, illustrating the differential phase, integrated phase, and the residual phase after removing the parabolic fit.

The reconstructed differential phases $\Phi ^{'}_x(x,y)$ and $\Phi ^{'}_y(x,y)$ for the horizontal and vertical directions are shown in Figs. 8(a) and 8(b), respectively. The integrated phase profile $\Phi (x,y)$ in Fig. 8(c) shows a large total phase which is approximately a factor of 18 larger than the single CRL phase shown in Figs. 4(a) and 4(b). After removing the ideal parabolic phase, the residual phase profile is shown in Fig. 8(d) with a peak-to-valley phase error of $6 \lambda$. The Zernike polynomial decomposition of the residual phase profile shown in Fig. 8(d) suggests the dominating trefoil aberrations (i.e., $Z^{-3}_3$ and $Z^3_3$) and astigmatism (i.e., $Z^2_2$).

Using the same laptop computer described earlier, the total calculation time of the WSVT-C2 analysis was $40$ min, which is 2.6 times faster than the optimized CSVT code. The success in such a strongly focusing case further proves that the WSVT method can be suggested as a general XSVT treatment and a possible replacement of CSVT.

## 4. Conclusion

A new wavelet-transform-based speckle vector tracking method is introduced, analyzed, and compared to the cross-correlation-based method. The results show that the WSVT method has the unique feature of controllable noise reduction through different levels of coefficient cutoff and with improved efficiency at the same time. The WSVT method is faster than the CSVT method, especially when dealing with large image sizes and phase shifts. The noise robustness of WSVT can allow measurements under low flux conditions, faster scans, and fewer scan steps without sacrificing the reconstruction quality. The WSVT method has a high angular dynamic range, which makes it feasible for probing strong focusing optics with high sensitivity. In addition, the absorption and dark-field signal of the sample can also be obtained with the same process as other speckle-based imaging methods [22]. In summary, the WSVT method can provide high angular sensitivity, strong noise-robustness, and fast processing capability, which can strength the application of the speckle-based methods to wave-sensing related fields such as at-wavelength metrology [41] and optical design with specification of optical components [44]; calculation of corrective optics [45] or beam shapers in general; beamline alignment [46] and general wavefront diagnostics [47,48]; and phase-contrast imaging applications [49].

## Funding

U.S. Department of Energy (No. DE-AC02-06CH11357).

## Acknowledgments

The authors would like to thank Dr. Michael Wojcik (Argonne National Laboratory) for the experimental support at the APS 1-BM beamline and Dr. Thomas Roth (European Synchrotron Radiation Facility) for valuable discussions.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **D. H. Bilderback, P. Elleaume, and E. Weckert, “Review of third and next generation synchrotron light sources,” J. Phys. B: At., Mol. Opt. Phys. **38**(9), S773–S797 (2005). [CrossRef]

**2. **E. Allaria, D. Castronovo, P. Cinquegrana, P. Craievich, M. Dal Forno, M. B. Danailov, G. D’Auria, A. Demidovich, G. De Ninno, S. Di Mitri, B. Diviacco, W. M. Fawley, M. Ferianis, E. Ferrari, L. Froehlich, G. Gaio, D. Gauthier, L. Giannessi, R. Ivanov, B. Mahieu, N. Mahne, I. Nikolov, F. Parmigiani, G. Penco, L. Raimondi, C. Scafuri, C. Serpico, P. Sigalotti, S. Spampinati, C. Spezzani, M. Svandrlik, C. Svetina, M. Trovo, M. Veronese, D. Zangrando, and M. Zangrando, “Two-stage seeded soft-X-ray free-electron laser,” Nat. Photonics **7**(11), 913–918 (2013). [CrossRef]

**3. **A. Stevenson, T. Gureyev, D. Paganin, S. Wilkins, T. Weitkamp, A. Snigirev, C. Rau, I. Snigireva, H. Youn, I. Dolbnya, W. Yun, B. Lai, R. Garrett, D. Cookson, K. Hyodo, and M. Ando, “Phase-contrast x-ray imaging with synchrotron radiation for materials science applications,” Nucl. Instrum. Methods Phys. Res., Sect. B **199**, 427–435 (2003). [CrossRef]

**4. **R. A. Lewis, “Medical phase contrast x-ray imaging: current status and future prospects,” Phys. Med. Biol. **49**(16), 3573–3583 (2004). [CrossRef]

**5. **S. Berujon, E. Ziegler, and P. Cloetens, “X-ray pulse wavefront metrology using speckle tracking,” J. Synchrotron Radiat. **22**(4), 886–894 (2015). [CrossRef]

**6. **F. Pfeiffer, “X-ray ptychography,” Nat. Photonics **12**(1), 9–17 (2018). [CrossRef]

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

**8. **C. Kottler, C. David, F. Pfeiffer, and O. Bunk, “A two-directional approach for grating based differential phase contrast imaging using hard x-rays,” Opt. Express **15**(3), 1175–1181 (2007). [CrossRef]

**9. **H. H. Wen, E. E. Bennett, R. Kopace, A. F. Stein, and V. Pai, “Single-shot x-ray differential phase-contrast and diffraction imaging using two-dimensional transmission gratings,” Opt. Lett. **35**(12), 1932–1934 (2010). [CrossRef]

**10. **P. Modregger, F. Scattarella, B. Pinzer, C. David, R. Bellotti, and M. Stampanoni, “Imaging the Ultrasmall-Angle X-Ray Scattering Distribution with Grating Interferometry,” Phys. Rev. Lett. **108**(4), 048101 (2012). [CrossRef]

**11. **P. Modregger, S. Rutishauser, J. Meiser, C. David, and M. Stampanoni, “Two-dimensional ultra-small angle X-ray scattering with grating interferometry,” Appl. Phys. Lett. **105**(2), 024102 (2014). [CrossRef]

**12. **S. Berujon, H. Wang, and K. Sawhney, “X-ray multimodal imaging using a random-phase object,” Phys. Rev. A **86**(6), 063813 (2012). [CrossRef]

**13. **S. Bérujon, E. Ziegler, R. Cerbino, and L. Peverini, “Two-dimensional x-ray beam phase sensing,” Phys. Rev. Lett. **108**(15), 158102 (2012). [CrossRef]

**14. **K. S. Morgan, D. M. Paganin, and K. K. Siu, “X-ray phase imaging with a paper analyzer,” Appl. Phys. Lett. **100**(12), 124102 (2012). [CrossRef]

**15. **L. Assoufid, X. Shi, S. Marathe, E. Benda, M. J. Wojcik, K. Lang, R. Xu, W. Liu, A. T. Macrander, and J. Z. Tischler, “Development and implementation of a portable grating interferometer system as a standard tool for testing optics at the Advanced Photon Source beamline 1-BM,” Rev. Sci. Instrum. **87**(5), 052004 (2016). [CrossRef]

**16. **S. Marathe, L. Assoufid, X. Xiao, K. Ham, W. W. Johnson, and L. G. Butler, “Improved algorithm for processing grating-based phase contrast interferometry image sets,” Rev. Sci. Instrum. **85**(1), 013704 (2014). [CrossRef]

**17. **T. Inoue, S. Matsuyama, S. Kawai, H. Yumoto, Y. Inubushi, T. Osaka, I. Inoue, T. Koyama, K. Tono, H. Ohashi, M. Yabashi, T. Ishikawa, and K. Yamauchi, “Systematic-error-free wavefront measurement using an X-ray single-grating interferometer,” Rev. Sci. Instrum. **89**(4), 043106 (2018). [CrossRef]

**18. **F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance x-ray sources,” Nat. Phys. **2**(4), 258–261 (2006). [CrossRef]

**19. **Y. Liu, M. Seaberg, D. Zhu, J. Krzywinski, F. Seiboth, C. Hardin, D. Cocco, A. Aquila, B. Nagler, H. J. Lee, S. Boutet, Y. Feng, Y. Ding, G. Marcus, and A. Sakdinawat, “High-accuracy wavefront sensing for x-ray free electron lasers,” Optica **5**(8), 967 (2018). [CrossRef]

**20. **I. Zanette, T. Zhou, A. Burvall, U. Lundström, D. H. Larsson, M. Zdora, P. Thibault, F. Pfeiffer, and H. M. Hertz, “Speckle-Based X-Ray Phase-Contrast and Dark-Field Imaging with a Laboratory Source,” Phys. Rev. Lett. **112**(25), 253903 (2014). [CrossRef]

**21. **T. Zhou, I. Zanette, M.-C. Zdora, U. Lundstrom, D. H. Larsson, H. M. Hertz, F. Pfeiffer, and A. Burvall, “Speckle-based x-ray phase-contrast imaging with a laboratory source and the scanning technique,” Opt. Lett. **40**(12), 2822 (2015). [CrossRef]

**22. **H. Wang, Y. Kashyap, and K. Sawhney, “From synchrotron radiation to lab source: advanced speckle-based X-ray imaging using abrasive paper,” Sci. Rep. **6**(1), 20476 (2016). [CrossRef]

**23. **M. Seaberg, R. Cojocaru, S. Berujon, E. Ziegler, A. Jaggi, J. Krempasky, F. Seiboth, A. Aquila, Y. Liu, A. Sakdinawat, H. Lee, U. Flechsig, L. Patthey, F. Koch, G. Seniutinas, C. David, D. Zhu, L. Mikeš, M. Makita, T. Koyama, A. P. Mancuso, H. N. Chapman, and P. Vagovič, “Wavefront sensing at X-ray free-electron lasers,” J. Synchrotron Radiat. **26**(4), 1115–1126 (2019). [CrossRef]

**24. **S. Berujon and E. Ziegler, “Near-field speckle-scanning-based x-ray imaging,” Phys. Rev. A **92**(1), 013837 (2015). [CrossRef]

**25. **M.-C. Zdora, P. Thibault, T. Zhou, F. J. Koch, J. Romell, S. Sala, A. Last, C. Rau, and I. Zanette, “X-ray Phase-Contrast Imaging and Metrology through Unified Modulated Pattern Analysis,” Phys. Rev. Lett. **118**(20), 203903 (2017). [CrossRef]

**26. **S. Berujon and E. Ziegler, “X-ray multimodal tomography using speckle-vector tracking,” Phys. Rev. Appl. **5**(4), 044014 (2016). [CrossRef]

**27. **D. M. Paganin, H. Labriet, E. Brun, and S. Berujon, “Single-image geometric-flow x-ray speckle tracking,” Phys. Rev. A **98**(5), 053813 (2018). [CrossRef]

**28. **K. M. Pavlov, H. T. Li, D. M. Paganin, S. Berujon, H. Rougé-Labriet, and E. Brun, “Single-shot x-ray speckle-based imaging of a single-material object,” Phys. Rev. Appl. **13**(5), 054023 (2020). [CrossRef]

**29. **S. Berujon, R. Cojocaru, P. Piault, R. Celestre, T. Roth, R. Barrett, and E. Ziegler, “X-ray optics and beam characterization using random modulation: theory,” J. Synchrotron Radiat. **27**(2), 284–292 (2020). [CrossRef]

**30. **C. Faloutsos, M. Ranganathan, and Y. Manolopoulos, “Fast subsequence matching in time-series databases,” ACM SIGMOD Rec. **23**(2), 419–429 (1994). [CrossRef]

**31. **D. Wu, A. Singh, D. Agrawal, A. El Abbadi, and T. R. Smith, “Efficient retrieval for browsing large image databases,” in Proceedings of the fifth international conference on Information and knowledge management, (1996), pp. 11–18.

**32. **K.-P. Chan and A. W.-C. Fu, “Efficient time series matching by wavelets,” Proceedings 15th International Conference on Data Engineering (Cat. No.99CB36337) pp. 126–133 (1999).

**33. **K. Fukunaga, “Introduction to statistical pattern recognition,” Comput. Sci. Sci. Comput. (1990).

**34. **I. Popivanov and R. J. Miller, “Similarity Search Over Time-Series Data using Wavelets,” Proceedings 18th International Conference on Data Engineering pp. 212–221 (2002).

**35. **R. Frankot and R. Chellappa, “A method for enforcing integrability in shape from shading algorithms,” IEEE Trans. Pattern Anal. Machine Intell **10**(4), 439–451 (1988). [CrossRef]

**36. **A. Agrawal, R. Raskar, and R. Chellappa, “What is the range of surface reconstructions from a gradient field?” in * Computer Vision–ECCV 2006*, (Springer, 2006), pp. 578–591.

**37. **M. Harker and P. O’Leary, “Matlab toolbox for the regularized surface reconstruction from gradients,” Proc. SPIE **9534**, 95341E (2015). [CrossRef]

**38. **C. E. Heil and D. F. Walnut, “Continuous and discrete wavelet transforms,” SIAM Rev. **31**(4), 628–666 (1989). [CrossRef]

**39. **M. J. Shensa, “The discrete wavelet transform: wedding the a trous and mallat algorithms,” IEEE Trans. Signal Process. **40**(10), 2464–2482 (1992). [CrossRef]

**40. **A. Macrander, M. Erdmann, N. Kujala, S. Stoupin, S. Marathe, X. Shi, M. Wojcik, D. Nocher, R. Conley, J. Sullivan, K. Goetze, J. Maser, and L. Assoufid, “X-ray optics testing beamline 1-BM at the advanced photon source,” AIP Conf. Proc. **1741**, 030030 (2016). [CrossRef]

**41. **S. Berujon, R. Cojocaru, P. Piault, R. Celestre, T. Roth, R. Barrett, and E. Ziegler, “X-ray optics and beam characterization using random modulation: experiments,” J. Synchrotron Radiat. **27**(2), 293–304 (2020). [CrossRef]

**42. **Z. Qiao and X. Shi, “Wavelet-transform based speckle vector tracking (WSVT),” https://github.com/APS-XSD-OPT-Group/WaveletSpeckleTracking.

**43. **A. Koch, C. Raven, P. Spanne, and A. Snigirev, “X-ray imaging with submicrometer resolution employing transparent luminescent screens,” J. Opt. Soc. Am. A **15**(7), 1940–1951 (1998). [CrossRef]

**44. **R. Celestre, S. Berujon, T. Roth, M. Sanchez del Rio, and R. Barrett, “Modelling phase imperfections in compound refractive lenses,” J. Synchrotron Radiat. **27**(2), 305–318 (2020). [CrossRef]

**45. **F. Seiboth, A. Schropp, M. Scholz, F. Wittwer, C. Rödel, M. Wünsche, T. Ullsperger, S. Nolte, J. Rahomäki, K. Parfeniukas, S. Giakoumidis, U. Vogt, U. Wagner, C. Rau, U. Boesenberg, J. Garrevoet, G. Falkenberg, E. C. Galtier, H. Ja Lee, B. Nagler, and C. G. Schroer, “Perfect X-ray focusing via fitting corrective glasses to aberrated optics,” Nat. Commun. **8**(1), 14623 (2017). [CrossRef]

**46. **T. Zhou, H. Wang, O. J. L. Fox, and K. J. S. Sawhney, “Optimized alignment of X-ray mirrors with an automated speckle-based metrology tool,” Rev. Sci. Instrum. **90**(2), 021706 (2019). [CrossRef]

**47. **S. Berujon, E. Ziegler, and P. Cloetens, “X-ray pulse wavefront metrology using speckle tracking,” J. Synchrotron Radiat. **22**(4), 886–894 (2015). [CrossRef]

**48. **M. Seaberg, R. Cojocaru, S. Berujon, E. Ziegler, A. Jaggi, J. Krempasky, F. Seiboth, A. Aquila, Y. Liu, A. Sakdinawat, H. J. Lee, U. Flechsig, L. Patthey, F. Koch, G. Seniutinas, C. David, D. Zhu, L. Mikeš, M. Makita, T. Koyama, A. P. Mancuso, H. N. Chapman, and P. Vagovič, “Wavefront sensing at X-ray free-electron lasers,” J. Synchrotron Radiat. **26**(4), 1115–1126 (2019). [CrossRef]

**49. **S. Berujon and E. Ziegler, “Near-field speckle-scanning-based x-ray tomography,” Phys. Rev. A **95**(6), 063822 (2017). [CrossRef]