Development of a Hybrid Solution to Replacement of Clouds and Shadows in Remote Sensing Images

and (cloud parameter and shadow parameter, in the order given), increasing the algorithm flexibility. The equation that describes this operation is:

$$$\begin{aligned} m(x,y) =\left\{ \begin{array}{ll} f(x,y) < sp \times f_{av-sd}, &{} f(x,y) \varepsilon 0, \\ f_{av-sd} < f(x,y) < f_{av}, &{} f(x,y) \varepsilon 1, \\ f_{m}< f(x,y) < cp \times f_{av+sd}, &{} f(x,y) \varepsilon 2, \\ f(x,y) > cp \times f_{av+sd}, &{} f(x,y) \varepsilon 3. \\ \end{array} \right. \end{aligned}$$” src=”/wp-content/uploads/2016/03/A320009_1_En_15_Chapter_Equ1.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(1)</DIV></DIV>where <SPAN id=IEq3 class=InlineEquation><IMG alt=$$f(x,y)$$ src=$ is the pixel value, $f_{av}$ is the average value of image pixels, $f_{av+sd}$ is the sum of the average value pixels and the standard deviation value of the image and $f_{av-sd}$ is the subtraction of the average value by the standard deviation pixels.

The region labeled as $$0$$

represents a shadow region, labeled as $$1$$

means region not affected by atmospheric interference, while regions labeled as $$2$$

represents thin clouds, and dense cloud are labeled as $$3$$

. For images with multiple bands these labels are assigned if and only if the rule is valid for all bands.

To complete this process, it is applied a morphological opening operation that aims to remove very small objects that can cause mistakes in following steps.

2.2 Image Decomposition

The method proposed by Vese and Osher [18] decomposes an image into two sub-images, each representing the components of structure or texture, thus, a better image redefinition can be made. On the structured part, should be applied the technique of inpainting based on DCT, whilst into texture portions and heterogeneous areas is suitable to use texture synthesis.

The generalized model is defined as:

$\begin{aligned} f = u + v, \end{aligned}$

(2)

where

is the input image, $$u$$

is the structure image and $$v$$

is the texture image. So, given these sub-images one can reconstruct the original image. However, in practice it is observed that the original image can be only approximately reconstructed. The goal of the method is to have a structure image $$u$$

that preserves all strong edges with smoothed internal regions, and an image $$v$$

that contains all the texture and noise information.

The method used to construct the structured image $$u$$

is based on the assumption that $$u$$

is a

function and in attempt to minimize this function in the space of all Bounded Variation functions (BV). Functions in BV space are functions whose total variation are limited by some constant value less than infinity. Minimizing $$u$$

in BV space ensures resulting in a stable image and without infinite values. It should be noted, however, that this space allows for functions which have very large derivatives (although non-infinite), thereby ensuring that strong edges are preserved.

Taking in mind the intuition described above, the minimization problem should logically have two terms. One of them will be the fidelity, responsible for maintaining the difference between $$f$$

and

small. This term ensure that data of the input image are kept on result. The other one imply a smoothing over $$u$$

, although not necessarily in all $$u$$

components. The minimization is computed as [14]:

$\begin{aligned} F(u) = \left[ (\int {|\nabla u|}) + (\lambda \int {|f - u|^2})\right] dxdy, \end{aligned}$

(3)

with $u \in BV$ and $\nabla$ representing the gradient operator. The second term is the data term and the first one is a regularization term to ensure a relatively smooth image. $\lambda$ is a tuning parameter. As can be seen, this seeks find the optimal $$u$$

and ignores the $$v$$

image. The reason is that, in Vese and Osher [18], the authors had considered the $$v$$

image to be noise, and therefore to be discarded.

There is an unique result to this optimization problem, and methods exist for finding the solution. Noting that $$v = f - u$$

it is possible to easily modify the above equation to incorporate $$v$$

$\begin{aligned} F(u) = \int {|\nabla u|} + \lambda \int {\parallel v\parallel ^2} dxdy, \end{aligned}$

(4)

(still $u \in BV$ ). Which yields the Euler-Lagrange equation $u = f + \frac{1}{2\lambda }div(\frac{\nabla u}{|\nabla u|})$ . Making the right manipulation $v = f - u = -\frac{1}{2\lambda }div(\frac{\nabla u}{|\nabla u|})$ . At this point it is useful to break $$v$$

into its

and

components respectively. It will be denoted as $$g_1$$

and

, where:

$\begin{aligned} g_1 = -\frac{1}{2\lambda }div(\frac{\nabla u_x}{|\nabla u|})\ \text {and}\ g_2 = -\frac{1}{2\lambda }div(\frac{\nabla u_y}{|\nabla u|}). \end{aligned}$

(5)

This allows us to write $$v$$

as: $v = div \varvec{g}$ where $\varvec{g} = (g_1,g_2)$ . It can be seen that $g_1^2 + g_2^2 = \frac{1}{2\lambda }$ , so that $\parallel \sqrt{g_1^2 + g_2^2}\parallel = \frac{1}{2\lambda }$ .

This allows us to rewrite $$v$$

as:

$\begin{aligned} v(x,y) = div \varvec{g} = {\partial }_x g_1(x,y) + \partial _y g_2(x,y). \end{aligned}$

(6)

And leads to the final minimization problem ( $u \in BV$ ):

$\begin{aligned} G(u,v1,v2)&=\displaystyle \int {|\nabla u|} + \lambda \int {|f - u - \partial _x g_1} \nonumber \\&\quad - \partial _y g_2|^2 dxdy + \mu \displaystyle \int {\sqrt{g_1^2 + g_w^2}dxdy}. \end{aligned}$

(7)

Solving the minimization problem (Eq. 7) yields the Euler-Lagrange equations:

$\begin{aligned} u&= f - \partial _x g_1 - \partial _y g_2 + \frac{1}{2\lambda }div(\frac{\nabla u}{|\nabla u|}) \end{aligned}$

(8)

$\begin{aligned} \mu \frac{g_1}{\sqrt{g_1^2 + g_2^2}}&= 2\lambda [\frac{\partial }{\partial _x}(u-f) + \partial _{xx}^2 g_1 + \partial _{xy}^2 g_2] \end{aligned}$

(9)

$\begin{aligned} \mu \frac{g_2}{\sqrt{g_1^2 + g_2^2}}&= 2\lambda [\frac{\partial }{\partial _y}(u-f) + \partial _{xy}^2 g_1 + \partial _{yy}^2 g_2] \end{aligned}$

(10)

2.3 Inpainting by Smoothing Based on Multidimensional DCT

This method was proposed by Garcia [7], and so as in Bertalmio et al. [1], is based on the information propagation by smoothing. The specificity of this approach is related to the use of the Discrete Cosine Transform (DCT) to simplify and to solve linear systems, to an efficient smoothing.

2.3.1 Smoothing by Penalized Least Squares Regression

In statistics and data analysis, smoothing is used to reduce experimental noise or information and keeping the most important marks of the data set. Considering the following model for the one-dimensional noisy signal $$y$$

from the Eq. 11.

$\begin{aligned} y = \hat{y} + \varepsilon , \end{aligned}$

(11)

where $\varepsilon$ represents a Gaussian noise with zero mean and unknown variance, and $\hat{y}$ is the so-called smoothing, i.e., has continuous derivatives up to some order (usually $\ge$ 2) throughout the domain. The smoothing of $$y$$

depends on the best estimate of $\hat{y}$ and this operation is usually performed by a parametric or nonparametric regression.

A classic approach to smooth is the Penalized Least Squares Regression. This technique minimizes a criterion that balances the data fidelity, measured by the Residual Sum-of-Squares ( $$RSS$$

) and by a penalty term ( $$P$$

), which reflects the robustness of the smoothed data. Another simple and straightforward approach to express the robustness is by using a Second-order Divided Difference ( $$SoDD$$

), which produces an one-dimensional array of data.

Now, using $$RSS$$

and the

, the minimization of $F(\hat{y})$ results in a linear system, expressed in Eq. 12, which allows the smoothed data determination.

$\begin{aligned} (I_n + s D^{T}D)\hat{y} = y, \end{aligned}$

(12)

where

is the identity matrix $n \times n$ , $$s$$

is a positive real scalar that controls the grade of smoothing, so that, as it increases, the degree of smoothing of $\hat{y}$ increases too; and $$D^T$$

represents the transpose of $$D$$

. Its important to note that $(I_n + s D^{T}D)$ is a penta-diagonal symmetric matrix, and the last equation can be solved numerically in a computationally efficient way.

2.3.2 Smoothing Equally Spaced Data

Equation 12 can be solved using the left division matrix applied to sparse matrices [7]. Solving this linear system, however, can be a lot of time expensive for a large amount of data. But, this algorithm can be simplified and accelerated, since the data are evenly spaced, in images where pixels are equally spaced, resulting in the following equation for multidimensional data:

$\begin{aligned} \hat{y} = IDCT(\varGamma DCT(y)), \end{aligned}$