Road Information Extraction from High-Resolution Remote Sensing Images Based on Road Reconstruction

Traditional road extraction algorithms, which focus on improving the accuracy of road surfaces, cannot overcome the interference of shelter caused by vegetation, buildings, and shadows. In this paper, we extract the roads via road centerline extraction, road width extraction, broken centerline connection, and road reconstruction. We use a multiscale segmentation algorithm to segment the images, and feature extraction to get the initial road. The fast marching method (FMM) algorithm is employed to obtain the boundary distance field and the source distance field, and the branch backing-tracking method is used to acquire the initial centerline. Road width of each initial centerline is calculated by combining the boundary distance fields, before a tensor field is applied for connecting the broken centerline to gain the final centerline. The final centerline is matched with its road width when the final road is reconstructed. Three experimental results show that the proposed method improves the accuracy of the centerline and solves the problem of broken centerline, and that the method reconstructing the roads is excellent for maintain their integrity.


Introduction
As the main body of modern transportation system, roads mean geographical, political, and economic significance, which are also the main recorded and identified object in GIS and maps.The technology of road extraction using images appeared in the middle of the 20th century, due to the need for intelligent transportation [1].With the rapid development of remote sensing in the recent 10 years, high-resolution remote sensing images provide the possibility of using remote sensing images for road extraction.Common available satellite sensors include the QuickBird, WorldView2, WorldView3, and IKONS launched by the United States; the SPOT by France; and the GF-2 by China.
Scholars put forward many algorithms and models, according to the different characteristics of image sources, image segmentation, rule sets of road extraction, and purposes of usage.The methods can be divided into four categories: (1) Road extraction method based on region segmentation; (2) Road extraction method based on template matching; (3) Road extraction method based on edge; and (4) Road extraction method based on multifeature fusion.
(1) Road extraction method based on region segmentation: this method uses threshold segmentation to obtain the road area as per the gray statistical feature of road image [2][3][4][5][6].(2) Road extraction method based on template matching: this method mainly segments and trims roads by local prior knowledge, and obtains road surfaces by the similarity between computer image and road template.The template matching method is more effective but semi-automatic, as its given template needs to be realized [7][8][9][10].(3) Road extraction method based on edge: this methods treats the essential features of roads as a set of parallel lines, based on which a lot of road extraction methods are proposed [11].(4) Road extraction method based on multifeature fusion: this method manifests that the road not only has spectral and texture characteristics, but also topological and contextual relations, and that the comprehensive use of road surface spectrum features, edge gradient features, and edge shape features allows better road extraction, and directs the development trend of road extraction algorithm [12][13][14].
In addition to direct road extraction, according to different extraction purposes, scholars presented the algorithm of road centerline extraction, using the road centerline to express the characteristics of road network.Zelang Miaoet et al., who were devoted to the research of road centerline extraction, proposed a geodesic method and the method combining Gaussian mixture model (GMM) with subspace constraint mean shift (SCMS) to extract the road centerline [15,16].Ruyi Liu et al. and Lipeng Gao et al. used the tensor voting algorithm to extract the road centerline to ensure the centerline integrity [17].Robert Van Uiterta et al. employed the fast marching method combined with the branch back-tracking method to obtain subpixel precise centerlines [18].A novel road centerline extraction method was proposed in [19] by integrating multiscale spectral structural features, support vector machines (SVMs), and morphological thinning algorithm.Guangliang Cheng et al. [20] used multiscale segmentation to extract roads, and the tensor voting algorithm to obtain centerlines.
The first step of road information extraction is often to segment the image.Traditional image segmentation algorithms, including the classical Otsu threshold segmentation [21] algorithm, FCM clustering algorithm [22], K-means algorithm [23], and watershed algorithm [24], can only segment grayscale images.However, high-resolution remote sensing images provide high spatial resolution and multiple bands with abundant spectral information, so the multiscale segmentation algorithm is used to segment the images [25].On the basis of studying the extraction method of road centerline, this paper selects the subvoxel precise skeleton extraction algorithm proposed by Robert Van Uiterta et al. to extract the road centerline, and the method is improved by the FMM algorithm, which extracts the centerline of organs efficiently and accurately.Level set methods evolve an isosurface in the direction of the surface normal, which was first proposed by Osher, Sethian [26] to determine the evolution of curves and surfaces through the level sets evolution of functions.The fast marching method (FMM) proposed by Sethian [27,28] is a fast level set algorithm, which can be understood as the solution of eikonal equation.There are faults for extracted centerlines, given the influence of algorithm accuracy and shadow occlusion, thus, this paper selects the tensor voting algorithm to connect the broken centerlines.The tensor voting, as a technique widely used in computer vision field, embodies outstanding performance in detecting geometrical features [17,20,29,30], which is a process of delivering each pixel's geometrical feature to its neighborhood.Thereby, the exchange of information at both ends of a fault centerline is enabled, so that the broken centerlines can be connected.The tensor voting algorithm has been widely accepted in road feature extraction methods, and has been demonstrated to be a good approach to repair broken roads, and extract road centerlines and intersection characteristics.
At present, although the previous research on road extraction based on high resolution remote sensing images has made great progress, the integrity of the extracted road is difficult to guarantee due to the influence of other landcovers.This paper aims to propose a road extraction method, which can overcome the problem of incomplete or fractured road due to shadow occlusion of buildings and buildings.The innovative work of this paper is as follows: (1) The statistical region merge (SRM) segmentation algorithm is used to guarantee the integrity of the image by choosing the proper segmentation scale [31].(2) Proposed a road width extraction method to reconstruct the road by using the road width and centerline.(3) Fast marching method algorithm (FMM) [18] is adopted and improved to obtain initial centerline and the final centerline, respectively.(4) Tensor voting is used to connect the broken centerline [17,32].

Methodology
The proposed algorithms include (1) The preprocessing of the experimental data; (2) The initial road obtained by using multiscale segmentation and feature extraction; (3) The FMM used to obtain the boundary distance field and the source distance field, and the branch backing-tracking method used to get the initial centerline; (4) The road width of each centerline obtained by utilizing the boundary distance field and the initial road; (5) The tensor field used to connect the broken centerlines to get the final centerline; (6) The final centerline matched with its path width, and the final road reconstructed.The flowchart of this algorithm is shown in Figure 1.

