🏗️ Fortune's Algorithm Fundamentals
🎯 What is Fortune's Algorithm?
Fortune's Algorithm is a sweep line algorithm for computing Voronoi diagrams in O(n log n) time. It was published by Steven Fortune in 1986 and is considered one of the most elegant algorithms in computational geometry.
Key concept: The algorithm uses a sweepline that moves across the plane, maintaining a beachline of parabolic arcs that represent the Voronoi diagram.
Deep Explanation
Fortune's algorithm is a landmark in computational geometry. It demonstrates how to compute Voronoi diagrams efficiently using a combination of sweepline technique, parabolic curves, and event-driven processing.
Why it matters: Voronoi diagrams have applications in computer graphics, robotics, geography, and computational biology. Fortune's algorithm makes them practical for large datasets.
# Fortune's algorithm (simplified structure)
class FortunesAlgorithm:
def __init__(self, sites):
self.sites = sites # List of (x, y) points
self.events = [] # Priority queue of events
self.beachline = Beachline()
self.voronoi = VoronoiDiagram()
def compute(self):
# Initialize events with all sites
for site in self.sites:
self.events.append(SiteEvent(site))
# Process events in order
while self.events:
event = heapq.heappop(self.events)
if isinstance(event, SiteEvent):
self.handle_site_event(event)
else:
self.handle_circle_event(event)
return self.voronoi
def handle_site_event(self, event):
# Add new parabolic arc to beachline
pass
def handle_circle_event(self, event):
# Remove arc and finalize Voronoi vertex
pass📐 Voronoi Diagrams
🎯 What is a Voronoi Diagram?
A Voronoi diagram partitions a plane into regions based on distance to a set of points (sites). Each region contains all points closer to its site than to any other.
Voronoi Properties
Cells: Each site has a convex polygon cell. Edges: Perpendicular bisectors between sites. Vertices: Intersections of edges. Delaunay: Dual graph of Voronoi diagram.
Generalization: Weighted Voronoi diagrams, higher-order Voronoi diagrams.
# Voronoi diagram representation
class VoronoiDiagram:
def __init__(self):
self.vertices = [] # List of (x, y) points
self.edges = [] # List of (start_idx, end_idx, site1, site2)
self.cells = {} # site_id -> list of edge indices
def add_vertex(self, x, y):
self.vertices.append((x, y))
return len(self.vertices) - 1
def add_edge(self, start, end, site1, site2):
edge = (start, end, site1, site2)
self.edges.append(edge)
# Add to cells
if site1 not in self.cells:
self.cells[site1] = []
if site2 not in self.cells:
self.cells[site2] = []
self.cells[site1].append(len(self.edges) - 1)
self.cells[site2].append(len(self.edges) - 1)
return len(self.edges) - 1
def get_cell_vertices(self, site_id):
# Get polygon vertices for a cell
edges = self.cells.get(site_id, [])
vertices = []
for edge_idx in edges:
start, end, _, _ = self.edges[edge_idx]
vertices.append(self.vertices[start])
return vertices📏 Sweepline Technique
🎯 The Sweepline Concept
A sweepline is a vertical line that moves from left to right across the plane. The algorithm maintains the state of the diagram above the sweepline.
Sweepline Mechanics
Invariant: The Voronoi diagram is complete for sites to the left of the sweepline. State: The beachline represents the boundary between completed and uncompleted regions. Events: Site events (new site encountered) and circle events (vertex formed).
Advantage: Reduces a 2D problem to a 1D sweep.
# Sweepline implementation
class Sweepline:
def __init__(self, x_position=0):
self.x = x_position
self.y = 0
def move_to(self, x):
self.x = x
def distance_to_site(self, site, point):
# Distance from point to site (Euclidean)
return math.sqrt((point[0] - site[0])**2 + (point[1] - site[1])**2)
def site_at_point(self, point, sites):
# Find closest site to point (for the sweepline)
min_dist = float('inf')
closest = None
for site in sites:
dist = self.distance_to_site(site, point)
if dist < min_dist:
min_dist = dist
closest = site
return closest🏖️ Beachline
🎯 The Beachline
The beachline is the set of parabolic arcs that form the boundary between the Voronoi regions above and below the sweepline. It changes as the sweepline moves.
Beachline Properties
Parabolic arcs: Each site contributes a parabolic arc that is equidistant from the site and the sweepline. Breakpoints: Intersections between adjacent arcs. Evolution: Arcs grow and shrink as the sweepline moves.
Data structure: A balanced binary tree of arcs.
# Beachline data structure
class Beachline:
def __init__(self):
self.arcs = None # Root of binary tree
self.breakpoints = {} # Dictionary mapping arcs to breakpoints
class Arc:
def __init__(self, site):
self.site = site
self.left = None
self.right = None
self.parent = None
self.is_leaf = True
def insert_arc(self, site, position):
# Insert a new parabolic arc into the beachline
pass
def remove_arc(self, arc):
# Remove an arc when it disappears
pass
def find_arc_at(self, x, y):
# Find the arc containing a point
pass
def get_breakpoint(self, arc1, arc2):
# Calculate the breakpoint between two arcs
pass🔄 Parabolas
🎯 Parabolic Arcs
Each site generates a parabolic arc that is the set of points equidistant from the site and the sweepline. The arc changes as the sweepline moves.
Parabola Geometry
Equation: For site (a,b) and sweepline at x = s, the parabola is: (x-a)² + (y-b)² = (x-s)². Focus: The site. Directrix: The sweepline. Shape: The arc's curvature changes with distance.
Intersection: Breakpoints are where two parabolas intersect.
# Parabola calculations
def parabola_y(site, sweepline_x, point_x):
# Calculate y on parabola for given x
# Parabola: (x - a)² + (y - b)² = (x - s)²
a, b = site
s = sweepline_x
# Solve for y: (y - b)² = (x - s)² - (x - a)²
# y = b ± sqrt((x - s)² - (x - a)²)
inner = (point_x - s)**2 - (point_x - a)**2
if inner < 0:
return None
return b + math.sqrt(inner)
def parabola_intersection(site1, site2, sweepline_x):
# Find intersection point of two parabolas
# This is where the breakpoint occurs
a1, b1 = site1
a2, b2 = site2
s = sweepline_x
# Solve for x where parabolas intersect
# (x - a1)² + (y - b1)² = (x - a2)² + (y - b2)²
# This reduces to a linear equation in x and y
# Use quadratic formula for exact solution
pass📊 Events
🎯 Event-Driven Processing
Fortune's algorithm is event-driven. There are two types of events: site events and circle events.
Event Types
Site events: Occur when the sweepline passes a site. A new arc is added to the beachline. Circle events: Occur when three arcs meet at a point. An arc is removed and a Voronoi vertex is created.
Event queue: Priority queue ordered by x-coordinate.
# Event classes
class Event:
def __init__(self, x, y):
self.x = x
self.y = y
self.is_valid = True
def __lt__(self, other):
# For priority queue ordering
if self.x != other.x:
return self.x < other.x
return self.y < other.y
class SiteEvent(Event):
def __init__(self, site):
super().__init__(site[0], site[1])
self.site = site
class CircleEvent(Event):
def __init__(self, x, y, arc):
super().__init__(x, y)
self.arc = arc
self.vertex = None
def get_center(self):
# Calculate the center of the circle
# This becomes a Voronoi vertex
pass📍 Site Events
🎯 Handling Site Events
When the sweepline reaches a site, a new parabolic arc is inserted into the beachline, splitting an existing arc.
Site Event Process
Step 1: Find the arc above the new site. Step 2: Split the arc into two parts. Step 3: Insert the new arc between them. Step 4: Add new circle events (if applicable).
Complexity: O(log n) per event.
# Site event handling
def handle_site_event(event, beachline, events_queue):
site = event.site
# Find the arc that contains the site
arc = beachline.find_arc_at(site[0], site[1])
if arc is None:
# First site: create initial arc
beachline.arcs = Beachline.Arc(site)
return
# Split the arc
left_arc = arc
right_arc = Beachline.Arc(arc.site)
new_arc = Beachline.Arc(site)
# Replace arc with left_arc -> new_arc -> right_arc
# (Implementation depends on tree structure)
# Check for new circle events
# A new circle event may occur with left_arc, new_arc, right_arc
# Also with arcs adjacent to left_arc and right_arc
# Add new circle events to queue
for triple in get_adjacent_triples(left_arc, new_arc, right_arc):
if triple_forms_circle(triple):
circle_event = create_circle_event(triple)
events_queue.push(circle_event)⭕ Circle Events
🎯 Circle Events
When three arcs converge at a point, a circle event occurs. The middle arc is removed and a Voronoi vertex is created.
Circle Event Process
Condition: Three consecutive arcs have a circle that touches the sweepline at the bottom. Action: Remove the middle arc. Result: A Voronoi vertex at the circle center.
Validation: Circle events can become invalid if the beachline changes before they occur.
# Circle event handling
def handle_circle_event(event, beachline, voronoi, events_queue):
# Check if event is still valid
if not event.is_valid:
return
arc = event.arc
left_arc = arc.left
right_arc = arc.right
# Create Voronoi vertex at circle center
center = event.get_center()
vertex_idx = voronoi.add_vertex(center[0], center[1])
# Add edges from this vertex to adjacent sites
voronoi.add_edge(vertex_idx, event.vertex, left_arc.site, arc.site)
voronoi.add_edge(vertex_idx, event.vertex, arc.site, right_arc.site)
# Remove the arc from beachline
beachline.remove_arc(arc)
# Add new circle events for adjacent triples
# left_arc may now form a circle with its left neighbor and right_arc
# right_arc may form a circle with left_arc and its right neighbor
for triple in get_adjacent_triples(left_arc, right_arc):
if triple_forms_circle(triple):
circle_event = create_circle_event(triple)
events_queue.push(circle_event)📚 Data Structures
🎯 Efficient Implementation
Fortune's algorithm requires efficient data structures for the beachline (balanced tree), event queue (priority queue), and Voronoi diagram (graph).
Data Structure Choices
Beachline: Red-black tree or AVL tree for O(log n) operations. Event queue: Heap for O(log n) insertion/extraction. Voronoi diagram: Adjacency list for vertices and edges.
Memory: O(n) space for all structures.
# Data structures for Fortune's algorithm
class DataStructures:
def __init__(self):
# Beachline: balanced binary search tree
self.beachline = Beachline()
# Event queue: min-heap
self.event_queue = []
# Voronoi diagram: graph representation
self.voronoi = VoronoiDiagram()
# Site lookup: dictionary by site id
self.sites = {}
# Arc references: map arc to site
self.arc_to_site = {}
def add_event(self, event):
heapq.heappush(self.event_queue, event)
def pop_event(self):
return heapq.heappop(self.event_queue)
def is_empty(self):
return len(self.event_queue) == 0
def update_site(self, site_id, arc):
self.sites[site_id] = arc
self.arc_to_site[arc] = site_id⏱️ Complexity
🎯 Time and Space Analysis
Fortune's algorithm runs in O(n log n) time and uses O(n) space, where n is the number of sites.
Complexity Analysis
Time: O(n log n) — each event is processed in O(log n) time. Space: O(n) — the beachline has at most O(n) arcs, and the event queue has O(n) events. Optimality: Lower bound of Ω(n log n) for Voronoi diagrams.
Constants: The algorithm has a small constant factor, making it practical.
# Complexity analysis
def analyze_complexity(n_sites, n_events):
# n_sites: number of input sites
# n_events: number of events processed
# Each event is processed in O(log n)
# The number of events is O(n)
time_complexity = n_events * math.log(n_sites)
# Space: beachline has O(n) arcs
space_complexity = n_sites
print(f"Time: O({n_events} * log({n_sites})) = O(n log n)")
print(f"Space: O({space_complexity}) = O(n)")
return {
'time': time_complexity,
'space': space_complexity,
'n_sites': n_sites,
'n_events': n_events
}💻 Implementation
🎯 Practical Implementation
Implementation considerations for Fortune's algorithm include numerical precision, edge cases, and performance optimization.
Implementation Tips
Numerical precision: Use double precision or rational arithmetic to avoid floating-point errors. Edge cases: Handle duplicate sites, sites with same x-coordinate, and degenerate cases. Optimization: Avoid unnecessary operations and use efficient data structures.
Testing: Validate with known Voronoi diagrams and random datasets.
# Implementation details
class FortunesAlgorithmImplementation:
def __init__(self, sites):
self.sites = sites
self.epsilon = 1e-10
def handle_numerical_issues(self, x, y):
# Handle floating-point precision
if abs(x) < self.epsilon:
x = 0.0
if abs(y) < self.epsilon:
y = 0.0
return x, y
def check_duplicate_sites(self):
# Remove duplicate sites
unique_sites = []
for site in self.sites:
if site not in unique_sites:
unique_sites.append(site)
return unique_sites
def handle_degenerate_cases(self):
# Handle sites with same x-coordinate
# Special handling for vertical sweepline
pass
def validate_result(self, voronoi):
# Validate Voronoi diagram properties
# Check that cells are convex and partition the plane
pass💼 Applications
🎯 Voronoi Diagram Applications
Voronoi diagrams computed by Fortune's algorithm are used in computer graphics, robotics, geography, and many other fields.
Application Domains
Computer Graphics: Voronoi textures, procedural generation. Robotics: Path planning, obstacle avoidance. Geography: Spatial analysis, nearest neighbor queries. Computational Biology: Protein structure analysis.
Machine Learning: k-NN classification, clustering.
# Application examples
class VoronoiApplications:
def nearest_neighbor(self, voronoi, point):
# Find nearest site using Voronoi diagram
cell = voronoi.find_cell(point)
return cell.site
def path_planning(self, voronoi, start, end):
# Use Voronoi diagram for path planning (Voronoi road map)
# Find path along Voronoi edges
start_cell = voronoi.find_cell(start)
end_cell = voronoi.find_cell(end)
path = self.find_path(start_cell, end_cell, voronoi)
return path
def procedural_generation(self, sites, bounds):
# Generate Voronoi-based procedural texture
voronoi = FortunesAlgorithm(sites).compute()
texture = self.rasterize_voronoi(voronoi, bounds)
return texture🔬 Variations
🎯 Generalized Voronoi Diagrams
Extensions of Fortune's algorithm handle weighted sites, higher-order Voronoi diagrams, and different distance metrics.
Algorithm Variations
Weighted Voronoi: Sites have weights affecting distance. Higher-order: k-th order Voronoi diagrams. Manhattan distance: L1 metric instead of Euclidean. Constrained: Voronoi diagrams with obstacles.
Implementation: Each variation requires modifications to the parabolic arcs and event handling.
# Weighted Voronoi
class WeightedSite:
def __init__(self, x, y, weight):
self.x = x
self.y = y
self.weight = weight
def weighted_parabola(site, sweepline_x, point_x):
# Weighted version: distance = Euclidean distance - weight
a, b = site.x, site.y
w = site.weight
s = sweepline_x
# (x - a)² + (y - b)² = (x - s)² + w²
# Solve for y
inner = (point_x - s)**2 - (point_x - a)**2 + w**2
if inner < 0:
return None
return b + math.sqrt(inner)
# Higher-order Voronoi: k nearest sites
def kth_order_parabolas(k, sites, sweepline_x, point_x):
# Need to consider combinations of k sites
pass🚀 Future Directions
🎯 Emerging Research
Research continues on parallel Voronoi algorithms, GPU implementations, and integration with machine learning.
Future Developments
Parallelization: CUDA and distributed implementations. Machine learning: Learning Voronoi diagrams from data. Real-time: Interactive applications with dynamic sites. 3D: Extending to 3D Voronoi diagrams.
Quantum: Quantum algorithms for computational geometry.
# Future: GPU-accelerated Voronoi
def gpu_voronoi(sites, resolution):
# Using CUDA/OpenCL for parallel computation
# Each pixel computes nearest site independently
# This is O(n) per pixel, but parallelized over pixels
# Pseudo-code for GPU kernel
# For each pixel (x, y):
# min_dist = INF
# for site in sites:
# dist = distance(pixel, site)
# if dist < min_dist:
# min_dist = dist
# pixel_color = site.color
pass
# Future: Learning Voronoi from data
def neural_voronoi(points, labels, model):
# Train a neural network to predict Voronoi cells
# Useful for soft Voronoi diagrams and uncertainty
pass