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
This commit is contained in:
Luc Maisonobe 2014-01-14 18:26:32 +00:00
parent 692142457a
commit 9b3a46e205
4 changed files with 422 additions and 60 deletions

View File

@ -665,6 +665,46 @@ public class BSPTree<S extends Space> {
}
/** Prune a tree around a cell.
* <p>
* 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.
* </p>
* @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<S> pruneAroundConvexCell(final Object cellAttribute,
final Object otherLeafsAttributes,
final Object internalAttributes) {
// build the current cell leaf
BSPTree<S> tree = new BSPTree<S>(cellAttribute);
// build the pruned tree bottom-up
for (BSPTree<S> current = this; current.parent != null; current = current.parent) {
final SubHyperplane<S> parentCut = current.parent.cut.copySelf();
final BSPTree<S> sibling = new BSPTree<S>(otherLeafsAttributes);
if (current == current.parent.plus) {
tree = new BSPTree<S>(parentCut, tree, sibling, internalAttributes);
} else {
tree = new BSPTree<S>(parentCut, sibling, tree, internalAttributes);
}
}
return tree;
}
/** Chop off parts of the tree.
* <p>The instance is modified in place, all the parts that are on
* the minus side of the chopping hyperplane are discarded, only the

View File

@ -520,8 +520,8 @@ public class ArcsSet extends AbstractRegion<Sphere1D, Sphere1D> 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

View File

@ -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<Sphere2D, Sphere1D> {
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.
* <p>The leaf nodes of the BSP tree <em>must</em> have a
* {@code Boolean} attribute representing the inside status of
@ -140,6 +155,27 @@ public class SphericalPolygonsSet extends AbstractRegion<Sphere2D, Sphere1D> {
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.
* <p>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<Sphere2D, Sphere1D> {
@Override
protected void computeGeometricalProperties() throws MathIllegalStateException {
final List<Vertex> boundary = getBoundaryLoops();
final BSPTree<Sphere2D> tree = getTree(true);
if (tree.getCut() == null) {
// the instance has a single cell without any boundaries
if (boundary.isEmpty()) {
final BSPTree<Sphere2D> 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<Sphere2D, Sphere1D> {
}
/** Visitor computing geometrical properties. */
private static class PropertiesComputer implements BSPTreeVisitor<Sphere2D> {
/** 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<Sphere2D> node) {
return Order.MINUS_SUB_PLUS;
}
/** {@inheritDoc} */
public void visitInternalNode(final BSPTree<Sphere2D> node) {
// nothing to do here
}
/** {@inheritDoc} */
public void visitLeafNode(final BSPTree<Sphere2D> 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<Vertex> 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);
}
}
}
}

View File

@ -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<SubHyperplane<Sphere2D>> boundary = new ArrayList<SubHyperplane<Sphere2D>>();
// 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<Sphere2D>().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<Sphere2D> factory = new RegionFactory<Sphere2D>();
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<Sphere1D> factory = new RegionFactory<Sphere1D>();
@ -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<Sphere2D> {
// 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<Sphere2D> node) {
// return Order.MINUS_PLUS_SUB;
// }
// public void visitInternalNode(BSPTree<Sphere2D> node) {
// }
// public void visitLeafNode(BSPTree<Sphere2D> 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("&");
// }
// }
// }
}