Computational Geometry - Tutorials

Computational geometry is a branch of computer science that studies algorithms which can be expressed in terms of geometry. It is essential for fields such as computer graphics, robotics, GIS, and game development. This tutorial provides a thorough, implementation‑focused guide for software engineers who want to master geometric algorithms and apply them in real‑world projects.

Fundamental Concepts
Core Data Structures
Classic Algorithms
Advanced Topics
Implementation Tips & Best Practices
FAQ
Quiz

Fundamental Concepts

Before diving into algorithms, understand the basic geometric primitives and their mathematical representation. The most common primitives are:

Point: a pair (x, y) in 2‑D or a triple (x, y, z) in 3‑D.
Vector: an ordered tuple representing direction and magnitude.
Segment: a line portion bounded by two endpoints.
Ray: a half‑line starting at an origin and extending infinitely.
Polygon: a closed chain of line segments.

📝 Note: All coordinates are typically stored as double (or float) for precision, but integer arithmetic is preferred when exactness is required (e.g., for integer lattice problems).

Coordinate Systems & Orientation Tests

The orientation (or turn) test determines the relative position of three points and is the backbone of many geometric predicates.

cpp python

inline long long orient(const Point &a, const Point &b, const Point &c) {
    // >0 : counter‑clockwise, <0 : clockwise, 0 : collinear
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

def orient(a, b, c):
    """Return >0 if a,b,c make a counter‑clockwise turn, <0 for clockwise, 0 if collinear."""
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)

⚠ Warning: When using floating‑point numbers, be aware of rounding errors. Apply an epsilon tolerance (e.g., 1e‑9) to treat near‑zero values as zero.

Core Data Structures

KD‑Tree

A KD‑Tree is a space‑partitioning binary tree for organizing points in k‑dimensional space. It enables fast nearest‑neighbor and range queries, which are crucial for many geometric problems.

cpp python

struct KDNode {
    Point pt;
    KDNode *left, *right;
    int axis; // 0 for x, 1 for y
    KDNode(const Point& p, int ax) : pt(p), left(nullptr), right(nullptr), axis(ax) {}
};

KDNode* buildKD(vector<Point>& pts, int depth = 0) {
    if (pts.empty()) return nullptr;
    int axis = depth % 2;
    sort(pts.begin(), pts.end(), [&](const Point& a, const Point& b){
        return axis ? a.y < b.y : a.x < b.x;
    });
    size_t mid = pts.size() / 2;
    KDNode* node = new KDNode(pts[mid], axis);
    vector<Point> left(pts.begin(), pts.begin() + mid);
    vector<Point> right(pts.begin() + mid + 1, pts.end());
    node->left  = buildKD(left, depth + 1);
    node->right = buildKD(right, depth + 1);
    return node;
}

class KDNode:
    def __init__(self, point, axis):
        self.point = point
        self.axis = axis
        self.left = None
        self.right = None

def build_kd(points, depth=0):
    if not points:
        return None
    axis = depth % 2
    points.sort(key=lambda p: p[axis])
    median = len(points) // 2
    node = KDNode(points[median], axis)
    node.left = build_kd(points[:median], depth + 1)
    node.right = build_kd(points[median + 1:], depth + 1)
    return node

💡 Tip: When the dataset is static (no insertions/deletions), building the KD‑Tree once and reusing it yields optimal query performance.

Segment Tree for 1‑D Intervals

Segment trees allow fast query of intervals intersecting a given point or range. They are frequently used for sweep‑line algorithms.

cpp

struct SegTree {
    vector<vector<Segment>> tree;
    int n;
    SegTree(int size) : n(size) { tree.assign(4*n, {}); }
    void add(int idx, int l, int r, const Segment& seg) {
        if (l == r) { tree[idx].push_back(seg); return; }
        int mid = (l+r)/2;
        if (seg.l <= mid) add(idx*2, l, mid, seg);
        if (seg.r >  mid) add(idx*2+1, mid+1, r, seg);
    }
    // query point x returns all segments covering x
    void query(int idx, int l, int r, int x, vector<Segment>& out) {
        out.insert(out.end(), tree[idx].begin(), tree[idx].end());
        if (l == r) return;
        int mid = (l+r)/2;
        if (x <= mid) query(idx*2, l, mid, x, out);
        else query(idx*2+1, mid+1, r, x, out);
    }
};

Classic Algorithms

Convex Hull

The convex hull of a set of points is the smallest convex polygon that contains all points. Two popular O(n log n) algorithms are Graham Scan and Andrew’s Monotone Chain.

Andrew’s Monotone Chain (C++)

cpp python

vector<Point> convexHull(vector<Point> pts) {
    sort(pts.begin(), pts.end(), [](const Point& a, const Point& b){
        return a.x < b.x || (a.x == b.x && a.y < b.y);
    });
    vector<Point> lower, upper;
    for (const auto& p : pts) {
        while (lower.size() >= 2 && orient(lower[lower.size()-2], lower.back(), p) <= 0)
            lower.pop_back();
        lower.push_back(p);
    }
    for (int i = (int)pts.size()-1; i >= 0; --i) {
        const auto& p = pts[i];
        while (upper.size() >= 2 && orient(upper[upper.size()-2], upper.back(), p) <= 0)
            upper.pop_back();
        upper.push_back(p);
    }
    lower.pop_back(); upper.pop_back();
    lower.insert(lower.end(), upper.begin(), upper.end());
    return lower; // points in CCW order
}

def convex_hull(points):
    points = sorted(set(points))
    if len(points) <= 1:
        return points
    def cross(o, a, b):
        return (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])
    lower = []
    for p in points:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        lower.append(p)
    upper = []
    for p in reversed(points):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    return lower[:-1] + upper[:-1]  # CCW hull

