TEXTURAL FEATURES FROM WAVELETS ON GRAPHS FOR VERY HIGH
RESOLUTION PANCHROMATIC PLEIADES IMAGE CLASSIFICATION

Minh-Tan Pham¹, Grégoire Mercier¹, Julien Michel²

1: TELECOM Bretagne, UMR CNRS 6285 Lab-STICC/CID, 2900 Brest CEDEX 2, France

2: CNES, DCT/SI/AP - BPI 1219, 31401 Toulouse CEDEX 09, France

Résumé

Dans cet article, nous proposons une méthode de caractérisation locale des textures des images de très haute résolution spatiale, dans lesquelles l'hypothèse de stationnarité est peu respectée. Une approche ponctuelle (i.e. non-dense) est d'abord introduite pour la représentation de l'image en utilisant un ensemble de pixels d'intérêt au lieu de la totalité des pixels de l'image. Un graphe pondéré est ensuite construit à partir de ces pixels représentatifs. Le signal de texture, porté sur ce graphe, est ensuite analysé à travers une transformée en ondelettes sur graphe. La classification en texture, implémentée ici de façon non-supervisée, est donc réalisée par la classification des coefficients d'ondelettes sur le graphe. Les expérimentations appliquées aux images panchromatiques Pléiades nous donnent des résultats très prometteurs avec une bonne précision de classification tout en gardant une compléxité intéressante.

**Mots clés :** Ondelettes sur graphes, Image Pléiades de très haute résolution spatiale, Représentation éparse d'image, Texture, Classification non-supervisée d'image

Abstract

The aim of this paper is to focus on local texture characterization for very high resolution remote sensing images. To prevent from the lack of stationarity hypothesis, the image is represented though a pointwise approach (i.e. non dense) by using sparse representative pixels. A weighted graph is then constructed from these representative pixels, and a wavelet transformed is applied on the graph to extract textural features. A classification algorithm, implemented here in a non-supervised way, of the wavelet coefficients achieves the texture-based classification. Experiments applied on VHR panchromatic Pléiades images provide very promising and encouraging classification results since they seem to be accurate while remaining of low complexity with low memory requirements.

**Keywords :** Wavelets on graphs, Very high spatial resolution Pléiades image, Sparse image representation, Texture, Unsupervised image classification

1. Introduction

Texture characterization for segmentation and classi-
fication is one of the most crucial tasks of image proces-
sing which have been tackled in many research studies.
However, with the emergence of very high spatial reso-
lution (VHR) imagery technique, this task becomes more
challenging when texture zones to be classified are about
to be too small or they do not respect the stationarity hy-
pothesis. Therefore, our motivation in the present work is
to develop a texture-based image classification algorithm
based on non-dense local feature characteristics. Hence,
a graph-based approach is investigated as it seems to
be probably relevant for this kind of image modeling and
since no stationary condition is required.

Signal processing on graphs has become an emer-
ging field with several applications in many diverse do-
mains (Shuman et al., 2013). One of the most outstan-
ding tools of this field is the spectral graph wavelet trans-
form (SGWT) which was introduced in (Hammond et al.,

2011). Recent works have studied and used the SGWT
to tackle a lot of problems from various areas such as
(Tremblay and Borgnat, 2013), (Leonardi and Ville, 2011).
In the domain of image processing, SGWT has also been
exploited to attack the problem of image denoising, pro-
posed in (Hammond et al., 2010), or to develop a fil-
terbank algorithm for an edge-aware methodology (Na-
rang et al., 2012). Here, we investigate the use of graph
wavelet transform to characterize local texture features
and then perform an unsupervised classification algo-
rithm applied for the VHR panchromatic Pléiades images.

The paper is organized as follows. Section 2 gives
a brief review of signal processing on graphs in which
some basic notions of graph and the SGWT approach
are going to be presented. Then, the main proposed me-
thodology is described. In Section 4, a short introduction
of VHR panchromatic Pléiades image datasets is given
before experimental results are performed for a validation
of the proposed algorithm. We finally draw a conclusion
from this work as well as consider some perspectives for
future work in Section 5.

<footnote>This work is supported by the CNES grant and the Région
Bretagne grant.</footnote>
## 2. Review of Signal Processing on Graphs

### 2.1. Weighted graph and characteristics

A weighted graph $G = \{V, E, w\}$ consists of $|V| = N$ vertices connected by a set of edges $E$ with associated weights $w$ involving the similarity between vertices. The adjacency matrix $A$ of $G$ is given as follows:

