## Abstract

We propose a new method for polarimetric synthetic aperture radar (PolSAR) imagery classification based on stochastic distances in the space of random matrices obeying complex Wishart distributions. Given a collection of prototypes and a stochastic distance *d*(.,.), we classify any random matrix *X* using two criteria in an iterative set-up. First, we associate *X* with the class which minimizes the weighted stochastic distance *w*_{m}*d*(*X*,*Z*_{m}), where the positive weights *w*_{m} are computed to max- imize the class discrimination power. Second, we improve the result by embedding the classification problem into a diffusion–reaction partial differential system where the diffusion term smooths the patches within the image, and the reaction term tends to move the pixel values towards the closest class prototype. In particular, the method inherits the benefits of speckle reduction by diffusion-like methods. Results on synthetic and real PolSAR data show the performance of the method.

## 1. Introduction

Classification is one of the most important techniques for image analysis. It aims at mapping each pixel into a class, so it transforms values into information.

Classification can be performed using a variety of sources, mostly the spectral information (the observed value in each pixel), spatial or contextual information and ancillary data (ground truth, for instance). The latter is usually available only in very restricted areas, from which training samples can be obtained.

The simplest available classification techniques rely only on pixel-wise information, i.e. on the observation in each coordinate: isodata, parallelepiped and pointwise maximum-likelihood (ML) are examples of these methods. Arguably, the most successful techniques exploit both the spectral information and the context. This is mostly due to the fact that images exhibit a great deal of spatial redundancy, i.e. spatially neighbouring pixels tend to be alike.

As an example of contextual classification one should mention techniques based on Markovian models. Geman & Geman [1] posed the classification process as an estimation problem and, as such, proposed a number of estimators and algorithms. These techniques rely on variations of the following idea: the class *m* in each coordinate should satisfy a criterion that, simultaneously, optimizes the pointwise spectral evidence and a contextual measure of smoothness, for example maximizing
1.1where *ξ*∈[0,1] is the relative weight of the spectral evidence over the context, *f*_{m}(*z*(*i*,*j*)) is the likelihood of the observation *z*(*i*,*j*) in coordinate (*i*,*j*) with respect to the model characterized by the probability density function *f*_{m}, *m*∈{1,…,*M*} is one of the *M* possible classes, and *N*(*m*,∂_{ij}) is a non-decreasing function on the number of neighbours of coordinate (*i*,*j*), denoted by ∂_{ij}, that have been classified as *m*. If *ξ*=0, only the context is relevant and the most frequent class in the neighbours is the optimal solution. If *ξ*=1, only the radiometric pointwise information is relevant, and the ML estimator (assuming independent observations) is the optimal solution. Specialized algorithms, such as the iterated conditional modes (ICM), maximum posterior modes (MPM) and simulated annealing, are required to obtain (often approximate) solutions for intermediate values of *ξ*; cf. [2]. Equation (1.1) makes it clear that the context, in this approach, is incorporated by means of the neighbouring classes.

The aim of this work is proposing a new classification procedure based on both context and radiometric information, but without the use of classes as descriptors of the former. For the latter, we tackle the problem of classifying polarimetric synthetic aperture radar (PolSAR) imagery.

According to Lee & Pottier [3], the early developments of imaging radar aimed at the characterization of aircraft targets. Spatial resolution was a fundamental issue with these images; it depends directly on the dimensions of the antenna, more precisely, on its *aperture*. Deployable antennas would yield unacceptable spatial resolutions of the order of hundreds of metres, even kilometres. Higher resolutions can be obtained by controlling the illumination, forming longer *synthetic* (as opposed to physical) antennas, yielding the broad area of *synthetic aperture radar* (SAR).

The main differences between SAR sensors and sensors that operate in the visible spectrum are

*Wavelength*: the former operate in the microwave section of the electromagnetic spectrum;*Activeness*: SAR sensors carry their own source of illumination;*Coherence*: SAR sensors are able to record the relative phase of the emitted and incidence signals.

These characteristics lead to a number of interesting and, oftentimes, challenging properties, among them

— SAR sensors are little affected by adverse weather conditions, such as fog, rain, smog, etc., and they provide images, regardless the presence of daylight;

— the return is mostly sensitive to the target dielectric properties, and to its geometry at the wavelength scale, which typically ranges from millimetres to metres;

— SAR images, after being processed, can be adequately described by a multiplicative rather than by additive model for the data.

