A Method for Road Network Extraction from High-Resolution SAR Imagery Using Direction Grouping and Curve Fitting

: Roads are an important recognition target in synthetic aperture radar (SAR) image interpretation. Although a considerable number of high-quality SAR images are now available, the method of road extraction is lagging. To extract the road network with low missed and false rates, this paper proposed a road network extraction approach which includes line detection, road segmentation, road network extraction and optimization. First, the linear feature response and direction map are obtained from the SAR intensity image using the multiplicative Duda operation. Then, the backscattering coe ﬃ cient and coe ﬃ cient of variation are combined using a support vector machine to eliminate the linear structures of non-roads, and the binary image of road candidates is subsequently achieved by morphological proﬁles of path openings. Next, with the obtained direction map, a novel thinning method based on binary image decomposition and curve ﬁtting is presented to obtain line segments of the road network. Finally, a series of measures which involve overlap, continuity, and junction optimization are proposed to construct the road network. In the experiments, the proposed method was applied to Radarsat-2 and TerraSAR-X high-resolution images. The experimental results showed that the proposed method had an excellent performance in terms of both completeness and correctness.


Introduction
Roads are typical man-made objects that play an essential role in modern transportation systems. Especially in an emergency, for instance, when serious natural disasters occur, real-time road information is vital for emergency rescue. Therefore, accurate and timely road information extraction is considered to be essential. Synthetic aperture radar (SAR), with its imaging capability in all-weather and all-time, has become a popular tool to obtain road information. With the development of SAR technology, SAR sensors have been able to provide daily high-quality images with different modes. Taking TerraSAR-X as an example, the resolution can be of the order of 20 cm in staring spotlight mode [1]. Nevertheless, complete road detection from SAR images is still a challenging task because of the unique imaging mechanism of SAR, the diversity of road types, and the variable complex situations surrounding the road.
In SAR images with high resolution, road targets are often characterized by elongated dark strips with specified widths, instead of thin lines as the case of low-resolution SAR images [2]. Considering these road features, many different kinds of methods have been presented to detect road from SAR images in the past decades. Among those approaches, a significant number are comprise two main steps:

Segmentation of Road Candidates
For illustrative purposes, the two-dimensional image in this paper is described as x M y N x y Z where ( , ) A x y is the value of the pixel located at row x and column y , M and N are the height and width of the image respectively. Considering the linear indexing technique, the two-dimensional image can be also expressed as To acquire road candidates from high-resolution SAR images, the proposed method is based on [20] which uses multiplicative Duda operators and path operators. More specifically, the multiplicative Duda operator [5] is applied to obtain the linear feature response (LFR) using local sliding windows with different sizes and orientations. Then, road candidates are segmented by constructing morphological profiles using path openings [21] on the LFR map. However, due to the fact that only features of local intensity ratio and target length are considered, some non-road targets may be interpreted as roads when SAR images contain shadows and dense buildings. In order to tackle the drawback, backscattering coefficient ( σ 0 ), and coefficient of variation (CV) are merged to reduce the linear structures of non-roads in the LFR map. In this work, the σ 0 map is obtained as

