Special Section on Perceptually Driven Visual Information Analysis

Review of multiscale geometric decompositions in a remote sensing context

[+] Author Affiliations
Mariem Zaouali, Sonia Bouzidi, Ezzeddine Zagrouba

Université Tunis El Manar, Institut Supérieur de l’Informatique, Research Team SIIVA, Laboratoire LIMTIC, 2 Street Abou Rayhane Bayrouni, Ariana 2080, Tunisia

J. Electron. Imaging. 25(6), 061617 (Dec 02, 2016). doi:10.1117/1.JEI.25.6.061617
History: Received May 3, 2016; Accepted November 7, 2016
Text Size: A A A

Open Access Open Access

Abstract.  The quest for optimal representations is considered a challenging goal in the field of image processing. This consists of reducing the processing’s complexity while ensuring an efficient reconstruction. An optimal representation should conserve the properties of the image pertaining to smooth content and contours. The multiscale geometric decompositions (MGD) were designed to reach this finality. They were used in many fields and for different purposes, such as feature extraction, detail enhancing, and change detection. A state-of-art of these decompositions is proposed in this paper. We present their theoretical definitions and how they capture the feature of the objects within an image. An overview table is elaborated where we summarize the methods, the data and the different criteria of assessment used in the studied cases. We are interested, particularly, in the use of MGD in a remote sensing (RS) context. Thus, some examples of their applications on RS images are studied. A discussion is presented based on the analyzed cases.

Since 25 years ago, the definition of the wavelet has led to revolutionary achievements in the development of efficient encoding of piecewise regular signals. Their ability to detect singularities much more efficiently than Fourier transform1 has contributed to their remarkable success. The wavelet provides an optimal sparse approximation of a wide set of signals and digitalizes the continuous domain transform through the fast algorithmic implementation.

Nevertheless, this latter fails when dealing with singularities along curves, i.e., C2 functions, since reconstructing images with these anisotropic features could produce undesirable artifacts. A wavelet has only diagonal, vertical, and horizontal directions, and this limited directionality sensitivity leads to inaccurate results when dealing with multidimensional data. To overcome this drawback, some forms of directionalities have been introduced, such as the steerable pyramid,2 the directional filter banks (DFBs),3 and the two-dimensional (2-D) directional wavelet.4 A complexifyied version of wavelet was also proposed in Ref. 5 to tackle the directionality limitation by adding more precision to C2 functions. In fact, complex wavelet offers shift-invariance and yields, at each scale, six directional sub-bands oriented at +15 and 15  deg, +45 and 45  deg, +75 and 75  deg (indiscriminately by ordinary wavelet).

However, it was proven through several applications that these aforementioned methods cannot represent effectively the image singularities and its anisotropic features. A real transition occurred once Candès and Donoho6 introduced the second generation of curvelet (CT) in 2004. Since then, other multiscale geometric decompositions (MGD) were introduced to counteract the CT’s weakness, namely redundancy. In this work, we study some of these MGD and represent the characteristics behind their success. We are also interested in their application in the field of remote sensing (RS).

This field has drawn a great interest with the development of multiple types of sensors, leading to a huge heterogeneous amount of data. To overcome the curse of dimensionality while providing a powerful tool for environment monitoring, planning, and decision-making, several techniques have been employed. Convolutional neural network or the deep learning79 represents a tendency in this field thanks to its ability to automatically discover relevant contextual features in image categorization problems.10 In the same context, sparse representation (SR),1113 total variation,14,15 and machine learning16,17 techniques, to name a few, are also investigated for different purposes, such as enhancing spatial resolution, generating explanations, and extracting knowledge from the images. MGD, the scope of our paper, have been extensively used in several domains, namely pattern recognition18,19 and computer vision.20,21 Their use was also extended to the RS field, where we need to locate edges of roads, building, rivers, and forest to detect a potential change. Adding to that, MGD provide an analytical treatment of a scene by decomposing it in high frequencies and low frequencies. This essentially helps in fusing information from different sensors and in revealing hidden characteristics indiscriminately using only one sensor image. Therefore, we aim, in this work, at drawing attention to the MGD importance and contributions in the RS field. To the authors’ knowledge, several reviews were elaborated to describe MGD, but few of them focused on their particular use in RS.

This paper is organized as follows: we present, in the first section, an overview of MGD. Then we describe their characteristics and how they have been used in the field of RS. We conclude this paper with a discussion and a conclusion.

The notion of scales was introduced before wavelet emergence to better localize objects in the image observation. Methods like windowing or scale introduction in Fourier transform,22,23 Laplacian pyramid (LP), and derivatives of Gaussian (DoG) were the cornerstones of ideas a few years later. The definition of the continuous wavelet was a kind of generalization of previous works in scale incorporation since the newborn transform has offered a localization in time and frequency.

For numerical computation needs, the wavelet was discretized using multiresolution analysis (MRA). The MRA consists of projecting a signal successively onto subspaces Vj (as jZ increases, the approximation becomes coarser) and yielding both approximation and detail information (detail information is calculated from two successive approximations). Otherwise, MRA stipulates that any signal, in our case an image, could be constructed iteratively by exhibiting different characteristics in every scale. Figure 1 shows how smooth regions are emphasized in finest scale while contours are more likely to be salient in the coarsest ones. Taking advantage of the MRA, wavelet transform can address point-like singularities and offer less redundant information in scales, compared to LP and DoG. Besides, wavelet is characterized by being separable, which means that its 2-D function atoms could be written as the product of two other 1-D functions. Precisely, a wavelet’s atom ψ(x)=2jψk(2jxn) is obtained by three tensor products 2-D wavelets as follows:

  • ψV(x)=ψ(x1)φ(x2)
  • ψH(x)=φ(x1)ψ(x2)
  • ψD(x)=ψ(x1)ψ(x2)

where x=(x1,x2)R2 refers to spatial coordinates of a given point x, 2j refers to dyadic scalings, jZ refers to scale, nN is translation factor, φ is 1-D orthogonal scaling (acts like low-pass filter), ψ is wavelet function (acts like band-pass filter), and k stands for the only three directions supported by wavelets (horizontal, vertical, and diagonal).

