Unveiling water dynamics in fuel cells from time-resolved tomographic microscopy data

X-ray dynamic tomographic microscopy offers new opportunities in the volumetric investigation of dynamic processes. Due to data complexity and their sheer amount, extraction of comprehensive quantitative information remains challenging due to the intensive manual interaction required. Particularly for dynamic investigations, these intensive manual requirements significantly extend the total data post-processing time, limiting possible dynamic analysis realistically to a few samples and time steps, hindering full exploitation of the new capabilities offered at dedicated time-resolved X-ray tomographic stations. In this paper, a fully automatized iterative tomographic reconstruction pipeline (rSIRT-PWC-DIFF) designed to reconstruct and segment dynamic processes within a static matrix is presented. The proposed algorithm includes automatic dynamic feature separation through difference sinograms, a virtual sinogram step for interior tomography datasets, time-regularization extended to small sub-regions for increased robustness and an automatic stopping criterion. We demonstrate the advantages of our approach on dynamic fuel cell data, for which the current data post-processing pipeline heavily relies on manual labor. The proposed approach reduces the post-processing time by at least a factor of 4 on limited computational resources. Full independence from manual interaction additionally allows straightforward up-scaling to efficiently process larger data, extensively boosting the possibilities in future dynamic X-ray tomographic investigations.

Scientific RepoRtS | (2020) 10:16388 | https://doi.org/10.1038/s41598-020-73036-w www.nature.com/scientificreports/ However, as for the other methods described above, several parameters need to be tweaked in order to push the reconstructed image quality to its maximum. While the rSIRT-PWC algorithm offers an interesting possibility to exploit the continuity information across the available dynamic time sequence through time-regularization, it still relies, in a first step, on prior segmentation of the solid matter to separate the stationary region from the areas allowed to change through time (dynamic), which are continuously estimated and updated throughout the reconstruction process. For datasets with limited SNR, such as operando X-ray tomographic microscopy of fuel cells, extracting these dynamic and stationary regions from the reconstructed volumes and aligning the corresponding masks in a nearly perfect manner requires significant amount of manual labor. Moreover, as these sub-second dynamic experiments deliver as many as 10 tomograms per second for extended time periods, processing each of these large datasets prior to the algorithm application is infeasible.
To overcome these limitations, we propose a rSIRT-PWC-DIFF algorithm extended for interior tomography datasets, designed to reconstruct and segmented noisy datasets in a highly automatized manner enabling efficient scalability to large data volumes. The proposed algorithm operates on difference sinograms to separate the static and dynamic regions of the sample in a fully automatized manner. As no prior information of the static pixel location is needed, the proposed algorithm can be applied to reconstruct several datasets without manual interaction, making it ideal for cases where prior segmentation of the static features is difficult. In addition, the robustness of the time-domain function fitting step for highly noisy datasets is increased by considering small sub-regions (5 × 5 pixels) with a Gaussian weighting scheme instead of single pixels. Though the proposed algorithm protocol has been developed to overcome the current restrictions specifically in fuel cell data reconstruction, its application is not limited to fuel cells but can be applied to any dataset where evolving dynamics are present.

