⚙️AlgorithmAdvanced

Voronoi Diagram and Delaunay

Key Points

•
Voronoi diagrams partition the plane so each region contains points closest to one site, while the Delaunay triangulation connects sites whose Voronoi cells touch.
•
Delaunay triangulations are characterized by the empty circumcircle property and they maximize the minimum angle among all triangulations, avoiding skinny triangles.
•
Fortune’s sweep-line algorithm constructs Voronoi diagrams in O(n $lo g$ n) time; its dual gives the Delaunay triangulation with the same bound.
•
In practice, robust implementations rely on exact predicates for orientation and in-circle tests to avoid floating-point errors.
•
A simple and practical way to compute Delaunay is the incremental Bowyer–Watson algorithm, which is easy to code and typically O( $n^{2}$ ) without point location.
•
You can build the Voronoi diagram from a Delaunay triangulation by connecting triangle circumcenters across shared edges and handling hull edges as rays.
•
Delaunay/Voronoi structures power nearest-neighbor search, facility location, path planning, mesh generation, and natural-neighbor interpolation.
•
Competitive programmers should master geometric predicates, degeneracy handling, and adjacency construction to apply these structures in problems around 2200–2600 CF.

Prerequisites

→2D computational geometry basics — You must know points, vectors, dot/cross products, and orientation tests to implement and reason about triangulations.
→Planar graphs and Euler’s formula — Understanding planarity and combinatorial bounds explains why Delaunay/Voronoi have linear size.
→Convex hull — Hull structure determines Voronoi rays and boundary handling; many algorithms rely on hull properties.
→Sweep-line technique — Fortune’s algorithm is a sweep; knowing event queues and balanced trees helps grasp its mechanics.
→Randomization and incremental algorithms — Random insertion improves expected performance and avoids worst-case patterns.
→Exact geometric predicates — Robust orientation and in-circle tests prevent numerical errors that can break the mesh.
→Graph adjacency and duality — Building Voronoi from Delaunay uses edge-to-face duality and adjacency mapping.
→Bounding-box clipping — To output finite Voronoi edges, you must clip rays/segments to a rectangle.
→Time/space complexity analysis — To choose between Fortune’s algorithm and incremental methods, you need complexity trade-offs.
→Point location in triangulations — Walking and search structures make queries and incremental updates efficient.

Detailed Explanation

Tap terms for definitions

01Overview

A Voronoi diagram splits the plane into regions around a set of sites (points) so that every location in a region is closer to its defining site than to any other. Its geometric dual, the Delaunay triangulation, connects two sites if their Voronoi cells share an edge. These structures capture nearest-neighbor relationships globally, turning distance queries and spatial reasoning into combinatorial graph problems. Delaunay triangulations are prized for good triangle shapes (they maximize the smallest angle) and for strong theoretical properties tied to convexity via a lifting map. Voronoi and Delaunay appear in algorithms for nearest neighbor search, facility placement, mesh generation, and interpolation. While Fortune’s sweep-line algorithm constructs Voronoi in O(n \log n), in contests developers often prefer incremental Delaunay (Bowyer–Watson) or edge-flip methods for simplicity. Once you have a Delaunay triangulation, the Voronoi diagram is easy to derive: triangle circumcenters become Voronoi vertices, and shared triangle edges induce Voronoi edges. Understanding robust orientation and in-circle predicates, degeneracy cases (collinear/cocircular), and graph adjacency is essential to implement these reliably.

02Intuition & Analogies

Imagine you and your friends each plant a flag on a huge field. Now, for every spot on the field, ask: whose flag is closest? If you color each spot by its nearest flag-owner, the field breaks into natural provinces—these are Voronoi cells. The borders between provinces are places where you’d be equally willing to walk to either of two flags, and where three provinces meet, you’d be equally close to three flags. Now draw a straight line between any two friends whose provinces actually touch; you get a mesh of triangles. That triangle mesh is the Delaunay triangulation—the friend graph that “respects” which provinces are neighbors. Why is Delaunay special? Think of stretching a tight rubber band around pins stuck at the flag locations, then pulling a thin fabric over a bowl (a paraboloid). The triangular pattern you see when looking from above corresponds to Delaunay triangles—the ones that avoid creating super-skinny slivers. Another intuition: blow a circle that just touches three flags. If no other flag is inside the bubble, those three flags deserve to be connected—they “own” the space around that circle jointly. If a flag sneaks inside, the triangle is suspect and should be flipped into two better ones. For Voronoi edges on the outer boundary, picture walking away from a single triangle’s center along the perpendicular bisector of its hull edge forever—that’s a Voronoi ray to infinity. These mental models lead directly to practical algorithms: add points one-by-one and carve out circles (Bowyer–Watson), or sweep a beach line that tracks where nearest sites change (Fortune’s algorithm).

03Formal Definition

Given a finite set of distinct sites S = {

s_{1}

s_{2}

\dots

s_{n}

}

\subset

R^{2}

, the Voronoi cell of

s_{i}

is V(

s_{i}

) = { x

\in

R^{2}

: \

∣ x - s_{i} ∥

\leq

∣ x - s_{j} ∥

\forall

j }. The collection {V(

s_{i}

)} forms a partition of

R^{2}

into convex (possibly unbounded) polygons whose boundaries are segments and rays of perpendicular bisectors of site pairs. The Delaunay triangulation D(S) is a planar straight-line graph on S where an edge (

s_{i}

s_{j}

) exists if and only if V(

s_{i}

) and V(

s_{j}

) share a Voronoi edge; faces are triangles (assuming general position, i.e., no four cocircular points). Equivalently, D(S) can be characterized by the empty circumcircle property: a triangle

△

s_{i}

s_{j}

s_{k}

is in D(S) if and only if there exists a circle through

s_{i}

s_{j}

s_{k}

that contains no other site in its interior. Among all triangulations of the convex hull of S, D(S) maximizes the minimum angle and minimizes the largest circumcircle radius, yielding favorable mesh quality. D(S) is also the projection of the lower convex hull of lifted points (x, y)

\mapsto

(x, y,

x^{2}

y^{2}

) on the paraboloid

z = x^{2} + y^{2}

. Under mild nondegeneracy, D(S) is a triangulation with O(n) edges and faces and can be computed in O(n

lo g

n).

04When to Use

