# Road Extraction from High-Resolution SAR Images via Automatic Local Detecting and Human-Guided Global Tracking

Authors: Jianghua Cheng; Wenxia Ding; Xishu Ku; Jixiang Sun
Journal: International Journal of Antennas and Propagation
 (2012)

---

## Abstract

Because of existence of various kinds of disturbances, layover effects, and shadowing, it is difficult to extract road from high-resolution SAR images. A new road center-point searching method is proposed by two alternant steps: local detection and global tracking. In local detection step, double window model is set, which consists of the outer fixed square window and the inner rotary rectangular one. The outer window is used to obtain the local road direction by using orientation histogram, based on the fact that the surrounding objects always range along with roads. The inner window rotates its orientation in accordance with the result of local road direction calculation and searches the center points of a road segment. In global tracking step, particle filter of variable-step is used to deal with the problem of tracking frequently broken by shelters along the roadside and obstacles on the road. Finally, the center-points are linked by quadratic curve fitting. In 1 m high-resolution airborne SAR image experiment, the results show that this method is effective.


---

## Body

## 1. Introduction

As a kind of microwave remote sensing system, synthetic aperture radar (SAR) data acquisition could operate during both day and night and is independent from the influence of sunlight and clouds. With the increasing availability of high resolution (HR) aerial and satellite SAR sensors (such as TerraSAR-X, AIRSAR, and RadarSAT-2), huge images are gotten everyday, which promote the urgent demand for analyzing and interpreting these data. Not only in photogrammetry and remote sensing field but also in computer vision and image understanding field; automatic target detection [1] and discrimination [2] from HR SAR images is a hot subject of research. As the typical man-made object in SAR images, road is one of the most important parts of modern transportation systems, which has important geographical, political, economic, and military values. Currently, road information are gaining popularity as a reference tool in locating important civilian and military targets and have an extensive application such as city planning, transportation controlling, GIS database updating, and assistant navigation. Road network extraction from SAR images is in the growing emphasis by researchers around the world.In low and medium resolution (>9 m) SAR images, road can be viewed as a line or curve with certain orientation, whose gray value is lower than other surrounding objects. Compared with low and medium resolution ones, high-resolution (<5 m) SAR images contain more details of the road. Road may be more precisely modeled as dark elongated areas surrounded by bright edges, which are due to double-bounce reflections by surrounding buildings or uniform backscattering by the vegetation [3, 4]. However, it is more difficult to extract road network from high-resolution SAR images, because the magnitude of noise is also magnified. The existence of various kinds of disturbances, layover effects, and shadowing are more vivid than in low and medium resolution SAR images.Depending on whether a human being operator interacts with the process or not, it is generally acknowledged that road network extraction methods can be usually divided into two categories: automatic and semiautomatic. Different approaches have been proposed in the literature to address the issue of automatic road extraction from high-resolution SAR images. Negri et al. [4] proposed a novel street candidate detection method, which was further specified into a multiscale feature fusion detector and a segment extractor, followed by optimization using a Markov random field description of the road network. Lisini et al. [5] thought that classical methods based on line detection do not use all the information available and presented a road extraction method which was based on the fusion of classification (statistical information) and line detection (structural information). Then, a modified Markovian approach was used for solving the global network optimization problem. Hedman et al. [6] combined two different approaches for automatic road extraction: one is more successful in rural areas and the other one is explicitly designed for urban areas. In order to get the best combination of both, a rapid mapping filter for discriminating rural and urban scenes was utilized. Finally, the results were fused on a feature level and connected by means of a network optimization. Amberg et al. [7] introduced context information (such as buildings, vehicles, etc.) detection result in the road extraction process to avoid false alarms in confusion areas and overcome difficulties of layover and bright objects.However, viewed from the progress of current research, even if automatic road extraction from optical remotely sensed images is still difficult, let alone the SAR images which are relatively harder with the visual senses compared with the optical ones. Existing automatic road extraction algorithms cannot achieve satisfactory results in robustness, accuracy, and applicability and often requires setting lots of parameters and manual processing later, which does not give a great deal of enhancing efficiency in SAR images interpretation. Some researchers believe that the automatic road extraction method is not realistic in the foreseeable future. Semiautomatic methods can combine the rapid calculation of machine with human interpretation technique effectively. Therefore, the study of semiautomatic methods of manual intervention for quick and accurate extraction of the road is a more practical choice.Numerous semiautomatic road extraction methods have been proposed from high-resolution remote sensing images [8–14], which are usually categorized into iterative two phases: local detection and global tracking. The local detection contents mainly include direction, width, and central point position in current local road region. Global tracking is an iterative line or center-point growing process, used to recover the network topology by connecting local detected outcomes and suppress false positives. Kim et al. [8] defined a rectangle window around a user-given input point, and then by matching the window along the orientation of the road under consideration. New match proceeded by shifting a matched target window further along road orientation once matching succeeded. By repeating the process above, a series of points were obtained. Zhou et al. [9] introduced a human-computer interaction system for robust and efficient road tracking. Extended Kalman filters and particle filters are used in conjunction with human inputs to estimate road axis points and update the tracking algorithms. Lin et al. [10] proposed a new combination strategy to extract the road networks. One or more suitable road trackers were designed or combined for each type of roads based on the initialization information, such as the interlaced window matching, the profile matching, the rectangular window matching, and the parallelepiped angular texture signature (PATS).The universal problems of these methods exist in that: (1) more than two human inputted initial seed points are needed to calculate the orientation and width of the road and (2) the tracking result is easy effected by the initial seed points (3) which are generally used in optical remote sensing images, seldom reports on their applicability on high resolution SAR images. Compared with optical remote sensing images, the signal-to-noise ratio (SNR) is lower in SAR images, which are disturbed by multiplicative speckle noise and often blocked by fences, crossovers, cars, and so forth on the road. The tracking methods mentioned above would be broken frequently and cannot be directly used in high resolution SAR images.This paper presents a semiautomatic road center-point tracking method from high resolution SAR images, which consists of two iterative steps: (1) local double window detection model and (2) global particle filter tracking. Local double window detection only needs one input seed point, which has no impact on tracking. Global particle filter tracking can effectively jump over obstacles. Also, calculation efficiency is improved because of less human-computer interaction. Then, quadratic curve fitting is used to link the extracted road center-points into line. The remaining sections of this paper are organized as follows. In Section2, a general scheme of road centerpoint tracking is introduced, and the local double window detection model is proposed. In Section 3, experiments are described and the performances of the road trackers are discussed. In Section 4, some conclusions are reached.