Methods
Background. Tomography model. The underlying tomography model is introduced for a 2D case assuming parallel beam geometry. Its extension to 3D is trivial.
The object of interest, represented on a pixel grid of N pixels, is denoted by a column vector x = x j ∈ R N . A vector p = p i ∈ R M is a collection of the log-corrected measured projection values or phase-retrieved projections of the object x from all measured angular positions, typically distributed homogeneously between 0 and 180 degrees. The projection data p and object x are connected by p = Wx where W = w ij ∈ R M×N is a collection of weights that models the contribution of each pixel j to the projected value at index i.
Due to the ill-posedness of the system, it is typically infeasible to directly solve x from the system of linear equations. Analytical reconstruction methods, such as filtered back-projection (FBP) 27 , are accurate when enough angular views of the object are available. Limited angular views or high sparsity in angular sampling lead instead to streak artefacts deteriorating the reconstruction quality. Iterative reconstruction algorithms address these limitations by exploiting prior information of the sought object, while minimizing the difference between the forward-projected estimated reconstruction and the measured projection data in an iterative manner.
Conventional SIRT. The Simultaneous Iterative Reconstruction Technique (SIRT) 28,29 is an algebraic iterative reconstruction algorithm where the object x is considered to consist of an array of unknowns, represented by algebraic equations. Starting from an initial reconstruction x (0) , typically a zero vector, the SIRT algorithm updates the reconstruction at iteration k by where R = r ij ∈ R M×M and C = c ij ∈ R N×N are the diagonal matrices with the inverse row and column sums of W , respectively. The algorithm is known to converge to a solution of Conventional rSIRT-PWC. Van Eydhoven et al. introduced in 2015 an iterative reconstruction algorithm developed to follow in 4D dynamic fluid flow through a solid matrix 26 . The region-based SIRT with intermediate piecewise constant function estimation (rSIRT-PWC) algorithm, based on conventional SIRT, exploits prior knowledge of the dynamic and static regions of the sample to improve the spatial resolution of the reconstructions, resulting in improved spatial and temporal resolution of the experiment. The static and dynamic regions of the object are separated through masking (often requiring segmentation). The static regions are reconstructed using all available data across the complete time sequence, so to increase the sampling frequency, and a static mask is applied to extract the static pixels from the reconstructed high-quality image. Simultaneously, each of the dynamic time frames is reconstructed separately and a dynamic mask is applied to extract the respective dynamic pixels.
To further improve the dynamic reconstructions, suffering from high noise level and undersampling artefacts, the algorithm exploits piecewise constant functions (PWC) to model each single dynamic pixel's attenuation curve over the complete time sequence, so incorporating time regularization within the reconstruction procedure. Through the PWC fitting, noise is effectively suppressed and the reconstruction accuracy of the dynamic region is significantly improved 26 . These dynamic and static reconstructions are finally merged for each time frame, obtaining a full 4D reconstruction of the dynamic object in high-quality.
In addition to masks, iteration parameters have to be manually selected. This algorithm foresees nested loops of rSIRT reconstructions and PWC time-regularization steps. It is therefore necessary to select the total Scientific RepoRtS | (2020) 10:16388 | https://doi.org/10.1038/s41598-020-73036-w www.nature.com/scientificreports/ number of iterations as well as the number of rSIRT iterations performed prior to the first PWC fitting and the number of rSIRT iterations performed between each of the following PWC steps in addition to the total number of PWC steps. Depending on the sample, tweaking these parameters can have a significant impact on the reconstruction quality while it is typically difficult to know a priori which set of parameters would lead to the optimal reconstruction. rSiRt-pWc-Diff. While the conventional rSIRT-PWC algorithm has shown superior reconstruction performance on a neutron tomography dataset compared to the standard analytical reconstruction (FBP), the requirement of static and dynamic masks is a major limiting factor for datasets for which such masks are difficult and/or time-consuming to obtain. In addition, several iteration parameters need to be specified prior to reconstruction. In case of many in situ time-resolved experiments as fuel cell imaging, an adaptation to accommodate interior tomography datasets is also necessary. To overcome these limitations, an adapted rSIRT-PWC-DIFF algorithm was developed. A detailed description of the proposed algorithm is presented in the following sub-sections.
Difference sinogram approach. In the proposed rSIRT-PWC-DIFF algorithm, dynamic and static masks are no longer required, as the dynamic and static pixels are separated through an automated sinogram subtraction step. If a static scan or a static time frame of the object is available, the corresponding static sinogram S S can be subtracted from the dynamic sinogram S Dn for each time step n, as The resulting difference sinograms S Diff n contain only dynamic changes and noise for each time frame n and can be further used in the rSIRT-PWC-DIFF algorithm to reconstruct dynamic changes for each time frame n . Simultaneously, the dry sinogram S S is considered to reconstruct the static regions of the object. If the number of angular views between the static sinogram S S and the dynamic sinograms S Dn is not equal, interpolation can be applied to the sinograms prior the subtraction step.
Sinogram alignment. Ideally, angular information for each projection image is available and can be directly used to align the static and dynamic sinograms prior to the subtraction step. To generalize the proposed algorithm for datasets without angular information, a cross-correlation-based sinogram alignment step was implemented. In the sinogram alignment step, for each time frame n maximum cross-correlation 30 ( * ) between the dynamic sinogram S Dn and the static sinogram S S is computed as The maximum correlation position is used to shift the dynamic sinograms S Dn so to align them with respect to the static sinogram S S . This alignment process is fully automatized within the proposed algorithm. If the angular information is available, it can be applied directly to align the sinograms.
The cross-correlation step identifies the main structural similarities between the sinograms from the superimposed data and among the tested registration approaches produces the best alignment between the sinograms. For other types of samples, alternative registration methods based for instance on mutual information 31 could also potentially be exploited. Any possible mismatch between the static structures of the dry and dynamic datasets, for example due to membrane swelling ("Effects of membrane swelling" section), translates as artefacts in the subtracted sinograms. When a standard analytical reconstruction algorithm, such as FBP, is applied to reconstruct these difference sinograms, misalignment artefacts become visible in the reconstructed images. However, when an iterative reconstruction algorithm, such as SIRT, is considered instead, these misalignment artefacts are partly suppressed during the iterative forward-and back-propagation steps and so, not significant in the reconstructed images.
Virtual sinogram for interior tomography. To increase the versatility of the proposed algorithm, further developments were made to extend it to interior tomography (INT) problems. In INT, the object support is fully or partly outside of the field of view (FOV) of the detector. This geometry often occurs in medical imaging, material science and biology applications where high spatial resolution information of a smaller region of interest (ROI) is required. If INT projections are reconstructed with standard iterative or analytical algorithms without special consideration, the reconstructions will suffer from strong artefacts 32,33 , compromising any further analysis of the reconstructed volume. To minimize these artefacts, a fully automatized virtual sinogram step 34 was incorporated within the reconstruction algorithm.
In the virtual sinogram step, the INT sinogram is first edge-padded and a standard FBP reconstruction is performed. The pixels outside of the resolution circle (the largest circle fitting into the reconstruction grid) are then suppressed to zero to create an artificial object boundary. This modified object, now possessing an artificial compact support, is further forward projected to create a virtual sinogram, simulating now a standard sinogram where the object and its boundary are completely within the FOV. This virtual sinogram, further edge-padded, can then be used for iterative reconstruction.
Curve fitting on sub-regions. The curve fitting (PWC) part of the conventional approach was further developed to improve its robustness for highly noisy and slightly misaligned datasets. Instead of the gray level value of single pixels, the 2D Gaussian weighted average value of small sub-regions of 5 × 5 pixels was used to emphasize the Scientific RepoRtS | (2020) 10:16388 | https://doi.org/10.1038/s41598-020-73036-w www.nature.com/scientificreports/ direct neighborhood of each pixel during the fitting procedure. The image borders were padded prior to weighting. While the sub-region size can be adjusted, it is important to note that larger regions may lead to reduction in size or vanishing of small structures. For datasets with a very limited SNR, this compromise might however be necessary for successful curve fitting and further data processing. Thanks to the Gaussian weighting scheme, in this work, the used 5 × 5 pixel sub-regions did though not negatively impact the spatial resolution.
Automatic stopping criterion. An automatic stopping criterion for selecting the optimal number of iterations was defined. For this purpose, the proposed algorithm was first run for 700 rSIRT-DIFF iterations without applying the function fitting step. After every 10 iterations, the current reconstruction estimate was compared to the reconstruction obtained 10 iterations earlier through the Euclidean L2-norm, defined for images x and y as The normalized L2-norm was plotted as a function of the number of iterations and its gradient was calculated. Several gradient slopes were selected for further analysis (Fig. 1) and the corresponding reconstructions were post-processed by applying a single PWC fitting step. The obtained reconstructed dynamics (here water) were compared to a ground truth segmentation through sensitivity (true-positive-rate), specificity (true-negativerate) 35 and the dice coefficient 36 . The metrics are defined as: Automatic stopping criterion evaluation. Euclidean L2-norm (red curve) between the current and previous reconstruction estimate of PEFC_1 sample ("Real data" in "Materials" section) was measured after every 10 iterations for a total of 700 iterations. The gradient of the L2-curve was estimated and reconstruction quality measured at selected gradient slope positions. The gradient slope (blue curve) within region [−0.001, −0.008] was identified to correspond to the highest reachable reconstruction quality. The droplet diameter is 100 pixels.