Use Voronoi diagrams when nearest-site equivalence regions drive your objective: facility location (assigning customers to warehouses), influence zones (cell towers), motion planning (medial-axis approximations), and interpolation (natural neighbor). Use Delaunay triangulation when you need a high-quality mesh for numerical simulation, terrain modeling (TINs), or to build efficient neighborhood graphs for proximity queries. In competitive programming, Delaunay is handy to answer nearest neighbor queries (via Voronoi duality or triangulation walking), to construct spanners, to detect skinny-triangle pathologies, and to turn geometric optimization into local flips. Prefer Fortune’s algorithm for optimal asymptotics or when you need an explicit Voronoi structure with guaranteed O(n \log n), though it is implementation-heavy. Prefer incremental Delaunay (Bowyer–Watson) or edge-flip Delaunayization when coding speed and clarity matter, as they are concise and robust enough for typical input sizes. When dealing with weighted distances or additively weighted sites, consider the power diagram and regular (weighted Delaunay) triangulation as generalizations.

⚠️Common Mistakes

Ignoring degeneracies: collinear points prevent triangulation; cocircular points create non-unique Delaunay edges. Resolve by symbolic perturbation, deterministic tie-breaking, or exact predicates that treat zero consistently.
Floating-point robustness: naive orientation and in-circle tests with doubles can flip signs near zero, causing nonplanar graphs or infinite loops. Use adaptive exact predicates (e.g., Shewchuk’s) or eps-guards and consistent tie-breaking.
Wrong edge orientation: many algorithms assume counterclockwise triangle orientation; mixing CW/CCW breaks adjacency and Voronoi construction. Normalize triangles when creating them.
Forgetting hull handling: Voronoi edges along the convex hull are rays, not segments. Either clip to a bounding box or output rays explicitly.
Complexity assumptions: incremental Bowyer–Watson without point location is O(n^2) on average; don’t assume O(n \log n) unless you add a search structure (e.g., walking locate or quadtree).
Building Voronoi directly from sites: without Delaunay adjacency, you might incorrectly connect circumcenters that don’t share a primal edge. Always use edge adjacency between triangles to connect circumcenters.
Unbounded growth of super-triangle influence: forgetting to remove triangles that touch super-triangle vertices yields incorrect output. Always filter them at the end.

Key Formulas

Orientation (2x2 determinant)

orient (a, b, c) = det (b_{x} - a_{x} c_{x} - a_{x} b_{y} - a_{y} c_{y} - a_{y}) = (b_{x} - a_{x}) (c_{y} - a_{y}) - (b_{y} - a_{y}) (c_{x} - a_{x})

Explanation: The sign indicates whether $a \to b \to c$ makes a left turn (>0), right turn (<0), or is collinear (=0). It is the core predicate for consistent triangle orientation and point-in-triangle tests.

In-circle (4x4 determinant)

InCircle_{△ ab c} (p) = det a_{x} b_{x} c_{x} p_{x} a_{y} b_{y} c_{y} p_{y} a_{x}^{2} + a_{y}^{2} b_{x}^{2} + b_{y}^{2} c_{x}^{2} + c_{y}^{2} p_{x}^{2} + p_{y}^{2} 1111

Explanation: For CCW triangle abc, the value is >0 if p is strictly inside the circumcircle, =0 if cocircular, and <0 otherwise. This decides cavity triangles and edge flips.

Empty circumcircle property

△ ab c \in Delaunay (S) ⟺ ∄ p \in S ∖ {a, b, c} with InCircle_{△ ab c} (p) > 0

Explanation: A triangle belongs to the Delaunay triangulation exactly when no other site lies inside its circumcircle. This is the defining geometric condition.

Voronoi cell definition

V (s_{i}) = {x : ∥ x - s_{i} ∥ \leq ∥ x - s_{j} ∥, \forall j}, \partial V (s_{i}) \subseteq j ⋃ {x : ∥ x - s_{i} ∥ = ∥ x - s_{j} ∥}

Explanation: A Voronoi cell consists of points at least as close to $s_{i}$ as to any other site, and its boundary lies on perpendicular bisectors of site pairs.

Paraboloid lifting

(x, y) \mapsto (x, y, x^{2} + y^{2}), Delaunay (S) = π_{x y} (lower hull of lifted S)

Explanation: Mapping points onto a paraboloid turns Delaunay triangulation into the projection of the lower convex hull. This explains convexity, duality, and complexity bounds.

Euler’s formula (planar)

V - E + F = 2

Explanation: For a connected planar embedding, the numbers of vertices, edges, and faces satisfy this relation. It bounds the size of planar structures like Delaunay/Voronoi to O(n).

Sweep-line complexity

Fortune’s algorithm: T (n) = O (n lo g n), S (n) = O (n)

Explanation: Maintaining a balanced search tree for the beach line and a priority queue for events gives O(n log n) time and linear space for Voronoi construction.

Max–min angle / min circumradius

T max θ \in T min θ = T min △ \in T max R_{△}

Explanation: Among all triangulations of a point set, Delaunay maximizes the smallest triangle angle and equivalently minimizes the largest circumradius, avoiding skinny triangles.

Power distance (weighted)

π_{w} (x, s) = ∥ x - s ∥^{2} - w_{s}

Explanation: Power diagrams replace Euclidean distance with power distance; their dual is the regular triangulation. Useful for weighted sites and additive biases.

Planar duality (edge correspondence)

Voronoi edge between s_{i}, s_{j} ⟺ (s_{i}, s_{j}) \in Delaunay edges

Explanation: Two sites share a Voronoi edge exactly when they are connected in the Delaunay triangulation. This bijection lets you derive Voronoi from Delaunay and vice versa.

Complexity Analysis

Fortune’s algorithm achieves O(n log n) time by sweeping a vertical line across sites while maintaining a balanced-tree representation of the beach line and a priority queue of site/circle events; each event is processed in O(log n) amortized time, and there are O(n) events. Space is O(n) for arcs, edges, and events. Implementing it robustly requires careful numerical handling and event invalidation. The incremental Bowyer–Watson algorithm is simpler: insert points one-by-one, delete triangles whose circumcircles contain the new point (the cavity), and stitch the boundary to the point. If triangles are searched naively, each insertion may inspect O(m) triangles (m current triangles), yielding O(

n^{2}

) time overall in practice. With randomized insertion and a point-location structure (e.g., walking from the last located triangle or a history DAG), expected time improves toward O(n log n). Space usage is O(n) for triangles and adjacency. Building the Voronoi diagram from Delaunay adjacency is linear in the number of triangles/edges: computing circumcenters is O(1) per triangle; connecting adjacent triangle circumcenters is O(1) per edge. Handling hull edges produces rays that can be clipped to a bounding box in O(1) each. Thus, after triangulation, Voronoi construction is O(n) time and space. For queries, walking point location on a well-shaped Delaunay triangulation typically visits O(

n

) triangles on random data, with worst-case O(n); space remains O(n) for the triangulation and O(1) per query. Edge-flip Delaunayization from an arbitrary triangulation completes in O(

n^{2}

) worst-case but often runs near-linear for well-behaved inputs.

