Fourier Domain Anomaly Detection and Spectral Fusion for Stripe Noise Removal of TIR Imagery

Abstract

Stripe noise is a common and unwelcome noise pattern in various thermal infrared (TIR) image data including conventional TIR images and remote sensing TIR spectral images. Most existing stripe noise removal (destriping) methods are often difficult to keep a good and robust efficacy in dealing with the real-life complex noise cases. In this paper, based on the intrinsic spectral properties of TIR images and stripe noise, we propose a novel two-stage transform domain destriping method called Fourier domain anomaly detection and spectral fusion (ADSF). Considering the principal frequencies polluted by stripe noise as outliers in the statistical spectrum of TIR images, our naive idea is first to detect the potential anomalies and then correct them effectively in the Fourier domain to reconstruct a desired destriping result. More specifically, anomaly detection for stripe frequencies is achieved through a regional comparison between the original spectrum and the expected spectrum that statistically follows a generalized Laplacian regression model, and then an anomaly weight map is generated accordingly. In the correction stage, we propose a guidance-image-based spectrum fusion strategy, which integrates the original spectrum and the spectrum of a guidance image via the anomaly weight map. The final reconstruction result not only has no stripe noise but also maintains image structures and details well. Extensive real experiments are performed on conventional TIR images and remote sensing spectral images, respectively. The qualitative and quantitative assessment results demonstrate the superior effectiveness and strong robustness of the proposed method.

1. Introduction

Thermal infrared (TIR) imagery has been an indispensable data resource in the remote sensing field because of the peculiarity of this spectrum [1]. Typically, satellite-based TIR images are widely used to investigate and monitor earth resources and environment by characterizing surface temperature and emissivity [2,3,4]. Conventional TIR cameras are increasingly mounted onto some low-cost aerial platforms like drones and airships to collect images for detection, identification and analysis missions [5,6,7]. In these available TIR data, however, raw images are frequently found to suffer from a distinctive noise pattern that manifests itself as vertical or horizontal stripes. This phenomenon not only results in the visual degradation of images but also affects the precision of the downstream processing.

1.1. Two Striped TIR Image Cases

In real TIR image data, there are two common cases with stripe noise: one is conventional TIR images that mostly come from uncooled infrared focal plane array (IRFPA)-based thermal imagers [8,9,10,11] and the other is remote sensing TIR spectral images that are captured by spaceborne/airborne imaging spectrometers [12,13,14,15]. For conventional TIR images, the presence of stripe noise is attributed to the nonuniform response of IRFPA detectors, which are typically designed with a column-shared readout circuit architecture [8,11]. From the sensor perspective, in practice, to eliminate the response errors between detector elements, the calibration-based nonuniformity correction (CNUC) technology could be applied in TIR imaging systems. But this does not necessarily ensure the complete disappearance of stripe noise in the image, such as when the CNUC operation is not performed well or its efficiency declines in the long term usage due to the variation of the response characteristics of detectors. As for remote sensing spectral images, the occurrence of the striping effect has a direct relationship with the work mechanism of imaging spectrometers, which record the one dimensional cross-track measurement by whiskbroom or pushbroom scanning while the observation in the other spatial dimension depends on the platform in-track motion [14]. Because of unavoidable interferences and errors during the scan imaging, stripes are generated in the cross-track direction of spectral images, although a series of radiometric calibration procedures are also performed beforehand [15].

Given the above two cases, the destriping task is essentially identical and necessary. So in this paper we focus on these striped TIR images and strive to provide an universal and efficient solution to improve their quality for further thermal applications.

1.2. Related Work

Over the past decades, plenty of stripe noise removal methods have been proposed, which fall into two broad categories: processing in the spatial domain and in the transform domain.Both processing technologies are developed by excavating and utilizing the inherent differences between stripe noise and image. And more attention is paid to the single image based destriping problem.

Spatial methods for destriping can be roughly divided into 1D-filtering-based methods [11,16,17], statistics-based methods [12,13,14,18,19,20], optimization-based methods [21,22,23,24,25,26,27,28,29,30], and deep-learning-based methods [31,32,33,34]. Currently, 1D-filtering-based methods are mainly performed on conventional TIR images, and their basic idea is to utilize 1D edge-preserving filters to progressively separate stripe noise from the contaminated image in horizontal and vertical directions. For example, Cao et al. successively used 1D row and column guided filters to achieve the estimation of stripe noise [16]. However, this method may produce some spurious artifacts near strong vertical edges, as has been demonstrated in [17]. Statistics-based methods, as one kind of classical correction methods, build on the similarity assumption of the data that global or local subscenes (or striping lines) should share the same statistics such as the histogram distribution [13,19], and the mean and standard deviation [12,14]. Such prerequisites are simple, efficient but idealistic, often resulting in an inadequate destriping effect in the complex real-world cases. For this reason, the combination of the statistics-based methods and other processing methods is becoming a practicable solution to balance the effectiveness and efficiency in real applications [15]. In recent years, optimization-based methods have been studied popularly owing to powerful model representation, excellent destriping result and flexible regularization strategy. In these optimization frameworks, the desired clean image is computed by minimizing an energy function with different regularization terms, which mathematically represent the prior knowledge of stripe noise and image. Bouali et al. considered the directional feature of stripes and proposed an unidirectional total variation (UTV) model to remove stripe noise while preserving image details [21]. Subsequently, many UTV variants are developed to further enhance the model’s adaptivity [22,23,24,25]. Other prior knowledge of stripes including low rank [26,27] and sparsity [28,29,30] are primely explored and applied in the relevant models as well. From a pragmatic perspective, these sophisticated optimization approaches tend to be multi-parametric and time-consuming, which makes them have lower utility. More recently, some deep-learning-based methods have emerged and created impressive destriping results, benefiting from the neural networks’ formidable function. In [22], Xiao et al. designed a nine-layered convolutional neural networks (CNNs) architecture for removing stripe noise from a single meteorological satellite infrared cloud image. He et al. added a polynomial simulation module of stripe noise into their CNNs training to improve networks’ discrimination ability of the noise [23]. So for these deep-learning-based methods, there may be two aspects worth considering to upgrade their performance; namely, enhancing the discrimination of stripe noise and the comprehension of image content [35]. The concrete challenges could include a realistic simulation model of stripe noise, a particular objective function (with some prior knowledge about the noise and image), an elaborate network architecture and so on.

