diff --git a/Cargo.lock b/Cargo.lock
index 5e72df25df..adf02bc804 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1141,6 +1141,15 @@ version = "0.3.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "be1e0bca6c3637f992fc1cc7cbc52a78c1ef6db076dbf1059c4323d6a2048376"
 
+[[package]]
+name = "delaunator"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3e1ee323c1275374f7e612d3724d12707079fb6a2117349fe144def656f5a880"
+dependencies = [
+ "robust",
+]
+
 [[package]]
 name = "deranged"
 version = "0.4.0"
@@ -4993,6 +5002,12 @@ dependencies = [
  "serde",
 ]
 
+[[package]]
+name = "robust"
+version = "1.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4e27ee8bb91ca0adcf0ecb116293afa12d393f9c2b9b9cd54d33e8078fe19839"
+
 [[package]]
 name = "ron"
 version = "0.12.0"
@@ -6576,6 +6591,7 @@ name = "vector-nodes"
 version = "0.1.0"
 dependencies = [
  "core-types",
+ "delaunator",
  "dyn-any",
  "futures",
  "glam",
diff --git a/Cargo.toml b/Cargo.toml
index 5685e06aba..c0f932ceb6 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -209,6 +209,7 @@ syn = { version = "2.0", default-features = false, features = [
 	"proc-macro",
 ] }
 lyon_geom = "1.0"
+delaunator = "1.0"
 petgraph = { version = "0.7", default-features = false, features = ["graphmap"] }
 half = { version = "2.4", default-features = false, features = ["bytemuck"] }
 tinyvec = { version = "1", features = ["std"] }
diff --git a/node-graph/libraries/vector-types/src/vector/vector_modification.rs b/node-graph/libraries/vector-types/src/vector/vector_modification.rs
index 2e5f0f4e39..eb67bc1368 100644
--- a/node-graph/libraries/vector-types/src/vector/vector_modification.rs
+++ b/node-graph/libraries/vector-types/src/vector/vector_modification.rs
@@ -629,6 +629,10 @@ where
 	deserializer.deserialize_seq(visitor)
 }
 
+/// Distance below which a path's closing point is treated as coincident with its start point.
+/// Matches Kurbo's default path accuracy, so points within an offset operation's own precision are not split into separate anchors.
+const CLOSE_POINT_TOLERANCE: f64 = 1e-6;
+
 pub struct AppendBezpath<'a> {
 	first_point: Option<Point>,
 	last_point: Option<Point>,