Polarimetric SAR (PolSAR) extends the imaging ability of single SAR imaging. A PolSAR sensor transmits a linearly polarized signal and records two orthogonal polarizations of the returned signal. Such multidimensional information allows deeper investigation of the scattering mechanisms of the targets under study. Fully PolSAR images record the four possible combinations of the signal according to the polarization of the transmitting and receiving antennas. The first practical fully PolSAR sensor, the L-band AIRSAR, was developed by the Jet Propulsion Laboratory in 1985.

On the one hand, the deterministic approach has been used with success in the classification of PolSAR imagery, usually at the expense of using additional information. For instance, Shi *et al.* [4] perform manifold learning on layers of features derived from polarimetric decompositions, whereas Negri *et al.* [5] use neighbouring information and support vector machines.

On the other hand, statistical approaches able to cope with the departure from the usual Gaussian additive data formation provide an attractive framework for PolSAR image classification. Among others, references [6,7,8] provide a comprehensive account of the models, and of the relations among them, that have been used to describe polarized and fully polarimetric SAR images.

The classification of PolSAR imagery is an important research topic. Among the ones that use stochastic models for the fully polarimetric information, one should mention the pioneering work of Lee *et al.* [9], who present the multivariate complex Wishart distribution and some of its properties for the remote sensing community. Other techniques have been proposed as, for instance, the minimization of stochastic distances between segments of data and prototypes [10], and the use of spectral and contextual information by means of the ICM algorithm under the Potts model as prior distribution [2]. These works use the data as available to the users, differently from another approach, namely the classification based on the *H*/*α* decomposition [11]. It is noteworthy that the latter work presents a deterministic and unsupervised technique, whereas the others rely on the availability of prototypes and stochastic models.

Diffusion–reaction equations have been successfully applied in many different contexts, e.g. chemistry and ecology. In the field of image processing, diffusion–reaction equations have been applied, for instance in quantization [12], filtering and segmentation [13], and in medical applications [14].

The general shape of a diffusion–reaction partial differential equation is given by
1.2where *u*(*t*,*x*,*y*) represents the solution of the equation at instant *t*, and coordinates (*x*,*y*), is the diffusion term given, typically, by an elliptic second-order differential operator, for instance the Laplacian operator. The reaction term is *f*(*u*), whose shape strongly depends on the intended application. Roughly speaking, the reaction term tends to move the solution *u*(*t*,*x*,*y*) of the equation towards the nearest stable equilibrium state associated with the ordinary differential equation *u*′=*f*(*u*) (*u*_{m} is a stable equilibrium state of *u*′=*f*(*u*) if *f*(*u*_{m})=0 and *f*′(*u*_{m})<0). Usually, *u*(0,*x*,*y*) represents the observed data and *u*(*t*,*x*,*y*) represents its evolution under the action of the differential equation where the diffusion and reaction terms are in competition. On the one hand, the diffusion term tends to smooth the solution taking into account the neighbours values and, on the other hand, the reaction term tends to move the solution towards some asymptotic values.

In this work, we address a difficult problem: the classification of PolSAR images. Our approach is based on weighted stochastic distances and a diffusion–reaction partial differential equations system that imposes a certain level of spatial smoothness. This work is an extension of Gomez *et al*. [15], where we proposed a diffusion–reaction system, but using the Euclidean distance to perform PolSAR image classification. We will show that the use of weighted stochastic distances, instead of the Euclidean distance, improves significantly the classification results.

The rest of the article is organized as follows: §2 deals with the statistical description of PolSAR data, including expressions of stochastic distances between models. Section 3 is devoted to the classification of these data using weighted stochastic distances between complex Wishart distributions. Section 4 deals with the discretization of the diffusion–reaction differential system. Section 5 presents results on the classification of simulated and real PolSAR data, and finally, in §6, some conclusions are drawn.

## 2. The scaled complex Wishart distribution

This section is based on the work by Nascimento *et al.* [16].

### (a) Characterization

The polarimetric coherent information associates, in each frequency of operation, to each pixel a 2×2 complex matrix with entries *S*_{VV}, *S*_{VH}, *S*_{HV} and *S*_{HH}, where *S*_{ij} is the backscattered signal for the *i*th transmission and *j*th reception linear polarization, *i*, *j*=H,V. Under mild conditions, *S*_{HV}=*S*_{VH}, and the scattering matrix can be simplified into the three-component complex vector , where ^{⊤} denotes vector transposition. This random vector can be modelled by the zero-mean multivariate complex Gaussian distribution [17].

Different targets are characterized by different variances and, thus, are hard to discern visually. In order make such difference notable, and to improve the signal-to-noise ratio, *L* ideally independent measurements of the same target are obtained while processing the raw data. This quantity is the number of looks. These observations are used to produce the ‘multilook sample covariance matrix format’
where the superscript asterisk represents the complex conjugate transpose of a vector, and **s**_{i}, *i*=1,2,…,*L*, are the *L* scattering vectors. Assuming that these vectors are independent, * Z* is an Hermitian positive definite matrix and it follows a scaled complex Wishart distribution [18]. Having

**and**

*Σ**L*as parameters, the scaled complex Wishart distribution is characterized by the following probability density function 2.1where for

*L*≥3,

*Γ*(⋅) is the gamma function, tr(⋅) represents the trace operator, |⋅| denotes the determinant operator,

**is the covariance matrix associated with**

*Σ**,*

**s***Σ*=E(

**), where E(⋅) is the expectation. The first moment of*

**s****s***satisfies E(*

**Z***)=*

**Z****. We denote to indicate that**

*Σ**follows the scaled complex Wishart distribution. As can be seen in reference [19], this distribution is able to accommodate an arbitrary number of polarimetric components.*

**Z**### (b) Parameter estimation

Let {**Z**_{1},**Z**_{2},…,**Z**_{N}} be a random sample from . The ML estimators for ** Σ** and

*L*, namely and , respectively, are quantities that maximize the log-likelihood function associated with the Wishart distribution, which is given by 2.2where ,

**=[**

*θ*

*σ*^{⊤}

*L*]

^{⊤},

**=vec(**

*σ***), and vec(⋅) is the column stacking vectorization operator. Frery**

*Σ**et al.*[20] showed that is given by the sample mean 2.3and that satisfies the following nonlinear equation 2.4where is the zeroth-order term of the

*v*th-order multivariate polygamma function given by where

*φ*

^{(v)}(⋅) is the ordinary polygamma function expressed by for

*v*≥0, and

*φ*

^{(0)}(⋅) is the digamma function.

Nascimento *et al.* [16] noted that is biased, and proposed a number of correction techniques obtained subtracting an estimator of the bias from the ML estimator. Using their results, we employ the Box–Snell estimator of the bias, given by
2.5We first compute by solving (2.4), then we evaluate with (2.5) and, finally, we use the first-order bias-corrected estimator .

### (c) Stochastic distances between Wishart laws

Frery *et al.* [19], while seeking for contrast measures between samples, derived expressions for the stochastic distances between pairs of Wishart laws with (possibly) different covariance matrices *Σ*_{1} and *Σ*_{2}, and same number of looks *L*. They obtained the Kullback–Leibler (KL), Hellinger, Bhattacharyya and Rényi (of order *β*) distances.

The two first are of interest in this work, and they are given, respectively, by 2.6and 2.7Note that both distances depend on two very simple operations: the determinant and the inverse.

The two other distances, namely Rényi and Bhattacharyya, were not used, because the former depends on an extra parameter (*β*), and the latter is a simple increasing transformation of the Hellinger distance (HD), namely , so it adds little to this study.

## 3. Complex Wishart matrices classification using weighted stochastic distances

In this paper, we address a multiclass classification problem in , the space (or cone) of the of Hermitian positive definite matrices. All elements in are assumed to follow a scaled complex Wishart distribution, as defined by the density given in equation (2.1).

We assume there are *M* classes and, for each class *m*∈{1,…,*M*}, we denote by *Z*_{m} its prototype, given by a random matrix in that follows the scaled complex Wishart distribution . For each point of a PolSAR image, we have a random matrix *X* which is assumed to belong to one of the predefined classes.

A basic classification procedure is to assign each to the class *m* for which *X* is closest to *Z*_{m}, in a stochastic distance sense. This pointwise procedure is prone to errors owing to the variability of the observations. In order to improve this basic classification procedure, we introduce a collection of weights {*w*_{m}} such that *X* is associated with the class satisfying
3.1i.e. among the *M* possible classes, is the class whose prototype *Z*_{m} is closest to *X* with respect to a distance *d* weighted by a factor that depends on the class. If for any class, that class is forbidden; if *w*_{m}=0, the classification into *m* is mandatory. In §4, we will propose a procedure to estimate {*w*_{m}} in order to maximize the discrimination power of the classes.

These multiplicative factors have the effect of modifying the number of looks *L*_{m} of each class either directly, as in the KL distance, equation (2.6), or as a first-order approximation in the HD, equation (2.7). Their purpose is to exert control on the relative importance each class has on the procedure.

Each of the *M* prototypes is built by selecting a region of interest (ROI) that will serve as ground truth or reference. Each ROI is split into two disjoint sets of equal size by simple random sampling without replacement. The first set is used as training data, whereas the second is employed as test data. For each class, the training samples are used to estimate (*Σ*_{m},*L*_{m}) by improved ML, respectively, as described in §2b. Denote by *M*_{m} the number of pixels in the training set of class *m*, 1≤*m*≤*M*, and its pixels values by , 1≤*k*≤*M*_{m}.

The weights *w*_{m} are computed such that they maximize the discrimination power of the stochastic distance. This is done by using an energy optimization strategy, with the aid of the following energy function
3.2where
and λ>0 is a parameter. In order to simplify the optimization problem, we assume that . Roughly speaking, by minimizing equation (3.2), we make the observations be as close as possible to *Z*_{m}, and as far as possible to any other class representative *Z*_{m′}. The minimum of the above energy function is found using a gradient descent-type algorithm. In all the experiments showed in this paper, *w*_{m} is initialized as 1/*M* and λ is fixed to 1.

## 4. Diffusion–reaction differential system

This section follows the ideas introduced in reference [15], but using weighted stochastic distances instead of the Euclidean distance. In order to introduce a smoothing procedure in the classification step, we embed the classification problem in a diffusion–reaction differential system.

We denote by *Σ*(0,*x*,*y*) the initial covariance matrix provided by the PolSAR image, and by *Σ*(*t*,*x*,*y*) its evolution under the action of the diffusion–reaction system. Let *Z*(*t*,*x*,*y*) be the family of Wishart random matrices associated with *Σ*(*t*,*x*,*y*). The diffusion–reaction system we propose is given by
4.1where Δ*Σ* represents the spatial Laplacian differential operator applied on each component of the matrix *Σ*(*t*,*x*,*y*), {*Z*_{m}} represents distributions that describe the classes present in the PolSAR image, and *α*≥0 is a weight parameter that balances the influence of the diffusion and reaction terms, and is obtained using equation (3.1). The diffusion term, given by the Laplacian operator, tends to smooth the covariance matrix *Σ*(*t*,*x*,*y*), and the reaction term, given by the remainder of the equation, tends to move *Σ*(*t*,*x*,*y*) towards the nearest covariance matrix *Σ*_{m} according to the weighted stochastic distance *w*_{m}*d*(*Z*(*t*,*x*,*y*),*Z*_{m}).

### (a) Diffusion–reaction system discretization

To discretize equation (4.1), we denote by the approximation *Z*(*nδt*,*ih*,*jh*), and by its associated covariance matrix, where *i*, *j*∈*N*, and *δt* and *h* are the discretization values. We propose the following two-step algorithm where, in the first step, we discretize the contribution of the diffusion part of the equation and, then, in the second step, we discretize the contribution of the reaction part of the equation taking into account that, locally, using a linear approximation of the reaction term, the solution of the equation is given by an exponential function:

*Step 1*:
*Step 2*:
where represents the covariance evolution using only the diffusion operator, and its associated random matrix. Note that the first step incorporates the contextual evidence.

We observe that if 1−4*αδt*/*h*^{2}≥0, then is a convex combination of , , , , , and, in particular, . Moreover, because is a convex combination of and , then .