## 2. Methodology

### 2.1. Road Characteristics

Ideally, roads appear as a continuous strip region with constant width and can be modeled as dark elongated areas surrounded by bright edges (see Figure1(a)) in high resolution SAR images. Actually, due to the interference of speckle noise and surrounding objects, the road shows the following characteristics in high-resolution SAR images.(i)
Obstacles such as iron fences and vehicles would destroy the continuity of the border of road and contaminate it with a large number of bright points or small bright areas. Road surfaces may not have adequate contrast with adjacent areas. Sometimes even human beings cannot distinguish road areas from other object areas (see Figure1(b)).(ii)
Road conglutinating with parking areas or covered by buildings shadow would make original bright bilateral edges change into single edge (see Figure1(c)). This makes the width inconsistent in the same road segment.(iii)
High buildings shelter and trees shadow, make parts of road segment even have no edge (see Figure1(d)).Road edges in high resolution SAR images. (a) Distinct bilateral edges. (b) Blurry road edges. (c) Single edge. (d) No edge.
(a)
(b)
(c)
(d)As the actual complex circumstances of high-resolution SAR images, it is difficult to use a classical SAR edge detection method to extract road.


### 2.2. Framework

The framework of our semiautomatic road extraction algorithm consists of two components (see Figure2).(1)
User. The user is the master of the whole progress of the algorithm. He (or she) inputs an initial seed point when new high resolution SAR image provided, or an orientation leading point when the tracking cannot be carried on.(2)
Computer. The computer extracts road center lines by three steps: local detection, particle filter tracking, and center-points fitting. Double window model is used to get the local road information such as orientation, width, and center point in local detection step. Particle filter tracking is used to predict and update the next road center point. When searching obstacle is met (such as end of road, border of image, tracking departure road area, crossroad, etc.), it hands over the control to the user.Figure 2
Semiautomatic road centerline extraction framework.


### 2.3. Local Detection

The type of window sometimes plays a very important role in road detection. Several windows have been used to detect local road information from high resolution remote sensing images, such as circular window [11], rectangular window [8, 12] and interlaced window [10]. Cheng et al. [11] proposed fixed size circular window to search the actual road cent-point perpendicular to the road direction (see Figure 3(a)). This method needs not to rotate in window matching calculation. Rectangular window is frequently used in road center-point searching, orientation calculation, and window matching, because of its simpleness and practicality. Figure 3(b) shows road centerlines tracking by rectangular window least squares correlation matching. Figure 3(c) shows several rectangular windows rotating around a road centerline point, in order to compute the road direction. Interlaced window is composed of two parts: some cross-section profiles (i.e., each is a typical intensity profile perpendicular and symmetrical to the road axis) and some rectangular windows of road markings (i.e., some intensity rectangles whose width is as wide as lane markings). Lin et al. [13] used interlaced window to predict next most possible position of the road axis.Local detection window. (a) Circle window. (b) Rectangular window. (c) Rotating rectangular window. (d) Interlaced window.
(a)
(b)
(c)
(d)However, the windows mentioned above cannot adapt to different width road detection because of fixed size. This paper presents a novel double window type, which is composed of two windows: outer square window and inner rectangular window (see Figure4(a)). The size of outer square window is fixed and chosen based on the resolution of the SAR sensor. Additionally, the outer square window cannot rotate anymore. It is used to calculate the local road direction. The inner rectangular window can shift (see Figure 4(b)), expand (see Figure 4(c)), and rotate (see Figure 4(d)) to adapt to different width and direction of road detection. It is used to search the center points and calculate the width of a road.Double window. (a) Definition. (b) Inner window shifting. (c) Inner window expanding. (d) Inner window rotating.
(a)
(b)
(c)
(d)


#### 2.3.1.  Outer Square Window Orientation Estimation

