Computational Geometry Part 1: Intro, Primitives, Area, Circles - - PowerPoint PPT Presentation

Computational Geometry

Part 1: Intro, Primitives, Area, Circles

Lucca Siaudzionis and Jack Spalding-Jamieson 2020/02/25

University of British Columbia

Announcements

• A3 is due Sunday!
• Project topics are due next Tuesday!

1

What is Computational Geometry?

Computational Geometry is a huge field, with its own conferences and labs dedicated to it. It includes all sorts of problems. It can include problems involving:

• Point sets
• Line segments
• Polygons
• Convexity

Some classic problem classes you may have heard of are:

• Art gallery problems
• Voronoi diagrams and Delauney triangulations
• Ray tracing
• Graph drawing (which is a field in its own right)

2

Primitives: Vectors

What is a vector? Not dynamic arrays! A d-dimensional vector is a list of exactly d numbers. Normally we work with 2-dimensional and 3-dimensional vectors. Vectors are used for everything in geometry! They can be used to represent:

• Points
• Displacement
• Velocity
• Direction

3

Primitives: Vectors - Notation and Implementation

A 2D vector v can be denoted as a tuple, i.e. v = (x, y), and likewise a 3D vector can appear as (x, y, z). In code, this can be represented many ways, such as:

1

typedef long double ld;

2

struct vec2 {

3

ld x, y;

4

vec2() : x(0), y(0) {}

5

};

6

struct vec3 {

7

ld x, y, z;

8

vec3() : x(0), y(0), z(0) {}

9

};

10

template <size_t N>

11

struct vec {

12

array<ld,N> a; };

4

Primitives: Vectors - Operations (1)

There are many operations we can do on a vector, such as:

• Addition/subtraction:

v1 ± v2 = (x1 ± x2, y1 ± y2)

• Scalar multiplication: a

v = (ax, ay) These are the simplest operations, and in general can be implemented as follows:

1

template <size_t N> struct vec {

2

array<ld,N> a;

3

vec operator+(const vec& o) const {

4

vec ret;

5

for (size_t i = 0; i < N; ++i) ret.a[i] = a[i] + o.a[i];

6

return ret; }

7

vec operator*(const int c) const {

8

vec ret;

9

for (size_t i = 0; i < N; ++i) ret.a[i] = c*a[i];

10

return ret; }

5

Primitives: Vectors - Operations (2)

We can also measure the ”length” of a vector in many different ways:

• L1 norm:

v1 = |x| + |y| (“Manhattan distance”)

• L2 norm:

v2 =

• x2 + y 2 (Euclidean distance, often written as |

v|)

• L∞ norm:

v∞ = max(|x|, |y|)

• Lk norm (general):

vk = (|x|k + |y|k)1/k Generally, this can be implemented as:

11

ld norm(ld k) const {

12

ld ret = 0;

13

for (size_t i = 0; i < N; ++i) ret += powl(abs(a[i]),k);

14

return powl(ret, 1/k); }

6

Primitives: Vectors - Operations (3)

Dot product

v1 · v2 = x1x2 + y1y2 = | v1|| v2| cos θ, where θ is the angle between v1 and v2

v · v = | v|2

• Interpretation: signed length of projection from

v1 to v2, scaled by length | v2|

v1 · v2 = 0 ⇔ v1 ⊥ v2 Generally, this can be implemented as:

15

friend vec dot(const vec& u, const vec& v) {

16

vec ret;

17

for (size_t i = 0; i < N; ++i) ret[i] = u.a[i] + v.a[i];

18

return ret; }

19

};

7

Primitives: Vectors - Operations (4)

2D cross-product:

v1 × v2 = x1y2 − x2y1 = | v1|| v2| sin θ

• Interpretation: signed distance from

v1 to its projection on v2, scaled by length | v2|

v1 × v2 = 0 ⇔ v1 v2

• Returns a scalar.

This can be implemented as:

1