That is, if the original covariance matrices provided by the original PolSAR image belong to , then for all *n*>0 and, therefore, the associated random matrices are well defined. Note that the above discretization scheme depends on three parameters: *α*, *δt* and *h*. However, looking at the shape of the scheme, we point out that it depends just on the values *αδt*/*h*^{2} (step 1), and on *δt* (step 2); therefore, without loss of generality, we fix *h*=1. The only parameters that have influence on the results are *α* and *δt*. As explained above, in order to have a well-defined procedure with solutions in , these parameters must satisfy the condition 1−4*αδt*≥0; *δt* represents the time discretization step in the evolution equation.

The smaller *δt* is, the better becomes the approximation of the above scheme to the actual continuous solution of the differential system given by equation (4.1). In practice, as far as *δt* is small enough, we do not expect a strong influence of *δt* in the results. In practice, the most important parameter is *α*. The larger *α* is, the stronger the smoothing effect of the differential system given in equation (4.1) becomes.

In all the experiments presented in this paper, the discretization parameters have been fixed to *δt*=0.01 and *α*=0.5. Different choices and combinations are possible, and a detailed comparison study will be performed.

## 5. Experimental set-up

We compare the performance of the proposed methodology by classifying both simulated and real PolSAR images with respect to the following techniques:

(i) ML under the Wishart model: each observation is assigned to the class whose density, as expressed in equation (2.1), is maximized;

(ii) Euclidean distance (ED): each observation is classified into the class which minimizes the Euclidean distance between the observation and its prototype;

(iii) HD: each observation receives the label of that class which minimizes this distance to its prototype; cf. equation (2.7);

(iv) KL distance: same as above, but using equation (2.6);

(v) KL distance and optimized weights (KL+OW): same as in (iv), but using the optimized weights computed minimizing the energy functional (3.2);

(vi) Diffusion–reaction equation with KL+OW (

*DR*+KL+OW+*n*): the reaction–diffusion equation is solved using*n*iterations of the iterative discretization scheme, then each resulting observation is classified using*KL*+OW.

The visualization of PolSAR images requires stipulating a projection of the elements of into the unitary cube, then mapping these three components to the red, green and blue components of a colour image (the RGB representation). Lee & Pottier [3] discuss a number of such projections. We adopt here a simple visualization technique that emphasizes the quality of the classifications obtained. We assign a unique colour to each class representative *Σ*_{m}, then, we assign a colour to any using a linear interpolation procedure based on the Euclidean distance between *Σ*_{i,j} and *Σ*_{m}. The smaller the Euclidean distance between *Σ*_{i,j} and *Σ*_{m}, the larger the weight of the colour component of class *m* is in the interpolation procedure.

### (a) Simulated image

In the first experiment, we use a 300×300 pixels phantom image with three classes generated using three Wishart distributions. The parameters used for the simulation are those estimated in the training samples used in the experiment discussed in §4b. These training samples come from the ocean, forest and urban classes identified in figure 2. Each deviate from a Wishart law is obtained by simulating entries of the complex matrix entries, which obey complex Gaussian observations, then forming the sample covariance matrix and adding as many of them as the (integer) number of looks. This requires stipulating the number of looks, which in our case was *L*=4 for all classes, and the variance and covariances of the these entries [17] (see the electronic supplementary material for the obtained simulated data).

Table 1 presents the overall accuracy of classifications produced by the techniques here considered. The optimized weights, obtained through the minimization of the energy function presented in equation (3.2), were *ω*=(0.23,0.50,0.28). In all cases, the equivalent number of looks was informed *L*=4 and fixed for the three classes. For each class, the technique which provides the worst classification is used as baseline for computing the improvement each technique provides. These improvements are informed as percentages between parentheses.

All techniques classified all points of class 1 properly. The worst classifications for classes 2 and 3 were produced, respectively, by the Euclidean and by the KL distances. Class 3 presents the largest variability, hence improvement, in the results. The HD is the second best technique with a marginal improvement of 35.6% over the baseline. The Euclidean distance improves 36.7%. Introducing optimized weights improves the results of using the KL distance by 74.2%. Nevertheless, these classification procedures are outperformed by the ML rule, which is 87.3% better than the baseline. Albeit remarkably good, this last result is overshadowed by our proposal that attains 100% of accuracy, thus improving the baseline 100%.

Figure 1 shows these results. We include two images that show the evolution of the solution of the diffusion–reaction equation for *n*=25 and *n*=50 iterations. Table 1 also shows the CPU time (in seconds) for each classification method. We use a simple `C++` implementation in an Intel Core i5 (2.50 GHz) architecture without parallelization or other optimization techniques. We note that, even with this basic implementation, the methods are quite fast and in particular, the proposed method (DR+*KL*+OW+50) takes just 3.63 s.

### (b) Real image

In the second experiment, we use a PolSAR image obtained by the AIRSAR sensor over the San Francisco Bay.^{1} The image size is 600×448 pixels. Four classes are readily identifiable: ocean, grass, park and urban areas. We selected a ROI for each class, whose observations were used to estimate the Wishart parameters according to the results presented in §2b; these ROIs are shown in figure 2. These data are also used as the reference to estimate the accuracy of the different classification strategies.

Table 2 presents the overall accuracy of classifying the four classes observed in the San Francisco image by the techniques here considered. With the exception of the ocean class, that is almost perfectly classified by all methods, the worst result is used a baseline. Class-wide improvements with respect to the baseline are reported in parentheses.

The results for the real image are consistent with the ones obtained for the simulated data: Euclidean and KL distances provide the worst results. HD classification is the next best technique, except in the urban area where the Euclidean distance is slightly better. ML clearly outperforms Euclidean, KL and HDs, and even the weighted KL distance except in the case of urban areas where the weighted KL distance is considerably better (53.40%) than ML (39.18%). The optimized weights for this problem were *w*=(0.44,0.38,0.12,0.06). As expected, the urban area presents a high variability which is consistent with the small weight value (0.06) obtained for this area in the optimization procedure.

The proposed technique using the diffusion–reaction equation consistently improves in a significant way the classification accuracy (90.4%, 82.2% and 77.8% for the grass, park and urban areas, respectively), and when a perfect classification is obtained, as in ocean, the iterative procedure does not make it worse. These improvements are not incremental with respect to both the baseline and to the use of optimized weights, and render classified images that can be considered excellent. In figures 2 and 3, we illustrate the results of the different classification techniques for the San Francisco Bay image. In table 2, we also show the CPU time for the different classification methods. The results are consistent with the ones obtained for the phantom image taking into account the size of the images. In particular, the proposed method (DR+KL+OW+50) takes just 7.67 s.

Figure 4 shows the evolution of the diffusion–reaction system described by equation (4.1) observed in the classification of the San Francisco Bay PolSAR image. The ordinates of the two sequences are shown in logarithmic scale. As explained above, the reaction term of the diffusion–reaction system tends to move each PolSAR image pixel value *Σ*_{i,j} towards its nearest class representative ; in particular, we expect that the weighted distance average of *Σ*_{i,j} to goes to zero when the data evolves according to the diffusion–reaction system. This behaviour is illustrated in the continuous grey curve. Initially, such distance average is equal to, approximately, 4.38; we observe that it goes very quickly, faster than exponentially, to zero. In a complementary fashion, the diffusion term tends to smooth the solution and, as a consequence, it makes some pixels change their classification. The circles show the evolution of the percentage of pixels which change their associated class in each iteration (i.e. the value is altered). We observe that 1.9% of pixels change their associated class in the first iteration, and that this rate decreases across the iterations also faster than exponentially.

## 6. Conclusion

In this paper, we address the problem of PolSAR image classification. The main contributions of the paper are, on the one hand, to introduce a stochastic weighted distance strategy to increase the class discrimination power, and on the other hand, to embed the classification in a diffusion–reaction differential system that aims to include smoothing constraints in the classification procedure. The weights used to improve the class discrimination power of the distance are estimated by minimizing a new energy functional. These weights take into account the PolSAR matrix variability in the class regions (the lower the weight, the larger the expected variability inside the class).

We show that introducing these weights is equivalent to modifying the number of looks of the Wishart model for each class and, as discussed by Torres *et al.* [21], this takes account of the variability owing to texture that is not embedded in the Wishart model. We present experiments with simulated and real PolSAR images. These experiments confirm that the introduction of the weighted distances improves considerably the classification results.

By embedding the classification problem in a diffusion–reaction differential system, and using the weighted distances, the classification obtained with simulated and real data outperforms, in a significant way, the classification score obtained with the usual classification strategies including ML, and Euclidean, Hellinger and KL distances.

Further comparisons, including computational costs, will be made with respect to other contextual techniques as, for instance, the simple mode and the ICM algorithm [2]. Also, techniques for estimating optimal values of the parameters that govern the process (*α* and *δt*) will be sought.

## Competing interests

The authors declare that they have no competing interests.

## Funding

A.C.F. acknowledges financial support from CNPq and CAPES.

## Authors' contributions

Luis Gomez participated in the design of the diffusion-reaction model, the classification analysis and performing experiments simulating Wishart matrices; Luis Alvarez participated in the conception and design of the diffusion-reaction model and the associated discretization scheme, in the drafting of the manuscript and performing experiments; Luis Mazorra participated in the conception and design of the diffusion-reaction model and the associated discretization scheme; Alejandro C. Frery participated in the use of stochastic distances, in drafting the manuscript and in the design of the classification experiment. All authors gave final approval for publication.

## Footnotes

One contribution of 13 to a theme issue ‘Topics on non-equilibrium statistical mechanics and nonlinear physics (II)’.

↵1 Details about these freely available data can be accessed at ESA's web page https://earth.esa.int/web/polsarpro/data-sources/sample-datasets.

- Accepted July 17, 2015.

- © 2015 The Author(s)