From 9b3a46e20540dd07cef3e996cc498e5fd6bd4996 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Tue, 14 Jan 2014 18:26:32 +0000 Subject: [PATCH] Fixed problem for spherical polygons with holes. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1558149 13f79535-47bb-0310-9956-ffa450edef68 --- .../math3/geometry/partitioning/BSPTree.java | 40 +++ .../geometry/spherical/oned/ArcsSet.java | 4 +- .../spherical/twod/SphericalPolygonsSet.java | 237 ++++++++++++++---- .../twod/SphericalPolygonsSetTest.java | 201 ++++++++++++++- 4 files changed, 422 insertions(+), 60 deletions(-) diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java index 040092a92..cc526f2d7 100644 --- a/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java @@ -665,6 +665,46 @@ public class BSPTree { } + /** Prune a tree around a cell. + *

+ * This method can be used to extract a convex cell from a tree. + * The original cell may either be a leaf node or an internal node. + * If it is an internal node, it's subtree will be ignored (i.e. the + * extracted cell will be a leaf node in all cases). The original + * tree to which the original cell belongs is not touched at all, + * a new independent tree will be built. + *

+ * @param cellAttribute attribute to set for the leaf node + * corresponding to the initial instance cell + * @param otherLeafsAttributes attribute to set for the other leaf + * nodes + * @param internalAttributes attribute to set for the internal nodes + * @return a new tree (the original tree is left untouched) containing + * a single branch with the cell as a leaf node, and other leaf nodes + * as the remnants of the pruned branches + */ + public BSPTree pruneAroundConvexCell(final Object cellAttribute, + final Object otherLeafsAttributes, + final Object internalAttributes) { + + // build the current cell leaf + BSPTree tree = new BSPTree(cellAttribute); + + // build the pruned tree bottom-up + for (BSPTree current = this; current.parent != null; current = current.parent) { + final SubHyperplane parentCut = current.parent.cut.copySelf(); + final BSPTree sibling = new BSPTree(otherLeafsAttributes); + if (current == current.parent.plus) { + tree = new BSPTree(parentCut, tree, sibling, internalAttributes); + } else { + tree = new BSPTree(parentCut, sibling, tree, internalAttributes); + } + } + + return tree; + + } + /** Chop off parts of the tree. *

The instance is modified in place, all the parts that are on * the minus side of the chopping hyperplane are discarded, only the diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java index 923d5a786..326bc1b5d 100644 --- a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java @@ -520,8 +520,8 @@ public class ArcsSet extends AbstractRegion implements Itera current = firstStart; if (firstStart == null) { - // the tree has a single node - if ((Boolean) getTree(false).getAttribute()) { + // all the leaf tree nodes share the same inside/outside status + if ((Boolean) getFirstLeaf(getTree(false)).getAttribute()) { // it is an inside node, it represents the full circle pending = new double[] { 0, MathUtils.TWO_PI 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 cf7dfc24c..2284c6601 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 @@ -25,7 +25,9 @@ import java.util.List; import java.util.Map; import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.MathInternalError; import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.geometry.euclidean.threed.Rotation; import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; import org.apache.commons.math3.geometry.partitioning.AbstractRegion; import org.apache.commons.math3.geometry.partitioning.BSPTree; @@ -67,6 +69,19 @@ public class SphericalPolygonsSet extends AbstractRegion { tolerance); } + /** Build a polygons set representing a regular polygon. + * @param center center of the polygon (the center is in the inside half) + * @param meridian point defining the reference meridian for first polygon vertex + * @param outsideRadius distance of the vertices to the center + * @param n number of sides of the polygon + * @param tolerance below which points are consider to be identical + */ + public SphericalPolygonsSet(final Vector3D center, final Vector3D meridian, + final double outsideRadius, final int n, + final double tolerance) { + this(tolerance, createRegularPolygonVertices(center, meridian, outsideRadius, n)); + } + /** Build a polygons set from a BSP tree. *

The leaf nodes of the BSP tree must have a * {@code Boolean} attribute representing the inside status of @@ -140,6 +155,27 @@ public class SphericalPolygonsSet extends AbstractRegion { super(verticesToTree(hyperplaneThickness, vertices), hyperplaneThickness); } + /** Build the vertices representing a regular polygon. + * @param center center of the polygon (the center is in the inside half) + * @param meridian point defining the reference meridian for first polygon vertex + * @param outsideRadius distance of the vertices to the center + * @param n number of sides of the polygon + * @return vertices array + */ + private static S2Point[] createRegularPolygonVertices(final Vector3D center, final Vector3D meridian, + final double outsideRadius, final int n) { + final S2Point[] array = new S2Point[n]; + final Rotation r0 = new Rotation(Vector3D.crossProduct(center, meridian), outsideRadius); + array[0] = new S2Point(r0.applyTo(center)); + + final Rotation r = new Rotation(center, MathUtils.TWO_PI / n); + for (int i = 1; i < n; ++i) { + array[i] = new S2Point(r.applyTo(array[i - 1].getVector())); + } + + return array; + } + /** Build the BSP tree of a polygons set from a simple list of vertices. *

The boundary is provided as a list of points considering to * represent the vertices of a simple loop. The interior part of the @@ -583,70 +619,28 @@ public class SphericalPolygonsSet extends AbstractRegion { @Override protected void computeGeometricalProperties() throws MathIllegalStateException { - final List boundary = getBoundaryLoops(); + final BSPTree tree = getTree(true); + + if (tree.getCut() == null) { + + // the instance has a single cell without any boundaries - if (boundary.isEmpty()) { - final BSPTree tree = getTree(false); if (tree.getCut() == null && (Boolean) tree.getAttribute()) { // the instance covers the whole space setSize(4 * FastMath.PI); + setBarycenter(new S2Point(0, 0)); } else { setSize(0); + setBarycenter(S2Point.NaN); } - setBarycenter(new S2Point(0, 0)); + } else { - // compute some integrals around the shape - double sumArea = 0; - Vector3D sumB = Vector3D.ZERO; - - for (final Vertex startVertex : boundary) { - - int n = 0; - double sum = 0; - Vector3D sumP = Vector3D.ZERO; - for (Edge edge = startVertex.getOutgoing(); - n == 0 || edge.getStart() != startVertex; - edge = edge.getEnd().getOutgoing()) { - final Vector3D middle = edge.getPointAt(0.5 * edge.getLength()); - sumP = new Vector3D(1, sumP, edge.getLength(), middle); - - // find path interior angle at vertex - final Vector3D previousPole = edge.getCircle().getPole(); - final Vector3D nextPole = edge.getEnd().getOutgoing().getCircle().getPole(); - final Vector3D point = edge.getEnd().getLocation().getVector(); - double alpha = FastMath.atan2(Vector3D.dotProduct(nextPole, Vector3D.crossProduct(point, previousPole)), - -Vector3D.dotProduct(nextPole, previousPole)); - if (alpha < 0) { - alpha += MathUtils.TWO_PI; - } - sum += alpha; - n++; - } - - // compute area using extended Girard theorem - // see Spherical Trigonometry: For the Use of Colleges and Schools by I. Todhunter - // article 99 in chapter VIII Area Of a Spherical Triangle. Spherical Excess. - // book available from project Gutenberg at http://www.gutenberg.org/ebooks/19770 - final double area = sum - (n - 2) * FastMath.PI; - sumArea += area; - - sumB = new Vector3D(1, sumB, area, sumP); - - } - - if (sumArea < 0) { - sumArea = 4 * FastMath.PI - sumArea; - sumB = sumB.negate(); - } - - setSize(sumArea); - final double norm = sumB.getNorm(); - if (norm == 0.0) { - setBarycenter(S2Point.NaN); - } else { - setBarycenter(new S2Point(new Vector3D(1.0 / norm, sumB))); - } + // the instance has a boundary + final PropertiesComputer pc = new PropertiesComputer(getTolerance()); + tree.visit(pc); + setSize(pc.getArea()); + setBarycenter(pc.getBarycenter()); } @@ -849,4 +843,137 @@ public class SphericalPolygonsSet extends AbstractRegion { } + /** Visitor computing geometrical properties. */ + private static class PropertiesComputer implements BSPTreeVisitor { + + /** Tolerance below which points are consider to be identical. */ + private final double tolerance; + + /** Summed area. */ + private double summedArea; + + /** Summed barycenter. */ + private Vector3D summedBarycenter; + + /** Simple constructor. + * @param tolerance below which points are consider to be identical + */ + public PropertiesComputer(final double tolerance) { + this.tolerance = tolerance; + this.summedArea = 0; + this.summedBarycenter = Vector3D.ZERO; + } + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree node) { + return Order.MINUS_SUB_PLUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree node) { + // nothing to do here + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree node) { + if ((Boolean) node.getAttribute()) { + + // transform this inside leaf cell into a simple convex polygon + final SphericalPolygonsSet convex = + new SphericalPolygonsSet(node.pruneAroundConvexCell(Boolean.TRUE, + Boolean.FALSE, + null), + tolerance); + + // extract the start of the single loop boundary of the convex cell + final List boundary = convex.getBoundaryLoops(); + if (boundary.size() != 1) { + // this should never happen + throw new MathInternalError(); + } + + // compute the geometrical properties of the convex cell + final double area = convexCellArea(boundary.get(0)); + final Vector3D barycenter = convexCellBarycenter(boundary.get(0)); + + // add the cell contribution to the global properties + summedArea += area; + summedBarycenter = new Vector3D(1, summedBarycenter, area, barycenter); + + } + } + + /** Compute convex cell area. + * @param start start vertex of the convex cell boundary + * @return area + */ + private double convexCellArea(final Vertex start) { + + int n = 0; + double sum = 0; + + // loop around the cell + for (Edge e = start.getOutgoing(); n == 0 || e.getStart() != start; e = e.getEnd().getOutgoing()) { + + // find path interior angle at vertex + final Vector3D previousPole = e.getCircle().getPole(); + final Vector3D nextPole = e.getEnd().getOutgoing().getCircle().getPole(); + final Vector3D point = e.getEnd().getLocation().getVector(); + double alpha = FastMath.atan2(Vector3D.dotProduct(nextPole, Vector3D.crossProduct(point, previousPole)), + -Vector3D.dotProduct(nextPole, previousPole)); + if (alpha < 0) { + alpha += MathUtils.TWO_PI; + } + sum += alpha; + n++; + } + + // compute area using extended Girard theorem + // see Spherical Trigonometry: For the Use of Colleges and Schools by I. Todhunter + // article 99 in chapter VIII Area Of a Spherical Triangle. Spherical Excess. + // book available from project Gutenberg at http://www.gutenberg.org/ebooks/19770 + return sum - (n - 2) * FastMath.PI; + + } + + /** Compute convex cell barycenter. + * @param start start vertex of the convex cell boundary + * @return barycenter + */ + private Vector3D convexCellBarycenter(final Vertex start) { + + int n = 0; + Vector3D sumB = Vector3D.ZERO; + + // loop around the cell + for (Edge e = start.getOutgoing(); n == 0 || e.getStart() != start; e = e.getEnd().getOutgoing()) { + final Vector3D middle = e.getPointAt(0.5 * e.getLength()); + sumB = new Vector3D(1, sumB, e.getLength(), middle); + n++; + } + + return sumB.normalize(); + + } + + /** Get the area. + * @return area + */ + public double getArea() { + return summedArea; + } + + /** Get the barycenter. + * @return barycenter + */ + public S2Point getBarycenter() { + if (summedBarycenter.getNormSq() == 0) { + return S2Point.NaN; + } else { + return new S2Point(summedBarycenter); + } + } + + } + } 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 c5af3973e..639549254 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 @@ -205,9 +205,6 @@ public class SphericalPolygonsSetTest { Assert.assertEquals(3, count); Assert.assertEquals(1.5 * FastMath.PI, threeOctants.getSize(), 1.0e-10); - Assert.assertEquals(0.0, - new Vector3D(-1, -1, 1).normalize().distance(((S2Point) threeOctants.getBarycenter()).getVector()), - 1.0e-10); } @@ -267,6 +264,116 @@ public class SphericalPolygonsSetTest { } + @Test + public void testSeveralParts() { + double tol = 0.01; + double sinTol = FastMath.sin(tol); + List> boundary = new ArrayList>(); + + // first part: +X, +Y, +Z octant + boundary.add(create(Vector3D.PLUS_J, Vector3D.PLUS_K, Vector3D.PLUS_I, tol, 0.0, 0.5 * FastMath.PI)); + boundary.add(create(Vector3D.PLUS_K, Vector3D.PLUS_I, Vector3D.PLUS_J, tol, 0.0, 0.5 * FastMath.PI)); + boundary.add(create(Vector3D.PLUS_I, Vector3D.PLUS_J, Vector3D.PLUS_K, tol, 0.0, 0.5 * FastMath.PI)); + + // first part: -X, -Y, -Z octant + boundary.add(create(Vector3D.MINUS_J, Vector3D.MINUS_I, Vector3D.MINUS_K, tol, 0.0, 0.5 * FastMath.PI)); + boundary.add(create(Vector3D.MINUS_I, Vector3D.MINUS_K, Vector3D.MINUS_J, tol, 0.0, 0.5 * FastMath.PI)); + boundary.add(create(Vector3D.MINUS_K, Vector3D.MINUS_J, Vector3D.MINUS_I, tol, 0.0, 0.5 * FastMath.PI)); + + SphericalPolygonsSet polygon = new SphericalPolygonsSet(boundary, tol); + + UnitSphereRandomVectorGenerator random = + new UnitSphereRandomVectorGenerator(3, new Well1024a(0xcc5ce49949e0d3ecl)); + 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, polygon.checkPoint(new S2Point(v))); + } else if ((v.getX() < sinTol) && (v.getY() < sinTol) && (v.getZ() < sinTol)) { + Assert.assertEquals(Location.BOUNDARY, polygon.checkPoint(new S2Point(v))); + } else if ((v.getX() > sinTol) && (v.getY() > sinTol) && (v.getZ() > sinTol)) { + Assert.assertEquals(Location.INSIDE, polygon.checkPoint(new S2Point(v))); + } else if ((v.getX() > -sinTol) && (v.getY() > -sinTol) && (v.getZ() > -sinTol)) { + Assert.assertEquals(Location.BOUNDARY, polygon.checkPoint(new S2Point(v))); + } else { + Assert.assertEquals(Location.OUTSIDE, polygon.checkPoint(new S2Point(v))); + } + } + + Assert.assertEquals(FastMath.PI, polygon.getSize(), 1.0e-10); + Assert.assertEquals(3 * FastMath.PI, polygon.getBoundarySize(), 1.0e-10); + + // there should be two separate boundary loops + Assert.assertEquals(2, polygon.getBoundaryLoops().size()); + + } + + @Test + public void testPartWithHole() { + double tol = 0.01; + double alpha = 0.7; + S2Point center = new S2Point(new Vector3D(1, 1, 1)); + SphericalPolygonsSet octant = new SphericalPolygonsSet(center.getVector(), Vector3D.PLUS_K, alpha, 6, tol); + SphericalPolygonsSet hole = new SphericalPolygonsSet(tol, + new S2Point(FastMath.PI / 6, FastMath.PI / 3), + new S2Point(FastMath.PI / 3, FastMath.PI / 3), + new S2Point(FastMath.PI / 4, FastMath.PI / 6)); + SphericalPolygonsSet octantWithHole = + (SphericalPolygonsSet) new RegionFactory().difference(octant, hole); + + for (double phi = center.getPhi() - alpha + 0.1; phi < center.getPhi() + alpha - 0.1; phi += 0.07) { + Location l = octantWithHole.checkPoint(new S2Point(FastMath.PI / 4, phi)); + if (phi < FastMath.PI / 6 || phi > FastMath.PI / 3) { + Assert.assertEquals(Location.INSIDE, l); + } else { + Assert.assertEquals(Location.OUTSIDE, l); + } + } + + // there should be two separate boundary loops + Assert.assertEquals(2, octantWithHole.getBoundaryLoops().size()); + + Assert.assertEquals(octant.getBoundarySize() + hole.getBoundarySize(), octantWithHole.getBoundarySize(), 1.0e-10); + Assert.assertEquals(octant.getSize() - hole.getSize(), octantWithHole.getSize(), 1.0e-10); + + } + + @Test + public void testConcentricSubParts() { + double tol = 0.001; + Vector3D center = new Vector3D(1, 1, 1); + SphericalPolygonsSet hexaOut = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.9, 6, tol); + SphericalPolygonsSet hexaIn = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.8, 6, tol); + SphericalPolygonsSet pentaOut = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.7, 5, tol); + SphericalPolygonsSet pentaIn = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.6, 5, tol); + SphericalPolygonsSet quadriOut = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.5, 4, tol); + SphericalPolygonsSet quadriIn = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.4, 4, tol); + SphericalPolygonsSet triOut = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.25, 3, tol); + SphericalPolygonsSet triIn = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.15, 3, tol); + + RegionFactory factory = new RegionFactory(); + SphericalPolygonsSet hexa = (SphericalPolygonsSet) factory.difference(hexaOut, hexaIn); + SphericalPolygonsSet penta = (SphericalPolygonsSet) factory.difference(pentaOut, pentaIn); + SphericalPolygonsSet quadri = (SphericalPolygonsSet) factory.difference(quadriOut, quadriIn); + SphericalPolygonsSet tri = (SphericalPolygonsSet) factory.difference(triOut, triIn); + SphericalPolygonsSet concentric = + (SphericalPolygonsSet) factory.union(factory.union(hexa, penta), factory.union(quadri, tri)); + + // there should be two separate boundary loops + Assert.assertEquals(8, concentric.getBoundaryLoops().size()); + + Assert.assertEquals(hexaOut.getBoundarySize() + hexaIn.getBoundarySize() + + pentaOut.getBoundarySize() + pentaIn.getBoundarySize() + + quadriOut.getBoundarySize() + quadriIn.getBoundarySize() + + triOut.getBoundarySize() + triIn.getBoundarySize(), + concentric.getBoundarySize(), 1.0e-10); + Assert.assertEquals(hexaOut.getSize() - hexaIn.getSize() + + pentaOut.getSize() - pentaIn.getSize() + + quadriOut.getSize() - quadriIn.getSize() + + triOut.getSize() - triIn.getSize(), + concentric.getSize(), 1.0e-10); + + } + private SubCircle create(Vector3D pole, Vector3D x, Vector3D y, double tolerance, double ... limits) { RegionFactory factory = new RegionFactory(); @@ -280,5 +387,93 @@ public class SphericalPolygonsSetTest { return new SubCircle(phased, set); } +// private static void displayCells(String name, SphericalPolygonsSet polygons, double maxStep, double grid) { +// TesselationDisplay td = new TesselationDisplay(name, maxStep, grid, polygons.getTolerance()); +// polygons.getTree(false).visit(td); +// td.close(); +// } +// +// private static class TesselationDisplay implements BSPTreeVisitor { +// private double maxStep; +// private double grid; +// private double tolerance; +// private PrintStream out; +// public TesselationDisplay(String name, double maxStep, double grid, double tolerance) { +// try { +// this.out = new PrintStream(new File(new File(System.getProperty("user.home")), name)); +// this.maxStep = maxStep; +// this.grid = grid; +// this.tolerance = tolerance; +// } catch (IOException ioe) { +// Assert.fail(ioe.getMessage()); +// } +// } +// public Order visitOrder(BSPTree node) { +// return Order.MINUS_PLUS_SUB; +// } +// public void visitInternalNode(BSPTree node) { +// } +// public void visitLeafNode(BSPTree node) { +// if ((Boolean) node.getAttribute()) { +// final SphericalPolygonsSet cell = +// new SphericalPolygonsSet(node.pruneAroundConvexCell(true, false, null), tolerance); +// display(out, cell, maxStep, grid); +// } +// } +// public void close() { +// out.close(); +// } +// } +// +// private static void display(String name, SphericalPolygonsSet polygons, double maxStep, double grid) { +// PrintStream out = null; +// try { +// out = new PrintStream(new File(new File(System.getProperty("user.home")), name)); +// display(out, polygons, maxStep, grid); +// } catch (IOException ioe) { +// Assert.fail(ioe.getMessage()); +// } finally { +// if (out != null) { +// out.close(); +// } +// } +// } +// +// private static void display(PrintStream out, SphericalPolygonsSet polygons, double maxStep, double grid) { +// for (final SphericalPolygonsSet.Vertex start : polygons.getBoundaryLoops()) { +// int n = 0; +// for (SphericalPolygonsSet.Edge edge = start.getOutgoing(); +// n == 0 || edge.getStart() != start; +// edge = edge.getEnd().getOutgoing()) { +// int m = (int) FastMath.ceil(edge.getLength() / maxStep); +// double step = edge.getLength() / m; +// for (int i = 0; i <= m; ++i) { +// final Vector3D p = edge.getPointAt(i * step); +// out.format(Locale.US, "%9.6f %9.6f %9.6f%n", p.getX(), p.getY(), p.getZ()); +// } +// n++; +// } +// out.println("&"); +// } +// if (grid > 0) { +// for (double phi = grid; phi < FastMath.PI - grid; phi += grid) { +// Location previous = Location.OUTSIDE; +// for (double theta = 0; theta < MathUtils.TWO_PI; theta += grid) { +// S2Point point = new S2Point(theta, phi); +// Location current = polygons.checkPoint(point); +// if (current == Location.OUTSIDE) { +// if (previous != Location.OUTSIDE) { +// out.println("&"); +// } +// } else { +// out.format(Locale.US, "%9.6f %9.6f %9.6f%n", +// point.getVector().getX(), point.getVector().getY(), point.getVector().getZ()); +// } +// previous = current; +// } +// out.println("&"); +// } +// } +// } }