diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java index 14a565639..3d914a3e9 100644 --- a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java @@ -211,7 +211,7 @@ public class SphericalPolygonsSet extends AbstractRegion { } /** Recursively build a tree by inserting cut sub-hyperplanes. - * @param hyperplaneThickness tolerance below which points are consider to + * @param hyperplaneThickness tolerance below which points are considered to * belong to the hyperplane (which is therefore more a slab) * @param node current tree node (it is a leaf node at the beginning * of the call) @@ -264,12 +264,12 @@ public class SphericalPolygonsSet extends AbstractRegion { if (!outsideList.isEmpty()) { insertEdges(hyperplaneThickness, node.getPlus(), outsideList); } else { - node.getMinus().setAttribute(Boolean.FALSE); + node.getPlus().setAttribute(Boolean.FALSE); } if (!insideList.isEmpty()) { insertEdges(hyperplaneThickness, node.getMinus(), insideList); } else { - node.getPlus().setAttribute(Boolean.TRUE); + node.getMinus().setAttribute(Boolean.TRUE); } } @@ -494,7 +494,7 @@ public class SphericalPolygonsSet extends AbstractRegion { private void split(final Circle splitCircle, final List outsideList, final List insideList) { - // get the inside chord, synchronizing its phase with the edge itself + // get the inside arc, synchronizing its phase with the edge itself final double edgeStart = circle.getPhase(start.getLocation().getVector()); final double arcRelativeStart = MathUtils.normalizeAngle(circle.getInsideArc(splitCircle).getInf(), edgeStart + FastMath.PI) - edgeStart; @@ -502,80 +502,93 @@ public class SphericalPolygonsSet extends AbstractRegion { final double unwrappedEnd = arcRelativeStart - FastMath.PI; // build the sub-edges - if (arcRelativeStart < length) { - if (unwrappedEnd > 0) { - // the edge starts inside the circle, then goes outside, then comes back inside + final double tolerance = circle.getTolerance(); + Vertex previousVertex = start; + if (unwrappedEnd >= length - tolerance) { - // create intermediate vertices - final Vertex vExit = new Vertex(new S2Point(circle.getPointAt(edgeStart + arcRelativeEnd))); - final Vertex vEnter = new Vertex(new S2Point(circle.getPointAt(edgeStart + arcRelativeStart))); - vExit.bindWith(splitCircle); - vEnter.bindWith(splitCircle); + // the edge is entirely contained inside the circle + // we don't split anything + insideList.add(this); - // create sub-edges - final Edge eStartIn = new Edge(start, vExit, unwrappedEnd, circle); - final Edge eMiddleOut = new Edge(vExit, vEnter, arcRelativeStart - unwrappedEnd, circle); - final Edge eEndIn = new Edge(vEnter, end, length - arcRelativeStart, circle); - eStartIn.setNode(node); - eMiddleOut.setNode(node); - eEndIn.setNode(node); - - // distribute the sub-edges in the appropriate lists - insideList.add(eStartIn); - insideList.add(eEndIn); - outsideList.add(eMiddleOut); - - } else { - // the edge starts outside of the circle, then comes inside - - // create intermediate vertices - final Vertex vEnter = new Vertex(new S2Point(circle.getPointAt(edgeStart + arcRelativeStart))); - vEnter.bindWith(splitCircle); - - // create sub-edges - final Edge eStartOut = new Edge(start, vEnter, arcRelativeStart, circle); - final Edge eEndIn = new Edge(vEnter, end, length - arcRelativeStart, circle); - eStartOut.setNode(node); - eEndIn.setNode(node); - - // distribute the sub-edges in the appropriate lists - outsideList.add(eStartOut); - insideList.add(eEndIn); - - } } else { - if (unwrappedEnd > 0) { - if (unwrappedEnd > length) { - // the edge is entirely contained inside the circle - // we don't split anything - insideList.add(this); + + // there are at least some parts of the edge that should be outside + // (even is they are later be filtered out as being too small) + double alreadyManagedLength = 0; + if (unwrappedEnd >= 0) { + // the start of the edge is inside the circle + previousVertex = addSubEdge(previousVertex, + new Vertex(new S2Point(circle.getPointAt(edgeStart + unwrappedEnd))), + unwrappedEnd, insideList, splitCircle); + alreadyManagedLength = unwrappedEnd; + } + + if (arcRelativeStart >= length - tolerance) { + // the edge ends while still outside of the circle + if (unwrappedEnd >= 0) { + previousVertex = addSubEdge(previousVertex, end, + length - alreadyManagedLength, outsideList, splitCircle); } else { - // the edge starts inside the circle, then goes outside - - // create intermediate vertices - final Vertex vExit = new Vertex(new S2Point(circle.getPointAt(edgeStart + arcRelativeEnd))); - vExit.bindWith(splitCircle); - - // create sub-edges - final Edge eStartIn = new Edge(start, vExit, arcRelativeEnd, circle); - final Edge eEndOut = new Edge(vExit, end, length - arcRelativeEnd, circle); - eStartIn.setNode(node); - eEndOut.setNode(node); - - // distribute the sub-edges in the appropriate lists - insideList.add(eStartIn); - outsideList.add(eEndOut); - + // the edge is entirely outside of the circle + // we don't split anything + outsideList.add(this); } } else { - // the edge is entirely outside of the circle - // we don't split anything - outsideList.add(this); + // the edge is long enough to enter inside the circle + previousVertex = addSubEdge(previousVertex, + new Vertex(new S2Point(circle.getPointAt(edgeStart + arcRelativeStart))), + arcRelativeStart - alreadyManagedLength, outsideList, splitCircle); + alreadyManagedLength = arcRelativeStart; + + if (arcRelativeEnd >= length - tolerance) { + // the edge ends while still inside of the circle + previousVertex = addSubEdge(previousVertex, end, + length - alreadyManagedLength, insideList, splitCircle); + } else { + // the edge is long enough to exit outside of the circle + previousVertex = addSubEdge(previousVertex, + new Vertex(new S2Point(circle.getPointAt(edgeStart + arcRelativeStart))), + arcRelativeStart - alreadyManagedLength, insideList, splitCircle); + alreadyManagedLength = arcRelativeStart; + previousVertex = addSubEdge(previousVertex, end, + length - alreadyManagedLength, outsideList, splitCircle); + } } + } } + /** Add a sub-edge to a list if long enough. + *

+ * If the length of the sub-edge to add is smaller than the {@link Circle#getTolerance()} + * tolerance of the support circle, it will be ignored. + *

+ * @param subStart start of the sub-edge + * @param subEnd end of the sub-edge + * @param subLength length of the sub-edge + * @param splitCircle circle splitting the edge in several parts + * @param list list where to put the sub-edge + * @return end vertex of the edge ({@code subEnd} if the edge was long enough and really + * added, {@code subStart} if the edge was too small and therefore ignored) + */ + private Vertex addSubEdge(final Vertex subStart, final Vertex subEnd, final double subLength, + final List list, final Circle splitCircle) { + + if (subLength <= circle.getTolerance()) { + // the edge is too short, we ignore it + return subStart; + } + + // really add the edge + subEnd.bindWith(splitCircle); + final Edge edge = new Edge(subStart, subEnd, subLength, circle); + edge.setNode(node); + list.add(edge); + return subEnd; + + } + } /** {@inheritDoc} */ diff --git a/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java b/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java index 7fd922039..b8f7a0174 100644 --- a/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java +++ b/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java @@ -18,6 +18,7 @@ package org.apache.commons.math3.geometry.spherical.twod; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.apache.commons.math3.geometry.partitioning.Region.Location; +import org.apache.commons.math3.geometry.partitioning.RegionFactory; import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; import org.apache.commons.math3.random.Well1024a; import org.apache.commons.math3.util.FastMath; @@ -35,6 +36,8 @@ public class SphericalPolygonsSetTest { Vector3D v = new Vector3D(random.nextVector()); Assert.assertEquals(Location.INSIDE, full.checkPoint(new S2Point(v))); } + Assert.assertEquals(4 * FastMath.PI, new SphericalPolygonsSet(0.01, new S2Point[0]).getSize(), 1.0e-10); + Assert.assertEquals(0, new SphericalPolygonsSet(0.01, new S2Point[0]).getBoundarySize(), 1.0e-10); } @Test @@ -43,7 +46,7 @@ public class SphericalPolygonsSetTest { double sinTol = FastMath.sin(tol); SphericalPolygonsSet south = new SphericalPolygonsSet(Vector3D.MINUS_K, tol); UnitSphereRandomVectorGenerator random = - new UnitSphereRandomVectorGenerator(3, new Well1024a(0x852fd2a0ed8d2f6dl)); + new UnitSphereRandomVectorGenerator(3, new Well1024a(0x6b9d4a6ad90d7b0bl)); for (int i = 0; i < 1000; ++i) { Vector3D v = new Vector3D(random.nextVector()); if (v.getZ() < -sinTol) { @@ -51,8 +54,50 @@ public class SphericalPolygonsSetTest { } else if (v.getZ() > sinTol) { Assert.assertEquals(Location.OUTSIDE, south.checkPoint(new S2Point(v))); } else { - Assert.assertEquals("" + v.getX() + " " + v.getY() + " " + v.getZ(), - Location.BOUNDARY, south.checkPoint(new S2Point(v))); + Assert.assertEquals(Location.BOUNDARY, south.checkPoint(new S2Point(v))); + } + } + } + + @Test + public void testPositiveOctantByIntersection() { + double tol = 0.01; + double sinTol = FastMath.sin(tol); + RegionFactory factory = new RegionFactory(); + SphericalPolygonsSet plusX = new SphericalPolygonsSet(Vector3D.PLUS_I, tol); + SphericalPolygonsSet plusY = new SphericalPolygonsSet(Vector3D.PLUS_J, tol); + SphericalPolygonsSet plusZ = new SphericalPolygonsSet(Vector3D.PLUS_K, tol); + SphericalPolygonsSet octant = + (SphericalPolygonsSet) factory.intersection(factory.intersection(plusX, plusY), plusZ); + UnitSphereRandomVectorGenerator random = + new UnitSphereRandomVectorGenerator(3, new Well1024a(0x9c9802fde3cbcf25l)); + for (int i = 0; i < 1000; ++i) { + Vector3D v = new Vector3D(random.nextVector()); + if ((v.getX() > sinTol) && (v.getY() > sinTol) && (v.getZ() > sinTol)) { + Assert.assertEquals(Location.INSIDE, octant.checkPoint(new S2Point(v))); + } else if ((v.getX() < -sinTol) || (v.getY() < -sinTol) || (v.getZ() < -sinTol)) { + Assert.assertEquals(Location.OUTSIDE, octant.checkPoint(new S2Point(v))); + } else { + Assert.assertEquals(Location.BOUNDARY, octant.checkPoint(new S2Point(v))); + } + } + } + + @Test + public void testPositiveOctantByVertices() { + double tol = 0.01; + double sinTol = FastMath.sin(tol); + SphericalPolygonsSet octant = new SphericalPolygonsSet(tol, S2Point.PLUS_I, S2Point.PLUS_J, S2Point.PLUS_K); + UnitSphereRandomVectorGenerator random = + new UnitSphereRandomVectorGenerator(3, new Well1024a(0xb8fc5acc91044308l)); + for (int i = 0; i < 1000; ++i) { + Vector3D v = new Vector3D(random.nextVector()); + if ((v.getX() > sinTol) && (v.getY() > sinTol) && (v.getZ() > sinTol)) { + Assert.assertEquals(Location.INSIDE, octant.checkPoint(new S2Point(v))); + } else if ((v.getX() < -sinTol) || (v.getY() < -sinTol) || (v.getZ() < -sinTol)) { + Assert.assertEquals(Location.OUTSIDE, octant.checkPoint(new S2Point(v))); + } else { + Assert.assertEquals(Location.BOUNDARY, octant.checkPoint(new S2Point(v))); } } }