THIS MANUSCRIPTIS A PREPRINT AND HAS BEEN SUBMITTED FOR PUBLICATION IN IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING 1 River Planform Extraction From High-Resolution SAR Images Via Generalised Gamma Distribution Superpixel Classification

The extraction of river planforms from remotely sensed satellite images is a task of crucial importance to many applications such as land planning, water resource monitoring or flood prediction. In this paper we present a novel framework for the extraction of rivers from Synthetic Aperture Radar (SAR) images, based on superpixel segmentation and subsequent classification. Superpixel segmentation is achieved via a modelling of the image pixels’ amplitudes and spatial coordinates as a finite mixture model, where the Generalised Gamma distribution is used to accurately model a variety of high-resolution SAR scenes. A number of features describing texture and statistics are extracted on a superpixel level, facilitating the identification of river superpixels planforms are then extracted via unsupervised, agglomerative clustering thus eliminating the need for labelled training data. We present results of our proposed method on ICEYE-X2 and SENTINEL-1 SAR data demonstrating its ability to produce pixel-accurate river masks.


I. INTRODUCTION
M ONITORING and analysis of river networks across the globe is a task of paramount importance, as they are an integral part of the Earth's water cycle and have the potential to severely impact ecosystems and human activities. Applications such as land use planning, landscape evolution [1] [2], sustainable water resource monitoring [3], flood monitoring & forecasting [4] [5] and even surveying compliance in multilateral environmental treaties [6] [7] are dependent on the accurate mapping of rivers and the monitoring of flooded areas, planform change and flow rates.
In-situ manual inspection of river networks can be extremely laborious even at small scales, and downright infeasible on a global scale as river gauge monitoring networks are far from ubiquitous and even when in place are typically sparsely distributed [8] [9]. However, technological advancements in the capabilities of remote sensing platforms and the ever increasing availability of remote sensing imagery create the potential for highly detailed, high temporal resolution largescale maps of river networks derived from such data [10] [11].
The wetted width of a river, its centerline, curvature and sinuosity can be inferred and/or measured from such data, in turn helping infer the morphodynamics of the river as well as estimate its discharge using methods such as at-many-stations hydraulic geometry [8]. Temporal monitoring of rivers via remote sensing data can allow for the study of bank erosion and/or accretion over time and provide meaningful insights on river channel evolution and migration mechanics. This ability to gain a better understanding of river evolution and the enhanced resource monitoring capabilities opened up via the use of remote sensing data has led to a demand for the development of algorithms for the automated extraction and interpretation of such data.
Algorithms and automated toolboxes that perform river segmentation and analysis have already started appearing in the literature [12]. These aim to produce river segmentation masks from remotely sensed images (most often optical and multispectral data), as well as produce a number of qualitative metrics such as river width estimates. RivMAP [13] is such a toolbox, capable of producing annual bankfull-resolving river masks from stacks of Landsat multispectral images, as well as quantify spatial and temporal changes over time. RivaMap [14] is another toolbox offering similar functionality, where a modified multiscale singularity index (MSI) [15] is used to detect rivers in Landsat images. RivWidth [16] is another popular software tool that uses spectral classification methods to produce river masks and that has also recently been redeveloped to integrate with the Google Earth Engine [17].
While multispectral data, especially those covering the Near-Infrared and Short-wave Infrared spectrum (NIR, SWIR), can be very convenient in that they inherently convey information that allows for pixel-level differentiation between land and water cover (and hence make river mask extraction an easier task via metrics such as the Normalised Difference Water Index [14]), they are also susceptible to adverse weather effects such as cloud cover, limited by the day/night cycle and tend to offer poor spatial resolution compared to other modalities. Synthetic Aperture Radar (SAR) is an alternative modality of particular interest to river monitoring [18]; its ability to operate regardless of meteorological conditions and throughout the day/night cycle, its high spatial resolution and the ability to differentiate between the surface roughness of land cover and water bodies make it a prime candidate for the task.
In this paper we are proposing an approach for extracting river masks in SAR images based on unsupervised superpixel classification. The use of superpixels is largely motivated by their ability to adhere well to object boundaries, congruent with our desire to accurately delineate the river's banks. A novel SAR-specific superpixel segmentation algorithm based on Generalised Gamma modelling of SAR images is utilised to oversegment the image into superpixels. We then propose a series of simple features that can be used to identify a superpixel as being a river superpixel based on its radiometric (statistical), structural and textural properties. River superpixels can then be merged together to produce a river planform mask; we demonstrate how this can be done accurately with a very simple, unsupervised agglomerative clustering method thus removing the need for labelled training data.
Our proposed method presents a number of novel features and improvements over other methods. The great descriptive power of the Generalised Gamma Distribution (GΓD) allows for improved superpixel segmentation performance, with better boundary adherence compared to previous methods utilising other statistical models. This is generally desirable for any superpixel-related application but even more crucial for applications such as river planform extraction, as the derived masks need to be as pixel-accurate and free of discontinuities as possible if they are to provide actionable information. This ability of the proposed superpixel segmentation algorithm to accurately delineate the river boundaries becomes particularly advantageous when considering other methods that rely on filtering operations which may inadvertently affect edge information in the image.
We also demonstrate the performance of our river mask algorithm over products from a number of SAR platforms, including state-of-the-art satellites such as the ICEYE-X2. These satellites, while they present furthered operational capabilities courtesy of their extremely high spatial resolution they also pose new challenges for the exact same reason. The increased resolution brings with it increased discriminatory power of individual scatterers, making for example for a more heterogeneous surface when imaging water surfaces [19] -the increasingly heavy-tailed nature of these state-of-the-art SAR products calls for ever more powerful heavy-tailed statistical models like the GΓD [20].
The remainder of this paper is structured as follows: Sections II discusses related work on river detection in SAR images and the associated challenges. Sections III & IV discuss preliminary material, the former including a discussion on superpixel segmentation and related algorithms and the latter discussing the Generalised Gamma Distribution (GΓD) and its use in SAR modelling. Section V describes the proposed Generalised Gamma Mixture Model superpixel algorithm. Superpixel feature extraction as well as clustering/classification are discussed in Sections VI and VII respectively. Finally, experimental results on ICEYE-X2 as well as SENTINEL-1 data are presented in Section VIII, with Section IX providing some closing remarks & conclusions, along with suggested avenues of future research.

