Transcript Slide 1

An Introduction to
Computational Geometry
Joseph S. B. Mitchell
Stony Brook University
Mini-Course Info



Joe Mitchell (Math 1-109)
[email protected] (631-632-8366)
http://www.ams.sunysb.edu/~jsbm/nasa/tutorial.html
Recommended Texts:
• Computational Geometry, by Mark de Berg, Marc van Kreveld,
Mark Overmars and Otfried Schwartzkopf (2nd edition)
• Computational Geometry in C, by Joe O’Rourke (2nd Edition)

Surveys/Reference books:
• CRC Handbook of Discrete and Computational Geometry,
editted by Goodman and O’Rourke (2nd edition)
• Handbook of Computational Geometry, editted by Sack and
Urrutia
Sessions and Topics

Session 1:
• Introduction
• Convex Hulls in 2D, 3D
• A variety of algorithms

Session 2:
• Searching for intersections: Plane sweep
• Halfspace intersections, low-dimensional
linear programming (LP)
• Range Search
3
Sessions and Topics

Session 3:
• Triangulation
• Proximity: Voronoi/Delaunay diagrams
• Point location search

Session 4:
• Arrangements
• Duality and applications
• Visibility graphs, shortest paths
4
Sessions and Topics

Session 5:
• Spatial partitioning, optimization, load
balancing
• Binary Space Partitions
• Clustering trajectories

Session 6:
• GeoSect: Overview and algorithms
• 3D partitioning
• Flow-respecting partitioning
5
Outline



Triangulation
Proximity: Voronoi/Delaunay
Point location search
6
Triangulation

Input: Set S of n points

Input: Other shapes
Simple polygon

Polygon with holes
Planar Straight-Line Graph (PSLG)
3D: Surfaces and solids (tetrahedralization)
7
Triangulation applet for simple polygons
Triangulation Theory in 2D

Thm: A simple polygon has a triangulation.
• Lem: An n-gon with n4 has a diagonal.




Also with holes
But, NOT true in 3D!
Thm: Any triangulation of a simple n-gon has n-3
diagonals, n-2 triangles.
Thm: The “dual” graph is a tree.
Thm: An n-gon with n4 has 2 “ears”.
Thm: The triangulation graph can be 3-colored.
Proofs: Induction on n.
8
Triangulating a Simple Polygon


Simple “ear-clipping” methods: O(n2 )
Cases with simple O(n) algorithms:
fan
• Convex polygons (trivial!)
• Monotone polygons, monotone mountains

General case (even with holes!):
• Sweep algorithm to decompose into
monotone mountains
• O(n log n)

Best theoretical results:
• Simple polygons: O(n) [Chazelle’90]
Not practical!
• Polygons with h holes: O(n+h log1+ h), (n+h log h) [BC]

Good practical method: FIST [Held], based on clever
9
2
methods of ear clipping (worst-case O(n ) )
FIST: Fast Industrial-Strength Triangulation
Based on
ear clipping
Simple polygon
FIST
Constrained Delaunay
Works nicely also for highly degenerate and “crazy” polygons
http://www.cosy.sbg.ac.at/~held/projects/triang/triang.html
3D cycles
Ear-Clipping Triangulation
Input: Simple polygon P
vi+1
vi-1
vi
pq is a diagonal, cutting off a
single triangle (the “ear”)
Naive: O(n3)
Smarter: Keep track of “ear tip
status” of each vi (initialize: O(n2) )
Each ear clip requires O(1) ear tip
tests ( @ O(n) per test )
Thus, O(n2) total, worst-case
Ear-clipping applet
Triangulating a Monotone
Mountain in O(n)

Ear clipping is easy!
• Testing if vi-1 vi+1 is a diagonal takes only
O(1) time
vi-1
vi is ear tip iff Left(vi+1 , vi , vi-1 )
Just traverse vertices from top
to bottom. Test/clip ears. If an
ear is clipped, re-test the earity
of the upper endpoint (vi-1 ) of the
diagonal just clipped.
vi
vi+1
monotone
12
Triangulation in O(n log n)

(1) Plane sweep to get horizontal trapezoidalization
L
Fire bullets
left/right
from each
vertex
SLS: left-to-right ordering of segments crossed by L (balanced binary tree)13
Events: L hits a vertex
Time: O(n log n)
Triangulation in O(n log n)