After an initial seed center point is given by a human being, the outer square window size and position are fixed. And the main orientation is estimated within the outer square window. There exist several orientation estimation methods, such as edge direction histogram (EDH) [15], histograms of oriented gradients (HOG) [16], and tensor voting [17]. These methods cannot be directly used in SAR images because of the influence of multiplicative speckle noise. We firstly introduce the speckle reduction process before orientation estimation.Many speckle reduction filters have been proposed in the literature for processing SAR images corrupted by multiplicative noise [18–20]. Among them, Frost filter [18] possesses better edge preserving properties and has shown superior performances for processing of speckle images.Frost et al. [18] has proposed an optimum minimum mean-square error (MMSE) estimation starting from the following observed SAR image model:
(1)I(x,y)=(r(x,y)·n(x,y))⊗h(x,y),
where h(x,y) is the system impulse response, ⊗ is the convolution operation, r(x,y) is the “ideal SAR image,” and n(x,y) is the multiplicative speckle noise. The impulse response m(x,y) of the MMSE filter that provides an estimate of  r(x,y) from  I(x,y) is obtained by minimizing the mean square error ε given by
(2)ε2=E[(r(x,y)-I(x,y)⊗m(x,y))2].
The minimization leads to the MMSE filter with the following impulse response:
(3)m(x,y)=k2(σII¯)2exp[-k(σII¯)2|D(x,y)|],
where k is a scale constant, I¯ and σI2 are the mean and variance of the observed SAR image, respectively, and D(x,y) is the distance from current point (x,y) to local center point.After speckle reduction by Frost filter, canny operator [21] is used to detect road edges. Let I(x,y) denote the SAR image, and G(x,y) is the 2-D Gaussian function. We have the following road edge intensity and orientation formulation, respectively:
(4)C(x,y)=|∇G⊗I(x,y)|,n→=∇G⊗I(x,y)|∇G⊗I(x,y)|,
where ∇G is the gradient of G, and ⊗ stands for the convolution operation.Then, the histogram of every pixel orientation is carried out within the outer square window. The maximum of the orientation histogram is the main road direction. Figures5 and 6 show the main road orientation estimation process within 100 × 100 the outer square window. Figures 5(d) and 6(d) are the orientation histogram of road edge.Orientation estimation 1. (a) Cropped SAR image. (b) Frost filter. (c) Canny edge detection. (d) Edge orientation histogram.
(a)
(b)
(c)
(d)Orientation estimation 2. (a) Cropped SAR image. (b) Frost filter. (c) Canny edge detection. (d) Edge orientation histogram.
(a)
(b)
(c)
(d)


#### 2.3.2. Inner Rectangular Window Road Width and Center-Point Calculation

The center-point position of a road is determined by its width. In previous semiautomatic road tracking methods, the road width was typically confirmed by two [9] or three [10] human inputted seed points. Whereas in our system, the road width is calculated automatically by inner rectangular window shifting, expanding, or rotating. The width and center-point calculation algorithm is detailed as follows.Step 1 (initialization).
Take the main direction of outer square window as the inner rectangular window rotating orientation, and take the outer square window’s center point as the initial inner rectangular window’s center point. Define the fixed length of the inner rectangular window asL. Define the initial width of the inner rectangular window as W, which is changeable.Step 2 (shifting).
Perpendicular to the main road direction, and shift the inner rectangular window with fixed step lengthd1 within the outer square window area, and compute the variations and means of series shifted inner rectangular windows.Step 3 (expanding).
Calculate the minimum variation and mean among the series shifted inner rectangular windows, and take corresponding inner rectangular window as the initial expanding window. Perpendicular to the main road direction, and expand the inner rectangular window with fixed step lengthd2 within the outer square window area, and compute the variations and means of series expanded inner rectangular windows.Step 4 (result).
Calculate the minimum variation and mean among the series shifted inner rectangular windows, and take corresponding inner rectangular window as the best fitted window. Take its width as the local road width, and take its center point as the local road center point.
Our road width and center-point calculation algorithm only need one initial seed input.


### 2.4. Particle Filter Tracking

Particle filtering, specifically thecondensation algorithm proposed by Isard et al. [22], is a recursive procedure to estimate the parameters of a dynamic system and has been successfully used in modeling nonlinear and non-Gaussian road tracking system [9, 23, 24]. It starts from a seed point  p0 of the road and then moves along its direction.


#### 2.4.1. State and Observation Model

The current stateSk is relative to the previous state and updated by the following nonlinear equation:
(5)sk=[xkykθk]=[xk-1+dk-1·sin(θk-1)+uk1yk-1+dk-1·cos(θk-1)+uk2θk-1+uk3],
where (xk,yk) is the current road center point, dk is the current step length, θk is the current road direction, and uki(i=1,2,3) is the noise came from state transform. The tracking procedure is represented by an ordered state sequence S0:k=(S0,…,Sk).The current observation vector is represented byZk as follows:
(6)zk=h(sk)+vk=[100010001]sk+vk,
where vk is the measure noise. The observation sequence Zk is formed from an ordered observed vectors Z0:k=(Z0,…,Zk).


#### 2.4.2. Particle Filtering

