Online dynamic flat-field correction for MHz microscopy data at European XFEL

The high pulse intensity and repetition rate of the European X-ray Free-Electron Laser (EuXFEL) provide superior temporal resolution compared with other X-ray sources. In combination with MHz X-ray microscopy techniques, it offers a unique opportunity to achieve superior contrast and spatial resolution in applications demanding high temporal resolution. In both live visualization and offline data analysis for microscopy experiments, baseline normalization is essential for further processing steps such as phase retrieval and modal decomposition. In addition, access to normalized projections during data acquisition can play an important role in decision-making and improve the quality of the data. However, the stochastic nature of X-ray free-electron laser sources hinders the use of standard flat-field normalization methods during MHz X-ray microscopy experiments. Here, an online (i.e. near real-time) dynamic flat-field correction method based on principal component analysis of dynamically evolving flat-field images is presented. The method is used for the normalization of individual X-ray projections and has been implemented as a near real-time analysis tool at the Single Particles, Clusters, and Biomolecules and Serial Femtosecond Crystallography (SPB/SFX) instrument of EuXFEL.

1. Introduction

Over the past few decades, X-ray microscopy has been established as an important tool in a broad range of scientific applications. In addition to spatial resolutions reaching the nanometre scale (Sakdinawat & Attwood, 2010), there have been significant advances in temporal resolution at third- and fourth-generation X-ray synchrotron sources (Fezzaa & Wang, 2008; Olbinado et al., 2017; Parab et al., 2018). Fast megahertz (MHz)-rate X-ray microscopy synchronized to individual X-ray pulses has been pioneered at the Advanced Photon Source, USA (Fezzaa & Wang, 2008), and at the European Synchrotron Radiation Facility, France (Olbinado et al., 2017), which made it possible to image fast stochastic phenomena. Recently, MHz X-ray microscopy with frame acquisition synchronized to individual X-ray pulses has been demonstrated at the European X-ray Free-Electron Laser (EuXFEL), Germany (Vagovič et al., 2019). This development was enabled by the unique properties of X-ray free-electron lasers (XFELs) and the availability of fast MHz detectors. In particular, EuXFEL is currently the only XFEL source providing intense pulses arriving at MHz rates (Decking et al., 2020), which opens up the potential for development of novel methods and applications, such as MHz X-ray multi-projection X-ray imaging (Villanueva-Perez et al., 2023).

EuXFEL provides three orders of magnitude more photons per pulse (∼10¹²–10¹³) than the most brilliant synchrotron source, i.e. ESRF. EuXFEL also provides high spatial coherence, with a source distance of 1 km and a maximum divergence of 4 µrad, with photon energy up to 24 keV (Sinn et al., 2011). The unique properties of EuXFEL come with unique issues, which complicate the acquisition of high-quality MHz microscopic data. An XFEL's X-ray pulses are generated through the self-amplified spontaneous emission (SASE) process (Milton et al., 2001), which leads to strongly stochastic spatial and temporal fluctuations.

The noise appearing in microscopy data is often referred to as fixed-pattern noise. This noise originates from the different responses of the detector's pixels and imperfections in the optic systems, such as scintillator screens, lenses, etc. Ideally, image aberrations should be corrected, so that only the signal variations originating from the X-rays' interaction with a sample are observed. Methods aiming to achieve this goal are often referred to as `flat-field correction' methods.

Flat-field correction methods rely on using reference images averaged over many realizations of acquired images. Reference images acquired with beam illumination are referred to as `flat-field images', while images acquired without beam illumination are referred to as `dark-field images'. However, as a stationary method, traditional flat-field correction, which assumes a flat-field that is stable in time, is unable to sufficiently normalize data produced by an XFEL probe which is stochastic and varies in a pulse-to-pulse manner.

