Lesson 48 of 50 16 minMaster

Computational geometry basics

The cross product as the one tool behind most geometry code: left or right turns, segment intersection, polygon area with the shoelace formula, point in polygon, and the convex hull in O(n log n).

Site plans, CAD drawings, maps, games and robots all need geometry: do two walls cross, how big is this plot, is this point inside the boundary, what is the shortest fence around these pegs? Almost all of it is built on one operation, the cross product, and on one rule: avoid floating-point division when you can, because rounding errors make “exactly on the line” impossible to test.

Left, right or straight: the cross product

For three points O, A, B, the cross product of the vectors OA and OB is (A.X − O.X)·(B.Y − O.Y) − (A.Y − O.Y)·(B.X − O.X). Its sign says which way you turn when walking O → A → B: positive means left (counter-clockwise), negative means right (clockwise), zero means the three points are on one line. Its size is twice the area of the triangle OAB. With integer coordinates it is exact.

Orientation.cs

1using Pt = (long X, long Y);2 3Pt a = (0, 0), b = (10, 0);                          // a wall from A to B4Console.WriteLine(Turn(a, b, (5, 4)));5Console.WriteLine(Turn(a, b, (5, -3)));6Console.WriteLine(Turn(a, b, (20, 0)));7 8static long Cross(Pt o, Pt a, Pt b) => (a.X - o.X) * (b.Y - o.Y) - (a.Y - o.Y) * (b.X - o.X);9 10static string Turn(Pt o, Pt a, Pt b) => Cross(o, a, b) switch11{12    > 0 => "left (counter-clockwise)",13    < 0 => "right (clockwise)",14    _ => "straight (collinear)"15};

Output

left (counter-clockwise)
right (clockwise)
straight (collinear)