II. RELATED WORK
Rivers in SAR images appear as dark longitudinal, curvilinear structures, with comparatively smooth surface texture. Bodies of water in general tend to appear as smooth regions in SAR images as the radar returns are much smaller in amplitude and relatively uniform compared to those caused by hard scatterers on land. These radiometric,structural and textural characteristics have been used in algorithms seeking to identify and segment rivers in SAR images, often borrowing ideas from similar applications such as road network detection, fingerprint analysis and biomedical vessel detection.
Thresholding methods are commonly applied for SAR river extraction given the radiometric properties of water surfaces in SAR images. An example of such an approach can be found in [21] where Otsu's thresholding is used as the starting step for a morphological process for river channel extraction. While thresholding can in principle be very efficient, in practice the performance of such methods is largely dependent on the amount of speckle noise present. Such methods are also susceptible to the presence of any other object with similar radiometric response as well as to the presence of discernible clutter on the water surface. These concerns become particularly poignant when dealing with high-resolution state of the art SAR sensors, as their increased resolution leads to higher appreciability of individual scatterers on the water surface thus lending water bodies a rougher texture in the SAR image.
Structural methods that look for longitudinally continuous curvilinear dark objects in SAR images have been proposed for a number of detection tasks, including that of rivers. Sun and Mao [22] segment rivers in SAR images based on their structural characteristics following a concept borrowed from fingerprint identification that involves performing edge detection in the wavelet domain followed by ridge tracing to enforce connectivity. Yang et al. [23] use Gabor filters to first enhance river cross-sections then followed by a pathopening morphological operation to enforce connectivity of longitudinal segments. Their method is demonstrated on Landsat panchromatic images but the principle could be further generalisable to SAR imagery.
Variants of active contours have also been proposed [24], utilising for example cross entropy [25] as well as adaptive global fitting energies [26]. Matgen et al [27] use an active contour model seeded by low backscatter pixels in the image. The use of active contour models can produce very good results, with accurate delineation of the river boundaries, however they do so with an important caveat -the appropriate selection of initial seeding points. Thresholding operations have been used by some of the above methods to automatically provide seed points for active contour models, however as discussed previously those can in practice suffer greatly from the presence of speckle.
The textural profile of river surfaces in SAR images has been the basis for the work of Senthilnath et al. [28] who use gray level cooccurrence matrix (GLCM) texture analysis and mean shift segmentation to extract river masks for the purpose of flood monitoring. Wavelet energy is used to analyse texture and differentiate between land and river in [29]. Textural information combined with structural information also form the basis for the river extraction algorithm of [30], which uses the standard deviation of the Structural Feature Set (SFS-SD) along with morphological path-opening operations to identify narrow and elongated areas of homogeneous texture.
One of the most important aspects of any river detection algorithm is its ability to accurately delineate the river bank in the image. River planform masks are useful not so much in their own right but in that they facilitate measurements and/or estimates of a number of quantities of interest such as river centerline & width, channel migration, area of accretion/erosion over time, discharge rates etc. SAR image products typically have a spatial resolution measured in metres (if not tens of) per pixel. Therefore any ambiguities in the actual boundaries of the river can significantly impact the quality of actionable information that can be deduced from the mask. Any smoothing effect as a result of a filtering operation (such as a decimated wavelet filter, a morphological operator etc) can negatively impact the usefulness of the produced mask. Our proposed method addresses this particular problem of accurate edge extraction by performing river detection via superpixel segmentation.

