EXTRACTION OF ROAD NETWORK USING A MODIFIED ACTIVE CONTOUR APPROACH

Said Mssedi1,2, Mohamed Ben Salah2, Riadh Abdelfattah3,4 and Amar Mitiche12

1 Ecole Polytechnique de Tunis, 2 Institut national de la recherche scientifique – Canada3 URISA, École Supérieure des Télécoms, Université de Carthage – Tunisie

4 Institut Telecom ; Telecom Bretagne – France.

ABSTRACT

Road network extraction from digitized map consists in par-titioning a map into two different classes (road and back-ground) distinguished with their intensity values and geomet-ric shape. In this paper, we propose a semi-automatic methodwhere user intervention is minimal. We adopt a variationalframework where a minimization energy is proposed. Roadsare tracked by active contours which evolve according to evo-lution equations derived from the energy minimization. Weconsider that the energy is composed of three main termswhich have never been combined and used in such applica-tion. The first term is the geodesic active contour (GAC)which is contour-based and which evolves toward high gra-dients in the image. The second term matches the intensitydistribution inside the contour with a model distribution. Thethird term is a geometric constraint which describes the thinstructure of routes. The user intervention is restricted to learn-ing the road intensity model.

Index Terms— Road network, variational framework, ac-tive contour, minimization energy.

1. INTRODUCTION

Road extraction from digitized maps is a fundamental prob-lem which has major practical importance for GIS update [1].The classic manual techniques are time consuming. Besides,manual techniques are known to be dependent on the operator[2]. Thus, the solution may change even for the same mapand is prone to human errors. Proposed solutions for roadextraction can be classified among maps features extractionmethods. However, all these methods almost use the samebasic features which are intensity, geometry and topology.The intensity information is used, especially, in color im-ages where roads have different color profile [3]. Geometrycharacteristics help the extraction process in the sense that,for example, roads are known to be elongated ribbons ratherthan other forms. However, the fact that roads are continueslines which form a network is rather topological information.Bases on these features, many automatic and semi-automaticmethods have been published [4].Regarding the amount of user intervention, the published

methods are often classified into automatic and semi-automaticmethods. A fully automatic method is theoretically supposedto take, as input, any kind of maps and provides, as output,the extracted road network. However, road characteristicsvary significantly with map type, road type, image resolu-tion, surrounding objects, presence of occluding objects, etc.Hence, a fully automatic method, as defined above, is quiteimpossible to achieve [3]. Moreover, the published automaticmethods are known to work for a given range of images withspecific conditions [5, 3].

In this paper, we propose a semi-automatic method whereuser intervention is minimal. We adopt a variational frame-work where a minimization energy is proposed. This frame-work allows embedding various constraints in the energy and,thus, solving many kinds of problems and various types ofmaps. Roads are tracked by active contours which evolve ac-cording to evolution equations derived from the minimizationenergy. We consider that the energy is composed of threemain terms which have never been combined and used insuch application. The first term is the geodesic active con-tour (GAC) [6] which is contour- based and which evolvestoward high gradients in the image. The second term matchesthe intensity distribution inside the contour with a model dis-tribution [7]. This term has been used mainly for tracking,but never joined to the GAC model for segmentation purposes[3, 2]. The third term is a geometric constraint relative to thethin structure of routes, and which was developed to satisfyto our need. Hence, the contour evolution is guided locallyby the GAC term, and globally by the distribution matchingterm and take finally the form of the geometric term. The userintervention is restricted to learning the road intensity model.Indeed, this task is done only one time for a given map, thenthe same model is used for maps of the same class. How-ever, other methods can be used for this task as clustering, orgraph-cut methods.The remainder of this paper is organized as follows: Section2 describes the minimization energy terms, and explains thelearning process for the distribution matching term. In section3 we details the developed constraint proposed to describe ge-ometrically the road network. Section 4 reports experimentalresults, and Section 5 contains concluding remarks.

2011 18th IEEE International Conference on Image Processing

978-1-4577-1303-3/11/$26.00 ©2011 IEEE 2929

(a) (b)

