Repository: commons-math Updated Branches: refs/heads/master e4a613ffe -> 3ea45970d
MATH-1442: fixing size calculation in PolyhedronsSet and updating the code to make the algorithm clearer Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/31e4a6df Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/31e4a6df Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/31e4a6df Branch: refs/heads/master Commit: 31e4a6df19cf4a13f8f0e9048198ce7488ad2614 Parents: 6d9bc1a Author: darkma773r <[email protected]> Authored: Tue Jan 23 22:11:30 2018 -0500 Committer: darkma773r <[email protected]> Committed: Tue Jan 23 22:11:30 2018 -0500 ---------------------------------------------------------------------- .../euclidean/threed/PolyhedronsSet.java | 94 ++++++++++----- .../euclidean/threed/PolyhedronsSetTest.java | 117 +++++++++++++++++++ 2 files changed, 181 insertions(+), 30 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/31e4a6df/src/main/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSet.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSet.java b/src/main/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSet.java index cf15f22..4459d34 100644 --- a/src/main/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSet.java +++ b/src/main/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSet.java @@ -353,30 +353,65 @@ public class PolyhedronsSet extends AbstractRegion<Euclidean3D, Euclidean2D> { /** {@inheritDoc} */ @Override protected void computeGeometricalProperties() { - - // compute the contribution of all boundary facets - getTree(true).visit(new FacetsContributionVisitor()); - - if (getSize() < 0) { - // the polyhedrons set as a finite outside - // surrounded by an infinite inside + // check simple cases first + if (isEmpty()) { + setSize(0.0); + setBarycenter((Point<Euclidean3D>) Cartesian3D.NaN); + } + else if (isFull()) { setSize(Double.POSITIVE_INFINITY); setBarycenter((Point<Euclidean3D>) Cartesian3D.NaN); - } else { - // the polyhedrons set is finite, apply the remaining scaling factors - setSize(getSize() / 3.0); - setBarycenter((Point<Euclidean3D>) new Cartesian3D(1.0 / (4 * getSize()), (Cartesian3D) getBarycenter())); } + else { + // not empty or full; compute the contribution of all boundary facets + final FacetsContributionVisitor contributionVisitor = new FacetsContributionVisitor(); + getTree(true).visit(contributionVisitor); + final double size = contributionVisitor.getSize(); + final Cartesian3D barycenter = contributionVisitor.getBarycenter(); + + if (size < 0) { + // the polyhedrons set is a finite outside surrounded by an infinite inside + setSize(Double.POSITIVE_INFINITY); + setBarycenter((Point<Euclidean3D>) Cartesian3D.NaN); + } else { + // the polyhedrons set is finite + setSize(size); + setBarycenter((Point<Euclidean3D>) barycenter); + } + } } - /** Visitor computing geometrical properties. */ - private class FacetsContributionVisitor implements BSPTreeVisitor<Euclidean3D> { + /** Visitor computing polyhedron geometrical properties. + * The volume of the polyhedron is computed using the equation + * <code>V = (1/3)*Σ<sub>F</sub>[(C<sub>F</sub>⋅N<sub>F</sub>)*area(F)]</code>, + * where <code>F</code> represents each face in the polyhedron, <code>C<sub>F</sub></code> + * represents the barycenter of the face, and <code>N<sub>F</sub></code> represents the + * normal of the face. (More details can be found in the article + * <a href="https://en.wikipedia.org/wiki/Polyhedron#Volume";>here</a>.) + * This essentially splits up the polyhedron into pyramids with a polyhedron + * face forming the base of each pyramid. + * The barycenter is computed in a similar way. The barycenter of each pyramid + * is calculated using the fact that it is located 3/4 of the way along the + * line from the apex to the base. The polyhedron barycenter then becomes + * the volume-weighted average of these pyramid centers. + */ + private static class FacetsContributionVisitor implements BSPTreeVisitor<Euclidean3D> { + + private double volumeSum; + private Cartesian3D barycenterSum = Cartesian3D.ZERO; + + public double getSize() { + // apply the 1/3 pyramid volume scaling factor + return volumeSum / 3.0; + } - /** Simple constructor. */ - FacetsContributionVisitor() { - setSize(0); - setBarycenter((Point<Euclidean3D>) new Cartesian3D(0, 0, 0)); + public Cartesian3D getBarycenter() { + // Since the volume we used when adding together the facet contributions + // was 3x the actual pyramid size, we'll multiply by 1/4 here instead + // of 3/4. This will make the overall polyhedron volume the weighted + // average of the pyramid barycenters. + return new Cartesian3D(1.0 / (4 * getSize()), barycenterSum); } /** {@inheritDoc} */ @@ -404,34 +439,33 @@ public class PolyhedronsSet extends AbstractRegion<Euclidean3D, Euclidean2D> { public void visitLeafNode(final BSPTree<Euclidean3D> node) { } - /** Add he contribution of a boundary facet. + /** Add the contribution of a boundary facet. * @param facet boundary facet * @param reversed if true, the facet has the inside on its plus side */ private void addContribution(final SubHyperplane<Euclidean3D> facet, final boolean reversed) { final Region<Euclidean2D> polygon = ((SubPlane) facet).getRemainingRegion(); - final double area = polygon.getSize(); + final double area = polygon.getSize(); if (Double.isInfinite(area)) { - setSize(Double.POSITIVE_INFINITY); - setBarycenter((Point<Euclidean3D>) Cartesian3D.NaN); + volumeSum = Double.POSITIVE_INFINITY; + barycenterSum = Cartesian3D.NaN; } else { + final Plane plane = (Plane) facet.getHyperplane(); + final Cartesian3D facetBarycenter = plane.toSpace(polygon.getBarycenter()); - final Plane plane = (Plane) facet.getHyperplane(); - final Cartesian3D facetB = plane.toSpace(polygon.getBarycenter()); - double scaled = area * facetB.dotProduct(plane.getNormal()); + // the volume here is actually 3x the actual pyramid volume; we'll apply + // the final scaling all at once at the end + double scaledVolume = area * facetBarycenter.dotProduct(plane.getNormal()); if (reversed) { - scaled = -scaled; + scaledVolume = -scaledVolume; } - setSize(getSize() + scaled); - setBarycenter((Point<Euclidean3D>) new Cartesian3D(1.0, (Cartesian3D) getBarycenter(), scaled, facetB)); - + volumeSum += scaledVolume; + barycenterSum = new Cartesian3D(1.0, barycenterSum, scaledVolume, facetBarycenter); } - } - } /** Get the first sub-hyperplane crossed by a semi-infinite line. http://git-wip-us.apache.org/repos/asf/commons-math/blob/31e4a6df/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSetTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSetTest.java b/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSetTest.java index 43553ea..1a64d2a 100644 --- a/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSetTest.java +++ b/src/test/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSetTest.java @@ -59,6 +59,76 @@ import org.junit.Test; public class PolyhedronsSetTest { @Test + public void testFull() { + // act + PolyhedronsSet tree = new PolyhedronsSet(1e-10); + + // assert + Assert.assertEquals(Double.POSITIVE_INFINITY, tree.getSize(), 1.0e-10); + Assert.assertEquals(0.0, tree.getBoundarySize(), 1.0e-10); + Cartesian3D center = (Cartesian3D) tree.getBarycenter(); + Assert.assertEquals(Double.NaN, center.getX(), 1e-10); + Assert.assertEquals(Double.NaN, center.getY(), 1e-10); + Assert.assertEquals(Double.NaN, center.getZ(), 1e-10); + Assert.assertEquals(true, tree.isFull()); + Assert.assertEquals(false, tree.isEmpty()); + + checkPoints(Region.Location.INSIDE, tree, new Cartesian3D[] { + new Cartesian3D(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY), + Cartesian3D.ZERO, + new Cartesian3D(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY) + }); + } + + @Test + public void testEmpty() { + // act + PolyhedronsSet tree = new PolyhedronsSet(new BSPTree<Euclidean3D>(Boolean.FALSE), 1e-10); + + // assert + Assert.assertEquals(0.0, tree.getSize(), 1.0e-10); + Assert.assertEquals(0.0, tree.getBoundarySize(), 1.0e-10); + Cartesian3D center = (Cartesian3D) tree.getBarycenter(); + Assert.assertEquals(Double.NaN, center.getX(), 1e-10); + Assert.assertEquals(Double.NaN, center.getY(), 1e-10); + Assert.assertEquals(Double.NaN, center.getZ(), 1e-10); + Assert.assertEquals(false, tree.isFull()); + Assert.assertEquals(true, tree.isEmpty()); + + checkPoints(Region.Location.OUTSIDE, tree, new Cartesian3D[] { + new Cartesian3D(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY), + Cartesian3D.ZERO, + new Cartesian3D(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY) + }); + } + + @Test + public void testSingleInfiniteBoundary() { + // arrange + Plane plane = new Plane(Cartesian3D.ZERO, Cartesian3D.PLUS_K, 1e-10); + + List<SubHyperplane<Euclidean3D>> boundaries = new ArrayList<>(); + boundaries.add(plane.wholeHyperplane()); + + // act + PolyhedronsSet tree = new PolyhedronsSet(boundaries, 1e-10); + + // assert + Assert.assertEquals(Double.POSITIVE_INFINITY, tree.getSize(), 1.0e-10); + Assert.assertEquals(Double.POSITIVE_INFINITY, tree.getBoundarySize(), 1.0e-10); + Cartesian3D center = (Cartesian3D) tree.getBarycenter(); + Assert.assertEquals(Double.NaN, center.getX(), 1e-10); + Assert.assertEquals(Double.NaN, center.getY(), 1e-10); + Assert.assertEquals(Double.NaN, center.getZ(), 1e-10); + Assert.assertEquals(false, tree.isFull()); + Assert.assertEquals(false, tree.isEmpty()); + + checkPoints(Region.Location.BOUNDARY, tree, new Cartesian3D[] { + Cartesian3D.ZERO + }); + } + + @Test public void testBox() { PolyhedronsSet tree = new PolyhedronsSet(0, 1, 0, 1, 0, 1, 1.0e-10); Assert.assertEquals(1.0, tree.getSize(), 1.0e-10); @@ -98,6 +168,53 @@ public class PolyhedronsSetTest { } @Test + public void testInvertedBox() { + // arrange + PolyhedronsSet tree = new PolyhedronsSet(0, 1, 0, 1, 0, 1, 1.0e-10); + + // act + tree = (PolyhedronsSet) new RegionFactory<Euclidean3D>().getComplement(tree); + + // assert + Assert.assertEquals(Double.POSITIVE_INFINITY, tree.getSize(), 1.0e-10); + Assert.assertEquals(6.0, tree.getBoundarySize(), 1.0e-10); + + Cartesian3D barycenter = (Cartesian3D) tree.getBarycenter(); + Assert.assertEquals(Double.NaN, barycenter.getX(), 1.0e-10); + Assert.assertEquals(Double.NaN, barycenter.getY(), 1.0e-10); + Assert.assertEquals(Double.NaN, barycenter.getZ(), 1.0e-10); + + for (double x = -0.25; x < 1.25; x += 0.1) { + boolean xOK = (x < 0.0) || (x > 1.0); + for (double y = -0.25; y < 1.25; y += 0.1) { + boolean yOK = (y < 0.0) || (y > 1.0); + for (double z = -0.25; z < 1.25; z += 0.1) { + boolean zOK = (z < 0.0) || (z > 1.0); + Region.Location expected = + (xOK || yOK || zOK) ? Region.Location.INSIDE : Region.Location.OUTSIDE; + Assert.assertEquals(expected, tree.checkPoint(new Cartesian3D(x, y, z))); + } + } + } + checkPoints(Region.Location.BOUNDARY, tree, new Cartesian3D[] { + new Cartesian3D(0.0, 0.5, 0.5), + new Cartesian3D(1.0, 0.5, 0.5), + new Cartesian3D(0.5, 0.0, 0.5), + new Cartesian3D(0.5, 1.0, 0.5), + new Cartesian3D(0.5, 0.5, 0.0), + new Cartesian3D(0.5, 0.5, 1.0) + }); + checkPoints(Region.Location.INSIDE, tree, new Cartesian3D[] { + new Cartesian3D(0.0, 1.2, 1.2), + new Cartesian3D(1.0, 1.2, 1.2), + new Cartesian3D(1.2, 0.0, 1.2), + new Cartesian3D(1.2, 1.0, 1.2), + new Cartesian3D(1.2, 1.2, 0.0), + new Cartesian3D(1.2, 1.2, 1.0) + }); + } + + @Test public void testTetrahedron() throws MathArithmeticException { Cartesian3D vertex1 = new Cartesian3D(1, 2, 3); Cartesian3D vertex2 = new Cartesian3D(2, 2, 4);
