Reconstruction of River Boundaries at Sub-Pixel Resolution : Estimation and Spatial Allocation of Water Fractions

Boundary pixels of rivers are subject to a spectral mixture that limits the accuracy of river areas extraction using conventional hard classifiers. To address this problem, unmixing and super-resolution mapping (SRM) are conducted in two steps, respectively, for estimation and then spatial allocation of water fractions within the mixed pixels. Optimal band analysis for the normalized difference water index (OBA-NDWI) is proposed for identifying the pair of bands for which the NDWI values yield the highest correlation with water fractions. The OBA-NDWI then incorporates the optimal NDWI as predictor of water fractions through a regression model. Water fractions obtained from the OBA-NDWI method are benchmarked against the results of simplex projection unmixing (SPU) algorithm. The pixel swapping (PS) algorithm and interpolation-based algorithms are also applied on water fractions for SRM. In addition, a simple modified binary PS (MBPS) algorithm is proposed to reduce the computational time of the original PS method. Water fractions obtained from the proposed OBA-NDWI method are demonstrated to be in good agreement with those of SPU algorithm (R2 = 0.9, RMSE = 7% for eight-band WorldView-3 (WV-3) image and R2 = 0.87, RMSE = 9% for GeoEye image). The spectral bands of WV-3 provide a wealth of choices through the proposed OBA-NDWI to estimate water fractions. The interpolation-based and MBPS methods lead to sub-pixel maps comparable with those obtained using the PS algorithm, while they are computationally more effective. SRM algorithms improve user/producer accuracies of river areas by about 10% with respect to conventional hard classification.


