2.3

CHANGELOG.md

@@ -5,6 +5,15 @@ All notable changes to PGS will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). Dates are *YYYY-MM-DD*.
 
+## **2.3** *(2026-xx-xx)*
+
+### Added
+* `fillPack()` to `PGS_CirclePacking`. Fills gaps in an existing circle packing using the Largest Empty Circle (LEC) algorithm, seeded with the existing circles as constraints.
+
+### Changes
+* Reimplemented `PGS_CirclePacking.frontChainPack()` to pack circles directly within the shape boundary, rather than packing the bounding envelope and filtering.
+* Reimplemented `PGS_Voronoi.additivelyWeightedVoronoi()`. About 10x faster than before. The `forceConforming` argument has been removed, as it's no longer necessary.
+
 ## **2.2** *(2026-04-06)*
 
 ### Added

README.md

@@ -750,10 +750,12 @@ Much of the functionality (but by no means all) is demonstrated below:
   </tr>
 
   <tr>
-    <td align="center" valign="center" colspan="2"><b>Obstacle</td>
+    <td align="center" valign="center"><b>Obstacle</td>
+    <td align="center" valign="center"><b>Fill</td>
   </tr>
   <tr>
     <td valign="top"><img src="resources/circle_packing/obstaclePack.gif"></td>
+    <td valign="top"><img src="resources/circle_packing/fillPack.gif"></td>
   </tr>
 </table>
 

pom.xml

@@ -354,7 +354,7 @@
 		<dependency>
 			<groupId>com.github.micycle1</groupId>
 			<artifactId>BetterBeziers</artifactId>
-			<version>1.0</version>
+			<version>1.1</version>
 		</dependency>
 		<dependency>
 			<groupId>com.github.micycle1</groupId>
@@ -385,8 +385,8 @@
 		</dependency>
 		<dependency>
 			<groupId>com.github.micycle1</groupId>
-			<artifactId>geoblitz</artifactId>
-			<version>0.9.6</version>
+			<artifactId>GeoBlitz</artifactId>
+			<version>0.9.9</version>
 		</dependency>
 		<dependency>
 			<groupId>com.github.micycle1</groupId>

src/main/java/micycle/pgs/PGS.java