$$A_{ij} = \begin{cases} w(i, j) & \text{if } (i, j) \in E \\ 0 & \text{otherwise} \end{cases} \quad (1)$$

The graph Laplacian matrix is defined as $L = D - A$ where $D$ denotes the degree matrix with $D_{ii} = \sum_j A_{ij}$. These matrices are squared with $N$-by-$N$ dimension, $D$ is a diagonal matrix while $A$ and $L$ are symmetric. Spectral graph analysis is performed based on the eigen decomposition of $L$, given by $\{L\chi_k = \lambda_k\chi_k\}_{k=0,...,N-1}$ (Chung, 1997). Since $L$ is symmetric, positive and semi-definite, $\{\lambda_k\}_{k=0,...,N-1}$ provides a non-negative eigenvalue set with $0 = \lambda_0 < \lambda_1 \le \lambda_2 \le \dots \le \lambda_{N-1}$ and $\{\chi_k\}_{k=0,...,N-1}$ gives an orthogonal eigenvector basis.

### 2.2. Graph wavelet transform

Using a function $g: \mathbb{R}^+ \to \mathbb{R}^+$ representing a transfer function of a band-pass filter in the spectral domain of graph, Hammond et al. (2011) generated the wavelet $\psi_{t,n}$ localized in vertex $n$ at scale $t$ as follows:

$$\psi_{t,n}(m) = \sum_{k=0}^{N-1} g(t\lambda_k)\chi_k^*(n)\chi_k(m) \quad (2)$$

Thus, wavelet coefficients of a function $f \in \mathbb{R}^N$, called function on vertices, at scale $t$ and vertex $n$ is defined as the projection of $f$ onto the wavelet basis:

$$W_f(t, n) = \langle \psi_{t,n}, f \rangle = \sum_{k=0}^{N-1} g(t\lambda_k) \hat{f}(k) \chi_k(n) \quad (3)$$

where $\hat{f}(k) = \langle f, \chi_k \rangle = \sum_{n=1}^{N} f(n)\chi_k^*(n)$ represents the defined graph Fourier transform based on the expansion of $f$ in terms of graph Laplacian eigenvectors.

Similarly, to be analogous to classical wavelet transform, scaling function $\phi_n$ and corresponding coefficients $S_f(n)$ can be also generated with the support of a low-pass filter $h$ in graph spectral domain as follows:

$$S_f(n) = \langle \phi_n, f \rangle = \sum_{k=0}^{N-1} h(\lambda_k) \hat{f}(k) \chi_k(n) \quad (4)$$

Hence, given a function $f \in \mathbb{R}^N$, SGWT generates a set of scaling function coefficients $S_f \in \mathbb{R}^N$ and $J$ sets of wavelet coefficients $W_f(t_j) \in \mathbb{R}^N$ at scale $t_j$, hence performs a multiscale analysis of signal on graphs.

### 2.3. Illustration

Figure 1 shows an example of multiscale decomposition for a region of interest (ROI) from a VHR image by graph wavelets performing on a dense graph constructed by setting each pixel as a graph vertex. Each vertex is connected to its 8 closest neighbors on the image plane.

<figure><img src="image_1.png" /><figcaption>FIGURE 1: Example of multi-scale image decomposition using wavelet transform on a dense graph</figcaption></figure>

Edge weights are then calculated based on the similarity between pixel intensities and can be performed by an exponential form $e^{-||I(i)-I(j)||}$, where $I(i)$ and $I(j)$ denote the intensities of vertices $i$ and $j$. Scaling function image (i.e. lowpass components) is found on the top right of the figure while two wavelet images (i.e. bandpass components) at 2 different scales $t_1$ and $t_2$ are found at the bottom.

<figure><img src="image_2.png" /><figcaption>FIGURE 2: Corresponding low-pass and band-pass filters for SGWT image decomposition in Figure 1</figcaption></figure>

We also display in Figure 2 the corresponding low-pass and band-pass filters in spectral domain which were used to generate those scaling function and wavelet coefficients according to Equations 4 and 3, respectively.

## 3. Methodology

Our proposed method for texture-based VHR panchromatic Pléiades image classification is found in Figure 3. Each of process blocks is addressed in this section.

### 3.1. Extraction of representative pixels