Methodology
The proposed algorithms include (1) The preprocessing of the experimental data;

Multiscale Segmentation and Feature Extraction
Though the high-resolution remote sensing image has rich band information, previous algorithms often made use of the image gray when using the high-resolution remote sensing image, which results in the rich spectral information not being fully exploited.This paper uses statistical region merging (SRM) [31] multiscale segmentation algorithm to segment the image, which not only fully utilizes the spectral features of high-resolution remote sensing images, but also guarantees the integrity of the image by choosing the proper segmentation scale.After the multiscale segmentation, according to the statistical analysis of road characteristics, we know that the road is generally continuous, and the surface radiation characteristics are consistent, and have a big length and width ratio value.On the basis of multiscale segmentation, a binary image is obtained by setting a specific threshold range on the multiscale segmentation result.The road area and aspect ratio feature are used to remove the non-road interference.Considering that the road is a continuous area and the road generally occupies a large proportion of the overall image, we can set the area parameters to eliminate the small area of non-road interference.In addition, according to the shape characteristics of the road, the road is generally slender, then we set the length-width

Multiscale Segmentation and Feature Extraction
Though the high-resolution remote sensing image has rich band information, previous algorithms often made use of the image gray when using the high-resolution remote sensing image, which results in the rich spectral information not being fully exploited.This paper uses statistical region merging (SRM) [31] multiscale segmentation algorithm to segment the image, which not only fully utilizes the spectral features of high-resolution remote sensing images, but also guarantees the integrity of the image by choosing the proper segmentation scale.After the multiscale segmentation, according to the statistical analysis of road characteristics, we know that the road is generally continuous, and the surface radiation characteristics are consistent, and have a big length and width ratio value.On the basis of multiscale segmentation, a binary image is obtained by setting a specific threshold range on the multiscale segmentation result.The road area and aspect ratio feature are used to remove the non-road interference.Considering that the road is a continuous area and the road generally occupies a large proportion of the overall image, we can set the area parameters to eliminate the small area of non-road interference.In addition, according to the shape characteristics of the road, the road is generally slender, then we set the length-width ratio of the road greater than 1.5 to remove the large non-road area interference.The result of initial road extraction is shown in Figure 2.   In multiple dimensions, the equation of evolution is:

The Extraction of Centerline
Skeleton extraction refers to the topological structure of a two-dimensional or three-dimensional object with a one-dimensional connected centerline, and the centerline of the road can be extracted by extracting the skeleton of the road [27].On the basis of the acquired initial road from Section 2.1, this paper uses FMM to obtain the boundary distance field and the source distance field, and uses the minimum cost path algorithm to obtain the initialization centerline, as shown in Figure 3.

The Extraction of Centerline
Skeleton extraction refers to the topological structure of a two-dimensional or three-dimensional object with a one-dimensional connected centerline, and the centerline of the road can be extracted by extracting the skeleton of the road [27].On the basis of the acquired initial road from Section 2.1, this paper uses FMM to obtain the boundary distance field and the source distance field, and uses the minimum cost path algorithm to obtain the initialization centerline， as shown in Figure 3.

Obtaining Distance Field Using FMM Algorithm
Fast marching method (FMM) was introduced by Sethian [33,34] for a special case of front evolution.In this paper, the FMM algorithm is used to obtain the boundary distance field and the source distance field.Considering an interface  (that is, the boundary of road or the maximum value point of the boundary in this paper),  can be a point, a curve or a surface.Consider the case of a front moving with speed , where F is always either positive or negative.Then, let T(x) be the time at which  crosses the point x.In the one-dimensional case, where distance = speed  time, the equation of motion can be expressed as In multiple dimensions, the equation of evolution is:

Obtaining Distance Field Using FMM Algorithm
Fast marching method (FMM) was introduced by Sethian [33,34] for a special case of front evolution.In this paper, the FMM algorithm is used to obtain the boundary distance field and the source distance field.Considering an interface Γ (that is, the boundary of road or the maximum value point of the boundary in this paper), Γ can be a point, a curve or a surface.Consider the case of a front moving with speed F(x), where F is always either positive or negative.Then, let T(x) be the time at which Γ crosses the point x.In the one-dimensional case, where distance = speed × time, the equation of motion can be expressed as In multiple dimensions, the equation of evolution is: where F(x) represents the speed when Γ passes by the point x, ∇T(x) represents the time gradient of Γ crossing the point x, and T(Γ 0 ) = 0 indicates that the initial arrival time of Γ is 0.
The process of solving the FMM algorithm is the process of obtaining T(x).According to Equation (2), the only changing parameter is the speed function F(x).In this way, we can get a different T(x) by setting different speed functions.The initial position of the evolution curve and the speed of the evolution curve are the initial conditions of the curve evolution.We can understand the construction process of the distance field as follows: Γ evolves from the initial position in accordance with speed function F(x), continually, through all the pixels in the image.The cumulative time map of Γ reached from all pixels is the desired distance field.a. Calculating the boundary distance field D 1 The boundary distance field D 1 of the road describes the relative spatial distances between the pixels and the boundary.The outermost pixels of the road are used as the road boundaries and the pixels between boundaries are used as the road, because the road boundary occupies a pixel width, the width of the road is smaller than the actual value.As shown in Figure 4b, with q = 3 structural elements to dilate the road boundary (Figure 4c), the pixels occupied by the newly acquired boundary is only used as a road boundary and not as part of the road.The corresponding value at the road boundary in the boundary distance field D 1 should be 0. When F(x) = 1, T(x) is the Euclidean distance from the boundary to the current position, so a uniform speed function is defined as where x is the spatial position of the pixel in the image, and F D 1 (x) is the speed at which the FMM evolves at that pixel.Substitute F D The process of solving the FMM algorithm is the process of obtaining T(x) .According to Equation (2), the only changing parameter is the speed function F(x).In this way, we can get a different T(x) by setting different speed functions.The initial position of the evolution curve and the speed of the evolution curve are the initial conditions of the curve evolution.We can understand the construction process of the distance field as follows: evolves from the initial position in accordance with speed function F(x), continually, through all the pixels in the image.The cumulative time map of  reached from all pixels is the desired distance field.
a. Calculating the boundary distance field D1 The boundary distance field D1 of the road describes the relative spatial distances between the pixels and the boundary.The outermost pixels of the road are used as the road boundaries and the pixels between boundaries are used as the road, because the road boundary occupies a pixel width, the width of the road is smaller than the actual value.As shown in Figure 4b, with q = 3 structural elements to dilate the road boundary (Figure 4c), the pixels occupied by the newly acquired boundary is only used as a road boundary and not as part of the road.The corresponding value at the road boundary in the boundary distance field D1 should be 0. When ( ) Euclidean distance from the boundary to the current position, so a uniform speed function is defined as where x is the spatial position of the pixel in the image, and where D 1 is boundary distance field, and D max is the maximum value of the boundary distance field.
According to the definition of boundary distance field, the nearer to the centerline, the larger D 1 is.Hence, F D2 (x) has a greater value while near the centerline, that is, the nearer to the centerline, the faster the curve evolves and the smaller ∇T(x) is. Figure 5 is a schematic diagram of the boundary distance field and source distance field, and we assume that Figure 5a is the initial road obtained by image segmentation.For the road boundary from Figure 5a, the evolution is performed at F D 1 (x), and the boundary distance field can be obtained by masking the area outside the road, as shown in Figure 5b.For the maximum distance point E of the boundary distance field, the evolution is performed at F D 1 (x), and the area outside the road is masked to obtain the source distance field D Y , as shown in Figure 5c.The source point E evolves with speed function F D2 (x), and the area outside the road is masked to obtain the source distance field D Y , as shown in Figure 5d.
where 1 D is boundary distance field, and Dmax is the maximum value of the boundary distance field.According to the definition of boundary distance field, the nearer to the centerline, the larger x has a greater value while near the centerline, that is, the nearer to the centerline, the faster the curve evolves and the smaller

