River Detection in Remotely Sensed Imagery Using Gabor Filtering and Path Opening

Detecting rivers from remotely sensed imagery is an initial yet important step in space-based river studies. This paper proposes an automatic approach to enhance and detect complete river networks. The main contribution of this work is the characterization of rivers according to their Gaussian-like cross-sections and longitudinal continuity. A Gabor filter was first employed to enhance river cross-sections. Rivers are better discerned from the image background after filtering but they can be easily corrupted owing to significant gray variations along river courses. Path opening, a flexible morphological operator, was then used to lengthen the river channel continuity and suppress noise. Rivers were consistently discerned from the image background after these two-step processes. Finally, a global threshold was automatically determined and applied to create binary river masks. River networks of the Yukon Basin and the Greenland Ice Sheet were successfully detected in two Landsat 8 OLI panchromatic images using the proposed method, yielding a high accuracy (~97.79%), a high true positive rate (~94.33%), and a low false positive rate (~1.76%). Furthermore, experimental tests validated the high capability of the proposed method to preserve river network continuity.


Introduction
The quantification of river features, such as length, width, tortuosity, branching pattern, and temporal variation, among others, can provide useful insights into surface water supply, transport, distribution, and dynamics [1].The study of these river characteristics is usually performed via field work [2].However, rivers are widely distributed all over the world [3].It is commonly time-consuming and expensive to conduct river measurement in the field.In particular, it is unsafe or even impossible to measure rivers in certain environments, such as ice sheets [4], tidal flats [5], and floodplains [6].Given this situation, at present, satellite observations of the earth are playing substantial roles in surface river studies [7].Rivers, which are open water bodies, have very low visible/near-infrared spectral reflectance and thereby appear as dark, continuous, and curvilinear features in the imagery [8,9].These spectral and shape properties make rivers unique features in optical visible/near-infrared remotely sensed imagery.Detecting rivers efficiently will greatly contribute to river studies and improve our understanding of surface water on earth.
The detection of linear or curvilinear features is a classic topic in image processing studies.A considerable amount of research has focused on detecting road networks, shorelines, and blood vessels from both optical and synthetic aperture radar (SAR) images [10][11][12][13][14][15].Detection of river features, in general, has received little attention.This is because the formation of river networks is mainly controlled by surface topography and consequently digital terrain models (DTMs) have become the primary data sources to determine drainage networks via contributing area methods, area-slope methods, and grid network ordering methods [16].Some global or regional river network products have been released based on several freely available DTM data [17][18][19].Optical visible/near-infrared remote sensing images, on the other hand, are usually employed as secondary or supplementary data to validate the DTM-derived drainage networks.However, the DTM-modeled drainage networks sometimes have significant differences with the actual observed river networks [20].At least four reasons can be identified for this: first, the accuracies and spatial resolutions of the DTM data are limited [21,22]; second, the snapshot DTM data fail to capture temporal topographic changes in highly dynamic environments (e.g., ice sheets) [23,24]; third, other factors, such as tectonic faults [25], crevasses [26], and moulins (draining surface meltwater) [4], may impact river network patterns besides topography; and fourth, it is commonly acknowledged that the patterns of the DTM-derived drainage networks depend on the threshold adapted when using certain DTM processing approaches, often leading to headwater rivers being under-or over-estimated [22,27].Moreover, apart from river lengths and orders, very limited quantification information can be generated from DTM-derived drainage networks.As such, optical remotely sensed imagery should play more important roles in river studies.It is critical to detect rivers from remotely sensed imagery and generate river masks, from which river widths can be directly measured [28], and subsequently some other key information (e.g., river depth, flow velocity, and discharge) may be further determined.
Recently, both remote sensing and hydrology research communities have begun to highlight the importance of studying rivers from space [29,30].A new research subject named fluvial remote sensing has grown at a rapid rate in recent years [31,32].A few published papers have noticed the limitation of DTM-derived drainage networks and attempted to detect rivers from optical and SAR images [33].Among these studies, manual interpretation [6], semi-automatic [34], and ancillary-data-supported [35,36] methods are the three primary detection approaches used.However, these methods may not be feasible for operational use, especially when applied to map rivers in large regions.Increasing the level of automation in river detection is, thus, crucial.An automated method is proposed in [37] to detect rivers in high-resolution SAR images, employing adaptively created training samples and supervised image classification.The method of supervised image classification is also used to delineate rivers in [38], focusing on the performance of a support vector machine (SVM) classifier.These automatic methods show good success for delineating wide rivers.However, much less attention has been paid to identifying thin rivers (width < 5 pixels) in images.These thin rivers are usually headwater channels in networks and represent very lower spectral contrast compared to the image background.Two studies independently noticed the limitations of studying rivers at mono-scale and proposed two multi-scale classification approaches to detect rivers [39,40].However, the spectral contrast along thin river courses is usually unstable owing to variations in river depth, sediment load, and channel bank material.As an example, several river networks distributed in the Yukon River Basin, northwestern North America, are shown in Figure 1a.In this study area, headwater river channels are difficult to discern from the image background even after multi-scale analysis.This low and unstable contrast that occludes rivers in image is a common problem that complicates river enhancement and renders subsequent river detection unreliable.To enhance thin rivers appropriately, it is crucial to understand the characteristics of rivers in remotely sensed imagery.
The main expected characteristics of rivers in optical visible/near-infrared remotely sensed imagery can be summarized as follows: first, rivers appear as curvilinear features instead of blob-like features; second, rivers exhibit Gaussian-shaped cross-sections; and third, rivers are continuous and not a series of disconnected curvilinear segments.Three-dimensional view of the example river networks presents these characteristics very well (Figure 1b).Utilizing the above river properties, this paper proposes an automatic approach to enhance and detect rivers in remotely sensed imagery.The main contribution of this work is to characterize rivers according to their Gaussian-like cross-sections and longitudinal continuity.Combining the cross-sectional and longitudinal river information in optical remotely sensed imagery, rivers are successfully discerned from noisy image background.
The rest of this paper is organized as follows.Section 2.1 describes the procedures for pre-processing river images.Section 2.2 describes the implementation of Gabor filtering to enhance river cross-sections.Section 2.3 describes the implementation of path opening to enhance river longitudinal continuity.Section 2.4 describes the method used to classify the enhanced image to obtain a binary river map using an automatically determined threshold.Section 3 presents the experimental results.Section 4 concludes the paper.