struct vec2 {

2

ld x, y;

3

vec2() : x(0), y(0) {}

4

friend ld cross(const vec2& u, const vec2& v) {

5

return u.x*v.y-v.x*u.y; }

6

};

8

Primitives: Vectors - Operations (5)

3D cross-product:

• v1 ×

v2 =    x1 y1 z1    ×    x2 y2 z2    =

• i

j k x1 y1 z1 x2 y2 z2

• =

   y1z2 − z1y2 z1x2 − x1z2 x1y2 − y1x2   

• Still true: |

v1 × v2| = | v1|| v2| sin θ

• Interpretation: vector that is perpendicular to both

v1 and v2

• Direction determined by left-hand rule
• NOT associative, NOT commutative
• Returns a vector.

9

Primitives: Vectors - Operations (6)

The 3D cross-product can be implemented succinctly as:

7

struct vec3 {

8

array<ld,3> a;

9

friend vec3 cross(const vec3& u, const vec3& v) {

10

vec3 ret;

11

for (size_t i = 0; i < 3; ++i)

12

ret[i] = u.a[(i+1)%3]*v.a[(i+2)%3]-u.a[(i+2)%3]*v.a[(i+1)%3];

13

return ret; }

14

};

10

Primitives: Turning Left or Right

How can we tell if the 2D points x0, x1, x2 in that order are oriented clockwise or counterclockwise?

• x0
• x1
• x2
• x0
• x1
• x2

Figure 1: Left: 3 points that are CCW; Right: 3 points that are CW

i.e. how do we differentiate left and right turns?

11

Primitives: Turning Left or Right (Solution)

( x1 − x0) × ( x2 − x1)        > 0 Counterclockwise = 0 Collinear < 0 Clockwise If you remember your first year physics courses, you can use the right-hand rule to derive this. Otherwise, take this as a primitive operation.

12

Complex numbers

Implementing all these linear algebra routines seems tedious. Luckily, in 2D we have a representation using complex numbers! Recall that complex numbers are defined as C = {a + bi | a, b ∈ R, i = √ −1} Then, we can represent a 2D point (x, y) by x + yi. We can add/subtract/multiply/divide complex numbers by multiplying the real and imaginary parts, then remembering that i2 = −1. Addition and subtraction is the same as with vectors (translation). What about multiplication/division?

13

Complex numbers

Recall that we can represent points in 2D using polar coordinates (r, θ). Euler’s formula: eiθ = cos θ + i sin θ Then (r, θ) = r cos θ + ir sin θ = reiθ, so we get (r1eiθ1)(r2eiθ2) = (r1r2)ei(θ1+θ2) Multiplication = multiply absolute values, add angles! In particular, if r2 = 1 then this is a rotation of a point about the origin!

14

Complex numbers

Complex numbers make geometry much easier!

• Vector addition: add complex numbers
• Multiply vector by scalar: multiply complex by real
• Dot/cross product: zw = (

z · w) + ( z × w)i

• if z = a + bi then z = a − bi. This is called the conjugate of z.
• Rotation by θ: multiply by eiθ
• Length: |z|

15

Primitives: Representations of 2D Vectors

There are many different ways we can implement our 2D vectors: ld x, y struct vec2 vec<2> complex<T> Speed Fast Fast Fast Slow Bug frequency High Low Very low Negligible Code length Very short Short Short Very short! Vector Arithmetic Hard-code Operator Operator Operator Rotation Hard-code Matrix Matrix (difficult) Multiply by eiθ Dot Product Hard-code Function Function Re(zw) Cross Product Hard-code Function Function Im(zw) Length/Distance Hard-code Function Function |z| In C++, complex numbers are implemented by complex<T>. In Python, complex numbers are builtin (e.g. 1 + 2j). In Java, Point2D.Double is fairly useless for lots of typing; just write your own vector/complex class.

16

Primitives: Representing Lines / Line Segments (1)

Representing lines