To grow such a sequence, particle filtering approximates the posterior densityp(Sk/Zk) by the particle set {Ski,wki,i=1,…,N} in each time step k, where wki is a weight used to characterize the probability of the particle Ski, N is the number of particles which is set by experience and selected as 100, and ∑i=1Nwki=1. The posterior density function (pdf) can be approximated by
(7)p(SkZ1:k)≈∑i=1Nwkiδ(Sk-Ski),
where δ is the Dirac delta function. When we have enough samples from a pdf, the set of samples can be used to closely represent that distribution. However, the pdf of p(Sk/Z1:k) is unknown in reality, and we can only generate samples from importance density function q(sk/z1:k) [24]. The choice of q(sk/z1:k) is an important issue in the design of the particle filtering. Generally, the transitional prior probability density function p(sk/sk-1i) is chosen as the importance density function. Based on this choice, the samples and weights update equation can be defined, respectively, as
(8)ski~p(sksk-1i),(9)ωki=ωk-1i×p(zkski).We assume that the observation is normally distributed, and the likelihood of the observation is(10)p(zkski)∝12πexp(-li22σ2),
where  li is the Euclidean distance between the position of particle Ski and its corresponding observation.The procedure of particle filter tracking is described as follows.Step 1 (initialization).
s
0
i
~
p
(
s
0
/
z
0
)
=
p
(
s
0
) and  w0i=1/N, where i=1,…,N.Step 2 (prediction).
SampleSki from the importance density function p(sk/sk-1i), and predict the next state vector Sk+1 by using (5).Step 3 (update).
Calculate unnormalized weightwki by using (9) and (10), and then normalize weight wki~=wki/∑i=1Nwki.Step 4 (estimate output).
S
k
^
=
∑
i
=
1
N
w
k
i
S
k
i.Step 5 (resample).
CalculateN^eff=(∑i=1N(ω~ki)2)-1. If N^eff<Nth,   {Ski,wki}=RESAMPLE[{Ski,wki}].Step 6 (return).
Go to Step2.


#### 2.4.3. Starting Point Selection and Stopping Criteria

Before road center-point tracking procedure, the initial seed point is inputted by a human being. Then, the outer square window is set based on this inputted center point, and local main road direction is calculated. Subsequently, the inner rectangular window rotates to the local main road direction and searches the actual road center point. And then, particle filter tracking is carried out.In real applications, the road characteristics are more complex in high resolution SAR images. The road center-point tracking procedure is broken frequently because of shelters along the road or obstacles on the road. For example, the outer square window cannot get the proper orientation because of shelters along the road, and the road width changes abruptly due to obstacles on the road. In this case, our particle filter tracking method makes another state update based on previous state, using a changeable step lengthdk. When the tracking procedure goes continually, the step length dk is a constant. When broken is met, dk is increased to jump over shelters or obstacles. When contiguous jumps occur, the particle filter tracking recognizes that there is a tracking failure and returns control back to the human being. Then, the human being restarts the tracking procedure by inputting another initial seed point.


### 2.5. Center-Points Fitting

Our particle filter tracking method can detect most center points despite small shelters or obstacles existing. When interference seriously exists, the extracted center points would deviate from the actual road central points. In order to lower the influence of interference, a method is applied to generate smooth road center lines by fitting all the extracted center points along one road segment.Generally, the variety of a road centerline curvature is very small. A conic is sufficient to describe the road centerline. Define the extracted center points of a road are{xi,yi}i=1Nr. The conic is y(x)=a0+a1x+a2x2. The error of fitting is
(11)E=∑i=1Nr(y(xi)-yi)2.The parameters ofa0,   a1,  a2 are calculated by means of least square error.


## 2.1. Road Characteristics

Ideally, roads appear as a continuous strip region with constant width and can be modeled as dark elongated areas surrounded by bright edges (see Figure1(a)) in high resolution SAR images. Actually, due to the interference of speckle noise and surrounding objects, the road shows the following characteristics in high-resolution SAR images.(i)
Obstacles such as iron fences and vehicles would destroy the continuity of the border of road and contaminate it with a large number of bright points or small bright areas. Road surfaces may not have adequate contrast with adjacent areas. Sometimes even human beings cannot distinguish road areas from other object areas (see Figure1(b)).(ii)
Road conglutinating with parking areas or covered by buildings shadow would make original bright bilateral edges change into single edge (see Figure1(c)). This makes the width inconsistent in the same road segment.(iii)
High buildings shelter and trees shadow, make parts of road segment even have no edge (see Figure1(d)).Road edges in high resolution SAR images. (a) Distinct bilateral edges. (b) Blurry road edges. (c) Single edge. (d) No edge.
(a)
(b)
(c)
(d)As the actual complex circumstances of high-resolution SAR images, it is difficult to use a classical SAR edge detection method to extract road.


## 2.2. Framework

The framework of our semiautomatic road extraction algorithm consists of two components (see Figure2).(1)
User. The user is the master of the whole progress of the algorithm. He (or she) inputs an initial seed point when new high resolution SAR image provided, or an orientation leading point when the tracking cannot be carried on.(2)
Computer. The computer extracts road center lines by three steps: local detection, particle filter tracking, and center-points fitting. Double window model is used to get the local road information such as orientation, width, and center point in local detection step. Particle filter tracking is used to predict and update the next road center point. When searching obstacle is met (such as end of road, border of image, tracking departure road area, crossroad, etc.), it hands over the control to the user.Figure 2
Semiautomatic road centerline extraction framework.


## 2.3. Local Detection