💡 Tip: For integer coordinates, use 64‑bit integer arithmetic in the orientation test to avoid overflow.

Line Segment Intersection (Sweep Line)

The Bentley‑Ottmann algorithm finds all intersections among n segments in O((n + k) log n) time, where k is the number of intersections.

The algorithm uses an event queue (ordered by x‑coordinate) and a balanced binary search tree representing the active sweep line.

Insert segment endpoints as events.
When a left endpoint is processed, insert the segment into the active set and check neighbours for intersections.
When a right endpoint is processed, remove the segment and re‑check neighbours.
When an intersection event is encountered, report it and swap the order of the two intersecting segments in the active set.

📝 Note: Implementing a robust sweep line requires careful handling of vertical segments and coincident endpoints. Libraries such as CGAL or Boost.Geometry already provide tested implementations.

Advanced Topics

Voronoi Diagrams & Delaunay Triangulation

A Voronoi diagram partitions the plane into cells where each cell contains the points closest to a particular site. Its dual graph is the Delaunay triangulation, which maximizes the minimum angle of all triangles.

cpp

// Using CGAL (C++)
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using DT = CGAL::Delaunay_triangulation_2<K>;

int main(){
    std::vector<K::Point_2> pts = { {0,0}, {1,0}, {0,1}, {1,1}, {0.5,0.5} };
    DT dt;
    dt.insert(pts.begin(), pts.end());
    for(auto f = dt.finite_faces_begin(); f != dt.finite_faces_end(); ++f){
        // each face is a Delaunay triangle
        std::cout << f->vertex(0)->point() << " "
                  << f->vertex(1)->point() << " "
                  << f->vertex(2)->point() << "\n";
    }
    return 0;
}

💡 Tip: If you need a pure‑Python solution, the scipy.spatial.Voronoi class provides a convenient wrapper around Qhull.

Binary Space Partition (BSP) Trees

BSP trees recursively subdivide space using hyperplanes. They are widely used in rendering pipelines and collision detection for static scenes.

cpp

// Minimal BSP node (C++)
struct BSPNode {
    Segment splitter; // line that splits the space
    BSPNode *front, *back; // subspaces
    std::vector<Segment> polygons; // coplanar polygons stored at leaf
    BSPNode(const Segment& s) : splitter(s), front(nullptr), back(nullptr) {}
};

⚠ Warning: Recursive BSP construction can lead to O(n²) depth on pathological input (e.g., all segments collinear). Apply heuristics like choosing median splitters to keep the tree balanced.

Implementation Tips & Best Practices

Prefer immutable geometric primitives (e.g., const Point) to avoid accidental modifications during complex traversals.
When possible, work with integer arithmetic and rational numbers to guarantee exact predicates; use libraries like Boost::multiprecision if overflow is a concern.
Cache expensive operations such as orientation or distance calculations if they are reused many times within a loop.
Use robust geometry libraries (CGAL, Boost.Geometry, Shapely) for production code— they handle edge cases, EPS tolerances, and provide tested data structures.

💡 Tip: Profile your geometric code with realistic data sets. Many algorithms appear O(n log n) theoretically but degrade to O(n²) due to poor cache locality or excessive dynamic allocations.

Summary

📘 Summary: Computational geometry equips software engineers with powerful tools to solve spatial problems efficiently. By mastering fundamental predicates, core data structures (KD‑Tree, segment tree, BSP), and classic algorithms (convex hull, sweep line, Delaunay triangulation), you can build robust, high‑performance geometry pipelines for graphics, GIS, robotics, and more.

Q: When should I choose a KD‑Tree over a Quadtree?
A: KD‑Trees work well for balanced point queries in any dimension and provide logarithmic nearest‑neighbor searches. Quadtrees are preferable when you need hierarchical spatial subdivision aligned with axis‑aligned regions (e.g., image tiles) or when the dataset is highly non‑uniform.

Q: Is floating‑point arithmetic ever safe for orientation tests?
A: Only if you apply a tolerance (epsilon) and ensure the magnitude of coordinates does not cause catastrophic cancellation. For critical applications (e.g., CAD), use exact arithmetic or symbolic predicates.

Q: Can I reuse a convex hull implementation for 3‑D points?
A: No. The 3‑D convex hull problem requires algorithms like Quickhull or incremental hull construction, which handle facets (triangles) instead of edges.

Q. What is the time complexity of Andrew’s Monotone Chain algorithm for n points?

O(n)
O(n log n)
O(n²)
O(log n)

Answer: O(n log n)
Sorting the points dominates the runtime (O(n log n)), while the hull construction itself is linear.

Q. Which data structure is most appropriate for answering orthogonal range queries in 2‑D?

KD‑Tree
Segment Tree
Fenwick Tree
Trie

Answer: KD‑Tree
KD‑Trees support axis‑aligned range queries with O(√n + k) average time, where k is the number of reported points.

References

Tags: computational-geometry algorithms software-engineering data-structures C++ Python

#ad

Table of Contents

Fundamental Concepts

Coordinate Systems & Orientation Tests

Core Data Structures

KD‑Tree

Segment Tree for 1‑D Intervals

Classic Algorithms

Convex Hull

Andrew’s Monotone Chain (C++)

Line Segment Intersection (Sweep Line)

Advanced Topics

Voronoi Diagrams & Delaunay Triangulation

Binary Space Partition (BSP) Trees

Implementation Tips & Best Practices

Summary

Q. What is the time complexity of Andrew’s Monotone Chain algorithm for n points?

Q. Which data structure is most appropriate for answering orthogonal range queries in 2‑D?

References

Related

Latest articles