Variational Destriping in Remote Sensing Imagery: Total Variation with L¹ Fidelity

Abstract

This paper introduces a variational method for destriping data acquired by pushbroom-type satellite imaging systems. The model leverages sparsity in signals and is based on current research in sparse optimization and compressed sensing. It is based on the basic principles of regularization and data fidelity with certain constraints using modern methods in variational optimization, namely, total variation (TV), $L^1$ fidelity, and the alternating direction method of multipliers (ADMM). The proposed algorithm, TV–$L^1$, uses sparsity-promoting energy functionals to achieve two important imaging effects. The TV term maintains boundary sharpness of the content in the underlying clean image, while the $L^1$ fidelity allows for the equitable removal of stripes without over- or under-penalization, providing a more accurate model of presumably independent sensors with an unspecified and unrestricted bias distribution. A comparison is made between the TV–$L^2$ model and the proposed TV–$L^1$ model to exemplify the qualitative efficacy of an $L^1$ striping penalty. The model makes use of novel minimization splittings and proximal mapping operators, successfully yielding more realistic destriped images in very few iterations.

1. Introduction

Image striping is a well-known phenomenon that arises in multi-detector imaging systems ranging from pushbroom-type instruments, such as the Airborne Multi-angle Spectro Polarimetric Imager (AirMSPI), to atomic force mircroscopy (AFM). Biases in lateral detection occur as a result of response variation in spatial detectors, such as in satellite imaging systems, or temporal changes, such as in raster scans. Although these systems are optimally precalibrated, post-processing, such as destriping, of data is a prerequisite for accurate and valid analyses. Striping removal has been traditionally performed using either statistically based methods [1,2] or low-pass filtering in the frequency domain [3,4,5,6,7]. These methods, however, do not remove stripes completely and have effects of blurring the image. More recently, wavelet-based filtering methods have been proposed [8,9,10]. However, such methods also blur the images and produce ringing effects in reconstruction.

In [11], the authors introduced a robust destriping algorithm with a spatially adaptive unidirectional total variation (TV) model. The authors developed a new destriping method that combines the TV–Stokes model and unidirectional TV model in [12]. In [13], the authors introduced a TV and framelet regularization model for destriping. The same authors proposed an anisotropic spectral–spatial TV regularization to enhance the smoothness of the solution along both the spectral and spatial dimension in [14]. In [13,14], the authors employed the $L^2$ norm as the fidelity term, and we refer to these models as TV–$L^2$ in this paper.

We follow the pedigree of variational and partial differential equation (PDE)-based methods applied to images [15,16] in order to construct a well-defined, optimizable model yielding fast and high-quality destriping. The focus of our paper is not on creating a sparse wavelet representation of the destriped image, but rather on how to remove the optimal striping mask while preserving high image fidelity. In this paper, we propose using the $L^1$ penalty for striping size.

Our research is robust to both isotropic and anisotropic versions of TV, whereas [13] argue that the anisotropic case is the only appropriate case. While it is true that the anisotropic case uses decoupled energy for the measure of smoothness and is therefore easier to minimize, isotropic TV is not selective in which direction smoothness is penalized. Image content smoothness (or lack thereof) is not known a priori, and thus no preference should immediately be given to certain directions for evaluating smoothness. Using the $L^1$ penalty, and depending on the data, the isotropic TV, which theoretically uses more local information, allows for qualitatively better, less invasive and more intelligent destriping. We compare the $L^1$ and $L^2$ penalties, ultimately favoring the $L^1$ as a result of a wider yet tighter distribution of the striping mask. We include detailed derivations and a motivated evolution of the optimization problem with pedagogy in mind so that the proposed TV–$L^1$ method, along with its alternating direction method of multipliers (ADMM) (split Bregman) optimization, can be accessible to all academic disciplines involved with image processing.

We construct a variational model that is well-defined, qualitatively motivated, and easily minimized. The constructed energy uses sparsity-promoting energy functionals, on the basis of TV and $L^1$ energy, to achieve minimally invasive destriping. Both isotropic and anisotropic TV, along with $L^1$ energy, are considered in our variational model. The ADMM is used in conjunction with nonlinear proximal operators to efficiently optimize the energy, yielding quick and high-quality results.

2. Materials and Methods

In this section, we describe the striping problem, provide the background and motivations, discuss tools from functional analysis to be used in describing the model, develop the proposed TV–$L^1$ model, and provide the algorithm.

2.1. Striping Structure

Let $U(x,y)$ be a stripe-free image of size R by C, and let $G(y)$ be a multiplicative stripe noise of length R. Then the observed image, F, can be written as

$$F(x,y)=G(x,y)U(x,y)$$

(1)

Taking logarithms of both sides of Equation (1) yields an additive structure more suitable for energy minimization methods. The model can now be written as

$$f(x,y)=g(y)+u(x,y)$$

where $f(x,y)=ln(F(x,y))$, $g(x,y)=ln(G(x,y))$, and $u(x,y)=ln(U(x,y))$.

