DETECTION OF ROADS IN SAR IMAGES USING PARTICLE FILTER 1 1

Yilun Chen, 2 Qiong Yang, 1 Yuantao Gu, 1 Jian Yang

Tsinghua University

ABSTRACT A novel method is presented to detect roads in Synthetic Aperture Radar (SAR) images. A multi-segmented poly-line model is introduced to provide a more accurate description of the road as well as to ensure the road curve’s smoothness in the model level. We then solve the road detection problem using the Bayesian tracking theory, where the Particle Filtering algorithm is adopted to provide a simple and consistent framework. The effectiveness and robustness of the proposed method is demonstrated by experimental results.

2

Microsoft Research Asia

solve the tradeoff of segment’s length. Then a poly-line tracker is applied to sequentially discover the consecutive segments using Particle Filtering algorithm [11]-[13], which is a sequential Monte Carlo method based on the importance sampling and Bayesian theory. As we would see in latter section, Particle Filtering algorithm provide a consistent framework for solving the road detection problem via integrating the local line structure and global smoothness together into the tracking procedure. 2. THE TRACKING FRAMEWORK FOR ROAD DETECTION

1. INTRODUCTION Automatic or semi-automatic extraction of linear features from synthetic aperture radar (SAR) images including roads, rivers and railroads is an important research area in recent years. It has various applications, such as multi-sensor image registration, cartology, and geomorphologic studies. In the past 20 years, many approaches for road detection have been proposed [1]-[10]. For instance, dynamic programming, which is first used by Fischler et al. [1] and then applied on SAR images [2, 3, 4] is used to minimize some global cost function. In [4] ,the optimal search is used to find the minimum cost path linking straight-line segments extracted by local Hough transformation. Global optimization methods are also widely used in this area. Tupin et al. [7] proposes a nearly automatic detection algorithm based on Markov random field, Jeon et al. [10] presents a road detection technique using a genetic algorithm. Recently a Gibbs point process framework is adopted by Stoica et al. to detect road networks in remotely sensed images [8], where the simulated annealing algorithm is used for optimizing the defined energy function. Energy minimizing curves, such as the Snake method [5, 6] are another possible solutions for road detection, which are usually applied in post-processing because of the need of a series of control points. Most of the previous methods model a road curve as a multisegmented straight line and extract the road in a detect-and-link approach, i.e., the whole extraction procedure consists of two steps: first, straight-line segments are extracted in a local detection step , and then a global grouping step is introduced to selectively link the detected segments. Due to the multi-segmented straight line model, the length of each segment has to be carefully selected. If the length of each segment is too long, the real road curve may not be covered well; On the contrary, the speckle noise of SAR images may lead to an undesirable extraction result, if the segments are too short. In this paper, a novel method of road detection is proposed in the framework of Bayesian tracking theory. Instead of the traditional straight line model, a poly-line model is adopted to provide a more accurate description of the segmented road curve, therefore This work was supported by the National Important Fundamental Research Plan of China (2001CB309401), by the National Natural Science Foundation of China (60402030) and by the Fundamental Research Foundation of Tsinghua University.

The novel multi-segmented poly-line model is proposed in this paper to approximate a real road curve (Fig.1(left) ). In this model, the road is divided into small consecutive segments, which are approximated by 3-order polynomials: yk (x) = a0 + a1 x + a2 x2 + a3 x3 ,

(1)

where k is the index of the individual segment. Please notice that the coordinates system of the current segment are defined by the endpoint of the previous segment. See Fig.1(left), the origin of the kth coordinates is determined by the location of the k-1 th segment’s endpoint. In addition, the directions of new x and y axes follow the tangent direction and normal direction on the endpoint of the k-1 th segment, respectively.

Fig. 1. (Left) Multi-segmented poly-line road model. (Right) Tracking based road detection. To ensure the adjacent segments smoothly connected, the first order derivative is restricted to be zero at the joint point, yk (x) |x=0 = 0, yk0 (x) x=0 = 0. (2) Using eq.(2) in eq.(1), and denoting the kth endpoint as in the polar coordinate, we get x3 tan θk 2 x + ck ( − x2 ), (3) rk cos θk rk cos θk where ck = ak rk cos θk . It can be proved that the maximum absolute value of the second derivative on the kth segment is yk (x) =

