Technische Universität Berlin Department of Geoinformation Science Geoinformation Technology: Voronoi Diagrams Alexandra Stadler Institute for Geodesy and Geoinformation Science Technische Universität Berlin Credits: This material is mostly an english translation of the course module no. 2 (‘Geoobjekte und ihre Modellierung‘) of the open e-content platform www. geoinformation. net. WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 2 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 3 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Introduction Department of Geoinformation Science Planning of a supermarket chain profitable? How do people decide where to do their shopping? Similar question: Where do people post their letters? “Post Office Problem” (Computational Geometry). To study this problem we make the following simple assumptions: The price is the same at every site The cost of transportation to the site is a linear function of the Euclidean distance Customers try to minimize the cost of acquiring the good or service 4 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 5 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Geometric Interpretation Department of Geoinformation Science Voronoi diagram = subdivision of the total area into regions (the trading areas of the sites) People who live in the same region all go to the same site The trading area for a given site consists of all those points for which that site is closer than any other site 6 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 7 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
History Department of Geoinformation Science 1644: Descartes was first to use Voronoi diagrams – he described the regions of gravitational influence of the sun and other stars However, they are named after Georgy Fedoseevich Voronoi René Descartes 1596 -1650 studies on positive quadratic forms (together with Peter Gustav Lejeune Dirichlet) Georgy Voronoi 1868 -1908 8 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 9 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Applications in Real World Department of Geoinformation Science 10 Astronomy -- Identify clusters of stars and clusters of galaxies Biology, Ecology, Forestry -- Model and analyze plant competition ("Area potentially available to a tree", "Plant polygons") Geography -- Analyzing patterns of urban settlements Geometric Modelling -- Finding "good" triangulations of 3 D surfaces Marketing -- Modelling of market areas Mathematics -- Study of positive definite quadratic forms ("Dirichlet tesselation", "Voronoi diagram") Meteorology -- Estimate regional rainfall averages ("Thiessen polygons") Pattern Recognition -- Find simple descriptors for shapes that extract 1 D characterizations from 2 D shapes ("Medial axis" or "skeleton" of a contour) Physiology -- Analysis of capillary distribution in cross-sections of muscle tissue to compute oxygen transport ("Capillary domains") Robotics -- Path planning in the presence of obstacles Zoology -- Model and analyze the territories of animals A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Applications within the Field of Computer Science Department of Geoinformation Science All Nearest Neighbors -- Given a set of points, find each point's nearest neighbors Euclidean Minimum Spanning Tree, Convex Hull -- as parts of the Delaunay triangulation Closest Pair -- Given a set of points, which two are closest together? 11 Knuth's Post Office Problem -- Given a set of locations for post offices, how do you determine the closest post office to a given house? (ignoring zip codes) Largest Empty Circle -- also known as the Toxic Waste Dump Problem A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 12 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Definition Department of Geoinformation Science Let P={p 1, p 2, …, pn} be a set of n distinct points in the plane; these are the sites. We define the Voronoi diagram of P as the subdivision of the plane into n cells, one for each site in P, with the property that a point q lies in the cell corresponding to a site pi if and only if dist(q, pi) < dist(q, pj) for each pj є P with j≠i. Voronoi edge Voronoi node i site 13 The cells represent the trading areas of the sites pi; they are called Voronoi regions. The boundary of a Voronoi region may consist of line segments, half-lines and infinite lines, which we call Voronoi edges. Voronoi region The end points of Voronoi edges are called Voronoi nodes. They belong to the boundary of 3 or more Voronoi regions. A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 14 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Construction of the Voronoi Diagram Department of Geoinformation Science ? 15 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Construction of the Voronoi Diagram Department of Geoinformation Science Definition: Bisector of p 1 and p 2 = perpendicular bisector of the line segment p 1 p 2 Splits the plane into two half-planes, each containing the points closer to p 1 or p 2 respectively Retrieve all the points closer to p 1 than to any other point in the plane Calculate n-1 bisectors between p 1 and all the other points Intersection of n-1 half-planes (possibly unbounded) open convex polygonal region bounded by at most n-1 vertices and at most n-1 edges 16 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 17 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 18 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation I: Intersection of Half-planes Department of Geoinformation Science Algorithm I (Naïve method)… Input: list of n sites P Output: Voronoi diagram for P Step 1: For each site generate n-1 half-planes and compute their common intersection Step 2: Return all intersections as regions of the Voronoi diagram 19 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation I: Intersection of Half-planes Department of Geoinformation Science Algorithmic Complexity – Algorithm I (Naïve method)… Input: list of n sites P Output: Voronoi diagram for P naïve more sophisticated Step 1: For each site generate n-1 half-planes and compute their common intersection Step 2: Return all intersections as regions of the Voronoi diagram O(n³) O(n) O(n³) 20 A. Stadler – Geoinformation Technology: lecture 11 O(n²logn) WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 21 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Divide and Conquer… one of the fundamental paradigms for designing efficient algorithms original problem is recursively divided into several simple subproblems of almost equal size solution of the original problem is obtained by merging the solutions of the subproblems 22 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithm II. 1 (Divide and conquer method)… Input: List of sites P in increasing order of the x coordinates Output: Voronoi diagram for P Step 1: if n<=3 Direct construction of the Voronoi diagram Step 2: otherwise: Divide P into PL and PR (using the median) Compute Voronoi diagrams for PL and PR by algorithm II. 1 (recursion!) Merge the Voronoi diagrams by algorithm II. 3 Step 3: 23 Return the Voronoi diagram for P A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithm II. 1 (Divide and conquer method)… Animation: 24 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Essential problem: Merging of two Voronoi diagrams Generate new Voronoi nodes Generate new Voronoi edges result from the interaction of left and right sites Delete superfluous portions from the input Voronoi diagrams The sets of sites PL and PR can be separated by a vertical line. To construct the new nodes and edges, we want to start at the top We need an algorithmic tool for finding the topmost edge/line such that It touches two disjoint convex polygons and All other points of the polygons lie on the same side of the line. 25 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithm II. 2 (Upper common support)… Input: 2 non-overlapping convex polygons UL and UR (convex hull of PL and PR) Output: Pair of nodes such that the line L(u, w) forms the upper common support Step 1: Find the node in UL with the largest x coordinate u and the node in UR with the smallest x coordinate w Step 2: Until no more nodes are detected: While the next node on UL (ccw) lies above L(u, w): replace u While the next node on UR (cw) lies above L(u, w): replace w Step 3: 26 Return the pair of nodes (u, w) A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithm II. 2 (Upper common support)… Animation: 27 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithm II. 3 (Merger of two Voronoi diagrams)… Input: 2 Voronoi diagrams for PL and PR Output: Voronoi diagram for PL PR (=P) Step 1: Compute the convex hulls of PL and of PR Step 2: Find the upper common support by algorithm II. 2 and mark the corresponding regions as active Step 3: Find a point at infinity upward on the bisector of the resulting pair of nodes Step 4: While the upper common support is not reached: Find the intersections of the bisector with the active Voronoi region Choose the intersection with the largest y coordinate and activate the next Voronoi region 28 Step 5: Add the polygonal line and delete superfluous portions from the input Voronoi diagrams A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithm II. 3 (Merger of two Voronoi diagrams)… Animation: 29 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Data Structure Doubly-connected edge-list Advantages: Direct access to neighboring regions Iterate the boundary in linear time But: Half-lines & full lines cannot be represented in this structure Bounding box 30 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithm II. 1 (Divide and conquer method)… Input: List of sites P in increasing order of the x coordinates Output: Voronoi diagram for P Step 1: if n<=3 Direct construction of the Voronoi diagram Step 2: otherwise: Divide P into PL and PR (using the median) Compute Voronoi diagrams for PL and PR by algorithm II. 1 (recursion!) Merge the Voronoi diagrams by algorithm II. 3 Step 3: 31 Return the Voronoi diagram for P A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithmic Complexity – Algorithm II. 1 (Divide and conquer method)… Input: List of sites P in increasing order of the x coordinates Output: Voronoi diagram for P Step 1: if n<=3 Direct construction of the Voronoi diagram Step 2: otherwise: Divide P into PL and PR (using the median) O(1) Compute Voronoi diagrams for PL and PR by algorithm II. 1 (recursion!) 2 T(n/2) Merge the Voronoi diagrams by algorithm II. 3 Step 3: 32 Return the Voronoi diagram for P A. Stadler – Geoinformation Technology: lecture 11 O(n) 2 T(n/2)+O(n) WS 2006/07
Computation II: Divide and Conquer Department of Geoinformation Science Algorithmic Complexity – Alg II. 1 (Divide and conquer method)… T(n) = 2 T(n/2)+O(n) Recursion tree method Total: T(n) = O(nlogn) optimal in terms of worst case complexity 33 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 34 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science Sweep Plane Method… Another fundamental paradigm for designing efficient geometric algorithms Strategy: Sweep a line – the sweep line – from top to bottom across the plane While sweeping, information is maintained regarding the structure one wants to determine Information only changes at certain special points, the so called event points 35 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science …sites above l …sweep line l Problem: …sites below l 36 The part of V above the sweep line l depends also on sites below l (l enters a Voronoi region before it meets its corresponding site!) A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science DON’T maintain intersections of V with the sweep line DO maintain information about the part of V that is fix corresponding to the sites above l What part of V cannot be changed anymore? Which points q above l have fixed nearest neighbors? Definition of an area consisting of points q that lie nearer to any known site (above l) than to l and the points below All sites below l cannot influence q Special geometric feature, which is defined as the set of points with the same distance of a point and a line PARABOLA 37 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science Every site above l defines a complete parabola Merging all areas limited by the parabolas Boundary = beach line The beach line is the function that (for each x coordinate) passes through the lowest point of all parabolas At the breakpoints, the distance to both sides must be the same as the distance to l While sweeping, the breakpoints of the beach line exactly trace out the Voronoi diagram We maintain the beach line (but not explicitly since it changes continuously) – more later… 38 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science Changes in the combinatorial structure of the beach line – Where and how? When a new parabolic arc appears on it (a). When a parabolic arc shrinks to a point and disappears (b). (a) 39 A. Stadler – Geoinformation Technology: lecture 11 (b) WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science When does a new parabolic arc appear? When the sweep line reaches a new site. a new parabola is formed (defined by the new site) The part of the new parabola below the old beach line is now part of the new beach line. degenerated parabola with zero width 40 A. Stadler – Geoinformation Technology: lecture 11 SITE EVENT WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science What happens to V at a site event? 2 breakpoints appear trace out the same edge (they coincide first and then move in opposite directions) = the only way for the appearing of new arcs (max 2 n-1 arcs). initially the edge is not connected to the rest of V 41 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science When does an existing arc shrink to a point and disappear? When α' disappears, all 3 consecutive arcs meet in a point q, which is equidistant from l and each of the 3 sites circle through the 3 sites, q as center, touching l No 3 α other sites can be inside the circle equidistant sites to q, no nearer sites possible node of V α' α'' q q CIRCLE EVENT 42 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science 2 possibilities for changes in the structure of the beach line SITE EVENT The sweep line reaches a new site A new arc appears A new edge starts to grow CIRCLE EVENT Circle through 3 sites, q as center, touching l An existing arc drops out 2 growing edges meet to form a Voronoi node in q 43 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science Representation of the beach line: balanced binary tree T Leaves = arcs of the beach line (ordered!) stored value: site defining the arc pointer to: node representing circle event when arc disappears Internal nodes = breakpoint between arcs stored value: ordered tuple of sites (defining left/right parabola) pointer to: edge that breakpoint traces out We can find the arc lying above a new site in O(logn) time 44 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science The event queue Q as a priority queue for handling events Ordered Stored by y coordinates priority events: SITE EVENT site is stored (known in advance) CIRCLE EVENT lowest point of the circle is stored pointer to the leaf in T representing the disappearing arc (have to be detected) 45 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science Algorithm III (Sweep plane method)… Input: A set P of sites in the plane Output: V for P Step 1: Initialize the event queue Q with all site events Step 2: WHILE Q is not empty DO Remove the event with largest y coordinate from Q pi IF the event is a site event THEN Handle. Site. Event(pi) ELSE Handle. Circle. Event(γ), where γ is the leaf of T representing the arc that will disappear Step 3: the remaining internal nodes in T represent half-infinite edges add them to V Step 4: Return V 46 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science Procedure Handle. Site. Event(pi) Step 1: IF T is empty, insert pi into it Step 2: ELSE Search in T for arc α vertically above pi IF α has pointer to circle event false alarm delete from Q Replace in T the leaf containing α by a subtree with the leaves pj, pi, pj and the internal nodes <pj, pi> and <pi, pj> rebalancing Create new record in V for half-edge separating V(pi) and V(pj) Check the triple of consecutive arcs (for pi… left/right arc) IF the breakpoints converge insert circle event into Q 47 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science Procedure Handle. Circle. Event(γ) Step 1: Delete the leaf γ from T representing the disappearing arc α update internal nodes, rebalancing Delete circle events involving α from Q Step 2: Add the center of the circle causing the event as a node to V Step 3: Create 2 half-edge records in V corresponding to the new breakpoint of the beach line Step 4: Check the new triples of consecutive arcs (where the former left/right neighbor of α is the middle arc) IF the breakpoints converge insert circle event into Q 48 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Computation III: Sweep Plane (Fortune) Department of Geoinformation Science Algorithmic Complexity - Algorithm III (Sweep plane method)… Primitive operations on the tree and the event queue: O(logn) Primitive operations on the doubly-connected edge list: O(1) Handling an event: constant number of such primitive operations O(logn) n site events O(nlogn) Linear number of circle events O(nlogn) Total: O(nlogn) optimal in terms of worst case complexity 49 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Voronoi Diagrams – An Overview Department of Geoinformation Science Introduction Geometric Interpretation History Applications Definition Construction Computing the Voronoi Diagram Intersection of Half-planes Divide and Conquer Sweep Plane (Fortune) Literature 50 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Literature Department of Geoinformation Science Okabe, Atsuyuki et al. : Spatial Tesselations: Concepts and Applications of Voronoi-Diagramms. John Wiley & Sons, 2000. de Berg, van Kreveld et al. : Computational Geometry – Algorithms and Applications. Springer, 2000. http: //home. dti. net/crispy/Voronoi. html Voronoi game http: //www. pi 6. fernuni-hagen. de/Geom. Lab/Voro. Glide/index. html. en Interactive generation of Voronoi diagram, Delaunay traingulation and Convex hull http: //www. snibbe. com/scott/bf/index. htm Voronoi art 51 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Outlook Department of Geoinformation Science So, maybe in future you will explore some astonishing analogies. . . 52 A. Stadler – Geoinformation Technology: lecture 11 WS 2006/07
Скачать презентацию Technische Universität Berlin Department of Geoinformation Science Geoinformation
5c5e92bdea3ffadb19294e39763dbc72.ppt