Segmentation of Road Candidates
For illustrative purposes, the two-dimensional image in this paper is described as A = {A(x, y)|1 ≤ x ≤ M; 1 ≤ y ≤ N; x, y ∈ Z} where A(x, y) is the value of the pixel located at row x and column y, M and N are the height and width of the image respectively. Considering the linear indexing technique, the two-dimensional image can be also expressed as A = {a i |a i = A(x, y), i = N × (x − 1) + y}. Then, we define that B = {b i |b i ∈ {0, 1}} denotes the binary image and T = {t i |t i ∈ {0, 1, 2, . . . , N t }} is the labeled image, where b i and t i represent the pixel values at the site i, N t denotes the total number of labels.
To acquire road candidates from high-resolution SAR images, the proposed method is based on [20] which uses multiplicative Duda operators and path operators. More specifically, the multiplicative Duda operator [5] is applied to obtain the linear feature response (LFR) using local sliding windows with different sizes and orientations. Then, road candidates are segmented by constructing morphological profiles using path openings [21] on the LFR map. However, due to the fact that only features of local intensity ratio and target length are considered, some non-road targets may be interpreted as roads when SAR images contain shadows and dense buildings. In order to tackle the drawback, Remote Sens. 2019, 11,2733 4 of 17 backscattering coefficient (σ 0 ), and coefficient of variation (CV) are merged to reduce the linear structures of non-roads in the LFR map. In this work, the σ 0 map is obtained as P = {σ i |σ i = 10 × log 10 [argmax µ C { f (i|V(θ, w, l))}]} (1) where f (i|V(θ, w, l)) is the response of multiplicative Duda operation with the window V(θ, w, l) at site i (See Figure 2), µ C is the mean intensity value of region C in the window V which has the maximum response. Since roads have low surface roughness and low permittivity [22], the value σ 0 of roads is very low and often comparable to the SAR sensor noise level. The CV measures the heterogeneity of a region and is defined as where γ C = s C /µ C , µ C and s C denote the mean and standard deviation of intensity values in the region C of window V which has the maximum response. For high-resolution SAR images, roads can be considered as homogeneous patches [23] and thus have low values of CV. Figure 1 shows the examples of LFR, σ 0 and CV maps obtained from the test image.
After the extraction of σ 0 and CV, the feature vector U i = (σ i , cv i ) is constructed and the SVM classifier is subsequently applied to obtain the non-road mask. In this study, all the features are normalized and a Gaussian RBF kernel is selected. The SVM is chosen because it has the ability to generalize well even with limited training samples. Besides, a considerable number of studies have demonstrate the performance of SVM in remote sensing applications [24]. Finally, the values of non-road pixels in the LFR map are reset to zero before morphological profiles with path operators are applied, which analyze object length and shape features to determine road candidates.  (1) where ( | ( , , )) f i V w l θ is the response of multiplicative Duda operation with the window θ ( , , ) V w l at site i (See Figure 2), μ C is the mean intensity value of region C in the window V which has the maximum response. Since roads have low surface roughness and low permittivity [22], the value σ 0 of roads is very low and often comparable to the SAR sensor noise level. The CV measures the heterogeneity of a region and is defined as where γ μ , μ C and C s denote the mean and standard deviation of intensity values in the region C of window V which has the maximum response. For high-resolution SAR images, roads can be considered as homogeneous patches [23] and thus have low values of CV. Figure 1 shows the examples of LFR, σ 0 and CV maps obtained from the test image.
After the extraction of σ 0 and CV, the feature vector ( , ) is constructed and the SVM classifier is subsequently applied to obtain the non-road mask. In this study, all the features are normalized and a Gaussian RBF kernel is selected. The SVM is chosen because it has the ability to generalize well even with limited training samples. Besides, a considerable number of studies have demonstrate the performance of SVM in remote sensing applications [24]. Finally, the values of nonroad pixels in the LFR map are reset to zero before morphological profiles with path operators are applied, which analyze object length and shape features to determine road candidates.

Generation of Road Network.
This part aims to extract the road network by global optimization using line segments as elements. The thinning road network is the foundation for embedding the extracted road information into a geographic information system (GIS) database [25]. In this article, the road network is creatively modeled as a sequence of connected curves which have determined equations. The proposed method is done with two main steps: binary image decomposition through direction grouping and thinning by curve fitting. To begin with, the following operations on images are defined.
Given two binary images 1 2 , B B , "AND (&) " and "OR (||) " operations for images are defined as performing logical "AND" and "OR" on each pixel at the same location. That is, For a labeled image T , we define a binary conversion to extract the pixels that belong to a label set G as follows:

Generation of Road Network
This part aims to extract the road network by global optimization using line segments as elements. The thinning road network is the foundation for embedding the extracted road information into a geographic information system (GIS) database [25]. In this article, the road network is creatively modeled as a sequence of connected curves which have determined equations. The proposed method is done with two main steps: binary image decomposition through direction grouping and thinning by curve fitting. To begin with, the following operations on images are defined.
Given two binary images B 1 , B 2 , "AND (&)" and "OR (||)" operations for images are defined as performing logical "AND" and "OR" on each pixel at the same location. That is, For a labeled image T, we define a binary conversion to extract the pixels that belong to a label set G as follows: where:

Binary Image Decomposition with Direction Feature Grouping
For a binary image B of road candidates, connected-component (CC) labeling [26] can be applied which scans the binary image and groups its pixels into components based on pixel connectivity. However, the shape of CC may be rather complex and irregular because pixels in the same CC may be oriented differently. This leads to a tricky problem of describing the CC using a curve equation. To tackle this problem, binary image decomposition is proposed.
Suppose that G = {0, π/n, 2π/n, . . . , (n − 1)π/n} is a set of values that represent different direction levels. Then, the direction feature map is estimated by where f (i|V(θ, w, l)) is the response of multiplicative Duda operation with the window V(θ, w, l) at site i, θ ∈ G. As shown in Figure 1, the estimated direction map are marked by different colors representing different directions. According to the proximity of directions, the direction set G is divided into K(K ≤ n) groups which are denoted by G k (k = 1, 2, . . . , K). For example, in the case of 8 direction levels (n = 8) and 4 groups (K = 4), it can be defined that G 1 = {0, π/8},G 2 = {π/4, 3π/8},G 3 = {π/2, 5π/8}, and G 4 = {3π/4, 7π/8}. Next, the binary image B is decomposed into K binary image layers where B T,G k = {b i |b i = δ(t i , G k ), t i ∈ T}. Figure 1 presents the example of binary image layers decomposed from road candidates. It is certain that all the road candidates on an image layer B k have similar directional characteristics that belongs to the set G k . Most importantly, this step makes it possible to construct a coordinate system and apply the method of curve fitting on each CC. This process can be considered as a decomposition process of binary image B because we have B = B 1 ||B 2 || . . . ||B K and

