Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'fix_MATH-1442' 

Closes #75
This commit is contained in:
Gilles 2018-01-26 01:00:45 +01:00
parent e4a613ffe0 32c5735c8a
commit 3ea45970db
2 changed files with 193 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

106
src/main/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSet.java
View File

@ -353,30 +353,77 @@ 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)*&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> {
/** Simple constructor. */
FacetsContributionVisitor() {
setSize(0);
setBarycenter((Point<Euclidean3D>) new Cartesian3D(0, 0, 0));
/** 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 the volume.
*/
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 the barycenter.
*/
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 to adjust for the actual barycenter position in each pyramid.
return new Cartesian3D(1.0 / (4 * getSize()), barycenterSum);
}
/** {@inheritDoc} */
@ -404,34 +451,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.

117
src/test/java/org/apache/commons/math4/geometry/euclidean/threed/PolyhedronsSetTest.java
View File

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