As an alternative for destriping, transform domain methods that are chiefly designed and realized in the Fourier domain and wavelet domain have also received considerable attention and research [35,36,37,38,39,40,41,42,43,44]. One significant starting point of these methods is that stripe noise has a concentrated energy distribution in both domains due to its directionality and global similarity. Perhaps the most representative one among the transform domain methods is the combined method of wavelet-Fourier filtering [36,37]. It takes full advantage of the noise characteristics in the wavelet and Fourier domains to jointly filter out stripes while maintaining the detail information of image. Before this method, some simple filtering methods in the Fourier domain and wavelet domain had been proposed and achieved initial results [38,39,40,41]. Recently, researchers have introduced some advanced spatial processing techniques into transform domains and thus developed a number of state-of-the-art transform-based methods for destriping. In [42], a multi-scale operation of guided filtering is performed on the noisy wavelet coefficients to adaptively estimate stripe noise from vertical high-frequency details. One preliminary study of ours on destriping is to correct the dominant Fourier coefficients that are contaminated by stripe noise with a reference spectrum [43]. Moreover, the idea of deep learning is being integrated into the wavelet domain to produce possibly better destriping results [35,44].

1.3. Motivation

In this paper, we focus on the destriping task in the Fourier domain for TIR images. Early Fourier domain methods remove stripe noise by using simple 1D filters [36,37], and some newly related methods design adaptive 2D filters to accomplish the goal [45,46]. However, as many literatures pointed out, on the one hand, these naive filtering approaches have limited ability to handle stripes when faced with manifold real cases and most of them are in fact only available for periodic stripes; on the other hand, the methods may cause the so-called ringing artifact owing to the high discontinuity of image intensity and the blur of useful structure details that have the same frequencies as stripes. These drawbacks restrict the development of the Fourier domain destriping technique to some degree and make it fail to show a good competitiveness in processing real striped images. Therefore, this paper attempts to make some contributions and provide a new line of thought for Fourier domain based destriping methods.

To be specific, we propose a Fourier domain anomaly detection and spectral fusion method to remove stripe noise of TIR images. In this two-step framework, the abnormal frequencies that are likely to be contaminated by stripe noise are detected intentionally so that an anomaly weight map is generated to represent the information; then a fusion strategy between the original spectrum and the spectrum of a stripe-free guidance image is adopted to obtain the new corrected spectrum. Finally, the resultant image has no stripe noise while keeping the original structures and details well. Experimental results on conventional TIR images and remote sensing TIR spectral images fully demonstrate the superior destriping performance of the proposed method by comparing with other state-of-the-art destriping methods in both fields. The main idea and contributions of the proposed method are summarized as follows:

• The respective spectral characteristics of TIR images and stripe noise are explored systematically. On this basis, a Fourier domain anomaly detection method is proposed to locate the abnormal frequencies that are likely to be contaminated by stripe noise. As a result, an anomaly weight map is obtained.
• With the anomaly weight map, a fusion strategy between the original spectrum and the spectrum of a stripe-free guidance image is adopted to generate the new Fourier spectrum and then to reconstruct the destriped image. In the implementation, the guidance image is estimated via an adequate de-texture filter.
• Extensive experiments on real striped TIR data including conventional TIR images and remote sensing spectral images are performed. The results demonstrate that the proposed method has better destriping ability and stronger robustness, compared with the state-of-the-art destriping methods in the two fields.

The remainder of this paper is organized as follows. In Section 2, the spectral characteristics of TIR images and stripe noise are investigated. In Section 3, the proposed Fourier domain anomaly detection and spectral fusion method is described. In Section 4, we present the experimental results to validate the algorithm performance. Finally, conclusions are drawn in Section 5.

2. Spectral Analysis

2.1. Brief Reminder and Notations

In the Fourier theory of 2D image, the discrete Fourier transform (DFT) and the inverse discrete Fourier transform (IDFT) are defined as

$$\begin{array}{c}\hfill \mathcal{F}(u,v)=\sum _{x=0}^{X-1}\sum _{y=0}^{Y-1}I(x,y)e^{-j2\pi (ux/X+vy/Y)}\end{array}$$

(1)

$$\begin{array}{c}\hfill I(x,y)=\frac{1}{XY}\sum _{u=0}^{X-1}\sum _{v=0}^{Y-1}\mathcal{F}(u,v)e^{j2\pi (ux/X+uv/Y)}\end{array}$$

(2)

where $I(x,y)$ is a $X\times Y$ sized image and $\mathcal{F}(u,v)$ is its Fourier spectrum with the same size; $(x,y)$ and $(u,v)$ are spatial and spectral coordinates and have integer values; j is the imaginary unit. In practice, $\mathcal{F}(u,v)$ is often shifted so that the origin $(u,v)=(0,0)$ is at the center of the whole spectrum. Using $(u,v)\in [-X/2,X/2-1]\times [-Y/2,Y/2-1]$, we denote the corresponding normalized frequencies by $f_u=u/X$ and $f_v=v/Y$ (units are cycles per pixel) in each direction and thus $(f_u,f_v)\in [-0.5,0.5)\times [-0.5,0.5)$. $\mathcal{P}={\left|\mathcal{F}\right|}^2$ is the power spectrum of I.

In this work, since the representation of $\mathcal{P}$ with the spatial frequency is involved, we also use the polar coordinate $(f,\theta )$, where $f=\sqrt{f_u^2+f_v^2}$ is the radial distance from $(0,0)$ and $\theta$ is the counter-clockwise angle from the $f_u$ or u-axis. So the spatial frequency is simply denoted by f without the directions. These notations will be unified in later parts.

2.2. Spectral Statistics of TIR Images

As is well known, the power spectrum of an image roughly obeys a power law $1/f^2$, which has been widely applied to natural images [47,48,49]. However, the law may not be perfectly suited to TIR images. Morris et al. in the study of statistics of TIR images found that the average power spectra obtained from their TIR image datasets have increased DC and low-frequency components because of the relative lack of textural details in TIR images [50]. Considering this fact, the authors came up with a generalized Laplace distribution to model the average power of TIR images. It is expressed as