III. SUPERPIXEL SEGMENTATION
Superpixel segmentation algorithms aim to oversegment an image into homogeneous groupings of perceptually similar pixels, with the aim of utilising these segments as meaningful primitives for further image processing tasks. Superpixels have in the past been used for object segmentation & classification, target detection & tracking and other tasks in a variety of image data types and applications.
Superpixel methods are often graph-based, as is one of the earliest proposed methods using normalised cuts [31] and its many variants that followed. One such state of the art variant is the Linear Spectral Clustering (LSC) algorithm of Chen et al. [32] which follows a normalised cuts-based problem formulation but also incorporates an RGB colour similarity & spatial proximity distance metric.
Another popular approach for superpixel generation is clustering; SLIC [33] is one popular algorithm that uses a modified k-means clustering algorithm to iteratively group pixels together based again on their spatial proximity and colour similarity (this time as expressed within in the CIELAB colorspace). For a general overview of some of the most commonly utilised superpixel segmentation methods the reader is directed to [33].

A. Superpixels for SAR images
Superpixel segmentation has also been studied specifically for the context of SAR images. SAR images differ significantly from natural images; colour information is (typically) absent, the noise follows a multiplicative rather than additive model and the intensity/amplitude statistics are generally far from Gaussian [34] [20]. Superpixel algorithms developed for natural images are therefore often far from optimal for segmenting SAR images.
Efforts to adapt existing superpixel methods such as SLIC for use with SAR images include examples such as [35], where the Euclidean distance between a pixel and its possible superpixel centre is replaced with the Amplitude Ratio Distance of [36]. In the case of polarimetric SAR data, the additional channel information can be used in lieu of colour information; this is the case in [37] which uses an iterative clustering technique similar to SLIC but with utilising polarimetric information instead of colour and the Wishart distance instead of the Euclidean. Finally, Hu et al. [38] present an edgedominated SLIC-inspired iterative clustering algorithm where a metric penalising the crossing of an edge (the edge map being derived via a degenerate filter) is included in the distance metric.

B. Mixture Model Superpixels
Of particular interest here is the Mixture Model Superpixel algorithm (MISP) first proposed by Arisoy and Kayabol in [39]. MISP performs segmentation based on a finite mixture model (FMM) which takes a pixel's amplitude and spatial coordinates as features to be modelled. Amplitudes are modelled according to the Nakagami distribution and spatial coordinates via the normal distribution.
Two pixels belonging to the same cluster are then assumed to have their amplitudes generated from the same Nakagami distribution (the same mixture component) and their spatial coordinates from the same normal distribution. An iterative process is then employed to assign superpixel labels to each pixel in the image and to infer the related statistical parameters for the FMM.
For completeness we note here that the Nakagami probability density function is given by where µ and v are here the scale and shape parameter of the Nakagami distribution respectively and Γ(·) denotes the Gamma function. The Nakagami, along with distributions like the Rayleigh and the Rician, are commonly used to model signals in situations where a receiver picks up signals arriving from multiple different paths due to scattering. While the Nakagami has been used in the past as a theoretical model for multilook SAR images, it is not quite a state of the art model fit for accurately modelling modern, high-resolution, highly heavytailed SAR data [34]. The Nakagami distribution's theoretical validity has also come under question in the case of high-resolution modern data as these exhibit non-Gaussian behaviours that undermine the Gaussian conjecture involved in its derivation, as well as the multiplicative model involved [20]. More descriptive, highly generalisable models that can more adequately model the highly heavy-tailed nature of modern SAR data have, instead, been proposed in the literature [20], [40], [41].
MISP serves as the basis for the Generalised Gamma Mixture Model superpixel algorithm proposed here, and the algorithm will therefore be presented in more detail in a later section. Our proposal revolves around the use of the Generalised Gamma distribution instead of the Nakagami for modelling SAR images within the MISP framework, as it provides a more capable fit for state-of-the-art, high-resolution SAR images [20].

IV. THE GENERALISED GAMMA DISTRIBUTION
The Generalised Gamma Distribution (GΓD) [42] has over the past years been established as one of the most flexible and accurate empirical statistical models for high resolution SAR images [20]. The great descriptive power of the GΓD allows it to model both amplitude and intensity SAR images, something not always feasible with previously used models. It is also more capable of dealing with the ever-increasingly heavy-tailed nature of modern SAR products, as well as modelling a large variety of land cover topologies [20].
The probability density function of the GΓD is defined 1 as: where v ∈ R =0 is the distribution's power parameter which determines tail behaviour, and κ, σ ∈ R + are the shape and scale parameters respectively. Γ(·) denotes the Gamma function.
While the Gamma distribution is characterised by its shape κ and scale parameter, σ, the additional parametrisation offered by the inclusion of the power parameter v allows the GΓD to more flexibly control the shape of the pdf, making it particularly capable at modelling heavy-tailed histograms as the value of v decreases [20]. Examples of how the shape of the pdf can vary under different combinations of parameter values can be seen in Figure 1.