@@ -657,7 +661,11 @@ impl<'a> AppendBezpath<'a> {
 	}
 
 	fn append_segment_and_close_path(&mut self, point: Point, handle: BezierHandles) {
-		let handle = if self.first_point.unwrap() != point {
+		// A path's final point may return to approximately (but not bit-exactly) its start before closing, e.g. a contour
+		// produced by Kurbo's path offsetting. Treat a near-coincident final point as already on the start so we don't
+		// introduce a redundant duplicate anchor and a zero-length closing segment.
+		let endpoints_coincide = (self.first_point.unwrap() - point).hypot2() <= CLOSE_POINT_TOLERANCE * CLOSE_POINT_TOLERANCE;
+		let handle = if !endpoints_coincide {
 			// If the first point is not the same as the last point of the path then we append the segment
 			// with given handle and point and then close the path with linear handle.
 			self.append_segment(point, handle);
diff --git a/node-graph/nodes/vector/Cargo.toml b/node-graph/nodes/vector/Cargo.toml
index 87594c253b..dba094bd4e 100644
--- a/node-graph/nodes/vector/Cargo.toml
+++ b/node-graph/nodes/vector/Cargo.toml
@@ -24,6 +24,7 @@ repeat-nodes = { workspace = true }
 dyn-any = { workspace = true }
 glam = { workspace = true }
 kurbo = { workspace = true }
+delaunator = { workspace = true }
 rand = { workspace = true }
 rustc-hash = { workspace = true }
 log = { workspace = true }
diff --git a/node-graph/nodes/vector/src/lib.rs b/node-graph/nodes/vector/src/lib.rs
index 28ff413ea1..8bdf7005ae 100644
--- a/node-graph/nodes/vector/src/lib.rs
+++ b/node-graph/nodes/vector/src/lib.rs
@@ -2,6 +2,7 @@ pub mod generator_nodes;
 pub mod merge_qr_squares;
 pub mod vector_modification_nodes;
 mod vector_nodes;
+mod voronoi;
 
 #[macro_use]
 extern crate log;
diff --git a/node-graph/nodes/vector/src/vector_nodes.rs b/node-graph/nodes/vector/src/vector_nodes.rs
index c2efc76a16..3d4ee37bc7 100644
--- a/node-graph/nodes/vector/src/vector_nodes.rs
+++ b/node-graph/nodes/vector/src/vector_nodes.rs
@@ -20,9 +20,9 @@ use kurbo::simplify::{SimplifyOptions, simplify_bezpath};
 use kurbo::{Affine, BezPath, DEFAULT_ACCURACY, Line, ParamCurve, ParamCurveArclen, PathEl, PathSeg, Shape};
 use rand::{Rng, SeedableRng};
 use std::collections::hash_map::DefaultHasher;
+use std::collections::{HashMap, HashSet};
 use vector_types::gradient::{build_transform_with_y_preservation, initial_gradient_transform_for_bounding_box};
 use vector_types::subpath::{BezierHandles, ManipulatorGroup};
-use vector_types::vector::PointDomain;
 use vector_types::vector::algorithms::bezpath_algorithms::{self, TValue, eval_pathseg_euclidean, evaluate_bezpath, split_bezpath, tangent_on_bezpath};
 use vector_types::vector::algorithms::merge_by_distance::MergeByDistanceExt;
 use vector_types::vector::algorithms::offset_subpath::offset_bezpath;
@@ -33,6 +33,7 @@ use vector_types::vector::misc::{
 };
 use vector_types::vector::style::{DashPattern, Gradient, PaintOrder, Stroke, StrokeAlign, StrokeCap, StrokeJoin};
 use vector_types::vector::{FillId, PointId, RegionId, SegmentDomain, SegmentId, StrokeId, VectorExt};
+use vector_types::vector::{PointDomain, RegionDomain};
 use vector_types::{GradientSpreadMethod, GradientType};
 
 /// Implemented for `List` types that contain vector items reachable via mutable access.
@@ -1190,6 +1191,168 @@ async fn points_to_polyline(_: impl Ctx, points: Item<Vector>, #[default(true)]
 	points
 }
 
+/// Evens out the distances between points by applying Lloyd's relaxation, moving every interior point toward the center of its Voronoi cell.
+#[node_macro::node(category("Vector: Modifier"), path(core_types::vector))]
+async fn relax_points(
+	_: impl Ctx,
+	/// A vector path or point cloud to relax.
+	source: Item<Vector>,
+	/// The number of relaxation steps to apply. A fractional value runs the whole steps and then blends partway toward one more step, so the amount of relaxation can be animated smoothly.
+	#[default(1.)]
+	#[hard(0..1000)]
+	iterations: Item<f64>,
+) -> Item<Vector> {
+	let mut source = source;
+	let iterations = *iterations.element();
+
+	let vector = source.element_mut();
+	let relaxed = crate::voronoi::relax_sites(vector.point_domain.positions(), iterations);
+	for ((_, position), new_position) in vector.point_domain.positions_mut().zip(relaxed) {
+		*position = new_position;
+	}
+
+	source
+}
+
+/// Builds a Voronoi diagram from the anchor points. Each point claims the region of space closest to it, and those regions tessellate the plane. Cells around the outside are clipped to the convex hull of the points so the diagram stays finite.
+///
+/// When Connect Cells is off, every cell becomes its own closed, fillable subpath. When on, the cells share their common points and segments, forming a single connected mesh with no fillable regions.
+#[node_macro::node(category("Vector"), path(core_types::vector))]
+async fn voronoi_cells(_: impl Ctx, source: Item<Vector>, connect_cells: Item<bool>) -> Item<Vector> {
+	let mut source = source;
+	let connect_cells = *connect_cells.element();
+
+	let vector = source.element_mut();
+	let sites = vector.point_domain.positions().to_vec();
+	let cells = crate::voronoi::voronoi_cells(&sites);
+	if !cells.is_empty() {
+		replace_with_polygons(vector, cells, connect_cells);
+	}
+
+	source
+}
+
+/// Builds a Delaunay triangulation connecting the anchor points. It is the geometric dual of the **Voronoi** node: a mesh of triangles in which no point lies inside any triangle's circumscribed circle.
+///
+/// When Connect Cells is off, every triangle becomes its own closed, fillable subpath. When on, the triangles share their common points and segments, forming a single connected mesh with no fillable regions.
+#[node_macro::node(category("Vector"), path(core_types::vector))]
+async fn triangulate(_: impl Ctx, source: Item<Vector>, connect_cells: Item<bool>) -> Item<Vector> {
+	let mut source = source;
+	let connect_cells = *connect_cells.element();
+
+	let vector = source.element_mut();
+	let sites = vector.point_domain.positions().to_vec();
+	let triangles = crate::voronoi::delaunay_triangles(&sites);
+	if !triangles.is_empty() {
+		// `delaunator` emits triangle vertices clockwise; reverse to `[a, c, b]` so triangles wind counter-clockwise to
+		// match the Voronoi cells and the rest of the framework's fill winding.
+		let polygons = triangles.iter().map(|&[a, b, c]| vec![sites[a], sites[c], sites[b]]).collect();
+		replace_with_polygons(vector, polygons, connect_cells);
+	}
+
+	source
+}
+
+/// Replaces a vector's geometry (points, segments, and regions) with the given closed polygons, preserving its style.
+///
+/// Without `connect_cells`, each polygon becomes its own closed subpath with a fillable region.
+/// With it, coincident vertices are welded and each shared edge is emitted once, producing a connected mesh with no regions.
+pub(crate) fn replace_with_polygons(vector: &mut Vector, polygons: Vec<Vec<DVec2>>, connect_cells: bool) {
+	let mut point_domain = PointDomain::new();
+	let mut segment_domain = SegmentDomain::new();
+	let mut region_domain = RegionDomain::new();
+	let mut next_point = PointId::ZERO;
+	let mut next_segment = SegmentId::ZERO;
+	let mut next_region = RegionId::ZERO;
+
+	if !connect_cells {
+		for polygon in &polygons {
+			if polygon.len() < 3 {
+				continue;
+			}
+
+			let base = point_domain.ids().len();
+			for &position in polygon {
+				point_domain.push(next_point.next_id(), position);
+			}
+
+			let count = polygon.len();
+			let mut first_segment = None;
+			let mut last_segment = None;
+			for i in 0..count {
+				let start = base + i;
+				let end = base + (i + 1) % count;
+				let id = next_segment.next_id();
+				first_segment.get_or_insert(id);
+				last_segment = Some(id);
+				segment_domain.push(id, start, end, BezierHandles::Linear, StrokeId::ZERO);
+			}
+
+			if let (Some(first), Some(last)) = (first_segment, last_segment) {
+				region_domain.push(next_region.next_id(), first..=last, FillId::ZERO);
+			}
+		}
+	} else {
+		// Weld vertices that fall in the same quantization cell so adjacent polygons share points,
+		// and emit each undirected edge only once so adjacent polygons share segments.
+		let tolerance = mesh_weld_tolerance(&polygons);
+		let mut vertex_lookup: HashMap<(i64, i64), usize> = HashMap::new();
+		let mut seen_edges: HashSet<(usize, usize)> = HashSet::new();
+
+		for polygon in &polygons {
+			if polygon.len() < 2 {
+				continue;
+			}
+
+			let indices: Vec<usize> = polygon
+				.iter()
+				.map(|&position| {
+					let key = ((position.x / tolerance).round() as i64, (position.y / tolerance).round() as i64);
+					*vertex_lookup.entry(key).or_insert_with(|| {
+						let index = point_domain.ids().len();
+						point_domain.push(next_point.next_id(), position);
+						index
+					})
+				})
+				.collect();
+
+			let count = indices.len();
+			for i in 0..count {
+				let start = indices[i];
+				let end = indices[(i + 1) % count];
+				if start == end {
+					continue;
+				}
+				let edge = if start < end { (start, end) } else { (end, start) };
+				if seen_edges.insert(edge) {
+					segment_domain.push(next_segment.next_id(), start, end, BezierHandles::Linear, StrokeId::ZERO);
+				}
+			}
+		}
+	}
+
+	vector.point_domain = point_domain;
+	vector.segment_domain = segment_domain;
+	vector.region_domain = region_domain;
+}
+
+/// The distance below which two mesh vertices are welded into one, scaled to the diagram's size so it tracks coordinate magnitude.
+fn mesh_weld_tolerance(polygons: &[Vec<DVec2>]) -> f64 {
+	let mut min = DVec2::splat(f64::MAX);
+	let mut max = DVec2::splat(f64::MIN);
+	for polygon in polygons {
+		for &position in polygon {
+			min = min.min(position);
+			max = max.max(position);
+		}
+	}
+
+	let diagonal = (max - min).length();
+	// Floor the tolerance so an extremely tiny diagram can't underflow `diagonal * 1e-6` to zero, which would divide by
+	// zero when quantizing vertices and weld everything into a single point.
+	if diagonal.is_finite() && diagonal > 0. { (diagonal * 1e-6).max(1e-12) } else { 1e-6 }
+}
+
 #[node_macro::node(category("Vector: Modifier"), path(core_types::vector), properties("offset_path_properties"))]
 async fn offset_path(_: impl Ctx, content: Item<Vector>, distance: Item<f64>, join: Item<StrokeJoin>, #[default(4.)] miter_limit: Item<f64>) -> Item<Vector> {
 	let mut content = content;
@@ -3231,6 +3394,153 @@ mod test {
 		Item::new_from_element(row).with_attribute(ATTR_TRANSFORM, transform)
 	}
 
+	fn item<T>(value: T) -> Item<T> {
+		Item::new_from_element(value)
+	}
+
+	fn vector_item_from_points(points: &[DVec2]) -> Item<Vector> {
+		let mut vector = Vector::default();
+		let mut next_point = PointId::ZERO;
+		for &position in points {
+			vector.point_domain.push(next_point.next_id(), position);
+		}
+		Item::new_from_element(vector)
+	}
+
+	const SQUARE_WITH_CENTER: [DVec2; 5] = [DVec2::new(0., 0.), DVec2::new(10., 0.), DVec2::new(10., 10.), DVec2::new(0., 10.), DVec2::new(5., 5.)];
+
+	#[tokio::test]
+	async fn offset_path_does_not_duplicate_closing_anchors() {
+		// Offsetting closed triangles must not leave each subpath with a redundant start/end anchor (a Kurbo offset
+		// contour returns to approximately, not exactly, its start; that near-coincident point must close, not duplicate).
+		let delaunay = super::triangulate((), vector_item_from_points(&SQUARE_WITH_CENTER), item(false)).await;
+		let offset = super::offset_path((), delaunay, item(0.5), item(StrokeJoin::Miter), item(4.)).await;
+		let result = offset.element();
+
+		let mut subpaths = 0;
+		for (group, closed) in result.stroke_manipulator_groups() {
+			subpaths += 1;
+			assert!(closed, "offset of a closed triangle should stay closed");
+			let first = group.first().unwrap().anchor;
+			let last = group.last().unwrap().anchor;
+			assert!(first.distance(last) > 1e-6, "closed subpath has a duplicated start/end anchor: {first:?} ~= {last:?}");
+		}
+		assert!(subpaths > 0);
+	}
+
+	#[tokio::test]
+	async fn delaunay_disconnected_cells_make_one_region_per_triangle() {
+		let result = super::triangulate((), vector_item_from_points(&SQUARE_WITH_CENTER), item(false)).await;
+		let vector = result.element();
+		// The square plus its center tessellates into four triangles, each its own closed subpath.
+		assert_eq!(vector.region_domain.ids().len(), 4);
+		assert_eq!(vector.segment_domain.ids().len(), 4 * 3);
+		assert_eq!(vector.point_domain.ids().len(), 4 * 3);
+	}
+
+	#[tokio::test]
+	async fn delaunay_and_voronoi_cells_share_winding() {
+		fn signed_area(anchors: &[DVec2]) -> f64 {
+			(0..anchors.len()).map(|i| anchors[i].perp_dot(anchors[(i + 1) % anchors.len()])).sum::<f64>() / 2.
+		}
+		fn subpath_winding_signs(vector: &Vector) -> Vec<f64> {
+			vector
+				.stroke_manipulator_groups()
+				.map(|(group, _)| signed_area(&group.iter().map(|g| g.anchor).collect::<Vec<_>>()).signum())
+				.collect()
+		}
+
+		// The Rectangle and Ellipse generators define the framework's fill winding convention; each is built from these
+		// subpath constructors (`Subpath::new_rectangle` / `Subpath::new_ellipse`), so their winding is the source of truth.
+		use vector_types::subpath::Subpath;
+		let rectangle = Vector::from_subpath(Subpath::new_rectangle(DVec2::new(-50., -50.), DVec2::new(50., 50.)));
+		let ellipse = Vector::from_subpath(Subpath::new_ellipse(DVec2::new(-50., -25.), DVec2::new(50., 25.)));
+		let expected = subpath_winding_signs(&rectangle)[0];
+		assert_eq!(subpath_winding_signs(&ellipse)[0], expected, "Rectangle and Ellipse should agree on winding");
+
+		// Delaunay and Voronoi must emit subpaths that wind the same way as those generators.
+		let delaunay = super::triangulate((), vector_item_from_points(&SQUARE_WITH_CENTER), item(false)).await;
+		let voronoi = super::voronoi_cells((), vector_item_from_points(&SQUARE_WITH_CENTER), item(false)).await;
+		for sign in subpath_winding_signs(delaunay.element()) {
+			assert_eq!(sign, expected, "Delaunay subpath winding should match the Rectangle/Ellipse generators");
+		}
+		for sign in subpath_winding_signs(voronoi.element()) {
+			assert_eq!(sign, expected, "Voronoi subpath winding should match the Rectangle/Ellipse generators");
+		}
+	}
+
+	#[tokio::test]
+	async fn delaunay_shared_mesh_welds_points_and_shares_edges() {
+		let result = super::triangulate((), vector_item_from_points(&SQUARE_WITH_CENTER), item(true)).await;
+		let vector = result.element();
+		// The connected mesh reuses the five input points and shares edges, with no fillable regions.
+		assert_eq!(vector.region_domain.ids().len(), 0);
+		assert_eq!(vector.point_domain.ids().len(), 5);
+		// Four hull edges plus four spokes to the center, each emitted once.
+		assert_eq!(vector.segment_domain.ids().len(), 8);
+	}
+
+	#[tokio::test]
+	async fn voronoi_disconnected_cells_make_a_region_per_cell() {
+		let result = super::voronoi_cells((), vector_item_from_points(&SQUARE_WITH_CENTER), item(false)).await;
+		let vector = result.element();
+		let regions = vector.region_domain.ids().len();
+		assert!(regions > 0, "expected at least one Voronoi region");
+		// Every region is a closed subpath, so segments and points come in matched per-region loops.
+		assert_eq!(vector.segment_domain.ids().len(), vector.point_domain.ids().len());
+
+		// Clipping to the convex hull keeps all cell vertices within the input bounds.
+		for &position in vector.point_domain.positions() {
+			assert!(position.x >= -1e-6 && position.x <= 10. + 1e-6);
+			assert!(position.y >= -1e-6 && position.y <= 10. + 1e-6);
+		}
+	}
+
+	#[tokio::test]
+	async fn voronoi_shared_mesh_has_no_regions() {
+		let result = super::voronoi_cells((), vector_item_from_points(&SQUARE_WITH_CENTER), item(true)).await;
+		let vector = result.element();
+		assert_eq!(vector.region_domain.ids().len(), 0);
+		assert!(vector.segment_domain.ids().len() > 0);
+	}
+
+	#[tokio::test]
+	async fn voronoi_leaves_degenerate_input_untouched() {
+		// Two points cannot form a diagram, so the element passes through unchanged.
+		let points = [DVec2::new(0., 0.), DVec2::new(1., 1.)];
+		let result = super::voronoi_cells((), vector_item_from_points(&points), item(false)).await;
+		let vector = result.element();
+		assert_eq!(vector.point_domain.ids().len(), 2);
+		assert_eq!(vector.segment_domain.ids().len(), 0);
+	}
+
+	#[tokio::test]
+	async fn relax_points_redistributes_anchors() {
+		// Four hull corners plus two off-center interior points.
+		let points = [
+			DVec2::new(0., 0.),
+			DVec2::new(10., 0.),
+			DVec2::new(10., 10.),
+			DVec2::new(0., 10.),
+			DVec2::new(3., 4.),
+			DVec2::new(7., 5.),
+		];
+		let result = super::relax_points((), vector_item_from_points(&points), item(2.)).await;
+		let vector = result.element();
+
+		// Relaxation preserves the point count but repositions the interior anchors within the hull.
+		assert_eq!(vector.point_domain.ids().len(), points.len());
+		assert_ne!(vector.point_domain.positions(), &points[..]);
+		// The convex-hull corners are pinned.
+		for i in 0..4 {
+			assert_eq!(vector.point_domain.positions()[i], points[i], "hull corner {i} should be pinned");
+		}
+		for &point in vector.point_domain.positions() {
+			assert!(point.x >= -1e-6 && point.x <= 10. + 1e-6);
+			assert!(point.y >= -1e-6 && point.y <= 10. + 1e-6);
+		}
+	}
+
 	#[tokio::test]
 	async fn bounding_box() {
 		let bounding_box = super::bounding_box((), Item::new_from_element(Vector::from_bezpath(Rect::new(-1., -1., 1., 1.).to_path(DEFAULT_ACCURACY)))).await;
diff --git a/node-graph/nodes/vector/src/voronoi.rs b/node-graph/nodes/vector/src/voronoi.rs
new file mode 100644
index 0000000000..05f6eca931
--- /dev/null
+++ b/node-graph/nodes/vector/src/voronoi.rs
@@ -0,0 +1,472 @@
+//! Geometry for the Voronoi and Delaunay nodes.
+//!
+//! Both diagrams are derived from a single Delaunay triangulation (computed by the `delaunator` crate). The Voronoi
+//! diagram is the geometric dual of that triangulation: each Voronoi vertex is the circumcenter of a Delaunay triangle,
+//! and each Voronoi edge connects the circumcenters of two triangles that share a Delaunay edge.
+//!
+//! Each function here reduces a diagram to a set of closed polygons (one per Delaunay triangle or per Voronoi cell). The
+//! nodes then assemble those polygons into vector geometry, either as separate filled subpaths or as a shared mesh of
+//! welded points and segments. Voronoi cells around the convex hull are unbounded, so they are clipped to the convex hull
+//! of the input sites, which also bounds the whole diagram to a finite region.
+
+use delaunator::{EMPTY, Point, triangulate};
+use glam::DVec2;
+
+/// Computes the Delaunay triangulation of `sites`, returning each triangle as a triple of indices into `sites`.
+///
+/// Returns an empty vector when there are fewer than three points or they are all colinear (no triangle exists).
+pub fn delaunay_triangles(sites: &[DVec2]) -> Vec<[usize; 3]> {
+	let points: Vec<Point> = sites.iter().map(|p| Point { x: p.x, y: p.y }).collect();
+	let triangulation = triangulate(&points);
+	triangulation.triangles.chunks_exact(3).map(|t| [t[0], t[1], t[2]]).collect()
+}
+
+/// Computes the Voronoi cell of every site, each clipped to the convex hull of `sites`.
+///
+/// Returns one closed polygon per site that produces a non-empty cell (degenerate or fully-clipped cells are omitted,
+/// so the result may be shorter than `sites`). Returns an empty vector when no triangulation exists (fewer than three points or all colinear).
+pub fn voronoi_cells(sites: &[DVec2]) -> Vec<Vec<DVec2>> {
+	voronoi_cells_per_site(sites).0.into_iter().flatten().collect()
+}
+
+/// Applies Lloyd's relaxation: each step moves every interior site to the centroid of its Voronoi cell, yielding a more
+/// even (centroidal) point distribution. A fractional `iterations` runs the whole-number steps and then blends each site
+/// partway toward the result of one more step, so the relaxation can be animated smoothly. Returns the sites unchanged when
+/// `iterations` is 0 or no diagram can be formed.
+///
+/// The convex-hull (perimeter) sites are pinned so the point cloud's outline is preserved. Otherwise, clipping the
+/// unbounded perimeter cells would drag those sites around (inward for the convex hull, or outward into the corners of a
+/// fixed bounding box), distorting the shape over successive iterations.
+pub fn relax_sites(sites: &[DVec2], iterations: f64) -> Vec<DVec2> {
+	const MAX_STEPS_FOR_SAFETY: f64 = 1000.;
+	let iterations = iterations.clamp(0., MAX_STEPS_FOR_SAFETY);
+	let whole_steps = iterations.floor();
+	let fraction = iterations - whole_steps;
+
+	let mut current = sites.to_vec();
+	for _ in 0..whole_steps as u32 {
+		current = relax_once(&current);
+	}
+
+	// Blend each site partway toward one further step for the fractional remainder.
+	if fraction > 0. {
+		let next = relax_once(&current);
+		for (point, target) in current.iter_mut().zip(next) {
+			*point = point.lerp(target, fraction);
+		}
+	}
+
+	current
+}
+
+/// Performs a single Lloyd relaxation step: moves every interior site to its Voronoi cell centroid,
+/// leaving the pinned convex-hull (perimeter) sites in place.
+fn relax_once(sites: &[DVec2]) -> Vec<DVec2> {
+	let (cells, is_hull) = voronoi_cells_per_site(sites);
+	let mut relaxed = sites.to_vec();
+	for ((site, cell), on_hull) in relaxed.iter_mut().zip(cells).zip(is_hull) {
+		if on_hull {
+			continue;
+		}
+		if let Some(centroid) = cell.as_deref().and_then(polygon_centroid) {
+			*site = centroid;
+		}
+	}
+	relaxed
+}
+
+/// The area-weighted centroid of a simple polygon, or `None` if it has fewer than three vertices or zero area.
+fn polygon_centroid(polygon: &[DVec2]) -> Option<DVec2> {
+	if polygon.len() < 3 {
+		return None;
+	}
+	let mut double_area = 0.;
+	let mut weighted = DVec2::ZERO;
+	for i in 0..polygon.len() {
+		let a = polygon[i];
+		let b = polygon[(i + 1) % polygon.len()];
+		let cross = a.perp_dot(b);
+		double_area += cross;
+		weighted += (a + b) * cross;
+	}
+	(double_area.abs() >= f64::EPSILON).then(|| weighted / (3. * double_area))
+}
+
+/// Computes each site's clipped Voronoi cell, aligned with `sites` (index `i` is the cell of `sites[i]`), together with a
+/// per-site flag marking the convex-hull (perimeter) sites. A cell is `None` when the site has no incident triangle
+/// (e.g. a coincident duplicate) or its cell vanishes after clipping.
+fn voronoi_cells_per_site(sites: &[DVec2]) -> (Vec<Option<Vec<DVec2>>>, Vec<bool>) {
+	let points: Vec<Point> = sites.iter().map(|p| Point { x: p.x, y: p.y }).collect();
+	let triangulation = triangulate(&points);
+	if triangulation.triangles.is_empty() {
+		return (vec![None; sites.len()], vec![false; sites.len()]);
+	}
+
+	let triangles = &triangulation.triangles;
+	let halfedges = &triangulation.halfedges;
+	let hull_indices = &triangulation.hull;
+
+	// Mark which sites lie on the convex hull (the diagram's perimeter).
+	let mut is_hull = vec![false; sites.len()];
+	for &index in hull_indices {
+		is_hull[index] = true;
+	}
+
+	// One Voronoi vertex per Delaunay triangle.
+	let circumcenters: Vec<DVec2> = triangles.chunks_exact(3).map(|t| circumcenter(sites[t[0]], sites[t[1]], sites[t[2]])).collect();
+
+	// The convex hull polygon, which clips the diagram to a finite region.
+	let hull: Vec<DVec2> = hull_indices.iter().map(|&i| sites[i]).collect();
+
+	// `inedges[p]` is a half-edge ending at site `p`, preferring a hull half-edge so a hull cell's walk starts on the boundary.
+	let mut inedges = vec![EMPTY; sites.len()];
+	for edge in 0..triangles.len() {
+		let endpoint = triangles[next_halfedge(edge)];
+		if halfedges[edge] == EMPTY || inedges[endpoint] == EMPTY {
+			inedges[endpoint] = edge;
+		}
+	}
+
+	// Outward ray directions for the two hull edges meeting at each hull site, used to project its unbounded cell outward.
+	// Both are zero for interior sites.
+	let mut ray_in = vec![DVec2::ZERO; sites.len()];
+	let mut ray_out = vec![DVec2::ZERO; sites.len()];
+	if let Some(&last) = hull_indices.last() {
+		let mut previous = last;
+		for &current in hull_indices {
+			let p0 = sites[previous];
+			let p1 = sites[current];
+			// Perpendicular to the hull edge `previous -> current`, pointing away from the hull interior.
+			let perpendicular = DVec2::new(p0.y - p1.y, p1.x - p0.x);
+			ray_out[previous] = perpendicular;
+			ray_in[current] = perpendicular;
+			previous = current;
+		}
+	}
+
+	// Length to extend unbounded cell rays so they reach past the hull before clipping trims them back to it.
+	let far = bounding_diagonal(&hull) * 10. + 1.;
+
+	let cells = (0..sites.len())
+		.map(|site| {
+			let mut polygon = cell_polygon(site, halfedges, &circumcenters, &inedges)?;
+
+			// A hull site's cell is unbounded; cap its open ends with far points along the outward hull-edge normals so the
+			// convex-hull clip below closes it off at the boundary.
+			let unbounded = ray_in[site] != DVec2::ZERO || ray_out[site] != DVec2::ZERO;
+			if unbounded {
+				if let Some(&first) = polygon.first() {
+					polygon.insert(0, first + ray_in[site].normalize_or_zero() * far);
+				}
+				if let Some(&last) = polygon.last() {
+					polygon.push(last + ray_out[site].normalize_or_zero() * far);
+				}
+			}
+
+			let clipped = clip_to_convex(&polygon, &hull);
+			(clipped.len() >= 3).then_some(clipped)
+		})
+		.collect();
+
+	(cells, is_hull)
+}
+
+/// The circumcenter of a triangle, computed relative to `a` for numerical stability. Falls back to the centroid for a
+/// degenerate (colinear) triangle.
+fn circumcenter(a: DVec2, b: DVec2, c: DVec2) -> DVec2 {
+	let d = b - a;
+	let e = c - a;
+	let determinant = d.x * e.y - d.y * e.x;
+	if determinant.abs() < f64::EPSILON {
+		return (a + b + c) / 3.;
+	}
+	let factor = 0.5 / determinant;
+	let bl = d.length_squared();
+	let cl = e.length_squared();
+	DVec2::new(a.x + (e.y * bl - d.y * cl) * factor, a.y + (d.x * cl - e.x * bl) * factor)
+}
+
+/// Walks the Delaunay triangles incident to `site` and collects their circumcenters in order, forming the site's
+/// Voronoi cell polygon. The polygon is closed for interior sites and open (a fan ending at the hull) for hull sites.
+/// Returns `None` for a site with no incident triangle (e.g. a coincident duplicate point).
+fn cell_polygon(site: usize, halfedges: &[usize], circumcenters: &[DVec2], inedges: &[usize]) -> Option<Vec<DVec2>> {
+	let start = inedges[site];
+	if start == EMPTY {
+		return None;
+	}
+
+	let mut polygon = Vec::new();
+	let mut edge = start;
+	loop {
+		polygon.push(circumcenters[edge / 3]);
+		edge = halfedges[next_halfedge(edge)];
+		if edge == EMPTY || edge == start {
+			break;
+		}
+	}
+
+	Some(polygon)
+}
+
+/// The next half-edge within the same triangle (triangles store three consecutive half-edges).
+fn next_halfedge(edge: usize) -> usize {
+	if edge % 3 == 2 { edge - 2 } else { edge + 1 }
+}
+
+/// Clips `subject` to the convex polygon `clip` using the Sutherland–Hodgman algorithm. The clip polygon may wind either way.
+/// (The subject doesn't need to be convex.) Returns the clipped polygon (empty if it lies entirely outside the clip region).
+fn clip_to_convex(subject: &[DVec2], clip: &[DVec2]) -> Vec<DVec2> {
+	if clip.len() < 3 {
+		return Vec::new();
+	}
+
+	// Normalize the clip polygon to counter-clockwise so "inside" is consistently to the left of each directed edge.
+	let mut clip = clip.to_vec();
+	if signed_area(&clip) < 0. {
+		clip.reverse();
+	}
+
+	let mut output = subject.to_vec();
+	for i in 0..clip.len() {
+		if output.is_empty() {
+			break;
+		}
+
+		let edge_start = clip[i];
+		let edge_end = clip[(i + 1) % clip.len()];
+		let edge = edge_end - edge_start;
+		let inside = |p: DVec2| edge.x * (p.y - edge_start.y) - edge.y * (p.x - edge_start.x) >= 0.;
+
+		let input = std::mem::take(&mut output);
+		for j in 0..input.len() {
+			let current = input[j];
+			let previous = input[(j + input.len() - 1) % input.len()];
+			let current_inside = inside(current);
+			let previous_inside = inside(previous);
+
+			if current_inside {
+				if !previous_inside && let Some(crossing) = line_intersection(previous, current, edge_start, edge_end) {
+					output.push(crossing);
+				}
+				output.push(current);
+			} else if previous_inside && let Some(crossing) = line_intersection(previous, current, edge_start, edge_end) {
+				output.push(crossing);
+			}
+		}
+	}
+
+	output
+}
+
+/// The signed area of a polygon (positive for counter-clockwise winding).
+fn signed_area(polygon: &[DVec2]) -> f64 {
+	let mut area = 0.;
+	for i in 0..polygon.len() {
+		let a = polygon[i];
+		let b = polygon[(i + 1) % polygon.len()];
+		area += a.x * b.y - b.x * a.y;
+	}
+	area / 2.
+}
+
+/// The intersection point of the segment `p1 -> p2` with the infinite line through `a` and `b`, or `None` if parallel.
+fn line_intersection(p1: DVec2, p2: DVec2, a: DVec2, b: DVec2) -> Option<DVec2> {
+	let r = p2 - p1;
+	let s = b - a;
+	let denominator = r.x * s.y - r.y * s.x;
+	if denominator.abs() < f64::EPSILON {
+		return None;
+	}
+	let t = ((a.x - p1.x) * s.y - (a.y - p1.y) * s.x) / denominator;
+	Some(p1 + r * t)
+}
+
+/// The diagonal length of the axis-aligned bounding box of `points`.
+fn bounding_diagonal(points: &[DVec2]) -> f64 {
+	let mut min = DVec2::splat(f64::MAX);
+	let mut max = DVec2::splat(f64::MIN);
+	for &p in points {
+		min = min.min(p);
+		max = max.max(p);
+	}
+	let diagonal = (max - min).length();
+	if diagonal.is_finite() && diagonal > 0. { diagonal } else { 0. }
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+
+	fn square_with_center() -> Vec<DVec2> {
+		vec![DVec2::new(0., 0.), DVec2::new(10., 0.), DVec2::new(10., 10.), DVec2::new(0., 10.), DVec2::new(5., 5.)]
+	}
+
+	#[test]
+	fn delaunay_triangles_wind_counter_clockwise() {
+		// `delaunator` returns clockwise triangles, so `delaunay_triangles` keeps that order, but `voronoi_cells` are
+		// counter-clockwise. This documents the raw orientation; the Delaunay node reverses it to match the cells.
+		let sites = square_with_center();
+		for t in delaunay_triangles(&sites) {
+			let poly = [sites[t[0]], sites[t[1]], sites[t[2]]];
+			assert!(signed_area(&poly) < 0., "delaunator triangles are expected to be clockwise");
+		}
+		for cell in voronoi_cells(&sites) {
+			assert!(signed_area(&cell) > 0., "voronoi cells are expected to be counter-clockwise");
+		}
+	}
+
+	#[test]
+	fn delaunay_triangulates_square() {
+		let triangles = delaunay_triangles(&square_with_center());
+		// Four corner-to-center triangles tessellate the square.
+		assert_eq!(triangles.len(), 4);
+		for triangle in triangles {
+			for index in triangle {
+				assert!(index < 5);
+			}
+		}
+	}
+
+	#[test]
+	fn delaunay_degenerate_inputs_produce_no_triangles() {
+		assert!(delaunay_triangles(&[]).is_empty());
+		assert!(delaunay_triangles(&[DVec2::new(1., 1.)]).is_empty());
+		assert!(delaunay_triangles(&[DVec2::new(0., 0.), DVec2::new(1., 1.)]).is_empty());
+		// Colinear points have no triangulation.
+		let colinear = vec![DVec2::new(0., 0.), DVec2::new(1., 1.), DVec2::new(2., 2.)];
+		assert!(delaunay_triangles(&colinear).is_empty());
+	}
+
+	#[test]
+	fn voronoi_cells_tile_the_hull() {
+		// The clipped cells partition the convex hull, so their (counter-clockwise, positive) areas sum to the hull's area
+		// (100 for the 10x10 square). If the outward projection direction were inverted, the boundary cells would collapse
+		// inward and the total would fall well short of 100.
+		let sites = square_with_center();
+		let total: f64 = voronoi_cells(&sites).iter().map(|cell| signed_area(cell)).sum();
+		assert!((total - 100.).abs() < 1e-6, "cells should tile the hull (area 100), got {total}");
+	}
+
+	#[test]
+	fn voronoi_cells_stay_within_the_hull() {
+		let sites = square_with_center();
+		let cells = voronoi_cells(&sites);
+		assert!(!cells.is_empty());
+		// Clipping to the hull keeps every vertex inside the input bounds (with a small tolerance for float error).
+		for cell in &cells {
+			assert!(cell.len() >= 3);
+			for &vertex in cell {
+				assert!(vertex.x >= -1e-6 && vertex.x <= 10. + 1e-6, "x out of bounds: {}", vertex.x);
+				assert!(vertex.y >= -1e-6 && vertex.y <= 10. + 1e-6, "y out of bounds: {}", vertex.y);
+			}
+		}
+	}
+
+	#[test]
+	fn relaxation_with_zero_iterations_is_identity() {
+		let sites = square_with_center();
+		assert_eq!(relax_sites(&sites, 0.), sites);
+	}
+
+	#[test]
+	fn relaxation_moves_points_and_keeps_them_in_the_hull() {
+		// Add a point clustered near the center; relaxation should redistribute the points without leaving the hull.
+		let mut sites = square_with_center();
+		sites.push(DVec2::new(5.5, 4.5));
+		let relaxed = relax_sites(&sites, 3.);
+
+		assert_eq!(relaxed.len(), sites.len());
+		assert_ne!(relaxed, sites, "relaxation should move the points");
+		for &point in &relaxed {
+			assert!(point.x >= -1e-6 && point.x <= 10. + 1e-6, "x out of bounds: {}", point.x);
+			assert!(point.y >= -1e-6 && point.y <= 10. + 1e-6, "y out of bounds: {}", point.y);
+		}
+	}
+
+	#[test]
+	fn relaxation_pins_the_convex_hull() {
+		// The four corners form the convex hull and must stay fixed; the interior points must relax.
+		let sites = vec![
+			DVec2::new(0., 0.),
+			DVec2::new(10., 0.),
+			DVec2::new(10., 10.),
+			DVec2::new(0., 10.),
+			DVec2::new(3., 3.),
+			DVec2::new(7., 4.),
+		];
+		let relaxed = relax_sites(&sites, 4.);
+
+		for i in 0..4 {
+			assert_eq!(relaxed[i], sites[i], "convex hull point {i} should be pinned");
+		}
+		assert!(relaxed[4] != sites[4] || relaxed[5] != sites[5], "interior points should relax");
+	}
+
+	#[test]
+	fn relaxation_interpolates_fractional_iterations() {
+		let sites = vec![
+			DVec2::new(0., 0.),
+			DVec2::new(10., 0.),
+			DVec2::new(10., 10.),
+			DVec2::new(0., 10.),
+			DVec2::new(3., 4.),
+			DVec2::new(7., 6.),
+		];
+
+		// A fractional count lands exactly midway between the two bracketing whole-step results, exercising both the
+		// pure-fraction path (0.5) and the whole-steps-then-fraction path (2.5).
+		for whole in [0., 2.] {
+			let lower = relax_sites(&sites, whole);
+			let upper = relax_sites(&sites, whole + 1.);
+			let half = relax_sites(&sites, whole + 0.5);
+			for i in 0..sites.len() {
+				let expected = (lower[i] + upper[i]) / 2.;
+				assert!((half[i] - expected).length() < 1e-9, "index {i} at {whole}.5: {half:?} vs {expected:?}", half = half[i]);
+			}
+		}
+	}
+
+	#[test]
+	fn relaxation_leaves_degenerate_input_unchanged() {
+		// Fewer than three points cannot form a diagram, so relaxation is a no-op.
+		let sites = vec![DVec2::new(0., 0.), DVec2::new(1., 1.)];
+		assert_eq!(relax_sites(&sites, 5.), sites);
+	}
+
+	#[test]
+	fn relaxation_clamps_extreme_iteration_counts() {
+		let sites = vec![
+			DVec2::new(0., 0.),
+			DVec2::new(10., 0.),
+			DVec2::new(10., 10.),
+			DVec2::new(0., 10.),
+			DVec2::new(3., 4.),
+			DVec2::new(7., 6.),
+		];
+		// A huge or infinite count must clamp to the converged result rather than hang on a billions-long loop.
+		let converged = relax_sites(&sites, 1000.);
+		assert_eq!(relax_sites(&sites, 1e9), converged);
+		assert_eq!(relax_sites(&sites, f64::INFINITY), converged);
+		// NaN and negative counts resolve to zero steps, leaving the sites unchanged.
+		assert_eq!(relax_sites(&sites, f64::NAN), sites);
+		assert_eq!(relax_sites(&sites, -5.), sites);
+	}
+
+	#[test]
+	fn voronoi_center_cell_is_bounded() {
+		// A ring of points around a center yields a finite cell for the center site.
+		let mut sites = vec![DVec2::new(0., 0.)];
+		for i in 0..6 {
+			let angle = i as f64 / 6. * std::f64::consts::TAU;
+			sites.push(DVec2::new(angle.cos() * 10., angle.sin() * 10.));
+		}
+		let cells = voronoi_cells(&sites);
+		assert!(!cells.is_empty());
+		// Every cell is a finite polygon with no runaway coordinates.
+		for cell in &cells {
+			for &vertex in cell {
+				assert!(vertex.length() < 100., "unbounded cell vertex: {vertex:?}");
+			}
+		}
+	}
+}