The using Pt = (long X, long Y); line (C# 12) gives the tuple type a short name. long keeps the products from overflowing for coordinates up to about a billion.

Do two segments cross?

Segments P1P2 and Q1Q2 cross properly when P1 and P2 lie on opposite sides of the line through Q1Q2, and Q1 and Q2 lie on opposite sides of the line through P1P2: four cross products with opposite signs. If a cross product is zero, a point lies on the other line, and a simple bounding-box check tells whether it is actually on the segment (touching at an end, or overlapping on one line).

Intersect.cs

1using Pt = (long X, long Y);2 3Console.WriteLine(Intersects((0, 0), (10, 10), (0, 10), (10, 0)));   // an X4Console.WriteLine(Intersects((0, 0), (10, 0), (0, 5), (10, 5)));     // parallel5Console.WriteLine(Intersects((0, 0), (10, 0), (10, 0), (10, 8)));    // touching at an end6Console.WriteLine(Intersects((0, 0), (4, 0), (6, 0), (10, 0)));      // same line, with a gap7 8static long Cross(Pt o, Pt a, Pt b) => (a.X - o.X) * (b.Y - o.Y) - (a.Y - o.Y) * (b.X - o.X);9 10static bool OnSegment(Pt p, Pt q, Pt r) =>           // r is collinear with pq: is it between them?11    Math.Min(p.X, q.X) <= r.X && r.X <= Math.Max(p.X, q.X) &&12    Math.Min(p.Y, q.Y) <= r.Y && r.Y <= Math.Max(p.Y, q.Y);13 14static bool Intersects(Pt p1, Pt p2, Pt q1, Pt q2)15{16    long d1 = Cross(q1, q2, p1), d2 = Cross(q1, q2, p2);17    long d3 = Cross(p1, p2, q1), d4 = Cross(p1, p2, q2);18    if (Math.Sign(d1) * Math.Sign(d2) < 0 && Math.Sign(d3) * Math.Sign(d4) < 0) return true;   // proper crossing19    return (d1 == 0 && OnSegment(q1, q2, p1)) || (d2 == 0 && OnSegment(q1, q2, p2))20        || (d3 == 0 && OnSegment(p1, p2, q1)) || (d4 == 0 && OnSegment(p1, p2, q2));21}

Output

True
False
True
False

Area of a polygon: the shoelace formula

List the corners of a polygon in order. For each edge from corner a to corner b add a.X·b.Y − b.X·a.Y (the cross product of the two position vectors). Half the absolute total is the area, for any simple polygon, convex or not. The sign of the total even tells you the direction: positive for counter-clockwise corners.

Area.cs

1using Pt = (long X, long Y);2 3Pt[] plot = [(0, 0), (40, 0), (40, 20), (20, 20), (20, 50), (0, 50)];   // an L-shaped plot, metres4Console.WriteLine($"Area: {Area(plot)} square metres");5 6static double Area(Pt[] poly)7{8    long twice = 0;9    for (int i = 0; i < poly.Length; i++)10    {11        Pt a = poly[i], b = poly[(i + 1) % poly.Length];12        twice += a.X * b.Y - b.X * a.Y;              // doubled signed area of the triangle (origin, a, b)13    }14    return Math.Abs(twice) / 2.0;15}

Output

Area: 1400 square metres

Is a point inside a polygon?

Cast a ray from the point to the right and count how many polygon edges it crosses. An odd count means inside, even means outside (every time the ray crosses the boundary it switches between inside and outside). This ray casting test is O(n) for n edges and works for concave shapes too.

PointInPolygon.cs

1using Pt = (double X, double Y);2 3Pt[] plot = [(0, 0), (40, 0), (40, 20), (20, 20), (20, 50), (0, 50)];4foreach (Pt p in new Pt[] { (10, 40), (30, 40), (30, 10) })5    Console.WriteLine($"{p}: {(Inside(plot, p) ? "inside" : "outside")}");6 7static bool Inside(Pt[] poly, Pt p)8{9    bool inside = false;10    for (int i = 0, j = poly.Length - 1; i < poly.Length; j = i++)11    {12        Pt a = poly[i], b = poly[j];13        if ((a.Y > p.Y) != (b.Y > p.Y) &&                          // the edge spans the ray's height14            p.X < (b.X - a.X) * (p.Y - a.Y) / (b.Y - a.Y) + a.X)    // and crosses it to the right of p15            inside = !inside;16    }17    return inside;18}

Output

(10, 40): inside
(30, 40): outside
(30, 10): inside

The convex hull

The convex hull of a set of points is the smallest convex polygon containing them all: the shape of a rubber band stretched around pegs. It is the shortest fence, the outline of a point cloud, and the first step of many other algorithms (farthest pair, collision shapes, minimum bounding box).

Andrew’s monotone chain sorts the points by x, then walks left to right building the lower edge. Each new point is checked against the last two points kept: if the turn is not a left turn, the last point lies inside the hull and is removed (perhaps several times). Then the same walk from right to left builds the upper edge.

ConvexHull.cs

1using Pt = (int X, int Y);2 3Pt[] pegs = [(30, 60), (70, 230), (120, 140), (150, 40), (160, 260), (200, 120),4             (230, 210), (260, 70), (270, 170), (100, 90), (190, 190), (60, 150)];5var hull = ConvexHull(pegs);6Console.WriteLine(string.Join(" ", hull));7Console.WriteLine($"{hull.Count} of {pegs.Length} pegs");8 9static long Cross(Pt o, Pt a, Pt b) => (long)(a.X - o.X) * (b.Y - o.Y) - (long)(a.Y - o.Y) * (b.X - o.X);10 11static List<Pt> ConvexHull(Pt[] points)12{13    var p = points.OrderBy(q => q.X).ThenBy(q => q.Y).ToArray();14    if (p.Length < 3) return p.ToList();15    var h = new List<Pt>();16    foreach (var q in p)                                         // lower edge, west to east17    {18        while (h.Count >= 2 && Cross(h[^2], h[^1], q) <= 0) h.RemoveAt(h.Count - 1);19        h.Add(q);20    }21    int lower = h.Count;22    for (int i = p.Length - 2; i >= 0; i--)                      // upper edge, east to west23    {24        while (h.Count > lower && Cross(h[^2], h[^1], p[i]) <= 0) h.RemoveAt(h.Count - 1);25        h.Add(p[i]);26    }27    h.RemoveAt(h.Count - 1);                                     // the first point came round again28    return h;29}

Output

(30, 60) (150, 40) (260, 70) (270, 170) (230, 210) (160, 260) (70, 230)
7 of 12 pegs

Sorting costs O(n log n). In the walk every point is added once and removed at most once, so it is O(n). The <= 0 removes collinear points too; use < 0 to keep them on the fence.

Problem	Tool	Time
Left or right turn	Cross product sign	O(1)
Do two segments cross?	Four cross products	O(1)
Polygon area	Shoelace formula	O(n)
Point in polygon	Ray casting	O(n)
Convex hull	Monotone chain	O(n log n)
Any two of n segments cross?	Sweep line (Bentley-Ottmann)	O((n + k) log n)
Closest pair of points	Divide and conquer, or a grid	O(n log n)

Quick check

The cross product for O → A → B is −12. What does it tell you?

Quick check

A ray from a point crosses the polygon boundary 4 times. Where is the point?

Practice

Write IsConvex for a polygon given in order: it is convex when all turns along the boundary go the same way. The square should print True, the L-shaped plot False.

Show a hint

Compute the cross product at every corner (corner i, i + 1, i + 2, wrapping around). Convex means you never see both a positive and a negative value.

Show the solution

Solution.cs

1using Pt = (long X, long Y);2 3Pt[] square = [(0, 0), (10, 0), (10, 10), (0, 10)];4Pt[] plot = [(0, 0), (40, 0), (40, 20), (20, 20), (20, 50), (0, 50)];5Console.WriteLine(IsConvex(square));6Console.WriteLine(IsConvex(plot));7 8static bool IsConvex(Pt[] poly)9{10    bool left = false, right = false;11    for (int i = 0; i < poly.Length; i++)12    {13        Pt a = poly[i], b = poly[(i + 1) % poly.Length], c = poly[(i + 2) % poly.Length];14        long cr = (b.X - a.X) * (c.Y - a.Y) - (b.Y - a.Y) * (c.X - a.X);15        if (cr > 0) left = true;16        if (cr < 0) right = true;17    }18    return !(left && right);                         // every turn the same way19}

Output

True
False

Remember

The cross product’s sign gives the turn direction (left, right, straight); its size is twice a triangle’s area.
Segment intersection, polygon area (shoelace) and point-in-polygon (ray casting) are all built on it.
Andrew’s monotone chain finds the convex hull in O(n log n): sort, then pop points that do not turn left.
Prefer integer coordinates; with doubles, compare against a tolerance instead of zero.

Foundations

Linear data structures

Searching and sorting

Hashing and problem-solving techniques

Trees

Graphs

Algorithm design

Master level

Computational geometry basics

Left, right or straight: the cross product

Do two segments cross?

Area of a polygon: the shoelace formula

Is a point inside a polygon?

The convex hull

Fencing the site: convex hull

Remember