Segmentation of Urban Areas Using Road Networks ∗

Nicholas Jing Yuan Yu Zheng Xing Xie

Microsoft Research Asia Microsoft Research Asia Microsoft Research Asia
[email protected] [email protected] [email protected]

ABSTRACT 10 5 km
Region-based analysis is fundamental and crucial in many geospatial-
related applications and research themes, such as trajectory analy-
sis, human mobility study and urban planning. In this paper, we
report on an image-processing-based approach to segment urban ar-
eas into regions by road networks. Here, each segmented region is
bounded by the high-level road segments, covering some neighbor-
hoods and low-level streets. Typically, road segments are classified
into different levels (e.g., highways and expressways are usually
high-level roads), providing us with a more natural and semantic
segmentation of urban spaces than the grid-based partition method.
We show that through simple morphological operators, an urban
road network can be efficiently segmented into regions. In addi-
tion, we present a case study in trajectory mining to demonstrate
the usability of the proposed segmentation method.
Figure 1: Level 0,1,2-road segments in the urban area of Beijing
Categories and Subject Descriptors
H.2.8 [Database Management]: spatial databases and GIS, data
Intrinsically, this segmentation is more natural than grid-based seg-
mining.
mentation due to the following reasons: First, people live in these
roads-segmented regions and POIs (points of interests) fall in these
General Terms regions instead of the main roads. Second, the road-segmented re-
gions as the origin and destination of a trip are the root cause of
Algorithms
human mobility. In short, people travel among the road-segmented
regions.
1. INTRODUCTION In a Geographical Information System (GIS), there are two ma-
In many geospatial-related applications, such as trip planning[18], jor models to represent spatial data: vector-based model and raster-
urban planning/urban computing[6][25] and traffic analysis [24][15], based model. Vector-based model uses geometric primitives such as
an urban area is often segmented into sub-regions for in-depth anal- points, lines and polygons to represent spatial objects referenced by
ysis or complexity reduction. Intuitively, a digital map can be seg- Cartesian coordinates, while raster-based model quantizes an area
mented into equal-sized grids, where each grid is a rectangle, as into small discrete grid-cells (cuboids for the 3D spatial objects)
shown in Figure 2(a). While the grid-based approach offers simplic- indexing all the spatial objects. For example, the vector model of
ity in implementation, it has deficiencies as well: The segmented Beijing road network (shown in Figure 1) stores a road segment as a
grids do not have semantic meanings, thus the results of the analysis polyline (polygone for a circuit road), where a polyline is consist of
based on these grids do not provide us with a natural understanding, a sequence of shape points, represented by coordinates. A road seg-
with respect to the road network. ment r is associated with a direction symbol (r.dir) and two termi-
Typically, the road segments of a road network are categorizes in- nal points (r.s, r.e). If r.dir=one-way, r can only be traveled from
to different levels according to their functions, e.g., Figure 1 presents r.s to r.e, otherwise, people can start from both terminal points,
the road network of Beijing with Level-0,1,2 road segments. Here, i.e., r.s → r.e or r.e → r.s. Instead, the raster model rasterizes the
the Level-0 and Level-1 road segments (red-colored) are freeways road network into grid-cells, where each road segment is represent-
and city expressways in Beijing, and Level-2 road segments (blue- ed by a sequence of cell-IDs. Figure 3(a) and Figure 3(a) reveal a
colored) are urban arterial roads. These main roads naturally seg- portion of the Beijing road network stored with vector model and
ment the urban area into sub-regions with varying sizes and shapes. raster model(726×726 grid-cells) respectively.
Both of the two models have advantages and disadvantages de-
∗
This technical report details the process of map segmentation uti- pending on the specific applications. For instance, vector-based
lized in our previous papers[15],[25] and [23] method is more powerful for precisely finding shortest-paths, where-
as it requires intensive computation when performing topological
analysis, such as the problem of map simplification [10], which
 (a) Grid-based (b) Hierarchy-based (c) Morphology-based

Figure 2: Different map segmentation methods

drivers by constructing a spatio-temporal profitability map with a

