shapes, Coloring
break;
case RLF_BRUTE_FORCE_4COLOR :
int iterations = 0;
- long seed = 1337;
+ long seed = SEED; // init as default
do {
coloring = new RLFColoring<>(graph, seed).getColoring();
- seed = ThreadLocalRandom.current().nextLong();
+ seed = ThreadLocalRandom.current().nextLong(); // randomise seed
iterations++;
} while (coloring.getNumberColors() > 4 && iterations < 250);
break;
diff --git a/src/main/java/micycle/pgs/PGS_Construction.java b/src/main/java/micycle/pgs/PGS_Construction.java
index 6b7d5fd3..f520b25c 100644
--- a/src/main/java/micycle/pgs/PGS_Construction.java
+++ b/src/main/java/micycle/pgs/PGS_Construction.java
@@ -796,7 +796,7 @@ public static PShape createRandomSFCurve(int nColumns, int nRows, double cellWid
* @param cellWidth visual/pixel width of each cell
* @param cellHeight visual/pixel width of each cell
* @param seed random seed
- * @return a stroked PATH PShape
+ * @return a mitered stroked PATH PShape
* @see #createRandomSFCurve(int, int, double, double)
* @since 1.4.0
*/
diff --git a/src/main/java/micycle/pgs/PGS_Contour.java b/src/main/java/micycle/pgs/PGS_Contour.java
index 10456f2e..2dfa23ac 100644
--- a/src/main/java/micycle/pgs/PGS_Contour.java
+++ b/src/main/java/micycle/pgs/PGS_Contour.java
@@ -14,9 +14,12 @@
import java.util.Map;
import java.util.Set;
import java.util.stream.Collectors;
+import java.util.stream.Stream;
import javax.vecmath.Point3d;
+import org.jgrapht.alg.interfaces.ShortestPathAlgorithm;
+import org.jgrapht.alg.shortestpath.BFSShortestPath;
import org.jgrapht.graph.DefaultEdge;
import org.jgrapht.graph.SimpleGraph;
import org.joml.Vector2d;
@@ -232,15 +235,26 @@ public static PShape chordalAxis(PShape shape) {
* consisting of straight-line segments only. Roughly, it is the geometric graph
* whose edges are the traces of vertices of shrinking mitered offset curves of
* the polygon.
- *
+ *
+ * For a single polygon, this method returns a GROUP PShape containing three
+ * children:
+ *
+ * - Child 0: GROUP PShape consisting of skeleton faces.
+ * - Child 1: LINES PShape representing branches, which are lines connecting
+ * the skeleton to the polygon's edge.
+ * - Child 2: LINES PShape composed of bones, depicting the pure straight
+ * skeleton of the polygon.
+ *
+ *
+ * For multi-polygons, the method returns a master GROUP PShape. This master
+ * shape includes multiple skeleton GROUP shapes, each corresponding to a single
+ * polygon and structured as described above.
+ *
* @param shape a single polygon (that can contain holes), or a multi polygon
* (whose polygons can contain holes)
- * @return when the input is a single polygon, returns a GROUP PShape containing
- * 3 children: child 1 = GROUP PShape of skeleton faces; child 2 = LINES
- * PShape of branches (lines that connect skeleton to edge); child 3 =
- * LINES PShape of bones (the pure straight skeleton). For
- * multi-polygons, a master GROUP shape of skeleton GROUP shapes
- * (described above) is returned.
+ *
+ * @return PShape based on the input polygon structure, either as a single or
+ * multi-polygon skeleton representation.
*/
public static PShape straightSkeleton(PShape shape) {
final Geometry g = fromPShape(shape);
@@ -255,25 +269,7 @@ public static PShape straightSkeleton(PShape shape) {
return shape;
}
- /**
- *
- * @param polygon a single polygon that can contain holes
- * @return
- */
private static PShape straightSkeleton(Polygon polygon) {
- /*
- * Kenzi implementation (since PGS 1.3.0) is much faster (~50x!) but can fail on
- * more complicated inputs. Therefore try Kenzi implementation first, but fall
- * back to Twak implementation if it fails.
- */
- try {
- return straightSkeletonKendzi(polygon);
- } catch (Exception e) {
- return straightSkeletonTwak(polygon);
- }
- }
-
- private static PShape straightSkeletonTwak(Polygon polygon) {
if (polygon.getCoordinates().length > 1000) {
polygon = (Polygon) DouglasPeuckerSimplifier.simplify(polygon, 2);
}
@@ -301,7 +297,7 @@ private static PShape straightSkeletonTwak(Polygon polygon) {
skeleton.skeleton(); // compute skeleton
skeleton.output.faces.values().forEach(f -> {
- final List vertices = f.getLoopL().iterator().next().asList();
+ List vertices = f.getLoopL().iterator().next().stream().toList();
List faceVertices = new ArrayList<>();
for (int i = 0; i < vertices.size(); i++) {
@@ -350,95 +346,20 @@ private static PShape straightSkeletonTwak(Polygon polygon) {
return lines;
}
- private static PShape straightSkeletonKendzi(Polygon polygon) {
- final LinearRing[] rings = new LinearRingIterator(polygon).getLinearRings();
- Set edgeCoordsSet = new HashSet<>();
- final List points = ringToVec(rings[0], edgeCoordsSet);
- final List> holes = new ArrayList<>();
- for (int i = 1; i < rings.length; i++) {
- holes.add(ringToVec(rings[i], edgeCoordsSet));
- }
-
- final SkeletonOutput so = kendzi.math.geometry.skeleton.Skeleton.skeleton(points, holes, new SkeletonConfiguration());
- final PShape skeleton = new PShape(PConstants.GROUP);
- final PShape faces = new PShape(PConstants.GROUP);
- /*
- * Create PEdges first to prevent lines being duplicated in output shapes since
- * faces share branches and bones.
- */
- final Set branchEdges = new HashSet<>();
- final Set boneEdges = new HashSet<>();
- so.getFaces().forEach(f -> {
- /*
- * q stores the index of second vertex of the face that is a shape vertex. This
- * is used to rotate f.getPoints() so that the vertices of every face PShape
- * begin at the shape edge.
- */
- int q = 0;
- for (int i = 0; i < f.getPoints().size(); i++) {
- final Vector2dc p1 = f.getPoints().get(i);
- final Vector2dc p2 = f.getPoints().get((i + 1) % f.getPoints().size());
- final boolean a = edgeCoordsSet.contains(p1);
- final boolean b = edgeCoordsSet.contains(p2);
- if (a ^ b) { // branch (xor)
- branchEdges.add(new PEdge(p1.x(), p1.y(), p2.x(), p2.y()));
- q = i;
- } else {
- if (!a) { // bone
- boneEdges.add(new PEdge(p1.x(), p1.y(), p2.x(), p2.y()));
- } else {
- q = i;
- }
- }
- }
-
- List faceVertices = new ArrayList<>(f.getPoints().size());
- Collections.rotate(f.getPoints(), -q + 1);
- f.getPoints().forEach(p -> faceVertices.add(new PVector((float) p.x(), (float) p.y())));
-
- PShape face = PGS_Conversion.fromPVector(faceVertices);
- face.setStroke(true);
- face.setStrokeWeight(2);
- face.setStroke(ColorUtils.composeColor(147, 112, 219));
- faces.addChild(face);
- });
-
- final PShape bones = prepareLinesPShape(null, null, 4);
- boneEdges.forEach(e -> {
- bones.vertex(e.a.x, e.a.y);
- bones.vertex(e.b.x, e.b.y);
- });
- bones.endShape();
-
- final PShape branches = prepareLinesPShape(ColorUtils.composeColor(40, 235, 180), null, null);
- branchEdges.forEach(e -> {
- branches.vertex(e.a.x, e.a.y);
- branches.vertex(e.b.x, e.b.y);
- });
- branches.endShape();
-
- skeleton.addChild(faces);
- skeleton.addChild(branches);
- skeleton.addChild(bones);
-
- return skeleton;
- }
-
/**
- * Generates a topographic-like isoline contour map from the shape's vertices.
- * The "elevation" (or z value) of points is the euclidean distance between a
- * point in the shape and the given "high" point.
+ * Generates a topographic-like isoline contour map from the shape's vertices
+ * and a given "high point". Isolines represent the "elevation", or euclidean
+ * distance, between a location in the shape and the "high point".
*
* Assigns each point feature a number equal to the distance between geometry's
* centroid and the point.
*
- * @param shape
+ * @param shape the bounds in which to draw isolines
* @param highPoint position of "high" point within the shape
* @param intervalSpacing distance between successive isolines
- * @return PShape containing isolines linework
+ * @return PShape containing isolines linework
*/
public static PShape isolines(PShape shape, PVector highPoint, double intervalSpacing) {
-
/*
* Also See:
* https://github.com/hageldave/JPlotter/blob/master/jplotter/src/main/java/
@@ -636,6 +557,55 @@ public static PShape distanceField(PShape shape, double spacing) {
return i;
}
+ /**
+ * Generates a tree structure representing the shortest paths from a given start
+ * point to all other vertices in the provided mesh. The paths are computed
+ * using the existing connectivity of the mesh edges, ensuring that the
+ * shortest-path tree respects the original mesh structure. The tree is
+ * constructed using a Breadth-First Search (BFS) algorithm.
+ *
+ * The shortest-path tree represents the minimal set of mesh edges required to
+ * connect the start point to all other vertices in the mesh, following the
+ * mesh's inherent connectivity. This ensures that the paths are constrained by
+ * the mesh's topology rather than creating arbitrary connections between
+ * vertices.
+ *
+ * If the provided start point does not exactly match a vertex in the mesh, the
+ * closest vertex in the mesh to the start point is used as the actual starting
+ * point for the shortest-path computation.
+ *
+ * @param mesh A GROUP shape representing a mesh from which the graph is
+ * constructed. The mesh defines the connectivity between
+ * vertices via its edges.
+ * @param startPoint The starting point from which the shortest paths are
+ * calculated. If this point does not exactly match a vertex
+ * in the mesh, the closest vertex in the mesh will be used as
+ * the starting point.
+ * @return A PShape object representing the tree of shortest paths from the
+ * start point to all other vertices in the mesh. The paths are
+ * visualized with a semi-transparent pink stroke and are constrained by
+ * the mesh's edge connectivity.
+ * @since 2.1
+ */
+ public static PShape distanceTree(PShape mesh, PVector source) {
+ var g = PGS_Conversion.toGraph(mesh);
+ ShortestPathAlgorithm spa = new BFSShortestPath<>(g);
+
+ final PVector sourceActual = PGS_Optimisation.closestPoint(g.vertexSet(), source);
+ var paths = spa.getPaths(sourceActual);
+
+// var edges = g.vertexSet().stream().filter(v -> !v.equals(sourceActual)) // Exclude the source vertex
+// .flatMap(v -> paths.getPath(v).getEdgeList().stream()) // Flatten the edge lists into a single stream
+// .collect(Collectors.toSet()); // Collect the edges into a Set to remove duplicates
+
+ var pathLines = g.vertexSet().stream() //
+ .filter(v -> v != sourceActual) // Exclude the source vertex
+ .map(v -> PGS_Conversion.fromPVector(paths.getPath(v).getVertexList())).toList();
+ var group = PGS_Conversion.flatten(pathLines);
+ PGS_Conversion.setAllStrokeColor(group, ColorUtils.setAlpha(Colors.PINK, 50), 2);
+ return group;
+ }
+
/**
* Calculates the longest center line passing through a given shape (using
* default straightness weighting and smoothing parameters).
diff --git a/src/main/java/micycle/pgs/PGS_Conversion.java b/src/main/java/micycle/pgs/PGS_Conversion.java
index 51ba0103..109e47e3 100644
--- a/src/main/java/micycle/pgs/PGS_Conversion.java
+++ b/src/main/java/micycle/pgs/PGS_Conversion.java
@@ -71,6 +71,7 @@
import micycle.pgs.color.Colors;
import micycle.pgs.commons.Nullable;
import micycle.pgs.commons.PEdge;
+import net.jafama.FastMath;
import processing.core.PConstants;
import processing.core.PMatrix;
import processing.core.PShape;
@@ -900,24 +901,25 @@ public static List toPVector(PShape shape) {
/**
* Transforms a given PShape into a simple graph representation. In this
* representation, the vertices of the graph correspond to the vertices of the
- * shape, and the edges of the graph correspond to the edges of the shape. This
- * transformation is specifically applicable to polygonal shapes where edges are
- * formed by adjacent vertices.
+ * shape, and the edges of the graph correspond to the edges of the shape.
*
- * The edge weights in the graph are set to the length of the corresponding edge
- * in the shape.
+ * The edge weights in the graph are set to the length (euclidean distance) of
+ * the corresponding geometric edge in the shape.
*
- * @param shape the PShape to convert into a graph
+ * @param shape the PShape to convert into a graph. LINES and polygonal shapes
+ * are accepted (and GROUP shapes thereof).
* @return A SimpleGraph object that represents the structure of the input shape
* @since 1.3.0
* @see #toDualGraph(PShape)
*/
public static SimpleGraph toGraph(PShape shape) {
final SimpleGraph graph = new SimpleWeightedGraph<>(PEdge.class);
- for (PShape face : getChildren(shape)) {
- for (int i = 0; i < face.getVertexCount() - (face.isClosed() ? 0 : 1); i++) {
- final PVector a = face.getVertex(i);
- final PVector b = face.getVertex((i + 1) % face.getVertexCount());
+ for (PShape child : getChildren(shape)) {
+ final int stride = child.getKind() == PShape.LINES ? 2 : 1;
+ // Handle other child shapes (e.g., faces)
+ for (int i = 0; i < child.getVertexCount() - (child.isClosed() ? 0 : 1); i += stride) {
+ final PVector a = child.getVertex(i);
+ final PVector b = child.getVertex((i + 1) % child.getVertexCount());
if (a.equals(b)) {
continue;
}
@@ -947,23 +949,27 @@ public static PShape fromGraph(SimpleGraph graph) {
}
/**
- * Takes as input a graph and computes a layout for the graph vertices using a
- * Force-Directed placement algorithm (not vertex coordinates, if any exist).
- * Vertices are joined by their edges.
+ * Computes a layout for the vertices of a graph using a Force-Directed
+ * placement algorithm. The algorithm generates vertex coordinates based on the
+ * graph's topology, preserving its structure (i.e., connectivity and
+ * relationships between vertices and edges). Existing vertex coordinates, if
+ * any, are ignored.
*
- * The output is a rather abstract representation of the input graph, and not a
- * geometric equivalent (unlike most other conversion methods in the class).
+ * The output is an abstract representation of the input graph, not a geometric
+ * equivalent (unlike most other conversion methods in this class). The layout
+ * is bounded by the specified dimensions and anchored at (0, 0).
*
- * @param any vertex type
- * @param any edge type
- * @param graph the graph whose edges and vertices to lay out
- * @param normalizationFactor normalization factor for the optimal distance,
- * between 0 and 1.
- * @param boundsX horizontal vertex bounds
- * @param boundsY vertical vertex bounds
- * @return a GROUP PShape consisting of 2 children; child 0 is the linework
- * (LINES) depicting edges and child 1 is the points (POINTS) depicting
- * vertices. The bounds of the layout are anchored at (0, 0);
+ * @param the type of vertices in the graph
+ * @param the type of edges in the graph
+ * @param graph the graph whose vertices and edges are to be laid
+ * out
+ * @param normalizationFactor the normalization factor for the optimal distance
+ * between vertices, clamped between 0.001 and 1
+ * @param boundsX the horizontal bounds for the layout
+ * @param boundsY the vertical bounds for the layout
+ * @return a GROUP PShape containing two children: child 0 represents the edges
+ * as linework (LINES), and child 1 represents the vertices as points
+ * (POINTS)
* @since 1.3.0
*/
public static PShape fromGraph(SimpleGraph graph, double normalizationFactor, double boundsX, double boundsY) {
@@ -1059,32 +1065,43 @@ public static SimpleGraph toCentroidDualGraph(PShape mesh) {
*/
static SimpleGraph toDualGraph(Collection meshFaces) {
final SimpleGraph graph = new SimpleGraph<>(DefaultEdge.class);
- // map of which edge belong to each face; used to detect half-edges
- final HashMap edgesMap = new HashMap<>(meshFaces.size() * 4);
+ final Map> edgeFacesMap = new HashMap<>();
+ // Phase 1: Collect edges and their associated faces
for (PShape face : meshFaces) {
- graph.addVertex(face); // always add child so disconnected shapes are colored
+ graph.addVertex(face);
for (int i = 0; i < face.getVertexCount(); i++) {
- final PVector a = face.getVertex(i);
- final PVector b = face.getVertex((i + 1) % face.getVertexCount());
+ PVector a = face.getVertex(i);
+ PVector b = face.getVertex((i + 1) % face.getVertexCount());
if (a.equals(b)) {
continue;
}
- final PEdge e = new PEdge(a, b);
- final PShape neighbour = edgesMap.get(e);
- if (neighbour != null) {
- // edge seen before, so faces must be adjacent; create edge between faces
- if (neighbour.equals(face)) { // probably bad input (3 edges the same)
- System.err.println("toDualGraph(): Bad input — saw the same edge 3 times.");
- continue; // continue to prevent self-loop in graph
- }
- graph.addEdge(neighbour, face);
- } else {
- edgesMap.put(e, face); // edge is new
- }
+ PEdge edge = new PEdge(a, b);
+ edgeFacesMap.computeIfAbsent(edge, k -> new ArrayList<>()).add(face);
}
}
+
+ // Phase 2: Process edges in sorted order for graph iteration consistency
+ edgeFacesMap.entrySet().stream().sorted(Comparator.comparing(e -> e.getKey())) // Sort edges to ensure deterministic processing
+ .forEach(entry -> {
+ List faces = entry.getValue();
+ if (faces.size() == 2) {
+ // If exactly two faces share this edge, connect them in the dual graph
+ PShape f1 = faces.get(0);
+ PShape f2 = faces.get(1);
+ if (!f1.equals(f2)) {
+ graph.addEdge(f1, f2); // Avoid self-loops
+ } else {
+ // Handle case where the same face is associated with the edge twice
+ System.err.println("toDualGraph(): Bad input — saw the same edge 3+ times for face: " + f1);
+ }
+ } else if (faces.size() > 2) {
+ // Handle edges shared by more than two faces
+ System.err.println("toDualGraph(): Bad input — edge shared by more than two faces: " + entry.getKey().toString());
+ }
+ });
+
return graph;
}
@@ -1449,7 +1466,7 @@ public static PShape fromContours(List shell, @Nullable List points = toPVector(shape); // CLOSE
+ List points = toPVector(shape); // unclosed
if (shape.isClosed() && keepClosed) {
points.add(points.get(0).copy()); // since toPVector returns unclosed view
}
@@ -1486,7 +1503,8 @@ public static PShape fromArray(double[][] shape, boolean close) {
/**
* Flattens a collection of PShapes into a single GROUP PShape which has the
- * input shapes as its children.
+ * input shapes as its children. If the collection contains only one shape, it
+ * is directly returned.
*
* @since 1.2.0
* @see #flatten(PShape...)
@@ -1494,6 +1512,9 @@ public static PShape fromArray(double[][] shape, boolean close) {
public static PShape flatten(Collection shapes) {
PShape group = new PShape(GROUP);
shapes.forEach(group::addChild);
+ if (group.getChildCount() == 1) {
+ return group.getChild(0);
+ }
return group;
}
@@ -1599,6 +1620,7 @@ public static PShape setAllFillColor(PShape shape, int color) {
*
* @param shape
* @return the input object (having now been mutated)
+ * @see #setAllStrokeColor(PShape, int, double, int)
* @see #setAllFillColor(PShape, int)
*/
public static PShape setAllStrokeColor(PShape shape, int color, double strokeWeight) {
@@ -1610,6 +1632,27 @@ public static PShape setAllStrokeColor(PShape shape, int color, double strokeWei
return shape;
}
+ /**
+ * Sets the stroke color and cap style for the PShape and all of its children
+ * recursively.
+ *
+ * @param strokeCap either SQUARE, PROJECT, or
+ * ROUND
+ * @return the input object (having now been mutated)
+ * @see #setAllStrokeColor(PShape, int, double)
+ * @see #setAllFillColor(PShape, int)
+ * @since 2.1
+ */
+ public static PShape setAllStrokeColor(PShape shape, int color, double strokeWeight, int strokeCap) {
+ getChildren(shape).forEach(child -> {
+ child.setStroke(true);
+ child.setStroke(color);
+ child.setStrokeWeight((float) strokeWeight);
+ child.setStrokeCap(strokeCap);
+ });
+ return shape;
+ }
+
/**
* Sets the stroke color equal to the fill color for the PShape and all of its
* descendent shapes individually (that is, each child shape belonging to the
@@ -1707,21 +1750,40 @@ public static PShape disableAllStroke(PShape shape) {
/**
* Rounds the x and y coordinates (to the closest int) of all vertices belonging
- * to the shape, mutating the shape. This can sometimes fix a visual
- * problem in Processing where narrow gaps can appear between otherwise flush
- * shapes.
+ * to the shape. This can sometimes fix a visual problem in Processing where
+ * narrow gaps can appear between otherwise flush shapes. If the shape is a
+ * GROUP, the rounding is applied to all child shapes.
*
- * @return the input object (having now been mutated)
+ * @param shape the PShape to round vertex coordinates for.
+ * @return a rounded copy of the input shape.
+ * @see #roundVertexCoords(PShape, int)
* @since 1.1.3
*/
public static PShape roundVertexCoords(PShape shape) {
- getChildren(shape).forEach(c -> {
+ return roundVertexCoords(shape, 0);
+ }
+
+ /**
+ * Rounds the x and y coordinates (to n decimal places) of all
+ * vertices belonging to the shape. This can sometimes fix a visual problem in
+ * Processing where narrow gaps can appear between otherwise flush shapes. If
+ * the shape is a GROUP, the rounding is applied to all child shapes.
+ *
+ * @param shape the PShape to round vertex coordinates for.
+ * @param n The number of decimal places to which the coordinates should be
+ * rounded.
+ * @return a rounded copy of the input shape.
+ * @since 2.1
+ */
+ public static PShape roundVertexCoords(PShape shape, int n) {
+ return PGS_Processing.transform(shape, s -> {
+ var c = copy(s);
for (int i = 0; i < c.getVertexCount(); i++) {
final PVector v = c.getVertex(i);
- c.setVertex(i, Math.round(v.x), Math.round(v.y));
+ c.setVertex(i, round(v.x, n), round(v.y, n));
}
+ return c;
});
- return shape;
}
/**
@@ -1943,6 +2005,12 @@ private static T[] reversedCopy(T[] original) {
return reversed;
}
+ private static float round(float x, float n) {
+ float m = (float) FastMath.pow(10, n);
+
+ return FastMath.floor(m * x + 0.5f) / m;
+ }
+
/**
* A utility class for storing and manipulating the visual properties of PShapes
* from the Processing library. It encapsulates the stroke, fill, stroke color,
@@ -1990,8 +2058,9 @@ public static class PShapeData {
* Apply this shapedata to a given PShape.
*
* @param other
+ * @return other (fluent interface)
*/
- public void applyTo(PShape other) {
+ public PShape applyTo(PShape other) {
if (other.getFamily() == GROUP) {
getChildren(other).forEach(c -> applyTo(c));
}
@@ -2000,6 +2069,8 @@ public void applyTo(PShape other) {
other.setStroke(stroke);
other.setStroke(strokeColor);
other.setStrokeWeight(strokeWeight);
+
+ return other;
}
@Override
diff --git a/src/main/java/micycle/pgs/PGS_Meshing.java b/src/main/java/micycle/pgs/PGS_Meshing.java
index a4d992a1..c963a737 100644
--- a/src/main/java/micycle/pgs/PGS_Meshing.java
+++ b/src/main/java/micycle/pgs/PGS_Meshing.java
@@ -511,6 +511,9 @@ public static PShape centroidQuadrangulation(final IIncrementalTin triangulation
*/
private static PShape removeHoles(PShape faces, IIncrementalTin triangulation) {
List holes = new ArrayList<>(triangulation.getConstraints()); // copy list
+ if (holes.size() <= 1) {
+ return faces;
+ }
holes = holes.subList(1, holes.size()); // slice off perimeter constraint (not a hole)
STRtree tree = new STRtree();
diff --git a/src/main/java/micycle/pgs/PGS_Morphology.java b/src/main/java/micycle/pgs/PGS_Morphology.java
index d933360a..d29435ea 100644
--- a/src/main/java/micycle/pgs/PGS_Morphology.java
+++ b/src/main/java/micycle/pgs/PGS_Morphology.java
@@ -37,10 +37,12 @@
import micycle.pgs.color.Colors;
import micycle.pgs.commons.ChaikinCut;
import micycle.pgs.commons.CornerRounding;
+import micycle.pgs.commons.CornerRounding.RoundingStyle;
import micycle.pgs.commons.DiscreteCurveEvolution;
import micycle.pgs.commons.DiscreteCurveEvolution.DCETerminationCallback;
import micycle.pgs.commons.EllipticFourierDesc;
import micycle.pgs.commons.GaussianLineSmoothing;
+import micycle.pgs.commons.LaneRiesenfeldSmoothing;
import micycle.pgs.commons.ShapeInterpolation;
import micycle.uniformnoise.UniformNoise;
import processing.core.PConstants;
@@ -544,21 +546,60 @@ public static PShape smoothEllipticFourier(PShape shape, int descriptors) {
}
}
+ /**
+ * Smooths a shape using Lane-Riesenfeld curve subdivision with 4-point
+ * refinement to reduce contraction.
+ *
+ * @param shape A shape having lineal geometries (polygons or
+ * linestrings). Can be a GROUP shape consiting of
+ * these.
+ * @param degree The degree of the LR algorithm. Higher degrees
+ * influence the placement of vertices and the
+ * overall shape of the curve, but only slightly
+ * increase the number of vertices generated.
+ * Increasing the degree also increases the
+ * contraction of the curve toward its control
+ * points. The degree does not directly control the
+ * smoothness of the curve. A value of 3 or 4 is
+ * usually sufficient for most applications.
+ * @param subdivisions The number of times the subdivision process is
+ * applied. More subdivisions result in finer
+ * refinement and visually smoother curves between
+ * vertices. A value of 3 or 4 is usually
+ * sufficient for most applications.
+ * @param antiContractionFactor The weight parameter for the 4-point refinement.
+ * Controls the interpolation strength. A value of
+ * 0 effectively disables the contraction
+ * reduction. Generally suitable values are in
+ * [0...0.1]. Larger values may create
+ * self-intersecting geometry.
+ * @return A Shape having same structure as the input, whose geometries are now
+ * smooth.
+ * @since 2.1
+ */
+ public static PShape smoothLaneRiesenfeld(PShape shape, int degree, int subdivisions, double antiContractionFactor) {
+ return PGS.applyToLinealGeometries(shape, lineal -> LaneRiesenfeldSmoothing.subdivide(lineal, degree, subdivisions, antiContractionFactor));
+ }
+
/**
* Modifies the corners of a specified shape by replacing each angular corner
- * with a smooth, circular arc. The radius of each arc is determined
- * proportionally to the shorter of the two lines forming the corner.
+ * with a smooth, circular arc.
*
- * @param shape The original PShape object whose corners are to be rounded.
- * @param extent Specifies the degree of corner rounding, with a range from 0 to
- * 1. A value of 0 corresponds to no rounding, whereas a value of
- * 1 yields maximum rounding while still maintaining the validity
- * of the shape. Values above 1 are accepted but may produce
- * unpredictable results.
+ * @param shape A polygonal PShape, or GROUP shape having polygonal children.
+ * @param radius The radius of the circular arc used to round each corner. This
+ * determines how much a circle of the given radius "cuts into"
+ * the corner. The effective radius is bounded by the lengths of
+ * the edges forming the corner: If the radius is larger than half
+ * the length of either edge, it is clamped to the smaller of the
+ * two half-lengths to prevent overlapping or invalid geometry.
* @return A new PShape object with corners rounded to the specified extent.
*/
- public static PShape round(PShape shape, double extent) {
- return CornerRounding.round(shape, extent);
+ public static PShape round(PShape shape, double radius) {
+ return PGS_Processing.transform(shape, s -> {
+ var styling = PGS_Conversion.getShapeStylingData(shape);
+ var t = CornerRounding.roundCorners(s, radius, RoundingStyle.CIRCLE);
+ return styling.applyTo(t);
+ });
}
/**
diff --git a/src/main/java/micycle/pgs/PGS_Optimisation.java b/src/main/java/micycle/pgs/PGS_Optimisation.java
index bd3f229b..8b11aac8 100644
--- a/src/main/java/micycle/pgs/PGS_Optimisation.java
+++ b/src/main/java/micycle/pgs/PGS_Optimisation.java
@@ -39,6 +39,7 @@
import micycle.pgs.commons.LargestEmptyCircles;
import micycle.pgs.commons.MaximumInscribedAARectangle;
import micycle.pgs.commons.MaximumInscribedRectangle;
+import micycle.pgs.commons.MaximumInscribedTriangle;
import micycle.pgs.commons.MinimumBoundingEllipse;
import micycle.pgs.commons.MinimumBoundingTriangle;
import micycle.pgs.commons.Nullable;
@@ -107,7 +108,7 @@ public static PShape maximumInscribedCircle(PShape shape, PVector centerPoint) {
* Finds an approximate largest area rectangle (of arbitrary orientation)
* contained within a polygon.
*
- * @param shape
+ * @param shape a polygonal shape
* @return a rectangle shape
* @see #maximumInscribedAARectangle(PShape, boolean)
* maximumInscribedAARectangle() - the largest axis-aligned rectangle
@@ -117,6 +118,20 @@ public static PShape maximumInscribedRectangle(PShape shape) {
MaximumInscribedRectangle mir = new MaximumInscribedRectangle(polygon);
return toPShape(mir.computeMIR());
}
+
+ /**
+ * Finds an approximate largest area triangle (of arbitrary orientation)
+ * contained within a polygon.
+ *
+ * @param shape a polygonal shape
+ * @return a triangular shape
+ * @since 2.1
+ */
+ public static PShape maximumInscribedTriangle(PShape shape) {
+ Polygon polygon = (Polygon) fromPShape(shape);
+ var mit = new MaximumInscribedTriangle(polygon);
+ return toPShape(mit.computeMIT());
+ }
/**
* Finds the rectangle with a maximum area whose sides are parallel to the
@@ -567,7 +582,8 @@ public static PShape binPack(List shapes, double binWidth, double binHei
*
* @param shape
* @param point
- * @return
+ * @return closest point in the shape to the query (not a reference, but having
+ * the same coordinates)
* @see #closestPoints(PShape, PVector)
*/
public static PVector closestPoint(PShape shape, PVector point) {
@@ -576,6 +592,35 @@ public static PVector closestPoint(PShape shape, PVector point) {
return new PVector((float) coord.x, (float) coord.y);
}
+ /**
+ * Finds the closest point in the collection to a specified point.
+ *
+ * @param points the collection of points to search within
+ * @param point the point to find the closest neighbor for
+ * @return the closest point from the collection to the specified point
+ * @since 2.1
+ */
+ public static PVector closestPoint(Collection points, PVector point) {
+ if (points == null || points.isEmpty()) {
+ return null; // Handle empty or null collection
+ }
+
+ PVector closest = null;
+ float minDistanceSq = Float.MAX_VALUE;
+
+ for (PVector p : points) {
+ float dx = p.x - point.x;
+ float dy = p.y - point.y;
+ float distanceSq = dx * dx + dy * dy;
+ if (distanceSq < minDistanceSq) {
+ minDistanceSq = distanceSq;
+ closest = p;
+ }
+ }
+
+ return closest;
+ }
+
/**
* Returns the nearest point for each "island" / separate polygon in the GROUP
* input shape.
diff --git a/src/main/java/micycle/pgs/PGS_PointSet.java b/src/main/java/micycle/pgs/PGS_PointSet.java
index f5a4d869..76b48372 100644
--- a/src/main/java/micycle/pgs/PGS_PointSet.java
+++ b/src/main/java/micycle/pgs/PGS_PointSet.java
@@ -239,8 +239,8 @@ public static List random(double xMin, double yMin, double xMax, double
final SplittableRandom random = new SplittableRandom(seed);
final List points = new ArrayList<>(n);
for (int i = 0; i < n; i++) {
- final float x = (float) (xMin + (xMax - xMin) * random.nextDouble());
- final float y = (float) (yMin + (yMax - yMin) * random.nextDouble());
+ final float x = (float) random.nextDouble(xMin, xMax);
+ final float y = (float) random.nextDouble(yMin, yMax);
points.add(new PVector(x, y));
}
return points;
diff --git a/src/main/java/micycle/pgs/PGS_Processing.java b/src/main/java/micycle/pgs/PGS_Processing.java
index 23021ab2..638899cd 100644
--- a/src/main/java/micycle/pgs/PGS_Processing.java
+++ b/src/main/java/micycle/pgs/PGS_Processing.java
@@ -19,6 +19,7 @@
import java.util.function.BiConsumer;
import java.util.function.BiFunction;
import java.util.function.Consumer;
+import java.util.function.Function;
import java.util.function.Predicate;
import java.util.function.UnaryOperator;
import java.util.stream.Collectors;
@@ -1335,6 +1336,86 @@ public static PShape transformWithIndex(PShape shape, BiFunction function.apply(i, children.get(i))).filter(Objects::nonNull).toList());
}
+ /**
+ * Applies a specified transformation function to each child of the given PShape
+ * and returns a list of results produced by the function.
+ *
+ * This method processes each child of the input shape using the provided
+ * function, which can transform the shape into any desired type T. The function
+ * can return null to exclude a shape from the result list.
+ *
+ * Unlike the {@link #transform(PShape, UnaryOperator)} method, this method does
+ * not flatten the results into a PShape. Instead, it returns a list of
+ * arbitrary objects (type T) produced by the transformation function. This
+ * makes it more flexible for use cases where the transformation does not
+ * necessarily produce PShape objects.
+ *
+ * The transformation function can:
+ *
+ * - Transform the shape into a new object of type T
+ * - Return null to exclude the shape from the result list
+ *
+ *
+ * Note: This method does not modify the original shape or its children. It only
+ * applies the transformation function to each child and collects the results.
+ *
+ * @param The type of the objects produced by the transformation
+ * function.
+ * @param shape The PShape whose children will be transformed.
+ * @param function A Function that takes a PShape as input and returns an object
+ * of type T. If the function returns null for a shape, that
+ * shape will be excluded from the result list.
+ * @return A list of objects of type T produced by applying the transformation
+ * function to each child of the input shape.
+ * @see #transform(PShape, UnaryOperator)
+ * @since 2.1
+ */
+ public static List forEachShape(PShape shape, Function function) {
+ return PGS_Conversion.getChildren(shape).stream().map(function).filter(Objects::nonNull).toList();
+ }
+
+ /**
+ * Applies a specified transformation function to each child of the given PShape
+ * along with its index and returns a list of results produced by the function.
+ *
+ * This method processes each child of the input shape using the provided
+ * function, which takes both the index of the child and the child itself as
+ * input. The function can transform the shape into any desired type T or return
+ * null to exclude the shape from the result list.
+ *
+ * Unlike the {@link #transformWithIndex(PShape, BiFunction)} method, this
+ * method does not flatten the results into a PShape. Instead, it returns a list
+ * of arbitrary objects (type T) produced by the transformation function. This
+ * makes it more flexible for use cases where the transformation does not
+ * necessarily produce PShape objects.
+ *
+ * The transformation function can:
+ *
+ * - Transform the shape into a new object of type T
+ * - Return null to exclude the shape from the result list
+ *
+ *
+ * Note: This method does not modify the original shape or its children. It only
+ * applies the transformation function to each child and collects the results.
+ *
+ * @param The type of the objects produced by the transformation
+ * function.
+ * @param shape The PShape whose children will be transformed.
+ * @param function A BiFunction that takes an integer index and a PShape as
+ * input and returns an object of type T. If the function
+ * returns null for a shape, that shape will be excluded from
+ * the result list.
+ * @return A list of objects of type T produced by applying the transformation
+ * function to each child of the input shape along with its index.
+ * @see #transformWithIndex(PShape, BiFunction)
+ * @see #forEachShape(PShape, Function)
+ * @since 2.1
+ */
+ public static List forEachShapeWithIndex(PShape shape, BiFunction function) {
+ List children = PGS_Conversion.getChildren(shape);
+ return IntStream.range(0, children.size()).mapToObj(i -> function.apply(i, children.get(i))).filter(Objects::nonNull).toList();
+ }
+
/**
* Applies a specified function to each child of the given PShape.
*
diff --git a/src/main/java/micycle/pgs/PGS_ShapeBoolean.java b/src/main/java/micycle/pgs/PGS_ShapeBoolean.java
index 54a393be..efdbbae0 100644
--- a/src/main/java/micycle/pgs/PGS_ShapeBoolean.java
+++ b/src/main/java/micycle/pgs/PGS_ShapeBoolean.java
@@ -10,6 +10,7 @@
import java.util.Map;
import java.util.stream.Collectors;
+import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.prep.PreparedGeometry;
import org.locationtech.jts.geom.prep.PreparedGeometryFactory;
@@ -85,21 +86,34 @@ public static PShape intersect(final PShape a, final PShape b) {
* @since 1.3.0
*/
public static PShape intersectMesh(final PShape mesh, final PShape area) {
- final Geometry g = fromPShape(area);
- final PreparedGeometry cache = PreparedGeometryFactory.prepare(g);
+ final Geometry areaGeometry = fromPShape(area);
+ final PreparedGeometry preparedArea = PreparedGeometryFactory.prepare(areaGeometry);
+ final Envelope areaEnvelope = areaGeometry.getEnvelopeInternal();
+
+ return PGS_Conversion.toPShape(PGS_Conversion.getChildren(mesh).parallelStream().filter(s -> {
+ // Quick envelope check before more expensive operations
+ Geometry face = PGS_Conversion.fromPShape(s);
+ return areaEnvelope.intersects(face.getEnvelopeInternal());
+ }).map(s -> {
+ final Geometry face = PGS_Conversion.fromPShape(s);
+
+ // First try contains test as it's faster than intersection
+ if (preparedArea.containsProperly(face)) {
+ return face;
+ }
- List faces = PGS_Conversion.getChildren(mesh).parallelStream().map(s -> {
- final Geometry f = PGS_Conversion.fromPShape(s);
- if (cache.containsProperly(f)) {
- return f;
- } else {
- // preserve the fill etc of the PShape during intersection
- Geometry boundaryIntersect = OverlayNG.overlay(f, g, OverlayNG.INTERSECTION);
- boundaryIntersect.setUserData(f.getUserData());
- return boundaryIntersect;
+ // Only perform intersection if faces actually overlap
+ if (preparedArea.intersects(face)) {
+ Geometry intersection = OverlayNG.overlay(face, areaGeometry, OverlayNG.INTERSECTION);
+ // Only return non-empty intersections
+ if (!intersection.isEmpty()) {
+ intersection.setUserData(face.getUserData());
+ return intersection;
+ }
}
- }).collect(Collectors.toList());
- return PGS_Conversion.toPShape(faces);
+
+ return null;
+ }).filter(g -> g != null).collect(Collectors.toList()));
}
/**
@@ -380,7 +394,9 @@ public static PShape complement(PShape shape, double width, double height) {
shapeFactory.setNumPoints(4);
shapeFactory.setWidth(width);
shapeFactory.setHeight(height);
- return toPShape(shapeFactory.createRectangle().difference(fromPShape(shape)));
+ // unioning difference shape helps robustness (when it comprises overlapping
+ // children)
+ return toPShape(shapeFactory.createRectangle().difference(fromPShape(shape).union()));
}
}
diff --git a/src/main/java/micycle/pgs/PGS_ShapePredicates.java b/src/main/java/micycle/pgs/PGS_ShapePredicates.java
index 21186de8..39383651 100644
--- a/src/main/java/micycle/pgs/PGS_ShapePredicates.java
+++ b/src/main/java/micycle/pgs/PGS_ShapePredicates.java
@@ -5,7 +5,9 @@
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
+import java.util.HashMap;
import java.util.List;
+import java.util.Map;
import javax.vecmath.Point3d;
import javax.vecmath.Point4d;
@@ -415,11 +417,14 @@ public static double sphericity(final PShape shape) {
}
/**
- * Measures the elongation of a shape; the ratio of a shape's bounding box
- * length to its width.
+ * Measures the elongation of a shape as the ratio of the difference between the
+ * bounding box's length and width to the maximum dimension. A value of 1
+ * indicates a highly elongated shape, while a value of 0 indicates a square or
+ * nearly square shape.
*
* @param shape
- * @return a value in [0, 1]
+ * @return a value in the range [0, 1], where 1 represents high elongation and 0
+ * represents no elongation
*/
public static double elongation(final PShape shape) {
Geometry obb = MinimumDiameter.getMinimumRectangle(fromPShape(shape));
@@ -429,11 +434,9 @@ public static double elongation(final PShape shape) {
Coordinate c2 = rect.getCoordinates()[2];
double l = c0.distance(c1);
double w = c1.distance(c2);
- if (l >= w) {
- return w / l;
- } else {
- return l / w;
- }
+ double max = Math.max(l, w);
+ double min = Math.min(l, w);
+ return 1 - (min / max);
}
/**
@@ -510,6 +513,45 @@ public static double maximumInteriorAngle(PShape shape) {
return maxAngle;
}
+ /**
+ * Calculates the interior angles of a polygon represented by a {@link PShape}.
+ * The method calculates the interior angle at each vertex.
+ *
+ * The vertices of the input {@code shape} are assumed to represent a simple
+ * polygon.
+ *
+ * @param shape The {@link PShape} representing the polygon for which to
+ * calculate interior angles. It's expected to be a polygon shape.
+ * @return A map where keys are {@link PVector} vertices of the polygon and
+ * values are their corresponding interior angles in radians as double
+ * values.
+ * @since 2.1
+ */
+ public static Map interiorAngles(PShape shape) {
+ Map anglesMap = new HashMap<>();
+
+ var vertices = PGS_Conversion.toPVector(shape); // unclosed
+ if (!PGS.isClockwise(vertices)) {
+ Collections.reverse(vertices);
+ }
+
+ int n = vertices.size();
+ for (int i = 0; i < n; i++) {
+ PVector currentVertex = vertices.get(i);
+ PVector previousVertex = vertices.get((i - 1 + n) % n); // Get previous vertex, wrapping around
+ PVector nextVertex = vertices.get((i + 1) % n); // Get next vertex, wrapping around
+
+ Coordinate p0 = new Coordinate(previousVertex.x, previousVertex.y);
+ Coordinate p1 = new Coordinate(currentVertex.x, currentVertex.y);
+ Coordinate p2 = new Coordinate(nextVertex.x, nextVertex.y);
+ double angleRadians = Angle.interiorAngle(p0, p1, p2); // CW
+ anglesMap.put(currentVertex, angleRadians);
+ }
+
+ return anglesMap;
+
+ }
+
/**
* Quantifies the similarity between two shapes, by using the pairwise euclidean
* distance between each shape's Elliptic Fourier Descriptors (EFD).
diff --git a/src/main/java/micycle/pgs/PGS_Voronoi.java b/src/main/java/micycle/pgs/PGS_Voronoi.java
index 6415904a..a5464b89 100644
--- a/src/main/java/micycle/pgs/PGS_Voronoi.java
+++ b/src/main/java/micycle/pgs/PGS_Voronoi.java
@@ -533,7 +533,6 @@ public static PShape multiplicativelyWeightedVoronoi(Collection sites,
* @since 2.0
*/
public static PShape multiplicativelyWeightedVoronoi(Collection sites, double[] bounds, boolean forceConforming) {
-
var faces = MultiplicativelyWeightedVoronoi.getMWVFromPVectors(sites.stream().toList(), bounds);
Geometry geoms = PGS.GEOM_FACTORY.createGeometryCollection(faces.toArray(new Geometry[] {}));
if (forceConforming) {
diff --git a/src/main/java/micycle/pgs/commons/CornerRounding.java b/src/main/java/micycle/pgs/commons/CornerRounding.java
index f4af4ce5..56e14152 100644
--- a/src/main/java/micycle/pgs/commons/CornerRounding.java
+++ b/src/main/java/micycle/pgs/commons/CornerRounding.java
@@ -1,157 +1,215 @@
package micycle.pgs.commons;
+import static java.lang.Math.PI;
+
+import java.util.ArrayList;
import java.util.List;
-import micycle.pgs.PGS_Conversion;
import micycle.pgs.color.Colors;
-import processing.core.PApplet;
+import net.jafama.FastMath;
import processing.core.PConstants;
import processing.core.PShape;
import processing.core.PVector;
/**
- * Corner rounding for PShape polygons. Replaces corners with arcs.
- *
- * @author Michael Carleton
+ * A utility class for rounding the corners of a polygon represented as a
+ * {@link PShape}.
+ *
+ * The implementation is based on the algorithm described in the ObservableHQ
+ * notebook: Rounding
+ * Polygon Corners. It calculates the necessary arc points for each corner
+ * and constructs a new {@link PShape} with rounded corners.
+ *
+ * The class supports different rounding styles, such as strictly geometric
+ * rounding, natural rounding, and freehand-style rounding, through the
+ * {@link RoundingStyle} enum.
*
+ * @author Mathieu Jouhet
+ * @author Michael Carleton
*/
-public final class CornerRounding {
+public class CornerRounding {
- // Inspired by https://observablehq.com/@daformat/rounding-polygon-corners
-
- private CornerRounding() {
+ /**
+ * An enum representing the rounding style for the corners of the polygon.
+ *
+ * - {@link #CIRCLE}: Strictly geometric rounding, using perfect circular
+ * arcs.
+ * - {@link #APPROX}: Natural rounding, approximating circular arcs with
+ * Bézier curves.
+ * - {@link #HAND}: Freehand-style rounding, providing a more organic
+ * look.
+ *
+ */
+ public enum RoundingStyle {
+ CIRCLE, // Strictly geometric rounding
+ APPROX, // Natural rounding
+ HAND // Freehand-style rounding
}
/**
- *
- * @param shape
- * @param extent 0...1
- * @return
+ * Rounds the corners of a given {@link PShape} using the specified radius and
+ * rounding style. The method generates circular arcs for each corner and
+ * constructs a new {@link PShape} with the rounded corners.
+ *
+ * @param original The original {@link PShape} whose corners are to be rounded.
+ * @param radius The radius of the circular arc used to round each corner.
+ * This determines how much a circle of the given radius "cuts
+ * into" the corner. The effective radius is bounded by the
+ * lengths of the edges forming the corner: If the radius is
+ * larger than half the length of either edge, it is clamped to
+ * the smaller of the two half-lengths to prevent overlapping or
+ * invalid geometry.
+ * @param style The rounding style to apply, as defined in the
+ * {@link RoundingStyle} enum.
+ * @return A new {@link PShape} with rounded corners. If the input shape is null
+ * or the radius is zero, the original shape is returned unchanged.
*/
- public static PShape round(PShape shape, double extent) {
- if (shape.getChildCount() > 1) {
- PShape groupRound = new PShape(PConstants.GROUP);
- for (PShape child : shape.getChildren()) {
- groupRound.addChild(roundPolygon(child, extent));
+ public static PShape roundCorners(PShape original, double radius, RoundingStyle style) {
+ if (original == null || radius == 0) {
+ return original;
+ }
+
+ List vertices = extractAndFilterVertices(original);
+
+ int numVertices = vertices.size();
+ if (numVertices == 0) {
+ return original;
+ }
+
+ List corners = new ArrayList<>();
+
+ for (int i = 0; i < numVertices; i++) {
+ PVector c1 = vertices.get(i);
+ PVector c2 = vertices.get((i + 1) % numVertices);
+ PVector c3 = vertices.get((i + 2) % numVertices);
+
+ PVector vC1c2 = PVector.sub(c1, c2);
+ PVector vC3c2 = PVector.sub(c3, c2);
+
+ float dx1 = vC1c2.x;
+ float dy1 = vC1c2.y;
+ float mag1 = vC1c2.mag();
+ if (mag1 == 0) {
+ continue;
+ }
+ float unitX1 = dx1 / mag1;
+ float unitY1 = dy1 / mag1;
+
+ float dx3 = vC3c2.x;
+ float dy3 = vC3c2.y;
+ float mag3 = vC3c2.mag();
+ if (mag3 == 0) {
+ continue;
+ }
+ float unitX3 = dx3 / mag3;
+ float unitY3 = dy3 / mag3;
+
+ float cross = dx1 * dy3 - dy1 * dx3;
+ float dot = dx1 * dx3 + dy1 * dy3;
+ final float angle = abs(atan2(cross, dot)); // == angleBetween(vC1c2, vC3c2)
+
+ float cornerLength = min((float) radius, mag1 / 2, mag3 / 2);
+ if (cornerLength <= 0) {
+ continue;
}
- return groupRound;
+
+ float a = cornerLength * tan(angle / 2);
+
+ float idealCPDistance = 0;
+ idealCPDistance = switch (style) {
+ case CIRCLE -> {
+ double numPointsCircle = (2 * PI) / (PI - angle);
+ yield (4.0f / 3) * tan(PI / (2 * numPointsCircle)) * a;
+ }
+ case APPROX -> {
+ float multiplier = angle < PI / 2 ? 1 + cos(angle) : 2 - sin(angle);
+ yield (4.0f / 3) * tan(PI / (2 * ((2 * PI) / angle))) * cornerLength * multiplier;
+ }
+ case HAND -> (4.0f / 3) * tan(PI / (2 * ((2 * PI) / angle))) * cornerLength * (2 + sin(angle));
+ default -> (4.0f / 3) * cornerLength;
+ };
+
+ float cpDistance = cornerLength - idealCPDistance;
+
+ PVector c1CurvePoint = new PVector(c2.x + unitX1 * cornerLength, c2.y + unitY1 * cornerLength);
+ PVector c1CP = new PVector(c2.x + unitX1 * cpDistance, c2.y + unitY1 * cpDistance);
+ PVector c3CurvePoint = new PVector(c2.x + unitX3 * cornerLength, c2.y + unitY3 * cornerLength);
+ PVector c3CP = new PVector(c2.x + unitX3 * cpDistance, c2.y + unitY3 * cpDistance);
+
+ corners.add(new CornerData(c1CurvePoint, c1CP, c3CP, c3CurvePoint));
}
- else {
- return roundPolygon(shape, extent);
+
+ if (corners.isEmpty()) {
+ return original;
}
- }
-
- private static PShape roundPolygon(PShape shape, double extent) {
- PShape rounded = new PShape(PShape.GEOMETRY);
- PGS_Conversion.setAllFillColor(rounded, Colors.PINK);
+
+ PVector startPoint = corners.get(corners.size() - 1).c3CurvePoint;
+ PShape rounded = new PShape(PShape.PATH);
+ rounded.setStroke(Colors.PINK);
+ rounded.setStroke(true);
+ rounded.setFill(true);
+ rounded.setFill(Colors.WHITE);
rounded.beginShape();
+ rounded.vertex(startPoint.x, startPoint.y);
- final List l = PGS_Conversion.toPVector(shape);
- final int size = l.size();
- for (int i = 0; i < l.size(); i++) {
- final PVector p1 = l.get(Math.floorMod(i - 1, size));
- final PVector p2 = l.get(i);
- final PVector p3 = l.get((i + 1) % size);
- roundCorner(p1, p2, p3, extent, rounded);
+ for (CornerData corner : corners) {
+ rounded.vertex(corner.c1CurvePoint.x, corner.c1CurvePoint.y);
+ rounded.bezierVertex(corner.c1CP.x, corner.c1CP.y, corner.c3CP.x, corner.c3CP.y, corner.c3CurvePoint.x, corner.c3CurvePoint.y);
}
+
rounded.endShape(PConstants.CLOSE);
return rounded;
}
- /**
- * Round a triplet of points.
- *
- * @param a
- * @param b middle/enclosed point
- * @param c
- * @param extent
- * @param shape
- */
- private static void roundCorner(PVector a, PVector b, PVector c, double extent, PShape shape) {
- if (clockwise(a, b, c)) {
- PVector temp = a;
- a = c;
- c = temp;
- }
-
- // line vectors
- PVector ab = PVector.sub(a, b);
- PVector cb = PVector.sub(c, b);
-
- float theta = PApplet.acos(ab.dot(cb) / (ab.mag() * cb.mag())); // same as a.angleBetween(a, c)
-
-// final float maxRadius = PApplet.min(ab.div(2).mag(), cb.div(2).mag());
-// extent = extent * maxRadius;
- final float extentF = (float) extent;
+ private static List extractAndFilterVertices(PShape shape) {
+ List vertices = new ArrayList<>();
+ int vertexCount = shape.getVertexCount();
+ if (vertexCount > 0) {
+ PVector firstVertex = new PVector(shape.getVertexX(0), shape.getVertexY(0));
+ vertices.add(firstVertex);
+ PVector prevVertex = firstVertex;
+
+ for (int i = 1; i < vertexCount; i++) {
+ PVector currentVertex = new PVector(shape.getVertexX(i), shape.getVertexY(i));
+ if (currentVertex.x != prevVertex.x || currentVertex.y != prevVertex.y) {
+ vertices.add(currentVertex);
+ prevVertex = currentVertex;
+ }
+ }
- final PVector A = PVector.add(b, ab.mult(extentF / ab.mag())); // where circle touches AB
- final PVector C = PVector.add(b, cb.mult(extentF / cb.mag())); // where circle touches CB
+ // Check if last vertex is the same as the first (for closed shapes)
+ if (vertexCount > 1 && firstVertex.x == prevVertex.x && firstVertex.y == prevVertex.y) {
+ vertices.remove(vertices.size() - 1);
+ }
+ }
+ return vertices;
+ }
- PVector vBC = PVector.sub(C, b);
+ private record CornerData(PVector c1CurvePoint, PVector c1CP, PVector c3CP, PVector c3CurvePoint) {
+ }
- final float lengthBC = vBC.mag();
- final float lengthBD = PApplet.cos(theta / 2); // length
- final float lengthFD = PApplet.sin(theta / 2) * lengthBD; // length
- final float lengthBF = lengthBD * lengthBD; // length
- final float r = lengthFD / (lengthBF / lengthBC); // radius of circle
- PVector tangent = ab.normalize().rotate(PConstants.HALF_PI).mult(r); // tangent to ab
- tangent.add(A); // find tangent to ab at point A -- this is circle center
+ private static float min(float a, float b, float c) {
+ return Math.min(a, Math.min(b, c));
+ }
- shape.vertex(C.x, C.y);
- final float a1 = angleBetween(C, tangent);
- final float a2 = angleBetween(A, tangent);
- sampleArc(tangent, r, a1, a2, shape);
- shape.vertex(A.x, A.y);
+ private static float abs(float val) {
+ return Math.abs(val);
}
- /**
- * Sub-sample an arc into invididual vertices.
- *
- * @param center
- * @param r arc radius
- * @param startAngle
- * @param endAngle
- * @param shape
- */
- private static void sampleArc(PVector center, float r, float startAngle, float endAngle, PShape shape) {
- if (startAngle > endAngle) {
- startAngle -= PConstants.TWO_PI;
- }
- final float n = 4; // every n degrees // TODO magic constant
- final float angleInc = (endAngle - startAngle) / (360 / n);
- float angle = startAngle;
- while (angle < endAngle) {
- shape.vertex(r * PApplet.cos(angle) + center.x, center.y + r * PApplet.sin(angle));
- angle += angleInc;
- }
+ private static float atan2(float y, float x) {
+ return (float) FastMath.atan2(y, x);
}
- /**
- * Are the three given points in clockwise order?
- *
- * @param p1
- * @param p2 middle point
- * @param p3
- * @return true if the points consitute a "right turn" or clockwise orientation
- */
- private static boolean clockwise(PVector p1, PVector p2, PVector p3) {
- float val = (p2.y - p1.y) * (p3.x - p2.x) - (p2.x - p1.x) * (p3.y - p2.y);
- return val < 0;
+ private static float tan(double angle) {
+ return (float) FastMath.tan(angle);
}
- /**
- * East = 0; North = -1/2PI; West = -PI; South = -3/2PI | 1/2PI
- *
- * @param tail PVector Coordinate 1.
- * @param head PVector Coordinate 2.
- * @return float θ in radians.
- */
- private static float angleBetween(PVector tail, PVector head) {
- float a = FastAtan2.atan2(tail.y - head.y, tail.x - head.x);
- if (a < 0) {
- a += PConstants.TWO_PI;
- }
- return a;
+ private static float cos(double angle) {
+ return (float) FastMath.cos(angle);
}
-}
+ private static float sin(double angle) {
+ return (float) FastMath.sin(angle);
+ }
+}
\ No newline at end of file
diff --git a/src/main/java/micycle/pgs/commons/LaneRiesenfeldSmoothing.java b/src/main/java/micycle/pgs/commons/LaneRiesenfeldSmoothing.java
new file mode 100644
index 00000000..55998831
--- /dev/null
+++ b/src/main/java/micycle/pgs/commons/LaneRiesenfeldSmoothing.java
@@ -0,0 +1,219 @@
+package micycle.pgs.commons;
+
+import org.locationtech.jts.geom.*;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+
+/**
+ * Geometry smoothing using Lane-Riesenfeld (LR) curve subdivision with 4-point
+ * refinement to reduce contraction.
+ *
+ * The LR algorithm is a generalization of the Chaikin subdivision which
+ * generates splines with variable continuity. The 4-point (Dyn-Levin-Gregory)
+ * refinement is a variant that interpolates the control points. A combination
+ * of LR and 4-point refinement can be used to reduce contraction when smoothing
+ * a geometry.
+ *
+ * This class provides a utility method to subdivide geometries using the LR
+ * algorithm with 4-point refinement. The algorithm can be applied to both open
+ * and closed geometries (e.g., LineString and LinearRing).
+ *
+ * @author Michael Carleton
+ */
+public class LaneRiesenfeldSmoothing {
+
+ // https://observablehq.com/@esperanc/lane-riesenfeld-subdivision
+ // https://tiborstanko.sk/teaching/geo-num-2017/tp5.html
+
+ /**
+ * Subdivides the input geometry using the Lane-Riesenfeld algorithm with
+ * 4-point refinement.
+ *
+ * @param geometry The lineal input geometry to subdivide.
+ * @param degree The degree of the LR algorithm. Higher degrees
+ * influence the placement of vertices and the
+ * overall shape of the curve, but only slightly
+ * increase the number of vertices generated.
+ * Increasing the degree also increases the
+ * contraction of the curve toward its control
+ * points. The degree does not directly control the
+ * smoothness of the curve. A value of 3 or 4 is
+ * usually sufficient for most applications.
+ * @param subdivisions The number of times the subdivision process is
+ * applied. More subdivisions result in finer
+ * refinement and visually smoother curves between
+ * vertices. A value of 3 or 4 is usually
+ * sufficient for most applications.
+ * @param antiContractionFactor The weight parameter for the 4-point refinement.
+ * Controls the interpolation strength. A value of
+ * 0 effectively disables the contraction
+ * reduction. Generally suitable values are in
+ * [0...0.1]. Larger values may create
+ * self-intersecting geometry.
+ * @return A new subdivided geometry (LineString or LinearRing).
+ */
+ public static LineString subdivide(LineString geometry, int degree, int subdivisions, double antiContractionFactor) {
+ Coordinate[] coords = geometry.getCoordinates();
+ boolean closed = geometry.isClosed();
+ if (closed) {
+ coords = Arrays.copyOf(coords, coords.length - 1); // Remove the last coordinate if closed
+ }
+ Coordinate[] subdivided = lr4(coords, degree, closed, antiContractionFactor, subdivisions);
+ GeometryFactory factory = geometry.getFactory();
+ return createGeometry(factory, subdivided, closed);
+ }
+
+ private static LineString createGeometry(GeometryFactory factory, Coordinate[] coords, boolean closed) {
+ if (closed) {
+ List coordList = new ArrayList<>(Arrays.asList(coords));
+ coordList.add(new Coordinate(coordList.get(0)));
+ return factory.createLinearRing(coordList.toArray(new Coordinate[0]));
+ } else {
+ return factory.createLineString(coords);
+ }
+ }
+
+ /**
+ * Applies the Lane-Riesenfeld algorithm with 4-point refinement to an array of
+ * coordinates. This replaces the vanilla midpoint averaging with four-point
+ * averaging.
+ *
+ * @param points The input array of coordinates.
+ * @param degree The degree of the Lane-Riesenfeld algorithm.
+ * @param closed Whether the geometry is closed.
+ * @param w The weight parameter for the 4-point refinement.
+ * @param subdivisions The number of times the subdivision process is applied.
+ * @return A new array of subdivided coordinates.
+ */
+ private static Coordinate[] lr4(Coordinate[] points, int degree, boolean closed, double w, int subdivisions) {
+ if (degree < 1 || subdivisions < 1) {
+ return Arrays.copyOf(points, points.length);
+ }
+
+ Coordinate[] v = points;
+
+ for (int s = 0; s < subdivisions; s++) {
+ v = fourPoint(v, closed, w);
+
+ for (int d = 1; d < degree; d++) {
+ int n = v.length;
+ List u = new ArrayList<>();
+
+ for (int i = 0; i < n; i++) {
+ int prevIndex = getPreviousIndex(i, n, closed);
+ int nextIndex = getNextIndex(i, n, closed);
+ int nextNextIndex = getNextNextIndex(i, n, closed);
+
+ Coordinate p0 = v[prevIndex];
+ Coordinate p1 = v[i];
+ Coordinate p2 = v[nextIndex];
+ Coordinate p3 = v[nextNextIndex];
+
+ double qx = computeQ(p0.x, p1.x, p2.x, p3.x, w);
+ double qy = computeQ(p0.y, p1.y, p2.y, p3.y, w);
+ u.add(new Coordinate(qx, qy));
+ }
+
+ if (closed) {
+ v = u.toArray(new Coordinate[0]);
+ } else {
+ List newV = new ArrayList<>();
+ if (v.length > 0) {
+ newV.add(v[0]);
+ }
+ if (!u.isEmpty()) {
+ newV.addAll(u.subList(0, u.size() - 1));
+ }
+ if (v.length > 0) {
+ newV.add(v[v.length - 1]);
+ }
+ v = newV.toArray(new Coordinate[0]);
+ }
+ }
+ }
+
+ return v;
+ }
+
+ /**
+ * Applies the 4-point (Dyn-Levin-Gregory) refinement to an array of
+ * coordinates.
+ *
+ * @param points The input array of coordinates.
+ * @param closed Whether the geometry is closed.
+ * @param w The weight parameter for the 4-point refinement.
+ * @return A new array of refined coordinates.
+ */
+ private static Coordinate[] fourPoint(Coordinate[] points, boolean closed, double w) {
+ int n = points.length;
+ if (n == 0) {
+ return new Coordinate[0];
+ }
+
+ List result = new ArrayList<>(points.length);
+
+ for (int i = 0; i < n; i++) {
+ int prevIndex = getPreviousIndex(i, n, closed);
+ int nextIndex = getNextIndex(i, n, closed);
+ int nextNextIndex = getNextNextIndex(i, n, closed);
+
+ Coordinate p0 = points[prevIndex];
+ Coordinate p1 = points[i];
+ Coordinate p2 = points[nextIndex];
+ Coordinate p3 = points[nextNextIndex];
+
+ double qx = computeQ(p0.x, p1.x, p2.x, p3.x, w);
+ double qy = computeQ(p0.y, p1.y, p2.y, p3.y, w);
+ Coordinate q = new Coordinate(qx, qy);
+
+ result.add(p1);
+ result.add(q);
+ }
+
+ if (!closed && !result.isEmpty()) {
+ result.remove(result.size() - 1);
+ }
+
+ return result.toArray(new Coordinate[0]);
+ }
+
+ /**
+ * Computes a new coordinate value using the 4-point refinement formula.
+ *
+ * @param p0 The first coordinate value.
+ * @param p1 The second coordinate value.
+ * @param p2 The third coordinate value.
+ * @param p3 The fourth coordinate value.
+ * @param w The weight parameter.
+ * @return The computed coordinate value.
+ */
+ private static double computeQ(double p0, double p1, double p2, double p3, double w) {
+ return -w * p0 + (0.5 + w) * p1 + (0.5 + w) * p2 - w * p3;
+ }
+
+ private static int getPreviousIndex(int i, int n, boolean closed) {
+ if (closed) {
+ return (i - 1 + n) % n;
+ } else {
+ return Math.max(0, i - 1);
+ }
+ }
+
+ private static int getNextIndex(int i, int n, boolean closed) {
+ if (closed) {
+ return (i + 1) % n;
+ } else {
+ return Math.min(n - 1, i + 1);
+ }
+ }
+
+ private static int getNextNextIndex(int i, int n, boolean closed) {
+ if (closed) {
+ return (i + 2) % n;
+ } else {
+ return Math.min(n - 1, i + 2);
+ }
+ }
+}
\ No newline at end of file
diff --git a/src/main/java/micycle/pgs/commons/MaximumInscribedRectangle.java b/src/main/java/micycle/pgs/commons/MaximumInscribedRectangle.java
index 952d0556..e330a4fd 100644
--- a/src/main/java/micycle/pgs/commons/MaximumInscribedRectangle.java
+++ b/src/main/java/micycle/pgs/commons/MaximumInscribedRectangle.java
@@ -4,8 +4,8 @@
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.geom.GeometryFactory;
import org.locationtech.jts.geom.prep.PreparedGeometry;
import org.locationtech.jts.geom.prep.PreparedGeometryFactory;
@@ -13,156 +13,239 @@
import net.metaopt.swarm.FitnessFunction;
import net.metaopt.swarm.pso.Particle;
import net.metaopt.swarm.pso.Swarm;
-import processing.core.PVector;
+
+import org.apache.commons.math3.analysis.MultivariateFunction;
+import org.apache.commons.math3.optim.InitialGuess;
+import org.apache.commons.math3.optim.MaxEval;
+import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
+import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
+import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction;
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex;
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer;
+import org.apache.commons.math3.optim.PointValuePair;
+
+import java.util.Random;
/**
* Finds an approximate largest area rectangle of arbitrary orientation in a
- * polygon via particle swarm optimisation.
+ * concave polygon via particle swarm optimisation.
*
* @author Michael Carleton
*
*/
public class MaximumInscribedRectangle {
- private static final int SWARM_SIZE = 2500; // Reduced from 2500
- private static final int GENERATIONS = 1000; // Reduced from 1000
- private static final double ASPECT_WEIGHT = 0.05; // Small bonus for non-square shapes
+ private static final int SWARM_SIZE = 2500; // swarm size is tunable
+ private static final int GENERATIONS = 1000; // maximum PSO generations
+ private static final double ASPECT_WEIGHT = 0.05; // bonus for non-square shapes
private static final GeometryFactory GEOM_FACTORY = new GeometryFactory();
+ private static final Random RANDOM = new Random();
+
private final Swarm swarm;
+ private final RectangleFitness fitnessFunction; // used by both PSO and local optimization
public MaximumInscribedRectangle(Polygon polygon) {
- final MaximumInscribedCircle mic = new MaximumInscribedCircle(polygon, 2);
- final double minSquare = mic.getRadiusLine().getLength() / Math.sqrt(2);
+ // We compute a lower bound based on the maximum inscribed circle.
+ MaximumInscribedCircle mic = new MaximumInscribedCircle(polygon, 2);
+ double minSquare = mic.getRadiusLine().getLength() / Math.sqrt(2);
- // Create multiple sub-swarms with different initial conditions
- final FitnessFunction fitnessFunction = new RectangleFitness(polygon);
- swarm = new Swarm(SWARM_SIZE, new RectangleCandidate(), fitnessFunction);
+ fitnessFunction = new RectangleFitness(polygon);
+ swarm = new ParallelSwarm(SWARM_SIZE, new RectangleCandidate(), fitnessFunction);
- final Envelope e = polygon.getEnvelopeInternal();
- final double[] maxPosition = new double[] { e.getMaxX(), e.getMaxY(), e.getWidth(), e.getHeight(), Math.PI };
- final double[] minPosition = new double[] { e.getMinX(), e.getMinY(), minSquare, minSquare, -Math.PI };
+ Envelope e = polygon.getEnvelopeInternal();
+ double[] maxPosition = new double[] { e.getMaxX(), e.getMaxY(), e.getWidth(), e.getHeight(), Math.PI };
+ double[] minPosition = new double[] { e.getMinX(), e.getMinY(), minSquare, minSquare, -Math.PI };
swarm.setMaxPosition(maxPosition);
swarm.setMinPosition(minPosition);
-
- // Initialize particles in promising regions
swarm.init();
initializeSwarm(swarm, e, minSquare);
}
private void initializeSwarm(Swarm swarm, Envelope e, double minSquare) {
Particle[] particles = swarm.getParticles();
- int particlesPerRegion = SWARM_SIZE / 4;
+ final int particlesPerRegion = SWARM_SIZE / 4; // Four groups
+ double minX = e.getMinX();
+ double minY = e.getMinY();
+ double width = e.getWidth();
+ double height = e.getHeight();
for (int i = 0; i < SWARM_SIZE; i++) {
double[] position = new double[5];
if (i < particlesPerRegion) {
- // Group 1: Initialize along envelope edges
- position[0] = e.getMinX() + Math.random() * e.getWidth();
- position[1] = Math.random() < 0.5 ? e.getMinY() : e.getMaxY();
- position[2] = e.getWidth() * (0.3 + Math.random() * 0.7);
- position[3] = e.getHeight() * (0.3 + Math.random() * 0.7);
- position[4] = Math.random() * Math.PI - Math.PI / 2;
+ // Group 1: Along envelope edges
+ position[0] = minX + RANDOM.nextDouble() * width;
+ position[1] = RANDOM.nextDouble() < 0.5 ? e.getMinY() : e.getMaxY();
+ position[2] = width * (0.3 + RANDOM.nextDouble() * 0.7);
+ position[3] = height * (0.3 + RANDOM.nextDouble() * 0.7);
+ position[4] = RANDOM.nextDouble() * Math.PI - Math.PI / 2;
} else if (i < 2 * particlesPerRegion) {
- // Group 2: Try vertical rectangles
- position[0] = e.getMinX() + Math.random() * e.getWidth();
- position[1] = e.getMinY() + Math.random() * e.getHeight();
- position[2] = e.getWidth() * (0.1 + Math.random() * 0.3);
- position[3] = e.getHeight() * (0.7 + Math.random() * 0.3);
- position[4] = Math.random() * Math.PI / 6 - Math.PI / 12;
+ // Group 2: Vertical rectangles
+ position[0] = minX + RANDOM.nextDouble() * width;
+ position[1] = minY + RANDOM.nextDouble() * height;
+ position[2] = width * (0.1 + RANDOM.nextDouble() * 0.3);
+ position[3] = height * (0.7 + RANDOM.nextDouble() * 0.3);
+ position[4] = RANDOM.nextDouble() * (Math.PI / 6) - Math.PI / 12;
} else if (i < 3 * particlesPerRegion) {
- // Group 3: Try horizontal rectangles
- position[0] = e.getMinX() + Math.random() * e.getWidth();
- position[1] = e.getMinY() + Math.random() * e.getHeight();
- position[2] = e.getWidth() * (0.7 + Math.random() * 0.3);
- position[3] = e.getHeight() * (0.1 + Math.random() * 0.3);
- position[4] = Math.PI / 2 + Math.random() * Math.PI / 6 - Math.PI / 12;
+ // Group 3: Horizontal rectangles
+ position[0] = minX + RANDOM.nextDouble() * width;
+ position[1] = minY + RANDOM.nextDouble() * height;
+ position[2] = width * (0.7 + RANDOM.nextDouble() * 0.3);
+ position[3] = height * (0.1 + RANDOM.nextDouble() * 0.3);
+ position[4] = Math.PI / 2 + RANDOM.nextDouble() * (Math.PI / 6) - Math.PI / 12;
} else {
// Group 4: Random exploration
- position[0] = e.getMinX() + Math.random() * e.getWidth();
- position[1] = e.getMinY() + Math.random() * e.getHeight();
- position[2] = minSquare + Math.random() * (e.getWidth() - minSquare);
- position[3] = minSquare + Math.random() * (e.getHeight() - minSquare);
- position[4] = Math.random() * Math.PI - Math.PI / 2;
+ position[0] = minX + RANDOM.nextDouble() * width;
+ position[1] = minY + RANDOM.nextDouble() * height;
+ position[2] = minSquare + RANDOM.nextDouble() * (width - minSquare);
+ position[3] = minSquare + RANDOM.nextDouble() * (height - minSquare);
+ position[4] = RANDOM.nextDouble() * Math.PI - Math.PI / 2;
}
-
particles[i].setPosition(position);
particles[i].setVelocity(new double[] { 0, 0, 0, 0, 0 });
particles[i].setBestPosition(position.clone());
}
}
+ /**
+ * This method runs the swarm evolution and then refines the best candidate
+ * using Apache Commons Math Nelder–Mead simplex optimization.
+ */
public Polygon computeMIR() {
- int i = 0;
- int same = 0;
+ int gen = 0, stableCount = 0;
double lastFitness = Double.MIN_VALUE;
- while (i++ < GENERATIONS) {
+ // Evolve the swarm until convergence criterion or maximum generation reached.
+ while (gen++ < GENERATIONS) {
swarm.evolve();
-
- // Periodically reinitialize worst performing particles
- if (i % 50 == 0) {
+ if (gen % 50 == 0) {
reinitializeWorstParticles(swarm);
}
-
- if (swarm.getBestFitness() == lastFitness) {
- if (same++ > 50) {
+ double bestFitness = swarm.getBestFitness();
+ if (bestFitness == lastFitness) {
+ if (stableCount++ > 50) {
break;
}
} else {
- same = 0;
+ stableCount = 0;
}
- lastFitness = swarm.getBestFitness();
+ lastFitness = bestFitness;
}
- return getBestRectangleResult(swarm);
- }
-
- private Polygon getBestRectangleResult(Swarm swarm) {
- return rectFromCoords(swarm.getBestPosition());
+ // Retrieve the best candidate from the swarm
+ double[] bestCandidate = swarm.getBestPosition().clone();
+ // Refine using Apache Commons Math
+ bestCandidate = refineCandidate(bestCandidate);
+ return rectFromCoords(bestCandidate);
}
private void reinitializeWorstParticles(Swarm swarm) {
Particle[] particles = swarm.getParticles();
double[] bestPos = swarm.getBestPosition();
-
- // Reinitialize bottom 10% of particles
int reinitCount = SWARM_SIZE / 10;
+ double[] maxPos = swarm.getMaxPosition();
+ double[] minPos = swarm.getMinPosition();
for (int i = 0; i < reinitCount; i++) {
double[] newPos = new double[5];
- // Explore around current best solution
for (int j = 0; j < 5; j++) {
- double range = (swarm.getMaxPosition()[j] - swarm.getMinPosition()[j]) * 0.1;
- newPos[j] = bestPos[j] + (Math.random() - 0.5) * range;
- // Ensure bounds
- newPos[j] = Math.max(swarm.getMinPosition()[j], Math.min(swarm.getMaxPosition()[j], newPos[j]));
+ double range = (maxPos[j] - minPos[j]) * 0.1;
+ newPos[j] = bestPos[j] + (RANDOM.nextDouble() - 0.5) * range;
+ // Clamp within bounds.
+ newPos[j] = Math.max(minPos[j], Math.min(maxPos[j], newPos[j]));
}
particles[i].setPosition(newPos);
particles[i].setVelocity(new double[] { 0, 0, 0, 0, 0 });
}
}
+ /**
+ * Uses Apache Commons Math optimizer (Nelder–Mead Simplex) to refine the best
+ * candidate. We define a MultivariateFunction that returns the negative fitness
+ * for minimization.
+ */
+ private double[] refineCandidate(double[] candidate) {
+ // We want to maximize the fitness, so we minimize negative fitness.
+ MultivariateFunction objective = new MultivariateFunction() {
+ public double value(double[] point) {
+ // Return the negative fitness
+ return -fitnessFunction.evaluate(point);
+ }
+ };
+
+ // No additional bounds mapping is used here since our candidate is already
+ // inside the feasible region.
+ // However, you can wrap this with a MultivariateFunctionMappingAdapter if
+ // needed.
+ MultivariateOptimizer optimizer = new SimplexOptimizer(1e-8, 1e-8);
+ NelderMeadSimplex simplex = new NelderMeadSimplex(candidate.length);
+ try {
+ PointValuePair result = optimizer.optimize(new MaxEval(1000), new ObjectiveFunction(objective), GoalType.MINIMIZE, new InitialGuess(candidate),
+ simplex);
+ return result.getPoint();
+ } catch (Exception e) {
+ // If the optimizer fails, return the original candidate.
+ return candidate;
+ }
+ }
+
+ /**
+ * Builds a rectangle Polygon based on a candidate parameter vector.
+ *
+ * The candidate parameters are: [x, y, width, height, angle] and the rectangle
+ * corners are computed from these values.
+ */
+ private static Polygon rectFromCoords(double[] position) {
+ double x = position[0];
+ double y = position[1];
+ double w = position[2];
+ double h = position[3];
+ double a = position[4];
+
+ double cosA = FastMath.cos(a);
+ double sinA = FastMath.sin(a);
+ double dx = cosA * w;
+ double dy = sinA * w;
+ double dx2 = -sinA * h;
+ double dy2 = cosA * h;
+
+ Coordinate[] coords = new Coordinate[5];
+ coords[0] = new Coordinate(x, y);
+ coords[1] = new Coordinate(x + dx, y + dy);
+ coords[2] = new Coordinate(x + dx + dx2, y + dy + dy2);
+ coords[3] = new Coordinate(x + dx2, y + dy2);
+ coords[4] = coords[0];
+
+ return GEOM_FACTORY.createPolygon(coords);
+ }
+
+ /**
+ * Evaluates each candidate rectangle. A candidate rectangle (constructed from
+ * position) is given a fitness equal to its area, with a small bonus for
+ * non-square aspect ratios. If the rectangle is not properly contained in the
+ * geometry, it returns 0.
+ */
private class RectangleFitness extends FitnessFunction {
- private PreparedGeometry geometry;
+ private final PreparedGeometry geometry;
- RectangleFitness(Geometry geometry) {
- this.geometry = PreparedGeometryFactory.prepare(geometry);
+ RectangleFitness(Geometry geom) {
+ this.geometry = PreparedGeometryFactory.prepare(geom);
}
@Override
public double evaluate(double[] position) {
- final double w = position[2];
- final double h = position[3];
- if (!geometry.containsProperly(rectFromCoords(position))) {
+ double w = position[2];
+ double h = position[3];
+ Polygon rect = rectFromCoords(position);
+ if (!geometry.containsProperly(rect)) {
return 0;
}
- // Add small bonus for non-square shapes to encourage elongated rectangles
- double aspectRatio = Math.max(w / h, h / w);
- return h * w * (1 + ASPECT_WEIGHT * (aspectRatio - 1));
+ double aspectRatio = (w > h) ? w / h : h / w;
+ return w * h * (1 + ASPECT_WEIGHT * (aspectRatio - 1));
}
}
+ // Particle candidate for the swarm; note the dimension is 5.
private class RectangleCandidate extends Particle {
public RectangleCandidate() {
super(5);
@@ -173,31 +256,4 @@ public Object selfFactory() {
return new RectangleCandidate();
}
}
-
- private static final Polygon rectFromCoords(double[] position) {
- final double x = position[0];
- final double y = position[1];
- final double w = position[2];
- final double h = position[3];
- final double a = position[4];
-
- Coordinate[] coords = new Coordinate[5];
- PVector base = new PVector((float) x, (float) y);
- PVector dir1 = new PVector((float) FastMath.cos(a), (float) FastMath.sin(a));
- PVector dir2 = dir1.copy().rotate((float) Math.PI * 0.5f);
- PVector base2 = base.copy().add(dir1.copy().mult((float) w));
- coords[0] = coordFromPVector(base);
- coords[1] = coordFromPVector(base2);
- dir2.mult((float) h);
- coords[2] = coordFromPVector(base2.add(dir2));
- coords[3] = coordFromPVector(base.add(dir2));
- coords[4] = coords[0];
-
- return GEOM_FACTORY.createPolygon(coords);
- }
-
- private static final Coordinate coordFromPVector(PVector p) {
- return new Coordinate(p.x, p.y);
- }
-
}
diff --git a/src/main/java/micycle/pgs/commons/MaximumInscribedTriangle.java b/src/main/java/micycle/pgs/commons/MaximumInscribedTriangle.java
new file mode 100644
index 00000000..21adb40f
--- /dev/null
+++ b/src/main/java/micycle/pgs/commons/MaximumInscribedTriangle.java
@@ -0,0 +1,405 @@
+package micycle.pgs.commons;
+
+import java.util.List;
+import java.util.Random;
+import java.util.concurrent.ThreadLocalRandom;
+
+import org.apache.commons.math3.analysis.MultivariateFunction;
+import org.apache.commons.math3.optim.InitialGuess;
+import org.apache.commons.math3.optim.MaxEval;
+import org.apache.commons.math3.optim.PointValuePair;
+import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
+import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction;
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex;
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer;
+import org.locationtech.jts.algorithm.construct.MaximumInscribedCircle;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.Envelope;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.LineSegment;
+import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.geom.prep.PreparedPolygon;
+import org.locationtech.jts.index.strtree.STRtree;
+
+import net.metaopt.swarm.FitnessFunction;
+import net.metaopt.swarm.pso.Particle;
+import net.metaopt.swarm.pso.Swarm;
+
+/**
+ * Finds an approximate largest area triangle of arbitrary orientation in a
+ * concave polygon via particle swarm optimisation.
+ *
+ * @author Michael Carleton
+ *
+ */
+public class MaximumInscribedTriangle {
+
+ private static final int SWARM_SIZE = 3000;
+ private static final int MAX_GENERATIONS = 1000;
+ private static final GeometryFactory GEOM_FACTORY = new GeometryFactory();
+ private static final Random RANDOM = ThreadLocalRandom.current();
+
+ private final Swarm swarm;
+ private final TriangleFitness fitnessFunction;
+
+ double minArea = 0;
+
+ public MaximumInscribedTriangle(Polygon polygon) {
+ fitnessFunction = new TriangleFitness(polygon);
+ // For triangle, our candidate vector is of dimension 6: [x1,y1,x2,y2,x3,y3].
+ swarm = new ParallelSwarm(SWARM_SIZE, new TriangleCandidate(), fitnessFunction);
+
+ // Use the envelope of the polygon as the search bounds.
+ Envelope e = polygon.getEnvelopeInternal();
+ double[] minPosition = new double[] { e.getMinX(), e.getMinY(), e.getMinX(), e.getMinY(), e.getMinX(), e.getMinY() };
+ double[] maxPosition = new double[] { e.getMaxX(), e.getMaxY(), e.getMaxX(), e.getMaxY(), e.getMaxX(), e.getMaxY() };
+
+ var r = MaximumInscribedCircle.getRadiusLine(polygon, 1).getLength();
+ /*
+ * Area of equilateral triangle inscribed in MIC. Sets a minimum bound to the
+ * MIT area.
+ */
+ minArea = (3 * Math.sqrt(3) / 4) * Math.pow(r, 2);
+
+ swarm.setMinPosition(minPosition);
+ swarm.setMaxPosition(maxPosition);
+ swarm.init();
+ initializeSwarm(swarm, e);
+ }
+
+ /**
+ * Initialize the swarm particles. Several groups use different heuristics: -
+ * Some choose points on polygon envelope edges. - Others choose random points
+ * inside the envelope.
+ */
+ private void initializeSwarm(Swarm swarm, Envelope e) {
+ Particle[] particles = swarm.getParticles();
+ int particlesPerGroup = SWARM_SIZE / 4;
+ double minX = e.getMinX();
+ double minY = e.getMinY();
+ double width = e.getWidth();
+ double height = e.getHeight();
+
+ for (int i = 0; i < SWARM_SIZE; i++) {
+ double[] position = new double[6];
+
+ if (i < particlesPerGroup) {
+ // Group 1: Pick vertices on the envelope boundary.
+ position[0] = minX + RANDOM.nextDouble() * width;
+ position[1] = RANDOM.nextBoolean() ? e.getMinY() : e.getMaxY();
+ position[2] = minX + RANDOM.nextDouble() * width;
+ position[3] = RANDOM.nextBoolean() ? e.getMinY() : e.getMaxY();
+ position[4] = minX + RANDOM.nextDouble() * width;
+ position[5] = RANDOM.nextBoolean() ? e.getMinY() : e.getMaxY();
+ } else if (i < 2 * particlesPerGroup) {
+ // Group 2: Random points in the envelope.
+ position[0] = minX + RANDOM.nextDouble() * width;
+ position[1] = minY + RANDOM.nextDouble() * height;
+ position[2] = minX + RANDOM.nextDouble() * width;
+ position[3] = minY + RANDOM.nextDouble() * height;
+ position[4] = minX + RANDOM.nextDouble() * width;
+ position[5] = minY + RANDOM.nextDouble() * height;
+ } else if (i < 3 * particlesPerGroup) {
+ // Group 3: Two points from one edge and one from the opposite edge.
+ position[0] = minX + RANDOM.nextDouble() * width;
+ position[1] = e.getMinY();
+ position[2] = minX + RANDOM.nextDouble() * width;
+ position[3] = e.getMinY();
+ position[4] = minX + RANDOM.nextDouble() * width;
+ position[5] = e.getMaxY();
+ } else {
+ // Group 4: Random exploration.
+ position[0] = minX + RANDOM.nextDouble() * width;
+ position[1] = minY + RANDOM.nextDouble() * height;
+ position[2] = minX + RANDOM.nextDouble() * width;
+ position[3] = minY + RANDOM.nextDouble() * height;
+ position[4] = minX + RANDOM.nextDouble() * width;
+ position[5] = minY + RANDOM.nextDouble() * height;
+ }
+
+ particles[i].setPosition(position);
+ particles[i].setVelocity(new double[] { 0, 0, 0, 0, 0, 0 });
+ particles[i].setBestPosition(position.clone());
+ }
+ }
+
+ /**
+ * The main entry method which runs the swarm evolution, then refines the best
+ * candidate using Apache Commons Math's Nelder–Mead optimizer.
+ */
+ public Polygon computeMIT() {
+ int gen = 0, stableCount = 0;
+ double lastFitness = Double.MIN_VALUE;
+
+ // Global search via swarm evolution.
+ while (gen++ < MAX_GENERATIONS) {
+ swarm.evolve();
+
+ // Occasionally reinitialize worst performing particles.
+ if (gen % 100 == 0) {
+ reinitializeWorstParticles(swarm);
+ }
+ double bestFitness = swarm.getBestFitness();
+ if (bestFitness == lastFitness) {
+ if (stableCount++ > 75) {
+ break;
+ }
+ } else {
+ stableCount = 0;
+ }
+ lastFitness = bestFitness;
+ }
+
+ double[] bestCandidate = swarm.getBestPosition().clone();
+ bestCandidate = refineCandidate(bestCandidate);
+ return triangleFromCoords(bestCandidate);
+ }
+
+ /**
+ * Reinitialize the bottom 10% of swarm particles around the best candidate.
+ */
+ private void reinitializeWorstParticles(Swarm swarm) {
+ Particle[] particles = swarm.getParticles();
+ double[] bestPos = swarm.getBestPosition();
+ int reinitCount = SWARM_SIZE / 10;
+ double[] maxPos = swarm.getMaxPosition();
+ double[] minPos = swarm.getMinPosition();
+ for (int i = 0; i < reinitCount; i++) {
+ double[] newPos = new double[6];
+ for (int j = 0; j < 6; j++) {
+ double range = (maxPos[j] - minPos[j]) * 0.1;
+ newPos[j] = bestPos[j] + (RANDOM.nextDouble() - 0.5) * range;
+ newPos[j] = Math.max(minPos[j], Math.min(maxPos[j], newPos[j]));
+ }
+ particles[i].setPosition(newPos);
+ particles[i].setVelocity(new double[] { 0, 0, 0, 0, 0, 0 });
+ }
+ }
+
+ /**
+ * Uses Nelder–Mead simplex optimization (via Apache Commons Math) to refine the
+ * best candidate. The objective function is defined as the negative fitness.
+ */
+ private double[] refineCandidate(double[] candidate) {
+ MultivariateFunction objective = point -> -fitnessFunction.evaluate(point);
+
+ SimplexOptimizer optimizer = new SimplexOptimizer(1e-8, 1e-8);
+ NelderMeadSimplex simplex = new NelderMeadSimplex(candidate.length);
+ try {
+ PointValuePair result = optimizer.optimize(new MaxEval(1000), new ObjectiveFunction(objective), GoalType.MINIMIZE, new InitialGuess(candidate),
+ simplex);
+ return result.getPoint();
+ } catch (Exception e) {
+ return candidate;
+ }
+ }
+
+ /**
+ * Constructs a triangle polygon from a candidate vector. Candidate parameters:
+ * [x1, y1, x2, y2, x3, y3]
+ */
+ private static Polygon triangleFromCoords(double[] position) {
+ // Create three coordinates for the triangle.
+ Coordinate p0 = new Coordinate(position[0], position[1]);
+ Coordinate p1 = new Coordinate(position[2], position[3]);
+ Coordinate p2 = new Coordinate(position[4], position[5]);
+ // Ensure the ring is closed.
+ Coordinate[] coords = new Coordinate[] { p0, p1, p2, p0 };
+ return GEOM_FACTORY.createPolygon(coords);
+ }
+
+ /**
+ * The fitness function for candidate triangles in a concave polygon. Instead of
+ * calling PreparedGeometry.containsProperly(…), we use a custom “fast–path”
+ * test: 1. All triangle vertices must be inside the polygon (using
+ * ray–casting). 2. None of the triangle edges may intersect any polygon
+ * boundary edge.
+ *
+ * If both tests pass, return the triangle’s absolute area.
+ */
+ private class TriangleFitness extends FitnessFunction {
+ // Cached polygon vertices (exterior ring) for fast point testing.
+ private double[] polyX;
+ private double[] polyY;
+ // Spatial index over the polygon’s edges.
+ private STRtree polygonEdgeIndex;
+ private final boolean hasHoles;
+ private PreparedPolygon cache;
+
+ TriangleFitness(Polygon polygon) {
+ hasHoles = polygon.getNumInteriorRing() > 0;
+ if (hasHoles) {
+ cache = new PreparedPolygon(polygon);
+ } else { // faster approach (but doesn't detect holes)
+ Coordinate[] coords = polygon.getExteriorRing().getCoordinates();
+ int n = coords.length - 1; // last coordinate is a duplicate of the first
+ polyX = new double[n];
+ polyY = new double[n];
+ for (int i = 0; i < n; i++) {
+ polyX[i] = coords[i].x;
+ polyY[i] = coords[i].y;
+ }
+ // Build an STRtree on polygon edges.
+ polygonEdgeIndex = new STRtree();
+ for (int i = 0; i < n; i++) {
+ Coordinate p0 = coords[i];
+ Coordinate p1 = coords[(i + 1) % n];
+ LineSegment seg = new LineSegment(p0, p1);
+ polygonEdgeIndex.insert(new Envelope(seg.p0, seg.p1), seg);
+ }
+ polygonEdgeIndex.build();
+ }
+ }
+
+ @Override
+ public double evaluate(double[] position) {
+ double x0 = position[0], y0 = position[1];
+ double x1 = position[2], y1 = position[3];
+ double x2 = position[4], y2 = position[5];
+ // Compute absolute area using the cross–product formula.
+ double area = Math.abs(0.5 * (x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1)));
+ if (area < minArea) {
+ return 0;
+ }
+
+ if (hasHoles) {
+ Polygon tri = triangleFromCoords(position);
+ if (cache.containsProperly(tri)) {
+ return area;
+ } else {
+ return 0;
+ }
+ } else {
+ // 1. Check that all three vertices are strictly inside the polygon.
+ if (!pointInPolygon(x0, y0, polyX, polyY) || !pointInPolygon(x1, y1, polyX, polyY)
+ || !pointInPolygon(x2, y2, polyX, polyY)) {
+ return 0;
+ }
+
+ // 2. Check for any intersection between triangle edges and polygon boundary.
+ // Define the triangle’s three edges.
+ if (edgeIntersectsPolygon(x0, y0, x1, y1) || edgeIntersectsPolygon(x1, y1, x2, y2) || edgeIntersectsPolygon(x2, y2, x0, y0)) {
+ return 0;
+ }
+
+ return area;
+ }
+ }
+
+ /**
+ * Returns true if the triangle edge from (ax,ay) to (bx,by) intersects a
+ * polygon edge.
+ */
+ private boolean edgeIntersectsPolygon(double ax, double ay, double bx, double by) {
+ Envelope edgeEnv = new Envelope(ax, bx, ay, by);
+ List candidates = polygonEdgeIndex.query(edgeEnv);
+ LineSegment triangleEdge = new LineSegment(new Coordinate(ax, ay), new Coordinate(bx, by));
+ for (Object obj : candidates) {
+ LineSegment polyEdge = (LineSegment) obj;
+ if (segmentsIntersect(triangleEdge, polyEdge)) {
+ return true;
+ }
+ }
+ return false;
+ }
+
+ /**
+ * Determines whether the two given non-collinear line segments intersect.
+ *
+ *
+ * This fast geometric method assumes that the bounding boxes (envelopes) of the
+ * segments already intersect, so no additional envelope intersection test is
+ * performed. It also assumes that no three endpoints are collinear, ensuring
+ * that none of the computed cross products are zero.
+ *
+ *
+ *
+ * The algorithm proceeds by computing the cross products to determine the
+ * relative orientations of the endpoints of each segment with respect to the
+ * line defined by the other segment. Specifically, let segment 1 be defined by
+ * points A and B, and segment 2 by points C and D. The method computes:
+ *
+ *
+ *
+ * - d1 = cross product of (D - C) and (A - C)
+ * - d2 = cross product of (D - C) and (B - C)
+ * - d3 = cross product of (B - A) and (C - A)
+ * - d4 = cross product of (B - A) and (D - A)
+ *
+ *
+ *
+ * If d1 and d2 have the same sign, then both A and B lie on the same side of
+ * the line through C and D, meaning segment 1 does not cross segment 2.
+ * Similarly, if d3 and d4 have the same sign, segment 2 does not cross segment
+ * 1. Therefore, the segments intersect if and only if A and B lie on opposite
+ * sides of the line through C and D, and C and D lie on opposite sides of the
+ * line through A and B.
+ *
+ *
+ * @param seg1 the first line segment
+ * @param seg2 the second line segment
+ * @return {@code true} if the segments intersect; {@code false} otherwise
+ */
+ private boolean segmentsIntersect(LineSegment seg1, LineSegment seg2) {
+ // Cache coordinates
+ double ax = seg1.p0.x, ay = seg1.p0.y, bx = seg1.p1.x, by = seg1.p1.y;
+ double cx = seg2.p0.x, cy = seg2.p0.y, dx = seg2.p1.x, dy = seg2.p1.y;
+
+ // Compute differences for segment 2 (t)
+ double rdx = dx - cx, rdy = dy - cy;
+ // Compute cross products for seg1 endpoints relative to seg2's line
+ double d1 = rdx * (ay - cy) - rdy * (ax - cx);
+ double d2 = rdx * (by - cy) - rdy * (bx - cx);
+
+ // If d1 and d2 are of the same sign, seg1 does not cross seg2.
+ if ((d1 > 0) == (d2 > 0)) {
+ return false;
+ }
+
+ // Compute differences for segment 1 (s)
+ double sdx = bx - ax, sdy = by - ay;
+ // Compute cross products for seg2 endpoints relative to seg1's line
+ double d3 = sdx * (cy - ay) - sdy * (cx - ax);
+ double d4 = sdx * (dy - ay) - sdy * (dx - ax);
+
+ // The segments intersect if and only if seg2's endpoints lie on opposite
+ // sides of seg1's line.
+ return (d3 > 0) != (d4 > 0);
+ }
+
+ /**
+ * A standard ray-casting point-in-polygon test.
+ *
+ * @param x the test point x coordinate
+ * @param y the test point y coordinate
+ * @param polyX the array of polygon vertex x coordinates
+ * @param polyY the array of polygon vertex y coordinates
+ * @return true if the point is inside
+ */
+ private boolean pointInPolygon(double x, double y, double[] polyX, double[] polyY) {
+ boolean inside = false;
+ int n = polyX.length;
+ for (int i = 0, j = n - 1; i < n; j = i++) {
+ // Check if point is between the y-interval of the edge.
+ if (((polyY[i] > y) != (polyY[j] > y)) && (x < (polyX[j] - polyX[i]) * (y - polyY[i]) / (polyY[j] - polyY[i]) + polyX[i])) {
+ inside = !inside;
+ }
+ }
+ return inside;
+ }
+ }
+
+ /**
+ * A Particle subclass for triangle candidates. The dimensionality is 6.
+ */
+ private class TriangleCandidate extends Particle {
+ public TriangleCandidate() {
+ super(6);
+ }
+
+ @Override
+ public Object selfFactory() {
+ return new TriangleCandidate();
+ }
+ }
+}
diff --git a/src/main/java/micycle/pgs/commons/MultiplicativelyWeightedVoronoi.java b/src/main/java/micycle/pgs/commons/MultiplicativelyWeightedVoronoi.java
index 1928320e..bf372a5e 100644
--- a/src/main/java/micycle/pgs/commons/MultiplicativelyWeightedVoronoi.java
+++ b/src/main/java/micycle/pgs/commons/MultiplicativelyWeightedVoronoi.java
@@ -14,6 +14,7 @@
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.operation.overlayng.OverlayNG;
+import org.locationtech.jts.operation.overlayng.RingClipper;
import org.locationtech.jts.operation.polygonize.Polygonizer;
import org.locationtech.jts.operation.union.UnaryUnionOp;
import org.tinspin.index.Index.PointEntryKnn;
@@ -90,9 +91,9 @@ private static List getMWVDFast(List sites, Envelope exten
// intersection(all ap_circles containing s)-union(all ap_circles not containing
// s)
- sites.sort((s1,s2) -> Double.compare(s1.z, s2.z));
-
+ sites.sort((s1, s2) -> Double.compare(s1.z, s2.z));
final Geometry extentGeometry = geometryFactory.toGeometry(extent);
+ final RingClipper rc = new RingClipper(extent);
return sites.parallelStream().map(site -> { // NOTE parallel
// s2 dominates s1 (hence s1 contained in apollo circle)
@@ -124,7 +125,7 @@ private static List getMWVDFast(List sites, Envelope exten
if (inCircleData.isEmpty()) {
// if no incircles, cell is simply defined by difference between plane and
// exCircles
- exCircleData.forEach(c -> exCircles.add(createCircle(c[0], c[1], c[2])));
+ exCircleData.forEach(c -> exCircles.add(createClippedCircle(c[0], c[1], c[2], rc)));
var outerDom = UnaryUnionOp.union(exCircles);
return extentGeometry.difference(outerDom);
}
@@ -149,7 +150,7 @@ private static List getMWVDFast(List sites, Envelope exten
});
}
- essentialCircles.forEach(c -> inCircles.add(createCircle(c[0], c[1], c[2])));
+ essentialCircles.forEach(c -> inCircles.add(createClippedCircle(c[0], c[1], c[2], rc)));
// intersect all inCircles to find the dominant region for this site
var localDominance = inCircles.stream().reduce((geom1, geom2) -> OverlayNG.overlay(geom1, geom2, OverlayNG.INTERSECTION)).get();
@@ -172,17 +173,17 @@ private static List getMWVDFast(List sites, Envelope exten
}).collect(Collectors.toList());
/*
- * NOTE optmisation, similar to inCircleData optimisation, but this time, remove
- * any SMALLER circles that are contained by another.
+ * NOTE optimisation, similar to inCircleData optimisation, but this time,
+ * remove any SMALLER circles that are contained by another.
*/
exCircleData = filterContainedCircles(exCircleData);
- exCircleData.forEach(c -> exCircles.add(createCircle(c[0], c[1], c[2])));
+ exCircleData.forEach(c -> exCircles.add(createClippedCircle(c[0], c[1], c[2], rc)));
if (exCircles.isEmpty()) {
- return localDominance.intersection(extentGeometry);
+ return localDominance;
} else {
var outerDom = UnaryUnionOp.union(exCircles);
- return localDominance.difference(outerDom).intersection(extentGeometry);
+ return localDominance.difference(outerDom);
}
}).toList();
}
@@ -310,7 +311,7 @@ private static Geometry apolloniusCircle(Coordinate s1, Coordinate s2, double w1
Coordinate center = new Coordinate(circle[0], circle[1]); // often lies outside bounds
double radius = circle[2];
double distance = s1.distance(center);
- localDominanceCircle = createCircle(center.x, center.y, radius);
+ localDominanceCircle = createClippedCircle(center.x, center.y, radius, null);
/*
* The circle will either enclose site 1 (i.e. it's dominated by site 2), or
* will bend away from site 1 (enclosing and dominating site 2).
@@ -384,11 +385,11 @@ private static double[] calculateWeightedApollonius(Coordinate s1, Coordinate s2
final double d = Math.sqrt(((s1x - s2x) * (s1x - s2x) + (s1y - s2y) * (s1y - s2y)));
// NOTE r can be huge (as circle may tend towards straight line)
final double r = Math.abs(w1 * w2 * d * den);
-
+
return new double[] { cx, cy, r };
}
- private static Polygon createCircle(double x, double y, double r) {
+ private static Polygon createClippedCircle(double x, double y, double r, RingClipper rc) {
final double maxDeviation = 0.49;
// Calculate the number of points based on the radius and maximum deviation.
int nPts = (int) Math.ceil(2 * Math.PI / Math.acos(1 - maxDeviation / r));
@@ -408,6 +409,8 @@ private static Polygon createCircle(double x, double y, double r) {
}
pts[nPts] = new Coordinate(pts[0]); // Close the circle
- return geometryFactory.createPolygon(pts);
+ // NOTE clip the circle to bounds now. slightly speeds up 2d boolean ops later
+ // on.
+ return geometryFactory.createPolygon(rc == null ? pts : rc.clip(pts));
}
}
\ No newline at end of file
diff --git a/src/main/java/micycle/pgs/commons/PEdge.java b/src/main/java/micycle/pgs/commons/PEdge.java
index 1853265f..2023a603 100644
--- a/src/main/java/micycle/pgs/commons/PEdge.java
+++ b/src/main/java/micycle/pgs/commons/PEdge.java
@@ -11,7 +11,7 @@
* @author Michael Carleton
*
*/
-public class PEdge {
+public class PEdge implements Comparable {
public final PVector a, b;
@@ -70,8 +70,7 @@ public PEdge slice(double from, double to) {
to = 1;
}
if (from < 0 || to > 1 || from > to) {
- throw new IllegalArgumentException(
- "Parameters 'from' and 'to' must be between 0 and 1, and 'from' must be less than or equal to 'to'.");
+ throw new IllegalArgumentException("Parameters 'from' and 'to' must be between 0 and 1, and 'from' must be less than or equal to 'to'.");
}
final PVector pointFrom;
@@ -123,4 +122,24 @@ public String toString() {
private static boolean equals(PVector a, PVector b) {
return a.x == b.x && a.y == b.y;
}
+
+ @Override
+ public int compareTo(PEdge other) {
+ PVector thisMidpoint = midpoint();
+ PVector otherMidpoint = other.midpoint();
+ return comparePVectors(thisMidpoint, otherMidpoint);
+ }
+
+ /**
+ * Helper method to compare two PVectors lexicographically.
+ */
+ private int comparePVectors(PVector v1, PVector v2) {
+ if (v1.x != v2.x) {
+ return Float.compare(v1.x, v2.x);
+ }
+ if (v1.y != v2.y) {
+ return Float.compare(v1.y, v2.y);
+ }
+ return Float.compare(v1.z, v2.z);
+ }
}
\ No newline at end of file
diff --git a/src/main/java/micycle/pgs/commons/PMesh.java b/src/main/java/micycle/pgs/commons/PMesh.java
index c0062cb8..375da2cc 100644
--- a/src/main/java/micycle/pgs/commons/PMesh.java
+++ b/src/main/java/micycle/pgs/commons/PMesh.java
@@ -372,7 +372,7 @@ private static class PMeshVertex {
final PVector originalVertex;
final PVector smoothedVertex;
- boolean onBoundary;
+ boolean onBoundary = false;
List neighbors;
public PMeshVertex(PVector v) {
diff --git a/src/main/java/micycle/pgs/commons/ParallelSwarm.java b/src/main/java/micycle/pgs/commons/ParallelSwarm.java
new file mode 100644
index 00000000..3e41ee56
--- /dev/null
+++ b/src/main/java/micycle/pgs/commons/ParallelSwarm.java
@@ -0,0 +1,62 @@
+package micycle.pgs.commons;
+
+import java.util.stream.IntStream;
+
+import net.metaopt.swarm.FitnessFunction;
+import net.metaopt.swarm.pso.Particle;
+import net.metaopt.swarm.pso.Swarm;
+
+/**
+ * A particle swarm that evaluates particle fitness in parallel for improved
+ * performance.
+ *
+ * @author Michael Carleton
+ */
+public class ParallelSwarm extends Swarm {
+
+ public ParallelSwarm(int numberOfParticles, Particle sampleParticle, FitnessFunction fitnessFunction) {
+ super(numberOfParticles, sampleParticle, fitnessFunction);
+ }
+
+ @Override
+ public void evaluate() {
+ // Initialize bestFitness on first run.
+ if (Double.isNaN(bestFitness)) {
+ bestFitness = (fitnessFunction.isMaximize() ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY);
+ bestParticleIndex = -1;
+ }
+
+ final int len = particles.length;
+ final double[] fitnesses = new double[len];
+
+ // Evaluate each particle's fitness in parallel.
+ IntStream.range(0, len).parallel().forEach(i -> {
+ double fit = fitnessFunction.evaluate(particles[i]);
+ fitnesses[i] = fit;
+ });
+
+ // Now process the results sequentially to update the global best and
+ // neighborhood.
+ for (int i = 0; i < len; i++) {
+ numEvaluations++;
+ double fit = fitnesses[i];
+
+ // Update 'best global' position if this particle is better.
+ if (fitnessFunction.isBetterThan(bestFitness, fit)) {
+ bestFitness = fit;
+ bestParticleIndex = i;
+ if (bestPosition == null) {
+ bestPosition = new double[sampleParticle.getDimension()];
+ }
+ particles[i].copyPosition(bestPosition);
+ }
+
+ // Update 'best neighborhood' information.
+ if (neighborhood != null) {
+ neighborhood.update(this, particles[i]);
+ }
+
+ }
+ }
+
+}
diff --git a/src/main/java/micycle/pgs/commons/TangencyPack.java b/src/main/java/micycle/pgs/commons/TangencyPack.java
index 1c62d494..c84220b2 100644
--- a/src/main/java/micycle/pgs/commons/TangencyPack.java
+++ b/src/main/java/micycle/pgs/commons/TangencyPack.java
@@ -6,7 +6,6 @@
import java.util.HashSet;
import java.util.List;
import java.util.Map;
-import java.util.Map.Entry;
import java.util.Set;
import org.apache.commons.math3.complex.Complex;
@@ -43,6 +42,7 @@
*
* @author Michael Carleton
*/
+
public class TangencyPack {
/*-
@@ -52,89 +52,45 @@ public class TangencyPack {
* https://github.com/kensmath/CirclePack/blob/CP-develop/src/rePack/EuclPacker.java
*/
- private static final double TOLERANCE = 1 + 1e-8;
private static final double TWO_PI = Math.PI * 2;
private final IIncrementalTin triangulation;
- /**
- * Maps a vertex to a list of it neighbouring vertices; the neighbour list is
- * ordered radially around the given vertex.
- */
- private Map> flowers;
- /**
- * The radius of each circle (including boundary circles).
- */
- private Object2DoubleOpenHashMap radii;
- private double[] boundaryRadii;
- private Map placements = new HashMap<>();
+
+ private int[][] flowersIds; // Each row contains neighbor IDs for a flower
+ private int[] flowerSizes; // Precomputed sizes of each flower
+ private double[] delArray; // Precomputed sin(π / n) for each flower
+ private double[] factorDelArray; // Precomputed (1 - del)/del for each flower
+
+ private double[] radiiArray;
+ private Complex[] placementsArray;
private List circles;
private Vertex centralVertex;
+ private double[] boundaryRadii;
+
+ private List allVertices;
+ private Map vertexToId;
+ private List interiorVertices;
+ private Map interiorVertexIndex;
+ private int[] interiorVertexIds;
- /**
- * Creates a circle packing using tangancies specified by a triangulation.
- *
- * @param triangulation Pattern of tangencies; vertices connected by an edge in
- * the triangulation represent tangent circles in the
- * packing
- * @param boundaryRadii Radii of the circles (same for every circle) associated
- * with the boundary/perimeter vertices of the
- * triangulation
- */
public TangencyPack(IIncrementalTin triangulation, double boundaryRadii) {
this.triangulation = triangulation;
this.boundaryRadii = new double[] { boundaryRadii };
init();
}
- /**
- * Creates a circle packing using tangancies specified by a triangulation.
- *
- * @param triangulation Pattern of tangencies; vertices connected by an edge in
- * the triangulation represent tangent circles in the
- * packing
- * @param boundaryRadii List of radii values of the circles associated with the
- * boundary/perimeter vertices of the triangulation. The
- * list may have fewer radii than the number of boundary
- * vertices; in this case, boundary radii will wrap around
- * the list
- */
public TangencyPack(IIncrementalTin triangulation, List boundaryRadii) {
this.triangulation = triangulation;
this.boundaryRadii = boundaryRadii.stream().mapToDouble(Double::doubleValue).toArray();
init();
}
- /**
- * Creates a circle packing using tangancies specified by a triangulation.
- *
- * @param triangulation Pattern of tangencies; vertices connected by an edge in
- * the triangulation represent tangent circles in the
- * packing
- * @param boundaryRadii Array of radii values of the circles associated with the
- * boundary/perimeter vertices of the triangulation. The
- * list may have fewer radii than the number of boundary
- * vertices; in this case, boundary radii will wrap around
- * the list
- */
public TangencyPack(IIncrementalTin triangulation, double[] boundaryRadii) {
this.triangulation = triangulation;
this.boundaryRadii = boundaryRadii;
init();
}
- /**
- * Computes and returns a circle packing for the configuration of tangencies
- * given by the triangulation.
- *
- * @return a list of PVectors, each representing one circle: (.x, .y) represent
- * the center point and .z represents radius.
- */
- public List pack() {
- computeRadii();
- computeCenters();
- return circles;
- }
-
private void init() {
Set perimeterVertices = new HashSet<>();
triangulation.getPerimeter().forEach(e -> {
@@ -142,328 +98,350 @@ private void init() {
perimeterVertices.add(e.getB());
});
- flowers = new HashMap<>();
- radii = new Object2DoubleOpenHashMap<>(triangulation.getVertices().size());
+ allVertices = new ArrayList<>(triangulation.getVertices());
+ vertexToId = new HashMap<>();
+ for (int i = 0; i < allVertices.size(); i++) {
+ vertexToId.put(allVertices.get(i), i);
+ }
+
+ radiiArray = new double[allVertices.size()];
+ flowersIds = new int[allVertices.size()][];
+ interiorVertices = new ArrayList<>();
+ interiorVertexIndex = new HashMap<>();
SimpleGraph graph = PGS_Triangulation.toTinfourGraph(triangulation);
NeighborCache neighbors = new NeighborCache<>(graph);
- final PVector meanVertexPos = new PVector();
- int index = 0;
- for (Vertex v : graph.vertexSet()) {
+ int boundaryIndex = 0;
+ PVector meanVertexPos = new PVector();
+ int vi = 0;
+ for (Vertex v : allVertices) {
if (perimeterVertices.contains(v)) {
- radii.put(v, boundaryRadii[index++ % boundaryRadii.length]);
+ radiiArray[vertexToId.get(v)] = boundaryRadii[boundaryIndex++ % boundaryRadii.length];
} else {
List flower = neighbors.neighborListOf(v);
- RadialComparator c = new RadialComparator(v);
- flower.sort(c);
- flowers.put(v, flower);
- radii.put(v, boundaryRadii[0] / 10);
+ flower.sort(new RadialComparator(v));
+ int[] flowerIds = new int[flower.size()];
+ for (int i = 0; i < flowerIds.length; i++) {
+ flowerIds[i] = vertexToId.get(flower.get(i));
+ }
+ flowersIds[vi++] = flowerIds;
+ interiorVertices.add(v);
+ int vid = vertexToId.get(v);
+ radiiArray[vid] = boundaryRadii[0] / 10;
meanVertexPos.add((float) v.x, (float) v.y);
}
}
- // pick a rather central vertex, so output is same on identical input
- meanVertexPos.div(flowers.size());
- double maxDist = Double.MAX_VALUE;
- for (Vertex v : flowers.keySet()) {
+ interiorVertexIds = new int[interiorVertices.size()];
+ for (int i = 0; i < interiorVertices.size(); i++) {
+ Vertex v = interiorVertices.get(i);
+ interiorVertexIds[i] = vertexToId.get(v);
+ interiorVertexIndex.put(v, i);
+ }
+
+ meanVertexPos.div(interiorVertices.size());
+ double minDist = Double.MAX_VALUE;
+ for (Vertex v : interiorVertices) {
double dist = v.getDistanceSq(meanVertexPos.x, meanVertexPos.y);
- if (dist < maxDist) {
- maxDist = dist;
+ if (dist < minDist) {
+ minDist = dist;
centralVertex = v;
}
}
- }
- /**
- * Find radii of circles using numerical relaxation. Circle radii converge
- * rapidly to a unique fixed point for which all flower angles are are within a
- * desired tolerance of 2π, at which point iteration stops and a packing is
- * found.
- *
- * @deprecated in favor of superstep solution
- */
- @Deprecated
- private void computeRadiiSimple() {
- double lastChange = TOLERANCE + 1;
- while (lastChange > TOLERANCE) {
- lastChange = 1.0;
- for (Vertex v : flowers.keySet()) {
- double theta = flower(v);
- lastChange = Math.max(lastChange, theta);
- }
+ // Precompute flower sizes and trigonometric terms
+ flowerSizes = new int[interiorVertices.size()];
+ delArray = new double[interiorVertices.size()];
+ factorDelArray = new double[interiorVertices.size()];
+
+ for (int i = 0; i < interiorVertices.size(); i++) {
+ int n = flowersIds[i].length;
+ flowerSizes[i] = n;
+ double del = FastMath.sin(Math.PI / n); // del = sin(π / n)
+ delArray[i] = del;
+ factorDelArray[i] = (1 - del) / del; // Precompute (1-del)/del
}
}
- /**
- * This method implements the super acceleration described in 'A circle packing
- * algorithm'.
- */
+ public List pack() {
+ computeRadiiSuperStep();
+ computeCenters();
+ return circles;
+ }
+
private void computeRadii() {
- final double ttoler = 3 * radii.size() * 1e-11;
- int key = 1; // initial superstep type
+ double ttoler;
+ // adapt tolerance based on input size. seems sufficient
+ if (interiorVertices.size() <= 100) {
+ ttoler = 1e-3; // Base case
+ } else {
+ double exponent = 3.5 + (interiorVertices.size()) / 100.0;
+ ttoler = Math.pow(10, -exponent);
+ }
+ int key = 1;
double accumErr2 = Double.MAX_VALUE;
- int localPasses = 1;
-
- while ((accumErr2 > ttoler && localPasses < 1000)) { // main loop
- Object2DoubleMap R1 = new Object2DoubleOpenHashMap<>(radii);
- double c1;
+ int localPasses = 0;
- double factor;
-
- do { // Make sure factor < 1.0
- c1 = computeAngleSums();
- c1 = Math.sqrt(c1);
+ while (accumErr2 > ttoler && localPasses < 3 * interiorVertices.size()) {
+ Object2DoubleMap R1 = new Object2DoubleOpenHashMap<>(interiorVertices.size());
+ for (Vertex v : interiorVertices) {
+ R1.put(v, radiiArray[vertexToId.get(v)]);
+ }
- factor = c1 / accumErr2;
- if (factor >= 1.0) {
- accumErr2 = c1;
- key = 1;
- }
- } while (factor >= 1.0);
+ double c1 = computeAngleSums();
+ c1 = Math.sqrt(c1);
+ if (c1 < ttoler) {
+ break;
+ }
- // ================= superstep calculation ====================
+ double factor = c1 / accumErr2;
+ if (factor >= 1.0) {
+ accumErr2 = c1;
+ key = 1;
+ localPasses++;
+ continue;
+ }
- Object2DoubleMap R2 = new Object2DoubleOpenHashMap<>(radii);
+ Object2DoubleMap R2 = new Object2DoubleOpenHashMap<>(interiorVertices.size());
+ for (Vertex v : interiorVertices) {
+ R2.put(v, radiiArray[vertexToId.get(v)]);
+ }
- // find maximum step one can safely take
double lmax = 10000;
- double fact0;
-
for (Vertex v : R1.keySet()) {
double r1 = R1.getDouble(v);
double r2 = R2.getDouble(v);
double rat = r2 - r1;
- double tr;
if (rat < 0) {
- lmax = (lmax < (tr = (-r2 / rat))) ? lmax : tr; // to keep R>0
+ double tr = -r2 / rat;
+ lmax = Math.min(lmax, tr);
}
}
- lmax = lmax / 2;
+ lmax /= 2;
- // do super step
- double m = 1;
- int sct = 1;
- int fct = 2;
double lambda;
- if (key == 1) { // type 1 SS
- lambda = m * factor;
- double mmax = 0.75 / (1 - factor); // upper limit on m
- double mm = 0.0;
- m = (mmax < (mm = (1 + 0.8 / (sct + 1)) * m)) ? mmax : mm;
- } else { // type 2 SS
- fact0 = 0.0;
- double ftol = 0.0;
- if (sct > fct && Math.abs(factor - fact0) < ftol) { // try SS-2
- lambda = factor / (1 - factor);
- sct = -1;
- } else {
- lambda = factor; // do something
- }
+ if (key == 1) {
+ lambda = Math.min(0.75 / (1 - factor), factor);
+ } else {
+ lambda = factor;
}
- lambda = (lambda > lmax) ? lmax : lambda;
+ lambda = Math.min(lambda, lmax);
- // interpolate new radii labels
- for (Vertex v : R1.keySet()) {
+ for (Vertex v : interiorVertices) {
+ int vid = vertexToId.get(v);
double r1 = R1.getDouble(v);
double r2 = R2.getDouble(v);
- double nwr = r2 + lambda * (r2 - r1);
- radii.put(v, nwr);
+ radiiArray[vid] = r2 + lambda * (r2 - r1);
}
- // end of superstep
-
- // do step/check superstep
- accumErr2 = computeAngleSums();
- accumErr2 = Math.sqrt(accumErr2);
-
- // check results
- double pred = FastMath.exp(lambda * FastMath.log(factor)); // predicted improvement
- double act = accumErr2 / c1; // actual improvement
- if (act < 1) { // did some good
- if (act > pred) { // not as good as expected: reset
- if (key == 1) {
- key = 2;
- }
- } // implied else: accept result
- } else { // reset to before superstep
- for (Vertex v : R1.keySet()) {
- double r2 = R2.getDouble(v);
- radii.put(v, r2);
- }
- accumErr2 = c1;
- if (key == 2) {
+ // NOTE probably faster to not call computeAngleSums() again. simply use sum
+ // from before
+// accumErr2 = computeAngleSums();
+// accumErr2 = Math.sqrt(accumErr2);
+ accumErr2 = c1;
+
+ localPasses++;
+ }
+ }
+
+ private void computeRadiiSuperStep() {
+ // Precompute tolerance based on problem size
+ final double ttoler = interiorVertices.size() <= 10 ? 1e-2 : Math.pow(10, -(3.5 + interiorVertices.size() / 100.0));
+
+ // Preallocate working arrays once
+ final double[] R1 = new double[radiiArray.length];
+ final double[] R2 = new double[radiiArray.length];
+
+ int key = 1;
+ double accumErr2 = Double.MAX_VALUE;
+ int localPasses = 1;
+ final int maxPasses = 3 * interiorVertices.size();
+
+ while (accumErr2 > ttoler && localPasses < maxPasses) {
+ // Phase 1: Standard iteration
+ System.arraycopy(radiiArray, 0, R1, 0, radiiArray.length);
+ double c1;
+ double factor;
+
+ // Single-pass factor calculation
+ do {
+ c1 = Math.sqrt(computeAngleSums());
+ factor = c1 / accumErr2;
+ if (factor >= 1.0) {
+ accumErr2 = c1;
key = 1;
}
+ } while (factor >= 1.0);
+
+ // Phase 2: Super-step preparation
+ System.arraycopy(radiiArray, 0, R2, 0, radiiArray.length);
+
+ // Lambda calculation with precomputed values
+ final double lambda = calculateLambda(R1, R2, factor, key);
+
+ // Vectorized radius update
+ updateRadii(R1, R2, lambda);
+
+ // Error calculation with early exit check
+ accumErr2 = Math.sqrt(computeAngleSums());
+
+ // Convergence monitoring
+ if (!updateState(R1, R2, c1, accumErr2, factor, key, lambda)) {
+ key = (key == 1) ? 2 : 1;
}
localPasses++;
}
}
+
+ private double calculateLambda(double[] R1, double[] R2, double factor, int key) {
+ double lmax = Double.MAX_VALUE;
+
+ // Parallel safe iteration (if needed)
+ for (int vid : interiorVertexIds) {
+ final double r1 = R1[vid];
+ final double r2 = R2[vid];
+ final double rat = r2 - r1;
+ if (rat < 0) {
+ lmax = Math.min(lmax, -r2 / rat);
+ }
+ }
+ lmax /= 2;
- /**
- * Determine the centers of the circles using radii of the interior circles.
- *
- * @return
- */
- private void computeCenters() {
- if (flowers.size() > 0) {
- placements = new HashMap<>();
-
- Vertex k1 = centralVertex; // pick one internal circle
- placements.put(k1, new Complex(0, 0)); // place it at the origin
+ if (key == 1) {
+ return Math.min(0.75 / (1 - factor), factor);
+ } else {
+ return Math.min(factor / (1 - factor), lmax);
+ }
+ }
- Vertex k2 = flowers.get(k1).get(0); // pick one of its neighbors
- placements.put(k2, new Complex(radii.getDouble(k1) + radii.getDouble(k2))); // place it on the real axis
- place(k1); // recursively place the rest
- place(k2);
+ private void updateRadii(double[] R1, double[] R2, double lambda) {
+ for (int vid : interiorVertexIds) {
+ final double delta = R2[vid] - R1[vid];
+ radiiArray[vid] = R2[vid] + lambda * delta;
}
+ }
+
+ private boolean updateState(double[] R1, double[] R2, double c1, double accumErr2, double factor, int key, double lambda) {
+ final double pred = FastMath.exp(lambda * FastMath.log(factor));
+ final double act = accumErr2 / c1;
- circles = new ArrayList<>(radii.size());
- placements.forEach(
- (v, pos) -> circles.add(new PVector((float) pos.getReal(), (float) pos.getImaginary(), (float) radii.getDouble(v))));
+ if (act >= 1) {
+ System.arraycopy(R2, 0, radiiArray, 0, radiiArray.length);
+ return false;
+ }
+ return act <= pred;
}
- /**
- * Compute the angle sum for every flower.
- *
- * @return sum of angle error (difference between 2PI) across all flowers
- */
private double computeAngleSums() {
double error = 0;
- for (Entry> entry : flowers.entrySet()) {
- final Vertex v = entry.getKey();
- final List flower = entry.getValue();
-
- final double ra = radii.getDouble(v);
- double angleSum = angleSum(ra, flower);
-
- final int N = 2 * flower.size();
- final double del = FastMath.sin(TWO_PI / N);
- final double bet = FastMath.sin(angleSum / N);
- final double r2 = ra * bet * (1 - del) / (del * (1 - bet));
+ for (int i = 0; i < interiorVertices.size(); i++) {
+ int vId = interiorVertexIds[i];
+ int[] flower = flowersIds[i];
+ double ra = radiiArray[vId];
+
+ // Compute angle sum for the flower
+ double angleSum = 0;
+ int n = flowerSizes[i];
+ // NOTE inlined angleSum
+ for (int j = 0; j < n; j++) {
+ int currentId = flower[j];
+ int nextId = (j + 1 < n) ? flower[j + 1] : flower[0];
+ double b = radiiArray[currentId];
+ double c = radiiArray[nextId];
+ double bc = b * c;
+ double denominator = ra * ra + ra * (b + c) + bc;
+ angleSum += FastMath.acos(1 - (2 * bc) / denominator);
+ }
- // alternative form
-// double hat = ra / (1.0 / FastMath.sin(angleSum / (2 * flower.size())) - 1);
-// double r2 = hat * (1.0 / FastMath.sin(Math.PI / flower.size()) - 1);
+ // Update radius using precomputed values
+ double factorDel = factorDelArray[i];
+ double bet = FastMath.sin(angleSum / (2 * n));
+ double r2 = ra * bet * factorDel / (1 - bet);
- radii.put(v, r2);
- angleSum -= TWO_PI;
- error += angleSum * angleSum; // accum abs error
+ radiiArray[vId] = r2;
+ error += (angleSum - TWO_PI) * (angleSum - TWO_PI);
}
return error;
}
- /**
- *
- * @param rc radius of center circle
- * @param center center circle
- * @param flower center circle's petals
- * @return
- */
- private double angleSum(final double rc, final List flower) {
+ private double angleSum(final double rc, final List flower) {
final int n = flower.size();
double sum = 0.0;
- for (int i = 0; i < n; i++) {
- int j = i + 1 == n ? 0 : i + 1;
- sum += tangentAngle(rc, radii.getDouble(flower.get(i)), radii.getDouble(flower.get(j)));
+ for (int j = 0; j < n; j++) {
+ int currentId = flower.get(j);
+ int nextId = flower.get((j + 1) % n);
+ sum += tangentAngle(rc, radiiArray[currentId], radiiArray[nextId]);
}
return sum;
}
- /**
- * Compute the angle sum for the petals surrounding the given vertex and update
- * the radius of the vertex such that the angle sum would equal 2π.
- *
- * @param center target vertex
- * @return a measure of the error (difference between target angle sum (2π) and
- * the actual angle sum
- * @deprecated used by {@link #computeRadiiSimple()}
- */
- @Deprecated
- private double flower(final Vertex center) {
- List flower = flowers.get(center);
- final int n = flower.size();
- final double rc = radii.getDouble(center);
- double sum = 0.0;
- for (int i = 0; i < n; i++) {
- int j = i + 1 == n ? 0 : i + 1;
- sum += tangentAngleFast(rc, radii.getDouble(flower.get(i)), radii.getDouble(flower.get(j)));
+ private void computeCenters() {
+ if (interiorVertices.isEmpty()) {
+ circles = new ArrayList<>();
+ return;
}
- double hat = rc / (1.0 / FastMath.sin(sum / (2 * n)) - 1);
- double newrad = hat * (1.0 / FastMath.sin(Math.PI / n) - 1);
- radii.put(center, newrad);
+ placementsArray = new Complex[allVertices.size()];
+ int centralId = vertexToId.get(centralVertex);
+ placementsArray[centralId] = new Complex(0, 0);
- return Math.max(newrad / rc, rc / newrad);
- }
+ int centralIndex = interiorVertexIndex.get(centralVertex);
+ int[] centralFlower = flowersIds[centralIndex];
+ int k2Id = centralFlower[0];
+ placementsArray[k2Id] = new Complex(radiiArray[centralId] + radiiArray[k2Id], 0);
- /**
- * Recursively determine centers of all circles surrounding a given vertex.
- */
- private void place(final Vertex centre) {
+ place(centralVertex);
+ place(allVertices.get(k2Id));
- if (!flowers.containsKey(centre)) {
- return; // boundary vertex
+ circles = new ArrayList<>(allVertices.size());
+ for (Vertex v : allVertices) {
+ int id = vertexToId.get(v);
+ Complex pos = placementsArray[id];
+ if (pos != null) {
+ circles.add(new PVector((float) pos.getReal(), (float) pos.getImaginary(), (float) radiiArray[id]));
+ }
}
+ }
- List flower = flowers.get(centre);
- final int nc = flower.size();
- final double rcentre = radii.getDouble(centre);
-
- final Complex minusI = new Complex(0.0, -1);
+ private void place(Vertex centre) {
+ int centreId = vertexToId.get(centre);
+ if (!interiorVertexIndex.containsKey(centre)) {
+ return;
+ }
- for (int i = -nc; i < nc - 1; i++) {
- int ks = i < 0 ? nc + i : i;
- Vertex s = flower.get(ks);
- double rs = radii.getDouble(s);
+ int centreIndex = interiorVertexIndex.get(centre);
+ int[] flower = flowersIds[centreIndex];
+ int nc = flower.length;
+ double rcentre = radiiArray[centreId];
+ Complex minusI = new Complex(0, -1);
- int kt = ks + 1 < nc ? ks + 1 : 0;
- Vertex t = flower.get(kt);
- double rt = radii.getDouble(t);
+ for (int i = 0; i < nc; i++) {
+ int sId = flower[i];
+ int tId = flower[(i + 1) % nc];
- if (placements.containsKey(s) && !placements.containsKey(t)) {
+ if (placementsArray[sId] != null && placementsArray[tId] == null) {
+ double rs = radiiArray[sId];
+ double rt = radiiArray[tId];
double theta = tangentAngle(rcentre, rs, rt);
- Complex offset = (placements.get(s).subtract(placements.get(centre))).divide(new Complex(rs + rcentre));
+
+ Complex offset = placementsArray[sId].subtract(placementsArray[centreId]).divide(new Complex(rs + rcentre));
offset = offset.multiply(minusI.multiply(theta).exp());
- placements.put(t, placements.get(centre).add(offset.multiply(rt + rcentre)));
+ placementsArray[tId] = placementsArray[centreId].add(offset.multiply(rt + rcentre));
- place(t);
+ place(allVertices.get(tId));
}
}
}
private static double tangentAngle(double a, double b, double c) {
- /*
- * Overall computation time is actually reduced by foregoing trig approximation
- * functions (such as tangentAngleFast()), because the slight inaccuracies cause
- * solution to converge more slowly, performing more iterations overall.
- */
- final double q = b * c;
- final double o = 1 - 2 * q / (a * a + a * (b + c) + q);
+ double q = b * c;
+ double o = 1 - 2 * q / (a * a + a * (b + c) + q);
return FastMath.acos(o);
}
- /**
- * Computes the angle that circles y and z make with circle x (angle yxz). The
- * circles are given by their radii and are mutually tangent.
- *
- * @param rx radius of circle x, the circle of interest
- * @param ry radius of circle y, a "petal" circle
- * @param rz radius of circle z, a "petal" circle
- * @return angle of yxz
- * @deprecated
- */
- @Deprecated
- private static double tangentAngleFast(final double rx, final double ry, final double rz) {
- final double x = (ry * rz) / ((rx + ry) * (rx + rz));
- // return 2 * Fixed64.ToDouble(Fixed64.Asin(Fixed64.FromDouble(Math.sqrt(x))));
- // return 2 * Fixed64.ToDouble(Fixed64.Atan(Fixed64.FromDouble(Math.sqrt(x) / (Math.sqrt(1 - x)))));
- return 0;
- }
-
private static class RadialComparator implements Comparator {
-
private Vertex origin;
public RadialComparator(Vertex origin) {
@@ -475,34 +453,15 @@ public int compare(Vertex o1, Vertex o2) {
return polarCompare(origin, o1, o2);
}
- /**
- * Given two points p and q compare them with respect to their radial ordering
- * about point o. First checks radial ordering.
- *
- * @param o the origin
- * @param p a point
- * @param q another point
- * @return -1, 0 or 1 depending on whether angle p is less than, equal to or
- * greater than angle q
- */
private static int polarCompare(Vertex o, Vertex p, Vertex q) {
double dxp = p.x - o.x;
double dyp = p.y - o.y;
double dxq = q.x - o.x;
double dyq = q.y - o.y;
- int result = 0;
- double alph = FastAtan2.atan2(dxp, dyp);
- double beta = FastAtan2.atan2(dxq, dyq);
- if (alph < beta) {
- result = -1;
- }
- if (alph > beta) {
- result = 1;
- }
- return result;
+ double alph = FastMath.atan2(dyp, dxp);
+ double beta = FastMath.atan2(dyq, dxq);
+ return Double.compare(alph, beta);
}
-
}
-
}
diff --git a/src/test/java/micycle/pgs/PGSTests.java b/src/test/java/micycle/pgs/PGSTests.java
index b08f13c2..19591213 100644
--- a/src/test/java/micycle/pgs/PGSTests.java
+++ b/src/test/java/micycle/pgs/PGSTests.java
@@ -1,6 +1,7 @@
package micycle.pgs;
import static org.junit.jupiter.api.Assertions.assertEquals;
+import static org.junit.jupiter.api.Assertions.assertTrue;
import java.util.ArrayList;
import java.util.Arrays;
@@ -32,7 +33,7 @@ void testFromEdges() {
for (int i = 0; i < 15; i++) {
edges.add(new PEdge(i, i, i + 1, i + 1));
}
- edges.add(new PEdge(15, 15, 0,0)); // close
+ edges.add(new PEdge(15, 15, 0, 0)); // close
Collections.shuffle(edges);
@@ -41,4 +42,35 @@ void testFromEdges() {
assertEquals(16, orderedVertices.size());
}
+ @Test
+ void testOrientation() {
+ /*
+ * NOTE the isClockwise method tests for orientation in a y-axis-up coordinate
+ * system. The lists below are defined in terms of visual orientation in
+ * Processing (which uses y-axis-down orientation). Hence the results are
+ * inverted to check the method is geometrically correct.
+ */
+
+ // @formatter:off
+ // (0,0) --------> (1,0) (x-axis increasing to the right)
+ // | |
+ // | |
+ // V V
+ // (0,1) <--------- (1,1) (y-axis increasing downwards)
+ // @formatter:on
+ List clockwisePoints = List.of(new PVector(0, 0), new PVector(1, 0), new PVector(1, 1), new PVector(0, 1));
+ assertTrue(!PGS.isClockwise(clockwisePoints)); // NOTE inverted
+
+ List counterClockwisePoints = List.of(new PVector(0, 0), new PVector(0, 1), new PVector(1, 1), new PVector(1, 0));
+ assertTrue(PGS.isClockwise(counterClockwisePoints)); // NOTE inverted
+
+ List clockwisePointsClosed = new ArrayList<>(
+ List.of(new PVector(0, 0), new PVector(1, 0), new PVector(1, 1), new PVector(0, 1), new PVector(0, 0)));
+ assertTrue(!PGS.isClockwise(clockwisePointsClosed)); // NOTE inverted
+
+ List counterClockwisePointsClosed = new ArrayList<>(
+ List.of(new PVector(0, 0), new PVector(0, 1), new PVector(1, 1), new PVector(1, 0), new PVector(0, 0)));
+ assertTrue(PGS.isClockwise(counterClockwisePointsClosed)); // NOTE inverted
+ }
+
}
diff --git a/src/test/java/micycle/pgs/PGS_ConversionTests.java b/src/test/java/micycle/pgs/PGS_ConversionTests.java
index 169fb201..4c2decaa 100644
--- a/src/test/java/micycle/pgs/PGS_ConversionTests.java
+++ b/src/test/java/micycle/pgs/PGS_ConversionTests.java
@@ -352,14 +352,14 @@ void testMultiContour() {
@Test
void testVertexRounding() {
- final PShape shape = new PShape(PShape.GEOMETRY);
+ PShape shape = new PShape(PShape.GEOMETRY);
shape.beginShape();
shape.vertex(12.4985f, -97.234f);
shape.vertex(10, -10);
shape.vertex(999.99f, 0.0001f);
shape.endShape(PConstants.CLOSE); // close affects rendering only -- does not append another vertex
- PGS_Conversion.roundVertexCoords(shape);
+ shape = PGS_Conversion.roundVertexCoords(shape);
assertEquals(12, shape.getVertex(0).x);
assertEquals(-97, shape.getVertex(0).y);
@@ -368,6 +368,25 @@ void testVertexRounding() {
assertEquals(1000, shape.getVertex(2).x);
assertEquals(0, shape.getVertex(2).y);
}
+
+ @Test
+ void testVertexRounding1DP() {
+ PShape shape = new PShape(PShape.GEOMETRY);
+ shape.beginShape();
+ shape.vertex(12.4985f, -97.234f);
+ shape.vertex(10, -10);
+ shape.vertex(999.34f, 0.049f);
+ shape.endShape(PConstants.CLOSE);
+
+ shape = PGS_Conversion.roundVertexCoords(shape, 1);
+
+ assertEquals(12.5, shape.getVertex(0).x, 1e-5);
+ assertEquals(-97.2, shape.getVertex(0).y, 1e-5);
+ assertEquals(10, shape.getVertex(1).x, 1e-5);
+ assertEquals(-10, shape.getVertex(1).y, 1e-5);
+ assertEquals(999.3, shape.getVertex(2).x, 2e-5);
+ assertEquals(0, shape.getVertex(2).y, 1e-5);
+ }
@Test
void testDuplicateVertices() {
diff --git a/src/test/java/micycle/pgs/PGS_ShapePredicatesTests.java b/src/test/java/micycle/pgs/PGS_ShapePredicatesTests.java
index 1278f37a..b7fa57e4 100644
--- a/src/test/java/micycle/pgs/PGS_ShapePredicatesTests.java
+++ b/src/test/java/micycle/pgs/PGS_ShapePredicatesTests.java
@@ -18,7 +18,7 @@ class PGS_ShapePredicatesTests {
private static final double EPSILON = 1E-4;
- static PShape square, triangle;
+ static PShape square, triangle, rect;
@BeforeAll
static void initShapes() {
@@ -29,6 +29,14 @@ static void initShapes() {
square.vertex(10, 10);
square.vertex(0, 10);
square.endShape(PConstants.CLOSE); // close affects rendering only -- does not append another vertex
+
+ rect = new PShape(PShape.GEOMETRY); // 10x10 rect
+ rect.beginShape();
+ rect.vertex(0, 0);
+ rect.vertex(10, 0);
+ rect.vertex(10, 20);
+ rect.vertex(0, 20);
+ rect.endShape(PConstants.CLOSE); // close affects rendering only -- does not append another vertex
float[] centroid = new float[] { 0, 0 };
float side_length = 10;
@@ -63,6 +71,37 @@ void testMaximumInteriorAngle() {
assertEquals(Math.PI / 3, PGS_ShapePredicates.maximumInteriorAngle(triangle), EPSILON);
}
+ @Test
+ void testInteriorAnglesSquare() {
+ var angles = PGS_ShapePredicates.interiorAngles(square);
+ assertEquals(4, angles.size(), "Square should have 4 angles");
+
+ double expectedAngleRadians = Math.PI / 2.0; // 90 degrees in radians
+ double expectedAngleSumRadians = Math.PI * 2; // 360 degrees for a square
+ double actualAngleSumRadians = 0;
+
+ for (double angle : angles.values()) {
+ assertEquals(expectedAngleRadians, angle, 1e-6, "Interior angle should be approximately 90 degrees");
+ actualAngleSumRadians += angle;
+ }
+ assertEquals(expectedAngleSumRadians, actualAngleSumRadians, 1e-6, "Sum of square interior angles should be approximately 360 degrees");
+
+ }
+
+ @Test
+ void testInteriorAnglesTriangle() {
+ var angles = PGS_ShapePredicates.interiorAngles(triangle);
+
+ assertEquals(3, angles.size(), "Triangle should have 3 angles");
+
+ double expectedAngleSumRadians = Math.PI; // 180 degrees for a triangle
+ double actualAngleSumRadians = 0;
+ for (double angle : angles.values()) {
+ actualAngleSumRadians += angle;
+ }
+ assertEquals(expectedAngleSumRadians, actualAngleSumRadians, 1e-6, "Sum of triangle interior angles should be approximately 180 degrees");
+ }
+
@Test
void testHoles() {
assertEquals(0, PGS_ShapePredicates.holes(square));
@@ -77,15 +116,21 @@ void testHoles() {
coverage.removeChild(0); // remove a mesh face; mesh no longer forms a hole
assertEquals(0, PGS_ShapePredicates.holes(coverage));
}
-
+
@Test
void testIsClockwise() {
assertTrue(PGS_ShapePredicates.isClockwise(square));
- List ccw = PGS_Conversion.toPVector(square);//.reversed();
+ List ccw = PGS_Conversion.toPVector(square);// .reversed();
Collections.reverse(ccw);
ccw.add(ccw.get(0)); // close
assertFalse(PGS_ShapePredicates.isClockwise(PGS_Conversion.fromPVector(ccw)));
-
+
+ }
+
+ @Test
+ void testElongation() {
+ assertEquals(0, PGS_ShapePredicates.elongation(square));
+ assertEquals(0.5, PGS_ShapePredicates.elongation(rect));
}
}