The type of window sometimes plays a very important role in road detection. Several windows have been used to detect local road information from high resolution remote sensing images, such as circular window [11], rectangular window [8, 12] and interlaced window [10]. Cheng et al. [11] proposed fixed size circular window to search the actual road cent-point perpendicular to the road direction (see Figure 3(a)). This method needs not to rotate in window matching calculation. Rectangular window is frequently used in road center-point searching, orientation calculation, and window matching, because of its simpleness and practicality. Figure 3(b) shows road centerlines tracking by rectangular window least squares correlation matching. Figure 3(c) shows several rectangular windows rotating around a road centerline point, in order to compute the road direction. Interlaced window is composed of two parts: some cross-section profiles (i.e., each is a typical intensity profile perpendicular and symmetrical to the road axis) and some rectangular windows of road markings (i.e., some intensity rectangles whose width is as wide as lane markings). Lin et al. [13] used interlaced window to predict next most possible position of the road axis.Local detection window. (a) Circle window. (b) Rectangular window. (c) Rotating rectangular window. (d) Interlaced window.
(a)
(b)
(c)
(d)However, the windows mentioned above cannot adapt to different width road detection because of fixed size. This paper presents a novel double window type, which is composed of two windows: outer square window and inner rectangular window (see Figure4(a)). The size of outer square window is fixed and chosen based on the resolution of the SAR sensor. Additionally, the outer square window cannot rotate anymore. It is used to calculate the local road direction. The inner rectangular window can shift (see Figure 4(b)), expand (see Figure 4(c)), and rotate (see Figure 4(d)) to adapt to different width and direction of road detection. It is used to search the center points and calculate the width of a road.Double window. (a) Definition. (b) Inner window shifting. (c) Inner window expanding. (d) Inner window rotating.
(a)
(b)
(c)
(d)


### 2.3.1.  Outer Square Window Orientation Estimation

After an initial seed center point is given by a human being, the outer square window size and position are fixed. And the main orientation is estimated within the outer square window. There exist several orientation estimation methods, such as edge direction histogram (EDH) [15], histograms of oriented gradients (HOG) [16], and tensor voting [17]. These methods cannot be directly used in SAR images because of the influence of multiplicative speckle noise. We firstly introduce the speckle reduction process before orientation estimation.Many speckle reduction filters have been proposed in the literature for processing SAR images corrupted by multiplicative noise [18–20]. Among them, Frost filter [18] possesses better edge preserving properties and has shown superior performances for processing of speckle images.Frost et al. [18] has proposed an optimum minimum mean-square error (MMSE) estimation starting from the following observed SAR image model:
(1)I(x,y)=(r(x,y)·n(x,y))⊗h(x,y),
where h(x,y) is the system impulse response, ⊗ is the convolution operation, r(x,y) is the “ideal SAR image,” and n(x,y) is the multiplicative speckle noise. The impulse response m(x,y) of the MMSE filter that provides an estimate of  r(x,y) from  I(x,y) is obtained by minimizing the mean square error ε given by
(2)ε2=E[(r(x,y)-I(x,y)⊗m(x,y))2].
The minimization leads to the MMSE filter with the following impulse response:
(3)m(x,y)=k2(σII¯)2exp[-k(σII¯)2|D(x,y)|],
where k is a scale constant, I¯ and σI2 are the mean and variance of the observed SAR image, respectively, and D(x,y) is the distance from current point (x,y) to local center point.After speckle reduction by Frost filter, canny operator [21] is used to detect road edges. Let I(x,y) denote the SAR image, and G(x,y) is the 2-D Gaussian function. We have the following road edge intensity and orientation formulation, respectively:
(4)C(x,y)=|∇G⊗I(x,y)|,n→=∇G⊗I(x,y)|∇G⊗I(x,y)|,
where ∇G is the gradient of G, and ⊗ stands for the convolution operation.Then, the histogram of every pixel orientation is carried out within the outer square window. The maximum of the orientation histogram is the main road direction. Figures5 and 6 show the main road orientation estimation process within 100 × 100 the outer square window. Figures 5(d) and 6(d) are the orientation histogram of road edge.Orientation estimation 1. (a) Cropped SAR image. (b) Frost filter. (c) Canny edge detection. (d) Edge orientation histogram.
(a)
(b)
(c)
(d)Orientation estimation 2. (a) Cropped SAR image. (b) Frost filter. (c) Canny edge detection. (d) Edge orientation histogram.
(a)
(b)
(c)
(d)


### 2.3.2. Inner Rectangular Window Road Width and Center-Point Calculation

The center-point position of a road is determined by its width. In previous semiautomatic road tracking methods, the road width was typically confirmed by two [9] or three [10] human inputted seed points. Whereas in our system, the road width is calculated automatically by inner rectangular window shifting, expanding, or rotating. The width and center-point calculation algorithm is detailed as follows.Step 1 (initialization).
Take the main direction of outer square window as the inner rectangular window rotating orientation, and take the outer square window’s center point as the initial inner rectangular window’s center point. Define the fixed length of the inner rectangular window asL. Define the initial width of the inner rectangular window as W, which is changeable.Step 2 (shifting).
Perpendicular to the main road direction, and shift the inner rectangular window with fixed step lengthd1 within the outer square window area, and compute the variations and means of series shifted inner rectangular windows.Step 3 (expanding).
Calculate the minimum variation and mean among the series shifted inner rectangular windows, and take corresponding inner rectangular window as the initial expanding window. Perpendicular to the main road direction, and expand the inner rectangular window with fixed step lengthd2 within the outer square window area, and compute the variations and means of series expanded inner rectangular windows.Step 4 (result).
Calculate the minimum variation and mean among the series shifted inner rectangular windows, and take corresponding inner rectangular window as the best fitted window. Take its width as the local road width, and take its center point as the local road center point.
Our road width and center-point calculation algorithm only need one initial seed input.