A. Parameter Estimation via the Mellin Transform and secondkind statistics
A number of parameter estimation methods for the GΓD have been proposed in the past, including the second-kind statistics method used in [20]. This approach is capable of producing simpler characteristic equations, moment-and cumulant-generating functions than the classical approach of analysing a pdf via the Fourier transform, particularly when dealing with complicated distributions. These in turn allow for explicit parameter estimators that are computationally efficient and easy to use [43].
The Mellin transform has been proposed [43] as a more effective alternative to the Fourier transform specifically for analysing a distribution defined over R + . In analogous fashion to how a characteristic function is defined as the Fourier transform of a probability density function, the second-kind first characteristic function φ 1 (s) is defined as the Mellin transform of the probability density function p(x) and further the second-kind second characteristic function φ 2 (s) is given by the natural logarithm of the first one: Continuing with the analogy, the i-th order second-kind momentsm i (also referred to as log-moments) can be deduced by derivation of the first characteristic function: and the i-th order second-kind cumulantsζ i (also referred to as log-cumulants) 2 can be deduced by derivation of the second characteristic function: For the case of the GΓD specifically, it can be shown (see [20] for proof) that the log-cumulant expression of (6) can be equated to:ζ where Φ 0 (x) and Φ 0 (x, y) here refer to the Digamma and Polygamma functions respectively. This allows for estimating the parameters κ, v and σ by solving the equations (7) and (8) against the calculated sample log-cumulants ζ 1 , ζ 2 and ζ 3 . Using a second-order approximation of the Polygamma function Li et al. [20] show that finally the shape parameter κ can be estimated via the following closed-form equation where the quantities P, Q are defined as: and λ i as: Finally, obtaining an estimate of κ allows for obtaining closed form estimates of v and σ from (7) and (8) aŝ

B. Generalised Gamma Mixture Models
The use of the GΓD in SAR image modelling has also been extended to finite mixture models (FMMs), with the Generalised Gamma Mixture Model (GΓMM) proposed in [44] aiming to model high-resolution SAR images with ever greater accuracy. This more recent work further supports the notion that SAR images can be very efficiently modelled using GΓMM, obtaining more accurate models than in the case of e.g. Nakagami mixture models.
While the image being modellable by an GΓMM is an underlying assumption of the GΓD-MISP algorithm, during the course of the algorithm only parameter estimation methods for the GΓD are required -this will become evident in the next section where the algorithm is presented in detail. As direct estimation of the parameters of a GΓMM is not necessary in the context of the superpixel algorithm proposed it will not be further discussed here -the reader is directed to [44] and [45] instead for relevant literature on GΓMM parameter estimation.

V. GENERALISED GAMMA MIXTURE MODEL SUPERPIXELS
The Mixture Model Superpixel algorithm uses a pixel's amplitude a n and spatial coordinates q n = [x n , y n ] T as features to be statistically modelled; we denote the feature vector of the nth pixel in the image as f n = [a n , x n , y n ]. Note that here n = [1, ..., N ] is the lexicographically ordered pixel index with N being the total number of pixels in the image. The image is to be segmented into a total number of K mutually exclusive superpixels, with the kth superpixel denoted as S k where k ∈ [1, ..., K].
Note that this value K is not explicitly defined by the user. The user sets a desired RegionSize value governing the superpixel size and the algorithm initialises the superpixel map as an equispaced grid of RegionSize * RegionSize tiles. The initialised value K then becomes ImageW idth/RegionSize * ImageLength/RegionSize .
The entirety of the image is assumed to be modelled by a finite mixture model (FMM), with the amplitudes (and conversely the spatial coordinates) of two pixels belonging in the same superpixel assumed to be generated from the same amplitude (and conversely spatial) distribution, i.e. the same mixture component. The number of components in the mixture is therefore the same as the generated number of superpixels. The amplitudes and spatial coordinates are assumed to be statistically independent.
In [39] the authors opt to model the amplitudes of pixels with the Nakagami distribution, which is a commonly employed multilook amplitude model for SAR images. Here however we opt for the GΓD (2) and hence a pixel's amplitude a n is modelled as where v k , κ k , σ k are the GΓD power, shape and scale parameters respectively that model the kth superpixel. The spatial distribution of pixels around a superpixel centroid is assumed to be normal (bivariate Gaussian) and is given as follows where the random variable q n as discussed corresponds to the (x, y) spatial coordinates of the nth pixel of interest, m k is the kth superpixel's centroid spatial coordinates and Σ k is the covariance matrix of that same superpixel.
The kth superpixel can thus be described by the parameter set θ k = {v k , κ k , σ k , m k , Σ k }, i.e. the parameters of the fitted amplitude (GΓD) and spatial (bivariate Gaussian) distributions for the pixel population included in the superpixel S k .
Under conditional probability rules, the joint density of a pixel's feature vector f n and label vector z n can be equated to p(f n , z n |θ 1:K , ω 1:K ) = p(f n |z n , θ 1:K )p(z n |ω 1:K ), and it can be further shown [39] that the following finite mixture density can be obtained where ω k corresponds to the mixture proportion of the superpixels. A conjugate Dirichlet prior for these proportions can then be defined as where α is the concentration parameter and B(α) = Γ K (α)/Γ(αK).