Thinning Based on Polynomial Curve Fitting
After the decomposition of the binary image, the road candidate regions need to be thinned to achieve line segments, which are the basic elements of the road network. The traditional morphology thinning algorithm is a commonly used solution, but it has poor performance when the edge is rather rough and the image is seriously disturbed by noise. In this study, a thinning based on polynomial curve fitting is proposed to extract the thinned road line segments. Algorithm 1 gives the implementation details of the proposed thinning procedure. Step 1: Input all the road binary image layers B k (k = 1, 2, . . . , K).
Step 2: For k = 1, 2, . . . , K: a) Apply the method of connected-component (CC) labeling on B k ; b) Compute the pixel area of each CC and remove CCs in which pixel area is smaller than T area ; c) Create coordinate system, for each CC: 1) Extract all pixel coordinates in the CC as the discrete data points; 2) Apply the polynomial curves fitting on data points; 3) Record the curve equation, two endpoints, and corresponding tangential direction vectors; 4) Generate the line segment as the thinned curve of CC using the curve equation.
Step 3: Combine all the line segments to form the thinning road map.
The example of extracting the fitted curves and endpoints is illustrated in Figure 3a. Supposing that Q = {(x n , y n )|x n , y n ∈ N + , n = 1, 2, . . . , N} represents a set of pixel coordinates in a CC where N denotes the total number of pixels, the polynomial curves fitting is applied on Q to derive curve equation f (x) = c 0 + c 1 x + . . . + c p x p where p is the equation order and c 0 , c 1 , . . . , c p are parameters determined by the least square method. Furthermore, two endpoints can be obtained by coordinates (x min , f (x min )) and (x max , f (x max )) where x min and x max are the minimum and maximum of x n in Q. The corresponding angles of the tangent lines at two endpoints are expressed as θ 1 = π + arctan f (x min ), θ 2 = arctan f (x max ) and then tangential direction vectors are s 1 = (cos θ 1 , sin θ 1 ) and s 2 = (cos θ 2 , sin θ 2 ). The endpoints and tangential direction vectors are both recorded here for the following road network optimization. It is worth mention that the definition of x-axis and y-axis of coordinate system should be exchanged for the case of binary images B k with a vertical direction feature. As shown in Figure 3b

Road Network Generation and Optimization
Although the thinned road map obtained in the previous section has initiated formation of a road network, it is coarse and need to be optimized further. On one hand, due to the fact that road regions are decomposed according to the direction feature grouping, multiple extractions for the same road may occur when extracting the line segments from the image layers. This particular case is manifested as the overlap of line segments (see Figure 3b). On the other hand, as presented in Figure 3b, undesired gaps may also occur in the thinning results. There are two main factors that may contribute to gaps. One is the SAR imaging principles which cause roads approximately in the azimuth direction to be easily affected by the layovers and shadows of tall buildings or vegetation, especially in the high-resolution SAR images [27]. The other factor is the extraction approaches, e.g., image filtering and segmentation, which may remove some true road components.
In order to resolve these problems, the optimization of overlap, continuity, and junction are proposed successively. These proposed approaches of road network optimization are based on perceptual grouping principles which include proximity, parallelism, colinearity, continuity, and so on [28]. In terms of the network generalization, the perceptual grouping concepts have been widely

Road Network Generation and Optimization
Although the thinned road map obtained in the previous section has initiated formation of a road network, it is coarse and need to be optimized further. On one hand, due to the fact that road regions are decomposed according to the direction feature grouping, multiple extractions for the same road may occur when extracting the line segments from the image layers. This particular case is manifested as the overlap of line segments (see Figure 3b). On the other hand, as presented in Figure 3b, undesired gaps may also occur in the thinning results. There are two main factors that may contribute to gaps. One is the SAR imaging principles which cause roads approximately in the azimuth direction to be easily affected by the layovers and shadows of tall buildings or vegetation, especially in the high-resolution SAR images [27]. The other factor is the extraction approaches, e.g., image filtering and segmentation, which may remove some true road components.
In order to resolve these problems, the optimization of overlap, continuity, and junction are proposed successively. These proposed approaches of road network optimization are based on perceptual grouping principles which include proximity, parallelism, colinearity, continuity, and so on [28]. In terms of the network generalization, the perceptual grouping concepts have been widely used in image interpretation such as for road extraction [29] and river network generation [28]. In this article, the line segments will be perceived as elements and grouped for the optimization processes.

