Denoising X-ray Diffraction two-dimensional patterns by a Lattice Boltzmann Method

1. Introduction

X-ray tomography is a versatile technique for investigating the internal structure of nano(bio)-materials with applications in energy research [1], life sciences [2], medicine [3,4,5] and many other fields, both at synchrotron radiation facilities and at microfocus laboratory sources [6].

In recent years, advancements in X-ray diffraction (XRD) image denoising have continued to evolve, driven by both traditional [7]1 and machine learning approaches, particularly (convolutional) neural networks and techniques alike [8,9,10].

In spite of the prominence recently gained by deep learning techniques, it is worth stressing that these models have to be trained on large datasets to learn to distinguish between noise and signal, and improving the clarity of XRD images. For instance, new benchmark datasets and evaluation metrics are being developed to standardize the assessment of denoising techniques and, to enhance model performance, researchers use data augmentation techniques: this involves artificially increasing the size of the training dataset through transformations like rotations, scaling, and noise addition, which helps in training more robust models [11]. On the other hand, the integration of advanced machine learning techniques with traditional denoising methods is significantly enhancing the quality of X-ray diffraction images, leading to more accurate analyses and interpretations in various scientific and industrial applications. This is particularly the case of real-time processing, where advances in hardware and software are making it feasible to perform denoising in real-time, which is particularly useful for high-throughput XRD applications. Also combining machine learning with traditional signal processing techniques is becoming more common (hybrid approaches). For example, pre-processing steps might use classical denoising methods, followed by machine learning techniques for further refinement.

Furthermore tackling the problem of image restoration solely from an algorithmic perspective may even include generative aspects, that can bias the ground truth: noise profiles are unknown and often include noise from multiple distinct sources, substantially reducing the applicability of simulation-based approaches. Therefore, techniques like wavelet transform and Physics inspired adaptive filtering are still being refined to handle varying types of noise and different XRD data characteristics. In order to accomplish the task of preserving important features while reducing noise, focusing on a genuinely Physics-grounded method is still worth.

2. A Fluidistic Approach

An XRD two-dimensional pattern is described as a fluid representing the signal and the noise; the denoising process can be thought of as the result of the motion of such a fluid towards a steady state (the original noisy image): due to internal friction (viscosity) and thermal conduction, the effect of energy dissipation, occurring during the motion of this fluid, on that motion itself, results in the thermodynamic irreversibility of the motion. The different diffusive dynamics of signal/noise finally smears the noise and this process is irreversible. In order to obtain the equations describing the motion of a viscous fluid, some terms must be added to the Euler’s equations modelling an ideal fluid.

The Navier–Stokes equations are fundamental mathematical equations in fluid mechanics that describe the motion of viscous fluid substances. They are derived from the principles of conservation of mass, momentum, and energy. These nonlinear partial differential equations account for various factors, including pressure, velocity, and external forces, making them essential for understanding fluid behavior in diverse applications, from weather patterns to aerospace engineering. Solving the Navier–Stokes equations provides insights into complex fluid dynamics, although they are known for their mathematical challenges and are central to the unsolved problem of turbulence in fluid flow [12].

The diffusion equation and the Navier-Stokes equations are both fundamental in fluid dynamics but describe different physical phenomena.

• Diffusion Equation. This equation models the process of substance spreading due to concentration gradients. It is typically written as ${\displaystyle \frac{\partial c}{\partial t}=D{\nabla }^{2}c}$, where c is the concentration of the substance, D is the diffusion coefficient, and ${\nabla }^2$ is the Laplacian operator. It captures how substances diffuse through a medium over time.
• Navier–Stokes Equations. These equations describe the motion of fluid substances and account for viscosity and external forces. They can be expressed as:
  ${\displaystyle \frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\nabla )\mathbf{u}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{u}+\mathbf{f}}$, where $\mathbf{u}$ is the fluid velocity, p is the pressure, $\nu$ is the kinematic viscosity, $\rho$ is the density, and $\mathbf{f}$ represents external forces.