A. Inference
Having described this probabilistic model structure, it is clear that there now needs to be a procedure for inferring a number of variables; these are the parameter sets θ 1:K describing the amplitude & spatial distributions for each superpixel, the mixture proportion ω 1:K of each superpixel and finally the superpixel label z 1:N of each pixel. These can be estimated using the block iterated conditional mode algorithm (ICM) [46], which allows for conditional densities of variables to be maximised iteratively.
As outlined in [39] this allows the variables to be updated along iterations, in the following order: where t is the iteration/pseudo-time index. From these the close-form parameter-specific update equations can then be worked out, with the mixture proportion worked out as and the bivariate Gaussian parameters for the spatial distribution as Parameter estimation for the GΓD is performed according to the second-kind cumulants method outlined in the previous section, with the closed forms shown in equations (9), (12) and (13). Note that in the case of the GΓD 3 parameters need to be estimated, rather than 2 for the Nakagami case -the computational efficiency and robustness of the estimator in [20] however helps mitigate this possible computational overhead.
The algorithm initialises on an equispaced grid and terminates after a fixed number of iterations, similarly to MISPthis was here experimentally set at 20 iterations, as further iterations showed no significant change in the produced segmentation map. Figure 2 shows the MSIP-GΓD segmentation produced after 5, 20 and 80 iterations of the algorithm. By the 5th iteration the superpixels are already largely adhering to the edges of strong features in the image (such as the river banks) while the overall map is still reminiscent of the initialised equispaced grid. The superpixel edges corresponding to the river edges undergo very little change after the 20 th iteration. This can of course be experimented on depending on the application/data but is not anticipated to require modification.
At the end of the process any generated segments that are smaller than a defined threshold (1/20 th of RegionSize) are merged with their closest adjacent neighbour, as is standard practice. This is primarily to avoid any erroneously generated single-pixel superpixels in the segmentation map and should not require application specific tuning.

VI. SUPERPIXEL FEATURE EXTRACTION
In order to extract river planforms an automated method for classifying superpixels as land cover or water cover is necessary, hence the need to identify features that can be used to discriminate between these two classes of interest. Features of interest here tend to fall in one of three categories, these being: statistical features that reflect the radiometric difference between water and land cover (i.e. river superpixels appearing darker than land superpixels); textural features that reflect the smoother texture of river superpixels compared to land superpixels; and finally features based on structural and morphological filters. Note how these are analogous to the approaches for river detection discussed previously in Section II.
Features here are meant to be extracted out of (and be descriptive of) the pixel population within a superpixel, and not a sliding window as is often the case. This leads to the rather unique position where features are meant to be extracted from sample populations for whom there is some a priori knowledge available [19]. Specifically, all pixels belonging to a superpixel can be generally assumed to belong to the same object, as per the object boundary adherence property and content homogeneity property of the superpixels. As a consequence, certain features may be more difficult to define and calculate over a spatially irregular superpixel (compared to rectangular window) while other features may be far more meaningful when calculated over a homogeneous superpixel population (rather than the unknown population of an ad hoc sliding window).

A. Statistical Measures
The statistics of water-cover superpixels differ greatly from those of land-cover superpixels due to the differences in the radiometric response of the water surface versus that of the land clutter. Rivers appear as near-uniformly dark features in SAR images, and therefore river superpixel histograms can be expected to consist predominantly of low intensity values.
Accordingly, measures of central tendency can be employed to identify such superpixels of predominantly low intensity values. The sample median is a simple yet clear indicator for distinguishing river superpixels on this basis. It is also more robust in the presence of speckle noise (extreme values) compared to other statistical measures describing central tendency such as the sample mean. Another statistical indicator of central tendency readily available, is the scale parameter σ of the GΓD, which is estimated via (13) for each superpixel in the course of the MISP-GΓD algorithm. A graphical representation of how the parameter σ relates to the central tendency of the GΓD can be seen in Figure 1 (c). Figure 3 shows the histograms and statistics of a sample river and land superpixels, illustrating the difference described above.

B. Texture Features
Texture analysis in SAR images, for the purpose of image segmentation, is a typical procedure with the use of many classic texture features & analysis methods already investigated in the relevant literature [47].
At a most basic level, texture can be described by firstorder (occurrence) metrics such as mean, variance and most commonly entropy. Entropy is typically discussed in this context as a metric calculated over a small sliding window across the image, aiming to capture texture variations in small local neighbourhoods. With the superpixel segmentation in place however entropy can be calculated over the entirety of the superpixel population and serve as a direct indicator of "disorder" within it. A perfectly uniform superpixel (in terms of gray-levels, i.e. one with a single-peak histogram) would have zero entropy. A river superpixel, with its almost smooth dark content, will thus have a much lower entropy value compared to land-cover superpixels.
For a random population X with possible values {x 1 , x 2 , ..., x i , ...x n } the entropy value S(X) is given by Texture can also be characterised by second-order (cooccurrence) metrics, the most common of those being the Gray Level Co-Occurrence Matrix (GLCM) features. GLCM is statistical method that seeks to define how often pixel pairs of a certain spatial and gray-level relationship appear in the image. As the GLCM is therefore a table of frequencies its usefulness actually lies in the ability to calculate a number of features from it, such as GLCM-specific entropy, homogeneity, and contrast metrics [48]. Such GLCM features have been used to analyse SAR texture with the purpose of identifying and delineating rivers and flood boundaries [28] [49].
While for the purposes of the present work information entropy was found to be an adequate descriptor of texture, we note that more complex features, such as GLCM, could also be employed for superpixel regions. Superpixel regions of interest should be accordingly eroded prior to the calculation of the GLCM, effectively excluding border pixels to avoid contamination with external pixels. This is particularly important as a superpixel's edge is highly likely to be in alignment with an object edge in the image, and therefore, comparing a superpixel's edge pixel with its immediate external neighbours can severely affect the texture descriptors.