Overlap Optimization
As shown in Figure 4, the proposed overlap optimization is motivated by the matching principle presented by Wiedemann [30]. Let L m and L n denote two different line segments. A buffer with constant width is then constructed around L m . The parts of L n located in the buffer are considered as the overlapping parts. The overlapping ratio can be computed by r o = l n1 /(l n1 + l n2 ) where l n1 , l n2 are the lengths of L n within and without the buffer respectively. If the overlapping ratio is large than the given threshold T o , L m and L n will be merged into one line segments L o by polynomial curve fitting using the data points on the line segments L m and L n . The overlap optimization is finally implemented by repeating the above steps and constantly updating the line segments until there are no two line segments that satisfy the condition r o > T o where T o is 0.5 in this article. In practice, the buffer width can be set according to the linear width feature map W estmiated from where f (i|V(θ, w, l)) is the response of multiplicative Duda operation with the window V(θ, w, l) at site i. used in image interpretation such as for road extraction [29] and river network generation [28]. In this article, the line segments will be perceived as elements and grouped for the optimization processes.

Overlap Optimization
As shown in Figure 4, the proposed overlap optimization is motivated by the matching principle presented by Wiedemann [30].

Continuity Optimization
The continuity optimization is aimed at filling the gaps between line segments which have similar direction features. As illustrated in Figure 5, this procedure is ruled by two parameters: 1) D : the spatial Euclidean distances between endpoints that belong to the two target line segments using image coordinates; 2) Ω : the angle between tangential direction vectors s and connected direction vectors ω at the endpoints where ω is defined as the vector from one endpoint to another endpoint. For two given line segments , m n L L from the same image layer, the optimization is conducted if one of geometric constraint conditions below is satisfied:

Continuity Optimization
The continuity optimization is aimed at filling the gaps between line segments which have similar direction features. As illustrated in Figure 5, this procedure is ruled by two parameters: 1)D: the spatial Euclidean distances between endpoints that belong to the two target line segments using image coordinates; 2)Ω: the angle between tangential direction vectors s and connected direction vectors ω at the endpoints where ω is defined as the vector from one endpoint to another endpoint. For two given line segments L m , L n from the same image layer, the optimization is conducted if one of geometric constraint conditions below is satisfied: where T D1 , T D2 are two Euclidean distance thresholds and T Ω represents the angle threshold. Then, the optimization is done by merging line segments L m , L n into one L o using polynomial curves fitting. Finally, the above steps are repeated until no line segment pair satisfies the conditions in the whole image. For the experiments in this paper, D is computed using pixel coordinates, the distance threshold T D1 is equal to 10 and T D2 is 50, and the angle threshold T Ω is set to π/8.

Junction Optimization
Road junctions, which usually appear as the shapes of "L", "T", and "十", are significant parts in the road network. To resolve the gaps and burrs illustrated in Figure 6, a junction optimization method is proposed which can work well in the both cases. Let

Junction Optimization
Road junctions, which usually appear as the shapes of "L", "T", and "十", are significant parts in the road network. To resolve the gaps and burrs illustrated in Figure 6, a junction optimization method is proposed which can work well in the both cases. Let L m , L n denote two line segments derived from different image layers, which means they have different direction features and may form a junction. In order to determine the cross point J mn of L m , L n , one direct method is solving the combined curve equations of L m , L n . However, the solution will be complex and time-consuming when the equation order is high. Thus, the cross point is determined by searching the pixels located on the target curve. Suppose that S = {(x m , y m )} is the pixel coordinate set of L m and f n (x) is the curve equation of L n . Then, a potential cross point set can be determined by J = {(x m , y m )| f n (x m ) − y m < 1}. If J is not an empty set, the cross point is then determined by To further determine whether L m , L n constitute a junction and satisfy the condition of optimization, the spatial Euclidean distance D between J mn and E n is calculated using pixel coordinates where E n is the nearest endpoint of L n to the cross point J mn . If D < T D is satisfied, replace endpoint E n by J mn and redrawn line segment L n using the curve equation f n (x) and the new endpoint. Finally, iterate all the line segment pairs based on the above rules and the junction optimization is then finished. As can be noticed, T D is the only parameter in this step and is set to 20 in this paper. Normally, T D should not be too large or it may bring false optimization.