Centerline Extraction Based on Branch Backing-Tracking Method
The maximum point S of the source distance field DY is selected as the starting point, and the maximum value point E of the boundary distance field is chosen as the end point, and the QSE is defined as all paths from the starting point S to end point E.Then, a branch backing-tracking

Centerline Extraction Based on Branch Backing-Tracking Method
The maximum point S of the source distance field D Y is selected as the starting point, and the maximum value point E of the boundary distance field is chosen as the end point, and the Q SE is defined as all paths from the starting point S to end point E.Then, a branch backing-tracking method [18,25] is used for back-tracking the minimum cost path from the starting point S to the end point E on the source distance field D T , which is the centerline, as shown in Figure 6.method [18,25] is used for back-tracking the minimum cost path from the starting point S to the end point E on the source distance field DT, which is the centerline, as shown in Figure 6.From the given start point S to the end point E, a minimum cost path Q(t) can be found: [0,  ]  R n , where the cumulative cost function U is the function of the pixel position x, and the minimum cumulative cost is defined as where L is the cumulative distance (arrival time) from the starting point S along Q to the end point

E;
Sx Q is the set of all possible paths from the starting point S to x; ( ) U x , as a cost function, is defined as ( ) U x has a lower cost nearer the centerline and a higher cost away from the centerline, from which we can find the centerline of the road by finding the minimum value.
The extraction of centerline is the solution of the ordinary differential Equation (6).
where Q(t)) is the centerline, and E is the farthest point.From the starting point S to the end point E, the minimum-cost back-tracking path is the centerline.For the current location Qx, the next location where Qx+1 is the next point; ( ) x T Q is the local gradient of the current position; and w is the step.
The smaller w is, the higher the accuracy of the centerline, here, 0.01 is often chosen as the step.
Then, the source distance field DY and DT are obtained by using the centerline Q(t) as the source.
After the acquisition of the new source distance field, the starting point S and end point E are selected as the same way above, then, the branch backing-tracking method is used to obtain the centerline and looped until all the centerline is found.In this way, the initialization centerline C1 is obtained, and the centerline is shown as in Figure 6b.From the given start point S to the end point E, a minimum cost path Q(t) can be found: [0,∞]→R n , where the cumulative cost function U is the function of the pixel position x, and the minimum cumulative cost is defined as

Road Width Extraction
where L is the cumulative distance (arrival time) from the starting point S along Q to the end point E; Q Sx is the set of all possible paths from the starting point S to x; U(x), as a cost function, is defined as x) ; and U(x) has a lower cost nearer the centerline and a higher cost away from the centerline, from which we can find the centerline of the road by finding the minimum value.
The extraction of centerline is the solution of the ordinary differential Equation (6).
where Q(t) is the centerline, and E is the farthest point.From the starting point S to the end point E, the minimum-cost back-tracking path is the centerline.For the current location Q x , the next location point where Q x+1 is the next point; T(Q x ) is the local gradient of the current position; and w is the step.The smaller w is, the higher the accuracy of the centerline, here, 0.01 is often chosen as the step.
Then, the source distance field D Y and D T are obtained by using the centerline Q(t) as the source.
After the acquisition of the new source distance field, the starting point S and end point E are selected as the same way above, then, the branch backing-tracking method is used to obtain the centerline and looped until all the centerline is found.In this way, the initialization centerline C 1 is obtained, and the centerline is shown as in Figure 6b.

Road Width Extraction
As one of the key parameters of road extraction, road width is often neglected in previous studies.Due to the influence caused by vegetation, building shadows near roads, and other disturbances, the completeness of roads will be missing.In addition, for the effect of the surface features whose materials are similar to the road, the extracted roads will be wider than the actual roads.Assuming that the missing part of the road and the added part of the road are stochastic, this paper uses the coordinates of the centerline and the Euclidean distance from the centerline to the road boundary, to calculate the width of the road for each point in the centerline.We know that the width of a road is consistent, so we calculate the average width of each road branch as the road width D r .
From the previous section, the initial centerline C 1 has been obtained in Section 2.2.2.As shown in Figure 7, the road is extracted from an image of 400 × 400 pixels, where Figure 7a is the boundary distance field, and the closer to the centerline position, the larger of the value of the boundary distance field is, and the maximum value acquired at the centerline position; Figure 7b is the contour line of Figure 7a.For any point on the centerline of the road, suppose the value of boundary distance field x w is the distance from this point to road boundary, then the road width is D r = 2 × x w .The definition of road width is shown in Figure 7c.
studies.Due to the influence caused by vegetation, building shadows near roads, and other disturbances, the completeness of roads will be missing.In addition, for the effect of the surface features whose materials are similar to the road, the extracted roads will be wider than the actual roads.Assuming that the missing part of the road and the added part of the road are stochastic, this paper uses the coordinates of the centerline and the Euclidean distance from the centerline to the road boundary, to calculate the width of the road for each point in the centerline.We know that the width of a road is consistent, so we calculate the average width of each road branch as the road width Dr.
From the previous section, the initial centerline C1 has been obtained in Section 2.2.2.As shown in Figure 7, the road is extracted from an image of 400  400 pixels, where Figure 7a