There are relationships between the two equations, namely viscous effects, transport phenomena and limits. Indeed, the diffusion equation can be seen as a simplified version of the Navier–Stokes equations when considering mass transport in a viscous fluid. In this context, the term representing diffusion can be derived from the viscous effects of the fluid motion. Furthermore, the Navier–Stokes equations govern fluid flow, while the diffusion equation describes how quantities like temperature, concentration, or species disperse within that flow2. Finally in low Reynolds number flows (where inertial forces are negligible compared to viscous forces), the Navier–Stokes equations can approximate diffusion-like behavior.

In summary, while both equations are foundational to understanding fluid behavior, the diffusion equation focuses on concentration changes, whereas the Navier–Stokes equations address the overall dynamics of fluid motion.

2.1. A Navier–Stokes Inspired Model

In the Field Theory, the relevant functions (i.e. Lagrangian $\mathcal{L}$ and Hamiltonian $\mathcal{H}$ densities) to describe the model of a dynamical field (u), having a mass ($\propto w^{-1}$) and interacting with an external field ($\rho$), read:

$$\begin{array}{ccc}\hfill \mathcal{L}& \stackrel{def.}{=}& {\displaystyle \frac{ı}{2}}\left(u^{★}{\partial }_{t}u-{\partial }_t{u}^{★}u\right)-\alpha \left\{{\beta |u-\rho |}^2+\right.\hfill \end{array}$$

$$\begin{array}{ccc}& +& \left.w\nabla u^{★}·\nabla u+{\lambda }^2\left(w-lnw\right)\right\};\hfill \\ \hfill \mathcal{H}& \stackrel{def.}{=}& {\left[\Pi {\partial }_{t}u+{\Pi }^{★}{\partial }_t{u}^{★}-\mathcal{L}\right]}_{{\Pi }^{(★)}\stackrel{def.}{=}{\displaystyle \frac{\partial \mathcal{L}}{\partial \left({\partial }_t{u}^{(★)}\right)}}={\left({\displaystyle \frac{ı}{2}}{u}^{★}\right)}^{(★)}}=\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& =& \alpha \left[{\beta |u-\rho |}^2+w\nabla u^{★}·\nabla u+{\lambda }^2\left(w-lnw\right)\right].\hfill \end{array}$$

(2)

The model time evolution follows from the minimization of the functionals energy ($\mathcal{E}$) and action ($\mathcal{S}$):

$$\begin{array}{ccc}\hfill \mathcal{E}\left[u^{(★)},\nabla u^{(★)},w\right]& =& \int \int d^2\overrightarrow{x}\mathcal{H}\left(u^{(★)},\nabla u^{(★)},w\right),\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \mathcal{S}\left[u^{(★)},\nabla u^{(★)},w\right]& =& \int dt\int \int d^2\overrightarrow{x}\mathcal{L}\left(u^{(★)},\nabla u^{(★)},w\right),\hfill \end{array}$$

(4)

where the functional dependence on the dynamical fields $u^{(★)}$ (together with their conjugate momentum densities, $\nabla u^{(★)}$) and on the non-dynamical3 field w is made explicit.

Finally the equations of motion for the field(s) promptly read:

