Skip to content

Commit

Permalink
Merge 31e4a6d into 6d9bc1a
Browse files Browse the repository at this point in the history
  • Loading branch information
@darkma773r
darkma773r committed Jan 24, 2018
2 parents 6d9bc1a + 31e4a6d commit 3768292
Show file tree
Hide file tree
Showing 2 changed files with 181 additions and 30 deletions.
94 changes: 64 additions & 30 deletions src/main/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSet.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -353,30 +353,65 @@ public PolyhedronsSet buildNew(final BSPTree<Euclidean3D> tree) {
/** {@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)*&Sigma;<sub>F</sub>[(C<sub>F</sub>&sdot;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} */
Expand Down Expand Up @@ -404,34 +439,33 @@ public void visitInternalNode(final BSPTree<Euclidean3D> node) {
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.
Expand Down
117 changes: 117 additions & 0 deletions src/test/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSetTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,76 @@

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);
Expand Down Expand Up @@ -97,6 +167,53 @@ public void testBox() {
});
}

@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);
Expand Down

0 comments on commit 3768292

Please sign in to comment.