Code Examples

Incremental Delaunay Triangulation (Bowyer–Watson)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct Point {
5     double x, y;
6 };
7 
8 static const double EPS = 1e-12;
9 
10 // Cross product (b - a) x (c - a)
11 double orient(const Point& a, const Point& b, const Point& c) {
12     return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
13 }
14 
15 // Squared Euclidean distance
16 static inline double dist2(const Point& a, const Point& b) {
17     double dx = a.x - b.x, dy = a.y - b.y;
18     return dx*dx + dy*dy;
19 }
20 
21 // Circumcenter of triangle abc (assumes non-collinear)
22 Point circumcenter(const Point& a, const Point& b, const Point& c) {
23     double d = 2.0 * orient(a, b, c);
24     // If nearly collinear, fall back to a large but finite center (avoid division by ~0)
25     if (fabs(d) < EPS) {
26         // Use midpoint of the longest edge as a degenerate proxy
27         double dab = dist2(a,b), dbc = dist2(b,c), dca = dist2(c,a);
28         if (dab >= dbc && dab >= dca) return Point{(a.x + b.x)/2.0, (a.y + b.y)/2.0};
29         else if (dbc >= dab && dbc >= dca) return Point{(b.x + c.x)/2.0, (b.y + c.y)/2.0};
30         else return Point{(c.x + a.x)/2.0, (c.y + a.y)/2.0};
31     }
32     double ax = a.x, ay = a.y, bx = b.x, by = b.y, cx = c.x, cy = c.y;
33     double a2 = ax*ax + ay*ay;
34     double b2 = bx*bx + by*by;
35     double c2 = cx*cx + cy*cy;
36     double ux = (a2*(by - cy) + b2*(cy - ay) + c2*(ay - by)) / d;
37     double uy = (a2*(cx - bx) + b2*(ax - cx) + c2*(bx - ax)) / d;
38     return Point{ux, uy};
39 }
40 
41 // In-circle test using orientation-consistent sign
42 // Returns >0 if p is inside circumcircle of CCW(a,b,c)
43 double inCircle(const Point& a, const Point& b, const Point& c, const Point& p) {
44     // Translate so p is the origin to improve stability
45     double ax = a.x - p.x, ay = a.y - p.y;
46     double bx = b.x - p.x, by = b.y - p.y;
47     double cx = c.x - p.x, cy = c.y - p.y;
48     double det = (ax*ax + ay*ay) * (bx*cy - by*cx)
49                - (bx*bx + by*by) * (ax*cy - ay*cx)
50                + (cx*cx + cy*cy) * (ax*by - ay*bx);
51     // Ensure (a,b,c) is CCW for the sign convention
52     double s = orient(a,b,c);
53     if (s < 0) det = -det;
54     return det; // >0 => inside, <0 => outside, ~0 => cocircular
55 }
56 
57 struct Tri { int a, b, c; }; // CCW indices into P
58 
59 // Normalize triangle orientation to CCW
60 void make_ccw(const vector<Point>& P, Tri& t) {
61     if (orient(P[t.a], P[t.b], P[t.c]) < 0) swap(t.b, t.c);
62 }
63 
64 // Delaunay triangulation via Bowyer-Watson; returns triangles over original points
65 vector<Tri> delaunay_bowyer_watson(vector<Point> P) {
66     int n = (int)P.size();
67     // Build super-triangle that contains all points
68     double minx = P[0].x, miny = P[0].y, maxx = P[0].x, maxy = P[0].y;
69     for (int i = 1; i < n; ++i) {
70         minx = min(minx, P[i].x); miny = min(miny, P[i].y);
71         maxx = max(maxx, P[i].x); maxy = max(maxy, P[i].y);
72     }
73     double dx = maxx - minx, dy = maxy - miny;
74     double delta = max(dx, dy);
75     double cx = (minx + maxx) / 2.0, cy = (miny + maxy) / 2.0;
76 
77     // Append 3 super-triangle points
78     Point s1{cx - 2*delta, cy - delta};
79     Point s2{cx,            cy + 2*delta};
80     Point s3{cx + 2*delta,  cy - delta};
81     int sidx1 = n, sidx2 = n+1, sidx3 = n+2;
82 
83     vector<Point> Q = P; // work on extended array
84     Q.push_back(s1); Q.push_back(s2); Q.push_back(s3);
85 
86     vector<Tri> tris;
87     tris.push_back({sidx1, sidx2, sidx3}); // initial super triangle
88 
89     // Randomize insertion order to improve behavior
90     vector<int> order(n);
91     iota(order.begin(), order.end(), 0);
92     mt19937 rng(712367);
93     shuffle(order.begin(), order.end(), rng);
94 
95     for (int pi : order) {
96         const Point& p = Q[pi];
97         // 1) Find all triangles whose circumcircle contains p
98         vector<char> bad(tris.size(), 0);
99         for (size_t t = 0; t < tris.size(); ++t) {
100             const Tri& T = tris[t];
101             double ic = inCircle(Q[T.a], Q[T.b], Q[T.c], p);
102             if (ic > EPS) bad[t] = 1; // inside => remove
103         }
104         // 2) Collect boundary edges of the cavity
105         struct DEdge { int u, v; };
106         vector<DEdge> directedEdges; directedEdges.reserve(tris.size()*3);
107         unordered_map<long long,int> cnt; cnt.reserve(tris.size()*3*2);
108         auto key = [&](int u, int v){ if(u>v) swap(u,v); return ( (long long)u<<32 ) | (unsigned)v; };
109 
110         for (size_t t = 0; t < tris.size(); ++t) if (bad[t]) {
111             Tri T = tris[t]; make_ccw(Q, T);
112             int v[3] = {T.a, T.b, T.c};
113             for (int e = 0; e < 3; ++e) {
114                 int u = v[e], w = v[(e+1)%3];
115                 directedEdges.push_back({u, w});
116                 cnt[key(u,w)]++;
117             }
118         }
119         // 3) Remove bad triangles
120         vector<Tri> keep; keep.reserve(tris.size());
121         for (size_t t = 0; t < tris.size(); ++t) if (!bad[t]) keep.push_back(tris[t]);
122         tris.swap(keep);
123         // 4) For each boundary edge (appears exactly once as undirected), form new triangle with p
124         for (auto e : directedEdges) {
125             if (cnt[key(e.u, e.v)] == 1) {
126                 Tri nt{e.u, e.v, pi};
127                 make_ccw(Q, nt);
128                 tris.push_back(nt);
129             }
130         }
131     }
132 
133     // Remove triangles that use any super vertex
134     vector<Tri> out; out.reserve(tris.size());
135     for (const Tri& T : tris) {
136         if (T.a < n && T.b < n && T.c < n) {
137             Tri TT = T; make_ccw(P, const_cast<Tri&>(TT));
138             out.push_back(TT);
139         }
140     }
141     return out;
142 }
143 
144 int main() {
145     ios::sync_with_stdio(false);
146     cin.tie(nullptr);
147 
148     int n; cin >> n;
149     vector<Point> P(n);
150     for (int i = 0; i < n; ++i) cin >> P[i].x >> P[i].y;
151 
152     auto tris = delaunay_bowyer_watson(P);
153 
154     cout << tris.size() << "\n";
155     for (auto &t : tris) {
156         cout << t.a << " " << t.b << " " << t.c << "\n"; // 0-based indices of input points
157     }
158     return 0;
159 }
160