The first proposition from this paper is to extract a set of pixels for a pointwise image representation which means to use a set of pixels of characteristics to represent image contents instead of using all image pixels like dense approaches. The problem of missing or losing information may be encountered, especially for describing and analyzing image textures. However, with an acceptable
<figure><img src="image_3.png" /><figcaption>FIGURE 3: Proposed workflow for texture-based panchromatic Pléiades unsupervised image classification using SGWT</figcaption></figure>

level of information loss, such a sparse approach becomes interesting and prospective in terms of memory cost and processing time when dealing with huge size data which are favored in our case of VHR Pléiades images.

Several approaches can be considered to extract representative pixels from an image. The simplest way could be to spatially downsample the image by a given level. Others can be considered are SIFT (Scale Invariant Feature Transform) and SURF (Speeded Up Robust Features) operators which have already achieved successful performance in the scope of image recognition and image matching. In this work, we consider that textural features could be sparsely described by using the local maximum and minimum pixels from the image.

<figure><img src="image_4.png" /><figcaption>FIGURE 4: Texture representation and distinction by local maximum and minimum pixels</figcaption></figure>

Figure 4 provides an example of different texture zones from a VHR image. Considering two distinctive textures localized in the red and green circles, 1D signal was obtained at the bottom by the yellow trace. From this 1D model, we compare the first texture with high variation to the second one which is smoother and less varied. Our

proposition is to characterize and distinct these textures using only the local maxima (red points) and local minima (black points) by taking into account their amplitudes (i.e. pixel intensities), their relative distance (i.e. pixel spatial position) and how they arrange to form the texture. This principle leads to our point-wise approach for texture representation and characterization.

A pixel is supposed to be a local maximum (or local mininum) if it has the greatest intensity (or lowest intensity, respectively) within a defined search window centered at it. Denote $S_{\omega \times \omega}^{\text{MAX}}$ the set of maximum pixels with search window size $\omega \times \omega$, we have:

$$ (i, j) \in S_{\omega \times \omega}^{\text{MAX}} \Leftrightarrow \left\{ I(i, j) = \max_{(k, \ell) \in N_{\omega \times \omega}(i, j)} I(k, \ell) \right\} \quad (5) $$

where $I(i, j)$ denotes the intensity value of pixel $(i, j)$ and $N_{\omega \times \omega}(i, j)$ represents a set of neighboring pixels of $(i, j)$ within the search window.

The set of local minimum pixels $S_{\omega \times \omega}^{\text{MIN}}$ could be also extracted using similar principle from Equation 5. These sets of local pixel extrema are then used for sparsely describing local features and constructing weighted graphs.

## 3.2. Graph construction

One of our motivations is to construct a weighted graph connecting representative points from the image in which edge weights involve the measure of vertex similarity in terms of local texture features. As it plays a significant role in the present approach, we propose dividing this stage into two small steps: graph vertex definition and description followed by graph weighted edge creation.

**Vertex definition and description.** Graph vertex set $V$ are formed by the extracted representative pixels. We propose that $V = S_{\omega \times \omega}^{\text{MAX}}$. As a consequence, graph vertex density in the processed image depends on the parameter $\omega$ which indicates the level of sparse representation to be performed. The smaller $\omega$ is, more pixels are extracted and the constructed graph thus becomes denser.

Next, for each vertex, a vector of signatures which is used for texture description and then for similarity measures between graph vertices is created. The idea is to describe local feature properties represented by a set of $K$ closest local maxima and a set of $K$ closest local minima around each vertex. To do that, we propose calculating different kinds of measures associating the variation of intensities, distances and directions formed by the vertex itself and each of two mentioned extrema sets.

For a given vertex $p(x, y, v)$ having its intensity value $v$ at image position $(x, y)$ and a set of $K$ closest maxima (or minima) $\Upsilon_K = \{p_k(x_k, y_k, v_k); k = 1, \dots, K\}$, following features are computed. Note that each pair $(p, p_k)$ creates a geometrical vector on the image plane. Corresponding length and angle of each geometrical vector are denoted by $d_k$ and $\theta_k$, respectively.

— Mean and variance of intensities:

$$ \mu_v = \frac{1}{K} \sum_{k=1}^{K} v_k, \quad \sigma_v^2 = \frac{1}{K} \sum_{k=1}^{K} (v_k - \mu_v)^2 $$
— Mean and variance of spatial distances:

$$ \mu_d = \frac{1}{K} \sum_{k=1}^{K} d_k, \quad \sigma_d^2 = \frac{1}{K} \sum_{k=1}^{K} (d_k - \mu_d)^2 $$

