A universal and effective variational method for destriping: application to light-sheet microscopy, FIB-SEM and remote sensing images

Stripe artifacts are a common problem for various imaging techniques such as light-sheet fluorescence microscopy (LSFM), electron microscopy and remote sensing. These artifacts are characterized by their elongated shapes, compromise image quality and impede further analysis. We present a robust, openly available [1] and universally applicable variational method for stripe removal. We demonstrate that our approach generalizes well across different data sets and give recommendations for its parametrization. The method is compared against previously published methods on both synthetic LSFM data and real images from LSFM, FIB-SEM and remote sensing. Visual inspection and quantitative metrics demonstrate that our method outperforms existing solutions with a better removal of artifacts. The method’s flexibility in handling variations in stripe orientation and thickness ensures its broad applicability across diverse imaging scenarios.


Introduction
Stripe artifacts are a common problem in several image acquisition settings including light-sheet (fluorescence) microscopy (LSFM), focused ion beam scanning electron microscopy (FIB-SEM) and remote sensing.In general, these artifacts are characterised by highly elongated shapes of low width which point in a common direction.Their removal is required not only to improve visual quality but also to enable further analysis and image processing on the data.The examples shown in Fig. 1 highlight the diversity of image structures and stripes encountered in practice.While stripes are periodic, thin, long and of low visibility in (c), much more severe artifacts of varying width, smaller length and higher intensity can be seen in (a).There exists a large body of literature dealing with stripe removal including the recently published article by Ricci and colleagues [2] which provides a comprehensive overview of previous research on stripe removal in LSFM.Possible approaches are prevention [3][4][5][6][7], post-processing [8][9][10][11][12] and hybrid solutions [13,14].Preventing the formation of artifacts would be preferable but typically requires specialized hardware equipment or can limit the performance [15][16][17].In contrast, post-processing is cheap and requires only computational power.Furthermore, corrupted images may already exist and specimens cannot be re-imaged, e.g., due to degradation and aging.For these reasons it is desirable to have powerful post-processing tools for stripe removal available.
In this article, we propose and explain a general solution for removing stripe artifacts.The results in this article quantitatively compare our method to previous approaches across images from LSFM, FIB-SEM and remote sensing.Our solution generalizes exceptionally well on all data tested.It is highly flexible thanks to intuitively adjustable parameters to accommodate for different appearances of structures and stripes.The source code is shared on GitHub [1] together with comprehensive guidance for application without further knowledge.The method which belongs to the class of variational methods is compared against selected state-of-the-art solutions from the same category and Fourier filtering methods.The performance is quantitatively evaluated with the help of synthetic LSFM data which is obtained by physically correct simulation of light transport with the Python package biobeam [18].We include simulation data since it provides ground truth images which are unavailable for experimentally acquired images.This enables an objective assessment of performance using common quality metrics such as peak signal-to-noise ratio or the multi-scale structural similarity index measure (MS-SSIM).The evaluation is supplemented by a measure of stripe corruptions called "curtaining" proposed by Roldán [19] and the difference in line profiles.In summary, we propose a solution to destriping, present results for different imaging methods and using several metrics compare with previously published solutions.

Methods
The majority of methods proposed for stripe removal belong to the two categories Fourier filtering and variational methods.We will concentrate on these categories, as they generalize well to several imaging methods, variations in image structures and stripe appearance.Other approaches include average filtering [20,21], histogram matching [22,23], spline interpolation [24] and recently neural networks [25][26][27][28].However, these are usually tailored for a specific appearance of images and stripes such that they are harder to transfer to other scenarios.

Variational Methods
Variational methods perform tasks such as denoising, segmentation, active contours [29] or destriping [30] by minimizing a task-specific convex objective function.It is constructed by penalizing unwanted features of the image and stripes.Therefore, a minimizer of the functional should have desirable properties, e.g., should be free of stripe artifacts.Under mild assertions on the function, a unique minimizer exists and can be approximated using well-studied optimization algorithms, see [31].
Several variational methods for removing stripes were proposed in the past.The contributions [10,32] represent stripes using elementary stripe patterns.Other propositions assume sparsity and a low-rank assumption on the stripe image [33][34][35][36].Alternatively, [37] includes additional information from the Fourier domain.In the following, we introduce, motivate and explain our objective function which builds on previous work by Fitschen et al. [30].A similar objective function to ours was also explored by Liu et al. [38] for remote sensing images.In contrast to previous implementations our objective function aims to effectively address stripes and retain image structures while providing a reasonable intensity profile.