max

0≤x≤rk cos θk

|yk00 (x) | = 2|

tan θk | + 4|ck |. rk cos θk

(4)

Therefore, ck can be regarded as a parameter controlling the smoothness of the kth segment. If rk is defined as step size and restricted to be constant r , the kth segment can be determined by sk = (θk , ck ) ,

i=1

(i) (i)

wk sk

ˆeff i. Calculate N ˆeff ≤ Nth , for i = 1, · · · , N , ii. If N P (i) (j) (j) sk ∼ N j=1 wk δ(sk − sk ), and (i)

wk = 1/N One may refer to [11] for detailed deduction and analysis. 4. THE PARTICLE FILTER BASED ROAD DETECTION METHOD As mentioned above, extraction of the road in image can be achieved via tracking solutions. In addition, particle filtering algorithm, which is derived from Bayesian theory and implemented by Monte Carlo simulation, is an efficient approach to deal with the nonlinear and non-Gaussian problems. Consequently, a new road extraction method based on particle filter is proposed in this paper. 4.1. State Model and Propagation According to the road model introduced in section 2, the state vector of the particle tracker is defined as

3. THE BAYESIAN FILTER AND ITS MONTE CARLO IMPLEMENTATION The Bayesian filter calculate p(sk |z1:k ) sequentially by the following two stages: Z p(sk |z1:k−1 ) = p(sk |sk−1 )p(sk−1 |z1:k−1 )dsk−1 (6) (7)

where α is a normalizing constant. The recurrence relations, eq.(6) and eq.(7), however, are only a conceptual solution, for the integration in whole space is intractable in practice. In the case of a linear model and the Gaussian noise, the recursive construction of the posterior distribution can be handled analytically yielding the Kalman filter. However, for most cases in image processing problems the linear Gaussian hypothesis does not hold. So the Monte Carlo implementation is a recommended solution. Particle filter is a technique for implementing a recursive Bayesian filter by Monte Carlo simulations. The key idea is to approximate the distribution via discrete random measures defined by random samples with associate weights, named n o particles. The particle set is often denoted as χ = s(i) , w(i) , where s(i) is the ith random i=1:N

PN

(d) Resample:

(5)

where sk is named as the state vector. Note that to determine a third order polynomial one need to specify four parameters. However, to guarantee the smoothness of two adjacent segment, the four parameters are constrained by the boundary condition in eq.(2). Thus, only two free parameters are needed to determine one segment. For the coordinate x, y of the k th segment is based on the k-1 th one, instead of determining the segment directly, sk controls how the kth segment “grows” from the previous one. Therefore, the problem of road extraction from an image is equivalent to sequantially estimate the mentioned state vector sk in such an image, using pixels around the supposed road region denoted as zk . From the Bayesian perspective, the tracking problem is to recursively calculate the belief of sk given the data z1:k = {z1 , ..., zk } (Fig.1(right)), i.e., to construct the p(sk |z1:k ). Based on the above description, the general Bayesian tracking framework, which will be explained briefly in the following section, can be adopted to solve the problem of road extraction.

p(sk |z1:k ) = αp(zk |sk )p(sk |z1:k−1 ),

(c) Estimate: ˆsk =

sample and w(i) is its weight. Based on the discrete approximation of the pdf p (sk |z1:k ), the complete procedure of particle filtering is described in the following paragraph.

sk = (θk , ck ) ,

(8)

where θk determines the endpoint of the kth segment and ck controls its smoothness. Based on the state model defined above, the particles are propagated by (1) (2) θk = u k , c k = u k , (9) (1)

(2)

where uk and uk are random perturbation which determine the next segment’s endpoint and smoothness, respectively. From the dynamic system view, the state equation is zero order where the smoothness is not shown explicitly. However, as mentioned in the above section, the current kth segment’s coordinate is determined by the previous one. This means that the state vector sk determines how the current segment“grows” from the previous one rather than determining the segment’s locations directly. In this way, the smoothness between adjacent segments is achieved. Furthermore, through limiting the range of ck , the smoothness within each segment can also be obtained. 4.2. Particle Weight