## 2.3.1.  Outer Square Window Orientation Estimation

After an initial seed center point is given by a human being, the outer square window size and position are fixed. And the main orientation is estimated within the outer square window. There exist several orientation estimation methods, such as edge direction histogram (EDH) [15], histograms of oriented gradients (HOG) [16], and tensor voting [17]. These methods cannot be directly used in SAR images because of the influence of multiplicative speckle noise. We firstly introduce the speckle reduction process before orientation estimation.Many speckle reduction filters have been proposed in the literature for processing SAR images corrupted by multiplicative noise [18–20]. Among them, Frost filter [18] possesses better edge preserving properties and has shown superior performances for processing of speckle images.Frost et al. [18] has proposed an optimum minimum mean-square error (MMSE) estimation starting from the following observed SAR image model:
(1)I(x,y)=(r(x,y)·n(x,y))⊗h(x,y),
where h(x,y) is the system impulse response, ⊗ is the convolution operation, r(x,y) is the “ideal SAR image,” and n(x,y) is the multiplicative speckle noise. The impulse response m(x,y) of the MMSE filter that provides an estimate of  r(x,y) from  I(x,y) is obtained by minimizing the mean square error ε given by
(2)ε2=E[(r(x,y)-I(x,y)⊗m(x,y))2].
The minimization leads to the MMSE filter with the following impulse response:
(3)m(x,y)=k2(σII¯)2exp[-k(σII¯)2|D(x,y)|],
where k is a scale constant, I¯ and σI2 are the mean and variance of the observed SAR image, respectively, and D(x,y) is the distance from current point (x,y) to local center point.After speckle reduction by Frost filter, canny operator [21] is used to detect road edges. Let I(x,y) denote the SAR image, and G(x,y) is the 2-D Gaussian function. We have the following road edge intensity and orientation formulation, respectively:
(4)C(x,y)=|∇G⊗I(x,y)|,n→=∇G⊗I(x,y)|∇G⊗I(x,y)|,
where ∇G is the gradient of G, and ⊗ stands for the convolution operation.Then, the histogram of every pixel orientation is carried out within the outer square window. The maximum of the orientation histogram is the main road direction. Figures5 and 6 show the main road orientation estimation process within 100 × 100 the outer square window. Figures 5(d) and 6(d) are the orientation histogram of road edge.Orientation estimation 1. (a) Cropped SAR image. (b) Frost filter. (c) Canny edge detection. (d) Edge orientation histogram.
(a)
(b)
(c)
(d)Orientation estimation 2. (a) Cropped SAR image. (b) Frost filter. (c) Canny edge detection. (d) Edge orientation histogram.
(a)
(b)
(c)
(d)


## 2.3.2. Inner Rectangular Window Road Width and Center-Point Calculation

The center-point position of a road is determined by its width. In previous semiautomatic road tracking methods, the road width was typically confirmed by two [9] or three [10] human inputted seed points. Whereas in our system, the road width is calculated automatically by inner rectangular window shifting, expanding, or rotating. The width and center-point calculation algorithm is detailed as follows.Step 1 (initialization).
Take the main direction of outer square window as the inner rectangular window rotating orientation, and take the outer square window’s center point as the initial inner rectangular window’s center point. Define the fixed length of the inner rectangular window asL. Define the initial width of the inner rectangular window as W, which is changeable.Step 2 (shifting).
Perpendicular to the main road direction, and shift the inner rectangular window with fixed step lengthd1 within the outer square window area, and compute the variations and means of series shifted inner rectangular windows.Step 3 (expanding).
Calculate the minimum variation and mean among the series shifted inner rectangular windows, and take corresponding inner rectangular window as the initial expanding window. Perpendicular to the main road direction, and expand the inner rectangular window with fixed step lengthd2 within the outer square window area, and compute the variations and means of series expanded inner rectangular windows.Step 4 (result).
Calculate the minimum variation and mean among the series shifted inner rectangular windows, and take corresponding inner rectangular window as the best fitted window. Take its width as the local road width, and take its center point as the local road center point.
Our road width and center-point calculation algorithm only need one initial seed input.


## 2.4. Particle Filter Tracking

Particle filtering, specifically thecondensation algorithm proposed by Isard et al. [22], is a recursive procedure to estimate the parameters of a dynamic system and has been successfully used in modeling nonlinear and non-Gaussian road tracking system [9, 23, 24]. It starts from a seed point  p0 of the road and then moves along its direction.


### 2.4.1. State and Observation Model

