maximumInscribedPack(PShape shape, double minRadius,
return out;
}
+ /**
+ * Fills the gaps in an existing circle packing with new circles, using the
+ * Largest Empty Circle (LEC) algorithm seeded with the existing packing.
+ *
+ * The existing circles are treated as constraints: new circles will not overlap
+ * them (nor each other, nor the shape boundary). Circles are found
+ * largest-first until the next circle would be smaller than {@code minRadius}.
+ *
+ * @param shape The shape within which circles will be packed.
+ * @param existingCircles An existing circle packing to fill in (perhaps
+ * produced by another method in this class). Each
+ * PVector represents one circle: (.x, .y) is the center
+ * and .z is the radius.
+ * @param minRadius The minimum allowed radius for the new circles.
+ * Filling stops once the largest remaining gap is
+ * smaller than this.
+ * @param tolerance The tolerance value to control the LEC algorithm's
+ * accuracy. Higher values yield faster results but lower
+ * accuracy. A value of 1 is a good starting point.
+ * @return A list of the new circles only (the input circles are not
+ * echoed back), each as a PVector: (.x, .y) represent the center point
+ * and .z represents the radius.
+ * @since 2.3
+ * @see #maximumInscribedPack(PShape, double, double)
+ */
+ public static List fillPack(PShape shape, Collection existingCircles, double minRadius, double tolerance) {
+ tolerance = Math.max(0.01, tolerance);
+ minRadius = Math.max(0.01, minRadius);
+
+ final List seeds = existingCircles.stream().map(c -> new Coordinate(c.x, c.y, c.z)).toList();
+ final LargestEmptyCircles lec = new LargestEmptyCircles(fromPShape(shape), null, seeds, tolerance);
+
+ final List out = new ArrayList<>();
+ double[] c;
+ while ((c = lec.findNextLEC())[2] >= minRadius) {
+ out.add(new PVector((float) c[0], (float) c[1], (float) c[2]));
+ }
+ return out;
+ }
+
+ /**
+ * Fills the gaps in an existing circle packing with exactly {@code n} new
+ * circles, using the Largest Empty Circle (LEC) algorithm seeded with the
+ * existing packing.
+ *
+ * The existing circles are treated as constraints: new circles will not overlap
+ * them (nor each other, nor the shape boundary).
+ *
+ * @param shape The shape within which circles will be packed.
+ * @param existingCircles An existing circle packing to fill in (perhaps
+ * produced by another method in this class). Each
+ * PVector represents one circle: (.x, .y) is the center
+ * and .z is the radius.
+ * @param n The number of new circles to find and pack.
+ * @param tolerance The tolerance value to control the LEC algorithm's
+ * accuracy. Higher values yield faster results but lower
+ * accuracy. A value of 1 is a good starting point.
+ * @return A list of the new circles only (the input circles are not
+ * echoed back), each as a PVector: (.x, .y) represent the center point
+ * and .z represents the radius.
+ * @since 2.3
+ * @see #fillPack(PShape, Collection, double, double)
+ */
+ public static List fillPack(PShape shape, Collection existingCircles, int n, double tolerance) {
+ tolerance = Math.max(0.01, tolerance);
+
+ final List seeds = existingCircles.stream().map(c -> new Coordinate(c.x, c.y, c.z)).toList();
+ final LargestEmptyCircles lec = new LargestEmptyCircles(fromPShape(shape), null, seeds, tolerance);
+
+ final List out = new ArrayList<>(n);
+ for (int i = 0; i < n; i++) {
+ double[] c = lec.findNextLEC();
+ out.add(new PVector((float) c[0], (float) c[1], (float) c[2]));
+ }
+ return out;
+ }
+
/**
* Generates a circle packing having a pattern of tangencies specified by a
* triangulation.
diff --git a/src/main/java/micycle/pgs/PGS_Conversion.java b/src/main/java/micycle/pgs/PGS_Conversion.java
index b08f82dc..5fcf456e 100644
--- a/src/main/java/micycle/pgs/PGS_Conversion.java
+++ b/src/main/java/micycle/pgs/PGS_Conversion.java
@@ -2121,7 +2121,7 @@ public static PShape fromQuadraticBezier(PVector start, PVector controlPoint, PV
*/
public static PShape fromCubicBezier(PVector start, PVector controlPoint1, PVector controlPoint2, PVector end) {
CubicBezier bezier = new CubicBezier(start.x, start.y, controlPoint1.x, controlPoint1.y, controlPoint2.x, controlPoint2.y, end.x, end.y);
- double[][] samples = bezier.sampleEquidistantPoints(BEZIER_SAMPLE_DISTANCE);
+ double[][] samples = bezier.sampleAdaptive(1.0);
final List coords = new ArrayList<>(samples.length);
for (double[] sample : samples) {
coords.add(new PVector((float) sample[0], (float) sample[1]));
diff --git a/src/main/java/micycle/pgs/PGS_ShapeBoolean.java b/src/main/java/micycle/pgs/PGS_ShapeBoolean.java
index 926dba63..3454f3be 100644
--- a/src/main/java/micycle/pgs/PGS_ShapeBoolean.java
+++ b/src/main/java/micycle/pgs/PGS_ShapeBoolean.java
@@ -27,7 +27,7 @@
import org.locationtech.jts.operation.union.UnaryUnionOp;
import org.locationtech.jts.util.GeometricShapeFactory;
-import com.github.micycle1.geoblitz.DiskUnion;
+import com.github.micycle1.geoblitz.CircleUnion;
import micycle.pgs.commons.FastOverlapRegions;
import micycle.pgs.commons.Nullable;
@@ -231,7 +231,7 @@ public static PShape union(PShape... shapes) {
*/
public static PShape unionCircles(Collection circles) {
var disks = circles.stream().map(c -> PGS.coordFromPVector(c)).toList();
- var union = DiskUnion.union(disks, PGS_Conversion.BEZIER_SAMPLE_DISTANCE);
+ var union = CircleUnion.union(disks, PGS_Conversion.BEZIER_SAMPLE_DISTANCE);
return toPShape(union);
}
diff --git a/src/main/java/micycle/pgs/PGS_Voronoi.java b/src/main/java/micycle/pgs/PGS_Voronoi.java
index 80f66953..d4c29bd0 100644
--- a/src/main/java/micycle/pgs/PGS_Voronoi.java
+++ b/src/main/java/micycle/pgs/PGS_Voronoi.java
@@ -483,43 +483,6 @@ public static PShape compoundVoronoi(PShape shape, double[] bounds) {
return voronoiCells;
}
- /**
- * Generates an additively weighted Voronoi diagram (AWVD) for a set of
- * weighted point sites, clipped to the provided bounding rectangle.
- *
- * AWVDs are a generalisation of standard Voronoi diagrams where each site has
- * an additive weight. Distances are compared using an adjusted metric of the
- * form:
- *
- *
- * d(p, s) = ||p - s|| - w
- *
- *
- * where {@code s} is the site location and {@code w} is its weight. Increasing
- * a site's weight tends to expand its cell; decreasing it tends to shrink the
- * cell. Unlike standard Voronoi diagrams, AWVD cell boundaries are generally
- * curved (hyperbolic arcs), and some sites may end up with empty cells
- * depending on weights and configuration.
- *
- * Each input {@link PVector} encodes one site where:
- *
- * - {@code (.x, .y)} is the site coordinate
- * - {@code .z} is the site's weight (in the same units as {@code x/y})
- *
- *
- * @param weightedSites a collection of weighted sites encoded as PVectors:
- * {@code (.x, .y)} position and {@code .z} weight
- * @param bounds an array of the form {@code [minX, minY, maxX, maxY]}
- * defining the clipping bounds of the diagram; must fully
- * contain all sites
- * @return a GROUP {@link PShape} where each child shape is a (possibly curved)
- * AWVD cell polygon clipped to {@code bounds}
- * @since 2.2
- */
- public static PShape additivelyWeightedVoronoi(Collection weightedSites, double[] bounds) {
- return additivelyWeightedVoronoi(weightedSites, bounds, false);
- }
-
/**
* Generates an additively weighted Voronoi diagram (AWVD) for a set of
* weighted point sites, clipped to the provided bounding rectangle.
@@ -551,38 +514,28 @@ public static PShape additivelyWeightedVoronoi(Collection weightedSites
* simplify the interior boundaries.
*
*
- * @param weightedSites a collection of weighted sites encoded as PVectors:
- * {@code (.x, .y)} position and {@code .z} weight
- * @param bounds an array of the form {@code [minX, minY, maxX, maxY]}
- * defining the clipping bounds of the diagram; must
- * fully contain all sites
- * @param forceConforming whether to apply additional processing to try to
- * ensure adjacent cells form a conforming coverage
- * (i.e., no tiny gaps between cells)
+ * @param weightedSites a collection of weighted sites encoded as PVectors:
+ * {@code (.x, .y)} position and {@code .z} weight
+ * @param bounds an array of the form {@code [minX, minY, maxX, maxY]}
+ * defining the clipping bounds of the diagram; must fully
+ * contain all sites
* @return a GROUP {@link PShape} where each child shape is a (possibly curved)
* AWVD cell polygon clipped to {@code bounds}
* @since 2.2
*/
- public static PShape additivelyWeightedVoronoi(Collection weightedSites, double[] bounds, boolean forceConforming) {
+ public static PShape additivelyWeightedVoronoi(Collection weightedSites, double[] bounds) {
var sites = weightedSites.stream().map(s -> PGS.coordFromPVector(s)).toList();
var e = new Envelope(bounds[0], bounds[2], bounds[1], bounds[3]); // x,x,y,y
- AdditivelyWeightedVoronoi vd = new AdditivelyWeightedVoronoi(GEOM_FACTORY, 0.2);
+ AdditivelyWeightedVoronoi vd = new AdditivelyWeightedVoronoi(GEOM_FACTORY, 0.25);
List cells = vd.computeCells(sites, e);
- if (forceConforming) {
- cells = PGS_Meshing.fixBreaks(cells, 1); // gapWidth==1 suits errTol==0.2
- var simple = CoverageSimplifier.simplifyInner(cells.toArray(Geometry[]::new), 1);
- cells = Arrays.asList(simple);
- } else {
- // Produces rather dense output, so simplify using conservative DCE relevance.
- final DCETerminationCallback dceCallback = (currentVertex, relevance, verticesRemaining) -> relevance >= 15;
-
- cells = cells.stream().map(cell -> {
- var ring = DiscreteCurveEvolution.process((LineString) cell.getBoundary(), dceCallback);
- return ring;
- }).toList();
- }
+ // Produces rather dense output, so simplify using conservative DCE relevance
+ final DCETerminationCallback dceCallback = (currentVertex, relevance, verticesRemaining) -> relevance >= 10;
+ cells = cells.stream().map(cell -> {
+ var ring = DiscreteCurveEvolution.process((LineString) cell.getBoundary(), dceCallback);
+ return ring;
+ }).toList();
var awvd = toPShape(cells);
PGS_Conversion.setAllFillColor(awvd, Colors.WHITE);
diff --git a/src/main/java/micycle/pgs/commons/AdditivelyWeightedVoronoi.java b/src/main/java/micycle/pgs/commons/AdditivelyWeightedVoronoi.java
index da396d13..8082607d 100644
--- a/src/main/java/micycle/pgs/commons/AdditivelyWeightedVoronoi.java
+++ b/src/main/java/micycle/pgs/commons/AdditivelyWeightedVoronoi.java
@@ -2,494 +2,1120 @@
import java.util.ArrayList;
import java.util.Arrays;
+import java.util.Collections;
+import java.util.HashMap;
import java.util.List;
-import java.util.Objects;
+import java.util.Map;
import java.util.stream.IntStream;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
+import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.index.ItemVisitor;
+import org.locationtech.jts.index.hprtree.HPRtree;
import net.jafama.FastMath;
/**
- * Approximates 2D additively-weighted Voronoi cells (Apollonius diagram) for
- * input sites encoded as {@link Coordinate}s: (x,y) is the site center and z is
- * the additive weight (radius).
+ * Computes 2D additively-weighted Voronoi cells (the Apollonius diagram)
+ * for input sites encoded as {@link Coordinate}s: (x,y) is the site centre and
+ * z is the additive weight (radius). The distance to site i is
+ * {@code d_i(p) = |p - p_i| - w_i}; cell i is the locus where no other
+ * site is closer.
*
- *
- * Each cell boundary is represented in polar form around the site center as a
- * "lifted radius" function t(θ). For each direction θ, competitors induce a
- * candidate radius; the boundary is the lower envelope (minimum) over all
- * competitors. The envelope is approximated by angular sampling with adaptive
- * refinement. This idea is based on: Hans-Martin Will, “Practical and
- * efficient computation of additively weighted Voronoi cells for applications
- * in molecular biology”, 1998.
+ *
Formulation
*
- *
- * Cells are bounded to a rectangular clip envelope. For each sampled direction,
- * the radius is upper-bounded by the ray's exit distance from the rectangle.
- * This produces a polygon that is already confined to the clip rectangle.
+ * Around site i, a competitor j bounds the cell in direction
+ * {@code u(t) = (cos t, sin t)} at radius
*
- *
- * Dominated (empty) cells are omitted from the returned list.
+ *
+ * t_j(t) = N / (2 (dw + <u, c>)), c = p_j - p_i, dw = w_j - w_i, N = |c|² - dw²
+ *
+ *
+ * and the cell boundary is the lower envelope {@code min_j t_j}. Rather than
+ * sampling that envelope, this implementation works with the reciprocal radius:
+ *
+ *
+ * rho_j(t) = 1 / t_j(t) = (dw + c_x cos t + c_y sin t) / (N/2)
+ *
+ *
+ * which is affine in (cos t, sin t), hence a sinusoid. The lower envelope of
+ * conics becomes an upper envelope of sinusoids,
+ * {@code t(t) = 1 / max_j rho_j(t)}.
*
- *
- * Notes:
*
- * - For best results, the clip envelope should contain all sites.
- * - This implementation returns polygon shells only (no holes).
+ * - Two sinusoids cross at most twice, so every envelope breakpoint is
+ * available in closed form and nothing is discovered by sampling.
+ * - A clip half-plane contributes the same kind of reciprocal-radius
+ * constraint.
+ * - For a non-dominated pair N is positive, and non-positive reciprocal
+ * radius is naturally ignored by the upper envelope.
*
- *
+ *
+ * Structure
+ *
+ *
+ * - Coincident sites are deduplicated.
+ * - An R-tree over disk bounding boxes supplies a certified superset of the
+ * competitors that can bind each cell.
+ * - The exact sinusoidal envelope determines the symbolic boundary.
+ * - Vertices are keyed by cyclically ordered triples of constraints.
+ * - Edges are keyed by their curve and symbolic endpoints, globally interned,
+ * and flattened once.
+ * - Incident cells consume the same flattened edge in opposite directions,
+ * making shared polygon boundaries watertight by construction.
+ *
+ *
+ *
+ * JTS stores the resulting polygons. If sites lie outside a supplied clipping
+ * envelope, cells are first built in an enclosing polar box and JTS performs
+ * the final exceptional trim against the requested envelope.
+ *
* @author Michael Carleton
*/
public class AdditivelyWeightedVoronoi {
- private static class AV2DOptions {
-
- /**
- * Adaptive sampling error tolerance (squared). Smaller values yield smoother
- * curves at the cost of more vertices.
- */
- public double errorToleranceSq = 0.01 * 0.01;
-
- /** Maximum recursion depth for adaptive sampling. */
- public int maxRecursionDepth = 10;
-
- /** Scale factor for the default clip envelope when none is provided. */
- public double boundsScale = 1.5;
-
- /**
- * Small epsilon for denominator/degeneracy checks. This is an absolute epsilon;
- * for very large coordinate magnitudes, consider increasing it.
- */
- public double eps = 1e-12;
- }
+ private static final int CLIP_XMAX = -1, CLIP_XMIN = -2, CLIP_YMAX = -3, CLIP_YMIN = -4;
+ private static final double TWO_PI = 2 * Math.PI;
+ private static final double ANG_EPS = 1e-12;
+ private static final double MAX_PATCH_SPAN = Math.PI / 3;
+ private static final double BOUNDS_SCALE = 1.5;
private final GeometryFactory gf;
- private final AV2DOptions opt;
+ private final double tol;
/**
- * Creates an instance using the given geometry factory and an adaptive sampling
- * tolerance.
- *
* @param gf geometry factory used to create output polygons
- * @param errorTolerance maximum allowed deviation (in coordinate units) used by
- * the adaptive sampler; smaller values produce smoother
- * cell boundaries at the cost of more vertices
+ * @param errorTolerance maximum sagitta of each emitted conic chord
*/
public AdditivelyWeightedVoronoi(GeometryFactory gf, double errorTolerance) {
+ if (gf == null) {
+ throw new IllegalArgumentException("GeometryFactory must not be null");
+ }
this.gf = gf;
- this.opt = new AV2DOptions();
- this.opt.errorToleranceSq = errorTolerance * errorTolerance;
+ tol = Double.isFinite(errorTolerance) && errorTolerance > 0 ? errorTolerance : 1e-3;
}
/**
- * Computes clipped polygonal approximations of additively-weighted Voronoi
- * cells.
- *
- *
- * Input sites are provided as {@link Coordinate}s where (x,y) is the site
- * location and z is its additive weight. Cells are bounded to the supplied
- * {@code bounds}, or to a default envelope derived from the data if
- * {@code bounds} is {@code null}.
- *
- *
- * Dominated (empty) cells are omitted from the returned list. Each returned
- * polygon has its {@linkplain Polygon#getUserData() userData} set to the
- * corresponding input site index.
+ * Computes clipped polygonal approximations of the weighted Voronoi cells.
*
- * @param sites input sites; {@code (x,y)} is the site center and {@code z} is
- * the additive weight (non-finite z is treated as 0)
- * @param bounds optional clip envelope; if {@code null}, a default is computed
- * from the input extent and weight magnitudes
- * @return list of polygons for non-empty cells; order is unspecified due to
- * parallel execution
+ * @param sites input sites; (x,y) is the centre and z the additive weight
+ * @param bounds optional clipping envelope
+ * @return non-empty cells in ascending site order, carrying the site index as
+ * user data
*/
public List computeCells(List sites, Envelope bounds) {
if (sites == null || sites.isEmpty()) {
return List.of();
}
+ return new Ctx(sites, bounds).compute();
+ }
+
+ // ------------------------------------------------------------------
+ // per-invocation state
+ // ------------------------------------------------------------------
+
+ private final class Ctx {
+
+ final int n;
+ final double[] x, y, w;
+ final boolean[] dead;
+ final Envelope clip;
+ final double bMinX, bMaxX, bMinY, bMaxY;
+ final boolean trim;
+ final double distEps, seedRadius, boxSpan;
+ final DiskIndex index;
+
+ Ctx(List sites, Envelope bounds) {
+ n = sites.size();
+ x = new double[n];
+ y = new double[n];
+ w = new double[n];
+
+ Envelope dataEnv = new Envelope();
+ double maxAbsW = 0, maxAbsC = 0;
+
+ for (int i = 0; i < n; i++) {
+ Coordinate c = sites.get(i);
+ if (c == null || !Double.isFinite(c.x) || !Double.isFinite(c.y)) {
+ throw new IllegalArgumentException("Site " + i + " has non-finite coordinates");
+ }
+ x[i] = c.x;
+ y[i] = c.y;
+ w[i] = Double.isFinite(c.getZ()) ? c.getZ() : 0;
+ maxAbsW = Math.max(maxAbsW, Math.abs(w[i]));
+ maxAbsC = Math.max(maxAbsC, Math.max(Math.abs(x[i]), Math.abs(y[i])));
+ dataEnv.expandToInclude(x[i], y[i]);
+ }
+
+ clip = bounds == null ? defaultClip(dataEnv, maxAbsW) : new Envelope(bounds);
+ if (clip.isNull() || !(clip.getWidth() > 0) || !(clip.getHeight() > 0) || !Double.isFinite(clip.getMinX()) || !Double.isFinite(clip.getMaxX())
+ || !Double.isFinite(clip.getMinY()) || !Double.isFinite(clip.getMaxY())) {
+ throw new IllegalArgumentException("Bounds must be finite and non-empty");
+ }
+
+ double scale = Math.max(Math.max(clip.getWidth(), clip.getHeight()), maxAbsC);
+ if (!(scale > 0) || !Double.isFinite(scale)) {
+ scale = 1;
+ }
+ distEps = 1e-12 * scale;
+
+ /*
+ * The polar form requires every site strictly inside its box. When the
+ * requested clip does not provide that, use an enclosing box and trim the
+ * finished cells afterward.
+ */
+ double margin = 1e-6 * scale + tol;
+ boolean outside = false;
+ for (int i = 0; i < n && !outside; i++) {
+ outside = x[i] <= clip.getMinX() + margin || x[i] >= clip.getMaxX() - margin || y[i] <= clip.getMinY() + margin
+ || y[i] >= clip.getMaxY() - margin;
+ }
+
+ Envelope box = new Envelope(clip);
+ if (outside) {
+ box.expandToInclude(dataEnv);
+ box.expandBy(margin);
+ }
+
+ trim = outside;
+ bMinX = box.getMinX();
+ bMaxX = box.getMaxX();
+ bMinY = box.getMinY();
+ bMaxY = box.getMaxY();
+ boxSpan = box.getWidth() + box.getHeight();
+ seedRadius = 1.5 * Math.sqrt(box.getWidth() * box.getHeight() / Math.max(1, n));
+ index = new DiskIndex(x, y, w);
+
+ /*
+ * Keep the heaviest coincident site, using the lowest index on ties.
+ */
+ dead = new boolean[n];
+ IntStream.range(0, n).parallel().forEach(i -> {
+ IntList near = new IntList(8);
+ index.query(x[i], y[i], 0, near);
+ for (int t = 0; t < near.size; t++) {
+ int j = near.get(t);
+ if (j != i && x[j] == x[i] && y[j] == y[i] && (w[j] > w[i] || w[j] == w[i] && j < i)) {
+ dead[i] = true;
+ return;
+ }
+ }
+ });
+ }
+
+ /**
+ * Builds symbolic cells in parallel, interns their shared geometry, flattens
+ * each distinct edge once, then assembles the polygons.
+ */
+ List compute() {
+ SCell[] cells = IntStream.range(0, n).parallel().mapToObj(this::symbolicCell).toArray(SCell[]::new);
+ Graph graph = new Graph();
+ for (SCell cell : cells) {
+ if (cell != null) {
+ graph.add(cell);
+ }
+ }
+ graph.validate();
+ graph.flatten();
+
+ return IntStream.range(0, n).parallel().mapToObj(i -> polygons(cells[i], graph)).flatMap(List::stream).toList();
+ }
+
+ private final class Graph {
+
+ final Map vertices = new HashMap<>();
+ final Map edges = new HashMap<>();
+
+ void add(SCell cell) {
+ for (Use use : cell.boundary) {
+ EKey key = use.key;
+ double[] p0 = vertices.computeIfAbsent(key.v0, Ctx.this::vertex);
+ double[] p1 = vertices.computeIfAbsent(key.v1, Ctx.this::vertex);
+ Edge edge = edges.computeIfAbsent(key, k -> new Edge(k, p0, p1));
+ edge.uses++;
+ if (use.forward) {
+ edge.forwardUses++;
+ }
+ }
+ }
+
+ void validate() {
+ for (Edge edge : edges.values()) {
+ if (edge.key.isClip()) {
+ if (edge.uses != 1) {
+ throw failure("Clip edge has " + edge.uses + " uses: " + edge.key);
+ }
+ } else if (edge.uses != 2 || edge.forwardUses != 1) {
+ throw failure("Inconsistent shared edge: " + edge.key);
+ }
+ }
+ }
+
+ void flatten() {
+ edges.values().parallelStream().forEach(e -> e.line = flattenEdge(e));
+ }
+ }
+
+ // --------------------------------------------------------------
+ // certified candidates and sinusoidal envelope
+ // --------------------------------------------------------------
+
+ private Cell buildCell(int i) {
+ IntList ids = new IntList(16);
+ seed(i, ids);
+
+ for (int round = 0; round < 64; round++) {
+ Constraint[] constraints = constraints(i, ids);
+ if (constraints == null) {
+ return null;
+ }
+
+ Arrays.sort(constraints);
+ Cell cell = envelope(constraints);
+ double radius = cell.maxRadius();
+
+ IntList query = new IntList(8);
+ query(i, radius, query);
+ boolean grew = false;
+
+ for (int t = 0; t < query.size; t++) {
+ int j = query.get(t);
+ if (!ids.contains(j)) {
+ ids.add(j);
+ grew = true;
+ }
+ }
+ if (!grew) {
+ return cell;
+ }
+ }
+ throw failure("Candidate certification did not converge for site " + i);
+ }
- final int n = sites.size();
- final double[] x = new double[n];
- final double[] y = new double[n];
- final double[] w = new double[n];
+ private Constraint[] constraints(int i, IntList ids) {
+ List out = new ArrayList<>(ids.size + 4);
+ out.add(clipConstraint(i, CLIP_XMAX));
+ out.add(clipConstraint(i, CLIP_XMIN));
+ out.add(clipConstraint(i, CLIP_YMAX));
+ out.add(clipConstraint(i, CLIP_YMIN));
+
+ for (int t = 0; t < ids.size; t++) {
+ int j = ids.get(t);
+ if (j == i || dead[j]) {
+ continue;
+ }
- final Envelope dataEnv = new Envelope();
- double maxAbsW = 0.0;
+ double dx = x[j] - x[i], dy = y[j] - y[i], dw = w[j] - w[i];
+ double d = FastMath.hypot(dx, dy);
+ if (dw >= d - distEps) {
+ return null;
+ }
- for (int i = 0; i < n; i++) {
- Coordinate c = sites.get(i);
- x[i] = c.getX();
- y[i] = c.getY();
- double zi = c.getZ();
- w[i] = Double.isFinite(zi) ? zi : 0.0;
- maxAbsW = Math.max(maxAbsW, Math.abs(w[i]));
- dataEnv.expandToInclude(x[i], y[i]);
+ double den = 0.5 * (d * d - dw * dw);
+ if (den > 0) {
+ out.add(new Constraint(j, dw, dx, dy, den));
+ }
+ }
+ return out.toArray(Constraint[]::new);
}
- final Envelope clipEnv = bounds != null ? new Envelope(bounds) : defaultClipEnvelope(dataEnv, maxAbsW);
+ /**
+ * Opens the certified loop with a cheap radius guess. Nothing here needs to be
+ * complete; buildCell certifies the final candidate set.
+ */
+ private void seed(int i, IntList out) {
+ for (double radius = seedRadius;; radius *= 2) {
+ out.clear();
+ index.query(x[i], y[i], radius, out);
+
+ int usable = 0;
+ for (int t = 0; t < out.size; t++) {
+ int j = out.get(t);
+ if (j != i && !dead[j]) {
+ usable++;
+ }
+ }
+ if (usable >= 4) {
+ return;
+ }
+ if (!Double.isFinite(radius) || radius >= boxSpan) {
+ out.clear();
+ index.all(out);
+ return;
+ }
+ }
+ }
- final double minX = clipEnv.getMinX(), maxX = clipEnv.getMaxX();
- final double minY = clipEnv.getMinY(), maxY = clipEnv.getMaxY();
+ /**
+ * All j that can bind within radius r: {@code (d - w_j + w_i)/2 < r}.
+ */
+ private void query(int i, double radius, IntList out) {
+ IntList raw = new IntList(16);
+ double expand = Double.isFinite(radius) ? Math.max(0, 2 * radius - w[i]) : Double.POSITIVE_INFINITY;
+
+ if (Double.isFinite(expand)) {
+ index.query(x[i], y[i], expand, raw);
+ } else {
+ index.all(raw);
+ }
- final Precompute pc = precomputePairs(x, y, w);
+ for (int t = 0; t < raw.size; t++) {
+ int j = raw.get(t);
+ if (j == i || dead[j]) {
+ continue;
+ }
+ double d = FastMath.hypot(x[j] - x[i], y[j] - y[i]);
+ if (0.5 * (d - w[j] + w[i]) < radius) {
+ out.add(j);
+ }
+ }
+ }
- List out = IntStream.range(0, n).parallel().mapToObj(i -> {
- if (pc.empty[i]) {
+ // --------------------------------------------------------------
+ // symbolic topology
+ // --------------------------------------------------------------
+
+ private SCell symbolicCell(int i) {
+ if (dead[i]) {
return null;
}
- // per-cell prune bound: farthest corner distance from site to clip rectangle
- final double Rsite = radiusToCoverEnvelopeFromSite(clipEnv, x[i], y[i]);
+ Cell cell = buildCell(i);
+ if (cell == null) {
+ return null;
+ }
+ if (cell.k < 2) {
+ throw failure("Degenerate envelope for site " + i);
+ }
- Polygon approx = approximateCellPolygon(i, x, y, w, pc.candidates[i], Rsite, minX, maxX, minY, maxY);
- approx.setUserData(i);
- return approx;
- }).filter(Objects::nonNull).toList();
+ List