grid-based segmentation, where the probabilities are calculated by
the historical data. Compared with grid-based segmentation, our
solution considering the high-level roads as the boundary are more
natural for studying the human mobility on a map, as described in
Section 1.
Gonzalez et al. [11] propose a novel road network partition ap-
proach based on the road hierarchy. Specifically, the road networks
are first divided into areas by high level roads, then the partition pro-
cess is recursively performed for each area. The partition process
is implemented by finding the strongly connected components after
(a) Vector model (b) Raster model (726×726) the removal of the intersection nodes connected to high level road-
s as well as the terminals of high level road segments themselves.
Figure 3: A portion of Beijing road network printed from vec- Figure 2(b) presents the results of this approach for a portion of
tor and raster models Beijing road network. However, this method does not work in our
scenario since 1) our desired region bounded by high level road seg-
ments may contain several strongly connected components, and 2)
is proved to be NP-complete [9]. On the other hand, raster-based we aim to segment the whole area instead of just the road nodes
model is more computational efficient and succinct for territorial into regions, i.e., we need a mapping from any locations represent-
analysis, but the accuracy is limited by the number of cells used for ed by latitudes and longitudes (within the bounding box of the road
discretizing the road networks. network) to the region IDs.
In our method, we choose the raster-based model to represent Morphology operators are widely used in Geographical informa-
the road network and utilize morphological image processing tech- tion systems as well as image processing. Saradjian and Amini [20]
niques to address the task of map segmentation. Specifically, a employs mathematical morphology for map simplification from re-
raster-based map is regarded as a binary image (e.g., 0 stands for mote sensing images by extracting skeletons from the image and
road segments and 1 stands for blank space). In order to remove the convert the structure into vectors. Similar works are presented in
unnecessary details, such as the lanes of a road and the overpass- [16],[17] and [8]. Different from the above methods which are
es, we first perform a dilation operation to thicken the roads. As a based on remote sensing images, we aim to segment the urban area
result, we can fill the small holes and smooth out the unnecessary represented by vector-based model into regions, instead of simpli-
details. Second, we obtain the skeleton of the road networks by per- fying the map or extracting structures from the map.
forming a thinning operation. This operation recovers the size of a
region which was reduced by the dilation operation, while keep-
ing the connectivity between regions. The last step is to perform a 3. MAP SEGMENTATION
connected component labeling (CCL) that finds individual regions
We display the vector-based road network in a plane by perform-
by clustering “1”-labeled grids. We detail the above operations in
ing map projection [22], which transforms the surface of a shpere
Section 3 and provide a case study in Section 4.
(i.e., the Earth) into a 2D plane (We use Mercator projection in our
approach [1, 4]). Then we convert the vector-based road network
2. RELATED WORK into the raster model by gridding the projected map [5]. Intuitively,
Grid-based map segmentation is extensively used in geospatial- each pixel of the projected map image can be regarded as a grid-cell
related analysis. Krumm and Horvitz [13] predict the drivers’ des- of the raster map. Consequently, the road network is converted to a
tination by mapping the historical trips into grid-cells and learning binary image, e.g., 1 stands for the road segments (termed as fore-
the destination probabilities with respect to each cell. Powell et al. ground) and 0 stands for the blank areas (termed as background).
[18] propose an approach to suggest the most profit grid for taxi This section introduces an mathematical morphological approach
 (a) Before dilation (b) After dilation (a) After 7 iterations (b) Until convergence

Figure 4: Dilation operator Figure 5: Thinning operator

