A Bandelet-Based Inpainting Technique for Clouds ... - IEEE Xplore

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 47, NO. 7, JULY 2009

2363

A Bandelet-Based Inpainting Technique for Clouds Removal From Remotely Sensed Images Aldo Maalouf, Philippe Carré, Bertrand Augereau, and Christine Fernandez-Maloigne

Abstract—It is well known that removing cloud-contaminated portions of a remotely sensed image and then filling in the missing data represent an important photo editing cumbersome task. In this paper, an efficient inpainting technique for the reconstruction of areas obscured by clouds or cloud shadows in remotely sensed images is presented. This technique is based on the Bandelet transform and the multiscale geometrical grouping. It consists of two steps. In the first step, the curves of geometric flow of different zones of the image are determined by using the Bandelet transform with multiscale grouping. This step allows an efficient representation of the multiscale geometry of the image’s structures. Having well represented this geometry, the information inside the cloud-contaminated zone is synthesized by propagating the geometrical flow curves inside that zone. This step is accomplished by minimizing a functional whose role is to reconstruct the missing or cloud contaminated zone independently of the size and topology of the inpainting domain. The proposed technique is illustrated with some examples on processing aerial images. The obtained results are compared with those obtained by other clouds removal techniques. Index Terms—Image wavelets.

reconstruction,

image

restoration,

I. I NTRODUCTION

O