Striping in images can be viewed as structured noise, of which variations are mainly concentrated along one axis. This can be mathematically encoded as $\parallel {\nabla }_{x}G\parallel \ll \parallel {\nabla }_{y}G\parallel$, or with the logarithmic terms, as $\parallel {\nabla }_{x}g\parallel \ll \parallel {\nabla }_{y}g\parallel$.

2.2. Tikhonov Minimization

A classical Tikhonov minimization problem would consist of a smoothness regularizer and a data fidelity term, both easily differentiable, with the striping constraint [17]:

$$\begin{array}{c}\underset{u}{min}\left\{{\parallel \nabla u(x,y)\parallel }_2^2+\frac{\lambda }{2}{\parallel u(x,y)-f(x,y)\parallel }_2^2\right\}\hfill \\ \mathrm{such}\mathrm{that}\parallel {\nabla }_{x}g\parallel \ll \parallel {\nabla }_{y}g\parallel \hfill \end{array}$$

This constraint can be simplified by the approximation that ${\nabla }_{x}G(x,y)=0\forall (x,y)$, which would make $G(x,y)=G(y)$, and $g(x,y)=g(y)$ functions of only one variable. Using the additive identity between f, g, and u along with the constraint approximation, the new unconstrained minimization problem is

$$\underset{g}{min}\left\{\parallel {\nabla }_y{(f(x,y)-g(x,y))\parallel }_2^2+\frac{\lambda }{2}{\parallel g(x,y)\parallel }_2^2\right\}$$

By taking the first variation of the energy and setting it to zero, the closed-form solution to this minimization problem is

$$\begin{array}{cc}\hfill g(x,y)& ={({\nabla }_y·{\nabla }_y+\lambda I)}^{-1}({\nabla }_y·{\nabla }_{y}f(x,y))\hfill \\ & ={(\frac{{\partial }^2}{\partial y^2}+\lambda I)}^{-1}(f_{yy}(x,y))\hfill \end{array}$$

However, this solution would cause g to become bivariate, in contradiction with the constraint. Instead, using the Cartesian regularity of our rectangular domain $\Omega =I_x\times I_y$ and using $g(x,y)=g(y)$, we can come to a solution that is in agreement with the constraint by integrating with respect to x:

$$\begin{array}{c}\hfill {\int }_{\Omega }g(x,y)dx={\int }_{I_x}g(y)dx=g(y){\int }_{I_x}dx=g(y)\mu (I_x)={\int }_{I_x}{(\frac{{\partial }^2}{\partial y^2}+\lambda I)}^{-1}(f_{yy}(x,y))dx\Rightarrow \end{array}$$

$$g(y)=\frac{1}{\mu (I_x)}{\int }_{I_x}{(\frac{{\partial }^2}{\partial y^2}+\lambda I)}^{-1}(f_{yy}(x,y))dx$$

2.3. Fourier Interpretation

Utilizing the Plancherel Fourier isometry, the solution can be interpreted in spectral form as

$$\widehat{g}({\omega }_y)=\frac{1}{\mu (I_x)}\left(\frac{{\omega }_y^2}{\lambda +{\omega }_y^2}\right)\widehat{f}=\left(\frac{{\omega }_y^2}{\lambda +{\omega }_y^2}\right)\overline{\widehat{f}}$$

For a specific x, the stripping $g(y)$ is constant and is of higher frequency, whereas the underlying clean image varies more slowly (has more low-frequency content) and for each x has somewhat different frequencies. Therefore, the average frequencies of the clean image are low in magnitude and are of lower frequency, while the average frequencies of the stripping are high in magnitude and are of higher frequency. Therefore, the average frequencies (averaged over ${\omega }_x$) of the cleaned image are simply the average frequencies of the original image multiplied by a one-dimensional low-pass filter $\frac{\lambda }{\lambda +z^2}$. Likewise, the striping mask on the spectral side, $\widehat{g}$, is obtained analogously to a one-dimensional high-pass filter $\frac{z^2}{\lambda +z^2}$.

Although this minimization problem is readily solvable in closed form and has a motivated physical interpretation, we must abandon the quadratic energy terms so that we may have less penalization for heavier striping and to allow for less smooth solutions. Although the differentiability of terms is easy to work with, enough optimization machinery has been developed that we may tread forward without it. We now investigate and outline some tools from signal processing in order to refine our model. Stripe and ring artifact removal from this frequency perspective has been accomplished using wavelet and Fourier filtering [18].

2.4. Tools from Functional Analysis

We now introduce the notations and tools from functional analysis that are used in Section 2.5.

2.4.1. Total Variation: Anisotropic vs. Isotropic

The idea of using TV as a regularizer and denoiser that promotes sparsity and piecewise constant smoothness can be traced back to Rudin, Osher, and Fatemi [19,20]. We begin with defining the notion of TV, which is used as a regularizer in the model.

Definition 1.

The TV of a function $f\in L^1(\Omega )$ is

$$V(f,\Omega ):=sup\left\{{\int }_{\Omega }f(x)\mathrm{div}\varphi (x)\mathrm{d}x:\varphi \in C_c^1(\Omega ,{\mathbb{R}}^n),{\parallel \varphi \parallel }_{L^{\infty }(\Omega )}\le 1\right\}$$