Obtaining the Tensor Field
Although all the centerlines in the initial road have been detected, there are still many broken places in the centerlines.From the original image, we can know this part of the fault is caused by shadow and the cover of vegetation.The phenomenon of shadow, vehicle, vegetation, and other occlusions in the road is ubiquitous, and it is impossible to compensate by the optimized segmentation algorithm, so the tensor field is introduced for the link of the broken centerline.
Tensor voting [35] is an algorithm that represents the saliency structure of an image, a method of reasoning its implied geometry from a large amount of unreliable information by passing (voting) tensor (differential geometric information) between adjacent points.In the algorithm, the neighborhood information is gathered together by the way of voting, and the orientation information is estimated by the calculation rules of the vector field.The flowchart for tensor voting is shown in Figure 8.

Obtaining the Tensor Field
Although all the centerlines in the initial road have been detected, there are still many broken places in the centerlines.From the original image, we can know this part of the fault is caused by shadow and the cover of vegetation.The phenomenon of shadow, vehicle, vegetation, and other occlusions in the road is ubiquitous, and it is impossible to compensate by the optimized segmentation algorithm, so the tensor field is introduced for the link of the broken centerline.
Tensor voting [35] is an algorithm that represents the saliency structure of an image, a method of reasoning its implied geometry from a large amount of unreliable information by passing (voting) tensor (differential geometric information) between adjacent points.In the algorithm, the neighborhood information is gathered together by the way of voting, and the orientation information is estimated by the calculation rules of the vector field.The flowchart for tensor voting is shown in Figure 8.  if the input pixel is a single pixel point, then the eigenvalues are l1 = l2 = 1,

263
the input pixel is a point on the curve and the normal vector of the point is , its eigenvalues are l1 = 1, l2 = 0, . According to the definition of matrix theory, (1) In the tensor encoding phase, the location information of the initialization centerline C 1 is encoded as the initialization tensor T. For the initial centerline of the input (binary image): if the input pixel is a single pixel point, then the eigenvalues are l 1 = l 2 = 1, T = 1 0 0 1 ; if the input pixel is a point on the curve and the normal vector of the point is n = (n x , n y ), its eigenvalues are l 1 = 1, l 2 = 0, T = n x n x n x n y n x n y n y n y .According to the definition of matrix theory, any one second-order semidefinite tensor can be decomposed into two eigenvalues and two eigenvectors.The tensor T can be decomposed as follows: where e 1 and e 2 are eigenvectors, l 1 and l 2 are eigenvalues, l 1 ≥ l 2 ≥ 0, e 1 e 1 T is the stick tensor, (l 1 − l 2 ) is the saliency of stick tensor; e 1 e 1 T + e 2 e 2 T is the ball tensor; and l 2 is the saliency of ball tensor.(2) Tensor voting is carried out on the initial tensor after encoding, and tensors communicate their information with each other in a neighborhood.A tensor vote is first made on the initialization tensor obtained from 1).The tensor field is expressed by the attenuation function DF.As shown in Figure 9a, we assume that O and P are two points in T, where O is the voting point and P is the receiving point, then, the voting intensity DF of P received from O point is as shown in Equation ( 9).DF(s, k, σ) = e −((s 2 +ck 2 )/σ 2 ) (9 where s is the osculating circle arc length of the voting center point to the voting point; k is the curvature of the osculating circle; σ is the size of the voting window; and c is a constant.The value of the voting window σ is determined by the size of the centerline fracture, and, if the fracture is large, then σ takes a relatively large value; here, it takes a value of 18.25.Where c is the constant associated with σ, the equation is c = −16 ln(0.1 × (σ − 1)/π 2 ).if the input pixel is a single pixel point, then the eigenvalues are l1 = l2 = 1, the input pixel is a point on the curve and the normal vector of the point is , its eigenvalues are l1 = 1, l2 = 0, According to the definition of matrix theory, any one second-order semidefinite tensor can be decomposed into two eigenvalues and two eigenvectors.The tensor T can be decomposed as follows: T e e e e e e e e e e (8) where e1 and e2 are eigenvectors, l1 and l2 are eigenvalues, e e is the stick tensor, (l1 − l2)is the saliency of stick tensor; e e e e is the ball tensor; and l2 is the saliency of ball tensor.
2) Tensor voting is carried out on the initial tensor after encoding, and tensors communicate their information with each other in a neighborhood.A tensor vote is first made on the initialization tensor obtained from 1).

Final Centerline Extraction and Road Reconstruction
From Section 2.2, we know that the boundary distance field is a Euclidean distance field, which is used as the input of the initial centerline extraction.However, the distance field is faulty because of the segmentation algorithm and the condition of the road itself.The D 3 obtained in Section 2.3 is similar to the boundary distance field D 1 , that is, the pixels closer the centerline have higher values.Thus, the boundary distance field D 3 , acquired in the Section 2.4, is used to extract the final centerline, and the speed function is redefined as where D 3 is the boundary distance field, and D 3max is the maximum value of D 3 .
After that, we look for the maximum point E 2 for D 3 , and redo the processes in Section 2.2.1(b) and Section 2.2.2.With the source distance field D Y and D T updating constantly and centerline merging continuously, we get the final centerline C 2 , which is connected in the broken.
As shown in Figure 10b,c, the distances from different branches of the initial centerline C 1 to C 2 are calculated, and the road of the branch centerline with the minimum distance is selected as the road width of C 2 .For all the final connected centerlines in the road network, according to the final centerline, the C 2 position matches the initial centerline C 1 position and the corresponding road width D r .Here, the broken road width the in C 1 obtained from the surrounding road width.Through widening final centerline C 2 to the corresponding road width D r , the complete road W is obtained, and the complete road information extraction is realized.
value of the voting window  is determined by the size of the centerline fracture, and, if the fracture is large, then  takes a relatively large value; here, it takes a value of 18.25.Where c is the constant associated with  , the equation is T e e e e e e (10)