(2) Join top vertex to bottom vertex in each trapezoid
Lemma: Resulting pieces are monotone mountains
14
Triangulation in O(n log n)

(3) Triangulate each monotone mountain
Triangulate each, in time O(ni ), for total time O(n)
15
Summary: O(n log n) to triangulate n points or a planar straight-line graph (PSLG)
Triangulation


Sweep gives O(n log n) for triangulating PSLGs
Best theoretical results:
• Simple polygons: O(n) [Chazelle’90]
• Polygons with h holes: O(n+h log1+ h), (n+h log h) [BC]


Good practical method: FIST [Held], based on clever
methods of ear clipping (worst-case O(n2 ) )
3D decompositions (into tetrahedra, etc) are much
harder!
• NP-complete to decide if it is even possible (without adding
Steiner points); may need (r2 ) Steiner points
r= # reflex edges


Huge field of quality mesh generation See J. Shewchuk webpage
Important special class of triangulations: Delaunay 16
(next!)
Outline



Triangulation
Proximity: Voronoi/Delaunay
Point location search
17
Voronoi Diagrams
Georgy Voronoi
1868-1908
Historical Origins and Diagrams in Nature
René Descartes
1596-1650
1644: Gravitational
Influence of stars
Dragonfly wing
Giraffe pigmentation
Honeycomb
Ack: Streinu&Brock
Constrained soap bubbles
19
Starbucks
Post Office Problem
Query point
Post offices
Voronoi Diagram
Partition the plane according to the equivalence
classes: V(Q) = { (x,y) : the closest sites of S
to (x,y) is exactly the set Q  S }
• |Q| = 1
• |Q| = 2
• |Q|  3
Voronoi cells (2-faces)
Voronoi edges (1-faces)
Voronoi vertices (0-faces)
Q cocircular
21
Example
p
p
Voronoi cell of p
22
Delaunay Diagram
Join pi to pj iff there is an “empty circle”
A witness to the
through pi and pj
Delaunayhood of
If no 4 points cocircular
(degenerate), then
Delaunay diagram is a
(very special)
triangulation.


pi
(pi , pj)
pj
Equivalent definition: “Dual” of the
Voronoi diagram
Applet [Chew]
23
Voronoi and Delaunay Properties





The planar dual of Voronoi, drawn with nodes at the
sites, edges straight segments, has no crossing
edges [Delaunay]. It is the Delaunay diagram, D(S)
(defined by empty circle property).
Combinatorial size:  2n-5 Voronoi/Delaunay edges;
 3n-6 Voronoi vertices (Delaunay faces) O(n)
A Voronoi cell is unbounded iff its site is on the
boundary of CH(S)
CH(S) = boundary of unbounded face of D(S)
Delaunay triangulation lexico-maximizes the angle
vector among all triangulations
•
In particular, maximizes the min angle
24
Voronoi and Delaunay Properties



In any partition of S, S=B  R, into blue/red points,
any blue-red pair that is shortest B-R bridge is a
Delaunay edge.
•
•
•
D(S) is connected
MST  D(S) (MST=“min spanning tree”)
NNG  D(S) (NNG=“nearest neighbor graph”)
•
•
•
Divide and conquer
Sweep
Randomized incremental
Voronoi/Delaunay can be built in time O(n log n)
VoroGlide applet
25
Delaunay: Edge Flip Algorithm
Lawson Edge Swap
“Legalize” [BKOS]



Assume: No 4 co-circular points, for simplicity.
d
Start with any triangulation
Keep a list (stack) of “illegal” edges:
•
•




ab is illegal if InCircle(a,c,b,d) Locally non-Delauany
iff the smallest of the 6 angles goes up if flip
b
a
c
Flip any illegal edge; update legality status of
neighboring edges
Continue until no illegal edges
Theorem: A triangulation is Delaunay iff there are
no illegal edges (i.e., it is “locally Delaunay”)
Only O(n2 ) flips needed.
26
Connection to Convex Hulls in 3D