• Bad: y = mx + b (what can possibly go wrong with this?)
• Okay: ax + by + c = 0
• Good:

L(t) = x0 + v0t Representing line segments

S(t) = x0 + ( x1 − x0)t where t ∈ [0, 1] Both lines and line segments can (and should) be represented by a pair of vectors!

17

Primitives: Representing Lines / Line Segments (2)

• x0
• x1
• v0

Figure 2: Parametric representation of a line

18

Primitives: Line Segment Intersection (Seg X Seg) (1)

Given two line segments L1 and L2, how do we tell if they intersect?

Intersects Does not Intersect Endpoint intersection Overlap (Colinear)

Figure 3: Types of intersections

19

Primitives: Line Segment Intersection (Seg X Seg) (2)

After choosing one of the segments as a ”primary” segment (with vector A), create some vectors:

• B
• C
• A

Notice how we have two kinds of turns: A → B is a left turn, and A → C is a right turn.

20

Primitives: Line Segment Intersection (Seg X Seg) (3)

What can we say about intersections, given the turns?

• If both turns are left turns, then the segments do not intersect.
• If both turns are right turns, then the segments do not intersect.
• We can try this with either segment as the “primary” segment.
• If we never satisfy the first two points with either primary, then the two segments do

intersect (except for some edge cases).

21

Primitives: Line Segment Intersection (Seg X Seg) (4)

What are the edge cases?

• Colinearity!
• Endpoint intersection (easy to handle).

The two segments are colinear if all the cross-products computed for the direction are zero. There are several different ways for the segments to be colinear.

22

Primitives: Line Segment Intersection (Seg X Seg) (5)

No overlap Interior overlap Shared endpoint Containment

b′ a′ b a b′ b a′ a b′ a′ b a a′ a b′ b

23

Primitives: Line Segment Intersection (Seg X Seg) (6)

How do we detect each of these cases? Since everything is colinear, it reduces to an interval intersection check.

24

Primitives: Line Segment Intersection (Seg X Seg) (7)

UBC Code Archive implementation:

1

bool seg_x_seg(pt a1, pt a2, pt b1, pt b2) {

2

//if (eq(a1,a2) || eq(b1,b2)) return false; // uncomment to exclude endpoints

3

ld za = abs(a2-a1), zb = abs(b2-b1); za=za>EPS?1/za:0; zb=zb>EPS?1/zb:0;

4

int s1 = sgn(cp(a2-a1, b1-a1)*za), s2 = sgn(cp(a2-a1, b2-a1)*za);

5

int s3 = sgn(cp(b2-b1, a1-b1)*zb), s4 = sgn(cp(b2-b1, a2-b1)*zb);

6

if(!s1 && !s2 && !s3) { // collinear

7

if (cmp_lex(a2, a1)) swap(a1, a2); if (cmp_lex(b2, b1)) swap(b1, b2);

8

//return cmp_lex(a1, b2) && cmp_lex(b1, a2);//uncomment to exclude endpoints

9

return !cmp_lex(b2, a1) && !cmp_lex(a2, b1);

10

} return s1*s2 <= 0 && s3*s4 <= 0; } //change to < to exclude endpoints

• There is an option to include endpoints (depends on problem).
• Heavy reliance on subroutines (e.g. cmp lex, sgn).
• Epsilon safety.
• Shortcut to direction checking: Multiply the cross-products.

25

Aside: Epsilon Safety

Recall that operations with floating points (and doubles) should often use an epsilon. Remember Aldssignment 1 question A? Typically, an epsilon of about 10−8 should be added before rounding a sum. Other examples of use-cases:

• To check if a value z is positive, evaluate z>EPS.
• To check if a value z is zero, evaluate abs(z)<EPS.
• To check if two values z1 and z2 are equal, evaluate abs(z1-z2)<EPS.

26

Area: Triangles (1)

How do we compute the area of a triangle (given by 3 points)?