There are several methods used for correcting both these effects (Nieuwenhove et al., 2015; Hagemann et al., 2021; Buakor et al., 2022; Hodge et al., 2022). The issue of normalization of a dynamically changing flat-field has been tackled by a method (Nieuwenhove et al., 2015) in which the reference flat-fields are represented in a latent lower-dimensional vector space and an `effective' flat-field is chosen per input image in the latent space. A latent space is a vector space resulting from an injective transformation of the original data representation, in which similar samples are positioned closer to one another, as measured by a certain criterium. In this paper, the latent space is built in such a way that strongly correlated data variations are agglomerated in the same dimension. This method, called dynamic flat-field correction, was introduced by Nieuwenhove et al. (2015) for synchrotron data. The application of a similar method was demonstrated on EuXFEL's full-field imaging data (Hagemann et al., 2021) of high magnification and in combination with deep-learning approaches (Buakor et al., 2022).

In this work, we present an online method used at EuXFEL to perform dynamic flat-field normalization for MHz microscopy data. By `online method' we mean a method used during data acquisition with near real-time output. EuXFEL users often require fast feedback when performing an experiment, such that they may adapt their configurations to improve the quality of the data taken. In such conditions, obtaining normalized images during the experiment may be the deciding factor between a successful experiment and failure to obtain reliable data. The online flat-field correction method described in this paper is optimized to provide output to the user as data arrive, minimizing the delay between data acquisition and analysis. The method proposed here also provides an indication of whether the flat-field images collected as references are sufficient, or whether a new dataset must be collected due to changes in the experimental conditions. This helps to enhance the method's performance, as a dataset of flat-fields that sufficiently describe the illumination at a given moment of measurement has a crucial impact on the resulting quality of the normalization when affected by the fluctuations of SASE illumination. Moreover, the normalized data provided during experiments can be further re-used as input to other processing tools, for which pre-processed data are essential. Our main objective is to ensure high-quality image data at high speeds. We are using a software tool provided by EuXFEL, enabling us to use online data and follow with the normalization procedure in near real-time. The resulting online normalization tool for MHz microscopy was developed and implemented at the Single Particles, Clusters, and Biomolecules and Serial Femtosecond Crystallography (SPB/SFX) instrument of EuXFEL. Our online implementation based on the dynamic method led to a major improvement in the quality of corrected images compared with results using the conventional flat-field-correction approach.

The manuscript is organized as follows. In Section 2, we describe flat-field correction methods. The normalization algorithm and implementation are described in Section 3. We discuss the results and the performance of implementation using test data in Section 4, followed by conclusion and discussion given in Section 5.

2. Flat-field correction

In this section we briefly review the conventional flat-field correction method. We then expand on the dynamic flat-field correction method, which is implemented in this work. While we closely follow Nieuwenhove et al. (2015) and Buakor et al. (2022), differences in the methods are pointed out as needed.

A stationary, conventional flat-field correction removes the effects of spatially uneven illumination functions. Such effects may be caused by variations in the detector response in different pixels or modules, or effects arising from the scintillator and optics configurations. One assumes that there has been a previous data acquisition of N_f flat-field images f_j and N_d dark-fields d_k, where j ∈ [1, N_f] and k ∈ [1, N_d] refer to the reference images' index. In a stationary flat-field correction, one takes advantage only of the average flat-field and dark-field images, = and = .

A simple procedure to correct the uneven detector response would be to correct each raw input image p_i by calculating

As mentioned previously, due to the character of SASE processes, noise is shot-to-shot-dependent, introducing a change in the intensity profile on the scintillator screen between shots. This behavior cannot be captured solely by the average flat-field image, causing the conventional method to fail for such conditions. To correct for the stochastic effects observed in SASE sources such as EuXFEL, a dynamic flat-field correction substitutes the average flat-field with an effective flat-field , that depends on the collected sample image p_i itself. With this alternative procedure, the corrected image is