The current stateSk is relative to the previous state and updated by the following nonlinear equation:
(5)sk=[xkykθk]=[xk-1+dk-1·sin(θk-1)+uk1yk-1+dk-1·cos(θk-1)+uk2θk-1+uk3],
where (xk,yk) is the current road center point, dk is the current step length, θk is the current road direction, and uki(i=1,2,3) is the noise came from state transform. The tracking procedure is represented by an ordered state sequence S0:k=(S0,…,Sk).The current observation vector is represented byZk as follows:
(6)zk=h(sk)+vk=[100010001]sk+vk,
where vk is the measure noise. The observation sequence Zk is formed from an ordered observed vectors Z0:k=(Z0,…,Zk).


### 2.4.2. Particle Filtering

To grow such a sequence, particle filtering approximates the posterior densityp(Sk/Zk) by the particle set {Ski,wki,i=1,…,N} in each time step k, where wki is a weight used to characterize the probability of the particle Ski, N is the number of particles which is set by experience and selected as 100, and ∑i=1Nwki=1. The posterior density function (pdf) can be approximated by
(7)p(SkZ1:k)≈∑i=1Nwkiδ(Sk-Ski),
where δ is the Dirac delta function. When we have enough samples from a pdf, the set of samples can be used to closely represent that distribution. However, the pdf of p(Sk/Z1:k) is unknown in reality, and we can only generate samples from importance density function q(sk/z1:k) [24]. The choice of q(sk/z1:k) is an important issue in the design of the particle filtering. Generally, the transitional prior probability density function p(sk/sk-1i) is chosen as the importance density function. Based on this choice, the samples and weights update equation can be defined, respectively, as
(8)ski~p(sksk-1i),(9)ωki=ωk-1i×p(zkski).We assume that the observation is normally distributed, and the likelihood of the observation is(10)p(zkski)∝12πexp(-li22σ2),
where  li is the Euclidean distance between the position of particle Ski and its corresponding observation.The procedure of particle filter tracking is described as follows.Step 1 (initialization).
s
0
i
~
p
(
s
0
/
z
0
)
=
p
(
s
0
) and  w0i=1/N, where i=1,…,N.Step 2 (prediction).
SampleSki from the importance density function p(sk/sk-1i), and predict the next state vector Sk+1 by using (5).Step 3 (update).
Calculate unnormalized weightwki by using (9) and (10), and then normalize weight wki~=wki/∑i=1Nwki.Step 4 (estimate output).
S
k
^
=
∑
i
=
1
N
w
k
i
S
k
i.Step 5 (resample).
CalculateN^eff=(∑i=1N(ω~ki)2)-1. If N^eff<Nth,   {Ski,wki}=RESAMPLE[{Ski,wki}].Step 6 (return).
Go to Step2.


### 2.4.3. Starting Point Selection and Stopping Criteria

Before road center-point tracking procedure, the initial seed point is inputted by a human being. Then, the outer square window is set based on this inputted center point, and local main road direction is calculated. Subsequently, the inner rectangular window rotates to the local main road direction and searches the actual road center point. And then, particle filter tracking is carried out.In real applications, the road characteristics are more complex in high resolution SAR images. The road center-point tracking procedure is broken frequently because of shelters along the road or obstacles on the road. For example, the outer square window cannot get the proper orientation because of shelters along the road, and the road width changes abruptly due to obstacles on the road. In this case, our particle filter tracking method makes another state update based on previous state, using a changeable step lengthdk. When the tracking procedure goes continually, the step length dk is a constant. When broken is met, dk is increased to jump over shelters or obstacles. When contiguous jumps occur, the particle filter tracking recognizes that there is a tracking failure and returns control back to the human being. Then, the human being restarts the tracking procedure by inputting another initial seed point.


## 2.4.1. State and Observation Model

The current stateSk is relative to the previous state and updated by the following nonlinear equation:
(5)sk=[xkykθk]=[xk-1+dk-1·sin(θk-1)+uk1yk-1+dk-1·cos(θk-1)+uk2θk-1+uk3],
where (xk,yk) is the current road center point, dk is the current step length, θk is the current road direction, and uki(i=1,2,3) is the noise came from state transform. The tracking procedure is represented by an ordered state sequence S0:k=(S0,…,Sk).The current observation vector is represented byZk as follows:
(6)zk=h(sk)+vk=[100010001]sk+vk,
where vk is the measure noise. The observation sequence Zk is formed from an ordered observed vectors Z0:k=(Z0,…,Zk).


## 2.4.2. Particle Filtering

To grow such a sequence, particle filtering approximates the posterior densityp(Sk/Zk) by the particle set {Ski,wki,i=1,…,N} in each time step k, where wki is a weight used to characterize the probability of the particle Ski, N is the number of particles which is set by experience and selected as 100, and ∑i=1Nwki=1. The posterior density function (pdf) can be approximated by
(7)p(SkZ1:k)≈∑i=1Nwkiδ(Sk-Ski),
where δ is the Dirac delta function. When we have enough samples from a pdf, the set of samples can be used to closely represent that distribution. However, the pdf of p(Sk/Z1:k) is unknown in reality, and we can only generate samples from importance density function q(sk/z1:k) [24]. The choice of q(sk/z1:k) is an important issue in the design of the particle filtering. Generally, the transitional prior probability density function p(sk/sk-1i) is chosen as the importance density function. Based on this choice, the samples and weights update equation can be defined, respectively, as
(8)ski~p(sksk-1i),(9)ωki=ωk-1i×p(zkski).We assume that the observation is normally distributed, and the likelihood of the observation is(10)p(zkski)∝12πexp(-li22σ2),
where  li is the Euclidean distance between the position of particle Ski and its corresponding observation.The procedure of particle filter tracking is described as follows.Step 1 (initialization).
s
0
i
~
p
(
s
0
/
z
0
)
=
p
(
s
0
) and  w0i=1/N, where i=1,…,N.Step 2 (prediction).
SampleSki from the importance density function p(sk/sk-1i), and predict the next state vector Sk+1 by using (5).Step 3 (update).
Calculate unnormalized weightwki by using (9) and (10), and then normalize weight wki~=wki/∑i=1Nwki.Step 4 (estimate output).
S
k
^
=
∑
i
=
1
N
w
k
i
S
k
i.Step 5 (resample).
CalculateN^eff=(∑i=1N(ω~ki)2)-1. If N^eff<Nth,   {Ski,wki}=RESAMPLE[{Ski,wki}].Step 6 (return).
Go to Step2.