Graphic Jump Location
Fig. 1
F1 :

(a) Smooth region revealed in finest scale (b) while contours/edges revealed in coarsest scale.

The MGD were introduced to further improve the aspect of directionality limitation, whether by using the wavelet itself combined with other treatments or by defining rules for wavelet functions. These methods can be divided into two different families:

  • adaptive family, which follows a Lagrangian approach, such as bandlet and
  • nonadaptive family, which follows the Eurelian approach, such as CT.

The first family proposes to build an adapted structure of the data. For example, in order to locate singularities, the bandlet decomposes the image using wavelets combined with dyadic decompositions. Each square of the resulted segmentation contains one direction, i.e., a unique piecewise singularity. However, the second family tries to ameliorate wavelet by defining additive mathematical rules on the wavelet functions. This brings more accuracy to detect image content. A debate was settled on the strength of those families and their abilities in image understanding: smooth regions, contours, and texture. According to Ref. 24, the best representation is left an open question since it is data dependent.

In this section, we represent some of the MGD and describe their main characteristics and differences. Particularly, in this paper, we are interested in some MGD that are applied in the RS field. Figure 2 illustrates a map of MGD elaborated based on Ref. 24. In fact, we have written the names of methods that will be discussed in bold, so that the reader can have a clear vision of what is treated in the coming sub/sections.

Graphic Jump Location
Fig. 2
F2 :

MGD map locating methods discussed in this paper.

Curvelet Transform

The first generation of CT was proposed by Donoho and Duncan.25 By that definition, this transform extended the ridgelets (built on the ground of 1-D wavelets to detect segment singularity)26 to switch from a global scale analysis to fine scale analysis. The idea was to decompose the image into a set of subimages and perform ridgelet transform on each one of them.

This CT transform has not shown great results because of the failure of ridgelets in diagnosing the C2 edges. The second generation of CT is then proposed and was less redundant and faster compared to its first generation.27 Given a function fL2(R2), the coefficients C(j,l,k) of the continuous CT transform are defined by Eq. (1): Display Formula

C(j,l,k)=f,ϕj,l,k=R2f(x)ϕj,l,k(x)dx,(1)
where ϕj,l,k is a CT atom, at scale j, θl the angle of the orientation l, and position k=(k1,k2)R2. In the frequency domain, CT coefficients are represented at discrete intervals formed by sampling the continuous curvelet transform at dyadic intervals Ij=2j, θl=2π2j2l, in an equispaced rotated anisotropic grid, as shown in Fig. 3.

Graphic Jump Location
Fig. 3
F3 :

CT frequency space tiling: to get a discrete representation, the rotated tiling is transformed into squared tiling. (a) Radial tiling, (b) angular tiling, and (c) the resulting frequency tiling of CT.

CT’s atoms are defined by a combination of a translation xk(j,l) and a rotation Rθl of an angle θl of an atom ϕj,l,k(x): Display Formula

ϕj,l,k(x)=ϕj[Rθl(xxk(j,l))].(2)
The rotation with an angle θl is defined as in Eq. (3): Display Formula
Rθl=(cosθlsinθlsinθlcosθl)et  Rθl1=(Rθl).(3)

The atom of CT ϕj is defined in the frequency domain by the mean of its Fourier transform ϕj(ω)^=Uj(r,θ), which is written in polar coordinates as follows: Display Formula

Uj(r,θ)=23j4Wj(2jr)Vj(2j2θ),(4)
where Wj is the radial window and Vj is the angular window at a given scale j.

To discretize the CT, the authors of this transform propose to switch from rotated tiling to angular tiling and from concentric circles to concentric squares, as illustrated in Fig. 3.

Contourlet Transform

Contourlet is considered as the extension of Candes and Donoho’s work to improve the isotropic criterion of wavelet representation.28 This transform aims at obtaining the same frequency space tiling as the CT, without the need to move from continuous to discrete domain. To do so, the authors propose a nonseparable decomposition scheme by applying a DFB, combined with LP [as shown in Fig. 4(a)] and they also propose to decrease the high redundancy information in its sub-bands.

Graphic Jump Location
Fig. 4
F4 :

(a) Contourlet filter bank, the first treatment is to decompose image in multiscale sub-bands coefficients due to LP. Second, a directional decomposition filter will be applied on these coefficients. (b) Contourlet frequency space tiling.28

The DFB is designed to capture anisotropic features in the high frequencies of the image by allowing different numbers of direction at each scale while the low frequencies are processed by LP. In fact, as presented in Ref. 28, LP is defined based on orthogonal filters and downsampling by 2 in each iteration, similar to wavelet. The low-pass synthesis filter G, used for LP processing, defines a unique scaling function φ(t)L2(R2), tZ that satisfies Eq. (5): Display Formula

φ(t)=2nZ2g[n]φ(2tn)andLet:φj,n=2jφt2jn2j,(5)
where jZ is the bandpass number, nN2 is the coordinates of a given point, and φj,n is a member of the family {φj,n}nZ2, which is an orthonormal basis for an approximation subspace Vj at the scale 2j. {Vj}jZ is a sequence of multiresolution nested subspaces.

Otherwise, the DFB offers localization and directionality properties due to the family, defined in Eq. (6): Display Formula