$$ \text{where } d_k = \sqrt{(x - x_k)^2 + (y - y_k)^2} $$

— Measures of α-dispersion and concentration of directions, presented in (Mardia and Jupp, 2000):

$$ D_{\alpha} = \frac{1}{K} \sum_{k=1}^{K} \{1 - \cos(\theta_k - \alpha)\}, \quad \bar{R} = \sqrt{\bar{C}^2 + \bar{S}^2} $$

$$ \text{where } \bar{C} = \frac{1}{K} \sum_{k=1}^{K} \cos \theta_k \text{ and } \bar{S} = \frac{1}{K} \sum_{k=1}^{K} \sin \theta_k \text{ and } \alpha = 0^{\circ} \text{ in our implementation.} $$

Hence, denote $\delta(i)$ is the computed vector of signatures for vertex $i$, $\delta(i)$ gives a description of local textural features of the environment around $i$. Edge weights can be then computed based on a similarity measure between these vectors.

**Graph weighted edge creation.** The computation of graph edge weights enables us to build the graph adjacency matrix. Given that each vertex has itself a description vector, edge weights associated to the vertex similarity can be defined by two general ways. One is to use a threshold to preserve only strong edges between vertices while the other is to connect each vertex to its k-closest neighbors in terms of distance measure between their description vectors (i.e. k-strongest edges). The second is chosen in this work since it seems to be simpler to control the number of edges from each vertex than to fix a threshold for eliminating weak edges. So, edge weights are computed as follows:

$$ w(i, j) = \begin{cases} e^{-\gamma[\text{dist}(\delta(i), \delta(j))]^2} & \text{if } j \in \mathcal{N}_k(i), \quad (6) \\ 0 & \text{otherwise} \end{cases} $$

where $\text{dist}(\delta(i), \delta(j))$ represents the distance measure between description vectors of two vertices $i$ and $j$. For the actual work, the Euclidean distance $||\delta(i), \delta(j)||_2$ is used in which these vectors $\delta(i), \delta(j)$ are normalized so that each measure of intensities, distances and directions has the normal distribution $N(0, 1)$. $N_k(i)$ represents the set of k-closest neighbors of $i$ in terms of description vector distances and $\gamma$ is a free parameter fixed to 1 in our work.

### 3.3. Multiscale analysis using SGWT

A weighted graph $G$ with associated adjacency matrix $A$ has been built from the previous section. Multiscale analysis using SGWT approach is now effectuated. As reviewed in section 2, graph wavelet transform is mainly performed in the spectral domain via the eigen decomposition into eigenvalues and eigenvectors of the Laplacian matrix $L$. Given that such a decomposition requires huge memory and time consumption, in particular when facing with large-size graphs, Hammond et al. (2011) thus proposed a fast computation using the truncated Chebyshev polynomial approximation. The authors also confirmed that this approximation is strongly efficient in case

of sparse Laplacian matrix which is appropriate to our previous method of graph construction.

Implementation of SGWT requires the design of low-pass and band-pass, or even high-pass, filters in graph spectral domain. In their works, Hammond et al. (2011) designed and used a set of cubic splines for wavelet kernels (i.e. band-pass filters) and an exponential waveform capturing low frequencies for scaling function. Example of these filters has been shown in Figure 2. They also allowed the reconstruction of input signal for the application such as signal synthesis or compression. Here, it should be noted that our application mainly involves graph vertex classification by manipulating wavelet coefficients. Hence, the choice of filter kernels is not bounded by any constraint of signal reconstruction. This advantage gives us more possibilities when selecting a successive of filters to capture expected frequencies along the axis and then produce desirable wavelet coefficients.

Now, given the adjacency matrix $A$ of graph and a number of scales $J$, following steps are activated for the implementation:

+ Calculate the graph Laplacian matrix $L = D - A$ and estimate the upper bound $\lambda_{max}$ of the graph spectral domain. Note that we do not need to compute all the set of eigenvalues and eigenvectors by direct decomposition of $L$.

+ Calculate the set of scales $\{t_j\}_{j=1}^J$ which will decide the levels of stretching filter kernels on frequency domain. Then, design the set of filters comprising one low-pass $h(.)$ and $J$ band-pass $g(t_j, \cdot)$ corresponding to each scale $t_j$.

+ Compute the truncated Chebyshev polynomial approximation for the scaling function and wavelets as describing in (Hammond et al., 2011).