Junction Optimization
Road junctions, which usually appear as the shapes of "L", "T", and "十", are significant parts in the road network. To resolve the gaps and burrs illustrated in Figure 6, a junction optimization method is proposed which can work well in the both cases. Let

Post-Processing of Road Network
Since the geometric relationships are the only characteristics used in performing the network optimization, some of the false detections of roads still cannot be removed. For example, some low backscattering objects, such as rivers and shadows, are difficult to distinguish from roads because of their similarity of features in SAR images. It should be noted that those road line segments could be considered as a refined input processed by the methods in [9,10] which are based on Markov random fields (MRF) to obtain a better result. However, although the MRF methods are not the focus of this paper, the proposed approach is suitable for further improvement in the correctness of the road network.
Let x max x min f m (x)dx. Thus, the length-width ratio of L m is inferred by r lw (L m ) = l(L m )/w(L m ). Finally, the road network after post-processing is derived from where T lw is the length-width ratio threshold. This process can improve the correctness of road extraction because true road segments usually have high length-width ratios. Empirically, threshold T lw is fixed to be 3 in this paper.

Dataset Description and Parameter Setting
The experiments were conducted on high-resolution SAR images from two different satellite sensors, which are Radarsat-2 and TerraSAR-X. Table 1 lists the main characteristics of SAR images in experimental study sites. Study site 1 was collected from Radarsat-2 using ultrafine acquisition mode. Study site 1 is located in Dujiangyan, Sichuan, China, and covers an area of factories. Study site 2 lies in Chongzhou, Sichuan, southwest of China. It was acquired by TerraSAR-X sensor with staring spotlight mode and mainly contains road network and some villages. Study site 3 is a larger dataset located in Denver, Colorado, USA which includes more complex scenes such as cloverleaf junction, lakes, and a dense residential area. The original SAR images and geographic information are presented in Figure 7a, Figure 8a, and Figure 9a.  Table 2 presents the main parameters used in the experiments. For the parameter settings, some can be fixed without considerable loss for the most of the applications. For example, the direction levels of sliding windows are fixed at 8, i.e., G = {0, π/8, 2π/8, . . . , 7π/8}. When performing binary image decomposition, the image is decomposed into four image layers which contain road segments with direction features of {0, π/8}, {π/4, 3π/8}, {π/2, 5π/8} and {3π/4, 7π/8} respectively. The polynomial order p is 3 in the process of thinning based on polynomial curve fitting. As presented in previous sections, constant parameters are also used for the step of optimization in this article, which are estimated empirically. Except from the above mentioned parameters, there are some other parameters that are considered significant. A parameter is considered to be significant when it has obvious influence on the final result. Table 2 presents these significant parameters which include sliding window sizes used in the multiplicative Duda operation, and median gray level (MGL) together with minimum length (L min ) for the road candidate segmentation. The definitions of MGL and L min can be found in [20,21]. More analysis about these parameters is introduced in Section 5.3. To evaluate the proposed method, the results obtained are analyzed visually and numerically. image decomposition, the image is decomposed into four image layers which contain road segments with direction features of π π π π π {0, / 8},{ / 4, 3 / 8},{ / 2, 5 / 8} and π π {3 / 4,7 / 8} respectively. The polynomial order p is 3 in the process of thinning based on polynomial curve fitting. As presented in previous sections, constant parameters are also used for the step of optimization in this article, which are estimated empirically. Except from the above mentioned parameters, there are some other parameters that are considered significant. A parameter is considered to be significant when it has obvious influence on the final result. Table 2 presents these significant parameters which include sliding window sizes used in the multiplicative Duda operation, and median gray level ( MGL ) together with minimum length ( min L ) for the road candidate segmentation. The definitions of MGL and min L can be found in [20,21]. More analysis about these parameters is introduced in Section 5.3. To evaluate the proposed method, the results obtained are analyzed visually and numerically. For the visual analysis, the three original SAR images of test regions are shown in Figure 7a  For the numerical analysis, the quantitative indexes of the obtained results are computed using the reference data based on the method proposed by Wiedemann [30]. When applying the evaluation method, a buffer is set to determine the degree that one road network matches with the other. In this paper, the buffer is set to 5 pixels and three quantitative indexes (completeness, correctness, quality) are used to measure the detection performance. The completeness (CP) is the percentage of the reference road network located in the buffer around the extracted road network, while the correctness (CR) represents the percentage of the extracted road network located in the buffer around the reference road network. The quality (QL) index is a general index combining completeness and correctness. Table 3 proposed the quantitative indexes of the extracted results in the three test sites. It is worth mentioning that these quantitative indexes have more relative than absolute meaning [10] because they are relative to the manually selected reference data and buffer. For the numerical analysis, the quantitative indexes of the obtained results are computed using the reference data based on the method proposed by Wiedemann [30]. When applying the evaluation method, a buffer is set to determine the degree that one road network matches with the other. In this paper, the buffer is set to 5 pixels and three quantitative indexes (completeness, correctness, quality) are used to measure the detection performance. The completeness (CP) is the percentage of the reference road network located in the buffer around the extracted road network, while the correctness (CR) represents the percentage of the extracted road network located in the buffer around the reference road network. The quality (QL) index is a general index combining completeness and correctness. Table 3 proposed the quantitative indexes of the extracted results in the three test sites. It is worth mentioning that these quantitative indexes have more relative than absolute meaning [10] because they are relative to the manually selected reference data and buffer.  For the numerical analysis, the quantitative indexes of the obtained results are computed using the reference data based on the method proposed by Wiedemann [30]. When applying the evaluation method, a buffer is set to determine the degree that one road network matches with the other. In this paper, the buffer is set to 5 pixels and three quantitative indexes (completeness, correctness, quality) are used to measure the detection performance. The completeness (CP) is the percentage of the reference road network located in the buffer around the extracted road network, while the correctness (CR) represents the percentage of the extracted road network located in the buffer around the reference road network. The quality (QL) index is a general index combining completeness and correctness. Table 3 proposed the quantitative indexes of the extracted results in the three test sites. It is worth mentioning that these quantitative indexes have more relative than absolute meaning [10] because they are relative to the manually selected reference data and buffer.