For a differentiable function $f\in \Omega$, with $\Omega \subseteq {\mathbb{R}}^n$, the TV of f can be written as

$$V(f,\Omega )=\underset{\Omega }{\int }\left|\nabla f(x)\right|\mathrm{d}x$$

The choice of vectorial norm inside the integral yields the two different types of TV.

Definition 2.

Isotropic TV: $|·|$ denotes the $l_2$-norm, in which case

$$V(f,\Omega )=\underset{\Omega }{\int }{(\sum _i^n{f}_{x_i}^2(x))}^{\frac{1}{2}}\mathrm{d}x$$

Definition 3.

Anisotropic TV: $|·|$ denotes the $l_1$-norm, in which case

$$V(f,\Omega )=\underset{\Omega }{\int }\sum _i^n|f_{x_i}(x)|\mathrm{d}x$$

The isotropic and anisotropic cases differ in terms of the geometries they each preserve. While the decoupled anisotropic TV preserves piecewise constant orthogonal structures, such as rectangular roofs, the coupled isotropic TV preserves piecewise constant radial structures, such as silos. Our model is robust with respect to either choice of TV, and dual derivations of variable updates are shown; however, the experiments and results are based on the anisotropic definition.

2.4.2. Shrinkage Proximal Operator

We introduce a splitting variable and quadratic penalty into the model. The solution to the $l_1$-regularized least-squares problem:

$$\underset{\overrightarrow{x}}{arg\; min}\mu \parallel \overrightarrow{x}{\parallel }_1+\frac{1}{2}{\parallel \overrightarrow{x}-\overrightarrow{y}\parallel }_2^2$$

is given by the soft threshold proximal mapping operator, shrinkage [21,22]:

Definition 4.

$$Shrink(\overrightarrow{x},\mu )=S_{\mu }(\overrightarrow{x})=\frac{\overrightarrow{x}}{|\overrightarrow{x}|}max\{|\overrightarrow{x}|-\mu ,0\}$$

If $\parallel \overrightarrow{x}{\parallel }_1=|x_1|+|x_2|$, as in the anisotropic case of TV, the shrinkage is decoupled and done component wise. On the other hand, if $\parallel \overrightarrow{x}{\parallel }_1=\sqrt{|x_1{|}^2+{|x_2|}^2}$, as in the isotropic case, the terms are coupled and both components are updated simultaneously. Both variants have their merits; while the former is computationally simpler, the latter has the advantage of using more local information and may be more conformant to certain image processing applications.

2.5. TV–$L^1$ ADMM Model

We now describe our model using the motivations and tools introduced earlier. The two energy components of the minimization problem are the smoothness regularizer and the data fidelity term. The energy of the data fidelity term, $\frac{\lambda }{2}{\parallel f(x,y)-u(x,y)\parallel }_2^2=\frac{\lambda }{2}{\parallel g(y)\parallel }_2^2$, can be interpreted as the magnitude of the striping mask. The $L^2$ fidelity overly penalizes stripes of large magnitude and likewise underexaggerates the significance of stripes of small magnitude. In areas of no striping, we intend our (logarithm of the) striping mask to be very close to zero, while in areas of heavy striping, we wish to remove such striping and thus yield a larger magnitude of our striping mask in that region. The $L^1$ fidelity gives a smaller striping mask in areas of no striping, leaving enough energy to remove the heavier striping in localized areas of the image. Because there is no prior knowledge of the distribution of the stripes and qualitatively we may wish to remove deep striping effects while preserving sharp geometry, we believe it is better to update the model with an $L^1$ striping penalty, ${\parallel g\parallel }_1$.

An $L^2$ gradient term would cause oversmoothing of the retrieved clean image $u(x,y)$. This could cause a loss in boundary sharpness of elements in the image (e.g., lakes, rooftops, etc.), which seems important in the pursuit and usage of destriped images. Implementing a TV-based regularizer would give performance similar to the $L^2$ gradient but maintains boundary sharpness more natural to the underlying image. Although these terms are not differentiable, impeding a closed-form solution, state-of-the-art nonlinear optimization algorithms are available for fast convergence to qualitatively meaningful minimizers. The unconstrained TV $L^1$ model (TV–$L^1$) is

$$\underset{u}{min}\left\{V(u(x,y),I_x\times I_y){+\lambda ||u(x,y)-f(x,y)\parallel }_1\right\}$$

or equivalently, minimizing with respect to the striping mask g,

$$\underset{g}{min}\left\{V(f(x,y)-g(y),I_x\times I_y)+\lambda {\parallel g(y)\parallel }_1\right\}$$

2.5.1. Discretization

For the purpose of application and computation, we now move the problem into a discrete setting. Let $\Omega =\{x_1,\dots ,x_C\}\times \{y_1,\dots ,y_R\}$ be an $R\times C$ image. First variations are approximated via forward differences, so that $\frac{\partial f}{\partial y}(x_i,y_j)\approx f(x_i,y_{j+1})-f(x_i,y_j):={\delta }_y{f}_{i,j}$ for $j=1,\dots ,R-1$ and analogously $\frac{\partial f}{\partial x}\approx f(x_{i+1},y_j)-f(x_i,y_j):={\delta }_x{f}_{i,j}$ for $i=1,\dots ,C-1$. We take Neumann boundary conditions, so that on the forward boundary ($i=C$ or $j=R$), the derivative is set to zero.