Scientific RepoRtS
| (2020) 10:16388 | https://doi.org/10.1038/s41598-020-73036-w www.nature.com/scientificreports/ where TP are true-positive pixels (water reconstructed as water), FP are false-positive pixels (air reconstructed as water), TN are true-negative pixels (air reconstructed as air) and FN are false-negative pixels (water reconstructed as air). Based on quantitative and qualitative analysis ( Fig. 1) on both simulated ("Simulated phantom" section) and real data ("Real data" section), a gradient slope within a region of [−0.001, −0.008] was identified as an optimal stopping point for the rSIRT-DIFF iterations and starting point for the PWC step. In this gradient slope region, the reconstructions recover the sample features successfully without emphasizing noise, leading to a desirable starting point for the time-regularization. Moreover, within this gradient slope region the reconstruction quality remained stable indicating only a low dependence on the iteration number around its optimum. If the iteration number is not high enough, the water features are instead not properly defined. Too many iterations lead to noisy reconstructions. Additional rSIRT-DIFF/PWC fitting step loops as proposed in 26 did not improve the reconstruction quality. In addition to the higher computational load, more loops tend to increase the level of noise in the reconstructed volumes. The described analysis was conducted on a simulated fuel cell phantom with high and moderate noise levels ("Simulated phantom" section) and on two real fuel cell datasets ("Real data" section). For all considered cases, the normalized L2-curve was behaving in a similar manner, only its relative position was shifted in y-direction. In all cases the optimal stopping point was identified to correspond to a gradient slope of [−0.001, −0.008] , which consistently with Fig. 1 was reached with 100-200 iterations. The same optimal gradient slope, obtained for both simulated and real fuel cell data, was found to be independent on image size and contrast. Therefore, for all fuel cell samples the number of rSIRT-DIFF iterations was set, to maximize computational performance, to 100, followed by a single PWC fitting step. For completely different types of samples (for instance sandstone or catalytic materials), the optimal value for the gradient slope has to be confirmed.
Algorithm protocol. To apply the proposed rSIRT-PWC-DIFF algorithm, we propose the following protocol described in Fig. 2: 1. Align and subtract sinograms to extract the dynamic changes ("Sinogram alignment" section).
2. If the dataset was acquired in interior tomography geometry, apply the virtual sinogram step ("Virtual sinogram for interior tomography" section).