The criteria to calculate particle weight, which is used to measure the confidence of each sample is essential to the tracking performance of particle filter. According to the above state model, the weight of (i) (i) 1. Initialize: s0 ∼ p (s0 ) and w0 = 1/N , where i = 1, · · · , N a particle should reflect the probability how the corresponding seg2. For k = 1, 2, · · · ment behaves like a part of a road in such SAR image, i.e., the locations, and the line structure of that segment. In this paper, a modified (a) Propagate: for i = 1, · · · , N ratio of average (ROA) detector is employed and the particle weight (i) (i) sample sk from a given samping function q(sk |sk−1 , zk ) is defined based on the response of this line detector. (b) Calculate weight: The ROA detector was firstly introduced in [7, 14]. Given a i. Calculate un-normalized weights: for i = 1, · · · , N piece of segment, index 1 denotes the central region, while indexes 2 (i) (i) and 3 denote its bilateral sides, respectively (Fig.2). The correspondcalculate wk for each particle sk ing detector is defined as ii. Normalize weight: for i = 1, · · · , N (i) w µ1 µ1 (i) wk = PN k (j) r(s(i) ) = max( , ), (10) µ2 µ3 j=1 wk

k+1 k- 1

2 1 3

k

Fig. 2. The regions of the modified ROA detector, region1 denotes the road region, region 2 and 3 denote the background region. The central line denotes the road segment.

where µj is the radiometric empirical mean value of a given region j = 1, 2, 3. Since the road region is always darker than its neighbors in SAR images, the smaller the response is, the more likely the segment belongs to a road. The response of the modified ratio line detector is then mapped to the particle weight utilizing w(i) = √

r 2 (s(i) ) 1 − 2σ 2 e , 2πσ

(11)

where σ is a constant parameter. 4.3. Starting Point Selection and Stopping Rule The starting point can be pre-selected or automatically detected by moving window line/edge detector [7, 14]. To avoid false-alarms, the window of the detector should be chosen large enough and the threshold should be chosen strictly. Our method works better with human inspected starting points. For a given particle filter, the tracking process can be terminated for any of the following criteria:

Fig. 3. Road extraction result on the area of Cascade Head, US. (Left) The SAR image. (Middle) Road extraction result from the proposed method. (Right) Road extraction result from the previous method in [4].

Fig.4 gives a more detailed demonstration of the proposed method. Fig.4(left) is the original SAR image, Fig.4 shows the extraction result. To look through the proposed method, we draw the locations of all the particles during the tracking process(Fig.4(right)). We can find that the size of the area covered by particles depends on the noise condition of this area. In the low-noised area such as area 1, the particles are centralized near the real road segment; In the highnoised area such as area 2, the particles distributes dispersedly. This hints that the proposed method provides an adaptive searching region: for the area where the road is easy to see, the method diminish the searching region to reduce the false-alarm from other road-like segments; for the area where the road is hard to find, the method enlarge the searching region trying to look for a most likely road segment.

1. The maximum weight of the particles is less than a threshold wth and tracking steps are more than a threshold Lth . 2. The current segment is out of the boundary of an image. 5. EXPERIMENTAL RESULTS In this section, we applied our algorithm to images obtained from NASA L band AIRSAR, where the span image is used for road detection. The experimental results are shown in Fig.3-5, including single roads and multiple roads across various regions such as terrains, mountains and etc. In Fig.3(left), the SAR image contains a single-track road in a mountain area, the extraction results obtained by proposed method is shown in Fig.3(middle), where we can see the road is precisely detected. For comparison, we also applied another tracking based method proposed in [4], which can be regarded as a typical “detectand-link” method. This method is based on dynamic programming over straight segments extracted by local Hough transformation in each block (see Fig.3(right)). Having obtained routing sequence from start node to end node, a complete road can be detected as a smooth continuous curve after post-processing. The result is shown in Fig.3(right) where we can see that in some highly textured area, the road segment is easily mis-detected which leads to a poor performance. However, the proposed method shows its robustness over noisy regions which benefits from the constraint of smoothness between the adjacent poly-line segments. Moreover, the multi-segmented poly-line model provides a more accurate extraction result.