In principle, the true illumination function may be any image. However, some simplifying assumptions allow for an approximate estimate of . Firstly, it is natural to assume that such a flat-field could be expressed as a function of the previously collected flat-fields f_j. While such a function may be extremely complex, a simple first-order approximation would be that is a linear combination of the collected flat-fields. For the purposes of an online flat-field correction, this initial assumption leads to a fast and reliable implementation, which may be supplemented by further corrections offline.

While one may approach the issue of the approximation of an effective flat-field with as many weights as there are collected flat-fields, one may further simplify such a description by assuming that only a few collected flat-fields are the main contributors to the sum. A data-reduction technique could, therefore, be of use to reduce the number of coefficients required. We assume that the flat-fields, viewed as a multidimensional random variable with each dimension represented by the pixel content, are samples of a multivariate normal distribution. Note that the assumption is not that the pixels are independent samples of a Gaussian random variable but that the entire image is sampled from a normal distribution, allowing us to consider the correlation between pixel values as appropriate. Under such conditions one may use principal components analysis (PCA) (Pearson, 1901; Hotelling, 1936) to rewrite the linear flat-field combination as a linear combination of the principal components and keep only the components with the largest variance. Discarded components would then contribute little to the linear combination. In general, it is feasible to use M components, where M   <<  N_f. In such a setting, we expand the effective flat-field as

where u_m is the mth principal component of the mean-subtracted flat-fields, , and w_im are free parameters to be identified online.

While a functional form is available to parameterize the effective flat-field , one must determine a procedure to select the weights w_im such that the effective flat-field is uniquely chosen for a given sample image p_i. A constraint must be imposed on the corrected image for a meaningful choice of weights.

A key constraint is to preserve details originating from the interaction of X-rays with the sample. Such details appear through sharp transitions or edges in the image, separating the signal from the background content. The free parameter w_im must, therefore, preserve such edge effects, while removing spurious variations in the images caused by the uneven detector response. Additionally, one would like to remove spurious variations in the images added by an uneven flat-field. The regularization technique of total variation (TV) denoising is often used in image processing to achieve exactly those goals, by minimizing the TV of the corrected image, while aiming to obtain an image as similar as possible to the uncorrected one (Chambolle, 2004). Details concerning the calculation of weights w_i can be found in the literature (Nieuwenhove et al., 2015; Chambolle, 2004).

The method described here consists of two separate steps. Initially, reference flat-fields and dark-fields are acquired and PCA is used to obtain the most relevant principal components of the flat-field dataset. During data acquisition with a sample, for each individual frame, weights appropriate for the effective flat-field are found. After such calculation, the final corrected frame is obtained by applying the weights in the linear combination in equation (2).

Such a procedure must be implemented in a way that it can be run online at adequately fast rates to keep up with the data-acquisition rate of EuXFEL. We expand on the technical implementation in the next section.

3. Implementation

This section provides a schematic overview of the two steps of the main online normalization algorithm outlined in the previous section. Selected details of the implementation of the algorithm and specific features for visualization applied at EuXFEL are detailed here.

A schematic overview of the first part of the algorithm, employing PCA, can be found in Fig. 1(a), whereas Fig. 1(b) summarizes the second part which estimates an individual flat-field for every sample image and visualizes the resulting normalized data.

Figure 1 Algorithm overviews of the first part with principal component analysis of the flat-field dataset (a) and the second part performing the online dynamic flat-field correction and visualization (b).

