Rural Road Extraction from High-Resolution Remote Sensing Images Based on Geometric Feature Inference

Road information as a type of basic geographic information is very important for services such as city planning and traffic navigation, as such there is an urgent need for updating road information in a timely manner. Scholars have proposed various methods of extracting roads from remote sensing images, but most of them are not applicable to rural roads with diverse materials, large curvature changes, and a severe shelter problem. In view of these problems, we propose a road extraction method based on geometric feature inference. In this method, we make full use of the linear characteristics of roads, and construct a geometric knowledge base of rural roads using information on selected sample road segments. Based on the knowledge base, we identify the parallel line pairs in images, and further conduct grouping and connection instructed by knowledge reasoning, and finally obtain complete rural roads. The case study in Xiangtan City of China’s Hunan Province validates the performance of the proposed method.


Introduction
Road information is important basic geographic information that is directly related to services such as city planning, traffic navigation, and disaster evacuation [1,2].With the rapid development of the economy, construction and alteration projects of roads are becoming increasingly frequent, especially in developing countries, so there is an urgent need for methods for rapid updating of road information [3].The emergence of aerial and satellite remote sensing technology promotes the automatic or semi-automatic extraction of roads from remote sensing images to be a research focus [4,5].
In view of images with different resolutions, different data sources or roads of different grades, scholars have proposed a variety of road extraction methods or algorithms.They can be divided into four general categories: (1) methods based on spectral and geometric features.Because the surface material of roads is characterized by typical spectral reflectance characteristics, researchers can distinguish roads from non-road regions.For geometric features, a road appears as an image object with an obvious linear shape or two boundary lines in the remote sensing images.Some scholars extract road regions in the images according to geometric features of roads (in the form of lines or ribbons) and spectral features (spectrum of typical surface material) [6][7][8][9][10]; (2) methods based on multi-scale features.Road regions are identified according to the spectral and geometric characteristics of roads on images at different scales.Heipke first proposed a multi-resolution approach for automatic road extraction from aerial images.He used two different resolutions of the same aerial image: extract bright lines from the coarse-resolution image while extract road edges at high-resolution image [11].The edges were grouped into parallel line pairs, and then rules were used to combine results at two resolutions.
Mayer et al. [12], Mayer and Steger [13] supplemented the multi-resolution method using the concept of scale abstraction, and discussed the influence of the analysis scale on the image characteristics of roads.Since then, Hinz et al. [14], Huang and Zhang [15] have extended the multi-resolution extraction methods; (3) the methods based on context.The context in the process of road extraction refers to the background characteristics associated with road objects, which provides information on the relationship between roads and other objects [16,17], and compensates for information that cannot be extracted from roads.For example, Baumgartner [18] extracted road objects based on context related to roads (i.e., roads and houses, roads and vehicles).Wang et al. [19] established a semantic representation model of roads from the perspective of algebraic linguistics, and extracted urban road networks using topologic and spatial relations.In view of disaster assessment, Ma et al. [20] proposed a method to extract intact road regions from post-disaster remote sensing images under the guidance of prior road vectors; (4) the integration of the above two or more methods.Due to the complexity of imaging process, relying solely on road features or contextual features is unable to achieve satisfactory results in some cases.Therefore, it is sometimes necessary to take full advantage of geometrical features, spectral features, elevation feature, various contextual information as well as multi-scale image features to implement road extraction, such as the research conducted by [21][22][23].
Although the above studies have accumulated a lot of road extraction methods based on remote sensing images, most of them have been designed for high-grade or urban roads, while few are applicable to rural roads.According to statistics from the Ministry of Transport of China, the total mileage of roads in China is about 4.7 million kilometers by the end of 2016, and in total, the rural road mileage is 3.96 million kilometers, accounting for about 84% [24].Therefore, obtaining the accurate distribution and changes of rural roads is useful for government agencies to set up policies related to roads.The following problems occur when the above methods are applied to rural roads: First, the pavement of rural roads is made of various kinds of materials: asphalt, cement, gravel, stone, etc., so rural roads have different spectral reflection characteristics.Especially on high-resolution images without abundant spectral information, methods based mainly on spectral features do not perform well.Secondly, rural roads are usually not wide, and some road segments can be completely sheltered by buildings and trees along roads, so the shelter problem in rural areas is more severe than in urban areas.Third, rural roads have larger curvature than urban roads.It is difficult for existing methods to extract complete roads.
In order to solve the abovementioned challenges, this paper proposes a rural road extraction method based on geometric feature inference.In this method, a geometric knowledge base for rural roads is constructed from sample road segments, and parallel line pairs conforming to the characteristics of road edges are further identified.Then the grouping and linking are implemented on parallel line pairs by knowledge reasoning to extract complete rural roads.
The next section presents the workflow of the methodology and detailed algorithms.Section 3 introduces the experimental dataset and demonstrates extraction results in different steps.Section 4 discusses the merits and demerits of the proposed method compared with other approaches, and points out some limitations.Section 5 summarizes the work.

Methodology
In high-resolution remote sensing images, rural roads have typical image characteristics of roads as other types of roads such as urban roads and highways, but also have unique characteristics including material diversity, large curvature, and a severe shelter problem.Therefore, for automatically identifying rural roads from high spatial resolution images, the geometric features especially linear features should be major consideration factors while the spectral features are secondary factors.Based on this idea, this paper constructs a geometric feature knowledge base of rural roads, and uses it to guide the process of road extraction.The overall workflow of road extraction is shown in Figure 1, mainly divided into six steps: (1) knowledge base construction of geometric features of rural roads.In the area to be processed, the sample road segments are selected.The typical geometric features are calculated by related algorithms, and on the basis of them, a knowledge base is constructed; (2) image feature enhancement.The Gaussian filter is used to smooth the original image to remove random noises and enhance road-related features; (3) edge detection.The edge detection algorithm is used to detect high-gradient regions in the image, and the edge image is generated; (4) candidate road segment identification.The Hough Transform is used to detect all the lines existing in the edge image, and the parallel lines that conform to the geometric features of rural roads are identified as candidate segments; (5) adjacent segment linking.Based on geometric knowledge base, the candidate road segments that are adjacent are linked together; (6) missing road segment inference.For the missing road segments caused by image quality, occlusion or other reasons, the Ribbon Snake method is adopted for inferring the missing parts, and ultimately the complete road vectors are extracted.The edge detection algorithm is used to detect high-gradient regions in the image, and the edge image is generated; (4) candidate road segment identification.The Hough Transform is used to detect all the lines existing in the edge image, and the parallel lines that conform to the geometric features of rural roads are identified as candidate segments; (5) adjacent segment linking.Based on geometric knowledge base, the candidate road segments that are adjacent are linked together; (6) missing road segment inference.For the missing road segments caused by image quality, occlusion or other reasons, the Ribbon Snake method is adopted for inferring the missing parts, and ultimately the complete road vectors are extracted.