Figure 2.
Algorithm protocol flowchart. The algorithm protocol begins by automatically pre-processing the sinograms to extract the dynamic components. If the required number of iterations is not yet set for the current sample type, the protocol moves automatically to the "Automatic initialization" module. Once the number of required iterations has been determined, the pipeline proceeds automatically to the "Automatic reconstruction and feature extraction" step. The set iteration number will be directly re-used for other datasets of similar sample types: the "Automatic initialization" step will be skipped and the protocol will automatically transition from the pre-processing to the reconstruction step. Step 3 has to be performed once for each sample type to determine the desired number of iterations corresponding to the specified gradient slope region. Once the stopping point has been automatically detected, the algorithm can be further applied in a faster manner by directly applying the chosen number of iterations without estimating the L2-norm. If the number of rSIRT-DIFF iterations specified by the gradient slope is used, additional loops of reconstruction and time-regularization steps were not found to improve the reconstruction quality further. Therefore, for optimized performance and computational load, one single cycle of reconstruction and time-regularization is recommended.

Results
In this section, reconstruction results of the rSIRT-PWC-DIFF algorithm are presented for a simulated fuel cell phantom ("Simulated phantom" section) and two dynamic X-ray tomographic microscopy datasets of a fuel cell ("Real data" section). The reconstructions were completed with MATLAB (2018b) using the ASTRA toolbox [37][38][39] and exploiting a Tesla V100 GPU card for accelerated computations.