Optimization Problem
We assume that we observe a corrupted image  0 ∈ [0, 1]  which can be decomposed as  0 =  + , where  ∈ [0, 1]  is the clean image and  ∈ [−1, 1]  are the stripes.The image sizes are  =   ×   ×   .Let ∇  , ∇  and ∇  denote the directional difference operators for the three coordinate directions and assume stripes in -direction.The objective function of variational methods often consists of a data term  () and noise term  () which independently penalize properties of the clean image and stripes.We propose and ∞ otherwise and (2) is known as total variation [39] and encodes that clean images contain few strong edges.The choice of ∇   1 is based on the assumption that stripes have large differences orthogonal and small differences parallel to the stripe direction.Additionally, ∥∥ 1 promotes sparsity in the stripe image and reflects that typically only a small part of the image is directly affected by artifacts.The indicator function  [0,1]  () ensures the output to live on the same range of values as the input.It prevents the need for brightness adjustments after processing and ensures non-negativity similar to natural intensities.The parameters  1 and  2 control the stripe removal, regularization and smoothness.They offer adjustability to make our method applicable to different settings ranging from severe corruptions to weak obscurities.For a discussion of the choice of  1 and  2 we refer to Supplement 1 Section 2. A basic intuition is given by: •  1 adjusts the strength of stripe removal.Larger values result in stronger removal but may affect stripe-like structures and yield smoothing.
•  2 affects the precaution towards retaining image structures.Increasing  2 enforces sparsity of the stripe image.Hence, larger values impair the removal, especially for non-stripe elements.
• The ratio of  1 and  2 is essential.Scaling both by an equal factor changes the amount of stripe removal while retaining the effect on image structures for the most part.
Optimization of ( 1) is done with the primal-dual gradient hybrid method with extrapolation of the dual variable (PDHGMp) [31].The sequence generated by this iterative algorithm converges to the solution [31,40].To generalize (1) we propose to replace ∇  with ∇ *  =   ∇  where   ∈ [0, 1].Setting   = 0 reflects a 2D case while 0 <   < 1 implies that information along the -direction is less reliable and coherent, e.g., due to lower resolution.For oblique stripes in direction  we replace ∇   by ∇   as suggested in [38].Furthermore, simultaneous multi-directional destriping is possible by including a sum of penalization terms ∇    1 with different stripe directions   .

Fourier Filtering
The Fourier transform is a mathematical decomposition of an image into constituent frequencies.If the fundamental assumptions on stripes to be parallel and thin is fulfilled, stripe information is mostly encoded in frequency coefficients in a small band orthogonal to the stripe direction.In most cases, natural structures are dominated by low frequencies and concentrate near central coefficients.This decisive difference is exploited by Fourier filtering approaches by using masked filters [41][42][43], a decision-based algorithm [11] or a prior structural decomposition [8,9] to restrict tampering to stripe related coefficients.The latter approach reduces the influence on image structures the most.While wavelet-Fourier filtering [8] is commonly referenced, the approach by Liang et al. [9] using the non-subsampled contourlet transform (NSCT) [44] improves stripe and structure separation resulting in less interference during removal.We compare our method against a slightly enhanced version of this approach abbreviated by MDSR + .The method is explained in Supplement 1 Section 3.

Results
In this section, we evaluate and compare our method against state-of-the-art solutions from the category of Fourier filtering and variational methods.In particular, we consider the MDSR + , see Section 2.2, and the variational stationary noise remover (VSNR) [10].For the MDSR + we use the NSCT with 8 directions and the maximum possible depth which depends on the image size and was either 4 or 5.The VSNR was used with three differently sized stripe patterns created from Gabor filters to be generally applicable.Our method (1) was used with   = 0 as MDSR + and VSNR were only available for 2D processing.The number of optimization steps taken for our method and VSNR were chosen to be 25000 to obtain a sufficient approximation of the optimum.The remaining parameters of the three stated methods were optimized via grid search, numerical and visual assessment.For the numerical assessment on synthetic data with ground truth we utilize the common performance metrics peak signal-to-noise ratio (PSNR) and the multi-scale structural similarity index measure (MS-SSIM) [45].Furthermore, we consider the difference in line profiles of a selected area to the ground truth using the 2-norm.Additionally, we use the curtaining metric proposed by Roldán [19].In contrast to the other metrics, it exclusively measures stripe corruptions based on the ideas of Fourier filtering methods, see Section 2.2.We report the absolute difference between the curtaining metrics of the processed image and the ground truth.This choice adjusts for naturally occurring stripe-like patterns in the image structures and treats over-and under-performance alike.

Synthetic Data
Fig. 2 shows results on a synthetic embryo-like structure inspired by [18] and obtained by physically accurate modeling of light propagation in LSFM imaging, see Supplement 1 Section 1 and [46].Major differences between the methods can be spotted.While MDSR + reduces the stripes significantly, larger artifacts and thin oblique stripes remain visible in the outcome.The latter become only visible at closer inspection.In comparison, VSNR is able to reduce wider artifacts better with only weak remnants.However, fine oblique artifacts are entirely unaffected by the removal.Only our method removes all stripes including the trails with only faint artifacts remaining in particular at the right edge of the body.The line plot supports our findings.MDSR + shrinks the profile but is unable to adjust wider artefacts correctly.On the other hand, the line profile of VSNR has less extremes but shows small perturbations reflecting the struggle with fine artifacts.Our method is the only one which has both a significant reduction in stripes and a flat profile.The numerical values presented in Table 1 further confirm prior observations with our method outperforming the others in all metrics.