Fig. 1. The position of the active contour after some iterationsof detection of these two objects: (a) original image; (b) thespace of level set translated with colors.

2. GEODESIC TERM AND BHATTACHARYYAMEASURE

Segmentation by active contour is based essentially on theoptimization of a functional, called also energy. The energy,which embeds constraints on the required solution, is min-imized in a way that evolves the contour towards image ob-jects of interest. Let I : Ω→ Rn be the image function wheren = 1 for gray-scale images and n = 3 for color images, andlet γ(s, t) : [0, 1]×R+ → Ω be a parameterized planar curvewhere Ω is the image domain. In our application, the intensityenergy functional is composed of two essential terms:

EI = E1 + αE2, (1)

where α > 0 is a real constant. E1 represents the local energyterm which guides the active contour towards areas of highgradient and E2 is the global energy which guides the activecontour to areas which have a similar distribution to a modellearned beforehand.

2.1. Local Energy

The first term of energy (1), the local term, is dependent onthe contour γ(s, t) as follow:

E1 =

∫γ

g(I)ds, (2)

where g is the edge indicator function. g is a monotonicallydecreasing function which is often referred to as stoppingfunction. A common choice of this function is:

g(I) =1

1 + |∇Gσ ∗ I|p, p = 1, 2,

where ∇Gσ ∗ I is a smooth version of the image I , and Gσis generally the Gaussian kernel Gσ = 1√

σexp(−|x

2+y2|4σ ).

The edge indicator function which depends solely on imagegradient, stops the active contour on the edges of objects. In-deed, it vanishes on objects borders because image gradientis high. The curve flow minimizing this energy is known asthe geodesic active contour (GAC) [6]. Although, efficient in

many applications where objects are very well characterizedby image gradient transitions, the GAC term is insufficient inour case. For instance, all road maps we consider contain var-ious objects other than roads. Although unwanted, these ob-jects attract the GAC because they have high image gradientsalong their boundaries. Hence, another term which biases theevolving contour towards only the object of interest is needed.

2.2. Global Energy

Let z be the photometric variable of interest. For example,z could be an intensity, color vector, or texture vector. LetY : Ω → Z be a mapping from the image domain to thespace of the photometric variable and x ∈ Ω. Thus, if Z isthe space of intensities, then Y (x) is just a gray-scale image;if Z is the space of color vectors, then Y (x) is a color im-age. Denote by R ⊂ Ω the region surrounded by the contouri.e. γ = ∂R. The sample distribution within region R of thevariable z is estimated as follow: d(z,R) =

∫RK(z−Y (x))dx

A(R) ,

where A(R) =∫Rdx is the area of the region R, and K is

a kernel function. Most common choices of K are the Diracdistribution and the Gaussian kernel function. The purpose ofthis part is to track the region R where the probability densityd(z,R) most closely matches a model distribution M learnedbeforehand. To do so, we measure the similarity between bothdistributions by the Bhattacharyya measure. However, othermeasures can be used such as the Kullback-Leibler distance[7]. This energy term is global because it references the inten-sity distribution in a given region and, as such, is not writtenin a pixel wise form. Hence, the global energy term is:

E2 = −B(d,M) = −∫Z

√d(z,R)M(z)dz. (3)

This measure varies between 0 and 1, where 0 indicates com-plete mismatch, and 1, a complete match. In the next section,we derive evolution flows from the energies defined above.

2.3. Evolution Flows

We derive, in this part, the evolution equations of the activecontour γ minimizing EI in (1). These equations manage themovement of the active contour from its initial position to theborders of all roads in the map road network. Using functionalderivative calculus, we obtain: FI = ∂γ

∂t = −∂EI

∂γ . Thus, thetotal evolution flow is expressed as follows:

FI =∂γ

∂t= −∂EI

∂γ= −∂E1

∂γ− α∂E2

∂γ= F1 + αF2, (4)

where F1 is the local force, called also the geodesic flow,which derives from the local energy and F2 is the global force,we refer to it here by the Bhattacharyya flow, which derivesfrom global energy: F1 = (gk − ∇g.∇−→n )−→n , k and −→n arerespectively the curvature and the normal to the contour γ. ;F2 = 1