Comparison of Methods for Road Segmentation
For road segmentation, the proposed method is compared with results from method in [21] and method of SVM which uses LFR, backscattering σ 0 , and CV as feature vector. Figure 10 presents the segmentation results using different methods in the three study sites. As shown in Table 3, the three method were also compared in terms of quantitative indexes. It should mentioned that all the quantitative indexes are computed using the same thinning method and parameters.
As can be observed from results in Figure 10, all three methods have the ability to detect the road candidates. However, a problem with the method in [21] can be clearly seen, namely that non-road dark regions with large area are interpreted as roads, which leads to a low correctness. This is mainly due to the fact that the method is based on the gray level value and object length. As for the method of SVM, a high CP but low CR can be noted in Table 3 because many isolated non-road pixels are detected, which is similar with the results of other pixel-based classifiers. In order to overcome the aforementioned defects, the proposed method improves the CR of results by using features of intensity ratio and linear length in spite of the loss of a small percentage of completeness. For the study site 3 with more complex scenes, it can be observed that the proposed method brings a greater decrease in CP than the other two sites. This fact occurs mainly because roads around the cloverleaf junction in study site 3 have a low intensity contrast, which causes the missed detection. Overall, the proposed method provides results with better QL, i.e., compromise results in terms of CP and CR.
For road segmentation, the proposed method is compared with results from method in [21] and method of SVM which uses LFR, backscattering 0 σ , and CV as feature vector. Figure 10 presents the segmentation results using different methods in the three study sites. As shown in Table 3, the three method were also compared in terms of quantitative indexes. It should mentioned that all the quantitative indexes are computed using the same thinning method and parameters.
As can be observed from results in Figure 10, all three methods have the ability to detect the road candidates. However, a problem with the method in [21] can be clearly seen, namely that nonroad dark regions with large area are interpreted as roads, which leads to a low correctness. This is mainly due to the fact that the method is based on the gray level value and object length. As for the method of SVM, a high CP but low CR can be noted in Table 3 because many isolated non-road pixels are detected, which is similar with the results of other pixel-based classifiers. In order to overcome the aforementioned defects, the proposed method improves the CR of results by using features of intensity ratio and linear length in spite of the loss of a small percentage of completeness. For the study site 3 with more complex scenes, it can be observed that the proposed method brings a greater decrease in CP than the other two sites. This fact occurs mainly because roads around the cloverleaf junction in study site 3 have a low intensity contrast, which causes the missed detection. Overall, the proposed method provides results with better QL, i.e., compromise results in terms of CP and CR.

Comparison of Methods for Network Generation
To assess the performance of proposed network generation method, a comparison was carried out with two state-of-the-art methods: (1) morphological thinning [31] and; (2) tensor voting algorithm [32]. Figure 11 presents the comparison results on two test binary images using different methods. With regard to the proposed thinning approach, it is based on curve fitting on connected Figure 10. Segmentation results of road candidates using method in [21], SVM, and proposed approach: (a-c) study site 1; (d-e) study site 2; (g-i) study site 3.