$$\begin{array}{c}\hfill D=\left[\begin{array}{cccccc}-1& 1& 0& \cdots & 0& 0\\ 0& -1& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & -1& 1\\ 0& 0& 0& \cdots & 0& 0\end{array}\right]\in M_{R\times R}(\mathbb{R})\\ \hfill \mathrm{so}\mathrm{that}Du(x_{i_0},y)={({\delta }_y{u}_{i_0,1},\cdots ,{\delta }_y{u}_{i_0,R-1},0)}^T\end{array}$$

The isotropic TV is

$$|{\nabla }^{I}f(x_i,y_j)|\approx \sqrt{{({\delta }_x{f}_{i,j})}^2+{({\delta }_y{f}_{i,j})}^2}$$

The anisotropic TV is

$$|{\nabla }^{AI}f(x_i,y_j)|\approx |{\delta }_x{f}_{i,j}|+|{\delta }_y{f}_{i,j}|$$

With these discrete operators defined, the discrete unconstrained TV–$L^1$ minimization problem is

$$\underset{g}{min}\left\{\sum _{i,j}\parallel \langle {\delta }_x{f}_{i,j},{\delta }_y(f_{i,j}-g_j)\rangle {\parallel }_1+\lambda {\parallel g(y)\parallel }_1\right\}$$

The two types of the minimization problem are the following:

Anisotropic:

$$\underset{g}{min}\left\{\sum _{i,j}|{\delta }_y(f_{i,j}-g_j)|+\lambda \parallel g(y){\parallel }_1\right\}$$

Isotropic:

$$\underset{g}{min}\left\{\sum _{i,j}\sqrt{{({\delta }_x{f}_{i,j})}^2+{({\delta }_y(f_{i,j}-g_j))}^2}+\lambda {\parallel g(y)\parallel }_1\right\}$$

2.5.2. Augmented Lagrangian

Anisotropic: With the discrete forward difference approximation matrix defined above, we can rewrite the minimization problem as follows:

$$\begin{array}{cc}\hfill \mathit{Point}\mathit{Form}:& \underset{\overrightarrow{g}}{min}\left\{\sum _{i,j}|{\delta }_y(f_{i,j}-g_j)|+\lambda \parallel g(y){\parallel }_1\right\}\hfill \\ \hfill \mathit{Vector}\mathit{Form}:& \underset{\overrightarrow{g}}{min}\left\{\sum _i\parallel D(f_i-\overrightarrow{g}){\parallel }_1+\lambda {\parallel g(y)\parallel }_1\right\}\hfill \\ \hfill \mathit{Matrix}\mathit{Form}:& \underset{\overrightarrow{g}}{min}\left\{\parallel D(f-\overrightarrow{g}\otimes \overrightarrow{1_C}){\parallel }_{1,1}+\lambda {\parallel g(y)\parallel }_1\right\}\hfill \end{array}$$

To render the constrained minimization problem unconstrained, we introduce auxiliary variables, Lagrangian multipliers (split Bregman), and quadratic penalty terms, so that the augmented Lagrangian is defined as follows:

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\alpha ,\lambda }(b_i,h,g,q_i,r)& =& \sum _i\left(\parallel b_i{\parallel }_1+\frac{\alpha }{2}{\parallel b_i-D(g-f_i)\parallel }_2^2+\langle q_i,b_i-D(g-f_i)\rangle \right)\hfill \\ & & +\lambda \left({\parallel h\parallel }_1+\frac{\alpha }{2}{\parallel h-g\parallel }_2^2+\langle r,h-g\rangle \right)\hfill \\ & =& \sum _i\left(\parallel b_i{\parallel }_1+\frac{\alpha }{2}{\parallel b_i-D(g-f_i)+\frac{q_i}{\alpha }\parallel }_2^2\right)\hfill \\ & & +\lambda \left({\parallel h\parallel }_1+\frac{\alpha }{2}{\parallel h-g+\frac{r}{\alpha }\parallel }_2^2\right)+\mathcal{O}(q_i^2,r^2)\hfill \end{array}$$

We now solve the unconstrained saddle point problem:

$$\underset{b_i,h,g}{min}\underset{q_i,r}{max}{\mathcal{L}}_{\alpha ,\lambda }(b_i,h,g,q_i,r)$$

The solution to the original constrained minimization problem is now found as the saddle point of the augmented Lagrangian $\mathcal{L}$ in a sequence of iterative suboptimizations called the ADMM [23,24,25,26,27,28].

The splitting variables $b_i$ and h are updated by the proximal mapping operator:

$$\begin{array}{cc}\hfill b_i^{k+1}& =\underset{b_i}{arg\; min}{\mathcal{L}}_{\alpha ,\lambda }(b_i,h^k,g^k,q_i^k,r^k)\hfill \\ & =\underset{b_i}{arg\; min}\left\{\parallel b_i{\parallel }_1+\frac{\alpha }{2}{\parallel b_i-D(g-f_i)+\frac{q_i}{\alpha }\parallel }_2^2\right\}\hfill \\ & =S_{\frac{1}{\alpha }}\left(D(g-f_i)-\frac{q_i}{\alpha }\right)\hfill \\ \hfill h^{k+1}& =\underset{h}{arg\; min}{\mathcal{L}}_{\alpha ,\lambda }(b_i^{k+1},h,g^k,q_i^k,r^k)\hfill \\ & =\underset{h_i}{arg\; min}\left\{{\parallel h\parallel }_1+\frac{\alpha }{2}{\parallel h-g\parallel }_2^2\right\}\hfill \\ & =S_{\frac{1}{\alpha }}\left(g-\frac{r}{\alpha }\right)\hfill \end{array}$$

Because of the introduction of the splitting variables, g is only contained in quadratic terms and thus is easily solved for

$$\begin{array}{cc}\hfill g^{k+1}=\underset{g}{arg\; min}& {\mathcal{L}}_{\alpha ,\lambda }(b_i^{k+1},h^{k+1},g,q_i^k,r^k)\hfill \\ & =\underset{g}{arg\; min}\hfill & \left\{\frac{\alpha }{2}\sum _i\parallel b_i-D(g-f_i)+\frac{q_i}{\alpha }{\parallel }_2^2+\frac{\lambda \alpha }{2}{\parallel h-g+\frac{r}{\alpha }\parallel }_2^2\right\}\Rightarrow \hfill \end{array}$$

$$\begin{array}{c}\hfill \frac{\delta \mathcal{L}}{\delta g}=\alpha \sum _i-D^T(b_i-D(g-f_i)+\frac{q_i}{\alpha })-\lambda \alpha (h-g+\frac{r}{\alpha })=\\ \hfill \alpha (C{D}^{T}D+\lambda I)g-\alpha D^T(\sum _i{b}_i+D{f}_i+\frac{q_i}{\alpha })-\lambda \alpha h+\lambda r=0\Rightarrow \\ \hfill g={(C{D}^{T}D+\lambda I)}^{-1}\left(D^T(\sum _i{b}_i+D{f}_i+\frac{q_i}{\alpha })+\lambda (h-\frac{r}{\alpha })\right)\end{array}$$

The Langrangian multipliers (split Bregman variables) are updated through gradient ascent:

$$\begin{array}{cc}\hfill q_i^{k+1}& =q_i^k+\tau \alpha (b_i^{k+1}-D(g^{k+1}-f_i))\hfill \\ \hfill r^{k+1}& =r^k+\tau \lambda \alpha (h^{k+1}-g^{k+1}))\hfill \end{array}$$

Isotropic: As a result of the coupling of the terms in this version of the minimization problem, we cannot compactly write the problem with matrices as above; however, the solution is just as readily available. Here the ⋄ denotes the Hadamard matrix power operator, which acts pointwise on the matrix.

Point Form:

$$\underset{\overrightarrow{g}}{min}\left\{\sum _{i,j}\sqrt{{({\delta }_x{f}_{i,j})}^2+{({\delta }_y(f_{i,j}-g_j))}^2}+\lambda {\parallel g(y)\parallel }_1\right\}$$

Matrix Form:

$$\underset{\overrightarrow{g}}{min}\left\{\parallel {[{(f{D}^T)}^{\diamond 2}+{(D(f-\overrightarrow{g}\otimes \overrightarrow{1_C}))}^{\diamond 2}]}^{\diamond \frac{1}{2}}{\parallel }_{1,1}+\lambda {\parallel g(y)\parallel }_1\right\}$$

Just as before, we introduce splitting variables and Lagrangian multipliers to form the augmented Lagrangian:

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\alpha ,\lambda }(a_{i,j},b_{i,j},h,g,p_{i,j},q_{i,j},r)& =& \sum _{i,j}\sqrt{|a_{i,j}{|}^2+{|b_{i,j}|}^2}+\frac{\alpha }{2}{\parallel a_{i,j}-{\delta }_x{f}_{i,j}\parallel }_2^2+\langle p_{i,j},a_{i,j}-{\delta }_x{f}_{i,j}\rangle \hfill \\ & & +\frac{\alpha }{2}{\parallel b_{i,j}-{\delta }_y(f_{i,j}-g_j)\parallel }_2^2+\langle q_{i,j},b_{i,j}-{\delta }_y(f_{i,j}-g_j)\rangle \hfill \\ & & +\lambda (\parallel h{\parallel }_1+\frac{\alpha }{2}\parallel h-g{\parallel }_2^2+\langle r,h-g\rangle )\hfill \\ & =& \sum _{i,j}\sqrt{|a_{i,j}{|}^2+{|b_{i,j}|}^2}+\frac{\alpha }{2}{\parallel a_{i,j}-{\delta }_x{f}_{i,j}+\frac{p_{i,j}}{\alpha }\parallel }_2^2\hfill \\ & & +\frac{\alpha }{2}{\parallel b_{i,j}-{\delta }_y(f_{i,j}-g_j)+\frac{q_{i,j}}{\alpha }\parallel }_2^2\hfill \\ & & +\lambda \left({\parallel h\parallel }_1+\frac{\alpha }{2}{\parallel h-g+\frac{r}{\alpha }\parallel }_2^2\right)\hfill \\ & & +\mathcal{O}(q_{i,j}^2,p_{i,j}^2,r^2)\hfill \end{array}$$