Final Centerline Extraction and Road Reconstruction
From Section 2.2, we know that the boundary distance field is a Euclidean distance field, which is used as the input of the initial centerline extraction.However, the distance field is faulty because of the segmentation algorithm and the condition of the road itself.The D3 obtained in Section 2.3 is similar to the boundary distance field D1, that is, the pixels closer the centerline have higher values.
Thus, the boundary distance field D3, acquired in the Section 2.3, is used to extract the final centerline, and the speed function is redefined as where D3 is the boundary distance field, and 3 max D is the maximum value of D3.
After that, we look for the maximum point E2 for D3, and redo the processes in Section 2.

Experiment
Three sets of data from difference part of urban are selected in order to demonstrate the validity of the proposed algorithm, which are all from the website http://www.cse.iitm.ac.in/~{}vplab/satellite. html.In order to make the proposed method better used in urban regional scale road extraction, this paper selects urban images under different circumstances to extract the road.The data are QuickBird images fused by panchromatic and multispectral images, whose sizes are all 512 × 512 pixels with a resolution of 0.61 m.The experimental results are shown in Figures 11-13, where (a) is the picture fused by panchromatic and multispectral images; (b) is the initial road, which is extracted by SRM multiscale segmentation and feature extraction; (c) is the initial centerline, in which the broken centerlines are marked with red rectangles; (d) is the tensor field obtained by the initial centerline; (e) is the complete centerline after the connection, with an unconnected broken centerline marked with red rectangle; and (f) is the extracted road, which is reconstructed by the final centerline and road width.

Experiment on Image 1
For image 1, the road is in the field, and the road widths are from various different sources.Several sites of the roads are affected by the overpass shadows, and some are obscured by vegetation.The roads and centerlines in Figure 11b,c are broken by the shadows and vegetation.By the utilization of tensor voting, the broken centerlines are well connected.As shown in Figure 11e, the breakage in the red rectangle is too large to be connected.
experimental results are shown in Figure 11

Experiment on Image 2
For the image 2, the road is in a residential area.The materials on the main road and several examples of landcover before the buildings are the same, which causes great confusion for road extraction.In addition, there are other road materials at the junctions, which easily leads to road breakage.As shown in Figure 12b, many non-road features have been misjudged as roads, due to the disturbance of the landcover with the same materials as the road.As shown in Figure 12c-e, the broken centerlines are all connected, which greatly improves the integrity of the road extraction.

Experiment on Image 3
For image 3, the road is located between a residential area's wastelands.On the way to the house, there are similar small paths which are characterized as the same as the road.This experiment treats these paths as roads.Many parts of the road are disturbed by the shading of vegetation and buildings, as shown in Figure 13a, which results in 12 breakages, nine of which are connected by tensor voting.
For the image 2, the road is in a residential area.The materials on the main road and several examples of landcover before the buildings are the same, which causes great confusion for road extraction.In addition, there are other road materials at the junctions, which easily leads to road breakage.As shown in Figure 12b, many non-road features have been misjudged as roads, due to the disturbance of the landcover with the same materials as the road.As shown in Figure 12c-e

Experiment on Image 3
For image 3, the road is located between a residential area's wastelands.On the way to the house, there are similar small paths which are characterized as the same as the road.This

370
The manual drawing method is used to obtain reference road and reference centerline in order 371 to evaluate the quality of the road extraction algorithm.The reference data are shown in Figure 14.

Analysis of Centerline Extraction Results
The manual drawing method is used to obtain reference road and reference centerline in order to evaluate the quality of the road extraction algorithm.The reference data are shown in Figure 14.

Analysis of Centerline Extraction Results
The manual drawing method is used to obtain reference road and reference centerline in order to evaluate the quality of the road extraction algorithm.The reference data are shown in Figure 14.For the centerline, the centerline appraisal method proposed by Wiedemmann et al. [36] is used to evaluate the centerline result.The equations are as follows: /( ) where Com, Cor, and Q represent completeness, correct, and quality, respectively.TP is the number of pixels both in the centerline and the reference centerline; FP is the number of pixels in the centerline but not in the reference centerline; and FN is the number of pixels in the reference centerline but not in the centerline.The number of TP, FP, and FN obtained from our method are exhibited in Table 1.Based on the experimental results of the four methods in Table 2, we can know the Com, Cor, and Q of our algorithm for road centerline extraction are all higher than those of the previous three algorithms.According to the data statistics in Table 1, FP is always greater than FN because the tensor voting can solve, well, the faults of the centerline, and greatly improve the precision of centerline extraction.However, the tensor voting algorithm also extends the centerline at the end of For the centerline, the centerline appraisal method proposed by Wiedemmann et al. [36] is used to evaluate the centerline result.The equations are as follows: Cor = TP/(TP + FN) where Com, Cor, and Q represent completeness, correct, and quality, respectively.TP is the number of pixels both in the centerline and the reference centerline; FP is the number of pixels in the centerline but not in the reference centerline; and FN is the number of pixels in the reference centerline but not in the centerline.The number of TP, FP, and FN obtained from our method are exhibited in Table 1.Based on the experimental results of the four methods in Table 2, we can know the Com, Cor, and Q of our algorithm for road centerline extraction are all higher than those of the previous three algorithms.According to the data statistics in Table 1, FP is always greater than FN because the tensor voting can solve, well, the faults of the centerline, and greatly improve the precision of centerline extraction.However, the tensor voting algorithm also extends the centerline at the end of road, allowing the centerline to be greater than its actual value at the end, which also results in the completeness always being less than the accuracy.Ten samples are randomly selected for road width statistics on the reference road, as shown in Figure 15, so as to verify the accuracy of the extracted road width.With the purpose of ensuring the representativeness of the selected road width samples, the road width of three positions are measured separately, and the average value of the three positions is utilized as the reference road width.Then, the reference road widths are compared with the road widths, as shown in Figure 11f, and the width data statistics are displayed in Table 3, where the mean absolute error of the road width is 4.324%, and the average accuracy of the road width is 95.676%.The effective extraction of road width demonstrates that the proposed algorithm can extract complete road information well.width.Then, the reference road widths are compared with the road widths, as shown in Figure 11f, and the width data statistics are displayed in Table 3, where the mean absolute error of the road width is 4.324%, and the average accuracy of the road width is 95.676%.The effective extraction of road width demonstrates that the proposed algorithm can extract complete road information well.The proposed research focuses on the accurate extraction of road centerline and road width.As a main feature in urban construction, roads, with their rational planning and construction, are of great significance to urban development.The traditional method of road information collection is mainly carried out by manual means, where the road information is updated slowly and consumes enormous human physical and financial resources.Nowadays, the rapid development of satellites and the birth of satellite remote sensing images, with a resolution better than 1 m, provide the possibility to use remote sensing images for road extraction, and the remote sensing images for road extraction has the advantages of timely information, wide range, and low cost, compared with traditional methods of road information collection.
For example, the City Prosperity Index (CPI) is an important data for sustainable development analysis [37][38][39].Here, roads are an important indicator for the CPI, especially in the following aspects [40]: Economic Density: Economic density is the city product [41] divided by the city's area, its definition is as follows: Econmic Density = City Product ($PPP) City square kilometers , Where city square kilometer includes: residential and non-residential buildings, major roads, railways, and sport facilities.Here, the proposed road extraction method can be used to calculate the area of the road within square kilometers effectively.