Simulated phantom. Standard setting.
To test and quantify the developed algorithm performance, a dedicated fuel cell phantom was created (Fig. 3, Section "Simulations" in "Materials"). The difference between Reconstructions obtained with the rSIRT-PWC-DIFF algorithm were compared to standard FBP (Parzen filter) and rSIRT reconstructions of a difference sinogram (FBP-DIFF and rSIRT-DIFF respectively) without the function fitting step. For rSIRT-DIFF and rSIRT-PWC-DIFF, a positivity constraint was enforced by setting negative pixel values to zero during the iterative process. For both iterative methods, a total of 100 iterations was considered based on the defined stopping criterion ("Automatic stopping criterion" section), followed by a single PWC-fitting for rSIRT-PWC-DIFF. The RRMSE metric was measured within three regions: full regionof-interest (ROI), static ROI and dynamic ROI, selected by applying corresponding masks. Results for all three reconstruction methods are presented in Table 1.
The error measures in Table 1 show superior reconstruction performance for all three ROIs when the rSIRT-PWC-DIFF algorithm was applied. In comparison to FBP-DIFF, the proposed algorithm improves the reconstruction of the full ROI by a factor of 6.7. The quality in the static ROI was improved by a factor of 8.8 while in the dynamic ROI by a factor of 3.1. Comparing the results for rSIRT-DIFF and rSIRT-PWC-DIFF, it is evident that including the PWC step improves the reconstruction quality in all three ROIs, and in particular of a factor of 1.13 in the dynamic (water) region. These improvements are also clearly visible in Fig. 3, presenting an example of the reconstructions obtained with the three algorithms and ground truth phantom for the merged and dynamic-only reconstructions. Increased sharpness in the static, fibrous structures (Fig. 3a,b) could be achieved by considering a high-resolution static scan, as typically done in fuel cell imaging experiments ("Real data" in "Materials" section).
Effects of membrane swelling. During fuel cell operation, the hydration of the polymer electrolyte membrane can change, which leads to swelling or shrinkage of the membrane, such that the corresponding bright area within the cell center (Fig. 3) changes its thickness accordingly. Due to this swelling, all dry structures within the GDLs are repositioned by the increasing membrane layer size, causing a mismatch between the dry structures of the dry and wet (dynamic) scans. Typically, 1 to 2 pixels (2.75-5.5 µm) relative shift between some of the dry and wet scans' dry structures can be observed in X-ray tomographic volumes.
To quantify the effects of membrane swelling on the rSIRT-PWC-DIFF reconstruction quality, membrane swelling of 0 to 8 pixels was simulated (Fig. 4a) for the dynamic fuel cell phantom. The water reconstruction accuracy was quantified through sensitivity (true-positive-rate), specificity (true-negative-rate) and the dice coefficient (Section Automatic stopping criterion). The rSIRT-PWC-DIFF reconstructions were obtained by completing 100 iterations of rSIRT-DIFF, chosen based on the stopping criterion (Section Automatic stopping criterion) and followed by a single PWC fitting step. The reconstruction quality was evaluated by comparing the dynamic reconstructions (water) without merging the static structures to the ground truth dynamic phantom (water only) between time steps 10 and 30. Results are presented in Fig. 4b,c.
As demonstrated in Fig. 4b, the sensitivity of the reconstruction remains at 90% or above for mild swelling of maximum 6 pixels, decreasing with increasing swelling. This is also supported by qualitative analysis: Fig. 4c reveals that increasing the applied swelling translates to amplified differences between the reconstructed and ground truth phantom as the misalignment artefacts of the moving static structure boundaries become visible and pronounced, leading to decreased reconstruction accuracy. These artefacts could be partially suppressed by aligning the dry scan independently for the anode and cathode sides of the cell prior subtraction, as it is typically done in practice. The dice coefficient (Fig. 4b) was found to remain between 91 and 88% for all applied swellings. The specificity was found to remain stable at 99% for all considered swelling scenarios.

Real data. The proposed algorithm performance was evaluated on two dynamic Polymer Electrolyte Fuel
Cell (PEFC) synchrotron tomographic datasets described in the "Materials" section. For PEFC_1, the rSIRT-PWC-DIFF reconstructions were computed for a time sequence of 30 consecutive time frames and for PEFC_2 a total of 10 consecutive time frames were considered, based on pre-evaluation of the water dynamics in the data Table 1. Average relative root mean squared error (RRMSE) estimates and their standard deviations for the dynamic fuel cell phantom reconstructions computed using a difference sinogram with FBP (FBP-DIFF), rSIRT (rSIRT-DIFF) and rSIRT-PWC-DIFF. The error metric was computed for the full ROI, the static ROI and the dynamic ROI. Static and dynamic ROIs were separated for evaluation using corresponding masks. The metrics were calculated for 100 fuel cell phantoms, each having randomized fiber and water content.   . For both samples, the manual segmentations were generated by aligning the dynamic Gridrec reconstructions to an additional high-quality dry Gridrec reconstruction of the sample (1000 projections with 1 ms exposure time) followed by an image subtraction step, resulting in images containing only dynamic changes (water) and noise 11 . These subtracted images were further post-processed to generate the final segmentations for each time frame 12 . In addition, for PEFC_1 a static mask containing GDL fiber structures has been applied to the manually segmented images as a post-processing step to exclude any fiber structures misclassified as water. For comparison, the same mask was applied as a post-processing step to the rSIRT-PWC-DIFF reconstructed images. The reconstructions were evaluated through sensitivity, specificity and dice coefficient measures ("Automatic stopping criterion" section), results are presented in Table 2.