2A(R)

[M1/2(Y (γ))

d1/2(Y (γ))−B(m, d)

]−→n .

2011 18th IEEE International Conference on Image Processing

2930

2.4. Level set implementation

Since its emergence, level set method has been widelyadopted to implement active contour based methods [8].It makes no assumption about the topology of the objects.We use here the fact that, for a given evolution equation ofthe form ∂−→γ

∂t = β−→n , where β is scalar function, the level setform is as follows ∂φ

∂t = β |∇φ|. The total level set flow isthus derived (FT = ∂φ

∂t ).

FT =

[1

2A

(M

12

d12

−B(m, d)

)+

(g.k −∇g. ∇φ

|∇φ|

)]|∇φ|

(5)

3. GEOMETRIC TERM

The extraction of road network from digitized map cannotbe reduced only to the recognition of objects boundaries androutes intensity distribution, but needs to specify their geom-etry. Take for example the presence of different objects onthe digitized map which have the same intensity distributionas the road network as shown in Fig.2 (a). For this reason,adding to the intensity energy (1) a term which describes thegeometric shape of the roads is necessary.

ET = EI + βE3 = E1 + αE2 + βE3, (6)

(a) (b)

Fig. 2. first test image: (a) digitized map, (b) result of theextraction with the developed method.

3.1. The distance between pixels and the active contour

In our application, the level set is managed by the distancefunction between pixel and the zero level of the level set.Therefore, the value of level set at each point x of the im-age Ω is the distance between x and the contour γ. To distin-guish between the region inside and outside the contour andto keep the concept of distance for all points of level set, weadded a minus sign for any measure outside the contour γ.So, the distance function will be signed in our application:d(x, γ) = d(x, γ) if x ∈ R and −d(x, γ) if x ∈ Rc. Ris area inside the contour γ and Rc is the complement of R

in the image domain Ω. So, we remark that the maximumdistance between pixels and the active contour corresponds tothe centers of thick and convex objects. Since roads are struc-turally thin and not convex, the distance function can be a toolwhich distinguishes roads from thick and convex objects suchas circles.

3.2. The developed distance function

This mathematical tool is used in the first step of segmen-tation with the calculation of gradient and measurement ofBhattacharyya. So after determining the initial position ofour active contour, the level set is generated automaticallyby the geodesic term, Bhattacharyya measure and distancefunction. However, after several iterations, the level set nolonger obeys to the distance function. Indeed several terms(geodesic, Bhattacharyya, smoothness) manage the new val-ues of level set at each iteration. So it is impossible that thelinear combination of these terms keep the distance notion inthe structure of the level set. For this reason, the first chal-lenge of this new method is to approach after each iterationthe level set to the distance function.Li et al [8] resolved the first issue. For more accurate com-putation involving the level set function, they need to regu-larize the level set function by penalizing its deviation froma signed distance function. They add to the overall energy aterm which can be characterized by the following function:

P (φ) =

∫Ω

1

2(|∇φ (x)| − 1)

2dx. (7)

Retaining the concept of distance in the level set at each iter-ation, we added an additional term to the total energy to de-scribe the structure of thin roads. This geometric term appearsin the overall energy in the form of a Heaviside function:

E3 = H(φ− 1

2m

), (8)

wherem is the maximum width of the road network. The usersets value m in the initialization. Using the Heaviside func-tion hides some objectives. Firstly, the functionH

(φ− 1

2m)

takes the value 1 for all points distant more than 12m from

the border γ and 0 otherwise. So the minimization of E3

closes the curve until the distance between all pixels in regionR is less than 1

2m (the middle of the widest road in the net-work). Secondly, H

(φ− 1

2m)

is discontinuous in a criticalvalue which represents the geometric threshold between largeand small objects. This discontinuity can significantly differ-entiate the two types of objects and avoid topological changeof curve in uncertainty area. Indeed, if we use a continuousfunction, accumulation of small changes in pixels close to thecritical value can change incorrectly the position of γ in thisarea. Finally, Dirac function δ, the derivative of the Heavi-side function is very simple in the implementation level. The