Comparison of Methods for Network Generation
To assess the performance of proposed network generation method, a comparison was carried out with two state-of-the-art methods: (1) morphological thinning [31] and; (2) tensor voting algorithm [32]. Figure 11 presents the comparison results on two test binary images using different methods. With regard to the proposed thinning approach, it is based on curve fitting on connected components (CC). Compared with traditional morphological thinning methods, the proposed method has several advantages. On the one hand, the thinning results of the proposed method are definitely curves with two endpoints and known curvilinear functions, while the morphology thinning method is sensitive to the contour of CC which may bring burrs and loops (See Figure 11a) instead of neat curves. For most practical applications, the desired result for road extraction is a smooth thinning curve. On the other hand, the proposed method can easily obtain the endpoints and their corresponding tangential extension direction (See Figure 3a) using the curve equation, which are useful for the subsequent network optimization. As shown in Figure 11c, the connected part that is emphasized is the result using the proposed network optimization. As can be seen from Figure 11b, the tensor voting method also generates smooth centerlines and has the ability to connect small gaps. However, the missing parts occur at junction areas which need to be improved in the post-processing. Although the proposed method has the aforementioned strengths, the input of additional direction information is mandatory which is achieved by the sliding window detection in this paper. If the direction feature is inaccurate or discontinuous, the road candidates in the decomposed image layers may be fragmented and the performance of presented work will decline. One consequence of this is the missing part shown in Figure 11c. Figure 11d,e shows one more example of the comparison. Experimental results show that the proposed method achieves an improved performance in network generation, especially for the case of SAR image in which extracted road segments have no smooth edges.
Remote Sens. 2019, 11,2733 14 of 18 may be fragmented and the performance of presented work will decline. One consequence of this is the missing part shown in Figure 11c. Figure 11d,e shows one more example of the comparison. Experimental results show that the proposed method achieves an improved performance in network generation, especially for the case of SAR image in which extracted road segments have no smooth edges.
Morphological thinning Tensor voting Proposed

Parameter Analysis and Discussion
Concerning the multiplicative Duda operation, window size is the key parameter. Usually, the parameter must be chosen to obtain responses with high contrast between lines and non-lines. As shown in Figure 12a, an image slice which is marked in red and contains line pixels is selected. Figure  12b shows the line feature responses of selected pixels under different window widths, while the other parameters remain the same. It can be inferred that window width has an obvious influence on