Delaunay diagram  lower convex hull of
the lifted sites, (xi , yi , xi 2 +yi 2 ), on the
parabaloid of revolution, z=x2 + y2
Upper hull  “furthest site Delaunay”
3D CH applet
27
Voronoi and Delaunay

Algorithms:
•
•
•
•
Divide and conquer (first O(n log n))
Sweep
Randomized incremental
Any algorithm that computes CH in R3 ,
Qhull website
e.g., QHull
28
Fortune’s Sweep Algorithm
29
Applet
Ack: Guibas&Stolfi
Parabolic Front
30
Applet
Ack: Guibas&Stolfi
Site Events
a)
b)
c)
31
Applet
Ack: Guibas&Stolfi
Circle Events
32
Applet
Ack: Guibas&Stolfi
Scheduling Circle Events
33
Applet
Ack: Guibas&Stolfi
Incremental Construction

Add sites one by one, modifying the
Delaunay (Voronoi) as we go:
• Join vi to 3 corners of triangle containing it
• Do edge flips to restore local Delaunayhood


If added in random order, simple
“Backwards analysis” shows expected time
O(n log n) [Guibas, Knuth, Sharir]
Applet (Shows Voronoi and Delaunay and empty circles)
34
Example
CS691G Computational Geometry - UMass
Amherst
35
Voronoi Extensions



Numerous ! [see Okabe, Boots, Sugihara, Chiu]
Different metrics
Higher dimensions:
• Delaunay in Rd  lower CH in Rd+1
• O( n log n + n (d+1)/2  )




Order-k
Other sites: Voronoi of polygons, “medial axis”
Additive/multiplicative weights; power diagram
-shapes: a powerful shape representation
36
GeoMagic, biogeometry at Duke
VRONI: Fast, robust Voronoi of polygonal domains
Incremental algorithm
37
38
Computing offsets with a Voronoi diagram
39
Alpha Shapes, Hulls
“Erase” with a
ball of radius 
to get -shape.
Straighten edges
to get -hull
40
41
-Shapes


Theorem: For each Delaunay edge,
e=(pi,pj), there exists min(e)>0 and
max(e)>0 such that e  -shape of S iff
min(e)    max(e).
Thus, every alpha-hull edge is in the
Delaunay, and every Delaunay edge is in
some alpha-shape.
Applet
Application of Triangulation:
Guarding a Polygon


Sensor placement, “guarding”
Determine a small set of “guards” that see all of a
given polygon P
Art Gallery Theorem
[Klee; Chvatal; Fisk]

For a simple n-gon, guards are always
sufficient and sometimes necessary.
Fisk proof that n/3 guards
are enough (even can use vertex guards)

• Triangulate
• 3-color the triangulation
• Put guards at vertices of the smallest color class
Chvatal comb: n/3
guards are necessary
45
Vertex Guarding a Simple Polygon

Vertex guarding applet
Outline



Triangulation
Proximity: Voronoi/Delaunay
Point location search
47
Point Location

Point-in-polygon test: Is point q inside
simple polygon P ?
Naïve: O(n) per test
CG: O(log n)
q
P n-gon
5 crossings  q inside
Point Location

More generally: Search for which face
of a planar map contains query point q
q

Example: Map of Manhattan Applet
Trapezoid under the mouse is highlighted
49
One Method: Partition the Plane into Slabs
Query time : O(log n)
binary search in x and then
binary search in y direction.
Storage space O(n2)
n/2
n/2
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
50
An Optimal Method: Kirkpatrick
Hierarchical Triangulation


Begin with full triangulation
Remove a “large” set of low-degree ( 9)
vertices that are independent (no two adjacent)
(such a set always exists in a planar graph, by Euler)