2011 18th IEEE International Conference on Image Processing

2931

derivation of E3 gives a new term in the evolution equation:

F3 = δ

(φ− 1

2m

)(9)

This new method extracts only the thin structure which resem-bles to roads independently of their shapes (straight, curved,...) and their position (far, close to large objects). Fig.1 whichshows the position of the contour γ few iterations after theextraction of objects having the same distribution of roads,explains the influence of F3(9) in evolution of active contour.This new term divides large objects into three parts. It elim-inates, at the distance of 1

2m from contour γ, a very narrowarea of the region R. After penetration of contour γ in theselarge objects, the geodesic term looks again for new areas ofhigh gradient until the complete disappearance of these ob-jects in the region R.

(a) (b)

Fig. 3. Second test image: (a) digitized map, (b) result of theextraction with the developed method.

4. EXPERIMENTAL RESULTS

To demonstrate the efficiency of the developed road networkextraction approach, we present some experimental resultstested on the digital maps depicted by Fig. 2 (a), Fig. 3 (a)and Fig. 4 (a). The result presented in Fig. 2 (b) containssome objects on the background with the same intensity dis-tribution as roads. We also tested the developed approach onmaps with road which are very curved and very close (stuck)with different widths which vary between 10 and 30 pixels(Fig. 3 (b)). Despite these difficulties, the new method wasable to extract the majority of road network with few errorsmainly due to the narrow width of the road as the differencebetween Fig. 3 (a) and (b). Moreover, considering the roadnetwork of a titled digitized map (Fig. 4 (a)), the extractionresult (Fig. 4 (b)) is perfect. This is despite the roads aresometimes cut to display their names or other information.These obtained results show the rigidity of this method whichextract accurately the entire road network and avoid two im-portant problems: Firstly, the extraction of road network onlyeven though it has the same intensity distribution with otherobjects. Secondly, the continuity of the road network eventhough it is cut by words and information panels.

(a) (b)

Fig. 4. Third test image: (a) digitized map, (b) result of theextraction with the developed method.

5. CONCLUSION

In this paper, we presented a method for extracting roadnetworks based on image segmentation by active contoursand level set. This method gathers the contour informationfrom the geodesic model, region information from distribu-tion matching model and geometric information from thedistance function. It improves the efficiency of road net-work extraction by the rapid detection of borders throughthe geodesic term, matches the segmented area with a modelwhich contains accurate information about roads through theBhattacharyya measure and distinguishes, by a geometricterm, objects which have the same intensity distribution asroads and differ in geometric level.

6. REFERENCES

[1] J.B. Mena, “State of the art on automatic road extractionfor GIS update: A Novel Classification,” Pattern Recog-nition Letters, vol: 24, No. 16, pp. 3037-3058, 2003.

[2] L. Bibo, K. Yongxia, H. Wang, “Automatic Road Ex-traction Method Based on Level Set and Shape Analy-sis,” ICICTA, vol:3, pp. 511-514, 2009.

[3] M. Rochery, I. H. Jermyn, J. Zerubia, “Gap closure inroad networks using higher-order active contours.” ICIP,Vol: 3, pp. 1879-1882, 2004.

[4] T. Chan and L. Vese, “Active Contours without Edges.”IEEE Trans. on Ima. Proc., Vol:10, pp. 266-277, 2001.

[5] C. R. Silva and J. A. S. Centeno, “Automatic extractionof main roads using aerial images.” Pattern Recognitionand Image Analysis, Vol: 20, No. 2, pp. 225-233, 2010.

[6] V. Caselles, R. Kimmel, and G. Sapiro, “Geodesic activecontours,” in Proc. of IEEE ICCV, pp. 694-699, 1995.

[7] D. Freedman and T. Zhang. “Active contours for track-ing distributions.” IEEE Trans. IP, 13(4): 518-526, 2004.

[8] C. Li, C. Xu, C. Gui, M. D. Fox, “Level set evolu-tion without re-initialization: A new variational formu-lation”, In IEEE, CVPR, vol:1, pp 430–436 (2005).

2011 18th IEEE International Conference on Image Processing

2932