Next, we list some of the specific functions and details used in the implementation of the algorithm. For the total variation minimization of the objective function a quasi-Newton algorithm is applied (Byrd et al., 1995). Python implementation of the minimization algorithm (Byrd et al., 1995) contains the `factr' parameter, which controls the accuracy and is referred to as the precision parameter here. In this work, the precision parameter was tuned to maximize the processing rates of the correction algorithm. Adjustments have been made to preserve the quality of the corrected images while increasing the processing speed.

The speed of the algorithm is further increased by down-sampling all the images entering the minimization procedure (Nieuwenhove et al., 2015; Buakor et al., 2022), causing also reduction of noise in the corrected images. Here, in order not to introduce more time-consuming steps, we omitted any additional algorithms aiming exclusively for noise reduction, which may be included at later stages for offline processing.

This work demonstrates a novel implementation of the dynamic flat-field correction algorithm to be used in EuXFEL's infrastructure, which can be used as a normalization technique during experiments. Moreover, considering the goal of user-friendliness of the online normalization tool, EuXFEL's software framework Karabo (Heisen et al., 2013) is used here for online visualization. One of the two implementations uses the metropc framework and the extra-metro package (Fangohr et al., 2018), which is integrated within the Karabo framework and its graphical user interface (GUI). The metropc framework allows for flexible adjusting of a script during an experiment. It allows visualization of the corrected incoming data, in addition to raw sample images. The second implementation of the online dynamic flat-field correction algorithm is done using Karabo bridge (Fangohr et al., 2018) to access online data in the Karabo pipeline and perform parallelized analysis and further visualization using the Qt toolkit (PyQT, 2012).

4. Results

We have described here the results obtained by our dynamic flat-field correction implementation at EuXFEL. We demonstrated our algorithm on a Venturi tube dataset (Soyama & Hoshino, 2016). Data were acquired in March 2021 at the SPB/SFX instrument (Mancuso et al., 2019) of EuXFEL with a photon energy of 9.3 keV. A fast-frame-rate Shimadzu HPV-X2 camera, coupled to a 8 µm-thick LYSO:Ce scintillator, was used as an imaging system with a 10× magnification (Vagovič et al., 2019). Datasets consisting of 68 dark-field and 67 flat-field trains were taken to test our method. Each train consisted of 128 images with a size of 250 × 400 pixels. The effective pixel size was 3.2 µm. The exposure time of each image was given by the X-ray pulse duration, approximately 25 fs (Sinn et al., 2011), translated into the latent image at the scintillator. The duration of the signal emitted by the scintillator was given by the scintillator emission time τ (LYSO:Ce, τ ≃ 40 ns). The frequency of acquired frames within a train was 1.128 MHz and trains were collected at frequencies from 0.08 to 0.1 Hz, limited by the camera idle time. The performance of the dynamic flat-field correction algorithm on the aforementioned dataset is demonstrated in this section.

The first step toward the normalization of online data is the processing of flat-field and dark-field data. Our aim was to acquire hundreds to thousands of flat-field images, which, at the current acquisition speed, requires up to 10 min of measurements. After the first step of the algorithm, as described in Fig. 1(a), a preview of the cumulative sum of explained variance ratio of principal components can be viewed. This sum reflects how much variance in the data is captured, i.e. explained, by a different number of principal components. A full set of principal components, whose size equals the number of flat-field images, leads to a value of one for the cumulative sum of explained variance ratio, thus capturing 100% of the variance in the flat-field dataset. This allows monitoring of the rate at which the selected number of principal components explain variance in the flat-field dataset. In addition, it enables the increase or decrease of the number of components before proceeding with the normalization of sample images.

The flat-field dataset consists of approximately 10000 flat-field images with a few instances shown in Fig. 2(a). Figure 2(b) shows the mean flat-field image calculated over the whole dataset. In Fig. 3, selected principal components are depicted for the index of principal component m = {1, 2, 3, 16}. The first principal component, labeled m = 1, is similar to the mean flat-field image in Fig. 2(b) and explains the largest fraction of the variance in the dataset. The two components m ∈ {2, 3} capture the intensity variation in the horizontal and vertical directions. Components m ≤ 4 show mostly horizontal changes in intensity with pronounced stripe features, as depicted on the component m = 16. The cumulative sum of explained variance ratio of the flat-field dataset is shown in Fig. 3 for the first 25 principal components.

Figure 2 Example instances (a) of the test flat-field dataset consisting of approximately 10000 images and the average flat-field image (b) calculated over the whole test dataset. All images shown here have size (250, 400) pixels leading to a 0.8 mm × 1.28 mm field of view.

Figure 3 An example of principal components um for m = {1, 2, 3, 16} (a), obtained from the test flat-field data shown in Fig. 2. The cumulative sum of explained variance ratio of the first 25 principal components is depicted in (b).

Figure 3(b) shows that the first three components explain approximately 83% of the variance. Adding more components does not increase it significantly, and after component 3 the increases in cumulative sum of explained variance are minor. As discussed in more detail later in this section, in order to correct for the majority of vertical features in sample images one needs to include also those components with seemingly low values of explained variance, as they obtain variations needed for a high-quality correction.

After identifying the most important components describing flat-field illumination, one can continue with the normalization procedure of sample images described in Section 3. To better visualize the performance of the dynamic flat-field correction on the test data, in Fig. 4 we plot the sample image (a), and its counterparts corrected by both conventional (b) and dynamic (c) methods. The default value for the number of principal components taken into consideration here is 20.

Figure 4 Comparison of a sample image (a) without any normalization, conventionally flat-field-corrected image (b) and dynamic flat-field-corrected data (c). Only the bottom section of the images [(50, 400) pixels and 0.16 mm × 1.28 mm field of view] is shown here from the originally sized images [(250, 400) pixels and 0.8 mm × 1.28 mm field of view].

Figure 4 illustrates that the dynamic method (c) more successfully removes vertical features than the conventional correction method (b). We can notice vertical features along the whole width of an image changing position from frame to frame, which is captured in the examples shown in Fig. 2(a). Moreover, the lack of vertical features on the average flat-field image, Fig. 2(b), leads to the inability of conventional methods to successfully remove such artifacts from the acquired images. Even fine stripes, appearing on an average flat-field, left side in Fig. 2(b), are not removed using the conventional normalization method, shown in Fig. 4(b), due to the slight change of their position during measurement. Another challenging detail to eliminate is the circular spot from the middle of the sample image (a). Similar to the fine stripes from the previous example, better, however not perfect, reduction of the spot is reached only with the dynamic method, Fig. 4(c).

The quality of the flat-field-corrected images was mostly assessed visually. To estimate the method's performance quantitatively, we calculated the pixel value spread (pvs) as a deviation of an image from its average pixel value,

where index l iterates over pixels of image p. We were using only `flat' parts of images, without any sample, in order to estimate the amount of noise present in an image. Higher pvs values are expected for uncorrected images with uneven intensity profile, and lower values for images corrected in a way that reduces most of the uneven structures and noisy pixels.