FBP-DIFF rSIRT-DIFF rSIRT-PWC-DIFF
As revealed by Fig. 5a, the rSIRT-PWC-DIFF reconstructions are visually in good agreement with the manual segmentations. The dynamic changes, most rapid at time steps 10 to 12 (Fig. 5b), are well captured while slight misalignment artefacts in the manual segmentations (bright curve at the right side of Fig. 5c) are successfully suppressed in the rSIRT-PWC-DIFF reconstructions. Moreover, Fig. 5c reveals small water structures disappearing and re-appearing in the manual segmentations, while the structures remain stable in the rSIRT-PWC-DIFF reconstructions. This qualitative performance is supported by the error metrics in Table 2: high sensitivity of 97% is reached across the time sequence while the specificity is 99% and the dice coefficient 94%. Even slightly For the PEFC_2 sample, the rSIRT-PWC-DIFF reconstructions are presented in Fig. 6a (middle column) for time frames 1, 5 and 10, cropped to the area containing the dynamic changes. For comparison, the manual segmentations (left column) and standard Gridrec reconstructions (right column) are presented. For the manual segmentations of PEFC_2 no additional post-processing step to exclude fibers has been considered and so, this step has been omitted for the rSIRT-PWC-DIFF reconstructions as well. The manual segmentations were found to suffer from misalignment artefacts (Fig. 6a) due to sub-optimal alignment of the dry and dynamic reconstructions prior the subtraction step. To avoid evaluating such artefacts in the iterative reconstructions, the reconstruction quality was measured only in chosen sub-regions (Fig. 6a) for which the segmentations could be visually confirmed successful based on the original Gridrec reconstructions. The sensitivity, specificity and dice coefficient measures ("Automatic stopping criterion" section) are reported in Table 2. Figure 6 shows high agreement between the manual segmentations and the rSIRT-PWC-DIFF reconstructions. While the dynamic process was found to proceed at moderate speed, the slight changes through the time sequence (Fig. 6b,c) are captured well while the stable water features are correctly reconstructed and remain static through the time sequence. Moreover, the misalignment artefacts (Fig. 6a) are found to be suppressed in the rSIRT-PWC-DIFF reconstructions (Fig. 6d). The qualitative comparison is supported by the quality metrics considered at chosen sub-regions (Fig. 6a) and presented in Table 2: the sensitivity of the reconstructions across the time sequence was 92%, with high specificity of 98% and dice coefficient of 93%.