2.
Population Density [42].High-density neighborhoods tend to decrease the costs of public services.So, a reasonable urban planning of high population density is the key to determine urban sustainable development, which includes: police and emergency response, school transport, roads, water and sewage, etc.

3.
Street Intersection Density.It means number of street intersections per one square kilometer of urban area.In order to calculate this index, we need to: Obtain the street network map of the urban area; Verify the topology: each street segment must be properly connected to other segments; Obtain the start and end point of each segment; Collect events from start and end points; Exclude points with less than 3 events; Count the remaining points and divide by the urban area in km 2 .For the above information, the method of getting road network is the key technique.Using the proposed road extraction, can efficiently obtain the urban road network information and road details, to provide data support for the street intersection density.4.
Street Density.Number of kilometers of urban streets per square kilometer of land which is defined as: Street density = Total length urban streets Total o f urban sur f ace .The proposed algorithm can obtain the high precision road centerline, we can obtain the city street length through computing centerline length.In this way, the street density can be obtained effectively.

5.
Land Allocated to Streets.Land allocated to streets means total area of urban surface allocated to streets.It was defined as: Street density = 100 Total sur f ace o f urban streets Total sur f ace o f urban area .Here, we propose a road extraction algorithm that can directly carry out Road area statistics and get total surface of urban streets.
Remote sensing image road extraction method is not only an important measure for UN-Habitat to measure the sustainable development of urban, it also plays a very important role in people's daily life.Construction of smart cities makes the efficient road extraction method an urgent need.According to the wide range of remote sensing images, we can use the high-resolution remote sensing images to: • Update the GPS system information timely.Road is one of the most important information of GPS, the use of road extraction from high-resolution remote sensing image can update GPS information timely.

Comparison of the Methods for Centerline Extraction
Previous studies [2][3][4][5][6] have shown that the high-resolution remote sensing imagery for road extraction still has many problems.Unlike the remote sensing imagery with rich spectral information, traditional image segmentation algorithms tend to only consider the information of single bands, which are difficult for guaranteeing the accuracy of image segmentation.In addition, due to the vegetation shading, building shadows, and other effects, it is easy to cause centerline fractures and incomplete road surfaces.In view of these problems, this study does its work from three aspects: using the multiscale segmentation algorithm to enhance the accuracy of image segmentation; employing the tensor voting algorithm to connect broken centerlines; and combining road width and centerline to reconstruct roads.
In Figure 10c-e, Figure 11c-e, and Figure 12c-e, the tensor vote algorithm was applied for connecting broken centerlines.However, for large centerline fractures, such as the red rectangles in Figures 10e and 12e, the algorithm was still unable to connect the broken centerlines.This study used Com, Cor, and Q to evaluate the extracted centerline.The Experiment 1 showed high values of Com, Cor, and Q, whereas Experiment 2 had mistaken some non-road areas as roads, resulting in Com less than Cor, but all three indicators were above 90%.The background of the Experiment 3 data was more complex, because, during the initial road extraction, there was an unextracted road, due to the disturbance of the vegetation shading, which resulted in lower Com and Q.
In the previous study of road extraction algorithms [10,17,29], many scholars had studied the extraction of road centerline.In order to fully verify the superiority of this algorithm, the centerline algorithm was compared with Ruyi Liu's algorithm [17], Chaudhuri's [10] algorithm, and Shi's [29] algorithm.In Chaudhuri's method, Chaudhuri segmented the enhanced image and obtained the road by the artificial template with small noise removed.Finally, the road centerline was obtained through morphological thinning.In Shi's method, Shi obtained the initial road based on the spectral-spatial classification result, and fused the road groups and homogeneous property to improve road-group accuracy.Then, the local linear kernel smoothing regression algorithm was used to obtain the centerline of road after a morphological processing.In Ruyi Liu's method, initial road regions were obtained by combining the shear transform algorithm with directional segmentation.Then, the road map based on Mahalanobis distance and thresholding was fused with the initial road regions to improve accuracy.Finally, the road centerline was extracted by the skeleton algorithm proposed by Uitert and Bitter.
The comparison results are shown in Figure 16, where (a) is the experimental data; (b) is the reference centerline extracted manually; (c) is the centerline extracted by the Chaudhuri's method; (d) is the centerline by the Shi's method; (e) is the centerline by the Ruyi Liu's method; and (f) is the centerline by the proposed algorithm.Compared with the other three algorithms, this method used the tensor voting algorithm to connect the broken centerline, as shown in Figure 15, which greatly improves the completeness of the centerline.The results of the evaluation are exhibited in Table 4 where the completeness, correctness, and quality of the proposed algorithm are much better than those of the previous three algorithms.Based on the accurate extraction of the centerline, the paper studied the road width extraction and received a good result.In this study, the distance from each point on the centerline to the edge of the road was calculated to find the road width.The road width was replaced by the average of each section of road width according to the continuity of the road.For the road width extraction, ten samples were selected randomly in Section 3.3, and the calculated road width and the road width on the reference road were compared and analyzed.As shown in Table 3, the average accuracy of road Compared with the other three algorithms, this method used the tensor voting algorithm to connect the broken centerline, as shown in Figure 15, which greatly improves the completeness of the centerline.The results of the evaluation are exhibited in Table 4 where the completeness, correctness, and quality of the proposed algorithm are much better than those of the previous three algorithms.