{dk(l)[nSk(l)m]}0k<2l,mZ2,(6)
which are obtained by translating the impulse responses of the equivalent synthesis filters dk(l) over the sampling lattices by Skl as follows: Display Formula
Skl={diag(2l1,2)for  0k<2l1diag(2,2l1)for  2l1<k2l,
where Skl corresponds to the mostly horizontal and mostly vertical sets of directions, respectively, and l refers to level tree-structured DFB, which yields 2l real wedge-shaped frequency bands (l=3; 23=8 real wedge frequency bands). The frequency tiling of this transform is shown in Fig. 4(b). The contourlet has preserved and improved the characteristics of CT. Besides, better than wavelet, the contourlet can represent edges and singularities along curves. However, the latter transform suffers from pseudo-Gibbs effect due to downsampling and upsampling. To overcome its weakness, a modified version, built on the ground of contourlet conception, was designed. The nonsubsampled contourlet transform (NSCT)29 discards the sampling step and is characterized by being shift-invariant and by preserving multiresolution criteria.

Shearlet Transform

Shearlets30 provide a rich mathematical structure and are optimally sparse in representing the edges within an image thanks to the fact that they form a tight Parseval frame at various scales and directions. They are distinguished from CTs by being directly constructed in discrete domain, which gives them the ability to provide a more efficient multiresolution representation of the geometry.31 Besides, shearlet outperforms contourlet by substituting the directional filter with shear filter, which helps in breaking the limitation of directionalities. In the continuous domain, in dimension 2, a shearlet represents an affine system with composite dilations in Eq. (7): Display Formula

ϕj,s,k(x)=a34ϕ[Us1Va1(xk)],(7)
where ϕj,s,k is a shearlet atom, j refers to the scale, s refers to the shear direction, and k to the translation parameter, where Va=(a0aa) is the anisotropic dilation matrix (controls the support of the shearlet function) and Us=(1101) is the shear matrix (controls the orientation). As is presented in Ref. 32, the discrete collection of shearlet transform could be defined by two window functions localized on a pair of trapezoids, as illustrated in Fig. 5(a).

Graphic Jump Location
Fig. 5
F5 :

The structure of the frequency tiling by the shearlet: (a) the tiling of the frequency plane induced by the shearlets. (b) The size of the frequency support of a shearlets. (c) Shearlets conception using filters: LP combined with directional filter (shear filter).33

In fact, for a given point in the frequency domain ξ=(ξ1,ξ2)R2^, j0, and l=2j,,2j1, let Display Formula

Wj,l(0)(ξ)={ψ2^(2jξ2ξ1l)χD0(ξ)+ψ2^(2jξ1ξ2l+1)χD1(ξ)if  l=2jψ2^(2jξ2ξ1l)χD0(ξ)+ψ2^(2jξ1ξ2l1)χD1(ξ)if  l=2j1ψ2^(2jξ2ξ1l)otherwise,
and Display Formula
Wj,l(1)(ξ)={ψ2^(2jξ2ξ1l+1)χD0(ξ)+ψ2^(2jξ1ξ2l)χD1(ξ)if  l=2jψ2^(2jξ2ξ1l1)χD0(ξ)+ψ2^(2jξ1ξ2l)χD1(ξ)if  l=2j1ψ2^(2jξ2ξ1l)otherwise,
where Wj,l(0)(ξ) and Wj,l(1)(ξ) are the mathematical definitions of vertical and horizontal cones, respectively, ψ2 is a function satisfying properties in Ref. 32, and χD0 and χD1 are, respectively, truncated vertical and horizontal cones. The authors of this transform have also proposed a fast algorithmic implementation and a filter bank architecture which is similar to the one of contourlet transform as shown in Fig. 5(c)].

Wedgelet Transform

Wedgelet is proposed by Donoho et al.34 This transform decomposes the image iteratively into piecewise constant functions (class of horizon functions). The first step consists of a dyadic recursive decomposition. In each quadtree leaf, wedgelet will search for an “edgel” in order to forward the leaf dyadic decomposition. Figure 6(a) represents a wedgelet decomposition of “cameraman” image and its dyadic decomposition.

Graphic Jump Location
Fig. 6
F6 :

(a) The wedgelet segmentation process. (b) (i) Decomposition based on piecewise constant function (Wedgelet) and (ii) Decomposition based on polynomial function (Surflet).35

The intuitive evolution of this transform is to switch the constant horizon functions by polynomial functions in order to catch more irregular singularities. In fact, this was proposed with the surflet.36Figure 6(b) illustrates how the dyadic decomposition is influenced when polynomial functions are used instead of piecewise constant functions.

Bandlet Transform

The first generation of the Bandlet transform was proposed by Le Pennec and Mallat.37 The main idea was to explore the wavelet transform and overcome its failure in detecting anisotropic regularities, such as C2 curvatures. According to Ref. 37, Bandlet first performs the grassfire algorithm38 to detect edge regions and then they perform a deformation to let it fit wavelet orientations (horizontal, vertical, or diagonal). This version was redundant and failed to ensure good results, mainly in compression applications.35

Thus, the second generation of bandlet is proposed. Respecting the same adaptive approach, this version is considered as an anisotropic wavelet wrapped along the geometry flows.39 This version differs from its previous version by being based on quadtree segmentation algorithm. The idea is first to estimate irregularities in the image as a vector field, generated using the edge polynomial function estimation. A wavelet transform is then applied, after that a dyadic segmentation is performed in wavelet domain. This way, wavelet coefficients are calculated along the optimally found direction in each square. The bandlet could also be seen as a polynomial approximation in an orthogonal Alpert basis. Alpert transform is simply a wrapped wavelet transform adapted to a given irregular form, as illustrated in Fig. 7. The coefficients of the bandlet are calculated by following Eq. (8): Display Formula

bj,l,nk=pal,n[p]ψj,pk(x),(8)
where al,n[p] is an Alpert coefficient of a given point pR2, ψj,pk is the wavelet mother function at p, j refers to the scale, and l refers to the factor defining the elongation of the bandlet function. Thus, bandlet coefficients are generated by inner products f,bj,l,nk of the image f with bandlet function bj,l,nk.

Graphic Jump Location
Fig. 7
F7 :

(a) The original image I, (b) wavelet coefficients of I, and (c) zoom on wavelet coefficients in a square S including edge and wrapping operation w of the geometric flow to align it horizontally or vertically.40

There is no doubt that wavelet was a genuine processing tool that has led to major advances in natural image representation and understanding. Since then, many transformations have been defined to overcome the shortage of wavelet capacities by adding directions, shift invariance criteria, and so on. In Table 1, we summarized the characteristics of different transforms mentioned in this review. We have separated the MGD into two families. Their differences reside on the fact that the adaptive family constructs a special data decomposition that can fit each dataset rather than using a predefined system like the nonadaptive families. However, this adaptive form is very intensive in terms of numerical computation.1

Table Grahic Jump Location
Table 1MGD characteristics.

But, in general, all these transforms exhibit interesting parameters, which helps to bring attention to different details within an image.

  • Scale: MGD locates, in each scale, a specific feature (edge in the coarsest scale, smooth content in the finest). Keeping some scales while withdrawing others could help remove undesirable information. Indeed, discarding disturbance that could be part of the finest scales helps in providing a better representation of the image
  • Direction: Likewise, directions reveal features in the image. We can detect, thanks to this parameter, the main axes by locating the accumulation of strong coefficients magnitude. Consequently, we can enhance some image structures. Figure 8 illustrates a rendering of the content of scales and supported directions of some MGD.
  • Magnitude: In transformation domain, the greatest magnitude values could be interpreted as an influent part in the processed image. This is why the magnitude is useful for feature extraction applications. In fact, since the MGD allow the zooming in and out, the content of the image, especially texture, will not be affected by the size of the neighborhood. Several first-order and second-order statistics could be calculated: mean, standard deviation (SD), energy, entropy, contrast, sum of mean, variance, and so on. Consequently, several descriptors could be extracted and could be dimensionally reduced without affecting the discriminative power.42 Another interesting use of magnitude is detecting spatiotemporal directions as in Ref. 43.

Graphic Jump Location
Fig. 8
F8 :

Reconstructing images keeping only some scales. (a and b) Second and fourth scale of the CT. (c and d) Second and fourth scale of the contourlet. (e and f) Second and fourth scale of the shearlet. (g and h) Second and fourth scale of the bandlet first generation.

  • Frames: CT, shearlet, and contourlet are a set of functions forming a tight frame. This means that for a given sequence (ϕi)iI in Hilbert space H, where I is a countable indexing set, i is an element of I, and for all x in H, there are two constants 0<AB: Display Formula
    Ax2iI|x,ϕi|2Bx2.(9)
    When A=B, this is called a tight-frame, and if A=B=1 this is called a Parseval frame. The advantage of having such a function system is to be able to represent a signal as a linear combination of the vectors within the frame in several ways and to reconstruct it by the use of inner products [Eq. (10)]. Let S=T*T, where T is the analysis operator and T* is the synthesis operator of a given frame. T and T* will be explained in “sparsity” paragraph: Display Formula
    S:HHxiIx,ϕiϕi.(10)

The authors of bandlet have also proposed its tight frame version, which is called grouplet.44 For further details on this frame version, we refer the interested reader to Ref. 45.

  • Redundancy/overcompleteness: Although the orthonormal/orthogonal bases has interesting properties, they have several weaknesses. They are translation and rotation sensitive, especially when it comes to dealing with multidimensional data. That is why the MGD, such as curvelets, shearlets, and contourlets, are redundant. This means that the vectors of their basis are linearly dependent and exceeded the number of their space dimension. Thus, a given vector in this kind of space would have an infinite set of representations. In fact, the authors of these transforms want to move “cautiously” to overcompleteness without losing the interesting features of the orthonormal/orthogonal basis.
  • Sparsity: The aforementioned MGD have two main operators: the analysis T and the synthesis T*. Given a signal x, MGD analyze x by decomposing it in a set of coefficients: Display Formula
    T:H2(I)x(x,ϕi)iI.(11)
    The dual operator, the synthesis, helps in recovering a given x: Display Formula
    T*:2(I)H[(ci)iI]iIciϕi.(12)
    Curvelets, contourlets and shearlets use fast algorithms to resolve Eq. (11). A suitable choice of a representation system enables a SR which means that a given signal x could be expressed in a low-dimensional space. This representation is widely used in several approaches such as compressed sensing.46Figure 9 shows the huge difference between representing an image using histograms of a gray level and representing it sparsely using curvelets (as an example of the effect of the use of MGD).

RS images represent an efficient tool to analyze climate change and dynamics of land cover. In the last decades, a great number of new satellites has been launched, allowing a tremendous data availability with improved spatial and spectral resolutions. This has helped in enhancing our understanding and control of our surroundings. Throughout the various RS applications (image fusion, enhancement, super resolution, and so on), we are capable of monitoring the earth’s surface, predicting changes, and preventing disasters.

Graphic Jump Location
Fig. 9
F9 :

Histogram of an image in (left) the original (pixel) domain and the CT domain (right).47

Whereas, the management of RS images represents a challenging task. As the amount of data is incredibly growing, it is getting more complex to extract knowledge from this type of images. RS data are not only different but also have a rapid velocity and generally need several complicated corrections. That is why researchers seek new approaches to replace the traditional data processing algorithms.

MGD were used in an attempt to solve some of the aforementioned problems. In fact, being able to represent the datasets in a sparse way while describing accurately, the objects in the image make the MGD an inspirational tool to be exploited in RS field. In this section, we are interested especially in the use of wavelet, CT, contourlet, shearlet, and bandlet.

Classification

The classification aims at recognizing the class label of a given study area with the aid of ground truth data. Conventional classification algorithms exploit essentially the spectral information of the images.48 Nevertheless, it has been proven that incorporating spatial information in the classification process, i.e., taking the contextual information into account, helps significantly in boosting the obtained accuracies. The spectral–spatial combination is possible using MGD, which not only provide a frequency representation but also help in establishing correlations between neighbors in the spatial domain.

From the studied cases, we can separate the MGD use in classification in two categories:

  • Decomposing the whole image and using its frequency details for classification
  • Combining the MGD frequency details with other features obtained using other methods.

We have studied the Ref. 49 where wavelet has been used for a feature extraction purpose. Combined with fuzzy hybridization, the wavelet coefficients are exploited for classification. The three multispectral (MS) RS images, two from Indian remote sensing and one from SPOT, are decomposed band-by-band. This technique has essentially helped in accounting for contextual information. The overall classification accuracy (Mean=79.15%) using ground truth of the three images has shown that the biorthogonal wavelet exceeded the other wavelet functions.

In Ref. 50, the authors proposed a CT-based approach to extract features from SAR images that can effectively identify the dynamic ice from the consolidated one. The representation of the dynamic ice contains curves in different locations with different widths and lengths and is considered larger than those of consolidated ice representation. The proposed approach consists of extracting patches from the transform domain sub-bands according to a specific size. This later increases according to the CT’s scale in order to capture more information about the curves. Once the size of patches is computed, a feature vector is calculated by Eq. (13): Display Formula

F(PL)=1θij=1θi(Fi,j),(13)
where PL is the patch with a size L, Fi,j denotes the mean of CT coefficients at scale i, and θiN is the total number of orientations at scale i.

Using SVM along with different features, the obtained experimentation results proved that the CT-based feature extraction is effective for classification since the dynamic ice area is more accurately classified. In fact, thanks to the probability density function (PDF), we can see clearly that the component representing the different types of ice could be distinguished and thus separated (Fig. 10).

Graphic Jump Location
Fig. 10
F10 :

Estimated PDF for different features: (a) the PDF of SD, (b) gray-level co-occurrence matrix, and (c) CT-based method. The extracted features overlap a lot, which means they will be less effective in classification except for the PDF of CT-based method, which shows that the different types of ices could be distinguished in classification.50

Nevertheless, this quality of recognition loses its precision in CT’s finer scale.

The contourlet is used in Ref. 48, where a performance comparison between wavelet and contourlet is discussed. Based on the fact that wavelet suffers from its shortage of directionality and the fact that contourlet provides directions only in high frequency coefficients, a wavelet-based contourlet transform (WBCT) is proposed and is applied on linear imaging self-scanner (LISS) II, III, and IV. Figure 11 illustrates the proposed transform.

Graphic Jump Location
Fig. 11
F11 :

Filter bank of wavelet-based contourlet transform proposed in Ref. 48.

After decomposing the image using the WBCT, wavelet, and contourlet, PCA is applied to the obtained features in order to reduce their dimensionality while removing redundancy and preserving the most discriminant ones. Then a mean vector and a covariance matrix are calculated. Finally, the Gaussian kernel fuzzy C-means classifier is applied, and the obtained overall classification accuracy proves that the proposed decomposition (overall accuracy of LISSIV=89.57%) is better than the wavelet-based (overall accuracy of LISSIV=87.62%) and contourlet-based feature extraction methods (overall accuracy of LISSIV=88.38%).

As we mentioned earlier, the second category combines MGD sub-bands as features describing edges or textures, with features obtained using other methods. For example, the authors of Ref. 51 suggest to use CT jointly with morphological component analysis (MCA) and to improve RS classification using hyperspectral AVIRIS and AirSAR images. The idea is about separating a given image into two components (14): Display Formula

y=ys+yt+n,(14)
where ys is the smooth component, yt is the texture component, and n is the noise. Here, the CT is used to construct a dictionary As representing the smooth component and a Gabor wavelet is used to construct a dictionary At representing the texture component. The authors propose to estimate the component ys and yt by solving Eq. (15) using SunSAL algorithm.52Display Formula
ys^,yt^=argminys,yt12yysyt22+λ1Tsys+λ2Ttyt,(15)
where Ts=(TsTAs)1AsT and Tt=(TtTAt)1AtT are the pseudoinverse of As and At, respectively.

This approach was extended in Ref. 53, where the authors combined several methods such as Gabor wavelet and horizontal filters to construct an MCA kernel for feature extraction. The use of such a composite kernel has given better results (overall accuracy of AVIRIS Indian pines=93.54%) in characterizing image content compared to minimum noise fraction (MNF) components and helped in enhancing the characterization of the image’s content. This is explained by the fact that the proposed approach combines several methods, such as wavelet and CT.

In a similar fashion, the authors of Ref. 54 proposed to combine multiple features pertaining to spectral, texture, and shape, and proposed a multiple feature combining (MFC) framework, as shown in Fig. 12. The spectral feature of a given pixel is elaborated by arranging its digital number in all of the l bands. The texture feature is obtained by applying 2-D Gabor wavelet filter and the shape feature is constructed due by the pixel shape index method.55 To calculate the feature vectors using the MFC framework, a single feature-based dimensionality reduction technique is conducted in order to generate the alignment matrix. Then a Lagrangian function is calculated in order to determine the low-dimensional feature space Y (16): Display Formula

L(ω,λ)=Σi=1mωirtr(YM(i)YT)λ(Σi=1m1),(16)
where MRN×N is the alignment matrix of input samples, m is the number of features (in this case 3) and ωi is a relaxation factor with r>1. After that, Y is linearized using an explicit linear projection matrix. This method outperformed several others like principal component analysis (PCA), MNF, locally linear embedding, local tangent space alignment, Laplacian eigenmaps, and achieved the optimal classification performance on HYDICE HI and ROSIS hyperspectral datasets. The criteria of assessment used in this classification were the computation of overall accuracy and the kappa index.

Graphic Jump Location
Fig. 12
F12 :

Multiple features of the airborne data over the mall in Washington DC dataset. (a–c) Spectral feature images in band 36, 52, and 65. (d–f) Gabor texture feature images, with direction=1 and scale=1, 3, and 5, respectively. (g–i) Shape feature images in direction 1, direction 8, and direction 16, respectively.54

Change Detection

The detection of a change in land cover or in the Earth’s surface is considered as one of the most important applications in RS. In fact, it helps in disaster management, vegetation development, deforestation detection, and urban growth tracking. This application needs a set of multitemporal satellite images. Throughout the studied cases in this section, MGD were extensively used to ensure an accurate detection of affected regions in multitemporal images by enhancing the image’s details, especially the difference images (DIs), where noise information can easily be interfered.

In Ref. 56, the contourlet is used to denoise the images of interest in order to enhance the change detection process. First, it is applied to each single temporal SAR image to preserve its features and edges. The authors also proposed to reduce speckle noise by performing hard thresholding on its high frequency sub-bands using Eq. (17): Display Formula

T(m)=δfi,j2(m)2lgL(m),(17)
where δfi,j2(m) is the variance magnitude, L(m) is the number of the decomposition coefficients at the scale m, g is the magnitude value of the pixel, and l refers to the direction. Then the best decomposition scale is calculated by finding the minimum of the local variance of river courses magnitude. Figure 13 illustrates how the river courses are represented in SAR images and their rendering after applying contourlet.

Graphic Jump Location
Fig. 13
F13 :

The river course in an SAR images: (a) original image and (b) after denoising.

After that, markers are extracted from the contourlet domain in order to find the potential course rivers and eliminate the false alarms (content detected as river courses while it is not). After being processed, the SAR images become smoother than the ones before applying the noise and the speckle reduction. In contrast, the edges are preserved and become more visible. That is why the authors affirmed that this contourlet-based method can achieve higher accuracy of river courses detection.

In Ref. 57, the CT is used jointly with MCA and is applied to a high-resolution airbone SAR image to detect changes in urban areas acquisition. The MCA decomposes the image into K different components, as specified in Eq. (18): Display Formula

minα1Kk=1Kαk11s.t.  yk=1Kψkαk2<σ,(18)
where ψk is a transform function, αk is the matrix coefficients of the image y in the transformed domain, and σ is the error (fidelity to data). In this work, the authors propose an approach using CT and wavelet-based MCA. The idea is to calculate two DIs. The first DI is obtained by computing the difference between the two SAR images. The second one is obtained by downsampling the first one. This technique is usually considered as the simplest noise suppression mechanism to use in a change detector. Then the two DIs are elaborated thanks to wavelet and CT-based MCA components as shown in Eq. (19): Display Formula
DIMCA=ψcurveletαcurvelet+ψswtdb4αswtdb4,DIdMCA=ψcurveletαcurveletd+ψswtdb8αswtdb8d,(19)
where the first component of DIMCA is ψcurveletαcurveletcalculated using CT, ψcurvelet is the CT function, and αcurvelet is the set of coefficients in CT domain. The second component of DIMCA is calculated using wavelet, where ψswtdb4 is the Daubechies wavelet function and αswtdb4 are the coefficients in wavelet domain. The DIdMCA corresponds to the downsampled DI by factor L and it is also built using CT and wavelet-based MCA component, but this time using a different vanishing moment for the wavelet.

To suppress the undesirable information, the authors applied a soft-threshold on the CT and wavelet coefficients. Then they calculated the change map, as illustrated by Eq. (20): Display Formula

G=FandFCurv,(20)
where F corresponds to a binary map obtained by Eq. (21): Display Formula
F=F1andF2=(|DIMCA|T1)and(|DIdMCA|L)T1),(21)
and FCurv corresponds to the CT component, as given in Eq. (22): Display Formula
Fcurv=(|ψcurvαcurv|T1)and(|ψcurvαcurvd|L)T1.(22)
The CT component, Fcurv, is used to characterize the geometry of contours and edges. It also contains noise due to speckle. This noise is reduced due to the threshold T1, as is mentioned in Eq. (22). By calculating the change map of DI illustrated in Eq. (20), the authors aimed at obtaining a change map free of speckle noise, while conserving the contours. The experimentations show that this method can withdraw the undesirable noise without removing the cars and buildings edges. But only the geometry that has less than 8 deg of incidence angle could be maintained.

In Ref. 58, the authors detected anomaly in hyperspectral images by, first, using the shearlet to decompose the images into several directional sub-bands at multiple scales. Then, in each sub-band, the background signal is reduced and the fourth-order partial differential equation is applied to brighten up the anomaly. Experimental results with HYDICE HI data show that the presented algorithm can suppress the background, detect the anomaly signal effectively, and outperform the original RX algorithm.59

Another approach is proposed using a detail-enhancing approach using NSCT Ref. 60, where the authors suggested to detect change from multitemporal optical images using a detail-enhancing approach. To do so, they studied and analyzed two change detection methods. The first61 is based on combining PCA and K-means, which was efficient in terms of computational time, but since the PCA does not consider multiscale processing, when fusing data, this has yielded false detections. The second method was proposed by Ref. 62 and introduces a multiscale framework for change detection in multitemporal SAR images. This method decomposes a DI into S scales by the undecimated discrete wavelet transform (DWT). Based on this decomposition, the authors extracted a feature vector by first computing the intrascale features (sampling local neighborhood using two methods) in each sub-band in a given scale, and then by computing the interscale features, which regroup all the intrascale vectors together.

The same authors improved this work in Ref. 63 by proposing the use of CT-DWT (complex wavelet) to capture further directions. The authors also used the Bayesian inference to calculate a threshold to decide whether the details to fuse correspond to changed or unchanged pixels.

In Ref. 60, both works Refs. 61 and 62, were discussed. The authors estimated that the use of K-means algorithm, in Ref. 61, can stick to local optima, which can produce false detections, especially when choosing an inappropriate initial centroid. Moreover, they also affirmed that the use of wavelet, in Ref. 62, failed to characterize the geometric details of the DI efficiently. To overcome this weakness, instead of wavelet, the authors used the NSCT to decompose the DI into S scales. Each scale yielded one low pass approximation sub-band DSA and L high-pass directional sub-bands D1,1HDS,LH. The details are extracted from directional sub-bands. The high-pass directional sub-bands serve to extract the details and are thresholded to decrease the amount of noise. The intrascale fusion is performed using max rule, as shown in Eq. (23): Display Formula

DsDetail(i,j)=max{abs[Ds,lH,T(i,j)]}l=1L,(23)
where T is a threshold, s is the scale, l is a specific orientation of L obtained directions from the NSCT filter bank, and H corresponds to high-pass. Then the interscale fusion is applied using the max rule fusion, as well in Eq. (24): Display Formula
DDetail(i,j)=max{abs[DsDetail(i,j)]}s=1S.(24)

Then an enhanced DI is calculated [Eq. (25)]: Display Formula

DEnhanced=DBase+βDDetail,(25)
where DBase is obtained from the finest scale approximation from NSCT decomposition and β is a weight to balance the emphasis between the base image and the detail image. Then the authors extract patches from DEnhanced, transform the patch into vector via lexicographic ordering and use PCA to produce principal components. Finally, a PCA-guided K-means64 is performed to calculate the change map. Compared to approaches based on EM, on Bayesian, on PCA, or on multiscale, the proposed approach gives better results since it conserves the geometric details due to the NSCT. The results of experimentations are illustrated in Fig. 14.

Graphic Jump Location
Fig. 14
F14 :

Qualitative change detection results by using different change detection methods. (a) Ground truth change detection map created by manually analyzing the corresponding input images, (b) EM-based approach, (c) Bayesian-based approach, (d) PCA-based approach, (e) multiscale-based approach, and (f) proposed method.60

In Ref. 65, the proposed method applies wavelet on MS imagery in an anisotropic diffusion aggregation. The proposed approach is composed of three steps. We only cite, in depth, the steps where wavelet is used:

  • A band selection is first conducted to choose where the texture is best captured. Then, 2-D DWT for multiscale-multidirectional texture extraction is applied.
  • Textural and spectral segmentation are performed by anisotropic diffusion in order to reduce noise without blurring inter-region edges as well as creating the desired multiscale low-level primitives.
  • Change detection is then performed in two steps:
    • classification of the image into forest/nonforest thematic maps through thresholding
    • comparison of the classification map, obtained from the previous step, using logical modeling in order to identify changes

The experiments are conducted on Landsat TM and Landsat ETM+ datasets dated 1986 and 2001, respectively. The criteria of assessment used in this work are the transformed divergence measure which statistically determines the adequate wavelet levels to keep, overall accuracy and kappa.

In Ref. 66, a CT-based change detection algorithm is proposed between two co-registered SAR images for natural disaster mapping. After applying the CT on the two images, the coefficients are weighted to suppress noise-like structures using the mathematical relation [Eq. (26)]: Display Formula

Lx,y=Σi=1nCi.ki,(26)
where Lx,y is the amplitude found at position (x,y) in the SAR image, Ci are the sum of CT coefficients, n is the number of CTs coefficients, and ki is a complex coefficient varying according to the image content. Finally, the change in radar amplitude Dx,y is calculated by Eq. (27): Display Formula
Dx,y=L2x,yL1x,y.(27)

Fusion

Image fusion represents one of the most important RS applications. It aims at combining two or more images to create another one with enhanced features. This increases the possibility of taking advantage of multisensors images and opens the door for several uses. Generally speaking, fusion can be categorized into four families based on when the fusion rule is applied: fusion at signal level, at pixel level, feature level, and decision level. The MGD belong to the feature level family. In fact, in order to conduct a fusion process with these latter decompositions, we should first convert the intrinsic image properties and get the frequency coefficients. Fusion techniques, based on MGD, suggest several injection models but tend basically to extract the spatial detail information from high spatial resolution images to inject them in the low spatial resolution ones. Compared to the largely used methods like intensity hue saturation (IHS) and PCA, MGD enhance the spatial characteristics of fused images. Thus, they look very clear, have sharp edges, and are free of spectral distortion.67 MGD-based fusion are applied abundantly in pansharpening techniques, multisensor fusion, and spatial enhancement.

In Ref. 68, the authors propose a technique to fuse MS satellite images (high spectral and low spatial resolution) with a panchromatic (PAN) satellite image with low spectral and high spatial resolutions. They introduced an improved method of image fusion based on the ARSIS concept using the CT transform. The ARSIS is a French abbreviation for enhancing spatial resolution by injecting structures. CT-based image fusion has been used to merge a Landsat enhanced thematic mapper plus, PAN and MS images. Based on experimental results using indicators of bias, the CT-based method provides better visual and quantitative results for RS fusion. The indicators are

  • SD which indicates the dispersion degree between the gray values and the gray mean values,
  • correlation coefficients (CCs) which indicate the degree of correlation between two images,
  • the relative average spectral error characterizes the average performance of a method in the considered spectral bands,
  • relative global dimensional synthesis error (ERGAS) which measures the spectral distortion between the reference image and the fused one.

Shearlet and nonsubsampled shearlet transform (NSST) are also used in this regard. In Ref. 69, the authors propose in the first step to register two RS images. Then the shearlet transform is applied. The number of directions is set to 6 and scales to 5. The low frequency coefficients are selected based on the average rule while the high frequency coefficients are chosen using the maximum absolute value rule as mentioned in Eq. (28): Display Formula

F(i,j)={A(i,j)DA(i,j)DB(i,j)B(i,j)DA(i,j)DB(i,j),(28)
where DX(i,j)=iM,jN|YX(i,j)|, X=A,B is the absolute value of high frequency coefficients in the neighborhood of a pixel value YX at a location (i,j). M and N are equal to 3 and correspond to the size of the neighborhood window. X denotes the two source images. To enhance the fusion rule for high frequencies, a decision map is produced by affecting 1 if DA(i,j)DB(i,j) and 0 otherwise. If some coefficients come from an image A while all its neighbors are from B, then the pixel will be extracted from B. The fusion results are illustrated in Figs. 15(a) and 15(b), which are two RS images with different band characteristics. Those two images have different physical properties of the sensor. Figures 15(c)15(i) represent the fused images with other methods. Using shearlet, the authors assert that the obtained fusion image has clearly preserved characteristics of the surface. They found also that fusion result using contourlet has higher sharpness and entropy values than shearlet.

Graphic Jump Location
Fig. 15
F15 :

Comparing shearlet fusion algorithm with other methods.69 (a) Eight-band RS image, (b) three-band RS image, (c) shearlet, (d) contourlet, (e) Haar, (f) Daubechies, (g) LP, (h) average, and (i) PCA.

In Ref. 70, the authors have used NSST.32 They have involved the SR71 paradigm in pansharpening application.72 The idea is about fusion the intensity component (I) of the MS image (the hue and saturation components are not used) with the PAN image. First, the authors decompose the images using NSST into low frequencies {LI,LPAN} and high frequencies {HI(s,d),HPAN(s,d)}, where s corresponds to scale and d to direction. For the low frequency coefficients, they propose to construct vectors using patches from both LI and LPAN. Those vectors are represented sparsely using a learned dictionary with K-SVD73 as mentioned in Eq. (29). Display Formula

αi=argminαα0subject to  viDα22ε,ε0,(29)
where vi is the i’th vector patch. The D corresponds to the dictionary obtained by different ways: DCT, NSST and trained by K-SVD algorithm and α is the sparse vector. The “absolute max” fusion rule is adopted to obtain a fused sparse coefficients (30): Display Formula
αFi={αIiif  |αIi||αPANi|αPANiotherwise.(30)

Then the new fused vector patch vFi is reconstructed using the dictionary D and the fused sparse coefficients are obtained as follows (31): Display Formula

vFi=DαFi.(31)

Otherwise, the high frequency coefficients are fused according to large local energy rule. This means that the energy of each scale s and direction d in the shearlet transform is calculated, as mentioned in Eq. (32) and then fused: Display Formula

HF(s,d)(i,j)={HI(s,d)(i,j)if  LEI(s,d)(i,j)HPAN(s,d)(i,j)HPAN(s,d)(i,j)otherwise.(32)

The experimental results show that the proposed method conserves better spatial and spectral information and is able to enhance the fused image better than Refs. 74 and 75. The render result is illustrated in Fig. 16. We can notice how well preserved the spatial details and color information are. Fused images have better visual accuracy compared to existing methods: AIHS,76 AWLP,75 SVT,74 and ATWT.77

Graphic Jump Location
Fig. 16
F16 :

Reference and fused images of IKONOS. (a) Original MS (R, G, and B) image; (b) original PAN image; and (c–h) the fused image using the AIHS, AWLP, SVT, ATWT, SR, and proposed method, respectively.70

Graphic Jump Location
Fig. 17
F17 :

PCNN and Shearlets framework for image classification.83

The contourlets were also used78 in the same context to merge SPOT and ALSAT-2A images. The authors advanced a conception of a new fusion scheme by combining the PCA and the NSCT in order to overcome the spectral distortion caused by the PCA. The aim of this work was to find a compromise between enhancing the spatial resolution and preserving the spectral information at the same time. After decomposing the MS using PCA, a histogram matching is applied to adapt the contrast between PAN and the resulting components. Then the obtained components and PAN image are decomposed using NSCT into approximation coefficients app and details coefficients det. The fusion rules applied to this approach are represented by both Eqs. (33) and (34), where authors extract the approximation coefficients of the fused image reconstim from the PAN image and the detail coefficients from the MS image.Display Formula

Imapp=apppan,(33)
Display Formula
Imdet=detpc.(34)

The resulting image is obtained using Eq. (35): Display Formula

reconstim=ImappImdet.(35)

Then the inverse of NSCT is calculated. The fused image obtained due to the proposed method, represented edges, contours of roads and buildings, and any structure shapes on the ground, better than other method, namely intensity hue saturation (IHS), PCA-IHS, and high pass filter (HPF). Moreover, the proposed method conserves the spectral information as well. Compared to PCA-NSCT based method, the resulting fusion image was not blurred.

Another injection model in pansharpening application is proposed in Ref. 79. Authors worked on QuickBird and IKONOS-2 imagery. The injection model is built on an adaptive cross gain, i.e., a ratio of local SD. Both images are decomposed using curvelets and then merged together by applying interband structure model. Compared to IHS, ATWT, and HPF, the proposed method exceeds them since it has the best rate according to the indicators: ERGAS, spectral angle mapper, and Q4. The resulted fused image is visually superior and succeeds in producing a tradeoff between different sensors.

The authors of Ref. 80 propose an approach for multisensor image fusion, based on beyond the wavelet transform domain (CT, bandlet, contourlet, and wedgelet). The approach consisted of the following steps: first, the authors decomposed the images into coefficients using beyond wavelet transform. Second, they selected from the two images the low frequency coefficients using maximum local energy (MLE) rule. It is calculated in a local 3×3 sliding window as shown in Eq. (36): Display Formula

LBEξl,k=E1*fξ(0)2+E2*fξ(0)2+EN*fξ(0)2,(36)
where fξ(0) is the beyond wavelet low frequency coefficients, E1,E2,,Ek are the filter operators in different directions, and l and k, respectively, are the scale and direction of the transform. They obtain the fused high-frequency coefficients using the sum modified Laplacian (SML) method [Eq. (37)]: Display Formula
SMLxl,k=i=MMj=NNML2f(i+p,j+q)for  ML2f(i,j)T,(37)
where ML2 is the modified Laplacian, xA or B are the source images, l and k are the scale and the direction of transform, respectively, and finally, M and N determine size of the chosen window. p and q are variables and T is a threshold. Finally, the fused image is obtained by performing an inverse beyond wavelet transform. Through a large number of experiments, they concluded that the fused images are best processed when using MLE-contourlet transform (the type of images was not mentioned).

The use of several transforms in an approach could increase the quality of fusion. Works like Zhong et al.81 and Zhanga et al.82 combined wavelet and curvelet together, to enhance wavelets’ abilities.

Neural networks have been also investigated in RS fusion. MGD were combined with this powerful tool in order to provide more directionalities. In fact, in Ref. 83, a fusion method of SAR images based on pulse couple neural network (PCNN)84 is proposed. The highlight of this algorithm is to use the global feature of source images and calculate the gradient information of them after being decomposed into several directions. Then wavelet is applied to decompose the images further while the high frequency coefficients are selected to be the input of PCNN to get a fire map. The max rule is applied to get the final fused image. The framework of the proposed method is detailed in Fig. 17.

Graphic Jump Location
Fig. 18
F18 :

Bandlet cloud removal: (a) image contaminated with clouds shadow, (b) mask of the area to be reconstructed, (c) proposed method,85 and (d) method Li et al.86

PCNN was also used in Ref. 87. The idea is to use the contourlet hidden Markov tree (CHMT) model88 to describe the statistical characteristics of contourlet coefficients of RS images. Besides, PCNN is used in this work in order to select the high-frequency directional sub-bands. First, contourlet transform is performed on the registered multisource SAR images and then the contourlet coefficients are trained by expectation maximization algorithm to calculate their edge PDFs. The highest magnitude in low-frequency sub-bands is selected. The high-frequency directional coefficients are updated by multiplying it by its edge PDF. Finally, a new clarity saliency measure is defined and used to fuse the high-frequency sub-bands. This method is compared to others: wavelets hidden Markov trees + PCNN and CT + PCNN. The results show that CHMT-PCNN can capture more directional information. Moreover, it needs less training time even with images abundant in textures. The criteria of assessment used in this work were the mutual information, weighted fusion (evaluates the visual effect), edge-dependent fusion quality index proposed by Ref. 89 (describes edge information in the image), and common information (represents the gradient information) proposed in Ref. 90.

Inpainting