This program builds a Delaunay triangulation using the Bowyer–Watson incremental algorithm. It starts with a super-triangle that encloses all points, inserts points in random order, removes triangles whose circumcircles contain the new point (the cavity), then creates new triangles from the cavity boundary to the point. Finally, it discards triangles that involve super-triangle vertices and outputs the remaining triangles as 0-based indices into the input points.

Time: O(n^2) expected with naive cavity search; O(n log n) expected if augmented with efficient point location.Space: O(n) triangles and auxiliary structures.

Build Voronoi Diagram (segments and rays) from Delaunay Triangulation

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct Point { double x, y; };
5 struct Tri { int a, b, c; }; // CCW
6 struct Seg { Point s, t; };   // Voronoi segment (finite)
7 struct Ray { Point s, d; };   // Voronoi ray: start s, direction d (unit-ish)
8 
9 static const double EPS = 1e-12;
10 
11 double orient(const Point& a, const Point& b, const Point& c) {
12     return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
13 }
14 
15 Point circumcenter(const Point& a, const Point& b, const Point& c) {
16     double d = 2.0 * orient(a, b, c);
17     if (fabs(d) < EPS) {
18         // Fallback: use midpoint of the longest edge
19         auto d2 = [&](const Point& p, const Point& q){ double dx=p.x-q.x, dy=p.y-q.y; return dx*dx+dy*dy; };
20         double dab = d2(a,b), dbc = d2(b,c), dca = d2(c,a);
21         if (dab >= dbc && dab >= dca) return Point{(a.x+b.x)/2.0, (a.y+b.y)/2.0};
22         else if (dbc >= dab && dbc >= dca) return Point{(b.x+c.x)/2.0, (b.y+c.y)/2.0};
23         else return Point{(c.x+a.x)/2.0, (c.y+a.y)/2.0};
24     }
25     double ax = a.x, ay = a.y, bx = b.x, by = b.y, cx = c.x, cy = c.y;
26     double a2 = ax*ax + ay*ay;
27     double b2 = bx*bx + by*by;
28     double c2 = cx*cx + cy*cy;
29     double ux = (a2*(by - cy) + b2*(cy - ay) + c2*(ay - by)) / d;
30     double uy = (a2*(cx - bx) + b2*(ax - cx) + c2*(bx - ax)) / d;
31     return Point{ux, uy};
32 }
33 
34 // Intersect a ray p + t*d (t >= 0) with axis-aligned bounding box [xmin,xmax]x[ymin,ymax]
35 // Returns the closest intersection point; if no intersection, returns p.
36 Point intersectRayAABB(Point p, Point d, double xmin, double ymin, double xmax, double ymax) {
37     // Normalize direction to avoid huge t; keep zeros if tiny
38     double len = hypot(d.x, d.y);
39     if (len < EPS) return p;
40     d.x /= len; d.y /= len;
41 
42     vector<Point> candidates;
43     auto push_if = [&](double t){ if (t >= 0) candidates.push_back(Point{p.x + t*d.x, p.y + t*d.y}); };
44 
45     // Left x = xmin
46     if (fabs(d.x) > EPS) {
47         double t = (xmin - p.x) / d.x; double y = p.y + t*d.y;
48         if (t >= 0 && y >= ymin - 1e-9 && y <= ymax + 1e-9) push_if(t);
49     }
50     // Right x = xmax
51     if (fabs(d.x) > EPS) {
52         double t = (xmax - p.x) / d.x; double y = p.y + t*d.y;
53         if (t >= 0 && y >= ymin - 1e-9 && y <= ymax + 1e-9) push_if(t);
54     }
55     // Bottom y = ymin
56     if (fabs(d.y) > EPS) {
57         double t = (ymin - p.y) / d.y; double x = p.x + t*d.x;
58         if (t >= 0 && x >= xmin - 1e-9 && x <= xmax + 1e-9) push_if(t);
59     }
60     // Top y = ymax
61     if (fabs(d.y) > EPS) {
62         double t = (ymax - p.y) / d.y; double x = p.x + t*d.x;
63         if (t >= 0 && x >= xmin - 1e-9 && x <= xmax + 1e-9) push_if(t);
64     }
65 
66     if (candidates.empty()) return p;
67     // Choose the nearest intersection
68     auto d2 = [&](const Point& a, const Point& b){ double dx=a.x-b.x, dy=a.y-b.y; return dx*dx+dy*dy; };
69     int best = 0; double bestd = d2(p, candidates[0]);
70     for (int i = 1; i < (int)candidates.size(); ++i) {
71         double di = d2(p, candidates[i]); if (di < bestd) bestd = di, best = i;
72     }
73     return candidates[best];
74 }
75 
76 // Helper to ensure CCW triangle
77 void make_ccw(const vector<Point>& P, Tri& t) {
78     if (orient(P[t.a], P[t.b], P[t.c]) < 0) swap(t.b, t.c);
79 }
80 
81 int main(){
82     ios::sync_with_stdio(false);
83     cin.tie(nullptr);
84 
85     int n; cin >> n;
86     vector<Point> P(n);
87     for (int i = 0; i < n; ++i) cin >> P[i].x >> P[i].y;
88 
89     int m; cin >> m; // number of Delaunay triangles
90     vector<Tri> T(m);
91     for (int i = 0; i < m; ++i) { cin >> T[i].a >> T[i].b >> T[i].c; make_ccw(P, T[i]); }
92 
93     // Compute circumcenters
94     vector<Point> CC(m);
95     for (int i = 0; i < m; ++i) {
96         CC[i] = circumcenter(P[T[i].a], P[T[i].b], P[T[i].c]);
97     }
98 
99     // Build edge -> incident triangle(s)
100     struct EdgeRec { int tri; int u, v; }; // directed as in its triangle
101     struct Key { int x, y; bool operator==(const Key& o) const { return x==o.x && y==o.y; } };
102     struct KeyHash { size_t operator()(const Key& k) const { return (size_t)k.x*1000003u ^ (size_t)k.y; } };
103     unordered_map<Key, vector<EdgeRec>, KeyHash> adj; adj.reserve(m*3*2);
104 
105     auto add_edge = [&](int tri, int u, int v){
106         Key k{min(u,v), max(u,v)}; adj[k].push_back({tri, u, v});
107     };
108 
109     for (int i = 0; i < m; ++i) {
110         add_edge(i, T[i].a, T[i].b);
111         add_edge(i, T[i].b, T[i].c);
112         add_edge(i, T[i].c, T[i].a);
113     }
114 
115     // Bounding box for clipping hull rays
116     double minx = P[0].x, miny = P[0].y, maxx = P[0].x, maxy = P[0].y;
117     for (int i = 1; i < n; ++i) { minx = min(minx, P[i].x); miny = min(miny, P[i].y); maxx = max(maxx, P[i].x); maxy = max(maxy, P[i].y); }
118     double pad = (max(maxx-minx, maxy-miny) + 1.0) * 0.2; // small padding
119     double X0 = minx - pad, Y0 = miny - pad, X1 = maxx + pad, Y1 = maxy + pad;
120 
121     vector<Seg> segs; segs.reserve(m*2);
122     vector<Seg> rays_as_segs; rays_as_segs.reserve(m);
123 
124     for (auto &kv : adj) {
125         auto &vec = kv.second;
126         if (vec.size() == 2) {
127             // Interior edge: Voronoi segment between circumcenters
128             int t1 = vec[0].tri, t2 = vec[1].tri;
129             segs.push_back({CC[t1], CC[t2]});
130         } else if (vec.size() == 1) {
131             // Hull edge: Voronoi ray from circumcenter along outward normal of the perpendicular bisector
132             int t = vec[0].tri; int u = vec[0].u, v = vec[0].v;
133             // Directed edge u->v is oriented as in CCW triangle t; inside is to the left.
134             Point a = P[u], b = P[v];
135             Point dir_edge{b.x - a.x, b.y - a.y};
136             // Left normal is (-dy, dx); outward (right) normal is its negative
137             Point right_normal{dir_edge.y, -dir_edge.x}; // rotate -90° equivalently
138             // The Voronoi edge lies along the perpendicular bisector; direction is right_normal
139             Point start = CC[t];
140             Point end = intersectRayAABB(start, right_normal, X0, Y0, X1, Y1);
141             rays_as_segs.push_back({start, end});
142         }
143     }
144 
145     // Output: first finite segments, then hull rays as segments to the bbox
146     cout << segs.size() << "\n";
147     cout.setf(std::ios::fixed); cout << setprecision(10);
148     for (auto &s : segs) {
149         cout << s.s.x << " " << s.s.y << " " << s.t.x << " " << s.t.y << "\n";
150     }
151     cout << rays_as_segs.size() << "\n";
152     for (auto &s : rays_as_segs) {
153         cout << s.s.x << " " << s.s.y << " " << s.t.x << " " << s.t.y << "\n";
154     }
155     return 0;
156 }
157