Pre-Process
Rivers in remotely sensed images often show gray variations, poor contrast, and noise.This, along with the fact that small rivers tend to lack distinction, means that pre-processing is required.The main idea is to ease the river detection process by both reducing noise and enhancing the dark curvilinear features (possible cracks).A 3 × 3 mean filter was first applied to reduce salt and pepper noise [41].
Local contrast enhancement was then achieved by applying contrast limited adaptive histogram equalization (CLAHE) [42].Next, shade correction was performed by subtracting an image approximating the background.This approximation was obtained by applying a mean filter with a 50 × 50 pixel size kernel.This large size was chosen to ensure river regions were well preserved.The gray values of the filtered image were then inversed to represent rivers as bright features, as shown in Figure 2b.This pre-processed image was used as the input for the subsequent enhancement and detection processes.A 3D view of the pre-processed image is shown in Figure 3b; it shows that river features are better discerned from the image background as compared to the original image shown in Figure 3a.

Gabor Filtering
The primary difficulty in detecting rivers from remotely sensed imagery is that thin rivers commonly exhibit low spectral contrast with the image background (Figure 3a).Therefore, linear or curvilinear features are required to be enhanced before performing detection.Many approaches have been developed to achieve this goal.For example, different linear detectors have been proposed to enhance roads, cracks, and vessels in images [43,44].A primary limitation of these detectors is that their designs are highly dependent on the geometry of the features to be enhanced and the image types, leading to limited application ability.A matched filter is a more flexible linear feature detector [45].This technique is simple yet effective and hence has been widely used for vessel detection.However, a matched filter will respond not only to linear features but also to edges.This will cause incorrect enhancement and frequently create false linear features [46].Recently, the Gabor filter has been shown to perform better than the matched filter and has been applied to enhance blood vessels in image [47].In addition, it enhances curvilinear features and suppresses noise and background.We note that the appearance of single channel rivers in remotely sensed imagery resembles, to some degree, the appearance of vessels in biomedical images.In particular, both structures are long curvilinear features with variable elongated intensities.Hence, river enhancement and detection problems stand to benefit from methods created for the detection and analysis of vessel structures.For this reason, the Gabor filter was employed in this study to enhance the appearance of rivers.
The Gabor filter is a sinusoidally modulated Gaussian function that has optimal localization in both the frequency and space domains [48].The Gabor filter is effective for enhancing rivers because it exploits the fact that the cross-sections of rivers in remotely sensed imagery are Gaussian shaped.This idea is adapted here for identifying rivers: although thin rivers contrast poorly with the background, they nevertheless exhibit Gaussian-shaped cross-sections in a pre-processed image (Figure 3b) and therefore, their cross-sections can be enhanced using the Gabor filter.The real Gabor filter kernel oriented at an angle θ π/2 may be formulated as: ( This filter rotation offers an ability to detect rivers with varied directions.The response of a Gabor filter to an input image f (x, y) can be expressed by the following equation: The final maximum response of all orientations satisfies: To detect the thin rivers in the images (width < 5 pixels), w was set to 2 in this study.In addition, the angular interval was set to 15°, and 12 different kernels were thus developed to span all the possible orientations.The Gabor filter was calculated by using a 5 × 5 sliding window and the final response was allocated to the center pixel of this window.Mono-scale (w = 2) was used instead of multi-scale (e.g., [40]) because path opening is subsequently applied to improve longitudinal continuity, discerning both wide and thin rivers from image background (see Section 2.3).
The resulting image after Gabor filtering is shown in Figure 2c, which shows a better contrast of the linear characteristics of rivers, especially the headwater channels, with the image background.The corresponding 3D view in Figure 3c better highlights this point, illustrating that rivers are better discerned from the background than in Figure 3b.

Path Opening
Owing to significant longitudinal spectral variations, thin elongated rivers can be easily corrupted even after enhancing cross-sections, resulting in disconnected river segments.Figure 3c shows that considerable intensity variations along river courses still exist after Gabor filtering.To generate complete river networks, a common process is to first classify the enhanced image using a global threshold and then close the gaps along the river courses [49].However, this process has four considerable drawbacks: first, most linking methods are only suitable for linear features with low sinuosity or short discontinuous segments, lacking an appropriate ability to connect highly sinuous river segments with long gaps; second, the classified image may miss very small rivers before connection, thus underestimating actual existing rivers; third, the linked gaps are only one pixel wide, and lose information on the river width [50]; and fourth, linking methods are usually computationally complex, including endpoint detection, direction determination, and some post-processing steps.Given these limitations, instead of linking broken river gaps, it is better to "stabilize" the pixel values of the river longitudinal profiles before classification.More specifically, the pixel values along an individual river segment should be more similar to each other to improve river course continuity and facilitate the subsequent river detection.To this end, a state-of-the-art morphological filter, the path opening filter, was adapted in this study.
The path opening, a flexible morphological operator, was proposed in the last few years [51].This operator retains oriented and local paths brighter than their immediate surroundings.This transformation does not filter the image by the size of the components only (as does the surface area opening) [52].Instead, it filters curvilinear paths with specified size and orientation.The search for the paths is carried out in the 0°, 45°, 90°, and 135° directions to identify neighboring pixels, and therefore, the paths in these directions are all enhanced.Path opening studies the gray level scale in the image, and assigns to each pixel the highest gray level, where a path with the required length Lmin is formed.A primary advantage of the path opening filter is that it can enhance curvilinear features with flexible shapes (e.g., highly sinuous rivers).The parameter Lmin can be interpreted as the structuring element size for the pathopening operator.This makes the length an important variable in the classification, and rivers should therefore be linked together as much as possible to be distinguished from the image background.The path opening filter explores all the potential paths in an image, leading to low efficiency when applied to large remotely sensed imagery.To address this problem, a newly proposed parsimonious path opening (PPO) was employed in this study [53].This operator only identifies a subset of the paths considered by the conventional path opening, thus achieving a substantial speed-up while obtaining very similar results as the complete path opening operator.
The resulting image after the path opening process (Lmin = 40 pixels) is shown in Figure 2d; it shows that the pixel values along each river segment are much more stable than in the Gabor filtered image shown in Figure 2c.A comparison of the 3D views in Figure 3c,d also better highlights this result.It is clear that both the mainstream and headwater channels are stabilized and preserved while the image background is suppressed as well.The length parameter, Lmin, is a key factor in the path opening process.As Lmin increases, the paths that pass this length threshold are enhanced and stabilized while those that fail to pass are suppressed (Figures 5 and 6).
An example of the PPO processed images under different length threshold Lmin together with the corresponding histograms are shown in Figure 5.When Lmin = 5, many short false segments are preserved while the gray intensities along many actual river segments vary significantly; as Lmin increases to 40, the majority of the short noise paths are suppressed well while the actual river segments become more continuous.However, if Lmin is set too high, some actual river segments may fail to pass the Lmin threshold and will, thus, be suppressed.For example, when Lmin = 120, the small river network at the upper-left corner of the image is suppressed; when Lmin = 280, only the largest river network in the middle of the image is preserved.When Lmin is even larger, this largest network starts to expand, seizing the adjacent pixels that do not belong to the actual river segments and creating a large number of false paths.The 3D views in Figure 6 further highlights the impact of Lmin on the longitudinal enhancement of rivers.
Moreover, the histograms in Figure 5 indicate that the number of gray levels in the processed images decreases correspondingly as Lmin increases, suggesting a stronger "grouping" or "stabilizing" performance by the path opening filter.On the other hand, fewer gray levels limit the ability to distinguish features in the enhanced image.In summary, we suggest that a practical Lmin parameter be used to achieve a good balance between actual river enhancement and false path suppression.Following this rule, Lmin was set to 40 pixels in this study based on experimental tests.

Thresholding
After Gabor filtering and path opening, the rivers show significant differences with the image background (Figures 2d and 3d).Therefore, a simple single-thresholding operation can effectively distinguish rivers from the image background and generate a binary river map.To test the performance of global thresholding, an example threshold was imposed on the original, pre-processed, Gabor-filtered, and final-enhanced images.As illustrated in Figure 7, rivers are erroneously widened in both the original and the preprocessed images, leading to many false positives.Meanwhile, the small headwater rivers are not detected in both the images (Figure 7a,b).After Gabor filtering, curvilinear rivers are successfully detected; however, a large number of rivers are disconnected and some short noisy paths are obtained as well (Figure 7c).Path opening addresses this problem and produces the final classified river mask, in which rivers are continuous while the small headwater channels are successfully detected (Figure 7d).It can be clearly seen that the rivers are better delineated when considering both cross-sectional and longitudinal information.
The global threshold, T, used to classify the enhanced image should be determined carefully.A low T will overestimate river pixels and cause many false positives, while a high T will underestimate river pixels and eliminate actual small rivers.Figure 8 shows the histogram of the final enhanced image and the binary river maps based on three different T values.The result illustrates that if the gray value of a pixel is smaller than the mean value of the enhanced image, that pixel belongs to non-river regions.Subsequently, if the pixel value is larger than the mean + one standard derivation (1 std) value of the image, it belongs to river regions with very high confidence; and finally, if the value is between the mean and mean + 1 std, it has the potential to be classified as a river pixel.For example, the small rivers in the three red circles are assumed to belong to this category, which are not detected when T = mean + 1 std, but successfully identified when T = mean.In the present study, we empirically found that the setting T = mean + 0.5 std achieved a good balance between river detection and noise suppression (Figure 8).The primary advantage of this approach is that this global threshold only depends on the final enhanced image and hence can be determined automatically.

Implementation
The proposed method is named the Gabor-PPO method and several river detection experiments were conducted to estimate the performance of this method.Pre-process and Gabor filtering steps were implemented in Matlab and the PPO coding was obtained from [48].The thickness of the Gabor detector, w, was set to 2. The length parameter of PPO, Lmin, was set to 40 pixels.The global threshold, T, was set to the value of mean + 0.5 std of the enhanced image.

Study Areas
River networks distributed in the Yukon River Basin and the Greenland Ice Sheet (GrIS) ablation zone were employed as examples (Figure 9).These two types of river networks distributed in the Arctic were selected because they exhibit complex morphology and have been experiencing large pattern variations in recent years.Importantly, they show the responses of surface river systems to global climate change.High drainage densities and numerous thin rivers can be used to comprehensively assess the performance of river detection methods.In this study, two subsets of Landsat 8 OLI (Operational Land Imager) panchromatic images from the above two areas were used as test data.The image sizes for the Yukon and GrIS study areas were 1095 × 1269 and 1439 × 1153 pixels.Landsat 8 was launched in 2013 and the images captured from this satellite show considerable improvements over those from its predecessors (Landsat 5 and 7).It is widely acknowledged that these types of images will play substantial roles in earth observations in the future.Moreover, panchromatic images offer a relative high resolution (~15 m) and wide spatial coverage (~185 km), and therefore facilitate detailed river studies.
Two river images are shown in Figure 10a,b as examples, which show that the river networks are well developed and organized in the two study areas and that they include numerous thin rivers.In addition, the Yukon networks flowed in two different directions (to the upper and lower, respectively), while the GrIS networks flowed in one primary direction (from left to right).Furthermore, several river networks of similar sizes are developed in the Yukon Basin, while a very large supraglacial river network is developed on the GrIS surface, surrounded by several considerably smaller networks.
To evaluate the performance of the automatic river detection method, three typical test sites were selected for each study area, as shown in Figure 10.For the Yukon networks, Site 1 represents a well-developed dendritic network; Site 2 represents a parallel network; and Site 3 represents a network with meandered main channels.For the GrIS networks, Site 4 includes an independent small network at low elevations; Site 5 includes a small part of the largest network in the middle of the image; and Site 6 includes a small network at high elevations.In summary, these test sites were selected in a way to thoroughly assess the performance of the proposed method to detect river networks with different drainage patterns at various locations.Note that the proposed method focused on linear feature detection and did not detect lakes.

Detection Accuracy Evaluation
The final binary river maps were generated using the proposed method.The detection results are shown in Figure 10c,d.A simple qualitative comparison illustrates that the proposed method detected the rivers from the images with good success.Thin rivers were successfully discerned from the image background and the network continuity was well preserved.The different drainage patterns in the two study areas are described well using the classified river masks.
To further evaluate the performance of the proposed method, binary river maps were manually digitized for the six test sites and used as the ground truth.The classification results obtained by thresholding the original panchromatic image, the PPO processed image, and the Gabor filtered image were used as comparison groups.The global threshold, T, was also set to the value of mean + 0.5 std for the enhanced images.The river detection results for the six test sites are shown in Figure 11.
The classified panchromatic and PPO processed images apparently widened the actual rivers, failed to detect a large number of thin rivers, and were very sensitive to image noise.For example, large patches of dark regions were erroneously classified as rivers at Sites 2, 5, and 6; in addition, many thin rivers were not detected at Sites 1, 3, and 4.This poor performance is mainly because both the panchromatic and the PPO-processed images lacked linear feature enhancing.As a result, river pixels were mixed with the noisy background in these two images.After applying the Gabor filters, the linear river characteristics were enhanced well and the detection results agreed much better with the ground truth maps.However, two problems exist in the Gabor filtered images: first, very short linear noise segments also yielded large responses to the Gabor filters, causing many erroneous short river segments (see the detection result of Site 5); second, many actual rivers were classified as discontinuous segments owing to the largely varied responses to the Gabor filters along the river courses (see the detection result of Site 2).After incorporating the path opening process, these two problems were addressed well.Short noise features failed to pass the length threshold and were, therefore, not classified as rivers.The weak Gabor-filtered responses along the river courses were successfully enhanced through path propagations, yielding continuous river features.
Furthermore, the good performance of the proposed method is demonstrated by the quantitative comparison and statistics.The outcomes of these river detections were pixel-based classification results.Any given pixel was classified either as representing an area of river or the background.Consequently, there were four possibilities: true positive (TP), true negative (TN), false negative (FN), or false positive (FP).In this study, the true positive rate (TPR = TP / (TP + FN)), the false positive rate (FPR = FP / (TN + FP)), and the corresponding accuracy (Acc = (TP + TN) / N), were calculated to evaluate the performance of the procedure.N was the pixel count of the image.The TPR represents the fraction of pixels correctly identified as river pixels.The FPR is the fraction of pixels erroneously identified as river pixels.The Acc was calculated as the ratio of the total number of correctly classified pixels to the number of pixels in the images.
The accuracy evaluation results are shown in Tables 1 and 2. The results show that the proposed method outperformed the other comparison groups in TPR, FPR, and Acc.The average Acc value of the proposed method for all the six test sites was 97.79%, which is larger than the 88.42%, 89.18%, and 94.65% obtained by thresholding the original panchromatic image, the PPO processed image, and the Gabor filtered image, respectively.From the perspective of the average TPR values, the proposed method achieved a high value of 94.33%, indicating that a majority of the actual river pixels were successfully detected; meanwhile, the average TPR values for the three comparison groups were only 77.76%, 75.02%, and 84.02%, implying a significant decrease in actual river detection ability.From the perspective of the average FPR values, the proposed method achieved a low value of 1.76%, indicating a good ability to eliminate noise; meanwhile, the original panchromatic image and the PPO processed image performed poorly to suppress noise, yielding very high FPR values (10.18% and 8.72%, respectively).The majority of false positives were mainly caused by some weak noise features that were darker than the image background and were also connected to actual rivers.These noise features yielded relatively large responses to the Gabor filters and were then further enhanced by the PPO operation owing to their connections to the actual long rivers.The red circles in Figure 11 show these erroneously detected features.It should be noticed that some of these features actually have certain river-like characteristics and imply channel formation processes.For example, three short segments in Site 4 were detected by the proposed method but were not accepted as rivers in the manual investigation; however, judging from their shape and spectral characteristics in the panchromatic image, they may belong to slush pathways on the ice sheet, which may evolve into open channel flows after sufficient meltwater supply.
Table 1.Performance evaluation of the classification images when one global threshold (value: mean + 0.5 std) is applied to the original panchromatic image, the PPO processed image, the Gabor filtered image, and the Gabor + PPO processed (the proposed method) image, compared with the manual digitization for the Sites 1-3 in the Yukon image.The actual rivers that the proposed method failed to detect mainly belonged to the following two cases: first, some river features were strongly confused with the image background and thereby were not detected (see the yellow circles in Site 1 and 2, Figure 11); and second, some standalone river segments failed to pass the path-opening length threshold and were not detected automatically (see the yellow circle in Site 6, Figure 11).

Image
In short, it can be clearly seen that the rivers were better detected when considering both cross-sectional and longitudinal information.The proposed method outperformed all the other comparison groups and detected rivers with good accuracy in both study areas.

Impacts of Global Thresholding
The binary river map obtained from a single global threshold, T, may not indicate the performance of the enhancing and detecting approaches comprehensively.Hence, it is necessary to analyze the impacts of varied T values on the detection accuracy.For this purpose, Receiver Operating Characteristics (ROC) graphs were obtained by plotting the derived TPRs against FPRs under different T values [54].Each point in the resultant ROC graphs represents the detection performance under one global threshold.Points in ROC space appear to be better as we proceed to the northwest (TPR is higher; FPR is lower, or both).As such, ROC graphs provide a good means to describe the impacts of global thresholding on river detection.
The ROC curves generated from the original panchromatic image, the PPO processed image, the Gabor filtered image, and the Gabor-PPO processed image are shown in Figure 12.The results show that the proposed method outperformed the three other comparison groups for six test sites.Their ROC curves were the most to the northwest of the graph spaces.The Gabor filtered image performed a little worse than the Gabor-PPO processed image, both performing much better than the original panchromatic image and the PPO processed image, especially at Sites 2, 5, 6, where large patches of false positives were obtained from the latter two images.
In the ROC graphs, the TPRs increased with increasing FPR simultaneously.When T was high, low FPR and TPR values were obtained; when T was low, high FPR and TPR values were obtained.Hence, the performance of one global threshold may be evaluated by considering the tradeoff between the benefits (true positives) and costs (false positives).The common indicator to evaluate this tradeoff is P = TPR -FPR.The corresponding global threshold Tp that yielded the largest P value was considered as the optimal threshold to classify images.For all the six test sites, Tp values were very close to the global threshold T that was employed to classify river images in Section 3.3.This quantitatively illustrated that setting T = mean + 0.5 std achieves good river detection results.Notably, as shown in Figures 8 and 12, the values of T greatly control the patterns of the river networks.A larger T only captured the apparent rivers in the images, while a smaller T managed to detect small rivers.Therefore, on the basis of the reliability of the detected rivers, different T values could be used to map rivers in remotely sensed imagery.Setting a large T was "conservative", only resulting in positive identifications with strong evidence (low FPRs and TPRs), while setting a small T was "unconservative," capturing more small rivers and possibly more noise (high FPRs and TPRs).

River Network Continuity Evaluation
Apart from detecting more actual river pixels and achieving simultaneous noise elimination (i.e., higher Acc, higher TPR, and lower FPR), a good detection method should retain a river network's continuity.This continuity greatly facilitates the subsequent river analysis (e.g., ordering, morphology analysis, etc.).Therefore, the ability of Gabor filtering and the proposed method to retain river network continuity was further evaluated here.
The detected river networks in the two study areas were labeled into individual networks, as shown in Figure 13.These labeled network maps indicates that the networks derived from the Gabor filtered image were very fragmented (Figure 13a,b).Numerous short river segments were labeled as independent networks; meanwhile, large networks were erroneously divided into several smaller networks.For example, the largest networks in the two study areas were labeled with the red color.It can be clearly seen that the largest Yukon network was separated into three primary networks, along with many small networks.Although the continuity of the largest GrIS network was preserved well, owing to its strong contrast with the image background, some headwater channels were still fragmented and labeled into independent networks in the Gabor-filtered image.Incorporating PPO processing, this fragmentary label problem was efficiently addressed.As shown in Figure 13c,d, the continuity of the large river networks was preserved well; meanwhile, most of the small river segments were successfully eliminated.
Table 3 quantitatively describes the river network continuity.The Gabor filtered image yielded many more river networks than the Gabor-PPO processed image.Given the low FPRs for both images (Tables 1 and 2), it is reasonable to assume that the majority of the detected pixels were actual river pixels.Therefore, better river network continuity is illustrated by fewer networks and more river pixels.Average pixel count per river network (Ave) was thus used as an indicator to evaluate the network continuity.The Ave value for the Gabor-PPO processed Yukon (GrIS) image was 814 (433) pixels, which is 25.2 (26.4) times more than that for the corresponding Gabor filtered image.This result further demonstrated the good ability of the proposed method to simultaneously preserve river network continuity and eliminate noise.
Figure 13 suggests that the majority of the labeled networks in the Gabor filtered image are short river segments.The use of a simple size threshold may eliminate the small networks and label these networks more accurately.Therefore, to further compare the ability of the Gabor filter and the Gabor-PPO processing, different size thresholds were employed to eliminate small river networks.Figure 14 shows the count of the networks and the river pixels under different size thresholds.The results suggest that the Gabor network counts decreased greatly with the increase of the size thresholds, while the decrease trends for the Gabor-PPO image were smooth (Figure 14a,c).Furthermore, although similar network counts were obtained when the size threshold was larger than 200 pixels, the number of the Gabor-PPO river pixels were considerably more than that of the Gabor river pixels (Figure 14b,d).Furthermore, the number of Gabor-PPO networks became more than that of the Gabor networks when the size threshold was larger than 800 pixels.This is because the continuity of the Gabor networks was relatively low and some networks failed to pass this large size threshold.These networks expanded to upstream channels or constituted larger networks in the Gabor-PPO processed image and, therefore, survived after size thresholding.

Discussion and Conclusion
The river detection procedures presented here have four areas for further studies.First, two parameters, the length threshold Lmin for the path opening and the global threshold T for the final river classification, were determined by experience in this study.For the areas where river networks are well developed and landscapes are uniform (e.g., ice sheet), we recommend using T = mean + 0.5 std (see Section 2.4).However, although these empirical parameters performed well in experimental tests, adaptive selections would be more useful.More specifically, these two parameters should be determined adaptively according to the characteristics of river images (e.g., river type, drainage pattern, drainage density, and temporal variation).Second, ROC graphs were used to analyze the overall sensitivity of T values.A possible limitation of ROC graphs is that they cannot show the impacts of T values at different locations of the tested image [54].If the tested river landscapes are not uniform and different locations have varied responses to river detection, adaptive sensitivity analysis should be conducted.Third, multiband images instead of single-band images can be used to detect rivers in future studies.Fusing of river enhancement results obtained from multiple bands and building band ratio map may facilitate more detailed river mapping.Fourth, some river-like features were erroneously enhanced and detected by the proposed methods, causing some false positives.These interference features should be further studied based on the physical environment that the remotely sensed imagery implies.
The well-developed thin rivers distributed in the Yukon River Basin and the Greenland Ice Sheet were selected to conduct river detection experiment in this study.In reality, rivers run through nearly every landscape on earth, and thus there are various types of river channels, ranging from braided channels through meandering and straight to anastomosing channels [40].In addition, the widths of rivers vary significantly along their courses, ranging from less than 1 m to more than 10 km in certain places [55].Furthermore, aside from rivers, lakes are also very important components of surface water features.Therefore, it is crucial to extend the proposed method to detect surface water features (including rivers and lakes) distributed at varied landscapes on earth.We suggest that the multi-scale analysis proposed in [39] and [40] can be adopted to extend the proposed method.More specifically, lakes and wide rivers can be identified using the normalized difference water index (NDWI) method and its variations [4,9]; in contrast, thin rivers can be detected by using the proposed method.Finally, surface water coverage can be obtained by combining these two detection results.This extended approach will be very useful when applied to detect river networks with a wide mainstream and numerous associated narrow tributaries, or to map surface water features in lake-rich hydrologic environments.
The proposed method outperforms the multi-scale approach described in [40].The proposed method can efficiently enhance and detect continuous thin rivers.However, discontinuous segments along the courses of thin rivers were obtained after multi-scale matched filtering.In [40], a multi-points fast marching method was employed to rejoin these discontinuous segments.This processing requires a connection propagation constraint, which was determined by experience in [40].Some short segments failed to be connected and were erroneously eliminated, consequently underestimating thin rivers and network continuity.By combining Gabor filtering and path opening, these difficulties are substantially overcome.
Detecting rivers from remotely sensed imagery is an initial yet critical step in river studies.The low spectral contrast between thin rivers and image background as well as the varied intensities along river courses pose huge challenges for automatic river detection.This paper proposes an automatic method to enhance and detect rivers from remotely sensed imagery.The underlying idea behind the proposed method is to integrate the cross-sectional (through Gabor filtering) and longitudinal (through path opening) river information in images to accomplish continuous river delineation.The results of the experiments show the advantages of the proposed method in detecting complete river networks.

Figure 1 .
Figure 1.(a) Landsat 8 OLI panchromatic image of example river networks distributed in the Yukon River Basin, northwestern North America, and (b) its corresponding three-dimensional (3D) view, in which pixel values in the original image are inversed for comparison convenience.The (0, 0) point in the 3D view corresponds to the top-left pixel of the image.

Figure 2 .
Figure 2. Proposed process steps: (a) Original Landsat 8 OLI panchromatic river image, (b) pre-processed river image, (c) image after Gabor filtering, and (d) image after path opening.The thickness of the Gabor filter, w, was set to 2 and the length threshold, Lmin, in the path opening was set to 40 pixels.

Figure 3 .
Figure 3. Three-dimensional view of the images in Figure 2: (a) the original image, (b) the pre-processed image, (c) the image after Gabor filtering, and (d) the image after path opening.

Figure 4 .
Figure 4. Gabor filter composition: (a) Real Gabor filter kernel oriented at an angle θ π/2; w is set to 20 for better display, and (b) the corresponding Gaussian filter kernel.

Figure 5 .
Figure 5. Path-opening results and their corresponding gray histograms under different Lmin settings.n is the total number of gray levels in each image.The intensities of all the images are rescaled to [0, 255] and showed in a blue-red color ramp for better display and comparison.Blue indicates lower intensities.

Figure 6 .
Figure 6.Three-dimensional views of the images in Figure 5 after path opening processing with different Lmin settings.

Figure 7 .
Figure 7. Classified rivers using a global threshold for: (a) original image, (b) pre-processed image, (c) image after Gabor filtering, and (d) image after path opening.

Figure 8 .
Figure 8. Histogram of the image after Gabor filtering and the path opening process and the corresponding river detection results under different global thresholds.The red circles show the thin rivers considered as potential rivers whose gray values are in the range of mean and mean + 1 std of the image.The value of mean + 0.5 std is considered as an optimal threshold to discern thin rivers from the image background.

Figure 9 .
Figure 9. Study area map.The Yukon River Basin and the Greenland Ice Sheet are selected as study areas.Both areas are located in the Arctic.River networks distributed in these two areas are well developed and thus ideal for river detection experiments.

Figure 10 .
Figure 10.Landsat 8 OLI panchromatic images for the study area of (a) the Yukon Basin and (b) the Greenland Ice Sheet ablation zone.The river detection results using the proposed method are shown in (c) and (d).The six small sites selected to test the performance of the proposed method are indicated by the white rectangles.

Figure 11 .
Figure 11.River detection results for Sites 1-3 in the Yukon image (rows 1-3) and Sites 4-6 in the Greenland Ice Sheet image (rows 4-6).The first column shows the site number; the second shows the original Landsat 8 OLI panchromatic images; the third to sixth columns show the river detection results when one global threshold (mean + 0.5 std) is applied to the original panchromatic image, the PPO-processed image, the Gabor-filtered image, and the Gabor + PPO processed (the proposed method) image, respectively; and the last column shows the corresponding ground-truth river maps.The red circles indicate the places where pixels were erroneously detected as rivers by the proposed method, while the yellow circles indicate the places where actual rivers were not correctly detected.

Figure 12 .
Figure 12.ROC graphs (TPRs against FPRs) of the river detection results when different global thresholds are applied to the original panchromatic image, the PPO-processed image, the Gabor-filtered image, and the Gabor + PPO processed (the proposed method) image.The result for Sites 1-6 is depicted by the curves in (a)-(f), respectively.

Figure 13 .
Figure 13.River network labeling results for (a) the Gabor filtered Yukon image; (b) the Gabor filtered GrIS image; (c) the Gabor-PPO processed Yukon image; and (d) the Gabor-PPO processed GrIS image.The river networks obtained from the Gabor filtered images are very fragmented, while those from the Gabor-PPO processed images are much more continuous, indicating that the proposed method can preserve the river network continuity very well.

Figure 14 .
Figure 14.The network and river pixel counts under different size thresholds.(a) The network count and (b) the corresponding river pixel count for the classified Yukon images.(c) The network count and (d) the river pixel count for the classified GrIS image.

Table 2 .
Performance evaluation of the classification images when one global threshold (value: mean + 0.5 std) is applied to the original panchromatic image, the PPO processed image, the Gabor filtered image, and the Gabor + PPO processed (the proposed method) image, compared with the manual digitization for Sites 4-6 in the Greenland Ice Sheet image.

Table 3 .
Summary statistics of the network count, the river pixel count, and the average pixel count per river network for the Gabor filtered and the Gabor-PPO processed images in the two study areas.