Discussion
We have introduced the rSIRT-PWC-DIFF algorithm, designed to reconstruct and segment, in an unsupervised manner, dynamic, low SNR interior tomography datasets consisting of static and dynamic structures without prior knowledge of the sample composition and its inner features. The proposed algorithm exploits difference sinograms between static and dynamic tomograms to directly extract dynamic features automatically without need for prior reconstruction and segmentation. Time regularization in terms of a piecewise constant function fitting (PWC) approach is applied to 5 × 5 pixel sub-regions with a Gaussian weighting scheme, so to effectively suppress noise within the dynamic reconstructions. Moreover, a stopping criterion was implemented to determine automatically the number of required rSIRT iterations and the starting point of the time-regularization step. Strive for automatization has been guiding this work: prior knowledge and input required by the algorithm has been strongly minimized. In this way, the algorithm can be applied unsupervised to large datasets with hundreds of time steps and to datasets of different samples.
Currently the algorithm works on single slices independently assuming that the sample is stable in the vertical direction, while drifts in the axial plane are automatically accounted for. To generalize the approach further, it is necessary in the future to include sinogram alignment in 3D. In this way any arbitrary drift in the vertical direction, not uncommon in evolving samples, could also be corrected during the alignment process, so to minimize any misalignment artefacts.
The current sinogram alignment is sufficient for aligning scans for which the relative position of the sample with respect to the rotation axis has remained stable. However, in cases where the sample position relative to the rotation center has changed significantly between scans, the current sample alignment will result in reconstruction artifacts arising from disagreeing sinogram shapes. In such cases we recommend to perform the current alignment procedure on a selection of lines (batches), the selection size depending on the amplitude of the relative sample position change, so to achieve an optimized alignment. Alternatively, an additional pre-processing step can be performed, in which the scans are first reconstructed (e.g. with FBP), coarsely aligned and after this, forward projected for final alignment followed by subtraction and the proposed iterative reconstruction. By default, the approach aligns sinograms assuming stable relative sample position for faster computation.
Currently the reconstruction, exploiting a GPU for the rSIRT iterations, is performed in approximately 2.5 min for one slice with 30 time steps (PEFC_1 dataset). To expand the possible applications for large-scale datasets, improvements in computational efficiency are foreseen. For accelerated reconstruction, an interesting possibility is to exploit approximated algebraic filters [40][41][42] to obtain reconstructions with comparable image quality to SIRT with the computational effort of FBP. Moreover, further improvement of the efficiency of the timeregularization step is envisaged through neural networks which have already demonstrated impressive results in improving the tomographic reconstruction quality when used as a post-processing tool 43,44 .
For datasets with highly limited SNR and complicated sample structures, common in time-resolved X-ray tomographic microscopy investigations, the typical data processing pipeline includes data reconstruction, registration, filtering and segmentation, each step requiring manual parameter tuning and optimization. Once Table 2. Sensitivity, specificity and dice coefficient between the rSIRT-PWC-DIFF reconstructions and manual segmentations for PEFC_1 and PEFC_2 datasets. For PEFC_1 each metric was computed over the whole image plane and averaged over the time sequence. For PEFC_2 the metrics were averaged over the time sequence at chosen sub-regions due to remaining misalignment artefacts in the segmented images (Fig. 6a). www.nature.com/scientificreports/ the optimal parameters have been identified, the pipeline can be typically applied to a time sequence with only minor parameter tweaking and manual intervention. Unfortunately, the parameters do though not generalize well, requiring a new parameter evaluation when experimental conditions or samples are modified. These major manual requirements considerably extend the total data post-processing time, limiting possible dynamic analysis realistically to a few samples and time steps. The proposed algorithm protocol offers a significant reduction in manual labor while delivering robust and high-quality reconstructions. This was demonstrated by automatically extracting liquid water dynamics from sub-second fuel cell X-ray tomographic microscopy datasets. The post-processing pipeline used so far to extract water dynamics from reconstructed images relies on subtracting a dry scan from dynamic scans. Manual tuning and supervision are unavoidable to optimally adjust the scan alignment and filtering parameters for these SNR-limited data to obtain segmented water. The parameters are dataset specific and need to be adjusted for each experiment individually. Despite careful alignment, remaining misalignment artefacts are inevitable, leading to the need for careful post-processing to avoid such artefacts appearing as falsely detected water in the final segmented images. For good quality data of a single cell, considering 100 full volume time steps, the complete process from raw data to segmented water requires in the best case 2 weeks of work, significantly dominated by the necessary manual interaction and not by the computation distributed on a maximum of 15 CPU nodes. In case of challenging image quality or completely new cell types and materials, additional manual tuning and exploration of the parameter space might be necessary, extending the processing time from 2 weeks up to 1 month. Thanks to the generalization of the stopping criterion, the proposed algorithm enables to perform reconstruction and segmentation of comparable data fully automatically www.nature.com/scientificreports/ in 1 week on a single GPU node. In addition, thanks to the subtraction step performed prior to back-projection and time-regularization, misalignment artefacts are optimally suppressed. As the algorithm is completely independent from manual intervention, it can be easily scaled up to efficiently process even larger data amounts by utilizing additional GPUs and parallelization. The reduced total post-processing time as well as the fully automatized protocol will open up new possibilities in dynamic X-ray investigations (e.g. exploration of a considerably larger experimental parameter space), for which the manual effort required by conventional pipelines would represent an unsurmountable obstacle. Though this algorithm has been designed to meet especially the needs and requirements of fuel cell reconstruction, it is not limited solely to fuel cell imaging applications but can be applied to any dynamic dataset for which either a static reference scan or time frame is available. Thanks to the high-quality reconstructions of limited SNR datasets delivered by the proposed rSIRT-PWC-DIFF in a highly automatized manner, it will be possible to push the temporal resolution of dynamic tomography experiments even further to the sub-second scan time region and beyond and efficiently deal with the related large data volumes.