The splitting variables $a_{i,j}$ and $b_{i,j}$ are updated by the vectorial proximal mapping operator:

$$\begin{array}{ccc}\hfill \langle a_{i,j},b_{i,j}\rangle & =& \underset{\langle a_{i,j},b_{i,j}\rangle }{arg\; min}{\mathcal{L}}_{\alpha ,\lambda }(a_{i,j},b_{i,j},h^k,g^k,p_i^k,q^k)\hfill \\ & =& \underset{\langle a_{i,j},b_{i,j}\rangle }{arg\; min}\left\{\sqrt{|a_{i,j}{|}^2+{|b_{i,j}|}^2}+\frac{\alpha }{2}{\parallel a_{i,j}-{\delta }_x{f}_{i,j}+\frac{p_{i,j}}{\alpha }\parallel }_2^2\right.\hfill \\ & & +\left.\frac{\alpha }{2}{\parallel b_{i,j}-{\delta }_y(f_{i,j}-g_j)+\frac{q_{i,j}}{\alpha }\parallel }_2^2\right\}\hfill \\ & =& \underset{\langle a_{i,j},b_{i,j}\rangle }{arg\; min}\{\parallel \langle a_{i,j},b_{i,j}\rangle \parallel +\left.\frac{\alpha }{2}{\parallel \langle a_{i,j},b_{i,j}\rangle -\langle {\delta }_x{f}_{i,j}+\frac{p_{i,j}}{\alpha },{\delta }_y(f_{i,j}-g_j)+\frac{q_{i,j}}{\alpha }\rangle \parallel }^2\right\}\hfill \\ & =& S_{\frac{1}{\alpha }}(\langle {\delta }_x{f}_{i,j}+\frac{p_{i,j}}{\alpha },{\delta }_y(f_{i,j}-g_j)+\frac{q_{i,j}}{\alpha }\rangle )\hfill \end{array}$$

Each component of the vector is updated via shrinkage as follows:

$$a_{i,j}=\frac{{\delta }_x{f}_{i,j}+\frac{p_{i,j}}{\alpha }}{s}·max(s-\frac{1}{\alpha },0)$$

$$b_{i,j}=\frac{{\delta }_y(f_{i,j}-g_j)+\frac{q_{i,j}}{\alpha }}{s}·max(s-\frac{1}{\alpha },0)$$

$$s=\sqrt{{({\delta }_x{f}_{i,j}+\frac{p_{i,j}}{\alpha })}^2+{({\delta }_y(f_{i,j}-g_j)+\frac{q_{i,j}}{\alpha })}^2}$$

The splitting variable h, the striping mask g, and the Lagrangian multipliers are updated as before due to the common structure between the two models:

$$\begin{array}{cc}\hfill p_{i,j}^{k+1}& =p_{i,j}^k+\tau \alpha ({\delta }_x{f}_{i,j}+\frac{p_{i,j}}{\alpha }))\hfill \end{array}$$

We now describe the algorithm for the TV–$L^1$ model (see Algorithm 1). The first algorithm (anisotropic) is presented in vector form. The second algorithm (isotropic) is presented in matrix form.

2.5.3. Algorithm

Algorithm 1: Complete ADMM optimization of TV–$L^1$.