C. Filter Output Features
Finally, we wish to address the use of a number of filterderived detection methods that have been employed for river detection (or similar tasks), and how these can be reconciled within the superpixel-based framework discussed here. By filter-derived we here refer to methods such as the modified multiscale singularity index (MSI) of [15], the Gabor filters and path-opening morphological filters of [23], or the wavelet energies of [29].
Unlike a statistical moment, a histogram or an information theoretic measure, which can simply be evaluated specifically over a superpixel population's intensity values, these methods fundamentally rely on neighbourhood operations (i.e. filtering operations) that can prove problematic within the irregular spatial structure of a superpixel. Instead, one can employ such methods as originally proposed over the entirety of the image, and then apply the superpixel segmentation map over the detection map produced by such methods. The average filter response within each superpixel can then serve as a superpixellevel feature.
This leads to the rather circular proposition of performing river detection via superpixel classification which is in turn facilitated via river detection algorithms. Namely, if one were to employ a river detection method such as MSI why not simply take its output as the river mask, which is the intended use of the method in the first place, and instead employ it as one of many features in this superpixel detection framework?
This latter approach offers a number of benefits. It allows for circumvention of one of the problems often inherent in filtering methods, that being any edge blurring or smoothing effects. The river mask boundaries are directly formed from the superpixel boundaries, which themselves are meant to adhere as accurately as possible to the actual river bank in the image and thus result in an ideally pixel-accurate river mask. This can be of particular importance if for example the river width is to be estimated from such a generated mask. In the MSI method of [15] for example the authors propose an estimate for the river width based on interpolation between the scale with the highest singularity index response and its neighbour scales, as direct pixel-wise calculation on the MSI output would not be accurate.
Within this superpixel framework however, one can employ such methods more as a confidence value for the presence of a river within a superpixel rather than a direct detector, without having to worry about the effects of any filtering operation on edge preservation. As an example, we employ the mean MSI response per superpixel as one such feature. For completeness, we briefly describe the modified multiscale singularity index as described in [15] for the detection of river networks. For each scale σ M SI the image is first debiased by convolving the image with a Gaussian filter G with standard deviation σ M SI and subtracting it from the original: The algorithm then computes the index at location q n of the debiased image as M SI(q n , σ M SI ) = |f 0,θ M SI ,σ M SI (q n )f 2,θ M SI ,σ M SI (q n )| 1 + |f 1,θ M SI ,L M SI σ M SI (q n )| 2 (27) where here θ M SI is the direction orthogonal to the curvilinear mass, estimated by the second order derivatives of the input image calculated along evenly spaced directions at the beginning of the algorithm, f 0,θ M SI ,σ M SI , f 1,θ M SI ,σ M SI and f 2,θ M SI ,σ M SI are respectively the zero-first-and second-order Gaussian derivatives at scale σ M SI and direction θ M SI (q n ) and finally the scalar L M SI is a scaling constant with value 1.7754 [50]. This index is evaluated over all n = [1, 2, ..., N ] scales, and with the filter window size for each scale set at 6σ . The maximum value across all scales for any given pixel location is set as the filter output. Figure 4 demonstrates the problem of edge preservation under filter-based detectors. As an example, the output of the MSI method is shown for two SAR images. This method is by design not aiming to produce pixel-accurate masks,focusing instead on the ability to detect river centerlines over a variety of scales, and thus detecting multi-channel river networks of highly-variable channel widths in a single image. While superpixel and filter based techniques are generally used for different purposes, when utilised within a superpixel framework filter methods can provide a useful confidence metric on the presence of a river within the image, effectively combining the detecting power of filter algorithms with the edge preservation benefits of the superpixel-based approach.
Note that the inclusion of such filter-derived features is by no means necessary; a superpixel's population alone can provide sufficient information for classification, for example via the statistical and textural measures previously discussed. We here opt to combine elements from all 3 categories mentioned above as an illustration of the breadth our method.