$$\begin{array}{c}\hfill \tilde{\mathcal{P}}\left(f\right)={cexp(-|f/a|}^b)\end{array}$$

(3)

where a, b and c are the model parameters.

To revalidate this finding, we conducted the statistical test on a remote sensing TIR image dataset from [31], which is intentionally chosen to differ from Morris’s conventional TIR datasets. In detail, we first computed the average power spectrum of the dataset and then used the given Laplace model to fit it via the least squares method. Figure 1 shows the survey results, where the dotted curve is the fitted Laplacian distribution and the red part is the average power distribution. It can be seen that the fitted curve matches pretty well with the real distribution. This indicates the effectiveness of the Laplace model expressed in Equation (3). Hence, it is feasible to utilize the generalized Laplace distribution to represent the statistical power spectra of TIR images.

2.3. Spectral Distribution of Stripe Noise

It has been found that structural image components within a fixed direction have a directionally concentrated energy distribution in the Fourier domain and the two directions are orthogonal to each other. The spectral distribution of stripe noise follows this objective law due to its spatial directionality [51]. To further analyse and illustrate it, Figure 2 shows an experimental example where several common types of vertical stripes are given, including the periodic stripes, the non-periodic sparse stripes and the non-periodic non-sparse stripes. We can easily see from the example that the dominant frequencies affected by these vertical stripes are distributed in the horizontal centerline of image spectrums. For the periodic stripes, theoretically, their spectrum should be an ideal comb-like spectrum with only a few fixed frequencies, but as shown in Figure 2b, spectral leakage may occur in reality because of non-ideal periodicity [45]. For the non-periodic (sparse or non-sparse) stripes, their spectral energy evenly spreads to all frequencies [36,52], forming an arresting line as shown in Figure 2c,d. Precisely based on such spectral distribution characteristics, Fourier domain destriping methods can work effectively by detecting and suppressing the dominant noise frequencies. In addition, it’s worth noting that the power spectrum of the Pentagon image is actually a good proof of the above spatial-spectral law, where there are five obvious energy bands within different orientations corresponding to the five-sided structure in the image. This fact also implies a potential problem with the conventional Fourier domain destriping methods, that the image structure which has the same direction as stripe noise could be severely blurred because of the similar spectral distribution. Therefore, as has been stated in the introduction, how to reduce the stripe frequencies while protecting the structurally similar spectral information is an important and realistic issue for Fourier domain based destriping methods.

3. Methodology

In this paper, we focus on the destriping problem of TIR images in the Fourier domain and put forward an effective two-step solution called Fourier domain anomaly detection and spectral fusion (ADSF). Figure 3 is a schematic overview of the proposed method. The first step, anomaly detection, aims at identifying the abnormal Fourier coefficients that are dominantly contaminated by stripe noise. As a result, a weight map is generated to represent the abnormality of Fourier spectrum. The second step is to obtain a stripe-free correction spectrum by manipulating the detected anomalies. For this purpose, a guidance image based spectrum fusion strategy is adopted, which organically integrates the spectral information of both the original image and the guidance image. The final ”fused” image has no stripe noise while holding the original structure and detail information. In the following subsections, we will detail the two parts and some practical considerations in the proposed method.

3.1. Anomaly Detection

Our idea of anomaly detection is based on spectral characteristics of TIR images and stripe noise. As discussed in Section 2, on the one hand, the power spectrum of TIR images in natural statistics can be represented by a generalized Laplace distribution; on the other hand, the frequencies contaminated by stripe noise are relatively anomalous in the whole spectrum because of having high energy and a specific distribution area. For the purpose of this article, that means that the stripe frequencies could be easily detected by comparing the original spectrum with an adaptively estimated spectrum that follows the given Laplace distribution statistically. Driven by such a rough idea, a stripe frequency anomaly detector is reasonably designed. To be concrete, there are two main substeps in the implementation of anomaly detection. The first step is to statistically estimate the adaptive power spectrum for the original TIR image by using the Laplace model. The second step is to detect the stripe-related outliers by making a regional comparison between the real spectrum and the estimated spectrum. Figure 4 illustrates the whole implementation process.

3.1.1. Spectral Estimation

Given a striped TIR image I with size $X\times Y$, we first extract the sliding subimages with size $N\times N$ from I. If the sliding step is $\delta$, the number of subimages is

$$\begin{array}{c}\hfill K=\lceil \frac{X-N}{\delta }\rceil \lceil \frac{Y-N}{\delta }\rceil \end{array}$$

(4)

where $\lceil ·\rceil$ the ceil function. In the implementation, we fix $N=100$ and $\delta =8$.

Statistically, the subimages denoted by ${\{I_i^{\left(s\right)}\}}_{i=1}^K$ constitute a sample set of I. After obtaining the corresponding power spectrums ${\{{\mathcal{P}}_i^{\left(s\right)}\}}_{i=1}^K$ through DFT and a logarithmic transformation, the average power spectrum is easily calculated by

$$\begin{array}{c}\hfill {\overline{\mathcal{P}}}^{\left(s\right)}=\frac{1}{K}\sum _{i=1}^K{\mathcal{P}}_i^{\left(s\right)}\end{array}$$

(5)

Regardless of few exceptional peaks due to stripe noise, the natural distribution for ${\overline{\mathcal{P}}}^{\left(s\right)}$ should behave as a generalized Laplace distribution. Based on this, ${\overline{\mathcal{P}}}^{\left(s\right)}$ is fitted with the function expressed in Equation (3) by least squares-based nonlinear regression to obtain the adaptive model parameters $\overline{a}$, $\overline{b}$ and $\overline{c}$. So the estimated power spectrum is as follows

$$\begin{array}{c}\hfill {\tilde{\mathcal{P}}}^{\left(s\right)}=\overline{c}exp(-|f/\overline{a}{|}^{\overline{b}})\end{array}$$

(6)

We consider ${\tilde{\mathcal{P}}}^{\left(s\right)}$ as an adaptive reference spectrum so that the difference spectrum between ${\overline{\mathcal{P}}}^{\left(s\right)}$ and ${\tilde{\mathcal{P}}}^{\left(s\right)}$ is essentially an indication of all statistical anomalies