Assessment of Road Width Extraction
Based on the accurate extraction of the centerline, the paper studied the road width extraction and received a good result.In this study, the distance from each point on the centerline to the edge of the road was calculated to find the road width.The road width was replaced by the average of each section of road width according to the continuity of the road.For the road width extraction, ten samples were selected randomly in Section 3.3, and the calculated road width and the road width on the reference road were compared and analyzed.As shown in Table 3, the average accuracy of road width extraction was as high as 95.676%, which indicates that the method of using centerline and road width for reconstruction is appropriate.
In general, the proposed method improves the accuracy of the centerline extraction, and the method can extract the road width accurately.Using the centerline and road width to reconstruct roads can solve, well, the problem of incomplete road extraction caused by shading.

Conclusions
A new method of road extraction is presented in this paper, where the multiscale segmentation algorithm makes full use of the spectral information in high-resolution remote sensing images.Using the tensor voting algorithm to connect fault centerlines greatly improves the integrity of centerline.The method of reconstructing roads with road width and centerline information solves the problem that the road integrity is affected by vegetation and shadows.
To verify the effectiveness of the proposed algorithm, three QuickBird images were tested.Parameters of completeness, correctness, and quality were used to evaluate the centerline extraction results.In addition, we compared the centerline extraction results with the other three centerline extraction algorithms.The comparison results show that the introduction of the tensor voting algorithm solves the problem of broken centerline effectively.For road width, ten samples were selected on the road, with the comparison between the extracted road width and the road width on the reference data.The analysis results demonstrate that the average extraction accuracy of the road width is 95.68%.
In future work, we will enhance the segmentation algorithm to improve segmentation accuracy, study shadow removal algorithms to reduce road extraction interference, and improve the algorithms of broken centerline connection.Characteristics of roads on multispectral and hyperspectral images will also be attached with great importance, by using spectral features to improve road extraction accuracy.

( 2 )FMMFigure 1 .
Figure 1.Flowchart of the proposed road information extraction based on road reconstruction.

Figure 1 .
Figure 1.Flowchart of the proposed road information extraction based on road reconstruction.

Figure 2 .
Figure 2. The extraction of initial road.(a) Original image; (b) Initial road after multiscale

135 2 . 2 . 1 .
Obtaining Distance Field Using FMM Algorithm 136 Fast marching method (FMM) was introduced by Sethian [33,34] for a special case of front 137 evolution.In this paper, the FMM algorithm is used to obtain the boundary distance field and the 138 source distance field.Considering an interface  (that is, the boundary of road or the maximum 139 value point of the boundary in this paper),  can be a point, a curve or a surface.Consider the 140 case of a front moving with speed , where F is always either positive or negative.Then, let 141 T(x) be the time at which  crosses the point x.In the one-dimensional case, where distance = speed 142  time, the equation of motion can be expressed as

Figure 2 .
Figure 2. The extraction of initial road.(a) Original image; (b) Initial road after multiscale segmentation and feature extraction.

Figure 2 .
Figure 2. The extraction of initial road.(a) Original image; (b) Initial road after multiscale segmentation and feature extraction.

Figure 4 .Figure 4 .
Figure 4.The extraction of boundary distance field.(a) The boundary of initial road; (b) Local magnification image of the boundary, which shows a width of a pixel occupied by the boundary; (c) The values of boundary distance from (b). b.Calculating the source distance field DY and DT for initial road The source distance field records the spatial position relation from each pixel to the source point.On the basis of the boundary distance field D1, the maximum position of the boundary distance field E1 with the value 1 max D is obtained, and the source distance fields DY and DT are obtained by using this point as the source point.The source point is the starting point of the source distance field calculation, which is the position of the zero-level set.The speed function of DY is defined the same as 1 ( ) D F x , so the DY records the Euclidean distance from each point to the source point.The speed function of DT is 2 ( ) D F x :

Figure 5
is a schematic diagram of the boundary distance field and source distance field, and we assume that Figure 5a is the initial road obtained by image segmentation.For the road boundary from Figure 5a, the evolution is performed at 1 ( ) D F x , and the boundary distance field can be obtained by masking the area outside the road, as shown in Figure 5b.For the maximum distance point E of the boundary distance field, the evolution is performed at 1 ( ) D F x , and the area outside the road is masked to obtain the source distance field DY, as shown in Figure 5c.The source point E evolves with speed function 2 ( ) D F x , and the area outside the road is masked to obtain the source distance field DY, as shown in Figure 5d.

Figure 5 .
Figure 5. (a) The initial road; (b) Boundary distance field of initial road, the point E has the maximum value of the boundary distance field with the value of 1 max D

Figure 5 .
Figure 5. (a) The initial road; (b) Boundary distance field of initial road, the point E has the maximum value of the boundary distance field with the value of D 1max and E is selected as the source point of D Y and D T ; (c) The source distance field D Y with the speed function F D 1 (x), the point S has the maximum value of D Y and S is chosen as the start point in Section 2.2.2;(d) The source distance field D T with the speed function F D2 (x).

Figure 6 .
Figure 6.Centerline extraction.(a) The schematic diagram [18] of the source distance field DT with the speed of

Figure 6 .
Figure 6.Centerline extraction.(a) The schematic diagram [18] of the source distance field D T with the speed of F D2 (x)-the darker pixels represent lower values, whereas the lighter indicates higher values; (b) The centerlines of different road branches shown in different colors.
is the boundary distance field, and the closer to the centerline position, the larger of the value of the boundary distance field is, and the maximum value acquired at the centerline position; Figure 7b is the contour line of Figure 7a.For any point on the centerline of the road, suppose the value of boundary distance field xw is the distance from this point to road boundary, then the road width is Dr = 2  xw.The definition of road width is shown in Figure 7c.