Real Data
Fig. 3 shows results on real data from LSFM, FIB-SEM and remote sensing.We observe similar behaviour to the synthetic data with VSNR leaving very thin remnants and MDSR + having insufficient reduction.The LSFM image highlights the capabilities of variational methods in general.While MDSR + only reduces stripes, both VSNR and our method show near perfect on par performance.However, the FIB-SEM and remote sensing images reveal short-comings of  VSNR which leaves very thin stripes in the result.MDSR + even introduces severe image artifacts in areas of high contrast in the remote sensing image despite performing well in low contrast regions.On the other hand, our method delivers visually perfect results for both FIB-SEM and remote sensing in the sense that stripes are entirely removed and image structures remain intact.

Discussion
The results shown in Section 3 demonstrate that our method produces consistently good results in terms of stripe removal quality and retention of image details while offering intuitive adjustability through its two weighting parameters  1 and  2 , see Supplement 1 Section 2. The consistency across all shown settings can be attributed to the non-restrictive formulation of its objective function with only rough assumptions on properties of the image and stripes.Therefore, deviations in stripe direction and thickness remain adequately penalized such that a general application for stripe removal is provided.However, limitations exist.In particular, the alignment of image structures with the stripe direction should be avoided at all cost since both will become affected during removal, see Supplement 1 Section 4.
The frequently occurring insufficient removal by MDSR + is due to a large overlap in area of stripe and image related frequency coefficients when structures and stripes live on similar scales or stripe directions deviate from the assumption.Depending on the choice of parameters, either insufficient stripe removal or severe modification of image structures are the result.We deem the prior to be preferential since the latter renders the outcome unusable for any analysis or visualisation which is why we only showed these results.The introduction of artifacts in the remote sensing image is surprising as the setting reflects an ideal scenario for this approach with thin, long and periodic stripes.However, the high contrast in some areas produces unforeseen problems by filling separating spaces with smooth stripe-like artifacts.
In contrast, VSNR performs overall well but has problems removing thin or oblique stripes.This behaviour arises since the considered stripe patterns do not sufficiently reflect these particular instances of stripes.While we see potential of VSNR to improve using more image specific patterns we argue that optimizing patterns and corresponding parameters requires considerable expertise, is time consuming and impractical for a general use.The results in Fig. 3 (top) showcase that using optimized patterns, good results can be achieved.However, results are expected to be worse with less optimized patterns.
The usage of synthetic LSFM data with physically accurate modeling of light propagation is a novel supplement to real data for comparing quality of stripe removal.It allows for the calculation of a ground truth and enables the use of quality metrics such as PSNR.This results in a more accurate assessment of quality compared to the classical approach of superimposing a clean image with independent stripes [9] or relying only on visual assessment.
Roldáns curtaining metric [19] complements previously established metrics well for assessing quality on the synthetic data.Since it is specifically designed to measure the amount of stripe corruptions in an image it is more robust against other difference between the image and its ground truth, e.g., different levels of brightness.One could argue that the metric should be applied directly instead of relying on the absolute difference with the ground truth.However, overly aggressive destriping with any of the presented methods almost surely gives an optimal curtaining value despite performing poorly in every other aspect.By adjusting for the naturally occurring imperfections in the structures using the ground truth, performance is assessed in line with visual inspection and the other metrics.
The results shown in this article were processed as 2D images to simplify visualization and compare with methods provided only in 2D.However, our method is capable of application on 3D images too, see Supplement 1 Section 4. Processing in 3D improves consistency across slices and prevents possible errors in image structures that result from ignoring relations in the stack direction.

Conclusion
In this work, we have demonstrated the potential of variational methods as a general solution for removing stripe artifacts from a variety of image structures and different imaging methods using our objective function.We provide insight into the intuitive choice of parameters to optimize destriping and preservation of image structures.We demonstrate that it is highly effective in removing stripes and exceeds the performance of previously published solutions.The generality of the approach is backed by using images from three different imaging methods: LSFM, FIB-SEM and remote sensing.A quantitative comparison on synthetic LSFM data using several quality metrics helped to better understand the capabilities and limitations of the stripe removal methods.It confirmed that our solution provides best destriping while retaining image structures which matches the visual assessment of real images.

Supplemental document
See Supplement 1 for supporting content.

Fig. 1 .
Fig. 1.Stripe corruptions in cultivated mouse intestine cells imaged by LSFM (a), tin bronze cast alloy imaged by FIB-SEM (b) and Terra MODIS data from remote sensing (c).

Fig. 2 .
Fig. 2. Destriping results on synthetic LSFM and corresponding line profiles for the marked area.

Fig. 3 .
Fig. 3. Stripe removal results on real data acquired by LSFM (top), FIB-SEM (middle) and remote sensing (bottom).The red arrows highlight areas of particular interest.

Table 1 .
Quantitative evaluation of results on synthetic data from Fig.2.Arrows indicate better performance.Percentage improvements from the input are given in green.