Geometric Knowledge Base Construction of Rural Roads
In order to obtain typical geometric features of rural roads as a guide of extracting roads from images, we first select representative sample road segments.The samples are visually interpreted based on remote sensing images and stored as vector format.It is better to select samples in the regions with different landforms, which would increase the generalizability of samples.For each type of material such as gravel, cement, asphalt and other materials, a certain number of samples should be chosen.The indicators that can reflect the geometric and spectral features of road segments are calculated including index values, value range and distribution.According to the extent of a sample road, the average gray values of each band within road regions are calculated, and the gradient value range at the edge of roads is also obtained.In addition, this paper takes the normal distribution as the theoretical distribution of the width of rural roads.In other words, the width of rural roads is expected to be around a value, and the width of most roads is close to the mean value.The distribution is represented by the function

Geometric Knowledge Base Construction of Rural Roads
In order to obtain typical geometric features of rural roads as a guide of extracting roads from images, we first select representative sample road segments.The samples are visually interpreted based on remote sensing images and stored as vector format.It is better to select samples in the regions with different landforms, which would increase the generalizability of samples.For each type of material such as gravel, cement, asphalt and other materials, a certain number of samples should be chosen.The indicators that can reflect the geometric and spectral features of road segments are calculated including index values, value range and distribution.According to the extent of a sample road, the average gray values of each band within road regions are calculated, and the gradient value range at the edge of roads is also obtained.In addition, this paper takes the normal distribution as the theoretical distribution of the width of rural roads.In other words, the width of rural roads is expected to be around a value, and the width of most roads is close to the mean value.The distribution is represented by the function pro(w) = dis(w).Based on the sample road segments, the maximum likelihood estimation method is used to calculate the two parameters (standard deviation σ and mathematical expectation µ) of road width distribution.After the distribution is obtained, the Kolmogorov-Smirnov test is used to verify if the observed width values accord with the fitted distribution model.If the K-S statistic D n is smaller than test statistic D α (α = 0.10), we think the quality of samples are satisfactory, otherwise, more samples are needed.

Image Feature Enhancement
In order to improve the detection of road edges, an interactive histogram stretching tool [25] can be used to enhance the gradient between roads and surrounding objects before edge detection.This is not a necessary step, but it is useful for an experienced practitioner.In addition, the Gaussian filter [26] is used for smoothing original image to remove random noises and enhance road features.The variance and neighborhood size of the Gaussian filter function can be set in advance, and then the Gaussian filter and the original image are convoluted to obtain a smoothed image.The Gaussian filter is calculated as follows: where σ is filter coefficient, and the larger it is, the stronger the ability of filtering noise is.While the Gaussian filter is removing noise, it would also blur the edges of roads, so it is necessary to select appropriate filter coefficients.

Edge Detection
Road is a typical striped geographic entity, and shows a strong edge feature on remote sensing images.As such, it is common to use edge information to identify possible road locations, and the detection of edges in images is a prerequisite for other analysis.We choose the edge detection method with embedded confidence [27] to generate edge images.

Line Detection
The Hough Transform is a commonly used straight line detection method for images [28,29].The main idea is to use mathematical transformation to convert a straight line in the image space into a point in the parameter space, and obtain the linear equation by detecting the position of points.As shown in Figure 2a, in the image space XOY, the straight line equations passing through the points (x i , y i ) and (x j , y j ) are expressed as y i = ax i + b and y j = ax j + b respectively, where a is the slope, and b is the intercept.The parameters a and b are treated as variables, to construct the parameter space AOB, and then the above-mentioned line equations are converted to b = −ax i + y i and b = −ax j + y j , that is the two lines passing through (a, b) in the parameter space, as shown in Figure 2b.The Hough Transform detects lines based on the global features in images, so it is robust in detecting the discontinuous boundaries of regions caused by noise disturbance.

Image Feature Enhancement
In order to improve the detection of road edges, an interactive histogram stretching tool [25] can be used to enhance the gradient between roads and surrounding objects before edge detection.This is not a necessary step, but it is useful for an experienced practitioner.In addition, the Gaussian filter [26] is used for smoothing original image to remove random noises and enhance road features.The variance and neighborhood size of the Gaussian filter function can be set in advance, and then the Gaussian filter and the original image are convoluted to obtain a smoothed image.The Gaussian filter is calculated as follows: where σ is filter coefficient, and the larger it is, the stronger the ability of filtering noise is.While the Gaussian filter is removing noise, it would also blur the edges of roads, so it is necessary to select appropriate filter coefficients.

Edge Detection
Road is a typical striped geographic entity, and shows a strong edge feature on remote sensing images.As such, it is common to use edge information to identify possible road locations, and the detection of edges in images is a prerequisite for other analysis.We choose the edge detection method with embedded confidence [27] to generate edge images.

Line Detection
The Hough Transform is a commonly used straight line detection method for images [28,29].The main idea is to use mathematical transformation to convert a straight line in the image space into a point in the parameter space, and obtain the linear equation by detecting the position of points.As shown in Figure 2a, in the image space XOY, the straight line equations passing through the points ( , ) i i x y and ( , )

Parallel Line Identification
According to the geometric features of rural roads, the two edges of roads usually have the following characteristics [7,30]: two edges are approximately parallel; the distance between two parallel lines is the width of the road; the edges consist of straight lines and curves, but straight lines