Introduction
The extraction of river areas is the primary task required for a wide range of remote sensing applications in fluvial systems spanning from hydrological, ecological, and morphological studies to mapping habitat suitability for different aquatic species [1][2][3][4][5][6].Thresholding on spectral bands, such as the near-infrared (NIR) band, or on water indices, such as the normalized difference water index (NDWI), as well as image classification (supervised or unsupervised), are the main techniques for delineation of water bodies using optical imagery [7,8].However, most of the techniques consider a hard labeling approach for producing the water mask.This means that mixtures within the pixels are considered very roughly so that each pixel can be assigned to only a single water/non-water class which represents the most abundant class within the pixel.The spectral mixture can occur at every spatial resolution, particularly in the boundary pixels [9,10].This point deserves more attention in terms of remote mapping of rivers, as the boundaries of river channels are inevitably subject to uncertainties concerned with mixture of water with surrounding land covers.Mixed boundary pixels can suppress the extraction of river area, geometric features, as well as the construction of cross-sections.Accordingly, accurate mapping of river areas and the construction of boundaries can play a decisive role in studying river morphodynamics, ecological restoration, estimation of hydraulic parameters (e.g., discharge) and management of water resources [7,11,12].
To address the problem of spectral mixture, a wide variety of unmixing and soft classification algorithms, including physics-based and data-driven techniques, are developed [13][14][15][16].These techniques estimate the fraction of each class within the pixels and are represented on a set of grayscale images.Although soft classifiers can reveal sub-pixel information, the spatial distribution of fractions still remains unknown.To tackle this problem, super resolution mapping (SRM) techniques, such as pixel swapping [17], are developed in order to spatially allocate the sub-pixels corresponding to fractions obtained from soft classification.Although some techniques have been applied to retrieve sub-pixel details in rivers, most of the previous research is limited to the estimation of water fractions regardless of sub-pixel mapping (e.g., [11,18]).Concerned with sub-pixel mapping, there have been some attempts to apply SRM techniques for shoreline and waterline detection [19][20][21], flood inundation mapping [22,23], as well as for river mapping [24].
This research aims at the estimation of water fractions within the mixed pixels (i.e., unmixing) and then the spatial allocation of corresponding sub-pixels (i.e., SRM) to map river boundaries at the sub-pixel level.To this end, NDWIs are leveraged for the estimation of water fractions.The original NDWIs are developed based on a specific pair of spectral bands (e.g., NIR and green bands) and with the primary aim of hard delineation of water bodies [25,26].However, recent studies show that different band combinations can be used to form NDWIs, particularly using the imagery with high spectral resolutions [8,[27][28][29].This research proposes a full search approach to identify the optimal pair of bands leading to the highest correlation of NDWI values with water fractions.The selected optimal NDWI is then used as a predictor for the estimation of water fractions through a regression model.The accuracy of the proposed method is compared against an advanced unmixing method, namely fully-constrained simplex projection unmixing (SPU) [15].A thorough investigation is then carried out on the performance of SRM techniques in the context of river mapping.Several SRM techniques are focused including spatial optimization techniques such as pixel swapping (PS) [17], as well as some interpolation-based algorithms.Furthermore, the PS algorithm is modified to speed up the binary water/non-water classification.Both semi-simulated and the fractions derived from real imagery are used for evaluation of SRM techniques.The first of these provides the possibility of accuracy assessment of the sole spatial allocation of sub-pixels task, while the latter also considers the uncertainties involved in estimation of water fractions.In addition, the effectiveness of current thresholding methods on NDWIs are examined for hard water/non-water classification.Small rivers have been the interest of this study and, accordingly, high-resolution satellite imagery (HRSI), including WorldView-3 (WV-3) and GeoEye imagery, are used to exercise the implementations.
The remainder of the paper is organized as follows: the following section provides a brief overview about the sub-pixel reconstruction of river boundaries.The case studies and dataset are presented in Section 2. Section 3 describes the proposed methodology encompassing estimation of water fractions, as well as SRM techniques.The results of the implementations are outlined in Section 4. Section 5 provides discussions about results and, finally, concluding remarks and suggestions for future work are addressed in Section 6.

Study Area and Imagery Set
Sarca and Noce, two Alpine rivers located in Northeast Italy, are considered as the case studies.The areas of interest are about 19 km and 8 km long reaches for Sarca and Noce rivers, respectively (Figure 1).HRSIs including an eight-band WV-3 image of the Sarca river acquired on 1 September 2015 and a GeoEye image of the Noce river acquired on 17 June 2012 are used in this research.These sensors provide sufficient spatial resolutions (about 2 m) for the case studies where the mean channel widths are not more than 20 m.The multispectral bands of WV-3 and GeoEye imagery are described in Table 1.Note that WV-3 also provides eight shortwave-infrared (SWIR) bands which are not acquired in this study.In this research, the effectiveness of additional multispectral bands of WV-3 is investigated compared to a conventional four-band high resolution sensor (i.e., GeoEye).This analysis would be valid for a wide range of high-resolution sensors providing the same or similar spectral and radiometric resolutions (Table 1).
Table 1.Multispectral bands of GeoEye and WV-3 sensors and list of the sensors providing same/similar spectral resolutions.

Sensor
Multispectral Bands

Methodology
As mentioned, the estimation of water fractions is the requirement and input for the SRM of river boundaries.Both real and semi-simulated water fractions are considered in order to isolate and examine the effect of fractions on final sub-pixel maps.The "semi-simulated" fractions are contrived inputs for SRM algorithms with known water fractions.As the real geometry of the river is considered for simulations, the term "semi-simulated" is preferred rather than "simulated".On the other hand, the "real fractions" refer to the fractions obtained from unmixing of real imagery.These two testing scenarios are highlighted in the general workflow of the proposed methodology in Figure 2. The proposed methodology is elaborated in detail in the following sub-sections emphasizing both unmixing and SRM steps.The multispectral bands of WV-3 and GeoEye imagery are described in Table 1.Note that WV-3 also provides eight shortwave-infrared (SWIR) bands which are not acquired in this study.In this research, the effectiveness of additional multispectral bands of WV-3 is investigated compared to a conventional four-band high resolution sensor (i.e., GeoEye).This analysis would be valid for a wide range of high-resolution sensors providing the same or similar spectral and radiometric resolutions (Table 1).

Methodology
As mentioned, the estimation of water fractions is the requirement and input for the SRM of river boundaries.Both real and semi-simulated water fractions are considered in order to isolate and examine the effect of fractions on final sub-pixel maps.The "semi-simulated" fractions are contrived inputs for SRM algorithms with known water fractions.As the real geometry of the river is considered for simulations, the term "semi-simulated" is preferred rather than "simulated".On the other hand, the "real fractions" refer to the fractions obtained from unmixing of real imagery.These two testing scenarios are highlighted in the general workflow of the proposed methodology in Figure 2. The proposed methodology is elaborated in detail in the following sub-sections emphasizing both unmixing and SRM steps.The procedure considered for the mapping and assessment of river boundaries at the sub-pixel resolution using both real and semi-simulated fractions; ZF represents the zoom factor.

Semi-Simulated Fractions
The accuracy of estimated water fractions can inescapably affect the output sub-pixel maps of SRM techniques.In order to evaluate the performance of SRM algorithms only through the spatial arrangement of sub-pixels, it is essential to supply the algorithms with free-of-error inputs (fractions).In support of this objective, semi-simulated fractions are contrived in such a way as to provide known fractions without any uncertainty.
The river area is initially masked out on the high resolution image using a hard classifier to take it as the reference map which has binary water/non-water labels.Then, a moving window with the dimensions equivalent to the user selected zoom factor (ZF) is applied on the reference map.ZF defines the number of sub-pixels that a single pixel can be split in through SRM process.The number of water pixels is counted in each window, which can be simply converted to the fractions in the [0, 1] range.The content of each window (i.e., the number of water pixels) is known from the high-resolution reference map; thereby the fractions can be considered known and with no uncertainty.For instance, semi-simulated fractions are calculated by moving a 5 × 5 window (ZF = 5) on the reference water map in Figure 3.The semi-simulated fractions can provide free-of-error inputs not only for SRM algorithms, but also for hard classifiers.This is because the most abundant class is known within the windows (coarse pixels).Therefore, a water/non-water label can be assigned without the usual uncertainties of hard classifiers concerned with spectral confusions and similarities existing among different land cover classes (Figure 3c).The procedure considered for the mapping and assessment of river boundaries at the sub-pixel resolution using both real and semi-simulated fractions; ZF represents the zoom factor.

Semi-Simulated Fractions
The accuracy of estimated water fractions can inescapably affect the output sub-pixel maps of SRM techniques.In order to evaluate the performance of SRM algorithms only through the spatial arrangement of sub-pixels, it is essential to supply the algorithms with free-of-error inputs (fractions).In support of this objective, semi-simulated fractions are contrived in such a way as to provide known fractions without any uncertainty.
The river area is initially masked out on the high resolution image using a hard classifier to take it as the reference map which has binary water/non-water labels.Then, a moving window with the dimensions equivalent to the user selected zoom factor (ZF) is applied on the reference map.ZF defines the number of sub-pixels that a single pixel can be split in through SRM process.The number of water pixels is counted in each window, which can be simply converted to the fractions in the [0, 1] range.The content of each window (i.e., the number of water pixels) is known from the high-resolution reference map; thereby the fractions can be considered known and with no uncertainty.For instance, semi-simulated fractions are calculated by moving a 5 × 5 window (ZF = 5) on the reference water map in Figure 3.The semi-simulated fractions can provide free-of-error inputs not only for SRM algorithms, but also for hard classifiers.This is because the most abundant class is known within the windows (coarse pixels).Therefore, a water/non-water label can be assigned without the usual uncertainties of hard classifiers concerned with spectral confusions and similarities existing among different land cover classes (Figure 3c).Following the simulation of water fractions, SRM algorithms and conventional hard classification can be applied and compared.A possible random arrangement of sub-pixels for the given example is presented in Figure 3d where the goal of SRM algorithms is to find the proper spatial distribution of sub-pixels.
water fractions, SRM algorithms and conventional hard classification can be applied and compared.A possible random arrangement of sub-pixels for the given example is presented in Figure 3d where the goal of SRM algorithms is to find the proper spatial distribution of sub-pixels.

Real Fractions
Real fractions are used to examine the accuracy of the entire procedure of sub-pixel mapping considering both unmixing and SRM stages.To this end, a new method based on water indices (i.e.OBA-NDWI) is developed for estimation of water fractions.The proposed method is compared with an advanced spectral unmixing method, namely the fully-constrained simplex projection unmixing (SPU) developed by Heylen et al. [15].
According to the described procedure (Figure 2), the river area and a buffer zone along the channel is initially masked by applying an unsupervised clustering followed by a morphological dilation.This is for ease of computations and selection of training samples.However, the unmixing procedure can be applied on the entire image scene.The unmixing methods are then applied on down-sampled imagery (with the same ZF of desired sub-pixel map) so that the hard classification of original high resolution images are considered as reference maps for accuracy assessment of sub-pixel maps (see Figure 2).

•
Proposed Optimal Bands Analysis for NDWI (OBA-NDWI) Making use of water indices has been a long-standing interest of researchers for delineating water features [8,25,26,30].McFeeters [25] developed the original normalized difference water index (NDWI) based on green (G) and near-infrared (NIR) bands of multispectral imagery (Equation (1)).Afterward, Xu [26] presented the modified NDWI (MNDWI) by replacing the NIR band with the SWIR band: The spectral channels of WV-3 provide a wealth of choices for developing different types of spectral indices (e.g., water, vegetation, and soil indices) compared to the conventional HRSI (e.g., GeoEye images).Different combinations of spectral bands can be considered for calculation of NDWI from WV-3 images; for instance, (G, NIR-2) and (coastal-blue (CB), NIR-2) are pairs of bands used in some previous studies [27][28][29].More recently, Xie et al. [8] considered the normalized ratio of a relatively high-reflective band (bh) and a low-reflective band (bl) for which different pairs of bands can be employed for the computation of NDWI according to the following equation: The outputs of NDWIs ranges from −1 to 1, where selection of an appropriate threshold is the main barrier to distinguish among water and non-water pixels [2,17].McFeeters [25] and Xu [26] considered zero as the threshold when they introduced NDWI and MNDWI.Ji et al. [31] have investigated the efficiency of NDWIs by considering different combination of spectral bands on coarse and moderate-resolution imagery (e.g., MODIS and Landsat ETM+).They have determined

Real Fractions
Real fractions are used to examine the accuracy of the entire procedure of sub-pixel mapping considering both unmixing and SRM stages.To this end, a new method based on water indices (i.e., OBA-NDWI) is developed for estimation of water fractions.The proposed method is compared with an advanced spectral unmixing method, namely the fully-constrained simplex projection unmixing (SPU) developed by Heylen et al. [15].
According to the described procedure (Figure 2), the river area and a buffer zone along the channel is initially masked by applying an unsupervised clustering followed by a morphological dilation.This is for ease of computations and selection of training samples.However, the unmixing procedure can be applied on the entire image scene.The unmixing methods are then applied on down-sampled imagery (with the same ZF of desired sub-pixel map) so that the hard classification of original high resolution images are considered as reference maps for accuracy assessment of sub-pixel maps (see Figure 2).
• Proposed Optimal Bands Analysis for NDWI (OBA-NDWI) Making use of water indices has been a long-standing interest of researchers for delineating water features [8,25,26,30].McFeeters [25] developed the original normalized difference water index (NDWI) based on green (G) and near-infrared (NIR) bands of multispectral imagery (Equation (1)).Afterward, Xu [26] presented the modified NDWI (MNDWI) by replacing the NIR band with the SWIR band: The spectral channels of WV-3 provide a wealth of choices for developing different types of spectral indices (e.g., water, vegetation, and soil indices) compared to the conventional HRSI (e.g., GeoEye images).Different combinations of spectral bands can be considered for calculation of NDWI from WV-3 images; for instance, (G, NIR-2) and (coastal-blue (CB), NIR-2) are pairs of bands used in some previous studies [27][28][29].More recently, Xie et al. [8] considered the normalized ratio of a relatively high-reflective band (b h ) and a low-reflective band (b l ) for which different pairs of bands can be employed for the computation of NDWI according to the following equation: The outputs of NDWIs ranges from −1 to 1, where selection of an appropriate threshold is the main barrier to distinguish among water and non-water pixels [2,17].McFeeters [25] and Xu [26] considered zero as the threshold when they introduced NDWI and MNDWI.Ji et al. [31] have investigated the efficiency of NDWIs by considering different combination of spectral bands on coarse and moderate-resolution imagery (e.g., MODIS and Landsat ETM+).They have determined threshold values based on synthetic mixture of the three dominant land-cover types (water, vegetation, and soil).Moreover, Otsu's thresholding method [32], which is based on the maximum between-class variance criterion, is used in several articles to mask out water pixels [30,33,34].However, the main aim of previous research is mainly binary (hard) classification of water features, whereas the estimation of water fractions received less attention.
This research introduces optimal band analysis for NDWI (OBA-NDWI) for the estimation of water fractions using NDWIs.First, image-derived spectra for the aforementioned major land-cover types (i.e., water, vegetation, and bare earth/soil) are mixed linearly with all the possible fractions of the classes with one percent intervals (Figure 4).This approach has been adapted from Ji et al. [31]; however, in the case of presence of an outstandingly different major class (endmember), the number of classes can be increased.Note that the number of classes cannot exceed the number of spectral bands to solve the common unmixing methods, such as linear spectral unmixing and SPU [9].This issue can limit the number of unmixing classes using imagery with low spectral resolution.However, the mentioned three major classes which can be derived from grouping the sub-classes with similar spectral characteristics are representative for a wide range of applications and particularly for the riparian zones.This assumption on the number and type of major classes (endmembers) is well-recognized in the previous research [31,35,36].
threshold values based on synthetic mixture of the three dominant land-cover types (water, vegetation, and soil).Moreover, Otsu's thresholding method [32], which is based on the maximum between-class variance criterion, is used in several articles to mask out water pixels [30,33,34].However, the main aim of previous research is mainly binary (hard) classification of water features, whereas the estimation of water fractions received less attention.
This research introduces optimal band analysis for NDWI (OBA-NDWI) for the estimation of water fractions using NDWIs.First, image-derived spectra for the aforementioned major land-cover types (i.e., water, vegetation, and bare earth/soil) are mixed linearly with all the possible fractions of the classes with one percent intervals (Figure 4).This approach has been adapted from Ji et al. [31]; however, in the case of presence of an outstandingly different major class (endmember), the number of classes can be increased.Note that the number of classes cannot exceed the number of spectral bands to solve the common unmixing methods, such as linear spectral unmixing and SPU [9].This issue can limit the number of unmixing classes using imagery with low spectral resolution.However, the mentioned three major classes which can be derived from grouping the sub-classes with similar spectral characteristics are representative for a wide range of applications and particularly for the riparian zones.This assumption on the number and type of major classes (endmembers) is well-recognized in the previous research [31,35,36].After producing the synthetically mixed spectra with all the possible fractions of the three classes, the NDWI values of the spectra are regressed against known water fractions.Then, OBA-NDWI identifies the pair of bands for which the corresponding NDWI values yield the highest correlation with the water fractions.All the possible combinations of spectral bands are considered to calculate the NDWI values (Equation (3)) in order to examine their potential for estimation of water fractions by assessing the coefficient of determination and RMSE of the regression models.
In Equation (3), the pair of bands used for estimation of NDWI is denoted by bi and bj for an n-band image that subscripts i and j can be assigned to different bands of an image considering the specified ranges in the equation.This provides a full search among all the possible options for choosing a pair of spectral bands which is also common in other applications such as bathymetry from optical imagery [37].The total number of possible pairs of spectral bands would be for an n-band image, which means 28 pairs for a WV-3 image.After identification of the optimal combination of bands to be used for calculation of NDWI, the corresponding regression model is used to predict the water fractions of the image pixels.The regression-based approaches are previously used for estimation of fractional vegetation coverage using vegetation indices [38][39][40] but have not been explored yet for the estimation of water fractions based on water indices.Figure 5 illustrates the step-by-step procedure for estimation of water fractions based on the proposed OBA-NDWI.After producing the synthetically mixed spectra with all the possible fractions of the three classes, the NDWI values of the spectra are regressed against known water fractions.Then, OBA-NDWI identifies the pair of bands for which the corresponding NDWI values yield the highest correlation with the water fractions.All the possible combinations of spectral bands are considered to calculate the NDWI values (Equation ( 3)) in order to examine their potential for estimation of water fractions by assessing the coefficient of determination and RMSE of the regression models.

NDW I_proposed
In Equation ( 3), the pair of bands used for estimation of NDWI is denoted by b i and b j for an n-band image that subscripts i and j can be assigned to different bands of an image considering the specified ranges in the equation.This provides a full search among all the possible options for choosing a pair of spectral bands which is also common in other applications such as bathymetry from optical imagery [37].The total number of possible pairs of spectral bands would be n !/((n − 2)! × 2 !) for an n-band image, which means 28 pairs for a WV-3 image.After identification of the optimal combination of bands to be used for calculation of NDWI, the corresponding regression model is used to predict the water fractions of the image pixels.The regression-based approaches are previously used for estimation of fractional vegetation coverage using vegetation indices [38][39][40] but have not been explored yet for the estimation of water fractions based on water indices.Figure 5 illustrates the step-by-step procedure for estimation of water fractions based on the proposed OBA-NDWI.Moreover, the effectiveness of common thresholding methods is examined for the hard classification of river area using HRSI.In this regard, the sensitivity of thresholding methods on different NDWIs is investigated and, accordingly, the relevant cut-off water fraction for each NDWI is determined.The determination of the minimum water fraction within the masked pixels can reveal the performance of hard classification built upon each of thresholding methods.This is because, considering the concept of hard classification for binary mapping, at least half of a pixel should be occupied by the desired class (water) to be assigned to that class.

• Simplex Projection Unmixing (SPU)
The recent and advanced algorithm of SPU is also applied on the imagery to evaluate the effectiveness of the proposed OBA-NDWI method for estimation of water fractions.The SPU is a technique ensuring sum-to-unity, as well as positivity, of fractions.These constraints are prerequisites for most SRM techniques, such as the PS algorithm.The SPU method is built upon the geometrical interpretation of the unmixing problem [41].This algorithm uses a sequence of orthogonal projections on sub-simplexes formed by the endmembers.The solution derived from the SPU algorithm minimizes the least squares error while firmly respecting the sum-to-unity and non-negativity constraints on the fractions (see [42] for further mathematical details).

Super Resolution Mapping (SRM)
SRM is required for spatial allocation of water fractions in order to reconstruct the river area at the sub-pixel level.Most of the SRM algorithms rely on maximizing the spatial proximity of sub-pixels with the same labels in their neighborhood [43,44].Several techniques are developed based on artificial intelligence optimization techniques, such as genetic algorithms [45,46] and neural networks [47,48], and also other methods based on Markov random fields [49].In addition to these methods, the PS algorithm [17] is a widely-used and efficient method for SRM purposes [50].On the other hand, interpolation-based approaches could be considered as alternatives for SRM.These algorithms are non-iterative, which can be an advantage when the computation time is considered.Although some interpolation methods (e.g., IDW and kriging) are tested for SRM using simulated imagery [51], a comparative analysis with respect to iterative techniques remains poorly investigated, particularly in terms of mapping water bodies.This research concentrates on the PS algorithm as an iterative technique and a couple of interpolation methods as non-iterative solutions for sub-pixel mapping of river boundaries.In addition, a non-iterative approach is proposed for improving the speed of PS algorithm.Moreover, the effectiveness of common thresholding methods is examined for the hard classification of river area using HRSI.In this regard, the sensitivity of thresholding methods on different NDWIs is investigated and, accordingly, the relevant cut-off water fraction for each NDWI is determined.The determination of the minimum water fraction within the masked pixels can reveal the performance of hard classification built upon each of thresholding methods.This is because, considering the concept of hard classification for binary mapping, at least half of a pixel should be occupied by the desired class (water) to be assigned to that class.

• Simplex Projection Unmixing (SPU)
The recent and advanced algorithm of SPU is also applied on the imagery to evaluate the effectiveness of the proposed OBA-NDWI method for estimation of water fractions.The SPU is a technique ensuring sum-to-unity, as well as positivity, of fractions.These constraints are prerequisites for most SRM techniques, such as the PS algorithm.The SPU method is built upon the geometrical interpretation of the unmixing problem [41].This algorithm uses a sequence of orthogonal projections on sub-simplexes formed by the endmembers.The solution derived from the SPU algorithm minimizes the least squares error while firmly respecting the sum-to-unity and non-negativity constraints on the fractions (see [42] for further mathematical details).

Super Resolution Mapping (SRM)
SRM is required for spatial allocation of water fractions in order to reconstruct the river area at the sub-pixel level.Most of the SRM algorithms rely on maximizing the spatial proximity of sub-pixels with the same labels in their neighborhood [43,44].Several techniques are developed based on artificial intelligence optimization techniques, such as genetic algorithms [45,46] and neural networks [47,48], and also other methods based on Markov random fields [49].In addition to these methods, the PS algorithm [17] is a widely-used and efficient method for SRM purposes [50].On the other hand, interpolation-based approaches could be considered as alternatives for SRM.These algorithms are non-iterative, which can be an advantage when the computation time is considered.Although some interpolation methods (e.g., IDW and kriging) are tested for SRM using simulated imagery [51], a comparative analysis with respect to iterative techniques remains poorly investigated, particularly in terms of mapping water bodies.This research concentrates on the PS algorithm as an iterative technique and a couple of interpolation methods as non-iterative solutions for sub-pixel mapping of river boundaries.In addition, a non-iterative approach is proposed for improving the speed of PS algorithm.

Pixel Swapping (PS)
The number of sub-pixels corresponding to the fraction of a given class within each pixel can be calculated based on the desired ZF according to the following equation: where N k and F k stand, respectively, for the number of sub-pixels and fraction of the class k within a pixel.For the binary water/non-water classification, the number of sub-pixels is calculated for the water class based on Equation ( 4), and then its subtraction from ZF 2 gives the number of sub-pixels for the non-water class.
The PS algorithm [17,52] allocates sub-pixel labels in random positions within each pixel.Then, the attractiveness of each sub-pixel with respect to a particular class is predicted as a distance-weighted function of its neighbors according to the following equation: where A k denotes the attractiveness of a sub-pixel with respect to the class k, F k i is the fraction of class k in the i-th neighbor pixel.The number of neighbor pixels is n, and d i is the distance of the i-th neighbor pixel from the sub-pixel for which the attractiveness is computed.
Considering the attractiveness values of the class k in a pixel, the least attractive sub-pixel location initially allocated to the desired class (e.g., water) should be identified, as well as the most attractive location initially allocated to the other class (e.g., non-water).If the attractiveness of the least attractive location is less than that of the most attractive location, then the classes are swapped, otherwise no change is made [43,52].After applying the swapping process on all image pixels, this process should be iterated until reaching a point that the algorithm is not able to perform anymore swaps.An example of the swapping process is illustrated in Figure 6 where the central pixel with 58% water fraction is divided to nine sub-pixels (ZF = 3) and the attractiveness of each sub-pixel is calculated based on the eight neighboring pixels.The number of sub-pixels corresponding to the fraction of a given class within each pixel can be calculated based on the desired ZF according to the following equation: where k N and k F stand, respectively, for the number of sub-pixels and fraction of the class k within a pixel.For the binary water/non-water classification, the number of sub-pixels is calculated for the water class based on Equation (4), and then its subtraction from ZF 2 gives the number of sub-pixels for the non-water class.
The PS algorithm [17,52] allocates sub-pixel labels in random positions within each pixel.Then, the attractiveness of each sub-pixel with respect to a particular class is predicted as a distance-weighted function of its neighbors according to the following equation: where A k denotes the attractiveness of a sub-pixel with respect to the class k, F is the fraction of class k in the i-th neighbor pixel.The number of neighbor pixels is n, and di is the distance of the i-th neighbor pixel from the sub-pixel for which the attractiveness is computed.
Considering the attractiveness values of the class k in a pixel, the least attractive sub-pixel location initially allocated to the desired class (e.g., water) should be identified, as well as the most attractive location initially allocated to the other class (e.g., non-water).If the attractiveness of the least attractive location is less than that of the most attractive location, then the classes are swapped, otherwise no change is made [43,52].After applying the swapping process on all image pixels, this process should be iterated until reaching a point that the algorithm is not able to perform anymore swaps.An example of the swapping process is illustrated in Figure 6 where the central pixel with 58% water fraction is divided to nine sub-pixels (ZF = 3) and the attractiveness of each sub-pixel is calculated based on the eight neighboring pixels.

Proposed Modified Binary Pixel Swapping (MBPS)
After calculation of the attractiveness values for sub-pixel locations, the original PS algorithm commences with a random allocation of sub-pixels and then maximizes the spatial dependency in an iterative manner.This can be a barrier in terms of computational time when applying the algorithm on large extents.To tackle this problem, a simple non-iterative solution is examined in this paper for the spatial allocation of binary classes (e.g., water and non-water).The proposed method suggests to simply allocate the sub-pixels of the desired class (water) in the N k locations with the highest attractiveness values toward that class.The remaining sub-pixel locations are then directly assigned to the other class (non-water).The proposed modified binary PS (MBPS) method is applied on the same example discussed in Figure 6.In this case, the resultant sub-pixel map is the same as that of the original PS algorithm (Figure 7).

Proposed Modified Binary Pixel Swapping (MBPS)
After calculation of the attractiveness values for sub-pixel locations, the original PS algorithm commences with a random allocation of sub-pixels and then maximizes the spatial dependency in an iterative manner.This can be a barrier in terms of computational time when applying the algorithm on large extents.To tackle this problem, a simple non-iterative solution is examined in this paper for the spatial allocation of binary classes (e.g., water and non-water).The proposed method suggests to simply allocate the sub-pixels of the desired class (water) in the N k locations with the highest attractiveness values toward that class.The remaining sub-pixel locations are then directly assigned to the other class (non-water).The proposed modified binary PS (MBPS) method is applied on the same example discussed in Figure 6.In this case, the resultant sub-pixel map is the same as that of the original PS algorithm (Figure 7).

Interpolation-Based SRM
The interpolation methods can be employed to develop non-iterative SRM.Assigning the fractions of classes to the center of the corresponding pixels, sub-pixel fractions can be estimated using common interpolation methods.In this research, several techniques including bilinear, bicubic, and lanczos3 are tested; an example using bilinear method is shown in Figure 8.A hard labeling process can be then applied on interpolated points (i.e., sub-pixels) in such a way to assign each sub-pixel to the class with the highest estimated fraction.The main advantage of this technique is that the labeling procedure is done in a single step which improves the computational efficiency.

Estimation of Water Fractions
Semi-simulated fractions are produced along the two rivers based on the methodology described in Section 3.1.As discussed in detail, this provides an input without any error for SRM methods, as well as for the hard classification.The sub-pixel detail neglected by hard classifiers and the performance of SRM methods through allocation of sub-pixels are investigated thoroughly (Section 4.2).
Toward estimation of real water fractions using the proposed method relying on OBA-NDWI, training samples for the three desired classes (i.e., water, vegetation, and soil) are selected from the imagery.The spectra of the endmembers are estimated based on the average values of the samples which are selected from different spots distributed along the rivers.For each endmember, about 30 pixels are selected by visual interpretation.Then the linear mixture is performed to synthetize all the possible fractions concerned with the three endmembers.NDWIs with all the possible combinations of spectral bands are applied on the synthetically mixed WV-3 and GeoEye spectra; the outputs are plotted against the known values of water fractions (some examples are shown in Figures 9 and 10).The regression model providing the highest correlation is then used to predict the water fractions of all image pixels.In addition, zero-thresholds and Otsu's thresholds are applied on NDWI values to assess these two conventional threshold methods for hard delineation of water bodies.

Interpolation-Based SRM
The interpolation methods can be employed to develop non-iterative SRM.Assigning the fractions of classes to the center of the corresponding pixels, sub-pixel fractions can be estimated using common interpolation methods.In this research, several techniques including bilinear, bicubic, and lanczos3 are tested; an example using bilinear method is shown in Figure 8.A hard labeling process can be then applied on interpolated points (i.e., sub-pixels) in such a way to assign each sub-pixel to the class with the highest estimated fraction.The main advantage of this technique is that the labeling procedure is done in a single step which improves the computational efficiency.

Interpolation-Based SRM
The interpolation methods can be employed to develop non-iterative SRM.Assigning the fractions of classes to the center of the corresponding pixels, sub-pixel fractions can be estimated using common interpolation methods.In this research, several techniques including bilinear, bicubic, and lanczos3 are tested; an example using bilinear method is shown in Figure 8.A hard labeling process can be then applied on interpolated points (i.e., sub-pixels) in such a way to assign each sub-pixel to the class with the highest estimated fraction.The main advantage of this technique is that the labeling procedure is done in a single step which improves the computational efficiency.

Estimation of Water Fractions
Semi-simulated fractions are produced along the two rivers based on the methodology described in Section 3.1.As discussed in detail, this provides an input without any error for SRM methods, as well as for the hard classification.The sub-pixel detail neglected by hard classifiers and the performance of SRM methods through allocation of sub-pixels are investigated thoroughly (Section 4.2).
Toward estimation of real water fractions using the proposed method relying on OBA-NDWI, training samples for the three desired classes (i.e., water, vegetation, and soil) are selected from the imagery.The spectra of the endmembers are estimated based on the average values of the samples which are selected from different spots distributed along the rivers.For each endmember, about 30 pixels are selected by visual interpretation.Then the linear mixture is performed to synthetize all the possible fractions concerned with the three endmembers.NDWIs with all the possible combinations of spectral bands are applied on the synthetically mixed WV-3 and GeoEye spectra; the outputs are plotted against the known values of water fractions (some examples are shown in Figures 9 and 10).The regression model providing the highest correlation is then used to predict the water fractions of all image pixels.In addition, zero-thresholds and Otsu's thresholds are applied on NDWI values to assess these two conventional threshold methods for hard delineation of water bodies.

Estimation of Water Fractions
Semi-simulated fractions are produced along the two rivers based on the methodology described in Section 3.1.As discussed in detail, this provides an input without any error for SRM methods, as well as for the hard classification.The sub-pixel detail neglected by hard classifiers and the performance of SRM methods through allocation of sub-pixels are investigated thoroughly (Section 4.2).
Toward estimation of real water fractions using the proposed method relying on OBA-NDWI, training samples for the three desired classes (i.e., water, vegetation, and soil) are selected from the imagery.The spectra of the endmembers are estimated based on the average values of the samples which are selected from different spots distributed along the rivers.For each endmember, about 30 pixels are selected by visual interpretation.Then the linear mixture is performed to synthetize all the possible fractions concerned with the three endmembers.NDWIs with all the possible combinations of spectral bands are applied on the synthetically mixed WV-3 and GeoEye spectra; the outputs are plotted against the known values of water fractions (some examples are shown in Figures 9 and 10).The regression model providing the highest correlation is then used to predict the water fractions of all image pixels.In addition, zero-thresholds and Otsu's thresholds are applied on NDWI values to assess these two conventional threshold methods for hard delineation of water bodies.As evident in Figures 9 and 10, the water fractions corresponding to zero-threshold and also Otsu's threshold are variable for different NDWIs and their difference can be outstanding.The minimum water fraction associated with zero-threshold on WV-3 image is significantly variable, from 87% to 35% for the NDWIs constructed, respectively, by two band combinations of (CB, NIR1) and (B, RE).This variation ranges from 91% to 67% for the Otsu's threshold with the same NDWIs.These fluctuations also exist for the GeoEye image for which the minimum water fraction for zero-threshold ranges from 57% to 75%, and for the Otsu's threshold from 75% to 85%, respectively, using NDWIs based on (B, NIR) and (G, NIR) pairs of bands (Figure 11).Therefore, the range of water fractions of masked out pixels based on common thresholding methods can be variable depending on the spectral bands used in the calculation of NDWI.With generalization of all the classes excluding the water class to a single non-water class, each pixel with a minimum 50% water content can be labeled as water.This is a binary hard classification scheme which can be considered for delineation of water bodies based on thresholding on NDWI values.Keeping this in mind, Figure 11 reveals that Otsu's thresholding most probably occurs above the 50% water fraction which leads  As evident in Figures 9 and 10, the water fractions corresponding to zero-threshold and also Otsu's threshold are variable for different NDWIs and their difference can be outstanding.The minimum water fraction associated with zero-threshold on WV-3 image is significantly variable, from 87% to 35% for the NDWIs constructed, respectively, by two band combinations of (CB, NIR1) and (B, RE).This variation ranges from 91% to 67% for the Otsu's threshold with the same NDWIs.These fluctuations also exist for the GeoEye image for which the minimum water fraction for zero-threshold ranges from 57% to 75%, and for the Otsu's threshold from 75% to 85%, respectively, using NDWIs based on (B, NIR) and (G, NIR) pairs of bands (Figure 11).Therefore, the range of water fractions of masked out pixels based on common thresholding methods can be variable depending on the spectral bands used in the calculation of NDWI.With generalization of all the classes excluding the water class to a single non-water class, each pixel with a minimum 50% water content can be labeled as water.This is a binary hard classification scheme which can be considered for delineation of water bodies based on thresholding on NDWI values.Keeping this in mind, Figure As evident in Figures 9 and 10, the water fractions corresponding to zero-threshold and also Otsu's threshold are variable for different NDWIs and their difference can be outstanding.The minimum water fraction associated with zero-threshold on WV-3 image is significantly variable, from 87% to 35% for the NDWIs constructed, respectively, by two band combinations of (CB, NIR1) and (B, RE).This variation ranges from 91% to 67% for the Otsu's threshold with the same NDWIs.These fluctuations also exist for the GeoEye image for which the minimum water fraction for zero-threshold ranges from 57% to 75%, and for the Otsu's threshold from 75% to 85%, respectively, using NDWIs based on (B, NIR) and (G, NIR) pairs of bands (Figure 11).Therefore, the range of water fractions of masked out pixels based on common thresholding methods can be variable depending on the spectral bands used in the calculation of NDWI.With generalization of all the classes excluding the water class to a single non-water class, each pixel with a minimum 50% water content can be labeled as water.This is a binary hard classification scheme which can be considered for delineation of water bodies based on thresholding on NDWI values.Keeping this in mind, Figure 11 reveals that Otsu's thresholding most probably occurs above the 50% water fraction which leads to underestimation of water pixels.For instance, the pixels with fractions ranging from 50% to 85% will be mislabeled as non-water pixels for the NDWI based on (B, NIR1) bands.On the other hand, zero-thresholding deals with both underestimation and overestimation of water pixels.For instance, every pixel with the minimum 35% water content will be labeled as water by applying the zero-threshold on NDWI obtained from (B, RE) bands.This kind of mislabeling at the pixel level can increase the uncertainty of the water mask in mixed boundary pixels.However, extraction of the water pixels in the main river channel would be less challenging.This is because they are pure or near pure pixels (tending to 100% water fraction) so that their fractions are much higher than the minimum water fraction detectable by thresholding methods using different band combinations (Figure 11).zero-thresholding deals with both underestimation and overestimation of water pixels.For instance, every pixel with the minimum 35% water content will be labeled as water by applying the zero-threshold on NDWI obtained from (B, RE) bands.This kind of mislabeling at the pixel level can increase the uncertainty of the water mask in mixed boundary pixels.However, extraction of the water pixels in the main river channel would be less challenging.This is because they are pure or near pure pixels (tending to 100% water fraction) so that their fractions are much higher than the minimum water fraction detectable by thresholding methods using different band combinations (Figure 11).The proposed OBA-NDWI method is applied on the synthetically mixed spectra in order to exploit the strongest relation between water fractions and NDWI values (the regression lines are illustrated on Figures 9 and 10).The OBA-NDWI performs second order regression of NDWI values versus water fractions for all possible combinations of spectral bands in order to identify the optimal pair of bands which their corresponding NDWI yields the highest coefficient of determination (R 2 ).Consideration of all the possible pair of bands is a systematic approach to identify the optimal structure of the NDWI that the results are in line with the assumption of Xie et al. [8] regarding the use of a relatively high-reflective band and a low-reflective band for calculation of the NDWI.As zero-thresholding deals with both underestimation and overestimation of water pixels.For instance, every pixel with the minimum 35% water content will be labeled as water by applying the zero-threshold on NDWI obtained from (B, RE) bands.This kind of mislabeling at the pixel level can increase the uncertainty of the water mask in mixed boundary pixels.However, extraction of the water pixels in the main river channel would be less challenging.This is because they are pure or near pure pixels (tending to 100% water fraction) so that their fractions are much higher than the minimum water fraction detectable by thresholding methods using different band combinations (Figure 11).The proposed OBA-NDWI method is applied on the synthetically mixed spectra in order to exploit the strongest relation between water fractions and NDWI values (the regression lines are illustrated on Figures 9 and 10).The OBA-NDWI performs second order regression of NDWI values versus water fractions for all possible combinations of spectral bands in order to identify the optimal pair of bands which their corresponding NDWI yields the highest coefficient of determination (R 2 ).Consideration of all the possible pair of bands is a systematic approach to identify the optimal structure of the NDWI that the results are in line with the assumption of Xie et al. [8] regarding the use of a relatively high-reflective band and a low-reflective band for calculation of the NDWI.As The proposed OBA-NDWI method is applied on the synthetically mixed spectra in order to exploit the strongest relation between water fractions and NDWI values (the regression lines are illustrated on Figures 9 and 10).The OBA-NDWI performs second order regression of NDWI values versus water fractions for all possible combinations of spectral bands in order to identify the optimal pair of bands which their corresponding NDWI yields the highest coefficient of determination (R 2 ).Consideration of all the possible pair of bands is a systematic approach to identify the optimal structure of the NDWI that the results are in line with the assumption of Xie et al. [8] regarding the use of a relatively high-reflective band and a low-reflective band for calculation of the NDWI.As illustrated in Figure 12, coastal-blue (CB), blue (B), and green (G) bands can be considered as optimal high-reflective bands (b h ), while the portion of the spectrum covering red-edge (RE), NIR1, and NIR2 could be effective as low-reflective bands (b l ) using the WV-3 image.In particular, CB and RE bands provide the strongest relation with an R 2 on the order of 0.97 and an RMSE of 2%.This is while the original NDWI with the (G, NIR1) band combination provides significantly less accuracy (R 2 of 0.85 and RMSE of 12%).Blue and NIR are the optimal pair of bands for the GeoEye image.Their corresponding NDWI yields an R 2 value of 0.92 and an RMSE value of 7% through a quadratic relation with water fractions (Figure 13).
ISPRS Int.J. Geo-Inf.2017, 6, 383 12 of 21 illustrated in Figure 12, coastal-blue (CB), blue (B), and green (G) bands can be considered as optimal high-reflective bands (bh), while the portion of the spectrum covering red-edge (RE), NIR1, and NIR2 could be effective as low-reflective bands (bl) using the WV-3 image.In particular, CB and RE bands provide the strongest relation with an R 2 on the order of 0.97 and an RMSE of 2%.This is while the original NDWI with the (G, NIR1) band combination provides significantly less accuracy (R 2 of 0.85 and RMSE of 12%).Blue and NIR are the optimal pair of bands for the GeoEye image.Their corresponding NDWI yields an R 2 value of 0.92 and an RMSE value of 7% through a quadratic relation with water fractions (Figure 13).The quadratic model obtained for the optimal pair of bands is used to predict the water fractions of all image pixels.Moreover, the same endmembers used for linear mixtures are introduced to the SPU algorithm and the resultant fractions are compared with those obtained from the proposed method (Figure 14).In general, a strong correlation (R 2 = 0.9, RMSE = 7% for the WV-3 image and R 2 = 0.87, RMSE = 9% for the GeoEye image) is observed between the estimated fractions of the proposed OBA-NDWI method and the SPU algorithm.

Super Resolution Mapping (SRM)
The semi-simulated, as well as the real water fractions are used as inputs for SRM algorithms.The super resolved maps obtained from the semi-simulated fractions are represented for a river segment in Figure 15.The sub-pixel mapping procedure of interpolation-based techniques is also illustrated in Figure 16 where the interpolated water fractions are represented (Figure 16d) along The quadratic model obtained for the optimal pair of bands is used to predict the water fractions of all image pixels.Moreover, the same endmembers used for linear mixtures are introduced to the SPU algorithm and the resultant fractions are compared with those obtained from the proposed method (Figure 14).In general, a strong correlation (R 2 = 0.9, RMSE = 7% for the WV-3 image and R 2 = 0.87, RMSE = 9% for the GeoEye image) is observed between the estimated fractions of the proposed OBA-NDWI method and the SPU algorithm.
ISPRS Int.J. Geo-Inf.2017, 6, 383 12 of 21 illustrated in Figure 12, coastal-blue (CB), blue (B), and green (G) bands can be considered as optimal high-reflective bands (bh), while the portion of the spectrum covering red-edge (RE), NIR1, and NIR2 could be effective as low-reflective bands (bl) using the WV-3 image.In particular, CB and RE bands provide the strongest relation with an R 2 on the order of 0.97 and an RMSE of 2%.This is while the original NDWI with the (G, NIR1) band combination provides significantly less accuracy (R 2 of 0.85 and RMSE of 12%).Blue and NIR are the optimal pair of bands for the GeoEye image.Their corresponding NDWI yields an R 2 value of 0.92 and an RMSE value of 7% through a quadratic relation with water fractions (Figure 13).The quadratic model obtained for the optimal pair of bands is used to predict the water fractions of all image pixels.Moreover, the same endmembers used for linear mixtures are introduced to the SPU algorithm and the resultant fractions are compared with those obtained from the proposed method (Figure 14).In general, a strong correlation (R 2 = 0.9, RMSE = 7% for the WV-3 image and R 2 = 0.87, RMSE = 9% for the GeoEye image) is observed between the estimated fractions of the proposed OBA-NDWI method and the SPU algorithm.

Super Resolution Mapping (SRM)
The semi-simulated, as well as the real water fractions are used as inputs for SRM algorithms.The super resolved maps obtained from the semi-simulated fractions are represented for a river segment in Figure 15.The sub-pixel mapping procedure of interpolation-based techniques is also

Super Resolution Mapping (SRM)
The semi-simulated, as well as the real water fractions are used as inputs for SRM algorithms.The super resolved maps obtained from the semi-simulated fractions are represented for a river segment in Figure 15.The sub-pixel mapping procedure of interpolation-based techniques is also illustrated in Figure 16 where the interpolated water fractions are represented (Figure 16d) along with the hard labels assigned to each sub-pixel (Figure 16e).As it is clear from the illustrations, hard classification is very rough on the river boundaries, while the SRM techniques reconstruct boundaries with sub-pixel details.User and producer accuracies of the river area extracted from the hard classification and SRM of semi-simulated fractions are presented in Figures 17 and 18 for the Sarca and Noce rivers, respectively.User and producer accuracies of the river area extracted from the hard classification and SRM of semi-simulated fractions are presented in Figures 17 and 18 for the Sarca and Noce rivers, respectively.The user/producer accuracies of hard classification are remarkably lower than that of sub-pixel maps.For instance, the error maps of hard classification and a sub-pixel map of the same area are represented in Figure 19, which shows that the hard classified map encompasses a large number of misclassified pixels.This issue gets worse with higher ZF where the difference between the accuracies of hard classification and SRM methods reaches above 10%.Interpolation-based SRM and MBPS lead to comparable accuracies with respect to the PS algorithm except at higher zoom factors (ZF = 10) where the producer accuracies of interpolation-based techniques are lower.The real water fractions obtained from the SPU algorithm, as well as those obtained from the proposed OBA-NDWI, are used as inputs for SRM algorithms.This provides the possibility of investigating the effect of uncertainties associated with water fractions on the final sub-pixel maps.The user/producer accuracies of hard classification are remarkably lower than that of sub-pixel maps.For instance, the error maps of hard classification and a sub-pixel map of the same area are represented in Figure 19, which shows that the hard classified map encompasses a large number of misclassified pixels.This issue gets worse with higher ZF where the difference between the accuracies of hard classification and SRM methods reaches above 10%.Interpolation-based SRM and MBPS lead to comparable accuracies with respect to the PS algorithm except at higher zoom factors (ZF = 10) where the producer accuracies of interpolation-based techniques are lower.The real water fractions obtained from the SPU algorithm, as well as those obtained from the proposed OBA-NDWI, are used as inputs for SRM algorithms.This provides the possibility of investigating the effect of uncertainties associated with water fractions on the final sub-pixel maps.The user/producer accuracies of hard classification are remarkably lower than that of sub-pixel maps.For instance, the error maps of hard classification and a sub-pixel map of the same area are represented in Figure 19, which shows that the hard classified map encompasses a large number of misclassified pixels.This issue gets worse with higher ZF where the difference between the accuracies of hard classification and SRM methods reaches above 10%.Interpolation-based SRM and MBPS lead to comparable accuracies with respect to the PS algorithm except at higher zoom factors (ZF = 10) where the producer accuracies of interpolation-based techniques are lower.The user/producer accuracies of hard classification are remarkably lower than that of sub-pixel maps.For instance, the error maps of hard classification and a sub-pixel map of the same area are represented in Figure 19, which shows that the hard classified map encompasses a large number of misclassified pixels.This issue gets worse with higher ZF where the difference between the accuracies of hard classification and SRM methods reaches above 10%.Interpolation-based SRM and MBPS lead to comparable accuracies with respect to the PS algorithm except at higher zoom factors (ZF = 10) where the producer accuracies of interpolation-based techniques are lower.The real water fractions obtained from the SPU algorithm, as well as those obtained from the proposed OBA-NDWI, are used as inputs for SRM algorithms.This provides the possibility of investigating the effect of uncertainties associated with water fractions on the final sub-pixel maps.The real water fractions obtained from the SPU algorithm, as well as those obtained from the proposed OBA-NDWI, are used as inputs for SRM algorithms.This provides the possibility of investigating the effect of uncertainties associated with water fractions on the final sub-pixel maps.Figure 20 illustrates the sub-pixel maps based on water fractions derived from the OBA-NDWI method.In Figure 21, outputs of the PS/MBPS algorithms are selected from parts of maps where some isolated groups of sub-pixels exist inside the river channel, which is due to errors of estimated water fractions.However, applying a majority filter on the resultant sub-pixel maps can reduce this drawback in order to produce homogenous maps (Figure 21e).These gaps, caused by the underestimation of water fractions by unmixing methods, either the OBA-NDWI or the SPU, are not obviously apparent on the maps of interpolation-based SRM methods due to interpolated values of fractions.Note that the maps are selected randomly from spatially distributed regions along the river channel.Figure 20 illustrates the sub-pixel maps based on water fractions derived from the OBA-NDWI method.In Figure 21, outputs of the PS/MBPS algorithms are selected from parts of maps where some isolated groups of sub-pixels exist inside the river channel, which is due to errors of estimated water fractions.However, applying a majority filter on the resultant sub-pixel maps can reduce this drawback in order to produce homogenous maps (Figure 21e).These gaps, caused by the underestimation of water fractions by unmixing methods, either the OBA-NDWI or the SPU, are not obviously apparent on the maps of interpolation-based SRM methods due to interpolated values of fractions.Note that the maps are selected randomly from spatially distributed regions along the river Figure 20 illustrates the sub-pixel maps based on water fractions derived from the OBA-NDWI method.In Figure 21, outputs of the PS/MBPS algorithms are selected from parts of maps where some isolated groups of sub-pixels exist inside the river channel, which is due to errors of estimated water fractions.However, applying a majority filter on the resultant sub-pixel maps can reduce this drawback in order to produce homogenous maps (Figure 21e).These gaps, caused by the underestimation of water fractions by unmixing methods, either the OBA-NDWI or the SPU, are not obviously apparent on the maps of interpolation-based SRM methods due to interpolated values of fractions.Note that the maps are selected randomly from spatially distributed regions along the river channel.User and producer accuracies of the sub-pixel river maps resultant from the SRM of real fractions associated with the SPU and the OBA-NDWI algorithms are illustrated, respectively, in Figures 22 and 23 for Sarca and Noce rivers, respectively.Obviously, the accuracies of sub-pixel User and producer accuracies of the sub-pixel river maps resultant from the SRM of real fractions associated with the SPU and the OBA-NDWI algorithms are illustrated, respectively, in Figures 22 and 23 for Sarca and Noce rivers, respectively.Obviously, the accuracies of sub-pixel User and producer accuracies of the sub-pixel river maps resultant from the SRM of real fractions associated with the SPU and the OBA-NDWI algorithms are illustrated, respectively, in Figures 22  and 23 for Sarca and Noce rivers, respectively.Obviously, the accuracies of sub-pixel maps are decreased in the presence of uncertainties related to unmixing.However, the sub-pixel maps are accurate up to 90% with ZF ≤ 5.The accuracies of sub-pixel maps associated with water fractions derived from the proposed OBA-NDWI are comparable with those obtained from the SPU.This is while the proposed method is simpler and easy to be numerically implemented and applied for water mapping issues.
As evident in Figures 12 and 13, original bands of NDWI show significantly lower correlation with water fractions than that of optimal pair of bands.This lead to about 6% lower user/producer accuracies at the final sub-pixel map using original pair of bands for the WV-3 image.Figure 24 compares the sub-pixel maps derived from OBA-NDWI using optimal pair of bands (i.e., CB and RE) against the original bands of NDWI (i.e., NIR1 and G bands).In the shown example, the sub-pixel map derived from the original pair of bands lead to omission of some water sub-pixels particularly at the boundaries.
ISPRS Int.J. Geo-Inf.2017, 6, 383 17 of 21 maps are decreased in the presence of uncertainties related to unmixing.However, the sub-pixel maps are accurate up to 90% with ZF ≤ 5.The accuracies of sub-pixel maps associated with water fractions derived from the proposed OBA-NDWI are comparable with those obtained from the SPU.This is while the proposed method is simpler and easy to be numerically implemented and applied for water mapping issues.
As evident in Figures 12 and 13, original bands of NDWI show significantly lower correlation with water fractions than that of optimal pair of bands.This lead to about 6% lower user/producer accuracies at the final sub-pixel map using original pair of bands for the WV-3 image.Figure 24 compares the sub-pixel maps derived from OBA-NDWI using optimal pair of bands (i.e., CB and RE) against the original bands of NDWI (i.e., NIR1 and G bands).In the shown example, the sub-pixel map derived from the original pair of bands lead to omission of some water sub-pixels particularly at the boundaries.Regarding the computational cost of SRM techniques, interpolation-based techniques and the MBPS algorithm are on average 20 times and three times faster than the original PS algorithm, respectively.This is because the interpolation-based and MBPS algorithms allocate sub-pixels in a non-iterative process and, also, the calculation of attractiveness for sub-pixel locations is not required for interpolation-based SRM.

Discussion
Results obtained from the SRM of semi-simulated fractions demonstrate that hard classification is suffering remarkably from the mixture in boundary pixels where the extracted border lines are very rough and are missing the detailed proximity with riparian zone.In this regard, the user and producer accuracies of the hard classified river map would be lower over than 10% compared to the SRM techniques.In general, interpolation-based techniques, as well as the proposed MBPS algorithm, produce comparable results with the PS algorithm.However, it seems that the producer accuracy of interpolation-based techniques decreases much faster than PS/MBPS from an increase in ZF.As another key point, increasing the ZF leads to a higher degree of uncertainty in allocating the fractions in sub-pixel locations.Thus, a trade-off between the accuracy of SRM algorithms and ZF should be taken into account.
The resultant real water fractions from the proposed OBA-NDWI are coherent and very close to the results of the SPU method, showing promising outputs especially for ZF ≤ 5.However, the computational complexity of estimating water fractions based on the OBA-NDWI method is significantly less than that of the SPU algorithm.Comparing the two testing approaches reveals that the spatial allocation process of sub-pixels using SRM techniques is very accurate in the context of mapping river boundaries (above 95% user/producer accuracy for ZF ≤ 6).This is while, considering the entire sub-pixel mapping process, the importance of the estimating water fractions (unmixing) seems to be more crucial than the spatial allocation of the sub-pixels.The PS and MBPS algorithms preserve the input water fractions while the fraction values cannot be respected through the interpolation-based SRM.In the discussed case studies, some isolated non-water sub-pixels appear Regarding the computational cost of SRM techniques, interpolation-based techniques and the MBPS algorithm are on average 20 times and three times faster than the original PS algorithm, respectively.This is because the interpolation-based and MBPS algorithms allocate sub-pixels in a non-iterative process and, also, the calculation of attractiveness for sub-pixel locations is not required for interpolation-based SRM.

Discussion
Results obtained from the SRM of semi-simulated fractions demonstrate that hard classification is suffering remarkably from the mixture in boundary pixels where the extracted border lines are very rough and are missing the detailed proximity with riparian zone.In this regard, the user and producer accuracies of the hard classified river map would be lower over than 10% compared to the SRM techniques.In general, interpolation-based techniques, as well as the proposed MBPS algorithm, produce comparable results with the PS algorithm.However, it seems that the producer accuracy of interpolation-based techniques decreases much faster than PS/MBPS from an increase in ZF.As another key point, increasing the ZF leads to a higher degree of uncertainty in allocating the fractions in sub-pixel locations.Thus, a trade-off between the accuracy of SRM algorithms and ZF should be taken into account.
The resultant real water fractions from the proposed OBA-NDWI are coherent and very close to the results of the SPU method, showing promising outputs especially for ZF ≤ 5.However, the computational complexity of estimating water fractions based on the OBA-NDWI method is significantly less than that of the SPU algorithm.Comparing the two testing approaches reveals that the spatial allocation process of sub-pixels using SRM techniques is very accurate in the context of mapping river boundaries (above 95% user/producer accuracy for ZF ≤ 6).This is while, considering the entire sub-pixel mapping process, the importance of the estimating water fractions (unmixing) seems to be more crucial than the spatial allocation of the sub-pixels.The PS and MBPS algorithms preserve the input water fractions while the fraction values cannot be respected through the interpolation-based SRM.In the discussed case studies, some isolated non-water sub-pixels appear on the PS and MBPS maps that applying a majority filter (mainly with a kernel size corresponding to the ZF) reduced this effect properly.The main advantage of interpolation-based and MBPS algorithms is their high computational efficiency which are, respectively, 20 times and three times faster than the PS algorithm.
The common thresholding methods, including Otsu's method and simple zero-thresholding, are applied through OBA-NDWI to investigate their effectiveness in the hard classification of water features.According to the results, Otsu's threshold value and corresponding cut-off water fraction significantly depend on the spectral bands used for the calculation of the NDWI.The minimum water fraction in Otsu's water mask is, in general, observed higher than the adequate amount (50%) for binary hard classification.The zero-thresholding suffers from both underestimation and overestimation of water pixels.This shortcoming of thresholding methods is that they can suppress the labeling of boundary mixed pixels.

Conclusions and Suggestions
Addressing the problem of mixed boundary pixels can enhance the accuracy of remote sensing applications in river sciences.Unmixing and SRM are key tools, respectively, for the spectral and spatial decomposition of the mixture in river boundaries.The former estimates the water fractions and the latter allocates the fractions in proper sub-pixel locations.Both of these steps have been the focus of this study in order to examine and develop the techniques for mapping of river boundaries at the sub-pixel resolution.Two different testing approaches are considered to survey the efficiency of SRM algorithms, while accounting for the absence and presence of uncertainty in the input data (i.e., fractions obtained from unmixing).Semi-simulated fractions are used as a contrived input with known fractions which provide a unique means of assessing the performance of the spatial allocation of sub-pixels.Moreover, the comparison between hard classification and SRM algorithms can be facilitated.This is because the hard labeling of pixels can also be applied on semi-simulated fractions without any uncertainty by deciding on the most abundant class (water/non-water) within the pixels.In the second testing approach, real water fractions are estimated based on a proposed method, namely OBA-NDWI, as well as using the fully-constrained algorithm of SPU.The rationality behind the proposed OBA-NDWI is to take advantages from the ease of use of water indices in order to estimate water fractions.The OBA-NDWI performs the NDWI with all the possible combinations of spectral bands to identify the pair of bands for which the NDWI values yield the highest correlation with water fractions.To analyze the relation of NDWI values and water fractions through OBA-NDWI, a linear mixture of image-derived spectra is used to simulate all the possible mixing fractions of selected endmembers.The proposed OBA-NDWI method permits a systematic approach to find the optimal combination of bands for calculation of NDWI which can also be instrumental for images with high spectral resolution.This benefits both the hard and soft classification of water features where the highly-correlated NDWI to water fractions can enhance the accuracy of the extraction of water features either at the pixel or sub-pixel level.The semi-simulated and real fractions are employed for SRM based on the PS algorithm, as well as interpolation-based techniques.Moreover, the MBPS algorithm is developed to provide a non-iterative alternative to the traditional PS.The accuracy and computational proficiency of the techniques are explored in two case studies using HRSI.
The proposed OBA-NDWI demonstrated that additional spectral bands of WV-3 imagery provide a couple of choices for selecting the proper pair of bands to form NDWI equation.Although the (CB, RE) pair is selected as the optimal combination of bands for the Sarca River, other combinations, such as (CB, NIR2), (B, RE), and a few other combinations also demonstrated a strong relation with water fractions.In this regard, making use of several NDWIs to establish a multiple regression for the prediction of water fractions can potentially lead to maximum benefit from spectral bands which can be an area of investigation for future studies.Furthermore, the number and quality of endmembers can affect the results of unmixing methods.The assumption regarding the three aforementioned endmembers is valid for the constitution of the major land cover types [31,36], particularly in the river boundaries.This is mainly because each of classes is representative of a major land cover which can consist of sub-classes with similar spectral responses (e.g., different vegetation types can be grouped into only one major vegetation class).However, in case of presence of an obviously different endmember in the riparian/buffer zone of the river, the number of endmembers can be increased through the unmixing process using OBA-NDWI.In this case, identification of the optimal bands would be performed locally.The OBA-NDWI method can be applied in a segment-based approach to identify the optimal bands for each desired segment of the channel.
The proposed OBA-NDWI method has the potential for selecting the proper threshold for hard classification using NDWI values and requires more investigation.OBA-NDWI can be investigated further for sub-pixel mapping of any other land/water classes, such as coastlines, road and farm boundaries.Although the focus of this study was on HRSI due to the small size of the rivers, the proposed methodology can be applied on any other sensor (e.g.Landsat-8, Sentinel-2, etc.).However, the OBA-NDWI can be applied to find the optimal pair of bands for the specific sensor.

Figure 1 .
Figure 1.True color composites of (a) WV-3 image of Sarca River and (b) the GeoEye image of Noce River; the river channels are highlighted by blue lines.

Figure 1 .
Figure 1.True color composites of (a) WV-3 image of Sarca River and (b) the GeoEye image of Noce River; the river channels are highlighted by blue lines.

Figure 2 .
Figure 2.The procedure considered for the mapping and assessment of river boundaries at the sub-pixel resolution using both real and semi-simulated fractions; ZF represents the zoom factor.

Figure 2 .
Figure 2.The procedure considered for the mapping and assessment of river boundaries at the sub-pixel resolution using both real and semi-simulated fractions; ZF represents the zoom factor.

Figure 4 .
Figure 4.The scheme for linear spectral mixture of three dominant endmembers of the land surface (i.e., water, vegetation, and soil).Each point inside the triangle represents a possible combination of fractions; adapted from Ji et al. [31].

Figure 4 .
Figure 4.The scheme for linear spectral mixture of three dominant endmembers of the land surface (i.e., water, vegetation, and soil).Each point inside the triangle represents a possible combination of fractions; adapted from Ji et al. [31].

Figure 5 .
Figure5.The step-by-step procedure for implementation of the proposed optimal band analysis for the normalized difference water index (OBA-NDWI) for the estimation of water fractions.

Figure 5 .
Figure5.The step-by-step procedure for implementation of the proposed optimal band analysis for the normalized difference water index (OBA-NDWI) for the estimation of water fractions.

Figure 6 .
Figure 6.Water fractions and the pixel swapping (PS) process for spatial allocation of sub-pixels for a given pixel; (a) water fractions, (b) random allocation of sub-pixels, and (c,d,e) swaps (candidate sub-pixels for swapping are highlighted by dash-lines).Values of sub-pixels represent their attractiveness toward the water class.

Figure 6 .
Figure 6.Water fractions and the pixel swapping (PS) process for spatial allocation of sub-pixels for a given pixel; (a) water fractions, (b) random allocation of sub-pixels, and (c,d,e) swaps (candidate sub-pixels for swapping are highlighted by dash-lines).Values of sub-pixels represent their attractiveness toward the water class.

Figure 7 .
Figure 7. Water fractions and the proposed modified binary pixel swapping (MBPS) for spatial allocation of sub-pixels for a given pixel; (a) water fractions, (b) identification of sub-pixel locations with highest attractiveness, highlighted by dash-lines, and (c) allocation of water and non-water sub-pixels.The values of sub-pixels represent their attractiveness toward the water class.

Figure 8 .
Figure 8.An example of bilinear interpolation of water fractions; (a) water fractions at the pixel level, and (b) interpolated water fractions at the sub-pixel level with ZF = 5.

Figure 7 .
Figure 7. Water fractions and the proposed modified binary pixel swapping (MBPS) for spatial allocation of sub-pixels for a given pixel; (a) water fractions, (b) identification of sub-pixel locations with highest attractiveness, highlighted by dash-lines, and (c) allocation of water and non-water sub-pixels.The values of sub-pixels represent their attractiveness toward the water class.

Figure 7 .
Figure 7. Water fractions and the proposed modified binary pixel swapping (MBPS) for spatial allocation of sub-pixels for a given pixel; (a) water fractions, (b) identification of sub-pixel locations with highest attractiveness, highlighted by dash-lines, and (c) allocation of water and non-water sub-pixels.The values of sub-pixels represent their attractiveness toward the water class.

Figure 8 .
Figure 8.An example of bilinear interpolation of water fractions; (a) water fractions at the pixel level, and (b) interpolated water fractions at the sub-pixel level with ZF = 5.

Figure 8 .
Figure 8.An example of bilinear interpolation of water fractions; (a) water fractions at the pixel level, and (b) interpolated water fractions at the sub-pixel level with ZF = 5.

21 Figure 9 .
Figure 9.The regressions of normalized difference water index (NDWI) values with different band combinations obtained from synthetic WV-3 spectra against known water fractions; the pair of bands used for calculation of NDWI is indicated on each graph (CB and RE stand for coastal-blue and red-edge bands, respectively); zero-threshold and Otsu's threshold are illustrated respectively with blue (dashed) and red (dot-dashed) lines.

Figure 10 .
Figure 10.The regressions of NDWI values with different band combinations obtained from synthetic GeoEye spectra against known water fractions; the pair of bands used for the calculation of NDWI is indicated on each graph; zero-threshold and Otsu's threshold are illustrated, respectively, with blue (dashed) and red (dot-dashed) lines.

Figure 9 . 21 Figure 9 .
Figure 9.The regressions of normalized difference water index (NDWI) values with different band combinations obtained from synthetic WV-3 spectra against known water fractions; the pair of bands used for calculation of NDWI is indicated on each graph (CB and RE stand for coastal-blue and red-edge bands, respectively); zero-threshold and Otsu's threshold are illustrated respectively with blue (dashed) and red (dot-dashed) lines.

Figure 10 .
Figure 10.The regressions of NDWI values with different band combinations obtained from synthetic GeoEye spectra against known water fractions; the pair of bands used for the calculation of NDWI is indicated on each graph; zero-threshold and Otsu's threshold are illustrated, respectively, with blue (dashed) and red (dot-dashed) lines.

Figure 10 .
Figure 10.The regressions of NDWI values with different band combinations obtained from synthetic GeoEye spectra against known water fractions; the pair of bands used for the calculation of NDWI is indicated on each graph; zero-threshold and Otsu's threshold are illustrated, respectively, with blue (dashed) and red (dot-dashed) lines.

Figure 11 .Figure 12 .
Figure 11.Minimum water fractions corresponding to zero-threshold and Otsu's threshold for several NDWIs using synthetic WV-3 and GeoEye spectra; CB and RE stand, respectively, for coastal-blue and red-edge bands.

Figure 11 .
Figure 11.Minimum water fractions corresponding to zero-threshold and Otsu's threshold for several NDWIs using synthetic WV-3 and GeoEye spectra; CB and RE stand, respectively, for coastal-blue and red-edge bands.

Figure 11 .Figure 12 .
Figure 11.Minimum water fractions corresponding to zero-threshold and Otsu's threshold for several NDWIs using synthetic WV-3 and GeoEye spectra; CB and RE stand, respectively, for coastal-blue and red-edge bands.

Figure 12 .
Figure 12.OBA-NDWI for the synthetic WV-3 spectra; all possible combinations of spectral bands are considered in the structure of NDWI to perform a second-order regression against water fractions.Values of (a) R 2 and (b) RMSE are represented by color bars.

Figure 13 .
Figure 13.OBA-NDWI for the synthetic GeoEye spectra; all possible combinations of spectral bands are considered in the structure of NDWI to perform a second order regression against water fractions.Values of (a) R 2 and (b) RMSE are represented by the color bars.

Figure 14 .
Figure 14.Water fractions of the WV-3 image obtained from (a) the proposed OBA-NDWI method and (b) the SPU algorithm.

Figure 13 .
Figure 13.OBA-NDWI for the synthetic GeoEye spectra; all possible combinations of spectral bands are considered in the structure of NDWI to perform a second order regression against water fractions.Values of (a) R 2 and (b) RMSE are represented by the color bars.

Figure 13 .
Figure 13.OBA-NDWI for the synthetic GeoEye spectra; all possible combinations of spectral bands are considered in the structure of NDWI to perform a second order regression against water fractions.Values of (a) R 2 and (b) RMSE are represented by the color bars.

Figure 14 .
Figure 14.Water fractions of the WV-3 image obtained from (a) the proposed OBA-NDWI method and (b) the SPU algorithm.

Figure 14 .
Figure 14.Water fractions of the WV-3 image obtained from (a) the proposed OBA-NDWI method and (b) the SPU algorithm.
ISPRS Int.J. Geo-Inf.2017, 6, 383 13 of 21classification is very rough on the river boundaries, while the SRM techniques reconstruct boundaries with sub-pixel details.

Figure 18 .
Figure 18.User and producer accuracies of hard classification and SRM algorithms using semi-simulated water fractions across a range of ZF for Noce River; (HC: hard classification, PS: pixel swapping, MBPS: modified binary PS, BL: bilinear, BC: bicubic, L3: lanczos3).

Figure 19 .
Figure 19.Error maps of (a) hard classification and (b) sub-pixel map obtained from MBPS for a segment of the Sarca River using the WV-3 image with ZF = 6; red and blue pixels show erroneously committed and omitted water pixels, respectively.

Figure 18 .
Figure 18.User and producer accuracies of hard classification and SRM algorithms using semi-simulated water fractions across a range of ZF for Noce River; (HC: hard classification, PS: pixel swapping, MBPS: modified binary PS, BL: bilinear, BC: bicubic, L3: lanczos3).

Figure 19 .
Figure 19.Error maps of (a) hard classification and (b) sub-pixel map obtained from MBPS for a segment of the Sarca River using the WV-3 image with ZF = 6; red and blue pixels show erroneously committed and omitted water pixels, respectively.

Figure 18 .
Figure 18.User and producer accuracies of hard classification and SRM algorithms using semi-simulated water fractions across a range of ZF for Noce River; (HC: hard classification, PS: pixel swapping, MBPS: modified binary PS, BL: bilinear, BC: bicubic, L3: lanczos3).

Figure 18 .
Figure 18.User and producer accuracies of hard classification and SRM algorithms using semi-simulated water fractions across a range of ZF for Noce River; (HC: hard classification, PS: pixel swapping, MBPS: modified binary PS, BL: bilinear, BC: bicubic, L3: lanczos3).

Figure 19 .
Figure 19.Error maps of (a) hard classification and (b) sub-pixel map obtained from MBPS for a segment of the Sarca River using the WV-3 image with ZF = 6; red and blue pixels show erroneously committed and omitted water pixels, respectively.

Figure 19 .
Figure 19.Error maps of (a) hard classification and (b) sub-pixel map obtained from MBPS for a segment of the Sarca River using the WV-3 image with ZF = 6; red and blue pixels show erroneously committed and omitted water pixels, respectively.

Figure 23 .
Figure 23.User and producer accuracies of SRM algorithms using real water fractions of: (a,b) the SPU algorithm, and (c,d) the OBA-NDWI algorithm across a range of ZF for Noce River; (PS: pixel swapping, MBPS: modified binary PS, BL: bilinear, BC: bicubic, L3: lanczos3).

Figure 24 .
Figure 24.Sub-pixel maps derived from OBA-NDWI using (a) optimal pair of bands (CB and RE) against (b) original pair of bands (NIR1, G) for the WV-3 image.

Figure 24 .
Figure 24.Sub-pixel maps derived from OBA-NDWI using (a) optimal pair of bands (CB and RE) against (b) original pair of bands (NIR1, G) for the WV-3 image.

Table 1 .
Multispectral bands of GeoEye and WV-3 sensors and list of the sensors providing same/similar spectral resolutions.