+ Compute all SGWT coefficients consisting of a set of scaling function coefficients $S_f(n)$ and $J$ sets of wavelet coefficients $W_f(t_j, n)$ by supposing that the function on vertices $f(n)$ is the intensity value at vertex $n$.

Hence, with a set of $N$ representative pixels linked by a weighted graph, SGWT analysis generates $(J + 1)$ set of coefficients which encapsulate the notions of local features as well as the similarity between graph vertices at different scales. The next stage employs these coefficients for classification task.

### 3.4. K-means clustering

We effectuated an unsupervised clustering algorithm by performing the K-means (Hartigan and Wong, 1979) on the normalized wavelet coefficients. We note that after the transform, wavelet coefficients at each scale vary in a different range from the others. In particular, scaling function coefficients often possess large values while high-frequency wavelet coefficients vary in a range close to zero. Therefore, the normalization is required so that all coefficient values comply with the distribution $N(0, 1)$.

One advantage of using K-means classifier is that only the number of classes is required. Therefore, during the clustering process, users only need to import their expected number of classes.
# 4. Data sets and experiments

## 4.1. Data sets

Process datasets in our implementation are panchromatic images acquired by the VHR Pléiades-1A Constellation (PHR-1A) using its high spatial resolution panchromatic sensor of 70cm at nadir. The images were then resampled at 50cm resolution. They are available from www.astrium-geo.com.

## 4.2. Experiments

We performed the experiment on a crop of 1500 × 2000 pixels from a panchromatic Pléiades image acquired in Melbourne, Australia in February 25, 2012 by the PHR-1A. The input image is shown in Figure 5.

Figure 6 shows the extracted local maximum pixels using a sliding research window with size 11 × 11 pixels. As a result, 14588 pixels were extracted and used as representative pixels. Consequently, a weighted graph with $|V| = 14588$ vertices is constructed to connect these pixels for a sparse image representation.

Next, for the graph construction (Section 3.2), following parameters were applied. We set a number of 20 closest local maxima and minima around each vertex. Then, only 200 strongest edges were kept for each one. SGWT analysis (Section 2.2) was implemented by exploiting the cubic spline filters with a number of $J = 3$ scales. Fast computation algorithm using the truncated Chebyshev polynomial approximation was also applied.

In order to evaluate the performance of our proposed method, we would like to perform, in addition, the unsupervised K-means algorithm on the extracted Haralick textural features using the Gray Level Cooccurrence Matrix (GLCM) technique (Haralick et al., 1973). This GLCM approach was implemented as follows:

+ For each representative pixel, compute 4 cooccurrence matrices at 4 different directions (0°, 90°, 45° and 135°). GLCM parameters: window size 65 × 65 pixels and coexisting pairwise pixel distance equal to 2.

+ For each matrix, extract 5 Haralick textural features including Contrast, Correlation, Energy, Homogeneity and Entropy. Thus, a description vector consisting of $4 \times 5 = 20$ Haralick features is formed.

+ Normalize these vectors so that each feature has the normal distribution $N(0, 1)$.

+ Perform the K-means algorithm on these normalized feature vectors.

By setting the number of classes of the K-means algorithm equal to 4 for both SGWT and GLCM approach, classification results are obtained and displayed in Figures 7 and 8, respectively, for a comparison.

Visually, our proposed method shows better performance in terms of characterizing and classifying different texture zones according to the chosen vertex set from the image. We find the thematic coherence and robustness dealing with the image illumination condition. On the contrary, the GLCM approach mixes the classes and is not capable of providing a general discrimination between different classes. For example in Figure 8, there

often exists an intermediate class (red) in the transition between golf zone (blue) and forest (green) but we could not define its proper nature. Also, the small lake which was well classified in Figure 7 (red) could not be distinguished in Figure 8 (a mix of yellow, green and red). Moreover, it should be also noted that the result of GLCM approach significantly depends on the choice of its parameters. Here, we tried to perform one of its good results among several of experimental tests.

In terms of time consumption, a comparison can be found regarding to two Tables 1 and 2. All the tests were performed by a desktop computer Xeon 3.6GHz using Matlab. Before comparing the two tables, we would like to emphasize that the GLCM method is usually performed for a dense approach in the literature. That means, Haralick feature vector will be calculated for each pixel. We had also tested this dense GLCM approach and since our image consists of 3 million pixels, it took approximately 2 hours for the implementation. Here, for a comparison, we only show computational time of the proposed method and the GLCM method performed only on selected pixels (14588). Our method requires more computational time than the GLCM approach due to the stage of weighted graph construction. However, as a total time of 251.81s was taken for a complete algorithm with better classification performance, the proposed method is still considered to be efficient in terms of time consuming.