Parallel Line Identification
According to the geometric features of rural roads, the two edges of roads usually have the following characteristics [7,30]: two edges are approximately parallel; the distance between two parallel lines is the width of the road; the edges consist of straight lines and curves, but straight lines account for a major proportion.Based on the above rules, we calculate the index that any two straight lines conform to the characteristics of road edges, and then select the parallel line pair with a maximum grouping index as a candidate road segment.We use pro(∆θ), pro(w), pro(len) to represent parallel index, width index and extension index, respectively.The grouping index for two straight lines can be expressed as: pro = pro(∆θ) × pro(w) × pro(len). ( After Hough Transform, the three parameters of a straight line can be obtained: radius ρ, polar angle θ and cumulative edge point number A(ρ, θ), and the index values are calculated by constructing the following models.
(1) Parallel index The closer to zero the direction difference ∆θ is, the closer to be parallel the two lines are, and the higher the parallel index is.According to this criterion, the parallel index of two straight lines is defined as: when ∆θ = 0, the index goes up to the highest value of 1, indicating that the two lines are accurately parallel; when ∆θ = ±π/2, the index is 0, indicating that the two lines are perpendicular to each other.The distribution of parallel probability is shown in Figure 3.
account for a major proportion.Based on the above rules, we calculate the index that any two straight lines conform to the characteristics of road edges, and then select the parallel line pair with a maximum grouping index as a candidate road segment.We use ( ) pro len to represent parallel index, width index and extension index, respectively.The grouping index for two straight lines can be expressed as: After Hough Transform, the three parameters of a straight line can be obtained: radius ρ , polar angle θ and cumulative edge point number ( , ) A ρ θ , and the index values are calculated by constructing the following models.
(1) Parallel index The closer to zero the direction difference θ Δ is, the closer to be parallel the two lines are, and the higher the parallel index is.According to this criterion, the parallel index of two straight lines is defined as: (2) Width index For two straight lines, the distance has a stable value only when they are approximately parallel, otherwise the distance is quite different at different positions.Therefore, we define the distance formula as follows: where 1 ρ and 2 ρ are the radius of two straight lines respectively, m θ is the maximum angle difference for calculating the distance between lines.After the width of road is obtained, the index distribution function of the width value

( ) ( ) pro w dis w =
in the geometric knowledge base is used to calculate the width index. (

3) Extension probability
The straight lines detected by the Hough Transform only reflect global characteristics of the edges, but cannot reflect the actual number of edge points on the straight line, or the reliability of the straight line.We take the average number of edge points of sample road segments ( ref length ) as the reference threshold, and calculate the extension index: (2) Width index For two straight lines, the distance has a stable value only when they are approximately parallel, otherwise the distance is quite different at different positions.Therefore, we define the distance formula as follows: where ρ 1 and ρ 2 are the radius of two straight lines respectively, θ m is the maximum angle difference for calculating the distance between lines.After the width of road is obtained, the index distribution function of the width value pro(w) = dis(w) in the geometric knowledge base is used to calculate the width index. (

3) Extension probability
The straight lines detected by the Hough Transform only reflect global characteristics of the edges, but cannot reflect the actual number of edge points on the straight line, or the reliability of the straight line.We take the average number of edge points of sample road segments (length re f ) as the reference threshold, and calculate the extension index: The distribution of the extension index is shown in Figure 4.
The distribution of the extension index is shown in Figure 4.The joint extension index of two lines is calculated as follows: For a set of straight lines detected by the Hough Transform, the index of any two parallel lines is obtained according to the grouping index formula, and the parallel line pair whose index is larger than a predefined threshold limit pro is selected as a candidate road segment.The joint extension index of two lines is calculated as follows: For a set of straight lines detected by the Hough Transform, the index of any two parallel lines is obtained according to the grouping index formula, and the parallel line pair whose index is larger than a predefined threshold pro limit is selected as a candidate road segment.pro limit is recommended to set to 0.3-0.4.The larger pro limit is, the higher the extraction accuracy is, but the lower the recall rate is. Figure 5 shows the detection process of a candidate road segment.Multiple straight lines are detected in different regions in Figure 5a, and candidate road segments are identified.The result is shown in Figure 5b.
The distribution of the extension index is shown in Figure 4.The joint extension index of two lines is calculated as follows: For a set of straight lines detected by the Hough Transform, the index of any two parallel lines is obtained according to the grouping index formula, and the parallel line pair whose index is larger than a predefined threshold limit pro is selected as a candidate road segment.Considering that roads show as smoothing ribbons on high-resolution images, connecting adjacent candidate segments mainly takes into account the following indicators: (1) Width similarity index ( ) pro w : represents the similarity between the width values of two candidate segments, is calculated as: where width and 2 width are the width of two candidate road segments.
(2) Direction consistence index ( ) pro θ : represents the degree of consistence of directions of two candidate segments, is calculated as: where 1 θ and 2 θ are the directions of center lines of two candidate segments.
(3) Adjacency index ( ) pro Dis : represents the degree of adjacency of the endpoints of two candidate segments.The distance is calculated as the average distance between the left and right endpoints, that is ( 1, 2) ( (LB ,LB ) (RB ,RB ))/2 Dis seg seg Dis Dis = + .On the basis of it, the adjacency index is calculated as: ( ) pro Dis e − ⋅ = .
Each index metric is given a weight according to their impact on connection, and the linking index of the two road segments can be calculated by the following formula: where i pro are three indexes, i R are their respective weights.The road segments would be linked when their linking index is larger than the predefined threshold lim it P .Considering that roads in rural areas have large curvature changes, we assign a relatively low weight to the direction consistence indicator, so the weights of the width similarity index, direction consistence index, and adjacency index are recommended to be set to 0.4, 0.2, and 0.4, respectively.lim it P is recommended to be set to around 0.55 based on a series of experiments on rural roads in China.

Missing Segment Inference
Some road segments might appear discontinuous on images due to the quality of image acquisition or the shelter of buildings and trees, which means candidate road segments cannot be detected in these areas.Thus, after candidate road segment extraction, a series of discontinuous parallel lines is generated.How to infer the missing parts and connect those discontinuous segments to form a complete road is a key procedure.Based on the spatial location and geometric characteristics of rural roads, a geometry inference model is established using the Ribbon Snake algorithm [30], which evaluates the connection index between two road segments and determines whether candidate segments should be connected.Considering that roads show as smoothing ribbons on high-resolution images, connecting adjacent candidate segments mainly takes into account the following indicators: (1) Width similarity index pro(w): represents the similarity between the width values of two candidate segments, is calculated as: where width 1 and width 2 are the width of two candidate road segments.(2) Direction consistence index pro(θ): represents the degree of consistence of directions of two candidate segments, is calculated as: where θ 1 and θ 2 are the directions of center lines of two candidate segments.(3) Adjacency index pro(Dis): represents the degree of adjacency of the endpoints of two candidate segments.The distance is calculated as the average distance between the left and right endpoints, that is Dis(seg1, seg2) = (Dis(LB 1 , LB 2 ) + Dis(RB 1 , RB 2 ))/2.On the basis of it, the adjacency index is calculated as: pro(Dis) = e −0.1•Dis(Seg1,Seg2) .
Each index metric is given a weight according to their impact on connection, and the linking index of the two road segments can be calculated by the following formula: where pro i are three indexes, R i are their respective weights.The road segments would be linked when their linking index is larger than the predefined threshold P limit .Considering that roads in rural areas have large curvature changes, we assign a relatively low weight to the direction consistence indicator, so the weights of the width similarity index, direction consistence index, and adjacency index are recommended to be set to 0.4, 0.2, and 0.4, respectively.P limit is recommended to be set to around 0.55 based on a series of experiments on rural roads in China.

Missing Segment Inference
Some road segments might appear discontinuous on images due to the quality of image acquisition or the shelter of buildings and trees, which means candidate road segments cannot be detected in these areas.Thus, after candidate road segment extraction, a series of discontinuous parallel lines is generated.How to infer the missing parts and connect those discontinuous segments to form a complete road is a key procedure.Based on the spatial location and geometric characteristics of rural roads, a geometry inference model is established using the Ribbon Snake algorithm [30], which evaluates the connection index between two road segments and determines whether candidate segments should be connected.
The Ribbon Snake algorithm is an extension of the traditional Snake algorithm, which adds the width variable into the traditional Snake algorithm.It is widely used in the research of road extraction.The mathematical formula of a ribbon snake is: → v (s, t) = (x(s, t), y(s, t), w(s, t)), (0 < s < 1), (10) where s is the ratio length that is proportional to the length of the ribbon, t is the current time, x(s), y(s) are point coordinates of the center line of the ribbon, and w(s) is the width of ribbon at the point x(s), y(s).For each center point, there exist two points v L (s) = (x(s), y(s), w(s)) and v R (s) = (x(s), y(s), w(s)) on the left and right boundary line.Taking the endpoints of the identified candidate road segments as initialized control points, we use the Ribbon Snake algorithm to infer the missing segments.In order to make sure the identified candidate road segments and inferred segments together form complete and smoothed roads, the following constraint conditions are defined.

Continuity Constraint of Road Edge
In order to avoid the problem that the traditional Snake model may cause the accumulation of control points during the convergence process, an energy part is defined to ensure the control points evenly distribute along road edges: where v i and v i−1 are two adjacent control points, d is the average distance among control points on road edges, which is calculated as: where v i and v i−1 are two adjacent control points, N is the number of control points.The continuity energy of a road with two boundary lines is expressed by the average energy of the two edges, that is

Curvature Constraint of Road Edge
We adjust the curvature of curves through bending the energy.For a single curve, the discrete bending energy is calculated by the second derivative: The result of minimizing the energy function is to straighten the curve.As shown in Figure 7, when the control point V n (i, j) is iterated to a minimum value, and the iteration extent is assumed to be a 3 × 3 neighborhood, the point V n (i, j) would move to (i + 1, j), and the polyline would become a straight line.For a Ribbon road, the total bending energy at the current position is defined as the average energy of two lines, that is: ISPRS Int.J. Geo-Inf.2017, 6, 314 8 of 23 The Ribbon Snake algorithm is an extension of the traditional Snake algorithm, which adds the width variable into the traditional Snake algorithm.It is widely used in the research of road extraction.The mathematical formula of a ribbon snake is: ( , ) ( ( , ), ( , ), ( , )) where s is the ratio length that is proportional to the length of the ribbon, t is the current time, ( ), ( ) x s y s are point coordinates of the center line of the ribbon, and ( ) w s is the width of ribbon at the point ( ), ( ) x s y s .For each center point, there exist two points v ( ) ( ( ), ( ), ( )) L s x s y s w s = and ( ) ( ( ), ( ), ( )) R v s x s y s w s = on the left and right boundary line.Taking the endpoints of the identified candidate road segments as initialized control points, we use the Ribbon Snake algorithm to infer the missing segments.In order to make sure the identified candidate road segments and inferred segments together form complete and smoothed roads, the following constraint conditions are defined.

Continuity Constraint of Road Edge
In order to avoid the problem that the traditional Snake model may cause the accumulation of control points during the convergence process, an energy part is defined to ensure the control points evenly distribute along road edges: where i v and 1 i v − are two adjacent control points, d is the average distance among control points on road edges, which is calculated as: where i v and 1 i v − are two adjacent control points, N is the number of control points.The continuity energy of a road with two boundary lines is expressed by the average energy of the two edges, that is

Curvature Constraint of Road Edge
We adjust the curvature of curves through bending the energy.For a single curve, the discrete bending energy is calculated by the second derivative: . The result of minimizing the energy function is to straighten the curve.As shown in Figure 7, when the control point n V i j is iterated to a minimum value, and the iteration extent is assumed to be a 3 × 3 neighborhood, the point ( , ) , and the polyline would become a straight line.For a Ribbon road, the total bending energy at the current position is defined as the average energy of two lines, that is:

Width Constraint of Road Edge
The width of a road is assumed to be the same at most positions, and the direction of the road usually changes in a continuous way.Thus, the two edges of a road should be parallel at a local scale.The distance between the control point pairs on two edges should be within a certain range.According to this rule, the width constraint energy of Snake is defined as: where v L i and v R i are the left and right control points on the two edge lines along the clockwise direction, w is the expected road width, and i is the subscript of control point.

Image Feature Constraint
A road usually corresponds to the long edge feature in the image, and Snake draws the curve to a significant image feature through the external energy.We define the length energy E length (v(s i )) that reflects the continuity of the edge, and meanwhile, enlarge its impact according to distance to form the distance potential E d (v(s i )), and finally composes them to form the image feature constraint E ext (v(s i )) for each Snake point: where η 1 and η 2 are weights of length energy and distance energy for edges, E length (s i ) is the negative image of length energy of edges, E d (s i ) is the distance potential function, which can be calculated by: In summary, the total energy function is as follows: where α, β, γ, λ 1 , λ 2 are the weight coefficients for each type of energy.
Based on the equations of the four constraints, the smaller the total energy is, the more complete and smoothed the inferred roads are.Therefore, the task of inferring missing segments transforms to the task of minimizing the total energy function.In other words, the ultimate goal is to find a control point sequence v i (i = 0, 1, • • • , n − 1) in the image space that can minimize the total energy of the Snake curve: There are mainly three convergence algorithms for Snake model: the variational method [31], the dynamic programming [32] and the greedy algorithm [33].Considering that the greedy algorithm owns all advantages of dynamic programming, is faster in solving process than the variational method, and requires smaller storage space, we use the greedy algorithm for optimization.The greedy algorithm [33] follows the problem solving heuristic of making the locally optimal choice at each stage; based on the identified candidate road segments, it can approximate a global optimal solution in a reasonable time.It is called greedy because it provides an immediate output, but does not consider the problem as a whole.

Experiment Data
We choose a rural area in Xiangtan City, Hunan Province, China as the study area, with a geographical range of 27 • 48 28 -27 • 50 7 N, 112 • 22 28 -112 • 24 35 E. This area is a hilly area with a large number of distributed villages, connected by roads of various grades.A high-resolution aerial remote sensing image in this area with a spatial resolution of 0.5 m is used as the study dataset, and the image size is 7000 × 6000 with three bands: a red band, a green band, and a blue band.An image is shown in Figure 8.In order to verify the detection accuracy of the proposed method, we also produce road vectors of the study area through artificial interpretation, and take them as the ground truth data.

Experiment Data
We choose a rural area in Xiangtan City, Hunan Province, China as the study area, with a geographical range of 27°48′28″-27°50′7″ N, 112°22′28″-112°24′35″ E. This area is a hilly area with a large number of distributed villages, connected by roads of various grades.A high-resolution aerial remote sensing image in this area with a spatial resolution of 0.5 m is used as the study dataset, and the image size is 7000 × 6000 with three bands: a red band, a green band, and a blue band.An image is shown in Figure 8.In order to verify the detection accuracy of the proposed method, we also produce road vectors of the study area through artificial interpretation, and take them as the ground truth data.Twenty-four sample road segments are selected to construct a geometric knowledge base of rural roads, as shown in Figure 9a.Sample roads are manually interpreted from the high-resolution image including all levels of roads, and the road edges are stored in the form of a vector.The length of sample road segments is between 50 and 200 m. Figure 9b-d  Twenty-four sample road segments are selected to construct a geometric knowledge base of rural roads, as shown in Figure 9a.Sample roads are manually interpreted from the high-resolution image including all levels of roads, and the road edges are stored in the form of a vector.The length of sample road segments is between 50 and 200 m. Figure 9b-d show the superimposed display effect of sample road segments and the high-resolution remote sensing image in three local regions.

Experiment Data
We choose a rural area in Xiangtan City, Hunan Province, China as the study area, with a geographical range of 27°48′28″-27°50′7″ N, 112°22′28″-112°24′35″ E. This area is a hilly area with a large number of distributed villages, connected by roads of various grades.A high-resolution aerial remote sensing image in this area with a spatial resolution of 0.5 m is used as the study dataset, and the image size is 7000 × 6000 with three bands: a red band, a green band, and a blue band.An image is shown in Figure 8.In order to verify the detection accuracy of the proposed method, we also produce road vectors of the study area through artificial interpretation, and take them as the ground truth data.Twenty-four sample road segments are selected to construct a geometric knowledge base of rural roads, as shown in Figure 9a.Sample roads are manually interpreted from the high-resolution image including all levels of roads, and the road edges are stored in the form of a vector.The length of sample road segments is between 50 and 200 m. Figure 9b-d

Road Extraction Results
The entire experiment is implemented on a machine with dual 3.2 GHz CPU, 8 GB RAM, and 64-bit operating system.The whole process, except building a geometric knowledge base (since it relies on manual interpretation of road segments), takes 72 s.The remote sensing image is processed by the enhancement method described in Section 2.2 and the edge detection with embedded confidence in Section 2.3.The edge detection result of the entire study area is shown by Figure 10a.Figure 10b-e show the results in four small areas in detail.The left side shows the original remote sensing image, while the right side shows the edge detection result.It can be seen that the algorithm accurately detects all obvious gradient changes in the image.The two edges of roads are accurately identified whenever the curvature is small or large, or there are sudden turns.Especially in Figure 10e, where there are two types of roads, a cement road and a gravel road, the edge detection results for both are good.

Road Extraction Results
The entire experiment is implemented on a machine with dual 3.2 GHz CPU, 8 GB RAM, and 64-bit operating system.The whole process, except building a geometric knowledge base (since it relies on manual interpretation of road segments), takes 72 s.The remote sensing image is processed by the enhancement method described in Section 2.2 and the edge detection with embedded confidence in Section 2.3.The edge detection result of the entire study area is shown by Figure 10a.Figure 10b-e show the results in four small areas in detail.The left side shows the original remote sensing image, while the right side shows the edge detection result.It can be seen that the algorithm accurately detects all obvious gradient changes in the image.The two edges of roads are accurately identified whenever the curvature is small or large, or there are sudden turns.Especially in Figure 10e, where there are two types of roads, a cement road and a gravel road, the edge detection results for both are good.

Road Extraction Results
The entire experiment is implemented on a machine with dual 3.2 GHz CPU, 8 GB RAM, and 64-bit operating system.The whole process, except building a geometric knowledge base (since it relies on manual interpretation of road segments), takes 72 s.The remote sensing image is processed by the enhancement method described in Section 2.2 and the edge detection with embedded confidence in Section 2.3.The edge detection result of the entire study area is shown by Figure 10a.Figure 10b-e show the results in four small areas in detail.The left side shows the original remote sensing image, while the right side shows the edge detection result.It can be seen that the algorithm accurately detects all obvious gradient changes in the image.The two edges of roads are accurately identified whenever the curvature is small or large, or there are sudden turns.Especially in Figure 10e, where there are two types of roads, a cement road and a gravel road, the edge detection results for both are good.For parallel line identification, the maximum angle difference θ m is set to 30 • , and the grouping probability threshold pro limit is set to 0.35.For road segment linking, the weights of width similarity probability, direction consistence probability and adjacency probability are set to 0.4, 0.2 and 0.4 respectively, and the linking probability threshold P limit is set to 0.55. Figure 11a-e show the entire road extraction process based on geometric feature inference.Figure 11a presents the edge detection result, Figure 11b shows the result of candidate road segment identification, Figure 11c is the linking result of the candidate segments, Figure 11d is the inference result of missing road segments, and Figure 11e is the overlay of remote sensing images and extracted roads.It can be found by comparing automatic extraction results with images that the results are overall good.Both cement roads and gravel roads are extracted, and even road segments with a large curvature have been accurately extracted.However, there are still a small number of wrong extractions.A regular striped farmland next to the road in the left part of image has been misinterpreted as roads because it has multiple parallel edges.In addition, the house adjacent to the road in the lower right corner is mistaken for a road, and this is probably because the shape of the house is similar to a road.
ISPRS Int.J. Geo-Inf.2017, 6, 314 13 of 23 For parallel line identification, the maximum angle difference m θ is set to 30°, and the grouping probability threshold limit pro is set to 0.35.For road segment linking, the weights of width similarity probability, direction consistence probability and adjacency probability are set to 0.4, 0.2 and 0.4 respectively, and the linking probability threshold lim P it is set to 0.55. Figure 11a-e show the entire road extraction process based on geometric feature inference.Figure 11a presents the edge detection result, Figure 11b shows the result of candidate road segment identification, Figure 11c is the linking result of the candidate segments, Figure 11d is the inference result of missing road segments, and Figure 11e is the overlay of remote sensing images and extracted roads.It can be found by comparing automatic extraction results with images that the results are overall good.Both cement roads and gravel roads are extracted, and even road segments with a large curvature have been accurately extracted.However, there are still a small number of wrong extractions.A regular striped farmland next to the road in the left part of image has been misinterpreted as roads because it has multiple parallel edges.In addition, the house adjacent to the road in the lower right corner is mistaken for a road, and this is probably because the shape of the house is similar to a road.
where TP is true positives, denoting the total length of correctly extracted roads; FN is false negatives, denoting the total length of undetected roads; FP is false positives, denoting the total length of incorrectly extracted roads.For the road extraction result, TP, FN, and FP are 74.86 km, 9.82 km, and 15.77 km, respectively.Therefore, the recall rate and accuracy rate are 88.4% and 82.6%, respectively [34].
A reason why the accuracy rate is comparably lower than recall rate is that most roads can be extracted by their geometric features, but some entities with geometric features similar to roads are mistaken for roads.The incorrectly extracted roads are mainly those roads with a large portion of shelter by dense trees or buildings, which leads to inaccurate inferences or linking of road segments.In addition, the extraction of high-grade roads with obvious image features is more accurate than low-grade roads with vague features.In general, this method can perform well in extracting rural roads with diverse materials and large curvature changes.
where TP is true positives, denoting the total length of correctly extracted roads; FN is false negatives, denoting the total length of undetected roads; FP is false positives, denoting the total length of incorrectly extracted roads.For the road extraction result, TP, FN, and FP are 74.86 km, 9.82 km, and 15.77 km, respectively.Therefore, the recall rate and accuracy rate are 88.4% and 82.6%, respectively [34].
A reason why the accuracy rate is comparably lower than recall rate is that most roads can be extracted by their geometric features, but some entities with geometric features similar to roads are mistaken for roads.The incorrectly extracted roads are mainly those roads with a large portion of shelter by dense trees or buildings, which leads to inaccurate inferences or linking of road segments.In addition, the extraction of high-grade roads with obvious image features is more accurate than low-grade roads with vague features.In general, this method can perform well in extracting rural roads with diverse materials and large curvature changes.In order to show the performance of the proposed method in more details, we select a flat area (Region A) and a hilly area (Region B), and compare the road extraction results and the corresponding ground truth in these two regions as shown in Figure 12.The two regions correspond to Figure 13a,b respectively.The figures on the left show the superposition of the aerial remote sensing image and the ground truth data (yellow line) that are derived from artificial interpretation, while the figures on the right show the extraction results.Region A is mainly flat farmland, mixed with some small hills and woods.It can be seen that the extraction results and the ground truth are basically the same.Any small differences are mainly due to the farmland having wide ridges or building edges.In addition, the roads in the top middle corner of the image are not detected because of a long row of dense trees by the roadside.Region B is covered by small hills and woods, and the extraction results are fairly good.There are also a certain number of discrepancies.For example, a few roads in the middle of the image are not detected because they are dirt roads without obvious boundaries and regular shape.Another problem appears mainly at the end of village roads, where roads connect to a housing community.The shape of the roads in these areas is relatively irregular and complicated, resulting in a small amount of wrong extractions or omissions, but it has little impact on the overall length and distribution of roads.In order to show the performance of the proposed method in more details, we select a flat area (Region A) and a hilly area (Region B), and compare the road extraction results and the corresponding ground truth in these two regions as shown in Figure 12.The two regions correspond to Figure 13a,b respectively.The figures on the left show the superposition of the aerial remote sensing image and the ground truth data (yellow line) that are derived from artificial interpretation, while the figures on the right show the extraction results.Region A is mainly flat farmland, mixed with some small hills and woods.It can be seen that the extraction results and the ground truth are basically the same.Any small differences are mainly due to the farmland having wide ridges or building edges.In addition, the roads in the top middle corner of the image are not detected because of a long row of dense trees by the roadside.Region B is covered by small hills and woods, and the extraction results are fairly good.There are also a certain number of discrepancies.For example, a few roads in the middle of the image are not detected because they are dirt roads without obvious boundaries and regular shape.Another problem appears mainly at the end of village roads, where roads connect to a housing community.The shape of the roads in these areas is relatively irregular and complicated, resulting in a small amount of wrong extractions or omissions, but it has little impact on the overall length and distribution of roads.In order to show the performance of the proposed method in more details, we select a flat area (Region A) and a hilly area (Region B), and compare the road extraction results and the corresponding ground truth in these two regions as shown in Figure 12.The two regions correspond to Figure 13a,b respectively.The figures on the left show the superposition of the aerial remote sensing image and the ground truth data (yellow line) that are derived from artificial interpretation, while the figures on the right show the extraction results.Region A is mainly flat farmland, mixed with some small hills and woods.It can be seen that the extraction results and the ground truth are basically the same.Any small differences are mainly due to the farmland having wide ridges or building edges.In addition, the roads in the top middle corner of the image are not detected because of a long row of dense trees by the roadside.Region B is covered by small hills and woods, and the extraction results are fairly good.There are also a certain number of discrepancies.For example, a few roads in the middle of the image are not detected because they are dirt roads without obvious boundaries and regular shape.Another problem appears mainly at the end of village roads, where roads connect to a housing community.The shape of the roads in these areas is relatively irregular and complicated, resulting in a small amount of wrong extractions or omissions, but it has little impact on the overall length and distribution of roads. (a)

Discussion
Many scholars such as Simler [35] and Wang et al. [36] proposed to extract roads from images based on object-oriented methods.However, the object-oriented methods do not perform well for rural roads due to the low classification accuracy caused by material diversity.We adopt the strategy of extracting roads through identifying parallel elements on images, which avoids the problem of object-oriented methods.Some researchers [37][38][39] also proposed extracting roads based on geometric features or mathematical morphology features that are not much influenced by material diversity.However, their methods cannot work well on a road with many shadows and occlusions in rural areas.Our method is able to infer the missing road segments by geometric feature inference, which is designed for serious occlusion problems.
The quality of edge detection result directly influences the accuracy of road detection since the proposed method is based mainly on geometric features.There are various kinds of edge detection methods such as Sobel, LOG, and Canny, and they are applied to process a test area shown in Figure 14a.The detection results of both Sobel and LOG operator are wide edges with low location accuracy, as shown in Figure 14b,c, so they are not suitable for road detection.The Canny operator is able to detect thin edges from images with high accuracy; however, it would generate significant background noise while detecting weak edges, as shown in Figure 14d.In order to solve this problem, Meer and Georgescu proposed a edge detection operator with embedded confidence by defining a confidence measure and integrating it into all steps in a gradient-based edge detection procedure [27].Therefore, considering the ability of their method in detecting weak road edges while eliminating background noises, we choose it for edge detection. (a)

Discussion
Many scholars such as Simler [35] and Wang et al. [36] proposed to extract roads from images based on object-oriented methods.However, the object-oriented methods do not perform well for rural roads due to the low classification accuracy caused by material diversity.We adopt the strategy of extracting roads through identifying parallel elements on images, which avoids the problem of object-oriented methods.Some researchers [37][38][39] also proposed extracting roads based on geometric features or mathematical morphology features that are not much influenced by material diversity.However, their methods cannot work well on a road with many shadows and occlusions in rural areas.Our method is able to infer the missing road segments by geometric feature inference, which is designed for serious occlusion problems.
The quality of edge detection result directly influences the accuracy of road detection since the proposed method is based mainly on geometric features.There are various kinds of edge detection methods such as Sobel, LOG, and Canny, and they are applied to process a test area shown in Figure 14a.The detection results of both Sobel and LOG operator are wide edges with low location accuracy, as shown in Figure 14b,c, so they are not suitable for road detection.The Canny operator is able to detect thin edges from images with high accuracy; however, it would generate significant background noise while detecting weak edges, as shown in Figure 14d.In order to solve this problem, Meer and Georgescu proposed a edge detection operator with embedded confidence by defining a confidence measure and integrating it into all steps in a gradient-based edge detection procedure [27].Therefore, considering the ability of their method in detecting weak road edges while eliminating background noises, we choose it for edge detection.

Discussion
Many scholars such as Simler [35] and Wang et al. [36] proposed to extract roads from images based on object-oriented methods.However, the object-oriented methods do not perform well for rural roads due to the low classification accuracy caused by material diversity.We adopt the strategy of extracting roads through identifying parallel elements on images, which avoids the problem of object-oriented methods.Some researchers [37][38][39] also proposed extracting roads based on geometric features or mathematical morphology features that are not much influenced by material diversity.However, their methods cannot work well on a road with many shadows and occlusions in rural areas.Our method is able to infer the missing road segments by geometric feature inference, which is designed for serious occlusion problems.
The quality of edge detection result directly influences the accuracy of road detection since the proposed method is based mainly on geometric features.There are various kinds of edge detection methods such as Sobel, LOG, and Canny, and they are applied to process a test area shown in Figure 14a.The detection results of both Sobel and LOG operator are wide edges with low location accuracy, as shown in Figure 14b,c, so they are not suitable for road detection.The Canny operator is able to detect thin edges from images with high accuracy; however, it would generate significant background noise while detecting weak edges, as shown in Figure 14d.In order to solve this problem, Meer and Georgescu proposed a edge detection operator with embedded confidence by defining a confidence measure and integrating it into all steps in a gradient-based edge detection procedure [27].Therefore, considering the ability of their method in detecting weak road edges while eliminating background noises, we choose it for edge detection.In order to compare the difference between the proposed method that is based on geometric features and other approaches in a quantitative way, we implemented the object-oriented classification method (OOC) [40] using eCognition and the pixel-based classification (PBC) method using Envi on the same remote sensing image of the study area.In the OOC method, the multi-resolution segmentation method [41] is first used to partition the entire image into image objects based on the common characteristics of pixels, with the scale parameter, color factor, and smoothness set to 25, 0.5, and 0.3, respectively.Then we selected 38 road samples and 95 non-road samples to train a classification model for extracting roads.Figure 15 shows the extraction results in which the red color regions are road segments.Generally, the OOC method also obtains the overall spatial distribution of roads in the study area.In the PBC method, the same training dataset as the OOC method is used to train a maximum likelihood classifier, and Figure 16 shows the extraction results.In order to compare the difference between the proposed method that is based on geometric features and other approaches in a quantitative way, we implemented the object-oriented classification method (OOC) [40] using eCognition and the pixel-based classification (PBC) method using Envi on the same remote sensing image of the study area.In the OOC method, the multi-resolution segmentation method [41] is first used to partition the entire image into image objects based on the common characteristics of pixels, with the scale parameter, color factor, and smoothness set to 25, 0.5, and 0.3, respectively.Then we selected 38 road samples and 95 non-road samples to train a classification model for extracting roads.Figure 15 shows the extraction results in which the red color regions are road segments.Generally, the OOC method also obtains the overall spatial distribution of roads in the study area.In the PBC method, the same training dataset as the OOC method is used to train a maximum likelihood classifier, and Figure 16 shows the extraction results.It can be seen that the PBC result is so fragmented that few extracted roads are complete.This is probably because the limited spectral information on a single pixel cannot reflect the characteristics of rural roads.Nevertheless, the OOC result is much better.Figure 17 shows the overlap of the extraction result of the proposed method based on geometric feature inference (GFI) and the result of the OOC method.The total length of the roads extracted by both methods is 51.26 km, the length of roads that have been detected by GFI but missed by OOC is 39.37 km, and the length of roads that have been detected by OOC but missed by GFI is 5.48 km.It can be seen that the PBC result is so fragmented that few extracted roads are complete.This is probably because the limited spectral information on a single pixel cannot reflect the characteristics of rural roads.Nevertheless, the OOC result is much better.Figure 17 shows the overlap of the extraction result of the proposed method based on geometric feature inference (GFI) and the result of the OOC method.The total length of the roads extracted by both methods is 51.26 km, the length of roads that have been detected by GFI but missed by OOC is 39.37 km, and the length of roads that have been detected by OOC but missed by GFI is 5.48 km.It can be seen that the PBC result is so fragmented that few extracted roads are complete.This is probably because the limited spectral information on a single pixel cannot reflect the characteristics of rural roads.Nevertheless, the OOC result is much better.Figure 17 shows the overlap of the extraction result of the proposed method based on geometric feature inference (GFI) and the result of the OOC method.The total length of the roads extracted by both methods is 51.26 km, the length of roads that have been detected by GFI but missed by OOC is 39.37 km, and the length of roads that have been detected by OOC but missed by GFI is 5.48 km.The merits and demerits of the GFI method are analyzed by comparing the result of the OOC method.We can learn the following aspects: (1) For the four major roads in the study area, both methods obtain results of relatively good quality, and they correlate well with each other.This is attributed to the clear spectral and shape features of roads: the four major roads show apparent boundary lines on the image, and their pavement materials are consistent in different segments.Therefore, the ability of the two methods in handling high-grade roads is similar.(2) The extraction result of GFI is generally more complete and accurate than that of OOC in terms of connectivity.In rural areas, buildings and trees along the road often cause discontinuities of road segments on remote sensing images.The OOC method extracts road regions based on the spectral and texture features of image objects, and therefore cannot identify those segments sheltered by buildings or trees.By contrast, the GFI method first identifies road segments with clear boundaries and then connects nearby segments by knowledge inference, so it overcomes the discontinuity problem. Figure 18 demonstrates the extraction results of both methods on a road with a large portion of shelter by trees.(3) By comparison, we can easily find that the result of OOC missed quite a few roads.It is probably because the pavement of rural roads is made of various kinds of materials: asphalt, cement, gravel, stone, etc.The material heterogeneity makes it hard to select appropriate samples and establish classification rules in the OOC method.On the other hand, the GFI method is good at handling roads of heterogeneous material since it is based mainly on geometric features.(4) Because the GFI method extracts road segments by identifying parallel segments that conform to road characteristics, it is sensitive to parallel elements in the remote sensing image and can lead to wrong extractions.Figure 19 shows the extraction results in a village.It can be seen from Figure 19b that a large farmland ridge is wrongly identified as a road by GFI because it has two obvious boundary lines, and its width is similar to a road's.However, the OOC method obtains the correct result since the farmland ridge has different spectral and texture features from a road.Spectral and texture information would be taken into account in our method in future to increase the extraction accuracy, especially for high-resolution images with more spectral bands.(5) Another aspect of the GFI method we need to pay attention to is that its extraction result is not very accurate when processing road segments with an irregular shape.The extraction result of The merits and demerits of the GFI method are analyzed by comparing the result of the OOC method.We can learn the following aspects: (1) For the four major roads in the study area, both methods obtain results of relatively good quality, and they correlate well with each other.This is attributed to the clear spectral and shape features of roads: the four major roads show apparent boundary lines on the image, and their pavement materials are consistent in different segments.Therefore, the ability of the two methods in handling high-grade roads is similar.(2) The extraction result of GFI is generally more complete and accurate than that of OOC in terms of connectivity.In rural areas, buildings and trees along the road often cause discontinuities of road segments on remote sensing images.The OOC method extracts road regions based on the spectral and texture features of image objects, and therefore cannot identify those segments sheltered by buildings or trees.By contrast, the GFI method first identifies road segments with clear boundaries and then connects nearby segments by knowledge inference, so it overcomes the discontinuity problem. Figure 18 demonstrates the extraction results of both methods on a road with a large portion of shelter by trees.(3) By comparison, we can easily find that the result of OOC missed quite a few roads.It is probably because the pavement of rural roads is made of various kinds of materials: asphalt, cement, gravel, stone, etc.The material heterogeneity makes it hard to select appropriate samples and establish classification rules in the OOC method.On the other hand, the GFI method is good at handling roads of heterogeneous material since it is based mainly on geometric features.(4) Because the GFI method extracts road segments by identifying parallel segments that conform to road characteristics, it is sensitive to parallel elements in the remote sensing image and can lead to wrong extractions.Figure 19 shows the extraction results in a village.It can be seen from Figure 19b that a large farmland ridge is wrongly identified as a road by GFI because it has two obvious boundary lines, and its width is similar to a road's.However, the OOC method obtains the correct result since the farmland ridge has different spectral and texture features from a road.Spectral and texture information would be taken into account in our method in future to increase the extraction accuracy, especially for high-resolution images with more spectral bands.
ISPRS Int.J. Geo-Inf.2017, 6, 314 20 of 23 (5) Another aspect of the GFI method we need to pay attention to is that its extraction result is not very accurate when processing road segments with an irregular shape.The extraction result of GFI is formed by combinations of parallel lines, so it cannot express a complicated intersection well.For example, the road width gradually changes at the crossroad in the upper left corner of Figure 18, but the extraction result of GFI is only a parallel line.In future research, we would recognize the shape of intersection regions by image segmentation algorithms, and combine them with the GFI results to extract accurate road intersections.
ISPRS Int.J. Geo-Inf.2017, 6, 314 20 of 23 GFI is formed by combinations of parallel lines, so it cannot express a complicated intersection well.For example, the road width gradually changes at the crossroad in the upper left corner of Figure 18, but the extraction result of GFI is only a parallel line.In future research, we would recognize the shape of intersection regions by image segmentation algorithms, and combine them with the GFI results to extract accurate road intersections.Note that the image we used has not been orthorectified, and the acquisition angle and topography can produce some errors of location and shape of road detection.However, this does not affect the procedure of rural road extraction and detection accuracy.Future research is needed to investigate the influence of unorthorectified images on the detection results.In addition, the proposed method is designed to detect double edge lines of roads; therefore, it is only applicable to relatively wide roads that appear as ribbons rather than lines on images.Finally, a prerequisite of GFI is formed by combinations of parallel lines, so it cannot express a complicated intersection well.For example, the road width gradually changes at the crossroad in the upper left corner of Figure 18, but the extraction result of GFI is only a parallel line.In future research, we would recognize the shape of intersection regions by image segmentation algorithms, and combine them with the GFI results to extract accurate road intersections.Note that the image we used has not been orthorectified, and the acquisition angle and topography can produce some errors of location and shape of road detection.However, this does not affect the procedure of rural road extraction and detection accuracy.Future research is needed to investigate the influence of unorthorectified images on the detection results.In addition, the proposed method is designed to detect double edge lines of roads; therefore, it is only applicable to Note that the image we used has not been orthorectified, and the acquisition angle and topography can produce some errors of location and shape of road detection.However, this does not affect the procedure of rural road extraction and detection accuracy.Future research is needed to investigate the influence of unorthorectified images on the detection results.In addition, the proposed method is designed to detect double edge lines of roads; therefore, it is only applicable to relatively wide roads that appear as ribbons rather than lines on images.Finally, a prerequisite of our method is that the width of a road is approximately constant.A slight change in road width would not influence the detection result too much since our method uses the Ribbon snake algorithm to link road segments.However, the extraction result would be not very accurate when processing road segments with an irregular shape (variable width).

Conclusions
The current road extraction methods based on remote sensing image have been widely studied, and scholars have made great progress for different types of regions, image data sources, and data types; however, little research has focused on rural road extraction.In view of the fact that rural roads have diverse materials, large curvature changes, and a severe shelter problem, we propose a road extraction method based on geometric feature inference.The method makes full use of the geometric characteristics of roads, and constructs a geometric knowledge base of rural roads using the information of selected sample road segments.Based on the knowledge base, it identifies the parallel line pairs conforming to the characteristics of road edges, conducts grouping and connection on parallel lines instructed by knowledge reasoning, and finally extracts rural road vectors.We applied our method to a rural area in Xiangtan City, Hunan Province, China.By comparing our extraction result with the ground truth dataset, both the recall and accuracy rates were over 80%.The comparison between our method and the object-oriented and pixel-based classification method proves the advantages of the proposed method in handling roads in rural areas even if there are a few deficiencies.The proposed method is useful for providing a technical basis for large-scale rural road statistics and survey, and can inform national traffic construction and management.
ISPRS Int.J. Geo-Inf.2017, 6, 314 3 of 23 knowledge base is constructed; (2) image feature enhancement.The Gaussian filter is used to smooth the original image to remove random noises and enhance road-related features; (3) edge detection.

.
Based on the sample road segments, the maximum likelihood estimation method is used to calculate the two parameters (standard deviation σ and mathematical expectation μ ) of road width distribution.After the distribution is obtained, the Kolmogorov-Smirnov test is used to verify if the observed width values accord with the fitted distribution model.If the K-S statistic n D is smaller than test statistic D α ( 0.10 α = ), we think the quality of samples are satisfactory, otherwise, more samples are needed.
a is the slope, and b is the intercept.The parameters a and b are treated as variables, to construct the parameter space AOB, and then the above-mentioned line equations are converted to is the two lines passing through ( , ) a b in the parameter space, as shown in Figure2b.The Hough Transform detects lines based on the global features in images, so it is robust in detecting the discontinuous boundaries of regions caused by noise disturbance.

Figure 2 .
Figure 2. Points and lines in (a) image space and (b) parameter space.

Figure 2 .
Figure 2. Points and lines in (a) image space and (b) parameter space.

Figure 3 .
Figure 3. Distribution of parallel index.

Figure 3 .
Figure 3. Distribution of parallel index.

Figure 4 .
Figure 4. Distribution of extension index.

Figure 4 .
Figure 4. Distribution of extension index.

Figure 4 .
Figure 4. Distribution of extension index.

2. 5 .
Adjacent Segment Linking Candidate road segments (parallel line pair) include left boundary line LB and right boundary line RB , as shown in Figure 6.Their polar angles and radius are L L θ ρ , and R R θ ρ , respectively, and then there exists a theoretical center line, whose direction and width are (

2. 5 .
Adjacent Segment Linking Candidate road segments (parallel line pair) include left boundary line LB and right boundary line RB, as shown in Figure 6.Their polar angles and radius are θ L , ρ L and θ R , ρ R respectively, and then there exists a theoretical center line, whose direction and width are θ C = (θ R + θ L )/2 and width = |ρ L − ρ R |.

Figure 6 .
Figure 6.Geometric relationship between two candidate segments.

Figure 6 .
Figure 6.Geometric relationship between two candidate segments.

Figure 8 .
Figure 8. High-resolution aerial image of the study area.

Figure 8 .
Figure 8. High-resolution aerial image of the study area.

Figure 8 .
Figure 8. High-resolution aerial image of the study area.

Figure 9 .
Figure 9. Sample road segments of the study area.(a) spatial distribution of all sample roads, (b) Sample road 1, (c) Sample road 2, (d) Sample road 3.

Figure 9 .
Figure 9. Sample road segments of the study area.(a) spatial distribution of all sample roads, (b) Sample road 1, (c) Sample road 2, (d) Sample road 3.

Figure 11 .
Figure 11.Road extraction process.(a) edge detection result, (b) candidate segment identification result, (c) candidate segment linking result, (d) inference of missing segments, (e) overlay of remote sensing images and extracted roads

Figure 12 .
Figure 12.Road extraction result in the study area.

Figure 12 .
Figure 12.Road extraction result in the study area.

23 Figure 12 .
Figure 12.Road extraction result in the study area.

Figure
Figure Comparison of edge detection operators: (a) original image, (b) Sobel operator, (c) LOG operator, (d) Canny operator, (e) edge detection with embedded confidence.

Figure 15 .
Figure 15.Road extraction result of the object-oriented classification method.

Figure 16 .
Figure 16.Road extraction result of the pixel-based classification method.

Figure 15 . 23 Figure 15 .
Figure 15.Road extraction result of the object-oriented classification method.

Figure 16 .
Figure 16.Road extraction result of the pixel-based classification method.

Figure 16 .
Figure 16.Road extraction result of the pixel-based classification method.

Figure 17 .
Figure 17.Overlap of the GFI result (green) and the OOC result (red).

Figure 17 .
Figure 17.Overlap of the GFI result (green) and the OOC result (red).