@@ -167,7 +167,7 @@ static final Coordinate[] toCoords(final Collection<PVector> points) {
 	}
 
 	static final PVector toPVector(Coordinate c) {
-		return new PVector((float) c.x, (float) c.y);
+		return new PVector((float) c.x, (float) c.y, (float) c.z);
 	}
 
 	/**

src/main/java/micycle/pgs/PGS_CirclePacking.java

@@ -14,15 +14,17 @@
 import org.locationtech.jts.geom.Envelope;
 import org.locationtech.jts.geom.Geometry;
 import org.locationtech.jts.geom.Location;
+import org.locationtech.jts.geom.Polygon;
 import org.locationtech.jts.operation.distance.IndexedFacetDistance;
 import org.tinfour.common.IIncrementalTin;
 import org.tinfour.common.SimpleTriangle;
 import org.tinfour.common.Vertex;
+
+import com.github.micycle1.geoblitz.CircleIndex;
 import com.github.micycle1.geoblitz.PointDistanceIndex;
 import com.github.micycle1.geoblitz.YStripesPointInAreaLocator;
 
-import micycle.pgs.commons.CircleCoverTree;
-import micycle.pgs.commons.FrontChainPacker;
+import micycle.pgs.commons.GeometryFrontChainPacker;
 import micycle.pgs.commons.LargestEmptyCircles;
 import micycle.pgs.commons.RepulsionCirclePack;
 import micycle.pgs.commons.TangencyPack;
@@ -191,7 +193,7 @@
 	 *         the center point, and .z represents the radius.
 	 */
 	public static List<PVector> stochasticPack(final PShape shape, final int points, final double minRadius, boolean triangulatePoints, long seed) {
-		CircleCoverTree<PVector> tree = new CircleCoverTree<>();
+		CircleIndex<PVector> tree = new CircleIndex<>();
 
 		List<PVector> steinerPoints = PGS_Processing.generateRandomPoints(shape, points, seed);
 		if (triangulatePoints) {
@@ -266,32 +268,12 @@
 	 *         the center point and .z represents radius.
 	 */
 	public static List<PVector> frontChainPack(PShape shape, double radiusMin, double radiusMax, long seed) {
-		radiusMin = Math.max(1f, Math.min(radiusMin, radiusMax)); // choose min and constrain
-		radiusMax = Math.max(1f, Math.max(radiusMin, radiusMax)); // choose max and constrain
-		final Geometry g = fromPShape(shape);
-		final Envelope e = g.getEnvelopeInternal();
-		YStripesPointInAreaLocator pointLocator;
-
-		final FrontChainPacker packer = new FrontChainPacker((float) e.getWidth(), (float) e.getHeight(), (float) radiusMin, (float) radiusMax,
-				(float) e.getMinX(), (float) e.getMinY(), seed);
-
-		if (radiusMin == radiusMax) {
-			// if every circle same radius, use faster contains check
-			pointLocator = new YStripesPointInAreaLocator(g.buffer(radiusMax));
-			packer.getCircles().removeIf(p -> pointLocator.locate(PGS.coordFromPVector(p)) == Location.EXTERIOR);
-		} else {
-			pointLocator = new YStripesPointInAreaLocator(g);
-			IndexedFacetDistance distance = new IndexedFacetDistance(g);
-			packer.getCircles().removeIf(p -> {
-				// first test whether shape contains circle center point (somewhat faster)
-				if (pointLocator.locate(PGS.coordFromPVector(p)) != Location.EXTERIOR) {
-					return false; // keep if interior
-				}
-				return !distance.isWithinDistance(PGS.pointFromPVector(p), p.z * (2 / 3d));
-			});
-		}
-
-		return packer.getCircles();
+		radiusMin = Math.max(1d, Math.min(radiusMin, radiusMax)); // choose min and constrain
+		radiusMax = Math.max(1d, Math.max(radiusMin, radiusMax)); // choose max and constrain
+
+		final var packer = new GeometryFrontChainPacker((Polygon) fromPShape(shape), radiusMin, radiusMax, seed);
+
+		return packer.getCircles().stream().map(c -> PGS.toPVector(c)).toList();
 	}
 
 	/**
@@ -356,6 +338,83 @@
 		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.
+	 * <p>
+	 * 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 <b>new</b> 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<PVector> fillPack(PShape shape, Collection<PVector> existingCircles, double minRadius, double tolerance) {
+		tolerance = Math.max(0.01, tolerance);
+		minRadius = Math.max(0.01, minRadius);
+
+		final List<Coordinate> 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<PVector> 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.
+	 * <p>
+	 * 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 <b>new</b> 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<PVector> fillPack(PShape shape, Collection<PVector> existingCircles, int n, double tolerance) {
+		tolerance = Math.max(0.01, tolerance);
+
+		final List<Coordinate> 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<PVector> 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.

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<PVector> coords = new ArrayList<>(samples.length);
 		for (double[] sample : samples) {
 			coords.add(new PVector((float) sample[0], (float) sample[1]));

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<PVector> 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);
 	}

src/main/java/micycle/pgs/PGS_Voronoi.java

@@ -6,14 +6,14 @@
 
 import java.awt.geom.Rectangle2D;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Collection;
 import java.util.HashMap;
 import java.util.List;
 import java.util.Objects;
 import java.util.stream.Collectors;
 
 import org.locationtech.jts.coverage.CoverageSimplifier;
 import org.locationtech.jts.coverage.CoverageUnion;
 import org.locationtech.jts.densify.Densifier;
 import org.locationtech.jts.geom.Coordinate;
@@ -483,43 +483,6 @@
 		return voronoiCells;
 	}
 
-	/**
-	 * Generates an <b>additively weighted Voronoi diagram</b> (AWVD) for a set of
-	 * weighted point sites, clipped to the provided bounding rectangle.
-	 * <p>
-	 * 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:
-	 * 
-	 * <pre>
-	 *   d(p, s) = ||p - s|| - w
-	 * </pre>
-	 * 
-	 * 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
-	 * <i>curved</i> (hyperbolic arcs), and some sites may end up with empty cells
-	 * depending on weights and configuration.
-	 * <p>
-	 * Each input {@link PVector} encodes one site where:
-	 * <ul>
-	 * <li>{@code (.x, .y)} is the site coordinate</li>
-	 * <li>{@code .z} is the site's weight (in the same units as {@code x/y})</li>
-	 * </ul>
-	 *
-	 * @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<PVector> weightedSites, double[] bounds) {
-		return additivelyWeightedVoronoi(weightedSites, bounds, false);
-	}
-
 	/**
 	 * Generates an <b>additively weighted Voronoi diagram</b> (AWVD) for a set of
 	 * weighted point sites, clipped to the provided bounding rectangle.
@@ -551,38 +514,28 @@
 	 * simplify the interior boundaries.</li>
 	 * </ul>
 	 *
-	 * @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<PVector> weightedSites, double[] bounds, boolean forceConforming) {
+	public static PShape additivelyWeightedVoronoi(Collection<PVector> 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<? extends Geometry> 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);