Parameter Analysis and Discussion
Concerning the multiplicative Duda operation, window size is the key parameter. Usually, the parameter must be chosen to obtain responses with high contrast between lines and non-lines. As shown in Figure 12a, an image slice which is marked in red and contains line pixels is selected. Figure 12b shows the line feature responses of selected pixels under different window widths, while the other parameters remain the same. It can be inferred that window width has an obvious influence on the detection results and the optimal should be equal to the actual pixel width of road in the image. Since road width is usually unknown in advance and it changes with different road types, windows with different widths should be applied and the maximum response value should be taken even though it may increase the cost of time. Additionally, parameter analysis of road segmentation was conducted. The performance of road segmentation is controlled by the pixel length threshold min L and the pixel value threshold MGL [20]. In this article, both of the threshold min L and MGL are selected empirically. In order to present the effects of the two parameters, the results using the proposed method with different values of min L and MGL are compared. The results are obtained using the image of study site 1 and presented in Figure 13. As illustrated in Figure 13(a), there is a decline of index CP with the increase of min L , which is opposite to the variation trend of CR. The index QL reaches peak probably at the intersection of two curve CP and CR. The index of completeness decreases because some short road segments are removed when the threshold min L increases. Besides, it can be seen that both of the CP and CR are above 60% when min L is in the range from 60 to 180 pixels. Moreover, Figure 13(b) shows the changes of quantitative indexes when MGL varies but min L is fixed at 100. The graph indicates that the completeness has the downward trend but the correctness keeps rising when MGL changes from 0.05 to 1. To make a compromise between completeness and correctness, the optimal MGL can be determined at the peak of quality curve. However, such a quality curve is hard to achieve in practical applications which need complete reference road data. Thus, according to the QL curve, the optimal MGL is estimated to be in the range of 0.2 to 0.4, which provides a guide for parameter selection in the other applications. Additionally, parameter analysis of road segmentation was conducted. The performance of road segmentation is controlled by the pixel length threshold L min and the pixel value threshold MGL [20]. In this article, both of the threshold L min and MGL are selected empirically. In order to present the effects of the two parameters, the results using the proposed method with different values of L min and MGL are compared. The results are obtained using the image of study site 1 and presented in Figure 13. As illustrated in Figure 13a, there is a decline of index CP with the increase of L min , which is opposite to the variation trend of CR. The index QL reaches peak probably at the intersection of two curve CP and CR. The index of completeness decreases because some short road segments are removed when the threshold L min increases. Besides, it can be seen that both of the CP and CR are above 60% when L min is in the range from 60 to 180 pixels. Moreover, Figure 13b shows the changes of quantitative indexes when MGL varies but L min is fixed at 100. The graph indicates that the completeness has the downward trend but the correctness keeps rising when MGL changes from 0.05 to 1. To make a compromise between completeness and correctness, the optimal MGL can be determined at the peak of quality curve. However, such a quality curve is hard to achieve in practical applications which need complete reference road data. Thus, according to the QL curve, the optimal MGL is estimated to be in the range of 0.2 to 0.4, which provides a guide for parameter selection in the other applications. the completeness has the downward trend but the correctness keeps rising when MGL changes from 0.05 to 1. To make a compromise between completeness and correctness, the optimal MGL can be determined at the peak of quality curve. However, such a quality curve is hard to achieve in practical applications which need complete reference road data. Thus, according to the QL curve, the optimal MGL is estimated to be in the range of 0.2 to 0.4, which provides a guide for parameter selection in the other applications. Moreover, the choice of order p also has impact on the thinning results. Figure 14a proposes the quantitative indexes of road extraction results with different orders of p in the case of study site 1. When order p = 1, all the fitting curves are straight lines which do not fit in the case of a winding road and thus it is not an optimal value. Generally, when the order p increases, the thining process tends to overfiting which may lead to poor robustness and accuracy for the estimation of tangential extension directions at endpoints. As illustrated in Figure 14b, with the increase of order p, the root-mean-square error (RMSE) in polynomial fitting process decreases but the tangential directions at endpoints tend to deviate from the actual direction. The inaccuracy of tangential directions may bring false network optimization which explains the slow downward trends of quantitative indexes when the order p is higher than 3 in Figure 14a. As a consequence, the value of p should not be a large number in practice. Moreover, the choice of order p also has impact on the thinning results. Figure 14a proposes the quantitative indexes of road extraction results with different orders of p in the case of study site 1. When order = 1 p , all the fitting curves are straight lines which do not fit in the case of a winding road and thus it is not an optimal value. Generally, when the order p increases, the thining process tends to overfiting which may lead to poor robustness and accuracy for the estimation of tangential extension directions at endpoints. As illustrated in Figure 14b, with the increase of order p , the rootmean-square error (RMSE) in polynomial fitting process decreases but the tangential directions at endpoints tend to deviate from the actual direction. The inaccuracy of tangential directions may bring false network optimization which explains the slow downward trends of quantitative indexes when the order p is higher than 3 in Figure 14a. As a consequence, the value of p should not be a large number in practice.

Conclusions
In this article, a road network extraction method from high-resolution SAR images is proposed. The major contributions of this work lie in constructing the road network using smooth cross-linked curves with determined functions. This mathematical description of road segments is useful for road regularization, network optimization and even data fusion. Based on this idea, the multiplicative Duda operation is introduced to achieved line feature responses and their corresponding directions. To decrease the false rate, the non-road detection using backscattering coefficient and coefficient of variation, and the filter of morphological profiles using path openings, are then presented to obtain the road candidates. Next, binary image decomposition and polynomial curve fitting are proposed

Conclusions
In this article, a road network extraction method from high-resolution SAR images is proposed. The major contributions of this work lie in constructing the road network using smooth cross-linked curves with determined functions. This mathematical description of road segments is useful for road regularization, network optimization and even data fusion. Based on this idea, the multiplicative Duda operation is introduced to achieved line feature responses and their corresponding directions. To decrease the false rate, the non-road detection using backscattering coefficient and coefficient of variation, and the filter of morphological profiles using path openings, are then presented to obtain the road candidates. Next, binary image decomposition and polynomial curve fitting are proposed to linearize the road segments. Finally, the road network is obtained by overlap, continuity, and junction optimization using the defined geometric constrain conditions. The experimental results and quantitative comparisons on three different SAR images demonstrated the effectiveness of the proposed method. Besides, the proposed road thinning method showed better performance than that attained by the traditional morphology thinning and tensor voting. Though the presented method provides satisfactory compromise results in completeness and correctness, there is still some room for the improvement of road network extraction from the high-resolution SAR images. Future work will be carried on road information fusion from multi-temporal and multi-sensor images.