Fig. 4. Road extraction result on the area of Walnut Gulch, US. (Left) The SAR image. (Middle) Road extraction result from the proposed method. (Right)All of the particles during the extraction. In Fig.5 we attempt to apply the proposed method on the multiroads extraction, where four roads intersect with each other. Multiple particle filters are employed from respective starting points. Different from the single road tracker, a new criteria is added to terminate the detection processing. When the current estimated segment of one particle filter is close enough to that of another filter estimated before, the former filter terminates because the two filters are tracking the same road. After all of the trackers are terminated, the combination of results from all particle filters constructs the final multi-roads detection result. Comparing to the global searching methods [7, 10], the proposed method also shows its efficiency in computation. For a image with 1024 ×1024 pixels and tracking with 200 particles, the whole method can be carried out within 1 min on a Pentium III 1.13G pro-

road network extraction,” IEEE Trans. Geosci. Remote Sensing, vol. 36, pp. 434-453, 1998. [8] RB. Stoica, XB. Descombes and JB. Zerubia, “A Gibbs Point Process for Road Extraction from Remotely Sensed Images,” International Journal of Computer Vision, vol. 57, pp. 121 - 136, 2004. [9] M. Barzohar and D. B. Cooper, “Automatic finding of main roads in aerial images by using geometric-stochastic models and estimation,” IEEE Trans. Pattern Anal. Machine Intell., vol. 18, pp. 707-721, 1996. [10] B. Jeon, J. Jang, and K. Hong, “Road Detection in Spaceborne SAR Images Using a Genetic Algorithm,” IEEE Trans. Geosci. Remote Sensing, vol. 40, pp. 22-29, Jan. 2002 Fig. 5. Span image (left) and the road detection result (right) on the area of Pelliston, US.

[11] A. Doucet, A. Freitas, and N. Gordon, Sequential Monte Carlo Methods in Practice, New York: Springer, 2001

cessor.

[12] A. Doucet, “On Sequential Simulation-Based Methods for Bayesian Filtering,” Technical report, University of Cambridge, Dept. of Engineering, 1998 6. SUMMARY

In this paper, a particle filter based road detection method has been proposed where a road curve is modelled as multi-segmented polylines. By introducing the curve’s smoothness into the road model we reduce the parameter number of a given segment. The particle are given meaning to how a segment ”grows” from the previous one and the particle weight reflects the local line structure of the corresponding road segment. The experimental results demonstrate that our method shows robustness results even in heavy-noised regions. The Particle Filtering algorithm provides a simple and consistence framework, where we show the effectiveness of the importance sampling on adjusting the searching region adaptively. The experiments are taken base on road detection in SAR images. However, it can be readily modified to river, railway or any continuous curve detection in the images observed from various sensors. 7. REFERENCES [1] M. A. Fischler, J. M. Tenenbaum, and H. C. Wolf, “Detection of roads and linear structures in low resolution aerial imagery using a multisource knowledge integration technique,” Comput. Graph. Image Processing, vol. 15, no. 3, pp. 201-223, 1981. [2] J. W. Wood, “Line finding algorithms for SAR,” in Royal Signals and Radar Establishment, Memo. 3841, 1985. [3] R. Samadani and J. F. Vesecky, “Finding curviline structures in speckled images,” IEEE Trans. Geosci. Remote Sensing, vol. 28, pp. 669-673, July 1990. [4] Y. Rianto, S. Kondo and T. Kim, “Detection of roads from satellite images using optimal search”, Int. J. Pattern Recognit. Artif. Intell., Vol. 14, no. 8, pp. 1008-1023. Dec. 2000. [5] M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active contours models,” International Journal of Computer Vision, vol. 1, no. 4, pp. 321-331, 1988. [6] L. Bentabet, S. Jodouin, D. Ziou, J. Vaillancourt , “Road vectors update using SAR imagery: a snake-based method,” IEEE Trans. Geosci. Remote Sensing, Vol. 41, pp. 1785-1803, 2003. [7] F. Tupin, H. Maitre, J. F. Margin, J. M. Nicolars, and E. Pechersky, “Detection of line structures in sar images: Application to

[13] J.S. Liu and R. Chen, “Sequential Monte Carlo Methods for Dynamic Systems,” J. Am. Statist. Ass., 1998 [14] R. Touzi, A. Lopes, and P. Bousquet, “A statistical and geometrical edge detector for SAR images,” IEEE Trans. Geosci. Remote Sensing, vol. 26, pp. 764-773, Nov. 1988.