<table><caption>TABLE 1: Computation time for unsupervised classification using the proposed texture-based SGWT method</caption><thead><tr><th></th><th>Computation time</th></tr></thead><tbody><tr><td>Extraction of representative pixels</td><td>17.59 s</td></tr><tr><td>Construction of weighted graph</td><td>230.90 s</td></tr><tr><td>Fast SGWT</td><td>3.27 s</td></tr><tr><td>K-means</td><td>0.05 s</td></tr><tr><td><strong>Total time</strong></td><td><strong>251.81 s</strong></td></tr></tbody></table>

<table><caption>TABLE 2: Computation time for unsupervised classification using the non-dense GLCM method</caption><thead><tr><th></th><th>Computation time</th></tr></thead><tbody><tr><td>Extraction of representative pixels</td><td>17.59 s</td></tr><tr><td>Extraction of Haralick features</td><td>35.16 s</td></tr><tr><td>K-means</td><td>0.17 s</td></tr><tr><td><strong>Total time</strong></td><td><strong>52.92 s</strong></td></tr></tbody></table>

# 5. Conclusion

A novel texture characterization method for very high resolution images is proposed by using the spectral graph wavelets. In the current paper, the algorithm is applied for unsupervised classification of panchromatic Pléiades images. Using such a method, we would like to emphasize the effectiveness of a point-wise approach by a graph model for texture representation and description. During the paper, the proposed approach has been performed
<figure><img src="image_5.png" /><figcaption>FIGURE 5: A crop from a VHR panchromatic Pléiades image, © CNES (2012), distribution Airbus DS/Spot Image</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 6: Vertex definition from the extracted representative pixels (i.e. local maxima)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 7: Unsupervised classification result using the proposed SGWT approach</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 8: Unsupervised classification result using the GLCM approach</figcaption></figure>

using the local maximum and local minimum pixels to characterize textural features, escaping classical dense approaches when dealing with texture description for very huge size images. Graph wavelets are generated on our own graph and their coefficients are then exploited. Experiments have performed very promising and competitive results in terms of good texture classification as well as low time consumption.

Future work can concentrate on the development of the algorithm in terms of improving each block of the processing chain, especially for the extraction of representative pixels and the definition of vertex description vectors for graph construction. Other perspectives could be to investigate the present method on other types of very high resolution image data such as Radar or multispectral images.

Références

Chung, F. R. K., 1997. Spectral graph theory. In: AMS.

Hammond, D., Raoaroor, K., Jacques, L., Vandergheynst, P., April 2010. Image denoising with nonlocal graph wavelets. In: SIAM Conference on Imaging Science.

Hammond, D. K., Vandergheynst, P., Gribonval, R., 2011. Wavelets on graphs via spectral graph theory. Appl. Compt. Harmon. Anal. 30 (2), 129–150.

Haralick, R. M., Shanmugam, K., Dinstein, I., November 1973. Textural features for image classification. IEEE Trans. on Systems, Man. and Cybernetics 3 (6), 610–621.

Hartigan, J. A., Wong, M. A., 1979. Algorithm as 136 : A k-means clustering algorithm. Journal of the Royal Statistical Society. Series C (Applied Statistics) (1), 100–108.

Leonardi, N., Ville, D. V. D., March 2011. Wavelet frames on graphs defined by fMRI functional connectivity. In: IEEE International Symposium on Biomedical Imaging.

Mardia, K. V., Jupp, P. E., 2000. Directional statistics. John Wiley and Sons Ltd.

Narang, S. K., Chao, Y. H., Ortega, A., 2012. Graph-wavelet filterbanks for edge-aware image processing. In: Statistical Signal Processing Workshop. pp. 141–144.

Shuman, D. I., Narang, S. K., Frossard, P., Ortega, A., Vandergheynst, P., May 2013. The emerging field of signal processing on graphs : Extending high-dimensional data analysis to networks and other irregular domains. IEEE Signal Processing Society 30 (3), 83–98.

Tremblay, N., Borgnat, P., September 2013. Multiscale community mining in networks using spectral graph wavelets. In: EUSIPCO.