1: Initialize: $A^0=({\overrightarrow{a}}_i)$, $B^0=({\overrightarrow{b}}_i)$, $Q^0=({\overrightarrow{q}}_i)\leftarrow 0\in {\mathbb{R}}^{R\times C}$ , $\overrightarrow{g^0}$, $\overrightarrow{h^0}$, $\overrightarrow{r^0}\leftarrow 0\in {\mathbb{R}}^{R\times 1}$, $n\leftarrow 0\in \mathbb{R}$ 2: $f=(\overrightarrow{f_i})\leftarrow ln(F)$, $D\leftarrow 0\in {\mathbb{R}}^{R\times R}$, ${(D)}_{i,i}=-1$, ${(D)}_{i,i+1}=1$ for $i=1,\dots ,R-1$ 3: repeat 4:     $n\leftarrow n+1$ 5:     case Anisotropic: 6:         for $i=1:C$ do 7:            Update splitting variable for smoothness regularizer term via shrinkage: $${\overrightarrow{b_i}}^{n+1}\leftarrow S_{\frac{1}{\alpha }}\left(D({\overrightarrow{g}}^n-{\overrightarrow{f}}_i)-\frac{{\overrightarrow{q}}_i^n}{\alpha }\right)$$ 8:            Update Lagrangian multiplier for regularizer term via dual ascent: $${\overrightarrow{q_i}}^{n+1}\leftarrow {\overrightarrow{q_i}}^n+\tau \alpha ({\overrightarrow{b_i}}^{n+1}-D({\overrightarrow{g}}^n-\overrightarrow{f_i}))$$ 9:         end for 10:     case Isotropic: 11:         for $i=1:C$, $j=1:R$ do 12:            Update splitting variables for smoothness regularizer term via vectorial shrinkage: $$\langle a_{i,j}^{n+1},b_{i,j}^{n+1}\rangle \leftarrow S_{\frac{1}{\alpha }}(\langle {\delta }_x{f}_{i,j}+\frac{p_{i,j}^n}{\alpha },{\delta }_y(f_{i,j}-g_j^n)+\frac{q_{i,j}^n}{\alpha }\rangle )$$ 13:            Update Lagrangian multipliers for regularizer term via dual ascent: $$p_{i,j}^{n+1}\leftarrow p_{i,j}^n+\tau \alpha ({\delta }_x{f}_{i,j}+\frac{p_{i,j}}{\alpha }))$$ $$q_{i,j}^{n+1}\leftarrow q_{i,j}^n+\tau \alpha (b_{i,j}^{n+1}-{\delta }_y(g_j^n-f_{i,j}))$$ 14:         end for 15:     Update splitting variable for data fidelity term via shrinkage: $${\overrightarrow{h}}^{n+1}\leftarrow S_{\frac{1}{\alpha }}\left({\overrightarrow{g}}^n-\frac{{\overrightarrow{r}}^n}{\alpha }\right)$$ 16:     Update striping mask: $${\overrightarrow{g}}^{n+1}\leftarrow {(C{D}^{T}D+\lambda I)}^{-1}\left(D^T(\sum _i{\overrightarrow{b_i}}^{n+1}+D\overrightarrow{f_i}+\frac{{\overrightarrow{q_i}}^{n+1}}{\alpha })+\lambda {\overrightarrow{h}}^{n+1}-\frac{\lambda }{\alpha }{\overrightarrow{r}}^n\right)$$ 17:     Update Lagrangian multiplier for data fidelity term via dual ascent: $${\overrightarrow{r}}^{n+1}\leftarrow {\overrightarrow{r}}^n+\tau \lambda \alpha ({\overrightarrow{h}}^{n+1}-{\overrightarrow{g}}^{n+1}))$$ 18:     Update energy terms: $$E_1^{n+1}\leftarrow \sum _i\parallel D(\overrightarrow{f_i}-{\overrightarrow{g}}^{n+1}){\parallel }_1,E_2^{n+1}\leftarrow \lambda {\parallel {\overrightarrow{g}}^{n+1}\parallel }_1$$ $$E^{n+1}\leftarrow E_1^{n+1}+E_2^{n+1}$$ 19: until convergence: $$\parallel {\overrightarrow{g}}^{n+1}-{\overrightarrow{g}}^n{\parallel }_2^2/\parallel {\overrightarrow{g}}^n{\parallel }_2^2<{ϵ}_g\mathrm{and}|E^{n+1}-E^n{|}^2/{|E^n|}^2<{ϵ}_E$$ 20: Retrieve clean image: $$u\leftarrow f-g^{n+1}\otimes [1_1,1_2,...,1_C]$$ $$U\leftarrow exp(u)$$

3. Results

In our experiments, we used data acquired by the AirMSPI. The AirMSPI is an airborne prototype instrument similar to that of the future satellite-borne MSPI instrument for obtaining multi-angle polarization imagery [29]. The instrument was built for NASA by the Jet Propulsion Laboratory in Pasadena, California and has been flying aboard the NASA ER-2 high altitude aircraft since October 2010.

The AirMSPI is an eight-band (355, 380, 445, 470, 555, 660, 865, and 935 nm) pushbroom camera, measuring polarization in the 470, 660, and 865 nm bands, mounted on a gimbal to acquire multi-angular observations over a $\pm {67}^{\circ }$ along-track range. Two principal observing modes are employed: step-and-stare, in which 11 km × 11 km targets are observed at a discrete set of view angles with a spatial resolution of ∼10 m; and continuous sweep, in which the camera slews back and forth along the flight track within $\pm {67}^{\circ }$ to acquire wide-area coverage (11 km swath at nadir; target length of 108 km) with a ∼25 m spatial resolution. Step-and-stare provides more angles, but continuous sweep gives greater coverage. Multiple observing modes can be programmed into the instrument and activated under cockpit control. Multi-angle radiance and polarization imagery from the AirMSPI will provide 3D scene context where clouds and aerosol plumes are present. It will also enable the retrieval of aerosol and cloud macrophysical properties (distribution and height), microphysical properties (size distribution, single scattering albedo, and shape), and optical depth.

