From 31e4a6df19cf4a13f8f0e9048198ce7488ad2614 Mon Sep 17 00:00:00 2001 From: darkma773r Date: Tue, 23 Jan 2018 22:11:30 -0500 Subject: [PATCH 1/3] MATH-1442: fixing size calculation in PolyhedronsSet and updating the code to make the algorithm clearer --- .../euclidean/threed/PolyhedronsSet.java | 94 +++++++++----- .../euclidean/threed/PolyhedronsSetTest.java | 117 ++++++++++++++++++ 2 files changed, 181 insertions(+), 30 deletions(-) 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 cf15f2259..4459d3408 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 { /** {@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) Cartesian3D.NaN); + } + else if (isFull()) { setSize(Double.POSITIVE_INFINITY); setBarycenter((Point) Cartesian3D.NaN); - } else { - // the polyhedrons set is finite, apply the remaining scaling factors - setSize(getSize() / 3.0); - setBarycenter((Point) 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) Cartesian3D.NaN); + } else { + // the polyhedrons set is finite + setSize(size); + setBarycenter((Point) barycenter); + } + } } - /** Visitor computing geometrical properties. */ - private class FacetsContributionVisitor implements BSPTreeVisitor { + /** Visitor computing polyhedron geometrical properties. + * The volume of the polyhedron is computed using the equation + * V = (1/3)*ΣF[(CF⋅NF)*area(F)], + * where F represents each face in the polyhedron, CF + * represents the barycenter of the face, and NF represents the + * normal of the face. (More details can be found in the article + * here.) + * 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 { - /** Simple constructor. */ - FacetsContributionVisitor() { - setSize(0); - setBarycenter((Point) new Cartesian3D(0, 0, 0)); + private double volumeSum; + private Cartesian3D barycenterSum = Cartesian3D.ZERO; + + public double getSize() { + // apply the 1/3 pyramid volume scaling factor + return volumeSum / 3.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 { public void visitLeafNode(final BSPTree 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 facet, final boolean reversed) { final Region polygon = ((SubPlane) facet).getRemainingRegion(); - final double area = polygon.getSize(); + final double area = polygon.getSize(); if (Double.isInfinite(area)) { - setSize(Double.POSITIVE_INFINITY); - setBarycenter((Point) 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) 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. 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 43553ea01..1a64d2a2d 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 @@ -58,6 +58,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(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> 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); @@ -97,6 +167,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().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); From dc0465e864d7eace765434834d19575a106dbf64 Mon Sep 17 00:00:00 2001 From: darkma773r Date: Wed, 24 Jan 2018 22:08:18 -0500 Subject: [PATCH 2/3] MATH-1442: adding javadocs to PolyhedronsSet.FacetContributionVisitor --- .../euclidean/threed/PolyhedronsSet.java | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) 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 4459d3408..d7f846a66 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 @@ -398,19 +398,31 @@ public class PolyhedronsSet extends AbstractRegion { */ private static class FacetsContributionVisitor implements BSPTreeVisitor { + /** Accumulator for facet volume contributions. */ private double volumeSum; + + /** Accumulator for barycenter contributions. */ private Cartesian3D barycenterSum = Cartesian3D.ZERO; + /** Returns the total computed size (ie, volume) of the polyhedron. + * This value will be negative if the polyhedron is "inside-out", meaning + * that it has a finite outside surrounded by an infinite inside. + * @return + */ public double getSize() { // apply the 1/3 pyramid volume scaling factor return volumeSum / 3.0; } + /** Returns the computed barycenter. This is the volume-weighted average + * of contributions from all facets. All coordinates will be NaN if the + * region is infinite. + * @return + */ 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. + // of 3/4 to adjust for the actual barycenter position in each pyramid. return new Cartesian3D(1.0 / (4 * getSize()), barycenterSum); } From 32c5735c8afbf7e9eb4df30cfcac0d8251dda28e Mon Sep 17 00:00:00 2001 From: Gilles Date: Fri, 26 Jan 2018 01:00:13 +0100 Subject: [PATCH 3/3] Javadoc. --- .../math4/geometry/euclidean/threed/PolyhedronsSet.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 d7f846a66..ab2728359 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 @@ -407,7 +407,7 @@ public class PolyhedronsSet extends AbstractRegion { /** Returns the total computed size (ie, volume) of the polyhedron. * This value will be negative if the polyhedron is "inside-out", meaning * that it has a finite outside surrounded by an infinite inside. - * @return + * @return the volume. */ public double getSize() { // apply the 1/3 pyramid volume scaling factor @@ -417,7 +417,7 @@ public class PolyhedronsSet extends AbstractRegion { /** Returns the computed barycenter. This is the volume-weighted average * of contributions from all facets. All coordinates will be NaN if the * region is infinite. - * @return + * @return the barycenter. */ public Cartesian3D getBarycenter() { // Since the volume we used when adding together the facet contributions