Improved polygons creation with a numerically more robust constructor.

This is only a partial fix for issue 880 because there are still some
glitches in the first polygon involved in the issue.

JIRA: MATH-880

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1400717 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2012-10-21 20:08:32 +00:00
parent c323480cf6
commit b25f448d2f
4 changed files with 449 additions and 7 deletions

View File

@ -52,13 +52,15 @@ If the output is not quite correct, check for invisible trailing spaces!
<body>
<release version="3.1" date="TBD" description="
">
<action dev="luc" type="fix" issue="MATH-880">
Improved construction of polygons with an additional constructor, more robust numerically.
</action>
<action dev="tn" type="add" issue="MATH-474" due-to="Dan Checkoway">
Added new methods "merge(Frequency)", "merge(Collection&lt;Frequency&gt;)",
"incrementValue(Comparable&lt;?&gt;, long)" and "entrySetIterator()" to the "Frequency" class.
</action>
<action dev="tn" type="fix" issue="MATH-778" due-to="Sébastien Brisard">
Allow unlimited input values for "Dfp#multiply(int)".
</action>
<action dev="luc" type="fix" issue="MATH-641" due-to="Curtis Jensen">
Added distance to point to 2D Line and Segment.
</action>

View File

@ -25,12 +25,13 @@ import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D;
import org.apache.commons.math3.geometry.euclidean.oned.Interval;
import org.apache.commons.math3.geometry.euclidean.oned.IntervalsSet;
import org.apache.commons.math3.geometry.euclidean.oned.Vector1D;
import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane;
import org.apache.commons.math3.geometry.partitioning.BSPTree;
import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
import org.apache.commons.math3.geometry.partitioning.Side;
import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
import org.apache.commons.math3.geometry.partitioning.utilities.AVLTree;
import org.apache.commons.math3.geometry.partitioning.utilities.OrderedTuple;
import org.apache.commons.math3.util.FastMath;
@ -98,6 +99,26 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
super(boxBoundary(xMin, xMax, yMin, yMax));
}
/** Build a polygon 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
* region is on the left side of this path and the exterior is on its
* right side.</p>
* <p>This constructor does not handle polygons with a boundary
* forming several disconnected paths (such as polygons with holes).</p>
* <p>For cases where this simple constructor applies, it is expected to
* be numerically more robust than the {@link #PolygonsSet(Collection) general
* constructor} using {@link SubHyperplane subhyperplanes}.</p>
* <p>If the list is empty, the region will represent the whole
* space.</p>
* @param hyperplaneThickness tolerance below which points are considered to
* belong to the hyperplane (which is therefore more a slab)
* @param vertices vertices of the simple loop boundary
*/
public PolygonsSet(final double hyperplaneThickness, final Vector2D ... vertices) {
super(verticesToTree(hyperplaneThickness, vertices));
}
/** Create a list of hyperplanes representing the boundary of a box.
* @param xMin low bound along the x direction
* @param xMax high bound along the x direction
@ -119,6 +140,249 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
};
}
/** 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
* region is on the left side of this path and the exterior is on its
* right side.</p>
* <p>This constructor does not handle polygons with a boundary
* forming several disconnected paths (such as polygons with holes).</p>
* <p>For cases where this simple constructor applies, it is expected to
* be numerically more robust than the {@link #PolygonsSet(Collection) general
* constructor} using {@link SubHyperplane subhyperplanes}.</p>
* @param hyperplaneThickness tolerance below which points are consider to
* belong to the hyperplane (which is therefore more a slab)
* @param vertices vertices of the simple loop boundary
*/
private static BSPTree<Euclidean2D> verticesToTree(final double hyperplaneThickness,
final Vector2D ... vertices) {
if (vertices.length == 0) {
// the tree represents the whole space
return new BSPTree<Euclidean2D>(Boolean.TRUE);
}
// at start, none of the edges have been processed
final BSPTree<Euclidean2D> tree = new BSPTree<Euclidean2D>();
List<Vertex> list = new ArrayList<PolygonsSet.Vertex>(vertices.length);
for (final Vector2D vertex : vertices) {
list.add(new Vertex(vertex));
}
// build the tree top-down
insertVertices(hyperplaneThickness, tree, list);
return tree;
}
/** Recursively build a tree by inserting cut sub-hyperplanes.
* @param hyperplaneThickness tolerance below which points are consider to
* belong to the hyperplane (which is therefore more a slab)
* @param node current tree node (it is a leaf node at the beginning
* of the call)
* @param vertices list of vertices belonging to the boundary of the
* cell defined by the node
*/
private static void insertVertices(final double hyperplaneThickness,
final BSPTree<Euclidean2D> node,
final List<Vertex> vertices) {
Vertex current = vertices.get(vertices.size() - 1);
int index = 0;
Line inserted = null;
while (inserted == null && index < vertices.size()) {
final Vertex previous = current;
current = vertices.get(index++);
if (previous.outgoingNeedsProcessing() && current.incomingNeedsProcessing()) {
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 {
inserted = null;
}
}
}
}
if (node.getCut() == null) {
final BSPTree<Euclidean2D> parent = node.getParent();
if (parent == null || node == parent.getMinus()) {
node.setAttribute(Boolean.TRUE);
} else {
node.setAttribute(Boolean.FALSE);
}
return;
}
// distribute the remaining vertices in the two sub-trees
Side currentSide = Side.HYPER;
final List<Vertex> plusList = new ArrayList<Vertex>();
plusList.add(current);
int plusCount = 0;
final List<Vertex> minusList = new ArrayList<Vertex>();
minusList.add(current);
int minusCount = 0;
while (index < vertices.size()) {
final Vertex previous = current;
final Side previousSide = currentSide;
current = vertices.get(index++);
final double currentOffset = inserted.getOffset(current.getLocation());
currentSide = (FastMath.abs(currentOffset) <= hyperplaneThickness) ?
Side.HYPER :
((currentOffset < 0) ? Side.MINUS : Side.PLUS);
switch (currentSide) {
case PLUS:
if (previousSide == Side.MINUS) {
// we need to insert a split point on the hyperplane
final Line line = new Line(previous.getLocation(), current.getLocation());
final Vertex splitPoint = new Vertex(inserted.intersection(line));
splitPoint.addNode(node);
minusList.add(splitPoint);
plusList.add(splitPoint);
}
plusList.add(current);
if (current.incomingNeedsProcessing() || current.outgoingNeedsProcessing()) {
++plusCount;
}
break;
case MINUS:
if (previousSide == Side.PLUS) {
// we need to insert a split point on the hyperplane
final Line line = new Line(previous.getLocation(), current.getLocation());
final Vertex splitPoint = new Vertex(inserted.intersection(line));
splitPoint.addNode(node);
minusList.add(splitPoint);
plusList.add(splitPoint);
}
minusList.add(current);
if (current.incomingNeedsProcessing() || current.outgoingNeedsProcessing()) {
++minusCount;
}
break;
default:
current.addNode(node);
plusList.add(current);
minusList.add(current);
break;
}
}
// recurse through lower levels
if (plusCount > 0) {
insertVertices(hyperplaneThickness, node.getPlus(), plusList);
} else {
node.getPlus().setAttribute(Boolean.FALSE);
}
if (minusCount > 0) {
insertVertices(hyperplaneThickness, node.getMinus(), minusList);
} else {
node.getMinus().setAttribute(Boolean.TRUE);
}
}
/** Internal class for holding vertices while they are processed to build a BSP tree. */
private static class Vertex {
/** Vertex location. */
private final Vector2D location;
/** Nodes associated with the hyperplane containing this vertex. */
private final List<BSPTree<Euclidean2D>> nodes;
/** Indicator for incoming edges that still need processing. */
private boolean incomingNeedsProcessing;
/** Indicator for outgoing edges that still need processing. */
private boolean outgoingNeedsProcessing;
/** Build a non-processed vertex not owned by any node yet.
* @param location vertex location
*/
public Vertex(final Vector2D location) {
this.location = location;
this.nodes = new ArrayList<BSPTree<Euclidean2D>>();
this.incomingNeedsProcessing = true;
this.outgoingNeedsProcessing = true;
}
/** Get Vertex location.
* @return vertex location
*/
public Vector2D getLocation() {
return location;
}
/** Check if the instance and another vertex share a node.
* <p>
* When two vertices share a node, this means they are already handled
* by the hyperplane of this node, so there is no need to create a cut
* hyperplane for them.
* </p>
* @param vertex other vertex to check instance against
* @return true if the instance and another vertex share a node
*/
public boolean shareNodeWith(final Vertex vertex) {
for (final BSPTree<Euclidean2D> node1 : nodes) {
for (final BSPTree<Euclidean2D> node2 : vertex.nodes) {
if (node1 == node2) {
return true;
}
}
}
return false;
}
/** Add a node whose hyperplane contains this vertex.
* @param node node whose hyperplane contains this vertex
*/
public void addNode(final BSPTree<Euclidean2D> node) {
nodes.add(node);
}
/** Check incoming edge processed indicator.
* @return true if incoming edge needs processing
*/
public boolean incomingNeedsProcessing() {
return incomingNeedsProcessing;
}
/** Check outgoing edge processed indicator.
* @return true if outgoing edge needs processing
*/
public boolean outgoingNeedsProcessing() {
return outgoingNeedsProcessing;
}
/** Mark the incoming edge as processed.
*/
public void incomingProcessed() {
incomingNeedsProcessing = false;
}
/** Mark the outgoing edge as processed.
*/
public void outgoingProcessed() {
outgoingNeedsProcessing = false;
}
}
/** {@inheritDoc} */
@Override
public PolygonsSet buildNew(final BSPTree<Euclidean2D> tree) {

View File

@ -153,7 +153,7 @@ public class BSPTree<S extends Space> {
}
final SubHyperplane<S> chopped = fitToCell(hyperplane.wholeHyperplane());
if (chopped.isEmpty()) {
if (chopped == null || chopped.isEmpty()) {
cut = null;
plus = null;
minus = null;
@ -285,7 +285,7 @@ public class BSPTree<S extends Space> {
* sub-hyperplane that lie outside of the cell using the
* cut-hyperplanes of the parent nodes of the instance.</p>
* @param sub sub-hyperplane to fit
* @return a new sub-hyperplane, gueranteed to have no part outside
* @return a new sub-hyperplane, guaranteed to have no part outside
* of the instance cell
*/
private SubHyperplane<S> fitToCell(final SubHyperplane<S> sub) {

View File

@ -16,17 +16,16 @@
*/
package org.apache.commons.math3.geometry.euclidean.twod;
import java.io.FileNotFoundException;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math3.geometry.euclidean.oned.Interval;
import org.apache.commons.math3.geometry.euclidean.oned.IntervalsSet;
import org.apache.commons.math3.geometry.euclidean.oned.Vector1D;
import org.apache.commons.math3.geometry.euclidean.twod.Line;
import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
import org.apache.commons.math3.geometry.euclidean.twod.PolygonsSet;
import org.apache.commons.math3.geometry.partitioning.BSPTree;
import org.apache.commons.math3.geometry.partitioning.Region;
import org.apache.commons.math3.geometry.partitioning.Region.Location;
import org.apache.commons.math3.geometry.partitioning.RegionFactory;
import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
import org.apache.commons.math3.util.FastMath;
@ -802,6 +801,183 @@ public class PolygonsSetTest {
}
@Test
public void testIssue880Simplified() {
Vector2D[] vertices1 = new Vector2D[] {
new Vector2D( 90.13595870833188, 38.33604606376991),
new Vector2D( 90.14047850603913, 38.34600084496253),
new Vector2D( 90.11045289492762, 38.36801537312368),
new Vector2D( 90.10871471476526, 38.36878044144294),
new Vector2D( 90.10424901707671, 38.374300101757),
new Vector2D( 90.0979455456843, 38.373578376172475),
new Vector2D( 90.09081227075944, 38.37526295920463),
new Vector2D( 90.09081378927135, 38.375193883266434)
};
PolygonsSet set1 = new PolygonsSet(1.0e-10, vertices1);
Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.12, 38.32)));
Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.135, 38.355)));
}
@Test
public void testIssue880Complete() throws FileNotFoundException {
Vector2D[] vertices1 = new Vector2D[] {
new Vector2D( 90.08714908223715, 38.370299337260235),
new Vector2D( 90.08709517675004, 38.3702895991413),
new Vector2D( 90.08401538704919, 38.368849330127944),
new Vector2D( 90.08258210430711, 38.367634558585564),
new Vector2D( 90.08251455106665, 38.36763409247078),
new Vector2D( 90.08106599752608, 38.36761621664249),
new Vector2D( 90.08249585300035, 38.36753627557965),
new Vector2D( 90.09075743352184, 38.35914647644972),
new Vector2D( 90.09099945896571, 38.35896264724079),
new Vector2D( 90.09269383800086, 38.34595756121246),
new Vector2D( 90.09638631543191, 38.3457988093121),
new Vector2D( 90.09666417351019, 38.34523360999418),
new Vector2D( 90.1297082145872, 38.337670454923625),
new Vector2D( 90.12971687748956, 38.337669827794684),
new Vector2D( 90.1240820219179, 38.34328502001131),
new Vector2D( 90.13084259656404, 38.34017811765017),
new Vector2D( 90.13378567942857, 38.33860579180606),
new Vector2D( 90.13519557833206, 38.33621054663689),
new Vector2D( 90.13545616732307, 38.33614965452864),
new Vector2D( 90.13553111202748, 38.33613962818305),
new Vector2D( 90.1356903436448, 38.33610227127048),
new Vector2D( 90.13576283227428, 38.33609255422783),
new Vector2D( 90.13595870833188, 38.33604606376991),
new Vector2D( 90.1361556630693, 38.3360024198866),
new Vector2D( 90.13622408795709, 38.335987048115726),
new Vector2D( 90.13696189099994, 38.33581914328681),
new Vector2D( 90.13746655304897, 38.33616706665265),
new Vector2D( 90.13845973716064, 38.33650776167099),
new Vector2D( 90.13950901827667, 38.3368469456463),
new Vector2D( 90.14393814424852, 38.337591835857495),
new Vector2D( 90.14483839716831, 38.337076122362475),
new Vector2D( 90.14565474433601, 38.33769000964429),
new Vector2D( 90.14569421179482, 38.3377117256905),
new Vector2D( 90.14577067124333, 38.33770883625908),
new Vector2D( 90.14600350631684, 38.337714326520995),
new Vector2D( 90.14600355139731, 38.33771435193319),
new Vector2D( 90.14600369112401, 38.33771443882085),
new Vector2D( 90.14600382486884, 38.33771453466096),
new Vector2D( 90.14600395205912, 38.33771463904344),
new Vector2D( 90.14600407214999, 38.337714751520764),
new Vector2D( 90.14600418462749, 38.337714871611695),
new Vector2D( 90.14600422249327, 38.337714915811034),
new Vector2D( 90.14867838361471, 38.34113888210675),
new Vector2D( 90.14923750157374, 38.341582537502575),
new Vector2D( 90.14877083250991, 38.34160685841391),
new Vector2D( 90.14816667319519, 38.34244232585684),
new Vector2D( 90.14797696744586, 38.34248455284745),
new Vector2D( 90.14484318014337, 38.34385573215269),
new Vector2D( 90.14477919958296, 38.3453797747614),
new Vector2D( 90.14202393306448, 38.34464324839456),
new Vector2D( 90.14198920640195, 38.344651155237216),
new Vector2D( 90.14155207025175, 38.34486424263724),
new Vector2D( 90.1415196143314, 38.344871730519),
new Vector2D( 90.14128611910814, 38.34500196593859),
new Vector2D( 90.14047850603913, 38.34600084496253),
new Vector2D( 90.14045907000337, 38.34601860032171),
new Vector2D( 90.14039496493928, 38.346223030432384),
new Vector2D( 90.14037626063737, 38.346240203360026),
new Vector2D( 90.14030005823724, 38.34646920000705),
new Vector2D( 90.13799164754806, 38.34903093011013),
new Vector2D( 90.11045289492762, 38.36801537312368),
new Vector2D( 90.10871471476526, 38.36878044144294),
new Vector2D( 90.10424901707671, 38.374300101757),
new Vector2D( 90.10263482039932, 38.37310041316073),
new Vector2D( 90.09834601753448, 38.373615053823414),
new Vector2D( 90.0979455456843, 38.373578376172475),
new Vector2D( 90.09086514328669, 38.37527884194668),
new Vector2D( 90.09084931407364, 38.37590801712463),
new Vector2D( 90.09081227075944, 38.37526295920463),
new Vector2D( 90.09081378927135, 38.375193883266434)
};
PolygonsSet set1 = new PolygonsSet(1.0e-10, vertices1);
Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.0905, 38.3755)));
Assert.assertEquals(Location.INSIDE, set1.checkPoint(new Vector2D(90.09084, 38.3755)));
// TODO: the following assertion fails and should not
// this is due to a small spurious triangle being included in the polygon
// Assert.assertEquals(Location.OUTSIDE, set1.checkPoint(new Vector2D(90.0913, 38.3755)));
Vector2D[] vertices2 = new Vector2D[] {
new Vector2D( 90.13067558880044, 38.36977255037573),
new Vector2D( 90.12907570488, 38.36817308242706),
new Vector2D( 90.1342774136516, 38.356886880294724),
new Vector2D( 90.13090330629757, 38.34664392676211),
new Vector2D( 90.13078571364593, 38.344904617518466),
new Vector2D( 90.1315602208914, 38.3447185040846),
new Vector2D( 90.1316336226821, 38.34470643148342),
new Vector2D( 90.134020944832, 38.340936644972885),
new Vector2D( 90.13912536387306, 38.335497255122334),
new Vector2D( 90.1396178806582, 38.334878075552126),
new Vector2D( 90.14083049696671, 38.33316530644106),
new Vector2D( 90.14145252901329, 38.33152722916191),
new Vector2D( 90.1404779335565, 38.32863516047786),
new Vector2D( 90.14282712131586, 38.327504432532066),
new Vector2D( 90.14616669875488, 38.3237354115015),
new Vector2D( 90.14860976050608, 38.315714862457924),
new Vector2D( 90.14999277782437, 38.3164932507504),
new Vector2D( 90.15005207194997, 38.316534677663356),
new Vector2D( 90.15508513859612, 38.31878731691609),
new Vector2D( 90.15919938519221, 38.31852743183782),
new Vector2D( 90.16093758658837, 38.31880662005153),
new Vector2D( 90.16099420184912, 38.318825953291594),
new Vector2D( 90.1665411125756, 38.31859497874757),
new Vector2D( 90.16999653861313, 38.32505772048029),
new Vector2D( 90.17475243391698, 38.32594398441148),
new Vector2D( 90.17940844844992, 38.327427213761325),
new Vector2D( 90.20951909541378, 38.330616833491774),
new Vector2D( 90.2155400467941, 38.331746223670336),
new Vector2D( 90.21559881391778, 38.33175551425302),
new Vector2D( 90.21916646426041, 38.332584299620805),
new Vector2D( 90.23863749852285, 38.34778978875795),
new Vector2D( 90.25459855175802, 38.357790570608984),
new Vector2D( 90.25964298227257, 38.356918010203174),
new Vector2D( 90.26024593994703, 38.361692743151366),
new Vector2D( 90.26146187570015, 38.36311080550837),
new Vector2D( 90.26614159359622, 38.36510808579902),
new Vector2D( 90.26621342936448, 38.36507942500333),
new Vector2D( 90.26652190211962, 38.36494042196722),
new Vector2D( 90.26621240678867, 38.365113172030874),
new Vector2D( 90.26614057102057, 38.365141832826794),
new Vector2D( 90.26380080055299, 38.3660381760273),
new Vector2D( 90.26315345241, 38.36670658276421),
new Vector2D( 90.26251574942881, 38.367490323488084),
new Vector2D( 90.26247873448426, 38.36755266444749),
new Vector2D( 90.26234628016698, 38.36787989125406),
new Vector2D( 90.26214559424784, 38.36945909356126),
new Vector2D( 90.25861728442555, 38.37200753430875),
new Vector2D( 90.23905557537864, 38.375405314295904),
new Vector2D( 90.22517251874075, 38.38984691662256),
new Vector2D( 90.22549955153215, 38.3911564273979),
new Vector2D( 90.22434386063355, 38.391476432092134),
new Vector2D( 90.22147729457276, 38.39134652252034),
new Vector2D( 90.22142070120117, 38.391349167741964),
new Vector2D( 90.20665060751588, 38.39475580900313),
new Vector2D( 90.20042268367109, 38.39842558622888),
new Vector2D( 90.17423771242085, 38.402727751805344),
new Vector2D( 90.16756796257476, 38.40913898597597),
new Vector2D( 90.16728283954308, 38.411255399912875),
new Vector2D( 90.16703538220418, 38.41136059866693),
new Vector2D( 90.16725865657685, 38.41013618805954),
new Vector2D( 90.16746107640665, 38.40902614307544),
new Vector2D( 90.16122795307462, 38.39773101873203)
};
PolygonsSet set2 = new PolygonsSet(1.0e-10, vertices2);
PolygonsSet set = (PolygonsSet) new
RegionFactory<Euclidean2D>().difference(set1.copySelf(),
set2.copySelf());
Vector2D[][] vertices = set.getVertices();
Assert.assertTrue(vertices[0][0] != null);
// TODO: the resulting polygon has two boundaries but should have only one
// 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) {
ArrayList<SubHyperplane<Euclidean2D>> edges = new ArrayList<SubHyperplane<Euclidean2D>>();
for (int i = 0; i < vertices.length; ++i) {