In order to assess the quality of correction of images of the same position within a train we calculated TV and pvs for each of 128 images acquired in one train and averaged over 15 trains. Average values for both variables and their standard deviations are shown in Fig. 5. While the pvs and TV values vary for each frame, one may see in Fig. 5 that the dynamic flat-field-corrected images have a consistently lower TV (a) and spread (b), suggesting a lowering of noise level in the chosen sample-free part of the images. In contrast, the conventional flat-field correction shows very little benefit, in comparison with no correction at all, as for TV and pvs values no statistically significant difference was found for uncorrected images and images corrected using the conventional method. This result was expected based on the visual comparison in Fig. 4, in which no visible removal of vertical features was observed. The standard deviation expressed as a colored area for uncorrected images and both normalization methods in Fig. 5 describes train-to-train changes in TV and pvs values for a given frame in the dataset. On average, the pvs value for a dynamically corrected image was found to be 0.069, while the spread for conventionally corrected images was 0.098. Uncorrected images have a pvs of 0.102.

Figure 5 Comparison of total variation (a) and pixel value spread (b) for uncorrected images and images corrected by conventional and dynamical methods. Values of TV and pvs are calculated depending on their position within a train, where the maximum frame number is 128. Both TV and pvs for each frame number are averaged over 15 trains and train-to-train changes are captured by their standard deviation expressed as the colored area.