NE OF THE major limitations of remote-sensing imaging could be its sensitivity to weather conditions during the image acquisition process. The resulting images are frequently contaminated by clouds and cloud shadows. Remote sensing has to cope with the so-called cloud removal problem which is an important difficulty affecting the observation of the Earth surface. Diverse techniques have been proposed including different methods to solve this problem. They can be grouped in two categories. The first category regroups the measurementbased approaches. These approaches consist of selecting the best measurement (the most cloud-free pixel) among a set of measurements acquired over a limited time period to represent the considered multitemporal pixel over that time period. Examples of these approaches are the recent work of Melagani and Benabdelkader [2], [6], [18] and the work of Liew et al. [10]. In [6], the author proposes a contextual prediction process in order to determine the spectrotemporal relationships between Manuscript received February 5, 2008; revised June 5, 2008 and October 14, 2008. First published March 24, 2009; current version published June 19, 2009. The authors are with the XLIM Laboratory, Unite Mixte de Recherche 6172, Centre National de la Recherche Scientifique, Signal-Image-Communication Department, University of Poitiers, 86034 Poitiers Cedex, France (e-mail: [email protected]; [email protected]; [email protected]; [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2008.2010454

the set of acquired images. These spectrotemporal relationships are deduced from cloud-free areas in the spatial neighborhood of the cloud-contaminated region over the available series of the temporal images. The contextual prediction process is implemented using linear predictors and support vector machines. The major drawback of this approach is the assumption that the spatial structure of the image is identical over the image sequence. For that, Benabdelkader et al. proposed in [18] and [2] to use a postreconstruction methodology that is based on a contextual spatiospectral postprediction system. The role of this postreconstruction is to incorporate the spatial information by taking advantage from the local properties between pixels in a predefined neighborhood system. However, the major limitation of this approach is that the acquisition dates should be close to each other, and the spatial dynamics of the geographical area under analysis should be slow compared to the total time interval of the sequence. In the other hand, the adaptive reconstruction technique proposed by Lee and Crawford [8] avoids this kind of problem. They assumed that the temporal signature of a given pixel is contaminated by residual effects caused by imperfect sensing of the target and by spatially autocorrelated noise due to atmospheric attenuation. However, this method showed a high computational complexity and is not easily applicable to the general case of nonstationary temporal image series. Another method was presented by Remund et al. [5] for recovering Advanced Very High Resolution Radiometer measurements that are modified by the effects not only of clouds but also of cloud shadows. This method is simple and effective, but presents the drawback of being limited to data acquired over vegetated areas. The second category regroups cloud-removal techniques that consist of filling in the clouds-contaminated region using traditional synthesis and image inpainting techniques. An example of these approaches is the recent work of Peng et al. [4]. The goal of the approaches of this category is to seamlessly synthesize a complete, visually plausible and coherent image. The main drawback of these methods is that most of the inpainting techniques are applicable for small-scale flaws like scratches or stains. Recently, there has been a major work for improving the inpainting techniques in order to achieve long region connection and texture recovering. An example of the inpainting techniques are the ones presented in [3], [9], [12], [19], and [20]. In this paper, the strategy of the approaches of the second category is pursued, and a new efficient inpainting technique for missing data synthesis is presented. The advantage of this reconstruction technique is that it is capable of filling large regions and synthesizing even highly textured missing regions. Thus, the drawbacks of the measurement-based approaches

0196-2892/$25.00 © 2009 IEEE

2364


as well as the limitations of the inpainting techniques of the approaches of the second category are circumvented. Therefore, an efficient cloud-removing scheme is achieved. The proposed reconstruction technique is based on the Bandelet transform and the multiscale geometrical grouping. It consists of two steps. In the first step, the curves of geometric flow of different zones of the image are determined by using the Bandelet transform with multiscale grouping. This step allows a better representation of the multiscale geometry of the image’s structures than other geometrical representations [11]. The Bandelet transform with multiscale grouping is obtained by applying a bandeletization process on the wavelet coefficients at each scale and orientation. This bandeletization with Multiscale grouping offer more flexibility in the sense that we can define geometry along flow lines that are not parallel. This geometrical representation converges toward singularities. Therefore, fine structures are well represented, and the propagation of the represented geometrical data in the cloud-contaminated region may guarantee an efficient reconstruction of the hidden data. Having well represented this geometry, the information inside the cloudcontaminated zone is synthesized by propagating the geometrical flow curves inside that zone. The latter step is accomplished by minimizing a functional whose role is to reconstruct the missing or damaged zone independently of the size and topology of the reconstruction or inpainting domain. As a result, the flow lines are well tied inside the cloud-contaminated zone, and the missing data are accurately synthesized. The rest of this paper is organized as follows. In Section II, a review of the bandelet transform and multiscale grouping is given. In Section III, the reconstruction technique is described. Section IV is devoted for experimental results, and Section V is a conclusion. II. R EVIEW OF B ANDELET T RANSFORM AND M ULTISCALE G ROUPING In this section, the bandelet transform with multiscale grouping is described. This transform is then used to obtain a geometrical flow of the image. The propagation of this flow into the cloud-contaminated regions synthesizes the missing data. A. Orthogonal Bandelets Approximations A sparse representation takes advantage of some kind of regularity of the function to approximate. Geometric regularity along edges in images is an anisotropic regularity. Although the image may be discontinuous across a contour, the image can be differentiable in a direction parallel to the tangent of the edge curve. The Bandelet transform exploits such an anisotropic regularity by constructing orthogonal vectors that are elongated in the direction where the function has a maximum of regularity. The first bandelet bases constructed by Le Pennec and Mallat [13], [14] have bring optimal approximation results for geometrically regular functions. Later works have enriched this construction by the use of multiscale geometry defined over the coefficients of a wavelet basis [16], [17]. These multiscale bandelet bases are described in the following section.

Fig. 1. (a) Wavelet coefficients and geometric flow. (b) Sampling position and geometric flow. (c) Warped sampling.

Fig. 2.

Example of quadtree segmentation.

B. Basis of Orthonormal Bandelets The bandelet approximation [17] is obtained by computing the polynomial approximation by a thresholding in an orthogonal Alpert basis [1]. The Alpert transform can be thought as a warped wavelet transform adapted to an irregular sampling, as shown in Fig. 1. It is obtained by orthogonalizing multiresolution space of polynomials defined on the irregular sampling grid. The bandeletization of wavelet coefficients using an Alpert transform defines a set of bandelet coefficients. These coefficients can be written as inner products f, bkj,l,m of the original image f with bandelet functions that are linear combinations of wavelet functions k bkj,l,n (x) = al,n [p]ψj,p (x). (1) p

The al,n [p] are the coefficients of the Alpert transform, which depends on the local geometric flow, since this flow defines the warped sampling locations shown in Fig. 1(c). The bandelet function is defined at location n and scale 2j . The indice k designates the wavelet orientation. The Alpert transform introduces a new scale factor 2l > 2j which defines the elongation of the bandelet function. The bandelet bkj,l,n (x) inherits the regularity of the wavelet ψjk (x). C. Segmented Geometric Flow The family of orthogonal bandelets depends on the local adapted flow defined for each scale 2j and orientation k. This parallel flow is characterized by an integral curve γ shown in red in the Fig. 1(a) (the dashed curve). This integral curve does not need to be strictly parallel to the contour. This is due to the bidimensional regularization introduced by the smoothing of fj = f ∗ ψjk with the wavelet ψjk . In order to approximate the geometry by a polynomial flow, we need to segment the set of wavelet coefficients in square S. For each scale 2j and orientation k of the wavelet transform, this segmentation is obtained using a recursive subdivision in dyadic squares of various sizes, as shown in Fig. 2.

MAALOUF et al.: INPAINTING TECHNIQUE FOR CLOUDS REMOVAL FROM REMOTELY SENSED IMAGES

Such a division defines a quadtree that specifies if a square S should be further subdivided in four subsquares of size twice smaller. If there is no specific direction of regularity inside a square, which is the case either in uniformly regular regions or at the vicinity of edge junctions, then there is no geometric directional regularity to exploit. Thus, it is not necessary to modify the wavelet basis. The geometry of orthogonal bandelets is described by a locally parallel flow over each square of the dyadic segmentation. Such a geometry is suitable for the approximation of the geometrically regular images, but lacks flexibility to represent the complex geometry of turbulent textures of remotely sensed images. Junctions are not explicitly modeled and require a fine recursive segmentation to be isolated from the remaining contours. Furthermore, the segmentation in small square areas forbid to take advantage of the long range regularity of fine elongated structures such as the boundaries of different terrains being imaged. To overcome these problems, our attention is directed toward the use of multiscale grouping which was introduced in [15] and is summarized in the next section. D. Bandelet Transform With Multiscale Grouping The multiscale grouping uses a multiscale association field in order to group together coefficients in the direction specified by the flow [11]. These recursive groupings allow one to take into account junctions and long range regularities of remotely sensed images. The multiscale grouping is first computed by applying the Haar transform over pairs of points that are neighbors according to an association field. The role of this field is to group together points that have similar neighborhoods in order to exploit the geometry of the signal. Then, a weighted mean and a weighted difference are computed between the values of the signal that are grouped together. The means and differences are consequently stored, and the process of computing the associated field and the Haar transform is repeated iteratively by doubling the scale at each step. This iterative process decomposes the original image in an orthogonal basis called grouping basis. The bandelet transform by multiscale grouping is finally obtained by applying the multiscale grouping to the set of coefficient of a multiscale representation. One applies the multiscale grouping over each scale 2j and orientation k of an orthogonal wavelet transform in order to get the decomposition of the image in an orthogonal basis of grouping bandelets. Similarly to the original orthogonal bandelet transform exposed in Section II-A, the cascade of the orthogonal wavelet transform and the orthogonal multiscale grouping defines an orthogonal bandelet transform by multiscale grouping. The bandelitization by grouping is more flexible than the bandelet transform since the association fields can deviate from the integral lines in order to converge to singularity points such as junction or crossings. Fine image structures are consequently well represented. Therefore, the propagation of the represented information in the cloud-contaminated regions of the image following the integral lines of the association field yields to a

2365

precise reconstruction of the missing data. We represent in the following section this reconstruction technique. III. R ECONSTRUCTION OF THE C LOUD -C ONTAMINATED R EGIONS The geometric reconstruction algorithms have motivated and inspired this paper mainly by the advantage of the geometric nature of the bandelet transform by multiscale grouping. Given that the image geometry is very well presented and characterized by a multiscale association field, we represent in this section a geometric reconstruction method based on the propagation of the geometric data from outside the cloudcontaminated region to fill-in the missing data. First, we represent some preliminaries definitions and then we describe our inpainting technique. A. Preliminaries Let x be an arbitrary point in the association geometric field. For every x ∈ Rn , we associate an n-dimensional vector space called the tangent space at x, which is denoted by Tx (Rn ). The use of the word “tangent” is motivated by the generalization to manifolds, for which the tangent spaces will be “tangent” to points on the manifold. The association geometric field is a none other a function that assigns a vector v ∈ Tx (Rn ) to every point x ∈ Rn . (x) at each point x ∈ Rn actually The association field V belongs to a different tangent space. The range of the function is therefore the union T (Rn ) = Tx (Rn ) (2) x∈Rn

which is called the tangent bundle on Rn . Even though the way we describe vectors from Tx (Rn ) may appear the same for any x ∈ Rn , each tangent space is assumed to produce distinct vectors. To maintain distinctness, a point in the tangent bundle can be expressed with 2n coordinates. The association field can therefore be expressed using n real-valued functions on Rn . Let fi (x1 , . . . , xn ) for i from 1 to n denote such functions. Using these, the association vector field is specified as f (x) = [f1 (x) f2 (x) . . .

fn (x)] .

(3)

In this case, it appears that a vector field is a function f from Rn to Rn . Therefore, standard function notation will be used from this point onward to denote a vector association field. Imagine a point that starts at some x0 ∈ Rn at time t = 0 and then moves according to the velocities expressed in the association field f . Its trajectory starting from x0 can be expressed as a function C : [0, ∞) → Rn , in which the domain is a time interval, [0, ∞). A trajectory represents an integral curve (or solution trajectory) of the differential equations with initial condition C(0) = x0 if dC (t) = f (C(t)) dt

(4)

2366


Fig. 3. (a) Integral lines of the association field. (b) Continuation of the integral lines of the flow.

for every time t ∈ [0, ∞). This could be expressed in integral form as t C(t) = x0 +

f (C(s)) ds

(5)

0

and is called a solution to the differential equations in the sense of Caratheodory. Intuitively, the integral curve starts at x0 and flows along the directions indicated by the velocity vectors [Fig. 3(a)]. The integral curve C is propagated into the cloudcontaminated region at a constant velocity. Let the image I be a real value function on a spatial domain Ω. We designate by ΩD the inpainting domain (ΩD ⊂ Ω). The domain ΩD regroups the pixels of the cloud-contaminated region, and ∂ΩD denotes the border of the inpainting zone [Fig. 3(b)]. The proposed model aims at recovering cloud areas of an image in such a way as to tie in the integral lines of the association field along those areas. To find the continuation of the integral lines, from outside to inside of the inpainting domain, following the direction of the association field, we propose a variational approach that is described in the next section. B. Reconstruction Technique

Fig. 4.

Grouping of an association field at scale 22 .

3) The filling in of information on the pixels border (belonging to ∂ΩD ) is performed in such a way as to satisfy the condition given by (6). In the numerical discretization, these conditions are reached by the following procedure. a) Compute the 2-D orthogonal wavelet transform of the image at a scale 2j . b) Compute the association field on the wavelet coefficients by first dividing each subband grid into two subgrids called “even subgrid” (Ωeven ) and “odd subgrid” (Ωodd ) and then associating to each wavelet coefficient modd in the odd subgrid a coefficient meven from the even subgrid that satisfies the so-called “best fit of radius P ” which is defined by meven = arg min m∈Ωeven

|a[m − n] − a[modd − n]|2 .

(7)

|n|