$$\left.\begin{array}{c}\delta \mathcal{S}=0\hfill \\ \tau =ıt\hfill \end{array}\right\}⟹\left\{\begin{array}{c}{\partial }_{\tau }{u}^{(★)}=-\alpha \left[\beta (\rho -u^{(★)})\right.+\hfill \\ +\left.\begin{array}{|c|}\hline \nabla \left(w\nabla u^{(★)}\right)\\ \hline\end{array}\right],\hfill \\ \\ \begin{array}{|c|}\hline {\displaystyle w={\left(1+\frac{\nabla u^{★}·\nabla u}{{\lambda }^2}\right)}^{-1}}\\ \hline\end{array},\hfill \end{array}\right.$$

(5)

where the boxed expressions represent the density (e.g. XRD image) gradient-dependent diffusive term. While this approach has been fruitfully applied to the subject in the past [7], its feasibility is undermined by the difficulty in separating the different time-scales, in particular the collision time-scale and the propagating one.

2.2. Lattice Boltzmann Method: Overview

The Lattice Boltzmann Methods (LBM) is a numerical approach used to simulate fluid dynamics [13,14,15]. Unlike traditional methods like the Navier–Stokes equations, LBM models fluid flow by simulating the movement of particles on a discrete lattice grid as follows.

• Discrete Lattice. The fluid is represented by particles moving on a fixed grid with discrete directions: the lattice dimensions (D) and the neighborhood (Q, the number of nearest neighbors, including the site itself) define the particular lattice Boltzmann method: e.g. D2Q9 means nine nearest neighbors for each site on a two-dimensional lattice ${\displaystyle \underset{\swarrow \downarrow \searrow }{\stackrel{\nwarrow \uparrow \nearrow }{\leftarrow ·\to }}}$.
• Collision Step. Particles collide and exchange momentum according to probabilistic rules.
• Streaming Step. After collisions, particles move to neighboring lattice sites.
• Macroscopic Variables. From the particle distribution functions, macroscopic quantities like density and velocity are computed.

While LBM is particularly useful for complex boundary conditions and multiphase flows due to its straightforward implementation and ability to handle irregular geometries, this method allows for better handling of noise at different scales by adapting to varying noise levels across the two-dimensional X-ray patterns.

2.3. Lattice Boltzmann Method: Application

Due to the equivalence between the hydrodynamic model (transport equation) and the fluid-mechanics (Hamiltonian), an approach is explored where the two dynamics can be explicitly and straightforwardly disentangled. In the hydrodynamic model, indeed, the (microscopic) distribution ${\displaystyle f=f(t,\overrightarrow{x},\overrightarrow{p})}$ satisfies the Boltzmann transport equation ${\displaystyle \frac{df}{d\tau }=C\left[f\right]}$, where the total time derivative reads ${\displaystyle \frac{d(·)}{dt}=\frac{\partial (·)}{\partial t}+\stackrel{·}{\overrightarrow{x}}·{\nabla }_{\overrightarrow{x}}+\stackrel{·}{\overrightarrow{p}}·{\nabla }_{\overrightarrow{p}}}$, while $C\left[f\right]$ is the collision integral accounting for the particle interactions4. The change to macroscopic equations is achieved by integrating out the degrees of freedom which are integrals of the motion between two collisions (in the absence of an external field), e.g. ${\displaystyle u(\overrightarrow{x},t)=\int \int d^2\overrightarrow{p}f(\overrightarrow{x},\overrightarrow{p},t)}$.

According to the Hamilton’s equations [18] ${\displaystyle \stackrel{·}{\overrightarrow{p}}={\overrightarrow{F}}_{ext}}$, being the latter one the external force acting on the fluid particle (the two-dimensional X-ray pattern image pixel).

$$\begin{array}{ccc}\hfill C\left[f\right]& \approx & -\omega \left(f-f_{eq}\right),\hfill \\ \hfill f_{eq}(\overrightarrow{x},\overrightarrow{p})& \stackrel{def.}{=}& {\displaystyle \frac{\rho \left(\overrightarrow{x}\right)}{2\pi }}exp\left(-{\displaystyle \frac{|\overrightarrow{p}{|}^2}{2}}\right),\hfill \end{array}$$

(6)

where $f_{eq}$ satisfies the transport equation with ${\displaystyle {\overrightarrow{F}}_{ext}=-\nabla U}$ and the external potential ${\displaystyle U=-ln\left(\rho \right)}$. The collision integral has been approximated by the Bhatnagar–Gross–Krook model [19], where the linear collision time-scale (${\omega }^{-1}$) has been spotted as the relaxation of the time dependent X-ray two-dimensional pattern ${\displaystyle f(t,\overrightarrow{x},\overrightarrow{p})}$ to the constant noisy one ${\displaystyle f_{eq}(\overrightarrow{x},\overrightarrow{p})\propto \rho \left(\overrightarrow{x}\right)}$. Such an evolution accounts for the different diffusion dynamics of noise/signal to improve the signal-to-noise ratio. The lattice (Boltzmann) method simulates the evolution through a cycle of propagation and collision on a regular lattice where the microscopic distribution function propagates to the nearest neighbor sites in precisely one time-step, once a suitable discrete momenta set has been chosen. By performing such a spatio-temporal discretization, the Boltzmann transport equation reads

$$f_i(\overrightarrow{x}+{\overrightarrow{p}}_i\delta t,t+\delta t)-f_i(\overrightarrow{x},t)=-{\tilde{\omega }}^{-1}\left[f_i(\overrightarrow{x},t)-f_{eqi}\left(\overrightarrow{x}\right)\right],$$

(7)

where ${\tilde{\omega }}^{-1}$ is the re-scaled relaxation time and the index i refers to the set of momenta $\left\{{\overrightarrow{p}}_i\right\}$ chosen to connect the lattice sites5. The anisotropic diffusive term (the boxed expressions in Equation (5)) is added by hand. The algorithm straightforwardly follows: updating a few number of floating-point arithmetic operations is computationally not demanding and the task can be easily parallelized6.

3. Results and Conclusions

The only parameters that are left in the numerical solution of Equation (7) are $\alpha$ and $\beta$ (the parameter $\lambda$, related to the non-dynamical field w, is proportional to the norm of the conjugate momentum density $|\nabla u|$ and it adiabatically changes [21]). While $\alpha =1.15$ is the highest value that can be used to guarantee stability without slowing down the convergence when the method is applied to an image affected by a gaussian noise (e.g. the classical Lenna image of Figure 1), $\alpha =0.85$ appears to be the optimal solution for an application of the technique to the XRD two-dimensional patterns (e.g. Figure 2, Figure 3 and Figure 4). In fact, one has to balance between the need to use the highest possible value of ${\tilde{\omega }}^{-1}$ of the Equation (7) to reduce the number of iterations to achieve the final reconstructed image and the stability requirement that fixes an upper bound for ${\tilde{\omega }}^{-1}$. The efficiency in image filtering is evaluated by measuring the peak-to-signal noise ratio (PSNR) defined as

$$\mathrm{PSNR}\stackrel{def.}{=}-20{log}_{10}\left({\displaystyle \frac{\sqrt{\langle |u-{\rho }_0{|}^2\rangle }}{max\left({\rho }_0\right)}}\right),$$

(8)

where ${\rho }_0$ is the original image without noise and $\langle ...\rangle$ denotes a sum over all the pixels of the images. This quantity is measured in decibels, and higher values correspond to a better denoising. The parameter $\beta$ was varied in the method to obtain the highest possible values of PSNR.

Beside the gaussian-noise affected images, the XRD experimental setups collect data on two-dimensional arrays, where a Poisson counting statistics accordingly needs a noise reduction.

Wide Angle X-ray Diffraction (WAXD) patterns are used to study the crystalline structure of materials, including fibers [6]. When fibers are analyzed using WAXD, the diffraction patterns provide information about the spacing and orientation of the crystalline regions within the fiber. In a typical WAXD pattern for fibers, one observes:

• diffraction peaks, that indicate the presence of crystalline regions while their position and intensity can help determine the fiber’s crystallinity and the arrangement of its molecular chains;
• intensity distribution of the peaks in the pattern, that can show how the crystalline regions are distributed and oriented relative to the fiber axis;
• azimuthal scans, that can be used to study the orientation of the crystallites around the fiber axis, revealing information about fiber alignment and texture.

By analyzing these patterns, one gains insight into the structural properties of the fiber, such as its crystallinity, orientation, and potential mechanical properties. Thus, in order to investigate the present approach in such noisy data, the model is applied to denoise XRD patterns collected on collagen molecules in tendon-derived collagens, where, under a number of biochemical conditions, a superorganization of tissue into triple helices (with a preferred orientation displayed by the $\pi$-symmetric partial arcs of Figure 2 replacing the fully $2\pi$-symmetric circles) and a high crystalline domain can trigger mechanical stiffness [5]. In Figure 2 the noisy XRD pattern (central panel) has been integrated on a collection time of 2400s, resulting in a maximum counts below 1000, while the noiseless one (left panel) required 327350s, with maxima reaching some 60000. The input image was processed for several values of the parameter $\beta$. The results for PSNR (initial and final values) for the optimal $\beta$ are illustrated in Table 1, while the algorithm outcome is shown in Figure 2 (right panel). The PSNR has trends similar to the case of Figure 1 with gaussian noise: it can be seen that the PSNR stays constant for the optimal $\beta$ value, after an initial stagger (essentially similar to the PSNR vs time plot, as a funtion of $\beta$, illustrated in [7], although the LBM-inspired algorithm outperforms the one used in the reference7). Another example of WAXD fiber pattern denoising is illustrated in Figure 3, where the close-up images on the upper right corner show the improvement of the diffraction arc/ring details. The improvement of the corresponding figure of merit (PSNR) is reported on Table 1.

The LBM-based denoising scheme can also be fruitfully applied to the XRD patterns collected at the synchrotron radiation facilities. The Pilatus II is a single photon counting pixel detector, e.g. used at the PSI facility [22], renowned for its high spatial resolution and efficiency in detecting low-intensity X-ray beams. It features an advanced hybrid pixel array that counts individual photons, offering exceptional precision and dynamic range. This detector is pivotal for experiments requiring detailed imaging and analysis, such as those in structural biology and materials science. Its performance is enhanced by its low noise and high count rate capabilities, making it a versatile tool for cutting-edge research in various scientific fields. However, low statistics two-dimensional images, sometimes affected by artifacts (e.g. experimental setup shadows), could benefit from a pre-processing denoising: in Figure 4 the application of the LBM-inspired denoising algorithm is illustrated in such a case, while the slight improvement of the PSNR is still appreciated in Table 1. In that respect, in order to stress the algorithmic efficiency, the values of another figure of merit are added, namely the structure similarity (SSIM), as defined for instance in [23] and implemented in MatLab: PSNR/SSIM values increase for all cases analysed in this paper; they are in agreement with the visual inspection of the images before/after applying the algorithm, as reported in Figure 1, Figure 2, Figure 3 and Figure 4.

In summary, while LBM-based image denoising is still an evolving field, it shows significant potential due to its ability to model complex noise processes and preserve the quality of important features in the images. Moreover, this study highlights its effectiveness in dealing with specific types of noise, such as gaussian noise or poissonian noise. Finally, advances in computational techniques and hardware have enhanced the efficiency of LBM-based X-ray diffraction image denoising: on one hand improved algorithms and parallel processing have made it feasible to apply these methods to high-resolution images and real-time processing scenarios, on the other hand its ability to model diffusion and transport processes leverages the inherent properties of LBM to effectively remove noise on a well established multiscale and adaptive Physics-grounded approach. Ongoing research and development are likely to further enhance its capabilities and broaden its applications. Studies have compared LBM-based denoising methods with traditional approaches such as gaussian filtering, wavelet denoising, and non-local means: LBM methods often show improved performance in preserving edge details and structural information while effectively reducing noise. Furthermore, there is growing interest in combining LBM with machine learning techniques for image denoising: machine learning models are used to optimize LBM parameters or to refine denoising results, leading to improved performance and more adaptive solutions. Finally this versatile method could be efficiently embedded in a suite of integrated programs for X-ray imaging, e.g. in [24].