Additionally, we studied the quality of normalized images estimated by TV values as a function of a varying number of principal components. Resulting corrected images were also assessed visually. Figure 6 describes changes in the TV value due to a varying number of principal components used in the normalization procedure. One can notice a decrease in the average value of TV for numbers of principal components 15 and 20, followed by an increase starting after the 20 components mark. For a visual comparison, Fig. 7 provides examples of flat-field-corrected images using a varying number of principal components in the normalization algorithm. Figure 7 shows that including two to ten principal components in the normalization procedure leads to non-optimal flat-field correction results. Fine line structures are still visible on the left part of the images for the two- and four-component cases. Vertical features along the whole width are not fully eliminated using ten principal components. Increasing the number of principal components up to ten results in decreased visibility of vertical features, which also supports the observation of decreasing TV values in Fig. 6. The circular spot in the middle is becoming less noticeable with a higher number of principal components, although can be detected on all corrected images up to the case using 100 components. Including around 20 principal components seems to remove the majority of the unwanted features and results for 40 and 100 components in Fig. 6 are not leading to visibly better results. TV values in Fig. 6 are growing after 20 components, which may be caused by overfitting. The results listed in Figs. 6 and 7 indicate that the optimal number of principal components is around 20, where most of the features appearing on flat fields are removed.

Figure 6 Dependence of TV of images corrected by the dynamical method using a different number of principal components. Values are calculated for one train and are averaged over 128 images with the gray colored area given by standard deviation of frame-to-frame changes within a train.

Figure 7 Comparison of the uncorrected sample image and flat-field-corrected images using the dynamical method with [2, 4, 10, 20, 40, 100] principal components. Only the bottom section of the images [(50, 400) pixels and 0.16 mm × 1.28 mm field of view] is shown here from originally sized images [(250, 400) pixels and 0.8 mm × 1.28 mm field of view].

Considering the idle time of the camera, our goal was to reach a correction before the next acquired train arrives. The idle time of the camera is around 15 s, meaning that the camera is able to acquire one in 150 of EuXFEL's trains. This consideration does not take into account the setup time spent for the analysis of flat- and dark-field datasets, as this process usually takes a few minutes and needs to be repeated on an hourly basis, whenever new flat- and dark-fields are taken. The proposed algorithm and both its implementations mentioned in Section 3 were able to be faster and keep up with the volume of incoming data. A parallelized version reached a correction time of approximately 1 s for a dataset consisting of one train with 128 images. Tests were performed at online cluster INTEL, Gold-6140 CPU @ 2.30 GHz, 768 GB. Moreover, the processing speeds should be able to sustain the load with anticipated updates, which aim to shorten the idle time of the camera. The current performance of the algorithm is shown to be suitable for use during experiments as a near real-time analysis tool.

5. Conclusion and discussion

We have demonstrated a method for dynamic flat-field correction on MHz-repetition-rate EuXFEL data. It has been shown that, for X-ray MHz microscopy data (Buakor et al., 2022), the dynamic method results in an improved correction compared with the conventional flat-field correction method, taking into account temporal variations in the flat-field dataset, which is characteristic of XFEL facilities. To compare the performance of conventional and dynamic methods applied to test data, we have calculated total variation and the pixel value spread for the area of images without any sample. Both quantitative comparison and visual assessment have shown improved results with the dynamic normalization method compared with the conventional method. The normalization algorithm has been implemented online to visualize data and flat-field-corrected data during experiments at the SPB/SFX instrument of EuXFEL. The implementation was able to correct and visualize incoming data in near real-time before the next data from the MHz-frame-rate camera arrived. In future, we plan to extend the use by implementing new data analysis tools such as data decomposition methods and adding more features to the GUI.