## 2.4.3. Starting Point Selection and Stopping Criteria

Before road center-point tracking procedure, the initial seed point is inputted by a human being. Then, the outer square window is set based on this inputted center point, and local main road direction is calculated. Subsequently, the inner rectangular window rotates to the local main road direction and searches the actual road center point. And then, particle filter tracking is carried out.In real applications, the road characteristics are more complex in high resolution SAR images. The road center-point tracking procedure is broken frequently because of shelters along the road or obstacles on the road. For example, the outer square window cannot get the proper orientation because of shelters along the road, and the road width changes abruptly due to obstacles on the road. In this case, our particle filter tracking method makes another state update based on previous state, using a changeable step lengthdk. When the tracking procedure goes continually, the step length dk is a constant. When broken is met, dk is increased to jump over shelters or obstacles. When contiguous jumps occur, the particle filter tracking recognizes that there is a tracking failure and returns control back to the human being. Then, the human being restarts the tracking procedure by inputting another initial seed point.


## 2.5. Center-Points Fitting

Our particle filter tracking method can detect most center points despite small shelters or obstacles existing. When interference seriously exists, the extracted center points would deviate from the actual road central points. In order to lower the influence of interference, a method is applied to generate smooth road center lines by fitting all the extracted center points along one road segment.Generally, the variety of a road centerline curvature is very small. A conic is sufficient to describe the road centerline. Define the extracted center points of a road are{xi,yi}i=1Nr. The conic is y(x)=a0+a1x+a2x2. The error of fitting is
(11)E=∑i=1Nr(y(xi)-yi)2.The parameters ofa0,   a1,  a2 are calculated by means of least square error.


## 3. Experimental Results and Analysis

To evaluate our proposed method, three different scale high resolution SAR images are chosen, which are taken by airborne in the region of Hefei, Anhui province, China. Figure7 is zoomed in from the image with a ground pixel size of about 1 m. Figures 8 and 9 are normal scale with a ground pixel size of about 1 m. All the experiments are run on P(R) dual-core 2.8 GHz CPU, with 2 GB SDRAM, and the software platform is MATLAB R2007b.Initial seed point automatic update. (a) Input seed point. (b) The procedure of update. (c) Updated center-point.
(a)
(b)
(c)Figure 8
Road center-points extraction methods comparison 1.Figure 9
Road center-points extraction methods comparison 2.Figure7(a) shows that an initial seed road center point is inputted by human being, which has a little bias from the true one. Red squares are outer window, and yellow rectangles are inner window. Our double windows can adjust automatically to the real position. In order to explore the effectiveness of our proposed method, the comparisons with other center-points extraction method are made like Kim’s least squares correlation matching [8], which are shown in Figures 8 and 9. Yellow points are road center points extracted by our tracking method and connected by white lines. Green points are road center points extracted by the Kim’s least squares correlation matching method and connected by green lines. In Figure 8, there exist shelters (Windows 3) along the roadside. Our tracking strategy can jump over it and come back to the normal road trajectory (Windows 5). But the center points extracted by the Kim’s matching method are deviated from the normal road trajectory. In Figure 9, strong obstacles (Windows 2 and Windows 3) are on the road. Our tracking strategy can jump over these obstacles in the same way. However, the center points extracted by the Kim’s matching method are deviated from the normal road trajectory too.Figure10 shows a large SAR scene also taken by airborne in the region of Hefei. Yellow points are road center points extracted by our tracking method and connected by white lines. Pink points are road center-points alter fitting and connected by pink lines. In Figure 10, we can see that the quadratic curve fitting method can link the discrete road center points into a line and smooth the abrupt changes. It is worth mentioning that our tracking method also can jump over the crossroad.Figure 10
Extracted main road center-points fitting.


## 4. Conclusion

We have proposed a new double windows tracking method for extracting center points and center lines. This algorithm consists of two iterative steps, namely, local detection and global tracking. The local detection step uses double windows, which can estimate the local road main direction and search the deviated center point automatically. The global tracking step makes use of particle filter, which can jump over shelters, obstacles, and crossroads. The presented method has been tested on three high resolution SAR images with 1 m spatial resolution. The experimental results validate the accuracy and efficiency of our method. The road center points are extracted by only using a single seed point to initialize the tracking. Moreover, this algorithm remains robust regardless of abrupt changes introduced by disturbs.



---

*Source: 989823-2012-12-13.xml*