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

Finalized fix for MATH-880. 

This version should be much more robust on polygons computation than the
earlier ones. Note that the hyperplaneThickness parameter is a key
tuning parameter in difficult cases like the one involved in the issue.
As shown in the corresponding unit test, the value had to be raised to
1.0e-8 in order for the test to pass correctly.

JIRA: MATH-880

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1401022 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2012-10-22 19:21:36 +00:00
parent c2e855d6a0
commit 2a9cbbab46
2 changed files with 275 additions and 144 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

361
src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
View File

@ -111,6 +111,20 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
* constructor} using {@link SubHyperplane subhyperplanes}.</p> * constructor} using {@link SubHyperplane subhyperplanes}.</p>
* <p>If the list is empty, the region will represent the whole * <p>If the list is empty, the region will represent the whole
* space.</p> * space.</p>
* <p>
* Polygons with thin pikes or dents are inherently difficult to handle because
* they involve lines with almost opposite directions at some vertices. Polygons
* whose vertices come from some physical measurement with noise are also
* difficult because an edge that should be straight may be broken in lots of
* different pieces with almost equal directions. In both cases, computing the
* lines intersections is not numerically robust due to the almost 0 or almost
* &pi; angle. Such cases need to carefully adjust the {@code hyperplaneThickness}
* parameter. A too small value would often lead to completely wrong polygons
* with large area wrongly identified as inside or outside. Large values are
* often much safer. As a rule of thumb, a value slightly below the size of the
* most accurate detail needed is a good value for the {@code hyperplaneThickness}
* parameter.
* </p>
* @param hyperplaneThickness tolerance below which points are considered to * @param hyperplaneThickness tolerance below which points are considered to
* belong to the hyperplane (which is therefore more a slab) * belong to the hyperplane (which is therefore more a slab)
* @param vertices vertices of the simple loop boundary * @param vertices vertices of the simple loop boundary
@ -157,20 +171,50 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
private static BSPTree<Euclidean2D> verticesToTree(final double hyperplaneThickness, private static BSPTree<Euclidean2D> verticesToTree(final double hyperplaneThickness,
final Vector2D ... vertices) { final Vector2D ... vertices) {
if (vertices.length == 0) { final int n = vertices.length;
if (n == 0) {
// the tree represents the whole space // the tree represents the whole space
return new BSPTree<Euclidean2D>(Boolean.TRUE); return new BSPTree<Euclidean2D>(Boolean.TRUE);
} }
// at start, none of the edges have been processed // build the vertices
final BSPTree<Euclidean2D> tree = new BSPTree<Euclidean2D>(); final Vertex[] vArray = new Vertex[n];
List<Vertex> list = new ArrayList<PolygonsSet.Vertex>(vertices.length); for (int i = 0; i < n; ++i) {
for (final Vector2D vertex : vertices) { vArray[i] = new Vertex(vertices[i]);
list.add(new Vertex(vertex)); }
// build the edges
List<Edge> edges = new ArrayList<Edge>();
for (int i = 0; i < n; ++i) {
// get the endpoints of the edge
final Vertex start = vArray[i];
final Vertex end = vArray[(i + 1) % n];
// get the line supporting the edge, taking care not to recreate it
// if it was already created earlier due to another edge being aligned
// with the current one
Line line = start.sharedLineWith(end);
if (line == null) {
line = new Line(start.getLocation(), end.getLocation());
}
// create the edge and store it
edges.add(new Edge(start, end, line));
// check if another vertex also happens to be on this line
for (final Vertex vertex : vArray) {
if (vertex != start && vertex != end &&
FastMath.abs(line.getOffset(vertex.getLocation())) <= hyperplaneThickness) {
vertex.bindWith(line);
}
}
} }
// build the tree top-down // build the tree top-down
insertVertices(hyperplaneThickness, tree, list); final BSPTree<Euclidean2D> tree = new BSPTree<Euclidean2D>();
insertEdges(hyperplaneThickness, tree, edges);
return tree; return tree;
@ -181,45 +225,32 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
* belong to the hyperplane (which is therefore more a slab) * belong to the hyperplane (which is therefore more a slab)
* @param node current tree node (it is a leaf node at the beginning * @param node current tree node (it is a leaf node at the beginning
* of the call) * of the call)
* @param vertices list of vertices belonging to the boundary of the * @param edges list of edges to insert in the cell defined by this node
* cell defined by the node * (excluding edges not belonging to the cell defined by this node)
*/ */
private static void insertVertices(final double hyperplaneThickness, private static void insertEdges(final double hyperplaneThickness,
final BSPTree<Euclidean2D> node, final BSPTree<Euclidean2D> node,
final List<Vertex> vertices) { final List<Edge> edges) {
Vertex current = vertices.get(vertices.size() - 1); // find an edge with an hyperplane that can be inserted in the node
int index = 0; int index = 0;
Line inserted = null; Edge inserted =null;
while (inserted == null && index < vertices.size()) { while (inserted == null && index < edges.size()) {
final Vertex previous = current; inserted = edges.get(index++);
current = vertices.get(index++); if (inserted.getNode() == null) {
if (previous.outgoingNeedsProcessing() && current.incomingNeedsProcessing()) { if (node.insertCut(inserted.getLine())) {
inserted.setNode(node);
if (previous.shareNodeWith(current)) {
// both vertices are already handled by an existing node,
// closer to the tree root, they were probably created
// when split points were introduced
inserted = null;
} else {
inserted = new Line(previous.getLocation(), current.getLocation());
if (node.insertCut(inserted)) {
previous.addNode(node);
previous.outgoingProcessed();
current.addNode(node);
current.incomingProcessed();
} else { } else {
inserted = null; inserted = null;
} }
} else {
} inserted = null;
} }
} }
if (node.getCut() == null) { if (inserted == null) {
// no suitable edge was found, the node remains a leaf node
// we need to set its inside/outside boolean indicator
final BSPTree<Euclidean2D> parent = node.getParent(); final BSPTree<Euclidean2D> parent = node.getParent();
if (parent == null || node == parent.getMinus()) { if (parent == null || node == parent.getMinus()) {
node.setAttribute(Boolean.TRUE); node.setAttribute(Boolean.TRUE);
@ -229,67 +260,58 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
return; return;
} }
// distribute the remaining vertices in the two sub-trees // we have split the node by inserted an edge as a cut sub-hyperplane
Side currentSide = Side.HYPER; // distribute the remaining edges in the two sub-trees
final List<Vertex> plusList = new ArrayList<Vertex>(); final List<Edge> plusList = new ArrayList<Edge>();
plusList.add(current); final List<Edge> minusList = new ArrayList<Edge>();
int plusCount = 0; for (final Edge edge : edges) {
final List<Vertex> minusList = new ArrayList<Vertex>(); if (edge != inserted) {
minusList.add(current); final double startOffset = inserted.getLine().getOffset(edge.getStart().getLocation());
int minusCount = 0; final double endOffset = inserted.getLine().getOffset(edge.getEnd().getLocation());
while (index < vertices.size()) { Side startSide = (FastMath.abs(startOffset) <= hyperplaneThickness) ?
final Vertex previous = current; Side.HYPER : ((startOffset < 0) ? Side.MINUS : Side.PLUS);
final Side previousSide = currentSide; Side endSide = (FastMath.abs(endOffset) <= hyperplaneThickness) ?
current = vertices.get(index++); Side.HYPER : ((endOffset < 0) ? Side.MINUS : Side.PLUS);
final double currentOffset = inserted.getOffset(current.getLocation()); switch (startSide) {
currentSide = (FastMath.abs(currentOffset) <= hyperplaneThickness) ?
Side.HYPER :
((currentOffset < 0) ? Side.MINUS : Side.PLUS);
switch (currentSide) {
case PLUS: case PLUS:
if (previousSide == Side.MINUS) { if (endSide == Side.MINUS) {
// we need to insert a split point on the hyperplane // we need to insert a split point on the hyperplane
final Line line = new Line(previous.getLocation(), current.getLocation()); final Vertex splitPoint = edge.split(inserted.getLine());
final Vertex splitPoint = new Vertex(inserted.intersection(line)); minusList.add(splitPoint.getOutgoing());
splitPoint.addNode(node); plusList.add(splitPoint.getIncoming());
minusList.add(splitPoint); } else {
plusList.add(splitPoint); plusList.add(edge);
}
plusList.add(current);
if (current.incomingNeedsProcessing() || current.outgoingNeedsProcessing()) {
++plusCount;
} }
break; break;
case MINUS: case MINUS:
if (previousSide == Side.PLUS) { if (endSide == Side.PLUS) {
// we need to insert a split point on the hyperplane // we need to insert a split point on the hyperplane
final Line line = new Line(previous.getLocation(), current.getLocation()); final Vertex splitPoint = edge.split(inserted.getLine());
final Vertex splitPoint = new Vertex(inserted.intersection(line)); minusList.add(splitPoint.getIncoming());
splitPoint.addNode(node); plusList.add(splitPoint.getOutgoing());
minusList.add(splitPoint); } else {
plusList.add(splitPoint); minusList.add(edge);
}
minusList.add(current);
if (current.incomingNeedsProcessing() || current.outgoingNeedsProcessing()) {
++minusCount;
} }
break; break;
default: default:
current.addNode(node); if (endSide == Side.PLUS) {
plusList.add(current); plusList.add(edge);
minusList.add(current); } else if (endSide == Side.MINUS) {
minusList.add(edge);
}
break; break;
} }
} }
}
// recurse through lower levels // recurse through lower levels
if (plusCount > 0) { if (!plusList.isEmpty()) {
insertVertices(hyperplaneThickness, node.getPlus(), plusList); insertEdges(hyperplaneThickness, node.getPlus(), plusList);
} else { } else {
node.getPlus().setAttribute(Boolean.FALSE); node.getPlus().setAttribute(Boolean.FALSE);
} }
if (minusCount > 0) { if (!minusList.isEmpty()) {
insertVertices(hyperplaneThickness, node.getMinus(), minusList); insertEdges(hyperplaneThickness, node.getMinus(), minusList);
} else { } else {
node.getMinus().setAttribute(Boolean.TRUE); node.getMinus().setAttribute(Boolean.TRUE);
} }
@ -302,23 +324,23 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
/** Vertex location. */ /** Vertex location. */
private final Vector2D location; private final Vector2D location;
/** Nodes associated with the hyperplane containing this vertex. */ /** Incoming edge. */
private final List<BSPTree<Euclidean2D>> nodes; private Edge incoming;
/** Indicator for incoming edges that still need processing. */ /** Outgoing edge. */
private boolean incomingNeedsProcessing; private Edge outgoing;
/** Indicator for outgoing edges that still need processing. */ /** Lines bound with this vertex. */
private boolean outgoingNeedsProcessing; private final List<Line> lines;
/** Build a non-processed vertex not owned by any node yet. /** Build a non-processed vertex not owned by any node yet.
* @param location vertex location * @param location vertex location
*/ */
public Vertex(final Vector2D location) { public Vertex(final Vector2D location) {
this.location = location; this.location = location;
this.nodes = new ArrayList<BSPTree<Euclidean2D>>(); this.incoming = null;
this.incomingNeedsProcessing = true; this.outgoing = null;
this.outgoingNeedsProcessing = true; this.lines = new ArrayList<Line>();
} }
/** Get Vertex location. /** Get Vertex location.
@ -328,57 +350,160 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
return location; return location;
} }
/** Check if the instance and another vertex share a node. /** Bind a line considered to contain this vertex.
* @param line line to bind with this vertex
*/
public void bindWith(final Line line) {
lines.add(line);
}
/** Get the common line bound with both the instance and another vertex, if any.
* <p> * <p>
* When two vertices share a node, this means they are already handled * When two vertices are both bound to the same line, this means they are
* by the hyperplane of this node, so there is no need to create a cut * already handled by node associated with this line, so there is no need
* hyperplane for them. * to create a cut hyperplane for them.
* </p> * </p>
* @param vertex other vertex to check instance against * @param vertex other vertex to check instance against
* @return true if the instance and another vertex share a node * @return line bound with both the instance and another vertex, or null if the
* two vertices do not share a line yet
*/ */
public boolean shareNodeWith(final Vertex vertex) { public Line sharedLineWith(final Vertex vertex) {
for (final BSPTree<Euclidean2D> node1 : nodes) { for (final Line line1 : lines) {
for (final BSPTree<Euclidean2D> node2 : vertex.nodes) { for (final Line line2 : vertex.lines) {
if (node1 == node2) { if (line1 == line2) {
return true; return line1;
} }
} }
} }
return false; return null;
} }
/** Add a node whose hyperplane contains this vertex. /** Set incoming edge.
* @param node node whose hyperplane contains this vertex * <p>
* The line supporting the incoming edge is automatically bound
* with the instance.
* </p>
* @param incoming incoming edge
*/ */
public void addNode(final BSPTree<Euclidean2D> node) { public void setIncoming(final Edge incoming) {
nodes.add(node); this.incoming = incoming;
bindWith(incoming.getLine());
} }
/** Check incoming edge processed indicator. /** Get incoming edge.
* @return true if incoming edge needs processing * @return incoming edge
*/ */
public boolean incomingNeedsProcessing() { public Edge getIncoming() {
return incomingNeedsProcessing; return incoming;
} }
/** Check outgoing edge processed indicator. /** Set outgoing edge.
* @return true if outgoing edge needs processing * <p>
* The line supporting the outgoing edge is automatically bound
* with the instance.
* </p>
* @param incoming outgoing edge
*/ */
public boolean outgoingNeedsProcessing() { public void setOutgoing(final Edge outgoing) {
return outgoingNeedsProcessing; this.outgoing = outgoing;
bindWith(outgoing.getLine());
} }
/** Mark the incoming edge as processed. /** Get outgoing edge.
* @return outgoing edge
*/ */
public void incomingProcessed() { public Edge getOutgoing() {
incomingNeedsProcessing = false; return outgoing;
} }
/** Mark the outgoing edge as processed. }
/** Internal class for holding edges while they are processed to build a BSP tree. */
private static class Edge {
/** Start vertex. */
private final Vertex start;
/** End vertex. */
private final Vertex end;
/** Line supporting the edge. */
private final Line line;
/** Node whose cut hyperplane contains this edge. */
private BSPTree<Euclidean2D> node;
/** Build an edge not contained in any node yet.
* @param start start vertex
* @param end end vertex
* @param line line supporting the edge
*/ */
public void outgoingProcessed() { public Edge(final Vertex start, final Vertex end, final Line line) {
outgoingNeedsProcessing = false;
this.start = start;
this.end = end;
this.line = line;
this.node = null;
// connect the vertices back to the edge
start.setOutgoing(this);
end.setIncoming(this);
}
/** Get start vertex.
* @return start vertex
*/
public Vertex getStart() {
return start;
}
/** Get end vertex.
* @return end vertex
*/
public Vertex getEnd() {
return end;
}
/** Get the line supporting this edge.
* @return line supporting this edge
*/
public Line getLine() {
return line;
}
/** Set the node whose cut hyperplane contains this edge.
* @param node node whose cut hyperplane contains this edge
*/
public void setNode(final BSPTree<Euclidean2D> node) {
this.node = node;
}
/** Get the node whose cut hyperplane contains this edge.
* @return node whose cut hyperplane contains this edge
* (null if edge has not yet been inserted into the BSP tree)
*/
public BSPTree<Euclidean2D> getNode() {
return node;
}
/** Split the edge.
* <p>
* Once split, this edge is not referenced anymore by the vertices,
* it is replaced by the two half-edges and an intermediate splitting
* vertex is introduced to connect these two halves.
* </p>
* @param splitLine line splitting the edge in two halves
* @return split vertex (its incoming and outgoing edges are the two halves)
*/
public Vertex split(final Line splitLine) {
final Vertex splitVertex = new Vertex(line.intersection(splitLine));
splitVertex.bindWith(splitLine);
final Edge startHalf = new Edge(start, splitVertex, line);
final Edge endHalf = new Edge(splitVertex, end, line);
startHalf.node = node;
endHalf.node = node;
return splitVertex;
} }
} }

32
src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
View File

@ -16,7 +16,7 @@
*/ */
package org.apache.commons.math3.geometry.euclidean.twod; package org.apache.commons.math3.geometry.euclidean.twod;
import java.io.FileNotFoundException; import java.io.IOException;
import java.util.ArrayList; import java.util.ArrayList;
import java.util.List; import java.util.List;
@ -801,6 +801,14 @@ public class PolygonsSetTest {
} }
@Test
public void testSqueezedHexa() {
PolygonsSet set = new PolygonsSet(1.0e-10,
new Vector2D(-6, -4), new Vector2D(-8, -8), new Vector2D( 8, -8),
new Vector2D( 6, -4), new Vector2D(10, 4), new Vector2D(-10, 4));
Assert.assertEquals(Location.OUTSIDE, set.checkPoint(new Vector2D(0, 6)));
}
@Test @Test
public void testIssue880Simplified() { public void testIssue880Simplified() {
@ -821,7 +829,7 @@ public class PolygonsSetTest {
} }
@Test @Test
public void testIssue880Complete() throws FileNotFoundException { public void testIssue880Complete() throws IOException {
Vector2D[] vertices1 = new Vector2D[] { Vector2D[] vertices1 = new Vector2D[] {
new Vector2D( 90.08714908223715, 38.370299337260235), new Vector2D( 90.08714908223715, 38.370299337260235),
new Vector2D( 90.08709517675004, 38.3702895991413), new Vector2D( 90.08709517675004, 38.3702895991413),
@ -894,12 +902,13 @@ public class PolygonsSetTest {
new Vector2D( 90.09081227075944, 38.37526295920463), new Vector2D( 90.09081227075944, 38.37526295920463),
new Vector2D( 90.09081378927135, 38.375193883266434) new Vector2D( 90.09081378927135, 38.375193883266434)
}; };
PolygonsSet set1 = new PolygonsSet(1.0e-10, vertices1); PolygonsSet set1 = new PolygonsSet(1.0e-8, vertices1);
Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.0905, 38.3755))); Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.0905, 38.3755)));
Assert.assertEquals(Location.INSIDE, set1.checkPoint(new Vector2D(90.09084, 38.3755))); Assert.assertEquals(Location.INSIDE, set1.checkPoint(new Vector2D(90.09084, 38.3755)));
// TODO: the following assertion fails and should not Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.0913, 38.3755)));
// this is due to a small spurious triangle being included in the polygon Assert.assertEquals(Location.INSIDE, set1.checkPoint(new Vector2D(90.1042, 38.3739)));
// Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.0913, 38.3755))); Assert.assertEquals(Location.INSIDE, set1.checkPoint(new Vector2D(90.1111, 38.3673)));
Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.0959, 38.3457)));
Vector2D[] vertices2 = new Vector2D[] { Vector2D[] vertices2 = new Vector2D[] {
new Vector2D( 90.13067558880044, 38.36977255037573), new Vector2D( 90.13067558880044, 38.36977255037573),
@ -965,17 +974,14 @@ public class PolygonsSetTest {
new Vector2D( 90.16746107640665, 38.40902614307544), new Vector2D( 90.16746107640665, 38.40902614307544),
new Vector2D( 90.16122795307462, 38.39773101873203) new Vector2D( 90.16122795307462, 38.39773101873203)
}; };
PolygonsSet set2 = new PolygonsSet(1.0e-10, vertices2); PolygonsSet set2 = new PolygonsSet(1.0e-8, vertices2);
PolygonsSet set = (PolygonsSet) new PolygonsSet set = (PolygonsSet) new
RegionFactory<Euclidean2D>().difference(set1.copySelf(), RegionFactory<Euclidean2D>().difference(set1.copySelf(),
set2.copySelf()); set2.copySelf());
Vector2D[][] vertices = set.getVertices(); Vector2D[][] verticies = set.getVertices();
Assert.assertTrue(vertices[0][0] != null); Assert.assertTrue(verticies[0][0] != null);
// TODO: the resulting polygon has two boundaries but should have only one Assert.assertEquals(1, verticies.length);
// this is because for an unknown reason the boundary has two infinitely close
// parallel paths near the top left of the polygon
Assert.assertEquals(2, vertices.length);
} }
private PolygonsSet buildSet(Vector2D[][] vertices) { private PolygonsSet buildSet(Vector2D[][] vertices) {