Retriangulate the resulting “holes” O(n)
Repeat until only the big outer triangle
remains “large” :  Cn  only O(log n) iterations
Build DAG for searching
Kirkpatrick Example
Full triangulation, T0
Triangles in T0
T1
T2
53
T3
T4
Final DAG hierarchy and search path
for query point p
T5
54
Better: Partitioning into Trapezoids
Assumption :
Segments are in “general
position“
R
Observation :
Every vertical edge has
one point in common with
a segment end.
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
55
Observations
f
R
f is convex
f is bounded
Every non - vertical side
of f is part of a segment
of S or an edge of R
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
56
Trapezoidal Decomposition
Lemma : Each face in a trapezoidal map of a set S of line segments
in general position has 1 or 2 vertical sides and exactly
two non-vertical sides
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
57
Left Edge of a Trapezoid
R
For every trapezoid
Δ T(S), except the
left most one, the left
vertical
edge of Δ is defined by
a segment endpoint p,
denoted by leftp(Δ) .
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
58
The 5 Cases: Left Edges of a Trapezoid
a)
b)
top()
leftp()
leftp()
top()
bottom()
bottom()
d)
c)
top()
top()
leftp()
bottom()
leftp()
bottom()
e) It is left edge of R. This case occurs for a single trapezoid of
T(S) only, namely the unique leftmost trapezoid of T(S)
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
59
Size of the Trapezoidal Map
Theorem: The trapezoidal map T(S) of a set of n line segments in
general position contains at most 6n + 4 vertices and at most 3n + 1
trapezoids.
Proof (1): A vertex of T(S) is either
- a vertex of R or
4
- an endpoint of a segment in S or
2n
- a point where the vertical extension starting
in an endpoint abuts on another segment
2 * (2n)
or on the boundary R.
6n + 4
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
60
Size of the Trapezoidal Map
Theorem: The trapezoidal map T(S) of a set of n line segments in
general position contains at most 6n + 4 vertices and at most 3n + 1
trapezoids.
Proof (2):Each trapezoid has a unique point leftp(), which is
- the lower left corner of R
1
- the left endpoint of a segment (can be
leftp() of at most two different trapezoids)
2n
- the right endpoint of a segment (can be leftp() of
n
atmost one trapezoid)
3n + 1
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
61
Adjacent Trapezoids
Two trapezoids  and ´ are adjacent if they meet along
a vertical edge.
1) Segments in general position :
A trapezoid has atmost four
adjacent trapezoids
5
1
2

4
3
1
2
3
4
5
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann

6
2) Segments not in general
position: A trapezoid can have
an arbitrary number of
adjacent trapezoids.
62
Vertical Neighbours: Upper, Lower Left
Trapezoid ' is (vertical) neighbor of 
top() = top(’) or bottom() = bottom(’)
In the first case ´ is upper left
neighbor of , in the second case ´
is lower left neighbor of .
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
63
Representing Trapezoidal Maps
There are records for all line segments and endpoints of S,
the structure contains records for trapezoids of T(S), but
not for vertices or edges of T(S).
The record for trapezoid  stores pointers to top(), and
bottom(), pointers to leftp() and rightp() and finally
pointers to its atmost 4 neighbors.
 is uniquely defined by top(), bottom(), leftp() and
rightp().
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
64
Search Structure: DAG
1l
B
A
1l
1
2l
C
1r
D
E
2
A
1r
G
1
2r
F
B
2r
2l
C
End points decide between left, right
Segments decide between below, above
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
65
2
2
D
E
F
G
A Randomized Incremental Algorithm
Input : A set S of n non-crossing line segments
Output : The trapezoidal map T(S) and a search structure
D(S) for T(S) in a bounding box.
Determine a bounding box R, initialize T and D
Compute a random permutation s1,s2, ..., sn of the elements of S
for i = 1 to n
do add si and change T(Si - 1) into T(Si) and D(Si -1) into D(Si)
Invariant :
In the step i T(Si) is correct trapezoidal map of Si
and D(Si) is an associated search structure.
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
66
A Randomized Incremental Algorithm
Input : A set of n non-crossing line segments
Output : The trapezoidal map T(S) and a search structure D
for T(S) in a bounding box.
Determine a bounding box R, initialize T and D
Compute a random permutation s1,s2, ..., sn of the elements of S
for i = 1 to n
do Find the set 0, 1, ... , k of trapezoids in T properly
intersected by si.
Remove 0, 1,..., k from T and replace them by new
trapezoids that appear because of the intersection of si.
Remove the leaves for 0, 1,..., k from D and create
leaves for the new Trapezoids.
Link the new leaves to the existing inner nodes by adding
some new inner nodes.
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
67
Questions
How can we find the intersecting trapezoids?
How can T and D be updated
a) if new segment intersects no previous trapezoid ?
b) if new segment intersects previous trapezoids ?
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
68
Finding the Intersecting Trapezoids