VII. SUPERPIXEL CLASSIFICATION
A river planform mask can now be produced by classifying the superpixels into river superpixels and land superpixels. As is typical with classification tasks, this can be achieved in a number of ways broadly described as either supervised classification, where an algorithm is trained via labelled data, or as unsupervised clustering, where no labelled training data are required and instead the algorithm tries to group instances samples together based on similarity.
The previous section made specific note to how features of interest were selected to be particularly descriptive of water surface superpixels and to help differentiate them from land and other cover. This was done in the interest of facilitating unsupervised classification. Besides removing the need for labelled data and classifier training, such clustering methods provide great performance with minimal complexity. They are also quite an intuitive fit to the task of putting together river superpixels to form a river mask based on their inherent similarity.
Specifically, we employ agglomerative hierarchical clustering to group image superpixels into two clusters representing land-cover and river superpixels. The algorithm iteratively pairs superpixels (and then groupings of superpixels) based on the similarity of their feature vectors, until all superpixels have been grouped into two classes. Choosing the pairs of clusters to be merged on each iteration is done according to Ward's minimum variance method [51] and the distance metric used is the Euclidean distance.

VIII. EXPERIMENTAL RESULTS AND DISCUSSION
In this section we demonstrate the performance of the proposed river planform extraction method. This is demonstrated over a number of scenes derived from two real SAR images described below.

A. Datasets
Image 1 has been acquired by the ICEYE-X2 SAR satellite. ICEYE-X2 is a state-of-the-art X-band SAR satellite, having been only launched in late 2018 and is capable of providing incredibly high spatial resolution (down to 1m for Spotlight products). The scenes shown here measure 300 x 300 pixels and are derived from a VV-polarised Stripmap product (Level 1B, Ground Range Detected, Orthorectified, Intensity image) with a spatial resolution of 3m. The product has been multilooked, with 3 looks in the azimuth and 2 in the range direction. The image shows part of the Cổ Chiên River and its tributaries near the city of Trà Vingh, Vietnam and was acquired on March 1st, 2019. Image 2 has been acquired by the SENTINEL-1B SAR satellite. The scene shown measures 600 x 900 pixels and is derived from the VH-polarised band of an Interferometric Wide-Swath (IW) product (Level 1, Ground Range Detected, Amplitude image) with a spatial resolution of 10m. The product has been multlooked, with 1 look in the azimuth and 5 in the range direction. The image shows part of the Moselle river north of Metz, France near the Luxembourg border and was acquired on October 27th, 2019.

B. Superpixel Segmentation
We begin the discussion by looking at the performance of the proposed MISP-GΓD superpixel algorithm. Figure 5 shows two scenes from an ICEYE-X2 image and the superpixel segmentation maps produced by the SLIC, LSC, MISP and MISP-GΓD algorithms. Both the MISP and MISP-GΓD algorithms here have the superpixel size (RegionSize parameter) set at 20x20 pixels and the concentration parameter α set at 10 6 . We found RegionSize parameter values of 15 to 30 pixels to work quite well for the data shown here. Smaller values run the risk of creating superpixels too small for robust estimation of GΓD parameters, while significantly larger superpixels do not bring any advantages, as a river may well be spanned by multiple superpixels if its width in the image is larger. Having the RegionSize parameter set at 20 does not mean the algorithm cannot delineate river with a width under 20 pixels, as the superpixels can evolve into more oblong, narrower shapes. To obtain comparable results the LSC algorithm was also given a superpixel size value of 20x20, while SLIC was set to have a desired output of 225 superpixels (to equate to an average superpixel size of 20x20) and the compactness value was set to 20.
As discussed previously, most algorithms developed for superpixel segmentation of natural images struggle when dealing with SAR data, partly due to the lack of colour information. The characteristic speckle noise and irregular texture is also problematic in such images. This is evident in the SLIC and LSC examples. While the algorithms produce reasonable groupings of uniform pixels in the river regions, the adherence to the boundary is not particularly accurate. The overall segmentation is also far from uniform, particularly over the land cover with both algorithms creating non-convex superpixels of significantly varying size. More importantly for our application, the boundary adherence of the superpixels along the river bank is not on par with MISP and especially the proposed MISP-GΓD algorithm. Figure 6 illustrates the above point on boundary adherence. While the original MISP algorithm has been shown to perform very well in a variety of scenarios, the improved descriptive power of the Generalised Gamma Distribution utilised in our variant allows for even better performance in challenging cases such as the ones shown in Fig. 6. These correspond to parts of the imaged river where noticeable mid-channel gravel bars are present and whose accurate segmentation can prove challenging in the case of the original MISP algorithm. This is an improvement that may seem marginal for a general application, but can be crucial for our particular application of river planform extraction. Considering that the spatial resolution of many wide-swath SAR image products is often in the tens of metres per pixel, a single pixel may create a significant change in the area or river width estimation. The underlying superpixel segmentation needs to be as accurate as possible if the produced planform masks are to be of any practical use in a hydrological context. Improper boundary adherence can furthermore cause discontinuities in the river mask, as mixed superpixels comprised of both land and water surface risk being misclassified as land superpixels. Figure 7 shows examples of the proposed river planform extraction method on ICEYE-X2 high resolution SAR images. As before the MISP-GΓD algorithm has the superpixel size (RegionSize parameter) set at 20x20 pixels and the concentration parameter α set at 10 6 . The features used for the following examples are the sample median, the GΓD scale parameter σ, the information entropy and the mean MSI filter response per superpixel.