• x0
• x1
• x2
• v1
• v0

27

Area: Triangles (2)

A useful fact about 2D vectors: For any two vectors v0 and v1, the value | v0 × v1| is equal to the area of the following parallelogram:

• v1
• v0
• v0
• v1

Using this, we can easily get the area of the triangle. It’s half the cross product!

• v0 =

x1 − x0

• v1 =

x2 − x0 A = 1 2 | v0 × v1|

• v1
• v0
• v0
• v1
• x0
• x1
• x2

28

Area: Polygons (Shoelace Formula)

Turns out computing the area of polygon is just as easy! If the vertices are x1, . . . xn in CCW

• rder, then

A = 1 2

• (

x1 × x2) + ( x2 × x3) + · · · + ( xn−1 × xn) + ( xn × x1)

• This works for any polygon! (Area will be signed if the polygon is self-intersecting.)
• x1
• x2
• x3
• x4
• x5
• x6

Figure 4: A polygon with vertices in CCW order

29

Pick’s Theorem

If polygon has vertices on lattice points (i.e. integer coordinates.) then the area is A = i + b 2 − 1, where i = # of interior lattice points, b = # of perimeter lattice points.

Figure 5: Illustration of Pick’s theorem, from Wikipedia

With some (non-trival) work, we can now count the lattice points inside a polygon.

30

Discussion Problem: Building a Convex Polygon

Given 1 ≤ n ≤ 106 2D points in the form of integer coordinates that represent the vertices of a convex polygon, find a sorted order of the points that proceeds CW around the polygon.

31

Discussion Problem: Building a Convex Polygon – Insight

After fixing any reference point, we can sort all the other points by their angle around that

• point. A good reference could be the lowest point in terms of y value, breaking ties by x value

(prioritizing the left side).

1 2 3 4 5

32

Aside: Integers

Everything we’ve done so far has used floating point numbers, but most operations and problems in geometry never need to use floating points (or rationals of any kind) at all! The following operations can all be done with purely integer values:

• addition/subtraction
• L1 norm, L∞ norm
• comparing L2 norms/distances (compare the norms squared!)
• comparing Ld norms
• dot and cross-products
• line segment intersection
• sorting by angle! (how do we do this?)

Many more operations can be done with integers too. Think about this when you see new

• perations or representations of data!

33

Circles: Representation and Primitives

A circle can be represented by two things: a centre-point c and a radius r. As you likely know, a circle has area πr 2 and radius 2πr.

1

struct circle {

2

pt cen;

3

ld r;

4

circle() : cen(), r(0) {}

5

ld area() const {

6

return M_PI * r * r;

7

}

8

ld perimeter() const {

9

return M_PI * 2 * r;

10

}

34

Circles: Intersection

Given two circles, can we check if they intersect? This is actually fairly easy. See if you can read the following code:

11

bool intersects(const circle& oth) const {

12

return (cen-oth.cen).norm2() <= r + oth.r;

13

}

14

};

35

Circles: Tangents (1)

What is the shortest path from A to B, with this circle as an obstacle? A B

36

Circles: Tangents (2)

Shortest path is 2 straight lines tangent to circle + 1 arc. Remember that circle tangents are perpendicular to radius. A B After some derivation, this can be computed in constant time.

37

Discussion Problem: Obstacles

What is the shortest path from a point A to a point B with up to 1 ≤ n ≤ 400 non-overlapping circles as obstacles in the way?

A B

38

Discussion Problem: Obstacles – Insight

We can construct a ”visibility” graph of the ”important” points. In addition to the points A and B, the important points fall into the following categories:

A or B → tangent point circle-circle parallel points

We end up with O(n2) different points. However, notice that any of these ”edges” we’ve drawn could actually intersect another circle, and so we need to check this for each circle (how do we do circle-line intersection?). Finally, run Dijkstra to find the shortest path. Overall runtime complexity: O(n3 + n2 log n) = O(n3)

39