0
j+1 is lower right neighbour of
j
1
si

2

3
In T(Si) exactly those trapezoids are changed,
which are intersected by si
if rightp(j) lies above si
then Let j+1 be the lower right neighbor of j.
else Let j+1 be the upper right neighbor of j
Clue :
0 can be found by a query in the search structure D(Si -1)
constructed in iteration stage i -1.
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
69
New Segment Completely Contained in Trapezoid
t
Decomposition
s
p
q
T(Si-1)
A C
D
T(Si)
Search structure
D(Si-1)
D(Si-1)
t
p
q
A
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
4 new trapezoids, -1 old
trapezoid,
search depth + 3
70
B
B
s
C
D
New Segment Intersects Previous Ones
Decomposition
p
s
t
2
t t
0 1
t
3
E
D
A
t
1
t
2
T(Si)
Search structure
D(Si)
s
t
3
6 new trapezoids, – 4 old,
depth + 2
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
F
C
B
T(Si - 1)
D(Si - 1)
t
0
q
71
s
s
B
D
A
C
q
s F
E
Estimation of the Depth of the Search Structure
Let S be a set of n segments in general position,
q be an arbitrary fixed query point.
Depth of D(S):
worst case : 3n,
average case : O(log n)
Consider the path traversed by the query for q in D
Let Xi = # of nodes on the search path for q created in
iteration step i.
Xi <= 3
Pi = probability that there exists node on the search
path for q that is created in iteration step i.
E[Xi] <= 3 Pi
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
72
Observation
Iteration step i contributes a node to the search path for q
exactly if q(Si - 1), the trapezoid containing q in T(Si - 1), is not the
same as q(Si) , the trapezoid containing q in T(Si)
Pi = Pr[q(Si)  q(Si - 1)].
If q(Si) is not same as q(Si - 1), then q(Si) must be one of the
trapezoids created in iteration i.
q(Si) does not depend on the order in which the segments in Si
have been inserted.
Backwards analysis :
We consider T(Si) and look at the probability that q(Si)
disappears from the trapezoidal map when we remove the
segment si.
q(Si) disappears if and only if one of top(q(Si)), bottom(q(Si)),
leftp(q(Si)), or right(q(Si)) disappears with removal of si .
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
73
Observation
Prob[top(q(Si))] = Prob[bottom(q(Si))] = 1/i.
Prob[leftp(q(Si))] disappears is at most 1/i.
Prob[rightp(q(Si))] disappears is at most 1/i.
Pi = Pr[q(Si)  q(Si - 1)] = Pr[q(Si)  T(Si - 1)]
 4/i
n
n
12
1
n
 n
E  Xi   3Pi    12  12Hn  log n
i 1 i
i 1 i
 i 1  i 1
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
74
Analysis: Size of Search Structure
Leaves in D are in one – to – one correspondence with the
trapezoids in , of which there are O(n).
The total number of nodes is bounded by :
n
n   (# of inner nodes created in iteration step i)
i 1
The worst case upper bound on the size of the structure
n
n    i   n² 
i 1
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
75
Analysis: Size of Search Structure
Theorem: The expected number of nodes of D is O(n).
Proof: The # of leaves is in O(n). Consider the internal
nodes: Xi = # of internal nodes created in iteration step i
n
 n
E  Xi   EXi
 i 1  i 1
1
 , s  : 
0
if  disappears from T(Si) when s is removed from Si
otherwise
There are at most four segments that cause a given trapezoid
to disappear
  , s   4 T S   i 
i
sSi T ( Si )
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
76
Analysis: Size of Search Structure
1
i 
Eki      , s  
 1
i sSi T ( Si )
i
The expected number of newly created trapezoids is O(1) in
every iteration of the algorithm, from which the O(n) bound
on the expected amount of storage follows.
n

 E  Xi  n 
 i 1 
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
77
Summary
Let S be a planar subdivision with n edges. In O(n log n)
expected time one can construct a data structure that
uses O(n) expected storage, such that for any query
point q, the expected time for a point location query is
O(log n).
Computational Geometry, WS
2007/08
Prof. Dr. Thomas Ottmann
78