Materials
Simulations. The fuel cell phantom was designed to mimic the static structures of an operando X-ray tomographic microscopy fuel cell setup 45 : a flow field with two gas channels was placed both on the anode and cathode sides of the cell together with a gas diffusion layer (GDL) consisting of thin, randomly distributed carbon fibers. A bright layer was added between the anode and cathode GDLs to mimic the polymer electrolyte membrane coated with Pt-based catalyst. A total of 100 phantoms were simulated with randomized GDL fiber structure on a 400 × 400 pixel grid, each having one larger droplet randomly positioned in one of the channels. In addition, 2 to 20 smaller droplets were randomly chosen and positioned within the GDL. All simulated droplets were evolving in time.
Poisson distributed noise, assuming an incoming beam intensity of 5 × 10 3 photons per pixel, was applied to the standard simulated data to simulate strongly noisy images and to avoid algorithm overfitting. The estimated SNR of the simulated FBP reconstructions was a factor of 3.7 lower than the SNR of the phase retrieved Gridrec reconstructions of the real data ("Real data" section), leading to a performance assessment for strongly noisy imaging conditions.
The membrane swelling was simulated by shifting all dry and dynamic structures in the wet simulated scans in vertical direction while the dry scan phantom remained unchanged. To ensure that most misclassifications are related to misalignment, a reduced Poisson distributed noise, compared to the previous experiment was applied. The chosen incoming beam intensity of 30 × 10 3 photons per pixel lead to a comparable SNR in a simulated FBP reconstruction as measured in phase retrieved Gridrec reconstruction of real data ("Real data" section).
Real data. Two dynamic Polymer Electrolyte Fuel Cell (PEFC) synchrotron tomographic datasets were collected at the beamline for TOmographic Microscopy and Coherent rAdiology experimenTs (TOMCAT) at the Swiss Light Source (SLS) at the Paul Scherrer Institut, Switzerland 46 . Both PEFCs had a diameter of approximately 5 mm and were mounted on the rotation stage. For both experiments, a high-numerical-aperture macroscope 47 providing a 4 × magnification was coupled with a 150 µm thick LuAG:Ce scintillator (Crytur, Turnov, Czech Republic) to convert the X-rays into visible light. The in-house developed GigaFRoST detector 48 was exploited to enable continuous data acquisition at the maximum rate of nearly 8 GB s −1 . The radiograms were acquired using filtered polychromatic radiation with a peak energy of approximately 20 keV. To match the required time resolution, the horizontal field of view was reduced to approximately 4 mm, leading to interior tomography datasets.
The first cell (PEFC_1) (flow field plates made of BMA5, SGL Technologies; membrane electrode assembly: Gore Primea A510.1/M815.15/C510.4 with a 15 µm thick membrane and anode/cathode Pt loadings of 0.1/0.4 mg/cm 2 ; GDL: Toray TGP-H-060, Toray Industries Inc.) was connected from above with two long tubes which provided the cell with continuous gas flows, enabling continuous cell rotation during operation. To avoid blurring artefacts caused by continuous sample rotation during data acquisition, 300 evenly distributed projections of the PEFC_1 sample were collected between 0 and 180 degrees. The exposure time was set to 0.3 ms for each projection image, leading to a total scan time of 0.1 s per tomogram. The collected images were 1440 × 1100 pixels with a pixel size of 2.75 µm. A single dry scan was also collected with the same experimental parameters. Five sets of tomographic datasets were acquired, each consisting of 60 continuous time steps. The data for the first 3 sets (3 × 60 scans) are available in TomoBank 49 .
For the second cell (PEFC_2) (flow field plates made of BMA5, SGL Technologies; membrane electrode assembly: Gore Primea A510.1/M815.15/C510.4 with 15 µm thick membrane and anode/cathode Pt loadings of 0.1/0.4 mg/cm 2 ; GDL: SGL 28BC, SGL Carbon), a heated fluid rotary union was developed to exploit continuous gas flow and sample rotation at elevated operation temperature during the experiment. A total of 120 tomograms were collected continuously, each time frame consisting of 299 homogenously distributed projections. The scan time was set to 0.1 s per tomogram, leading to an exposure time of 0.33 ms for each projection. The collected images were 1440 × 1000 pixels with a pixel size of 2.75 µm. A single dry scan was collected with the same experimental parameters.
Prior to reconstruction, all projection images were dark and flat-field corrected, followed by phase retrieval 50 using the processing pipeline available at TOMCAT 51 . Paganin phase retrieval was applied together with a deconvolution step to enhance the contrast between different materials while aiming to maintain the spatial resolution, typically compromised by the phase retrieval procedure 52 .

Data availability
The PEFC_1 dataset is available in TomoBank (https ://tomob ank.readt hedoc s.io/en/lates t/sourc e/data/docs. data.dynam ic.html#foam-data). The PEFC_2 dataset and the code are available upon request.