Given input points and a Delaunay triangulation (e.g., from Example 1), this program computes triangle circumcenters and then builds Voronoi edges. Each shared Delaunay edge corresponds to a Voronoi segment connecting the two incident triangles’ circumcenters. Each hull edge corresponds to a Voronoi ray, which we clip to a padded bounding box for finite output. The output lists finite Voronoi segments followed by clipped hull rays.

Time: O(m) after triangulation, where m = number of Delaunay triangles (overall O(n)).Space: O(n) for circumcenters and edge adjacency.

Point Location and Nearest Neighbor via Delaunay Walking

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct Point { double x, y; };
5 struct Tri { int v[3]; }; // CCW
6 
7 static const double EPS = 1e-12;
8 
9 double orient(const Point& a, const Point& b, const Point& c) {
10     return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
11 }
12 
13 bool pointInTri(const vector<Point>& P, const Tri& T, const Point& q) {
14     const Point& a = P[T.v[0]], &b = P[T.v[1]], &c = P[T.v[2]];
15     double o1 = orient(a,b,q), o2 = orient(b,c,q), o3 = orient(c,a,q);
16     return (o1 >= -EPS && o2 >= -EPS && o3 >= -EPS);
17 }
18 
19 // Build adjacency: for each triangle, neighbor across each edge (or -1 if hull)
20 vector<array<int,3>> buildNeighbors(const vector<Point>& P, const vector<Tri>& T) {
21     int m = (int)T.size();
22     struct Key { int a,b; bool operator==(const Key& o) const { return a==o.a && b==o.b; } };
23     struct Hash { size_t operator()(const Key& k) const { return (size_t)k.a*1000003u ^ (size_t)k.b; } };
24     unordered_map<Key, pair<int,int>, Hash> mp; mp.reserve(m*3*2);
25     auto norm = [&](int u, int v){ if(u>v) swap(u,v); return Key{u,v}; };
26 
27     for (int i = 0; i < m; ++i) {
28         for (int e = 0; e < 3; ++e) {
29             int u = T[i].v[e], v = T[i].v[(e+1)%3];
30             Key k = norm(u,v);
31             if (!mp.count(k)) mp[k] = {i, -1}; else mp[k].second = i;
32         }
33     }
34     vector<array<int,3>> nei(m); for (auto &a:nei) a = {-1,-1,-1};
35     for (auto &kv : mp) {
36         int t1 = kv.second.first, t2 = kv.second.second;
37         if (t2 == -1) continue; // hull edge
38         // For t1, find which edge matches, set neighbor to t2; and vice versa
39         auto setN = [&](int ti, int tj){
40             for (int e = 0; e < 3; ++e) {
41                 int u = T[ti].v[e], v = T[ti].v[(e+1)%3];
42                 int a = kv.first.a, b = kv.first.b; // undirected
43                 if ((u==a && v==b) || (u==b && v==a)) { nei[ti][e] = tj; break; }
44             }
45         };
46         setN(t1, t2); setN(t2, t1);
47     }
48     return nei;
49 }
50 
51 // Greedy walking point location: start from startTri, walk across edges that the query is to the right of.
52 // Returns index of containing triangle if found, else -1 (outside hull).
53 int locateTriangle(const vector<Point>& P, const vector<Tri>& T, const vector<array<int,3>>& nei, const Point& q, int startTri) {
54     if (startTri < 0 || startTri >= (int)T.size()) startTri = 0;
55     int cur = startTri;
56     int steps = 0, maxSteps = (int)T.size() * 2 + 10;
57     while (steps++ < maxSteps) {
58         const Tri& tr = T[cur];
59         bool inside = true; int edgeToCross = -1;
60         for (int e = 0; e < 3; ++e) {
61             int u = tr.v[e], v = tr.v[(e+1)%3];
62             if (orient(P[u], P[v], q) < -EPS) { inside = false; edgeToCross = e; break; }
63         }
64         if (inside) return cur;
65         int nxt = nei[cur][edgeToCross];
66         if (nxt == -1) return -1; // outside the convex hull
67         cur = nxt;
68     }
69     return -1; // fallback
70 }
71 
72 int main(){
73     ios::sync_with_stdio(false);
74     cin.tie(nullptr);
75 
76     int n; cin >> n;
77     vector<Point> P(n);
78     for (int i = 0; i < n; ++i) cin >> P[i].x >> P[i].y;
79 
80     int m; cin >> m;
81     vector<Tri> T(m);
82     for (int i = 0; i < m; ++i) cin >> T[i].v[0] >> T[i].v[1] >> T[i].v[2];
83 
84     auto nei = buildNeighbors(P, T);
85 
86     int qn; cin >> qn;
87     cout.setf(std::ios::fixed); cout << setprecision(10);
88 
89     int start = 0; // reuse last triangle for temporal coherence
90     for (int qi = 0; qi < qn; ++qi) {
91         Point q; cin >> q.x >> q.y;
92         int t = locateTriangle(P, T, nei, q, start);
93         int bestIdx = -1; double bestD2 = 1e300;
94         if (t != -1) {
95             // If inside a triangle, the nearest site is one of its vertices (Voronoi duality)
96             for (int k = 0; k < 3; ++k) {
97                 int vi = T[t].v[k];
98                 double dx = q.x - P[vi].x, dy = q.y - P[vi].y; double d2 = dx*dx + dy*dy;
99                 if (d2 < bestD2) { bestD2 = d2; bestIdx = vi; }
100             }
101             start = t; // warm-start next query
102         } else {
103             // Outside hull: fallback to checking convex hull vertices only (simple, robust)
104             // Build hull via Andrew's monotone chain
105             vector<int> id(n); iota(id.begin(), id.end(), 0);
106             sort(id.begin(), id.end(), [&](int a, int b){ if (P[a].x != P[b].x) return P[a].x < P[b].x; return P[a].y < P[b].y; });
107             vector<int> H;
108             for (int phase = 0; phase < 2; ++phase) {
109                 size_t startH = H.size();
110                 for (int idx : id) {
111                     while (H.size() >= startH + 2) {
112                         Point a = P[H[H.size()-2]], b = P[H[H.size()-1]], c = P[idx];
113                         if (orient(a,b,c) <= 0) H.pop_back(); else break;
114                     }
115                     H.push_back(idx);
116                 }
117                 H.pop_back(); reverse(id.begin(), id.end());
118             }
119             for (int vi : H) {
120                 double dx = q.x - P[vi].x, dy = q.y - P[vi].y; double d2 = dx*dx + dy*dy;
121                 if (d2 < bestD2) { bestD2 = d2; bestIdx = vi; }
122             }
123         }
124         cout << bestIdx << "\n"; // nearest site index
125     }
126     return 0;
127 }
128