$$\begin{array}{c}\hfill {\mathcal{D}}^{\left(s\right)}=max({\overline{\mathcal{P}}}^{\left(s\right)}-{\tilde{\mathcal{P}}}^{\left(s\right)},0)\end{array}$$

(7)

3.1.2. Regional Detection

In order to further detect the anomalies caused by stripes in ${\mathcal{D}}^{\left(s\right)}$, a symmetric tapered slice is designated as the detection region of interest, according to the distribution rule of stripe frequencies. As shown in Figure 5, the size of the detection region $\mathsf{\Omega }$ depends on the taper angle $\alpha$ having a small degree. Such design is based on a practical consideration that high-frequency abnormal components tend to be more than low-frequency abnormal components in most real cases. Without taking the DC component into account, $\mathsf{\Omega }$ can be defined as

$$\begin{array}{c}\hfill \mathsf{\Omega }:|\frac{f_v}{f_u}|\le tan\frac{\alpha }{2}\end{array}$$

(8)

Then, a thresholding operation is performed within $\mathsf{\Omega }$ to locate the stripe-related anomalies in ${\mathcal{D}}^{\left(s\right)}$. The initial anomaly map is obtained by

$$\begin{array}{c}\hfill A^{\left(s\right)}=\left\{\begin{array}{cc}1,\hfill & \mathrm{if} {\mathcal{D}}_{\mathsf{\Omega }}^{\left(s\right)}>\tau \hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

(9)

where $\tau$ stands for the threshold. It is determined via a simple rule $\tau =t{\mu }_f$ where ${\mu }_f$ is the mean of ${\mathcal{D}}^{\left(s\right)}$ at the spatial frequency f and t is set to 3 in the implementation.

The final anomaly weight map W is obtained from $A^{\left(s\right)}$ by two steps. First, we resize $A^{\left(s\right)}$ from $N\times N$ to $X\times Y$ through the bilinear interpolation method. Subsequently, Gaussian filtering is carried out on the resized map to generate an ideal fusion weight map W. Like many conventional fusion methods, here the object of filtering is to repair some possible omissions around the anomalies in the initial weight map. The whole treatment process is formulated as

$$\begin{array}{c}\hfill W=G\otimes (\uparrow A^{\left(s\right)})\end{array}$$

(10)

where ↑ and ⊗ are the interpolation and convolution operations, respectively. G is the Gaussian kernel and its size and standard deviation are set to $5\times 5$ and 2.

Although the dominant frequencies affected by stripe noise are effectively identified by anomaly detection, this does not mean that filtering out the anomalies simply is sufficient to achieve a good destriping result. Because some significant structure information also might lie behind these abnormal frequencies, as we discussed before. A better and more reasonable way is to repair the anomalies and obtain an adaptive corrected Fourier spectrum to reconstruct satisfactory destriping results. Therefore, a novel fusion-based spectrum correction approach is presented in the next section.

3.2. Guidance Image Based Spectrum Fusion

To adaptively correct the abnormality of Fourier spectrum due to stripe noise, we introduce the fusion principle into Fourier domain and propose a guidance image based spectrum fusion scheme. This approach not only achieves good and robust destriping results but also preserves more detail information, especially of some possible structures that are similar to stripes. Figure 6 shows the guidance image based spectrum fusion process. Its main idea is to fuse the Fourier spectrum of the original image with that of a guidance image via the anomaly weight map to obtain the corrected spectrum, thereby reconstructing a stripe-free but detail-preserved result. In the implementation, there are two basic links, generating a desired guidance image and then performing the spectrum fusion.

3.2.1. Guidance Image Generation

For the purpose of spectrum fusion, the guidance image is expected to be a smoothed version that has no stripe noise but keeps the original image information. Here we introduce an existing structure-texture decomposition technology called interval gradient filtering (IGF) [53] to obtain the desired guidance image. On the one hand, the IGF method decomposes images by adopting row-wise and column-wise filtering operations. The 1D filtering based mechanism is very effectual for stripe noise removal, as introduced in Section 1.2. On the other hand, since TIR images usually have more smooth homogeneous components and lower textural details, the IGF method is a perfectly suitable option for the estimation of the guidance image. Figure 7 gives an example of structure-texture decomposition via IGF on the striped image. We can see that the structural image used as a guidance image nicely retains the basic image content, but without any stripe noise. In the following, the fundamentals and implementation procedure of 1D IGF are concisely described.

For pixel k in a row signal R, the interval gradient is defined to measure the difference between the Gaussian weighted averages of the left and right parts in a k-centered window $w_k$.

$$\begin{array}{c}\hfill I{G}_k=g_{\sigma }^r\left(R_k\right)-g_{\sigma }^l\left(R_k\right)\end{array}$$

(11)

where $g_{\sigma }^r$ and $g_{\sigma }^l$, respectively, stand for the right and left 1D Gaussian filtering with the standard deviation $\sigma$ in $w_k$. Compared with the common gradient $C{G}_k=R_{k+1}-R_k$, the interval gradient has an extremely important characteristic that $|I{G}_k|\ge |C{G}_k|$ if the window $w_k$ is texture-free. So the indicator ${\mathsf{\Lambda }}_k=\frac{|I{G}_k|+ϵ}{|C{G}_k|+ϵ}$ (where $ϵ$ is a small constant to prevent numerical instability) behaves that ${\mathsf{\Lambda }}_k\ge 1$ for structural regions and ${\mathsf{\Lambda }}_k<1$ for textural regions. With this feature, an adaptive gradient rescaling rule is given, $C{G}_k^{\prime }=C{G}_k·min(1,{\mathsf{\Lambda }}_k)$, to suppress textural gradients while retaining structural gradients.

Based on above concepts, the procedure of 1D IGF mainly comprises four steps: (1) calculating $CG$ and $IG$ of R; (2) rescaling $CG$ to $C{G}^{\prime }$ via the aforementioned rule; (3) reconstructing an initial filtered result $\widehat{R}$ by simply accumulating $C{G}^{\prime }$; (4) Obtaining the final filtered result $R^{\prime }$ by filtering R with the guide of $\widehat{R}$. Accordingly, 2D images are processed by alternately performing the 1D filtering operations in row and column directions. For more details about the IGF algorithm, please refer to the original paper.

By leveraging the IGF method on the original image I, the desired guidance image $I^{\prime }$ is generated.

3.2.2. Fourier Spectrum Fusion

With the guidance image $I^{\prime }$ and the anomaly weight map W, the corrected Fourier spectrum $\tilde{\mathcal{F}}$ can be obtained by a fusion process

$$\begin{array}{c}\hfill \tilde{\mathcal{F}}=(1-W)\odot \mathcal{F}+W\odot {\mathcal{F}}^{\prime }\end{array}$$

(12)

where ⊙ is the operator of element-wise multiplication; $\mathcal{F}$ and ${\mathcal{F}}^{\prime }$ are Fourier spectrums of I and $I^{\prime }$, respectively.

Finally, the destriping result $\tilde{I}$ is reconstructed with $\tilde{\mathcal{F}}$ via IDFT.

3.3. Practical Considerations

Because of some potentially negative effects in the Fourier domain processing of 2D image, which mainly include the so-called edge effect and ringing artifact [38], practical considerations are needed in the implementation of our algorithm to deal with these issues effectively.

The edge effect inevitably occurs when the Fourier transform is performed on discrete images, owing to strong discontinuities of image borders. As a result, a well-known cross structure made of high energy coefficients along the axes is formed in the Fourier spectrum. This could bring adverse impacts to the performance of the proposed method. To reduce the edge effect and ensure good reconstruction results, we adopt a dual preprocessing step: first, we expand the input image I to a larger size via an amount of padding in each dimension; and then a special method called periodic plus smooth decomposition is employed to decompose the expanded image into a periodic component (that is, an image resembling the original image but whose DFT does not suffer from the edge effect) and a smooth component that has very slow variations in the image domain [54]. So the resulting periodic image is used as the actual input image in our detection and fusion process.

In destriping results via DFT, the ringing artifacts resulting from the discontinuity of intensity in the image sometimes appear along the highly-contrasted edges [38]. As for our method, despite ideal Fourier coefficients from the guidance image to correct the spectrum, slight ringing artifacts might still be observed because of a mismatch between Fourier coefficients in the spectrum fusion. Thus, as an optional coping strategy, we introduce a spectrum interpolation method based on total variation minimization to mitigate these possible artifacts, which is originally applied in Fourier filtering methods [49,55]. Here we shall not detail the method, but to make it workable within the proposed framework, we define the spectrum region to be interpolated, ${\mathsf{\Omega }}^{\prime }:|\tilde{\mathcal{F}}|./|\mathcal{F}|<0.99$, where the dominant Fourier coefficients in the fused spectrum are ruled out. More implementation details about the method can be found in the note of [55].

4. Experimental Results and Discussion

In this section, the performance of the proposed method is comprehensively evaluated. We focus on performing validations on two types of real-world TIR images we mentioned in Section 1.1, i.e., conventional TIR images and remote sensing TIR spectral images. To impartially demonstrate the effectiveness and robustness of our approach, abundant experiments are carried out on several publicly available datasets; at the same time, considering the differences of these two images in the research, we use the relevant contrast methods and common evaluation indexes for each case to compare and analyze the experimental results, respectively.

4.1. Experiments on Conventional TIR Images

For testing the proposed method on conventional TIR images, we utilize the open data coming from Tendero’s dataset (https://perso.telecom-paristech.fr/ytendero/) and Morris’s dataset (https://www.dgp.toronto.edu/~nmorris/IR/) and five state-of-the-art comparison methods including MIRE [19] (Midway Infrared Equalization), 1D-GF [16] (1D Guided Filtering), TwSF [52] (Two-Stage Filtering), DLSNUC [32] (Deep Learning-based Stripe Nonuniformity Correction), and WDNN [44] (Wavelet Deep Neural Network). Here three experiments are reported to illustrate the destriping performance of our method and their test images are shown in Figure 8. The images used in Experiments 1 and 2 are from Tendero’s dataset and record outdoor scenes at a resolution of 288 × 384, while the image used in Experiment 3 is an indoor image with size 288 × 384 and from Morris’s dataset. We can see that these real TIR images are obviously corrupted by vertical stripe noise.

Figure 9, Figure 10 and Figure 11 show the destriping results of three experiments. From the first two groups of results, some slight stripes remain in the results of MIRE, 1D-GF and TwSF, such as in the labeled regions with red ovals in Figure 9c,d and Figure 10b,c. So these results look somewhat rough. Besides, in 1D-GF’s and TwSF’s results in Experiment 2, vertical structure edges (highlighted by red arrows in Figure 10c,d) are smoothed to some extent. The results produced by DLSNUC, WDNN and our method have a globally smooth and fine appearance with edge preservation. But from a local viewpoint, the stripe noise is not yet eliminated completely in certain areas of the WDNN’s result (see the red ovals in Figure 9f and Figure 10f). From the third group of results, we can find that there are some undesirable artifacts in the results of MIRE, GF, and TwSF, as marked with red arrows in Figure 11b–d. Some slight stripes are also observed in the WDNN’s result. By comparison, our method and DLSNUC achieve superior visual effects.

Further, as shown in Figure 12, we extract the pixels in a single row from each set of experimental results to comparatively analyse the performance of test methods in suppressing the noise as well as preserving the sharp changes of image details. The red arrows indicate some notable segments. Here we put focus on the comparison of DLSNUC, WDNN and the proposed method. From these results, DLSNUC has a stronger ability of smoothing the noise at the expense of some detailed information. WDNN is inclined to keep more original details, but this also means that the stripe noise may not be completely suppressed. In contrast, the proposed method is moderate, which can not only efficiently remove stripe noise but also preferably maintain image details.

To corroborate our observations in the qualitative comparison, an objective assessment is also made with four quantitative indicators and they are the roughness index $\rho$ [10], the smoothing difference index $D_{S_T}^{S_F}$ [11], the image content metric Q [56], and the root-mean-square error (RMSE), respectively. In TIR image quality assessment, $\rho$ is usually used to measure the smoothness of the whole image, defined as

$$\begin{array}{c}\hfill \rho =\frac{{\parallel H\parallel }_1+{\parallel V\parallel }_1}{{\parallel I\parallel }_1}\end{array}$$

(13)

where H and V are horizontal and vertical gradients of the image I, respectively; ${\parallel ·\parallel }_1$ is the $L_1$ norm. A smaller $\rho$ value signifies a better smoothing effect. Another indicator $D_{s_T}^{s_F}$ is devoted to evaluating the algorithm’s ability in edge preservation and noise suppression by comparing the horizontal gradient smoothing difference between in structure and non-structure regions of TIR images, as follows

$$\begin{array}{c}\hfill D_{s_T}^{s_F}=\frac{\begin{array}{c}{\sum }_{k\in s_T}{H}^{\prime }\end{array}}{\begin{array}{c}{\sum }_{k\in s_T}H\end{array}}-\frac{\begin{array}{c}{\sum }_{k\in s_F}{H}^{\prime }\end{array}}{\begin{array}{c}{\sum }_{k\in s_F}H\end{array}}\end{array}$$

(14)

where H and $H^{\prime }$ correspond to horizontal gradients of images before and after destriping treatment; $s_T$ and $s_F$ are pixels in structure and non-structure regions, which are classified by the measurement of 1D horizontal differential statistics. $D_{s_T}^{s_F}$ generally yields positive values, with higher values indicating better edge-preserving effect while destriping. But if the value is less than 0, it suggests that some potential exceptions occur in the result. Since $\rho$ and $D_{s_T}^{s_F}$ essentially assess the image from a global perspective, we introduce the image patch based no-reference metric Q to locally evaluate the image quality. Concretely, the Q value of an image patch (with size $8\times 8$ by default) is computed by

$$\begin{array}{c}\hfill Q=s{v}_1\frac{s{v}_1-s{v}_2}{s{v}_1+s{v}_2}\end{array}$$

(15)

where $s{v}_1$ and $s{v}_2$ are the first and second singular values of gradient matrix of the image patch. The mean Q value of all anisotropic patches is defined as the metric for the whole image. The higher the Q value is, the better quality the image has. What’s more, the error measurement RMSE is adopted to reflect the pixel accuracy of the destriped image. Its expression is

$$\begin{array}{c}\hfill \mathrm{RMSE}=\sqrt{\frac{1}{XY}\sum _{k=1}^{XY}{(I_k-I_k^{\prime })}^2}\end{array}$$

(16)

The quantitative evaluation results of three experiments for conventional TIR images are listed in

Table 1. Quantitative evaluation results of three experiments for conventional TIR images.

Method	Experiment 1				Experiment 2				Experiment 3
	$\mathit{\rho }$	${\mathit{D}}_{{\mathit{s}}_{\mathit{T}}}^{{\mathit{s}}_{\mathit{F}}}$	Q	RMSE	$\mathit{\rho }$	${\mathit{D}}_{{\mathit{s}}_{\mathit{T}}}^{{\mathit{s}}_{\mathit{F}}}$	Q	RMSE	$\mathit{\rho }$	${\mathit{D}}_{{\mathit{s}}_{\mathit{T}}}^{{\mathit{s}}_{\mathit{F}}}$	Q	RMSE
MIRE	0.138	0.160	9.82	0.0217	0.125	0.328	8.66	0.0286	0.177	−0.470	17.08	0.0059
1D-GF	0.131	0.139	9.11	0.0209	0.120	0.223	7.83	0.0282	0.164	−0.284	11.89	0.0068
TwSF	0.131	0.149	9.09	0.0231	0.119	0.235	7.81	0.0318	0.166	−0.382	12.46	0.0056
DLSNUC	0.100	0.499	11.05	0.0220	0.088	0.567	10.71	0.0292	0.126	0.255	18.07	0.0057
WDNN	0.129	0.256	9.34	0.0207	0.118	0.421	8.55	0.0281	0.166	−0.021	15.70	0.0047
Ours	0.126	0.526	11.47	0.0200	0.114	0.618	11.36	0.0271	0.162	0.283	19.32	0.0042

. In terms of $\rho$, DLSNUC has the lowest values, followed by our method and WDNN. This is in line with the conclusion in the qualitative analysis that our method does produce a moderately smooth result in between DLSNUC and WDNN. In respects of $D_{s_T}^{s_F}$ and Q, our method obtain the highest scores, illustrating its advantage in balancing between stripe noise removal and detail preservation. Moreover, the lowest RMSE values demonstrate that our method has better capability to hold the pixel accuracy of image when reducing stripe noise.

In our performance test, the extreme case that conventional TIR images are corrupted by very serious stripe noise is taken into consideration as well. Figure 13 provides one such example, where the test image is captured by Mars Orbiter satellite [44]. We can observe from the figure that there remain more or less striping effects in the other five results, compared to our result. This further reflects the practicability and superiority of the proposed method.

In summary, the proposed method exhibits the highly effective and robust destriping performance in a variety of real experiments on conventional TIR images.

4.2. Experiments on Remote Sensing Spectral Images

To demonstrate that our destriping approach is also applicable to remote sensing spectral images and has the consistently excellent performance, three sets of real striped multispectral data (see Figure 14), taken from Terra MODIS Level-1B products that are freely available on the website https://ladsweb.modaps.eosdis.nasa.gov, are tested with four pertinent methods for comparison. In Experiments 1 and 2, we separately selected a sequence of atmosphere images and land images within six TIR bands (Band 21, 24, 27, 33, 34 and 35), and they are extracted from the MOD021KM product and saved as 8-bit images with a cropped size of 500 × 500. From Figure 14, it can be seen that there are multiple types of stripes in these real spectral images: Band 21 and Band 24 images are degraded by periodically even stripes; Band 27 and Band 34 images suffer severely uneven stripes, while Band 33 and Band 35 images have mildly uneven stripes. The test data used in Experiment 3, including Band 27, 34 and 36 images, are collected from the MOD02SSH product and have a size of 270 × 406. In this group, the striping effect is remarkably inhomogeneous and mixed with some random noise, especially at Band 36. Additionally, the four single-image-based contrast methods for destriping of remote sensing images are as follows: WFAF [37] (Wavelet Fourier Adaptive Filtering), SLD [18] (Statistical Linear Destriping), LRSID [26] (Low Rank Single Image Decomposition), DSM-$l_0$ [30] (Directional $l_0$ Sparse Model).

Figure 15, Figure 16 and Figure 17 show the destriping results of three experiments. Meanwhile for a better observation and comparison, we deliberately picked out a couple of image parts from the results of Experiments 1 and 2 to magnify and highlight visual differences, which are shown in Figure 18 and Figure 19. According to these results, some conclusions can be drawn. WFAF and SLD restrain stripes to certain extent, but their overall destriping effects are still inferior and many residual stripes are visible in their results. Also seeing from Experiment 2, WFAF could result in some light or shaded ribbonlike artifacts when there are a short, sharp variation in image contrast along the stripe direction. LRSID performs well in the case that the stripe in one column of image is steady and even, but its efficiency will deteriorate if the stripe case is less than ideal. For example, in dealing with the Band 27, 34 and 36 images, the results of LRSID are unsatisfactory because of fragmented and uneven stripes that are not perfectly fitted with the low rank prior, as the paper pointed out. DSM-$l_0$ has a good robust performance in the multifarious stripe cases, and its results are globally smooth without any stripes. However, such processing also causes the loss of a number of image detail information. The most notable cost is that many subtle gray changes are erased in image contents. By comparison, our approach unfurls a consistently impressive ability to work with various stripe cases, where not only are stripes effectively removed but also image details are well protected. It is especially noteworthy that, when the original image is damaged with a mix of stripes and random noise, the proposed method is still competent and achieves a better visual result than those of the other four algorithms, like the results of Band 36 showed in Figure 17.

Also, we employ the mean cross-track profile and the mean row power spectrum to further verify the algorithm performance, which are two commonly qualitative analysis indicators for destriping of remote sensing images. Here the Band 21 image from Experiment 1 and the Band 27 image from from Experiment 2 are taken as example, and their mean cross-track profiles and mean row power spectrums are shown in Figure 20, Figure 21, Figure 22 and Figure 23. Generally, the mean cross-track profile of a good destriped result should be a smoothed version of that of the original image, with optimally preserving primary structure trends. But we can notice from Figure 20 and Figure 21 that the WFAF’s profile remains many obtrusive burrs, indicating a poor destriping effect. Conversely, the DSM-$l_0$’s profile is over-smoothed and fails to keep some underlying image details well. The profiles of SLD, LRSID and the proposed approach have roughly semblable trends, but our result is more smooth than that of SLD and retains more finely structural fluctuations than that of LRSID. So in terms of the mean cross-track profile, the proposed method obtains the desired result both on smoothing stripes and maintaining image details. On the other hand, from the perspective of power spectrum (see Figure 22 and Figure 23), noise impulses are not reduced completely in the results of WFAF and SLD, while LRSID, DSM-$l_0$ and the proposed method produce similar and good results. Especially from Figure 23, after removing fragmentary and uneven stripes, the power spectrum of our result is more well-distributed with power preservation than other four methods.

To make a quantitative analysis of the above results, three frequently-used no-reference indices are introduced. Specifically, we adopt the inverse coefficient of variation (ICV) [13,14] and the mean relative deviation (MRD) [15] to assess the ability of destriping methods on balancing noise cancellation and detail preservation. The two indexes are defined as

$$\begin{array}{c}\mathrm{ICV}=\frac{H{R}_m}{H{R}_s}\hfill \end{array}$$

(17)

$$\begin{array}{c}\mathrm{MRD}=\frac{1}{M}\sum _k\frac{|S{R}_k^{\prime }-S{R}_k|}{S{R}_k}\hfill \end{array}$$

(18)

where $H{R}_m$ and $H{R}_s$ are the mean and standard deviation of pixel values in a homogeneous region; $S{R}_k$ and $S{R}_k^{\prime }$ are the k-th pixel pair before and after destriping in a sharp region, and M is the total number of pixels. The ICV index is designed to measure the smoothness of a homogeneous region, while the MRD index reflects the relative distortion of a sharp region. So a good destriping performance should be the combination of a large ICV value and a small MRD value. Apart from evaluating in the dimension of image, we also give consideration to the influence of the destriping operation on the spectral fidelity. Here the common metric in the spectral imaging filed, spectral angle mapper (SAM) [57,58], is utilized to numerically evaluate the distortion between the original and destriped spectral signatures. Its definition is as follows

$$\begin{array}{c}\hfill \mathrm{SAM}={cos}^{-1}\frac{{\mathit{s}}_k^T·{\mathit{s}}_k^{\prime }}{∥{\mathit{s}}_k∥·∥{\mathit{s}}_k^{\prime }∥}\end{array}$$

(19)

where ${\mathit{s}}_k$ and ${\mathit{s}}_k^{\prime }$ stand for the k-th spectral signatures of the original and destriped spectral images, respectively. In general, the smaller the SAM value is, the less the distortion is between ${\mathit{s}}_k$ and ${\mathit{s}}_k^{\prime }$.

Table 2. Quantitative evaluation of destriping methods on remote sensing spectral images.

Image	Index	Original	WFAF	SLD	LRSID	DSM-${\mathit{l}}_0$	Ours
Band 21	ICV1	21.8910	28.0554	28.0749	31.8895	36.5370	31.4056
	ICV2	18.5833	19.6632	21.8711	24.0064	29.5712	24.6683
	MRD1	-	0.0167	0.0182	0.0253	0.0376	0.0248
	MRD2	-	0.0131	0.0139	0.0157	0.0320	0.0135
Band 24	ICV1	3.8014	3.8561	3.8729	3.9096	4.4591	4.1827
	ICV2	87.8540	122.4623	125.2533	151.5480	154.6837	149.2865
	MRD1	-	0.0224	0.0268	0.0350	0.0980	0.0333
	MRD2	-	0.0106	0.0087	0.0110	0.0212	0.0101
Band 27	ICV1	14.5392	17.2013	16.5835	17.8278	27.7905	18.5027
	ICV2	38.7568	54.0097	54.5366	65.7343	125.2244	73.0345
	MRD1	-	0.0204	0.0211	0.0286	0.0356	0.0258
	MRD2	-	0.0259	0.0302	0.0402	0.0544	0.0372
Band 33	ICV1	5.5825	5.4789	5.5875	5.7910	8.7910	5.6949
	ICV2	108.4053	164.1888	166.2775	231.5580	253.2472	212.7835
	MRD1	-	0.0102	0.0067	0.0136	0.0169	0.0124
	MRD2	-	0.0044	0.0034	0.0074	0.0115	0.0065
Band 34	ICV1	14.4918	15.4177	15.7350	15.9210	18.2459	16.3604
	ICV2	40.4405	54.1095	58.4728	77.3437	127.7196	90.1767
	MRD1	-	0.0319	0.0482	0.0651	0.1450	0.0585
	MRD2	-	0.0102	0.0110	0.0131	0.0147	0.0126
Band 35	ICV1	7.8907	8.2272	8.3232	9.5049	13.2531	10.3291
	ICV2	60.8430	70.1713	75.3580	92.3648	153.1217	92.1730
	MRD1	-	0.0115	0.0151	0.0289	0.0344	0.0252
	MRD2	-	0.0084	0.0099	0.0179	0.0414	0.0175
	SAM1	-	0.0261	0.0218	0.0293	0.0533	0.0214
	SAM2	-	0.0149	0.0124	0.0143	0.0178	0.0121

lists the ICV, MRD and SAM values of all the destriped spectral bands from Experiments 1 and 2. In terms of ICV and MRD indexes, it can be found that, DSM-$l_0$ achieves not only all the largest ICV values but also the largest MRD values, while WFAF and SLD are in the opposite situation, corresponding to the smallest ICV and MRD values alternately. Both results are unacceptable because they actually indicate the oversmoothing effect of DSM-$l_0$ and the undersmoothing effect of WFAF and SLD, combined with the visual results. By comparison, LRSID and the proposed method obtain modest, balanced ICV and MRD values for six bands. But the proposed method has higher ICV values and lower MRD values than LRSID in Band 27 and Band 34 where the stripe noise is severe and irregular. It shows the better coping capacity of the proposed method in these two cases. So judging from the results of ICV and MRD, our method reaches a good and robust trade-off between stripe noise suppression and detail preservation. In addition, from the SAM perspective, the proposed method gets the lowest values, which means that the spectral distortion caused by the method is slightest, compared with other destriping methods.

Overall, for remote sensing spectral images, the proposed method also can maintain the stable and competitive destriping performance with minimum loss of the spatial-spectral information.

4.3. Discussion

4.3.1. Parameter Analysis

In the implementation of the proposed method, there are two tunable parameters that significantly affect performance, the taper angle $\alpha$ and the standard deviation $\sigma$. Specifically, $\alpha$ determines the size of the tapered detection region $\mathsf{\Omega }$ in the stage of anomaly detection. So when $\alpha$ is tuned from a small degree to a relatively large degree, it means that more potential high-frequency anomalies may be detected and restrained. To illustrate the influence of the parameter $\alpha$ on the destriping performance, Figure 24 gives an example, where the original image has even stripes as well as some uneven fragmented stripes, such as in the area marked with a red square. From the destriping results with different $\alpha$ values, we can find that, when $\alpha$ has a small degree ($2^{\circ }$∼6${}^{\circ }$), the even regular stripes are removed well but the uneven fragmented stripes still remain in local areas; all stripes are suppressed when $\alpha$ is set to a relatively large degree (${10}^{\circ }$∼16${}^{\circ }$). Thus a large $\alpha$ is favorable for the given case. On the contrary, for most general cases, a small $\alpha$ is sufficient and suggested because more detailed information could be saved. In our implementation, $\alpha$ is fixed as ${10}^{\circ }$ by default.

The other parameter $\sigma$ regulates the generation of the guidance image in the stage of spectrum fusion. Actually, $\sigma$ is involved in the IGF method [53] and it determines the degree of filtering, that is, the filtered image becomes more smooth as $\sigma$ increases. Since the guidance image is supposed to be a stripe-free base image, $\sigma$ has to be properly adjusted with the level of actual stripe noise in the implementation of our approach. For most general cases, $\sigma$ is set in the range of 0.5∼3. It should be noted that, on the premise of ensuring the guidance image without stripes, a small $\sigma$ is suggested so that the guidance image keeps more fine details, which could be fused into the result.

4.3.2. Running Time

We also consider running times of different methods to examine their computation complexity. In this paper, all the experiments are conducted on a laptop personal computer with an Intel Core 2.4-GHz CPU and 8-GB RAM and the codes of algorithms are implemented in Matlab expect WDNN in Python.

Table 3. Average running time (s) of different methods for conventional TIR images with size 288 × 384 and remote sensing spectral images with size 500 × 500.

Image Size	MIRE	1D-GF	TwSF	DLSNUC	Ours	Image Size	WFAF	SLD	LRSID	DSM-${\mathit{l}}_0$	Ours
288 × 384	0.07	0.06	0.01	0.46	0.38	500 × 500	0.25	0.04	16.37	34.52	1.06

records the average running times of different methods in experiments on conventional TIR images with size 288 × 384 and remote sensing spectral images with size 500 × 500. It can be concluded that statistics-based methods like MIRE and SLD and filtering-based methods such as 1D-GF, TwSF and WFAF are relatively simple and have the low cost of computation; while optimization-based methods like LRSID and DSM-$l_0$ and deep-learning-based methods including DLSNUC and WDNN are usually complicated and require a large quantity of operation. By comparison, the time consumption of our method is modest and acceptable.

5. Conclusions

In this paper, we have proposed a novel Fourier domain anomaly detection and spectral fusion approach for the destriping of TIR imagery. First, based on spectral characteristics of both TIR images and stripe noise, the abnormal frequencies that are likely to be caused by stripe noise are detected to generate an anomaly weight map. With this map, the original spectrum is then fused with the spectrum of a guidance image, which is estimated by a proper de-texture filter. As a result, in the ”fused” image, not only is the stripe noise removed effectively, but original structures and details are preserved well. Lots of experiments are performed on two types of real striped TIR data including conventional TIR images and remote sensing spectral images. Both qualitative and quantitative results demonstrate the excellent destriping performance of the proposed method with less loss of the original information, compared with the other state-of-the-art destriping methods in the two fields.