C. River Planform Extraction
The algorithm is capable of identifying all river superpixels in the image as belonging to the same class, thus producing a river planform mask for the scene. Of particular note is the lack of discontinuities and the ability to delineate river boundaries accurately, courtesy of the accuracy afforded by the MISP-GΓD superpixel segmentation.
Note how the scene in Fig. 7(j) includes two river channels of noticeably different widths. The maleability of the MISP-GΓD algorithm allows it to generate superpixels of appropriate size for each of the two channel structures within the image even though the superpixel region size parameter remains constant. The desired superpixel size should be set by the user upon some consideration of the spatial and pixel resolution of the image in conjunction with the size of the desired features to be detected. However, as this example illustrates, this does not need to be exact as the algorithm is capable of modifying the shape of the generated superpixels to conform to river channels of varying width.
False positives can occur in the masks when there are superpixel regions in the image that appear uniformly dark and hence get clustered along with the river superpixels. These could be small lakes, partly visible river channels (as in Figure  7 (i)) or dark patches of land. Removal of these disjointed false positives can be relatively straightforward, using standard connected component analysis methods. Component size, major axis length and circularity are properties that can be used to distinguish between isolated superpixels and the river channel in the mask. In the future we aim to expand on more detailed, clever heuristics for this task. Figure 8 shows a further example on SENTINEL-1 data. The superpixel segmentation parameters and classification features are as of Figure 7. To illustrate the removal of false positives we include a scene where many lakes are visible. These false positives are filtered out from the final river mask by thresholding based on a connected component's major axis length, i.e. by thresholding out components that do not exhibit the high aspect ratio expected of river segments.
Note also that a particularly thin meander channel fails to be picked up completely, as it becomes too fine to be accurately segmented by the superpixel algorithm. The width of the river at that section is around 3 pixels. While the superpixels are capable of squeezing into narrow channels of such width (as seen at the edges of this particular channel that have been detected) they may run into problems over prolonged stretches of narrow channels. Such narrow river segments may (as in this case) be absorbed into the nearest land-cover dominated superpixel and thus fail to be picked up in the classification stage.
It should finally also be noted that if large scenes pose a computational challenge then the segmentation and characterisation could be performed in batches. Doing so may artificially introduce false positives, as the algorithm will produce detections over each tile even when a river channel is not present in that part of the image (as for example in the upper right corner of the image of Figure 8).

D. Performance Comparison
Finally, we demonstrate the problem of accurate edge extraction that is associated with filter-based methods. Figures 4  and 9 shows examples of the filtering operations involved in other river detection methods. Examples of the output of the MSI filter [15] were shown previously in Figure 4 (b) & (e). This kind of approach is not able to produce pixel-accurate masks as far as the river edges are concerned -its strength is the ability to detect river network channels of highly-variable width in a single image. When utilised within a superpixel framework, it can provide a useful confidence metric on the presence of a river (or river-like structure) within the image without any edge degradation/blurring effects affecting the end product mask.
Sghaier et al. [30] base their method on texture measurement via the Standard Deviation feature of the Structural Feature Set. The output of this SFS-SD measurement and the binary Fig. 9. Examples of river detection using the SFS-SD method [30] and the Gabor filtering/Path Opening algorithm [23]. The response of the SFS is highly dependent on the manual setting of parameters controlling spectral and spatial thresholds -this affects how accurately delineated the river banks will be in the end image. The filter can also be quite sensitive to the texture within the river when dealing with high-resolution data where clutter is discernible over the water surface. Their proposed method later employs morphological operations (specifically path openings) for the enforcement of connectivity and the removal of false positives, a process that can again affect edge information in the image. Figure 9 further shows an example of the river detection method described in [23]. A combination of spectral filtering, Gabor filter banks and path-opening/closing operations is used to detect rivers, an approach that can be applied to a great variety of remote sensing image types. These filtering operations again introduce significant edge blurring, which may not be a major issue when dealing with large scenes of low spatial resolution like Landsat data but would be a problem in the case of highly-detailed SAR data.

IX. CONCLUSION
In this paper we have proposed a method for producing accurate river planform masks out of SAR images using superpixel characterisation and subsequent unsupervised clustering. Central to our proposed method is MISP-GΓD, a novel superpixel segmentation algorithm that utilises the Generalised Gamma distribution to model SAR images. We show how MISP-GΓD achieves improved accuracy of segmentation compared to other methods, thereby leading to more accurate delineation of the river planform in the final mask. Generated superpixels are then clustered based on features describing their statistical and textural properties. We demonstrate the results of our method on SENTINEL-1 as well as ICEYE-X2 SAR products. A current MATLAB implementation can be found at: https : //github.com/odispap/SARSuperpixelRivers.
We aim to expand on this body of work in future by further enhancing the accuracy of the generated river planform masks and by incorporating more elaborate heuristic methods for the identification of river channels and the removal of false positives, making for example the algorithm more capable of handling the presence of bridges. We also aim to focus on providing methods for reliably extracting measurements from the river masks; we hope to make a toolbox providing such functionality available in the near future.