for segmenting the binary-image-based road network into region- Here, we employ the subfields-based parallel thinning algorithm
s. In general, a morphological operator calculates the output image proposed in [12]. This method first divide the binary image space
given the input binary image and a structure element, whose size into two disjointed subfields in a checkerboard pattern. Each itera-
and shape are pre-defined. tion is consist of two sub-iterations in these two subfields:
3.1 Dilation • In the first sub-iteration, check every pixel p in the first sub-
Dilation is one of the basic morphological operators. Let A be a field, delete p iff Condition 1, 2 and 3 are all satisfied.
binary image and B be the structure element, the dilation of A by
B is defined as: • In the second sub-iteration, check every pixel p in the second
[ subfield, delte p iff Condition 1, 2 and 4 are all satisfied.
A⊕B = Ab , (1)
b∈B C ONDITION 1. XH (p) = 1, where
where Ab = {a + b|a ∈ A}, i.e., the translation of A by the vector 4
X
b. The dilation operator is commutative, which can also be defined XH (p) = bi
by: i=1
[ (
A⊕B =B⊕A= Ba . (2) 1 if x2i−1 = 0 and (x2i = 1 or x2i+1 = 1)
a∈A
bi =
0 otherwise
Actually, for any p in A, after the dilation, p = 1 iff the intersection
between A and B, centred at p, is not empty. C ONDITION 2. 2 ≤ min{n1 (p), n2 (p)} ≤ 3, where
The purpose of the dilation operation is to remove the unneces- 4
sary details for map segmentation, avoiding the small connected ar- X
n1 (p) = x2k−1 ∨ x2k
eas induced by these unnecessary details such as bridges and lanes.
k=1
Figure 4(a) plots a portion of road network before the dilation op- 4
erator. As is shown, the small holes between the lanes and viaducts X
n2 (p) = x2k ∨ x2k+1
are filled, where we use a 3×3 matrix with all values one as the
k=1
structure element B.
C ONDITION 3. (x2 ∨ x3 ∨ x̄8 ) ∧ x1 = 0
3.2 Thinning
As a consequence of the previous dilation operator, the road seg- C ONDITION 4. (x6 ∨ x7 ∨ x̄4 ) ∧ x5 = 0
ments are turgidly thickened. In this step, we aim to extract the
skeleton of the road segments while keeping the topology structure Figure 5(a) and Figure 5(b) are the results after 7 iterations and
(such as the Euler number) of the original binary image. Specif- until convergence (no pixel are deleted any more) respectively.
ically, the thinning operator is performed to remove certain fore-
ground pixels from the input binary image. For a given pixel in the 3.3 Connected Component Labeling
input image, whether it should be removed depends on its neigh- The connected component labeling operation finds the connect-
bouring pixels. Lam et al. [14] present an in-depth survey of exist- ed 0 pixels (the blank area) in the binary image, after the thinning
ing thinning algorithms, most of which are implemented in a recur- operation. For a given pixel x, the neighbouring 4 pixels shown in
sive manner. According to the way an algorithm check the pixels, Figure 6(a) are called the 4-neighbours of x. Similarly, the 8 neigh-
the thinning algorithms can be categorized into parallel thinning al- bouring pixels shown in Figure 6(b) are called the 8-neighbours
gorithms and sequential thinning algorithms. The sequential algo- of x. We call the sequence y1 , y2 , . . . , yn an 8-path (4-path), if
rithms check each pixel according to a fixed order in each iteration, ∀i = 1, 2, . . . , n − 1, yi+1 is an 8-neighbour (4-neighbour) of yi .
and the deletion of pixel p in the nth iteration is not only related We say a region Q in a binary image is 8-connected (4-connected)
to the results in the n − 1th iteration, but also relate to the pixels iff all the pixels in Q have the same value and for any two pixels
examined before p in the nth iteration; In contrast, for the paral- in Q, there exist an 8-path (4-path) connecting the two pixels. Note
lel algorithms, the deletion of p only depends on the results in the that in the previous thinning operation, connectivity paradox may
n − 1th iteration, therefore, can be performed independently in a be induced if we keep the same type of connectivity (4-connected
parallel manner. or 8-connected) for both the foreground and the background[19].
 O O
x2 x5 x4 x3 (b)

R1 p1
x3 x x1 x6 x x2 R3
p2 (c)
p4 p5
x4 x7 x8 x1 p3
p5
(a) 4-neighbours (b) 8-neighbours R2 R4 (d)
D R D
5
Figure 6: The 4-neighbours and 8-neighbours of x (a)

Figure 8: Trajectory mapping

we can determine which region it belongs to (by examining the re-

gion ID of this GPS point’s occupied pixel). However, in many
applications related to trajectories [23][15], instead of mapping iso-
lated points, we need to transform a whole trajectory into a region
ID sequence. For example, T : S → p1 → p2 · · · → p5 → D
in Figure 8 is a trajectory originated from region R1 to region R5 .
Note that some points may located along the roads (such as p1 )
while some points actually enter the non-road regions (such as p5 ).
(a) Local result (b) Global result In order to obtain a meaningful mapping from the raw trajectory
to the region IDs, for each point x of the trajectory, we apply the
Figure 7: Segmented regions after connected component label- following rules:
ing
• If x is the start/destination point, we first retrieve its k nearest
neighbouring pixels N (x). Then we assign x with the most
Since it is desired for the road segments (foreground) to have unit non-zero (not on road segments) frequent label (region ID)
width, typically, we preserve the 8-connectivity of the foreground in N (x), e.g., the origin O is assigned in region R1 and the
(i.e., the 8-connectivity does not change before and after the thin- destination D is labeled in region R5 , as shown in Figure 8
ning process for the road segments) and the 4-connectivity of the (b) and 8 (d) respectively.
background [14].
There exist many algorithms for connected component labeling. • If x is not the start/destination point, if all its k nearest neigh-
Here, we employ the classical two-pass algorithm [21], the main bouring pixels N (x) have an identical label (region ID), we
steps are as follows: assign x with this label, as depicted in Figure 8 (c); otherwise,
we label x with 0 (along the road).
• Set the initial label l to 1 (note the label here has nothing to
do with the 0-1 value in the binary image). In the first pass, As a result, T is mapped to the region ID sequence R1 → 0 →
scan each pixel from left to right and from top to bottom. If a 0 → 0 → 0 → R4 → R5 .
pixel x is a backgroud pixel (with value 0), Based on the above approach, we map a large number of GPS
trajectories of taxicabs during 3 month [23] to the segmented re-
– Retrieve all the 4-neighbours of x, label x with the s- gions, and then extract the Origin-Destination (OD) pairs. Figure 9
mallest label among its 4-neighbours; if there are no presents the top-200 frequent OD pairs for the time window 8am–
labeled 4-neighbours of x, label x with the new label 10am on an average of weekdays, where a darker color indicates a
(l ← l + 1). higher frequency. Figure 10 further shows the number of transitions
– Store the equivalence of x’s label and its 4-neighbours’
labels.
3558

• Perform the second pass: scan from left to right and from top
to bottom, for each pixel x, label x with the smallest label
among all the labels in the equivalent class of x’s label. 1955

Applying this algorithm to the binary image 5(b), we obtain the 1083

segmented regions as shown in Figure 7(a). Figure 7(b) presents

the result for the whole road network of Beijing. 609

4. A CASE STUDY 351

After the connected component labeling, we obtain the segment-

ed regions as well as their labels (termed as region IDs) for the
whole map. Here, all the road segments form a region labeled as 0. Figure 9: Top-200 frequent OD pairs on weekdays during 8am
Therefore, given a GPS point represented by latitude and longitude, to 10am
 4
10
[5] http://www.mathworks.cn/help/toolbox/map/
f7-10852.html.
10
3
[6] http://research.microsoft.com/en-us/projects/
num of transitions
urbancomputing/.
[7] http://www.mathworks.cn/help/toolbox/comm/ref/
2
10
vec2mat.html.
[8] M. ANSOULT, P. SOILLE, and J. LOODTS. Mathematical
10
1
morphology- a tool for automated gis data acquisition from scanned
thematic maps. Photogrammetric engineering and remote sensing, 56
(9):1263–1271, 1990.
10
0

10
0 1
10
2
10
3
10
4
10
5
10
[9] R. Estkowski. No steiner point subdivision simplification is np-
k complete. In Proc. 10th Canadian Conf. Computational Geometry.
Citeseer, 1998.
Figure 10: Number of transitions w.r.t. k [10] R. Estkowski and J. Mitchell. Simplifying a polygonal subdivision
while keeping it simple. In Proceedings of the seventeenth annual
(OD pairs) changing over k (log-log plot), which suggests a typical symposium on Computational geometry, pages 40–49. ACM, 2001.
Zipf distribution. [11] H. Gonzalez, J. Han, X. Li, M. Myslinska, and J. P. Sondag. Adaptive
fastest path computation on a road network: a traffic mining approach.
5. CONCLUSION In Proceedings of the 33rd international conference on Very large data
bases, VLDB ’07, pages 794–805, 2007. ISBN 978-1-59593-649-3.
In this technical report, we present a morphological map segmen-
tation method based on high level road segments of urban road net- [12] Z. Guo and R. Hall. Parallel thinning with two-subiteration algorithms.
Communications of the ACM, 32(3):359–373, 1989.
works. Compared with traditional segmentation approaches, such
as the grid-based method and the hierarchy-based method, our method [13] J. Krumm and E. Horvitz. Predestination: Where do you want to go
produces a more meaningful and semantic segmentation, where the today? Computer, 40(4):105–107, 2007.
high-level road segments are utilized as natural boundaries for the [14] L. Lam, S. Lee, and C. Suen. Thinning methodologies-a comprehen-
segmented regions. sive survey. IEEE Transactions on pattern analysis and machine in-
telligence, 14(9):869–885, 1992.
APPENDIX [15] W. Liu, Y. Zheng, S. Chawla, J. Yuan, and X. Xing. Discovering
spatio-temporal causal interactions in traffic data streams. In Proceed-
A. PSEUDO-CODE ings of the 17th ACM SIGKDD international conference on Knowl-
edge discovery and data mining, pages 1010–1018. ACM, 2011.
Input: Vector-based road network consisting of high-level road
segments G [16] E. López-Ornelas. High resolution images: segmenting, extracting in-
Output: Matrix M ; formation and gis integration. World Academy of Science, Engineering
/* the numbers in the matrix indicate the and Technology, 54:172–177, 2009.
identities of segmented regions */
1 Convert G into a binary image B; [17] C. Martel, G. Flouzat, A. Souriau, and F. Safa. A morphological
/* construct a raster-based model from a method of geometric analysis of images: Application to the gravity
vector-based model */ anomalies in the indian ocean. Journal of Geophysical Research, 94
2 B ← Dilation(B); (B2):1715–1726, 1989.
/* Perform dialation */ [18] J. Powell, Y. Huang, F. Bastani, and M. Ji. Towards reducing taxicab
3 B ← Thinning(B);
cruising time using spatio-temporal profitability maps. In Proceed-
/* Perform thinning until the convergence */ ings of the 12th International Symposium on Advances in Spatial and
4 M ← ConnectedComponentLabeling(B);
Temporal Databases, SSTD ’11, 2011.
/* Consider 4-connectivity for the regions and
8-connectivity for the road segments */ [19] A. Rosenfeld and J. Pfaltz. Sequential operations in digital picture
5 return M processing. Journal of the ACM (JACM), 13(4):471–494, 1966.
[20] M. Saradjian and J. Amini. Image map simplification using mathe-
B. USEFUL LINKS matical morphology. International Archives of Photogrammetry and
For the input format of the vector-based road network and the Remote Sensing, 33:36–43, 2000.
vector-to-raster conversion, please refer to [7]; For the morpholog- [21] L. Shapiro and G. Stockman. Computer Vision. 2001. Prentice Hall,
ical operators, such as dilation and thinning, refer to [3]; For con- 2001.
nected component labeling, refer to [2].
[22] J. Snyder. Map projections–A working manual. Number 1395. USG-
PO, 1987.
References
[1] http://msdn.microsoft.com/en-us/library/ [23] J. Yuan, Y. Zheng, and X. Xie. Discovering regions of different func-
aa940991.aspx. tions in a city using human mobility and pois. In Proceedings of the
18th ACM SIGKDD international conference on Knowledge discovery
[2] http://www.mathworks.cn/help/toolbox/images/ and data mining. ACM, 2012.
ref/bwlabel.html.
[24] Y. Zheng and X. Zhou. Computing with spatial trajectories. Springer-
[3] http://www.mathworks.cn/help/toolbox/images/ Verlag New York Inc, 2011.
ref/bwmorph.html.
[25] Y. Zheng, Y. Liu, J. Yuan, and X. Xie. Urban computing with taxicabs.
[4] http://msdn.microsoft.com/en-us/library/ In Proceedings of the 13th International Conference on Ubiquitous
bb429590.aspx. Computing Ubicomp 2011.