This program supports nearest-neighbor queries using the Delaunay triangulation. It first builds triangle adjacency, then locates a query point by walking across edges where the point lies to the right, which typically reaches the containing triangle quickly. If the query lies inside a triangle, the nearest site is guaranteed to be one of the triangle’s vertices. If the query lies outside the convex hull, it falls back to checking hull vertices (simple and safe). It reuses the last located triangle as a warm start for subsequent queries.

Time: Amortized O(√n) per query on typical inputs (walking), worst-case O(n); fallback hull check is O(h) where h is hull size.Space: O(n) for triangulation and adjacency; O(1) per query.

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct Point {
5	double x, y;
6	};
7
8	static const double EPS = 1e-12;
9
10	// Cross product (b - a) x (c - a)
11	double orient(const Point& a, const Point& b, const Point& c) {
12	return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
13	}
14
15	// Squared Euclidean distance
16	static inline double dist2(const Point& a, const Point& b) {
17	double dx = a.x - b.x, dy = a.y - b.y;
18	return dxdx + dydy;
19	}
20
21	// Circumcenter of triangle abc (assumes non-collinear)
22	Point circumcenter(const Point& a, const Point& b, const Point& c) {
23	double d = 2.0 * orient(a, b, c);
24	// If nearly collinear, fall back to a large but finite center (avoid division by ~0)
25	if (fabs(d) < EPS) {
26	// Use midpoint of the longest edge as a degenerate proxy
27	double dab = dist2(a,b), dbc = dist2(b,c), dca = dist2(c,a);
28	if (dab >= dbc && dab >= dca) return Point{(a.x + b.x)/2.0, (a.y + b.y)/2.0};
29	else if (dbc >= dab && dbc >= dca) return Point{(b.x + c.x)/2.0, (b.y + c.y)/2.0};
30	else return Point{(c.x + a.x)/2.0, (c.y + a.y)/2.0};
31	}
32	double ax = a.x, ay = a.y, bx = b.x, by = b.y, cx = c.x, cy = c.y;
33	double a2 = axax + ayay;
34	double b2 = bxbx + byby;
35	double c2 = cxcx + cycy;
36	double ux = (a2(by - cy) + b2(cy - ay) + c2*(ay - by)) / d;
37	double uy = (a2(cx - bx) + b2(ax - cx) + c2*(bx - ax)) / d;
38	return Point{ux, uy};
39	}
40
41	// In-circle test using orientation-consistent sign
42	// Returns >0 if p is inside circumcircle of CCW(a,b,c)
43	double inCircle(const Point& a, const Point& b, const Point& c, const Point& p) {
44	// Translate so p is the origin to improve stability
45	double ax = a.x - p.x, ay = a.y - p.y;
46	double bx = b.x - p.x, by = b.y - p.y;
47	double cx = c.x - p.x, cy = c.y - p.y;
48	double det = (axax + ayay) * (bxcy - bycx)
49	- (bxbx + byby) * (axcy - aycx)
50	+ (cxcx + cycy) * (axby - aybx);
51	// Ensure (a,b,c) is CCW for the sign convention
52	double s = orient(a,b,c);
53	if (s < 0) det = -det;
54	return det; // >0 => inside, <0 => outside, ~0 => cocircular
55	}
56
57	struct Tri { int a, b, c; }; // CCW indices into P
58
59	// Normalize triangle orientation to CCW
60	void make_ccw(const vector<Point>& P, Tri& t) {
61	if (orient(P[t.a], P[t.b], P[t.c]) < 0) swap(t.b, t.c);
62	}
63
64	// Delaunay triangulation via Bowyer-Watson; returns triangles over original points
65	vector<Tri> delaunay_bowyer_watson(vector<Point> P) {
66	int n = (int)P.size();
67	// Build super-triangle that contains all points
68	double minx = P[0].x, miny = P[0].y, maxx = P[0].x, maxy = P[0].y;
69	for (int i = 1; i < n; ++i) {
70	minx = min(minx, P[i].x); miny = min(miny, P[i].y);
71	maxx = max(maxx, P[i].x); maxy = max(maxy, P[i].y);
72	}
73	double dx = maxx - minx, dy = maxy - miny;
74	double delta = max(dx, dy);
75	double cx = (minx + maxx) / 2.0, cy = (miny + maxy) / 2.0;
76
77	// Append 3 super-triangle points
78	Point s1{cx - 2*delta, cy - delta};
79	Point s2{cx, cy + 2*delta};
80	Point s3{cx + 2*delta, cy - delta};
81	int sidx1 = n, sidx2 = n+1, sidx3 = n+2;
82
83	vector<Point> Q = P; // work on extended array
84	Q.push_back(s1); Q.push_back(s2); Q.push_back(s3);
85
86	vector<Tri> tris;
87	tris.push_back({sidx1, sidx2, sidx3}); // initial super triangle
88
89	// Randomize insertion order to improve behavior
90	vector<int> order(n);
91	iota(order.begin(), order.end(), 0);
92	mt19937 rng(712367);
93	shuffle(order.begin(), order.end(), rng);
94
95	for (int pi : order) {
96	const Point& p = Q[pi];
97	// 1) Find all triangles whose circumcircle contains p
98	vector<char> bad(tris.size(), 0);
99	for (size_t t = 0; t < tris.size(); ++t) {
100	const Tri& T = tris[t];
101	double ic = inCircle(Q[T.a], Q[T.b], Q[T.c], p);
102	if (ic > EPS) bad[t] = 1; // inside => remove
103	}
104	// 2) Collect boundary edges of the cavity
105	struct DEdge { int u, v; };
106	vector<DEdge> directedEdges; directedEdges.reserve(tris.size()*3);
107	unordered_map<long long,int> cnt; cnt.reserve(tris.size()32);
108	auto key = [&](int u, int v){ if(u>v) swap(u,v); return ( (long long)u<<32 ) \| (unsigned)v; };
109
110	for (size_t t = 0; t < tris.size(); ++t) if (bad[t]) {
111	Tri T = tris[t]; make_ccw(Q, T);
112	int v[3] = {T.a, T.b, T.c};
113	for (int e = 0; e < 3; ++e) {
114	int u = v[e], w = v[(e+1)%3];
115	directedEdges.push_back({u, w});
116	cnt[key(u,w)]++;
117	}
118	}
119	// 3) Remove bad triangles
120	vector<Tri> keep; keep.reserve(tris.size());
121	for (size_t t = 0; t < tris.size(); ++t) if (!bad[t]) keep.push_back(tris[t]);
122	tris.swap(keep);
123	// 4) For each boundary edge (appears exactly once as undirected), form new triangle with p
124	for (auto e : directedEdges) {
125	if (cnt[key(e.u, e.v)] == 1) {
126	Tri nt{e.u, e.v, pi};
127	make_ccw(Q, nt);
128	tris.push_back(nt);
129	}
130	}
131	}
132
133	// Remove triangles that use any super vertex
134	vector<Tri> out; out.reserve(tris.size());
135	for (const Tri& T : tris) {
136	if (T.a < n && T.b < n && T.c < n) {
137	Tri TT = T; make_ccw(P, const_cast<Tri&>(TT));
138	out.push_back(TT);
139	}
140	}
141	return out;
142	}
143
144	int main() {
145	ios::sync_with_stdio(false);
146	cin.tie(nullptr);
147
148	int n; cin >> n;
149	vector<Point> P(n);
150	for (int i = 0; i < n; ++i) cin >> P[i].x >> P[i].y;
151
152	auto tris = delaunay_bowyer_watson(P);
153
154	cout << tris.size() << "\n";
155	for (auto &t : tris) {
156	cout << t.a << " " << t.b << " " << t.c << "\n"; // 0-based indices of input points
157	}
158	return 0;
159	}
160

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct Point { double x, y; };
5	struct Tri { int a, b, c; }; // CCW
6	struct Seg { Point s, t; }; // Voronoi segment (finite)
7	struct Ray { Point s, d; }; // Voronoi ray: start s, direction d (unit-ish)
8
9	static const double EPS = 1e-12;
10
11	double orient(const Point& a, const Point& b, const Point& c) {
12	return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
13	}
14
15	Point circumcenter(const Point& a, const Point& b, const Point& c) {
16	double d = 2.0 * orient(a, b, c);
17	if (fabs(d) < EPS) {
18	// Fallback: use midpoint of the longest edge
19	auto d2 = [&](const Point& p, const Point& q){ double dx=p.x-q.x, dy=p.y-q.y; return dxdx+dydy; };
20	double dab = d2(a,b), dbc = d2(b,c), dca = d2(c,a);
21	if (dab >= dbc && dab >= dca) return Point{(a.x+b.x)/2.0, (a.y+b.y)/2.0};
22	else if (dbc >= dab && dbc >= dca) return Point{(b.x+c.x)/2.0, (b.y+c.y)/2.0};
23	else return Point{(c.x+a.x)/2.0, (c.y+a.y)/2.0};
24	}
25	double ax = a.x, ay = a.y, bx = b.x, by = b.y, cx = c.x, cy = c.y;
26	double a2 = axax + ayay;
27	double b2 = bxbx + byby;
28	double c2 = cxcx + cycy;
29	double ux = (a2(by - cy) + b2(cy - ay) + c2*(ay - by)) / d;
30	double uy = (a2(cx - bx) + b2(ax - cx) + c2*(bx - ax)) / d;
31	return Point{ux, uy};
32	}
33
34	// Intersect a ray p + t*d (t >= 0) with axis-aligned bounding box [xmin,xmax]x[ymin,ymax]
35	// Returns the closest intersection point; if no intersection, returns p.
36	Point intersectRayAABB(Point p, Point d, double xmin, double ymin, double xmax, double ymax) {
37	// Normalize direction to avoid huge t; keep zeros if tiny
38	double len = hypot(d.x, d.y);
39	if (len < EPS) return p;
40	d.x /= len; d.y /= len;
41
42	vector<Point> candidates;
43	auto push_if = [&](double t){ if (t >= 0) candidates.push_back(Point{p.x + td.x, p.y + td.y}); };
44
45	// Left x = xmin
46	if (fabs(d.x) > EPS) {
47	double t = (xmin - p.x) / d.x; double y = p.y + t*d.y;
48	if (t >= 0 && y >= ymin - 1e-9 && y <= ymax + 1e-9) push_if(t);
49	}
50	// Right x = xmax
51	if (fabs(d.x) > EPS) {
52	double t = (xmax - p.x) / d.x; double y = p.y + t*d.y;
53	if (t >= 0 && y >= ymin - 1e-9 && y <= ymax + 1e-9) push_if(t);
54	}
55	// Bottom y = ymin
56	if (fabs(d.y) > EPS) {
57	double t = (ymin - p.y) / d.y; double x = p.x + t*d.x;
58	if (t >= 0 && x >= xmin - 1e-9 && x <= xmax + 1e-9) push_if(t);
59	}
60	// Top y = ymax
61	if (fabs(d.y) > EPS) {
62	double t = (ymax - p.y) / d.y; double x = p.x + t*d.x;
63	if (t >= 0 && x >= xmin - 1e-9 && x <= xmax + 1e-9) push_if(t);
64	}
65
66	if (candidates.empty()) return p;
67	// Choose the nearest intersection
68	auto d2 = [&](const Point& a, const Point& b){ double dx=a.x-b.x, dy=a.y-b.y; return dxdx+dydy; };
69	int best = 0; double bestd = d2(p, candidates[0]);
70	for (int i = 1; i < (int)candidates.size(); ++i) {
71	double di = d2(p, candidates[i]); if (di < bestd) bestd = di, best = i;
72	}
73	return candidates[best];
74	}
75
76	// Helper to ensure CCW triangle
77	void make_ccw(const vector<Point>& P, Tri& t) {
78	if (orient(P[t.a], P[t.b], P[t.c]) < 0) swap(t.b, t.c);
79	}
80
81	int main(){
82	ios::sync_with_stdio(false);
83	cin.tie(nullptr);
84
85	int n; cin >> n;
86	vector<Point> P(n);
87	for (int i = 0; i < n; ++i) cin >> P[i].x >> P[i].y;
88
89	int m; cin >> m; // number of Delaunay triangles
90	vector<Tri> T(m);
91	for (int i = 0; i < m; ++i) { cin >> T[i].a >> T[i].b >> T[i].c; make_ccw(P, T[i]); }
92
93	// Compute circumcenters
94	vector<Point> CC(m);
95	for (int i = 0; i < m; ++i) {
96	CC[i] = circumcenter(P[T[i].a], P[T[i].b], P[T[i].c]);
97	}
98
99	// Build edge -> incident triangle(s)
100	struct EdgeRec { int tri; int u, v; }; // directed as in its triangle
101	struct Key { int x, y; bool operator==(const Key& o) const { return x==o.x && y==o.y; } };
102	struct KeyHash { size_t operator()(const Key& k) const { return (size_t)k.x*1000003u ^ (size_t)k.y; } };
103	unordered_map<Key, vector<EdgeRec>, KeyHash> adj; adj.reserve(m32);
104
105	auto add_edge = [&](int tri, int u, int v){
106	Key k{min(u,v), max(u,v)}; adj[k].push_back({tri, u, v});
107	};
108
109	for (int i = 0; i < m; ++i) {
110	add_edge(i, T[i].a, T[i].b);
111	add_edge(i, T[i].b, T[i].c);
112	add_edge(i, T[i].c, T[i].a);
113	}
114
115	// Bounding box for clipping hull rays
116	double minx = P[0].x, miny = P[0].y, maxx = P[0].x, maxy = P[0].y;
117	for (int i = 1; i < n; ++i) { minx = min(minx, P[i].x); miny = min(miny, P[i].y); maxx = max(maxx, P[i].x); maxy = max(maxy, P[i].y); }
118	double pad = (max(maxx-minx, maxy-miny) + 1.0) * 0.2; // small padding
119	double X0 = minx - pad, Y0 = miny - pad, X1 = maxx + pad, Y1 = maxy + pad;
120
121	vector<Seg> segs; segs.reserve(m*2);
122	vector<Seg> rays_as_segs; rays_as_segs.reserve(m);
123
124	for (auto &kv : adj) {
125	auto &vec = kv.second;
126	if (vec.size() == 2) {
127	// Interior edge: Voronoi segment between circumcenters
128	int t1 = vec[0].tri, t2 = vec[1].tri;
129	segs.push_back({CC[t1], CC[t2]});
130	} else if (vec.size() == 1) {
131	// Hull edge: Voronoi ray from circumcenter along outward normal of the perpendicular bisector
132	int t = vec[0].tri; int u = vec[0].u, v = vec[0].v;
133	// Directed edge u->v is oriented as in CCW triangle t; inside is to the left.
134	Point a = P[u], b = P[v];
135	Point dir_edge{b.x - a.x, b.y - a.y};
136	// Left normal is (-dy, dx); outward (right) normal is its negative
137	Point right_normal{dir_edge.y, -dir_edge.x}; // rotate -90° equivalently
138	// The Voronoi edge lies along the perpendicular bisector; direction is right_normal
139	Point start = CC[t];
140	Point end = intersectRayAABB(start, right_normal, X0, Y0, X1, Y1);
141	rays_as_segs.push_back({start, end});
142	}
143	}
144
145	// Output: first finite segments, then hull rays as segments to the bbox
146	cout << segs.size() << "\n";
147	cout.setf(std::ios::fixed); cout << setprecision(10);
148	for (auto &s : segs) {
149	cout << s.s.x << " " << s.s.y << " " << s.t.x << " " << s.t.y << "\n";
150	}
151	cout << rays_as_segs.size() << "\n";
152	for (auto &s : rays_as_segs) {
153	cout << s.s.x << " " << s.s.y << " " << s.t.x << " " << s.t.y << "\n";
154	}
155	return 0;
156	}
157