We first compare destriping results generated using TV–$L^2$ and the proposed TV–$L^1$ models. Figure 1 shows the 355 nm ultraviolet (UV) channel image with stripes captured by the AirMSPI instrument from the nadir angle at Mojave, California. The image has been destriped using the TV–$L^2$ and TV–$L^1$ models. As we can see from the destriped images and corresponding differences between the captured and destriped images, the TV–$L^2$ model does not preserve radiometric intensities in regions where no stripes are present. Figure 2 (left) shows plots of the recovered function G for the TV–$L^2$ and TV–$L^1$ destriped images. The TV–$L^1$ recovered function G is closer to the identity, particularly for the rows with no stripes, suggesting the TV–$L^1$ reconstruction is more accurate than the TV–$L^2$ reconstruction. Figure 2 (right) shows plots of sums over all rows of the original image with stripes (from Figure 1), as well as sums of rows for the TV–$L^2$ and TV–$L^1$ destriped images. These plots indicate that the TV–$L^1$ model preserves radiometric intensities of the original images better, and the TV–$L^2$ model produces more artificial smoothing throughout the image. Figure 3 shows histograms of the function G for the TV–$L^2$ and TV–$L^1$ reconstructions. The TV–$L^1$ reconstruction is pointier than the TV–$L^2$ reconstruction. It is also centered at 1, as opposed to the TV–$L^2$ reconstruction, which further indicates better accuracy of the TV–$L^1$ model. We note that the true stripes, at around $G\approx 0.95$, are represented in the histograms by small bumps.

Figure 4, Figure 5 and Figure 6 display similar results to those shown in Figure 1, Figure 2 and Figure 3 for the 355 nm channel image with stripes depicting clouds over the Pacific Ocean captured by the AirMSPI instrument from the ${66.0}^{\circ }$F angle.

Because the TV–$L^1$ model generates results with better accuracy than the TV–$L^2$ model, we focus on the TV–$L^1$ model in the following examples. Figure 7, Figure 8 and Figure 9 show more examples of the TV–$L^1$ reconstruction of images captured using the continuous-sweep observing mode. Figure 10 and Figure 11 display images captured using the step-and-stare observing mode as well as destriped results using the TV–$L^1$ model.

4. Discussion

We have provided detailed derivations and a motivated evolution of the optimization problem so that the proposed TV–$L^1$ method can be accessible to all academic disciplines involved with image processing. We have presented a variational model that is well-defined, qualitatively motivated, and easily minimized. The constructed energy uses sparsity-promoting energy functionals, on the basis of TV and $L^1$ energy, to achieve minimally invasive destriping. Both isotropic and anisotropic TV, along with $L^1$ energy, are considered in our variational model. The ADMM is used in conjunction with nonlinear proximal operators to efficiently optimize the energy, yielding quick and high-quality results.

The $L^2$ fidelity overly penalizes stripes of large magnitude and likewise underexaggerates the significance of stripes of small magnitude. In areas of no striping, we intend our striping mask to be very close to the identity, while in areas of heavy striping, we wish to remove such striping and thus yield a larger magnitude of our striping mask in that region. The $L^1$ fidelity gives a smaller striping mask in areas of no striping, leaving enough energy to remove the heavier striping in localized areas of the image. Because there is no prior knowledge of the distribution of the stripes and qualitatively we may wish to remove deep striping effects while preserving sharp geometry, we conclude that the model with an $L^1$ striping penalty is more effective.

The TV–$L^1$ model preserves radiometric intensities of the original images better, and the TV–$L^2$ model produces more artificial smoothing throughout the image. TV–$L^1$ recovers the function G that is closer to the identity, particularly for the rows with no stripes, suggesting that the TV–$L^1$ reconstruction is more accurate than the TV–$L^2$ reconstruction. The histograms of the function G for the TV–$L^1$ reconstruction are more compact around 1 than those for the TV–$L^2$ reconstruction. They are also centered at 1, as opposed to the TV–$L^2$ reconstruction, which further indicates the better accuracy of the TV–$L^1$ model.

5. Conclusions

In this paper, we have presented a novel variational method for image destriping through fast minimization techniques of appropriately modeled energy functionals, namely TV and the $L^1$ data fidelity term. Our destriping model solves the inverse problem as follows: it minimally removes a univariate multiplicative striping mask from the data, such that the clean image is somewhat smooth and the removed stripe has low energy. We assess the smoothness of the clean image using the TV, which maintains sharp image features and preserves definition and contrast. We address both isotropic and anisotropic TV in this paper, each having their respective strengths and weaknesses. We use $L^1$ and, for comparison, $L^2$ energy to measure the removed striping, ensuring minimal data removal and thus a clean image of high fidelity.

The variational problem is solved efficiently using an ADMM approach: it introduces splitting variables and quadratic penalties for deviations from the splitting variables to allow for efficient optimization via proximal shrinkage operators, explicit quadratic solutions, and simple gradient ascent for the Lagrangian multipliers. In our experiments, we have shown that the proposed TV–$L^1$ method yields qualitatively good results and removes minimal masking, and it does so quickly in terms of both iterations and time. From the histogram distributions of G, we observe a narrower spread around 1, yet a wider, more equidistributed support, suggesting that most of the time, there is minimal masking removal (multiplier close to 1); however, in areas of heavy striping, the destriping effect is more prevalent and is of greater magnitude.

Applications of this algorithm are not limited to satellite imagery, and they may be analogized to other fields, such as raster scans in microscopy. Any scientific measurements (of images) made mostly along a curve—parameterizable by a single dimension—may be susceptible to such striping biases and may be a candidate for similar destriping. Future work will expand this model to multi-modal images and color images, and it may incorporate other specific priors on the data.