Figure 7 .
Figure 7. (a) The boundary distance field of road; (b) The contour line of road; (c) The definition of road width; (d) The branch with difference road width.

Figure 7 .
Figure 7. (a) The boundary distance field of road; (b) The contour line of road; (c) The definition of road width; (d) The branch with difference road width.

260 1 )
In the tensor encoding phase, the location information of the initialization centerline C1 is 261 encoded as the initialization tensor T. For the initial centerline of the input (binary image):

Figure 8 . 1 )
Figure 8. Flowchart of tensor voting.1) In the tensor encoding phase, the location information of the initialization centerline C1 is encoded as the initialization tensor T. For the initial centerline of the input (binary image):

Figure 9 .( 3 )
Figure 9. Tensor voting.(a) Schematic diagram of tensor voting; (b) Broken centerline; (c) The result of tensor voting.(3) Each point in the encoded initial centerline C 1 votes to the point neighborhood within its tensor field, and accepts other points to vote for itself, getting a new tensor field T S .(4) After the voting, each point on the centerline has a new tensor field, and Equation (10) is used to decompose the tensor field T S .After decomposition, we can get the stick tensor saliency map l S1 − l S2 and the ball tensor saliency map l 2 .Here, a stick tensor saliency map is used as boundary distance field in the centerline extraction process, and the new boundary distance field is called D 3 .In consideration of the fact that D 3 is decaying without boundary, here, we set the value of the positions in D 3 less than 0.3 to 0, and the value of about 0.3 indicates that the tensor field of the point is very small and negligible.The effect of tensor voting is shown in Figure 9c.T S = (l S1 − l S2 )e S1 e S2 T + l 2 (e S1 e S1T + e S2 e S2 T )(10)

. 3 ) 4 )
Each point in the encoded initial centerline C1 votes to the point neighborhood within its tensor field, and accepts other points to vote for itself, getting a new tensor field TS.After the voting, each point on the centerline has a new tensor field, and Equation (10) is used to decompose the tensor field TS.After decomposition, we can get the stick tensor saliency map tensor saliency map 2 l .Here, a stick tensor saliency map is used as boundary distance field in the centerline extraction process, and the new boundary distance field is called D3.In consideration of the fact that D3 is decaying without boundary, here, we set the value of the positions in D3 less than 0.3 to 0, and the value of about 0.3 indicates that the tensor field of the point is very small and negligible.The effect of tensor voting is shown in Figure9c.
2.1 (b) and Section 2.2.2.With the source distance field DY and DT updating constantly and centerline merging continuously, we get the final centerline C2, which is connected in the broken.

Figure 10 .
Figure 10.Matching road width.(a) Initial centerline; (b) Initial centerline shown in the initial road, where C1a, C1b, and C1c are three different branches of initial road, and Dr1, Dr2, and Dr3 are the corresponding widths of the three branches; (c) The connected centerline, (c) The reconstructed road.

Figure 10 .
Figure 10.Matching road width.(a) Initial centerline; (b) Initial centerline shown in the initial road, where C 1a , C 1b , and C 1c are three different branches of initial road, and D r1 , D r2 , and D r3 are the corresponding widths of the three branches; (c) The connected centerline, (c) The reconstructed road.
-13, where (a) is the picture fused by panchromatic and multispectral images; (b) is the initial road, which is extracted by SRM multiscale segmentation and feature extraction; (c) is the initial centerline, in which the broken centerlines are marked with red rectangles; (d) is the tensor field obtained by the initial centerline; (e) is the complete centerline after the connection, with an unconnected broken centerline marked with red rectangle; and (f) is the extracted road, which is reconstructed by the final centerline and road width.
, the broken centerlines are all connected, which greatly improves the integrity of the road extraction.

Figure 14 .
Figure 14.Reference road and reference centerline.(a) Reference road of experiment 1; (b) Reference centerline of experiment 1; (c) Reference road of experiment 2; (d) Reference centerline of experiment 2; (e) Reference road of experiment 3; (f) Reference centerline of experiment 3.

Figure 14 .
Figure 14.Reference road and reference centerline.(a) Reference road of experiment 1; (b) Reference centerline of experiment 1; (c) Reference road of experiment 2; (d) Reference centerline of experiment 2; (e) Reference road of experiment 3; (f) Reference centerline of experiment 3.

Figure 15 .
Figure 15.The samples of reference road width.

Figure 15 .
Figure 15.The samples of reference road width.

Figure 16 .
Figure 16.Comparison of four methods.(a) Experiment data; (b) Reference centerline; (c) Results of Chaudhuri's method; (d) Results of Shi's method; (e) Results of Ruyi's method; and (f) Results of our method.

Table 1 .
Statistics of TP, FP, and FN.

Table 2 .
The result of the comparison of four methods.

Table 1 .
Statistics of TP, FP, and FN.

Table 2 .
The result of the comparison of four methods.

Table 3 .
The assessment of road width extraction.
4.1.Significance of Road Extraction Based on Remote Sensing Image

Table 3 .
The assessment of road width extraction.

•
Conduct traffic situation assessment.According to the distribution of road and vehicle to figure out the congestion of road, not only to provide owners with the best travel program but also greatly improve the efficiency of traffic management.•Preparestatisticsforurbanroad distribution and land consumption rates, so as to facilitate urban planning.Road and building are the two most important components of urban.Recently, the global population is growing rapidly, taking up more and more resources on the earth, and remote sensing image road extraction can be used in statistics for the distribution of urban roads, and contribute to the rational planning and sustainable development of cities.•Evaluate the urban population density and urban economic growth based on road development.The intensity of road network is positively correlated with the degree of urban population and economic development.The use of road statistics can establish population density and economic development models, used to assess the comprehensive strength of urban development.• Implement urban change monitoring.Remote sensing images keep the historical data of the city well, analyze the road information obtained from remote sensing images at different times, and can analyze the changing trends of urban development effectively.

Table 4 .
The result comparison of four method.

Table 4 .
The result comparison of four method.