1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.math3.geometry.euclidean.twod; 18 19 import java.util.ArrayList; 20 import java.util.Collection; 21 import java.util.List; 22 23 import org.apache.commons.math3.geometry.Point; 24 import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D; 25 import org.apache.commons.math3.geometry.euclidean.oned.Interval; 26 import org.apache.commons.math3.geometry.euclidean.oned.IntervalsSet; 27 import org.apache.commons.math3.geometry.euclidean.oned.Vector1D; 28 import org.apache.commons.math3.geometry.partitioning.AbstractRegion; 29 import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane; 30 import org.apache.commons.math3.geometry.partitioning.BSPTree; 31 import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor; 32 import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute; 33 import org.apache.commons.math3.geometry.partitioning.Hyperplane; 34 import org.apache.commons.math3.geometry.partitioning.Side; 35 import org.apache.commons.math3.geometry.partitioning.SubHyperplane; 36 import org.apache.commons.math3.util.FastMath; 37 import org.apache.commons.math3.util.Precision; 38 39 /** This class represents a 2D region: a set of polygons. 40 * @since 3.0 41 */ 42 public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> { 43 44 /** Default value for tolerance. */ 45 private static final double DEFAULT_TOLERANCE = 1.0e-10; 46 47 /** Vertices organized as boundary loops. */ 48 private Vector2D[][] vertices; 49 50 /** Build a polygons set representing the whole plane. 51 * @param tolerance tolerance below which points are considered identical 52 * @since 3.3 53 */ PolygonsSet(final double tolerance)54 public PolygonsSet(final double tolerance) { 55 super(tolerance); 56 } 57 58 /** Build a polygons set from a BSP tree. 59 * <p>The leaf nodes of the BSP tree <em>must</em> have a 60 * {@code Boolean} attribute representing the inside status of 61 * the corresponding cell (true for inside cells, false for outside 62 * cells). In order to avoid building too many small objects, it is 63 * recommended to use the predefined constants 64 * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p> 65 * <p> 66 * This constructor is aimed at expert use, as building the tree may 67 * be a difficult task. It is not intended for general use and for 68 * performances reasons does not check thoroughly its input, as this would 69 * require walking the full tree each time. Failing to provide a tree with 70 * the proper attributes, <em>will</em> therefore generate problems like 71 * {@link NullPointerException} or {@link ClassCastException} only later on. 72 * This limitation is known and explains why this constructor is for expert 73 * use only. The caller does have the responsibility to provided correct arguments. 74 * </p> 75 * @param tree inside/outside BSP tree representing the region 76 * @param tolerance tolerance below which points are considered identical 77 * @since 3.3 78 */ PolygonsSet(final BSPTree<Euclidean2D> tree, final double tolerance)79 public PolygonsSet(final BSPTree<Euclidean2D> tree, final double tolerance) { 80 super(tree, tolerance); 81 } 82 83 /** Build a polygons set from a Boundary REPresentation (B-rep). 84 * <p>The boundary is provided as a collection of {@link 85 * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the 86 * interior part of the region on its minus side and the exterior on 87 * its plus side.</p> 88 * <p>The boundary elements can be in any order, and can form 89 * several non-connected sets (like for example polygons with holes 90 * or a set of disjoint polygons considered as a whole). In 91 * fact, the elements do not even need to be connected together 92 * (their topological connections are not used here). However, if the 93 * boundary does not really separate an inside open from an outside 94 * open (open having here its topological meaning), then subsequent 95 * calls to the {@link 96 * org.apache.commons.math3.geometry.partitioning.Region#checkPoint(org.apache.commons.math3.geometry.Point) 97 * checkPoint} method will not be meaningful anymore.</p> 98 * <p>If the boundary is empty, the region will represent the whole 99 * space.</p> 100 * @param boundary collection of boundary elements, as a 101 * collection of {@link SubHyperplane SubHyperplane} objects 102 * @param tolerance tolerance below which points are considered identical 103 * @since 3.3 104 */ PolygonsSet(final Collection<SubHyperplane<Euclidean2D>> boundary, final double tolerance)105 public PolygonsSet(final Collection<SubHyperplane<Euclidean2D>> boundary, final double tolerance) { 106 super(boundary, tolerance); 107 } 108 109 /** Build a parallellepipedic box. 110 * @param xMin low bound along the x direction 111 * @param xMax high bound along the x direction 112 * @param yMin low bound along the y direction 113 * @param yMax high bound along the y direction 114 * @param tolerance tolerance below which points are considered identical 115 * @since 3.3 116 */ PolygonsSet(final double xMin, final double xMax, final double yMin, final double yMax, final double tolerance)117 public PolygonsSet(final double xMin, final double xMax, 118 final double yMin, final double yMax, 119 final double tolerance) { 120 super(boxBoundary(xMin, xMax, yMin, yMax, tolerance), tolerance); 121 } 122 123 /** Build a polygon from a simple list of vertices. 124 * <p>The boundary is provided as a list of points considering to 125 * represent the vertices of a simple loop. The interior part of the 126 * region is on the left side of this path and the exterior is on its 127 * right side.</p> 128 * <p>This constructor does not handle polygons with a boundary 129 * forming several disconnected paths (such as polygons with holes).</p> 130 * <p>For cases where this simple constructor applies, it is expected to 131 * be numerically more robust than the {@link #PolygonsSet(Collection) general 132 * constructor} using {@link SubHyperplane subhyperplanes}.</p> 133 * <p>If the list is empty, the region will represent the whole 134 * space.</p> 135 * <p> 136 * Polygons with thin pikes or dents are inherently difficult to handle because 137 * they involve lines with almost opposite directions at some vertices. Polygons 138 * whose vertices come from some physical measurement with noise are also 139 * difficult because an edge that should be straight may be broken in lots of 140 * different pieces with almost equal directions. In both cases, computing the 141 * lines intersections is not numerically robust due to the almost 0 or almost 142 * π angle. Such cases need to carefully adjust the {@code hyperplaneThickness} 143 * parameter. A too small value would often lead to completely wrong polygons 144 * with large area wrongly identified as inside or outside. Large values are 145 * often much safer. As a rule of thumb, a value slightly below the size of the 146 * most accurate detail needed is a good value for the {@code hyperplaneThickness} 147 * parameter. 148 * </p> 149 * @param hyperplaneThickness tolerance below which points are considered to 150 * belong to the hyperplane (which is therefore more a slab) 151 * @param vertices vertices of the simple loop boundary 152 */ PolygonsSet(final double hyperplaneThickness, final Vector2D ... vertices)153 public PolygonsSet(final double hyperplaneThickness, final Vector2D ... vertices) { 154 super(verticesToTree(hyperplaneThickness, vertices), hyperplaneThickness); 155 } 156 157 /** Build a polygons set representing the whole real line. 158 * @deprecated as of 3.3, replaced with {@link #PolygonsSet(double)} 159 */ 160 @Deprecated PolygonsSet()161 public PolygonsSet() { 162 this(DEFAULT_TOLERANCE); 163 } 164 165 /** Build a polygons set from a BSP tree. 166 * <p>The leaf nodes of the BSP tree <em>must</em> have a 167 * {@code Boolean} attribute representing the inside status of 168 * the corresponding cell (true for inside cells, false for outside 169 * cells). In order to avoid building too many small objects, it is 170 * recommended to use the predefined constants 171 * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p> 172 * @param tree inside/outside BSP tree representing the region 173 * @deprecated as of 3.3, replaced with {@link #PolygonsSet(BSPTree, double)} 174 */ 175 @Deprecated PolygonsSet(final BSPTree<Euclidean2D> tree)176 public PolygonsSet(final BSPTree<Euclidean2D> tree) { 177 this(tree, DEFAULT_TOLERANCE); 178 } 179 180 /** Build a polygons set from a Boundary REPresentation (B-rep). 181 * <p>The boundary is provided as a collection of {@link 182 * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the 183 * interior part of the region on its minus side and the exterior on 184 * its plus side.</p> 185 * <p>The boundary elements can be in any order, and can form 186 * several non-connected sets (like for example polygons with holes 187 * or a set of disjoint polygons considered as a whole). In 188 * fact, the elements do not even need to be connected together 189 * (their topological connections are not used here). However, if the 190 * boundary does not really separate an inside open from an outside 191 * open (open having here its topological meaning), then subsequent 192 * calls to the {@link 193 * org.apache.commons.math3.geometry.partitioning.Region#checkPoint(org.apache.commons.math3.geometry.Point) 194 * checkPoint} method will not be meaningful anymore.</p> 195 * <p>If the boundary is empty, the region will represent the whole 196 * space.</p> 197 * @param boundary collection of boundary elements, as a 198 * collection of {@link SubHyperplane SubHyperplane} objects 199 * @deprecated as of 3.3, replaced with {@link #PolygonsSet(Collection, double)} 200 */ 201 @Deprecated PolygonsSet(final Collection<SubHyperplane<Euclidean2D>> boundary)202 public PolygonsSet(final Collection<SubHyperplane<Euclidean2D>> boundary) { 203 this(boundary, DEFAULT_TOLERANCE); 204 } 205 206 /** Build a parallellepipedic box. 207 * @param xMin low bound along the x direction 208 * @param xMax high bound along the x direction 209 * @param yMin low bound along the y direction 210 * @param yMax high bound along the y direction 211 * @deprecated as of 3.3, replaced with {@link #PolygonsSet(double, double, double, double, double)} 212 */ 213 @Deprecated PolygonsSet(final double xMin, final double xMax, final double yMin, final double yMax)214 public PolygonsSet(final double xMin, final double xMax, 215 final double yMin, final double yMax) { 216 this(xMin, xMax, yMin, yMax, DEFAULT_TOLERANCE); 217 } 218 219 /** Create a list of hyperplanes representing the boundary of a box. 220 * @param xMin low bound along the x direction 221 * @param xMax high bound along the x direction 222 * @param yMin low bound along the y direction 223 * @param yMax high bound along the y direction 224 * @param tolerance tolerance below which points are considered identical 225 * @return boundary of the box 226 */ boxBoundary(final double xMin, final double xMax, final double yMin, final double yMax, final double tolerance)227 private static Line[] boxBoundary(final double xMin, final double xMax, 228 final double yMin, final double yMax, 229 final double tolerance) { 230 if ((xMin >= xMax - tolerance) || (yMin >= yMax - tolerance)) { 231 // too thin box, build an empty polygons set 232 return null; 233 } 234 final Vector2D minMin = new Vector2D(xMin, yMin); 235 final Vector2D minMax = new Vector2D(xMin, yMax); 236 final Vector2D maxMin = new Vector2D(xMax, yMin); 237 final Vector2D maxMax = new Vector2D(xMax, yMax); 238 return new Line[] { 239 new Line(minMin, maxMin, tolerance), 240 new Line(maxMin, maxMax, tolerance), 241 new Line(maxMax, minMax, tolerance), 242 new Line(minMax, minMin, tolerance) 243 }; 244 } 245 246 /** Build the BSP tree of a polygons set from a simple list of vertices. 247 * <p>The boundary is provided as a list of points considering to 248 * represent the vertices of a simple loop. The interior part of the 249 * region is on the left side of this path and the exterior is on its 250 * right side.</p> 251 * <p>This constructor does not handle polygons with a boundary 252 * forming several disconnected paths (such as polygons with holes).</p> 253 * <p>For cases where this simple constructor applies, it is expected to 254 * be numerically more robust than the {@link #PolygonsSet(Collection) general 255 * constructor} using {@link SubHyperplane subhyperplanes}.</p> 256 * @param hyperplaneThickness tolerance below which points are consider to 257 * belong to the hyperplane (which is therefore more a slab) 258 * @param vertices vertices of the simple loop boundary 259 * @return the BSP tree of the input vertices 260 */ verticesToTree(final double hyperplaneThickness, final Vector2D ... vertices)261 private static BSPTree<Euclidean2D> verticesToTree(final double hyperplaneThickness, 262 final Vector2D ... vertices) { 263 264 final int n = vertices.length; 265 if (n == 0) { 266 // the tree represents the whole space 267 return new BSPTree<Euclidean2D>(Boolean.TRUE); 268 } 269 270 // build the vertices 271 final Vertex[] vArray = new Vertex[n]; 272 for (int i = 0; i < n; ++i) { 273 vArray[i] = new Vertex(vertices[i]); 274 } 275 276 // build the edges 277 List<Edge> edges = new ArrayList<Edge>(n); 278 for (int i = 0; i < n; ++i) { 279 280 // get the endpoints of the edge 281 final Vertex start = vArray[i]; 282 final Vertex end = vArray[(i + 1) % n]; 283 284 // get the line supporting the edge, taking care not to recreate it 285 // if it was already created earlier due to another edge being aligned 286 // with the current one 287 Line line = start.sharedLineWith(end); 288 if (line == null) { 289 line = new Line(start.getLocation(), end.getLocation(), hyperplaneThickness); 290 } 291 292 // create the edge and store it 293 edges.add(new Edge(start, end, line)); 294 295 // check if another vertex also happens to be on this line 296 for (final Vertex vertex : vArray) { 297 if (vertex != start && vertex != end && 298 FastMath.abs(line.getOffset((Point<Euclidean2D>) vertex.getLocation())) <= hyperplaneThickness) { 299 vertex.bindWith(line); 300 } 301 } 302 303 } 304 305 // build the tree top-down 306 final BSPTree<Euclidean2D> tree = new BSPTree<Euclidean2D>(); 307 insertEdges(hyperplaneThickness, tree, edges); 308 309 return tree; 310 311 } 312 313 /** Recursively build a tree by inserting cut sub-hyperplanes. 314 * @param hyperplaneThickness tolerance below which points are consider to 315 * belong to the hyperplane (which is therefore more a slab) 316 * @param node current tree node (it is a leaf node at the beginning 317 * of the call) 318 * @param edges list of edges to insert in the cell defined by this node 319 * (excluding edges not belonging to the cell defined by this node) 320 */ insertEdges(final double hyperplaneThickness, final BSPTree<Euclidean2D> node, final List<Edge> edges)321 private static void insertEdges(final double hyperplaneThickness, 322 final BSPTree<Euclidean2D> node, 323 final List<Edge> edges) { 324 325 // find an edge with an hyperplane that can be inserted in the node 326 int index = 0; 327 Edge inserted =null; 328 while (inserted == null && index < edges.size()) { 329 inserted = edges.get(index++); 330 if (inserted.getNode() == null) { 331 if (node.insertCut(inserted.getLine())) { 332 inserted.setNode(node); 333 } else { 334 inserted = null; 335 } 336 } else { 337 inserted = null; 338 } 339 } 340 341 if (inserted == null) { 342 // no suitable edge was found, the node remains a leaf node 343 // we need to set its inside/outside boolean indicator 344 final BSPTree<Euclidean2D> parent = node.getParent(); 345 if (parent == null || node == parent.getMinus()) { 346 node.setAttribute(Boolean.TRUE); 347 } else { 348 node.setAttribute(Boolean.FALSE); 349 } 350 return; 351 } 352 353 // we have split the node by inserting an edge as a cut sub-hyperplane 354 // distribute the remaining edges in the two sub-trees 355 final List<Edge> plusList = new ArrayList<Edge>(); 356 final List<Edge> minusList = new ArrayList<Edge>(); 357 for (final Edge edge : edges) { 358 if (edge != inserted) { 359 final double startOffset = inserted.getLine().getOffset((Point<Euclidean2D>) edge.getStart().getLocation()); 360 final double endOffset = inserted.getLine().getOffset((Point<Euclidean2D>) edge.getEnd().getLocation()); 361 Side startSide = (FastMath.abs(startOffset) <= hyperplaneThickness) ? 362 Side.HYPER : ((startOffset < 0) ? Side.MINUS : Side.PLUS); 363 Side endSide = (FastMath.abs(endOffset) <= hyperplaneThickness) ? 364 Side.HYPER : ((endOffset < 0) ? Side.MINUS : Side.PLUS); 365 switch (startSide) { 366 case PLUS: 367 if (endSide == Side.MINUS) { 368 // we need to insert a split point on the hyperplane 369 final Vertex splitPoint = edge.split(inserted.getLine()); 370 minusList.add(splitPoint.getOutgoing()); 371 plusList.add(splitPoint.getIncoming()); 372 } else { 373 plusList.add(edge); 374 } 375 break; 376 case MINUS: 377 if (endSide == Side.PLUS) { 378 // we need to insert a split point on the hyperplane 379 final Vertex splitPoint = edge.split(inserted.getLine()); 380 minusList.add(splitPoint.getIncoming()); 381 plusList.add(splitPoint.getOutgoing()); 382 } else { 383 minusList.add(edge); 384 } 385 break; 386 default: 387 if (endSide == Side.PLUS) { 388 plusList.add(edge); 389 } else if (endSide == Side.MINUS) { 390 minusList.add(edge); 391 } 392 break; 393 } 394 } 395 } 396 397 // recurse through lower levels 398 if (!plusList.isEmpty()) { 399 insertEdges(hyperplaneThickness, node.getPlus(), plusList); 400 } else { 401 node.getPlus().setAttribute(Boolean.FALSE); 402 } 403 if (!minusList.isEmpty()) { 404 insertEdges(hyperplaneThickness, node.getMinus(), minusList); 405 } else { 406 node.getMinus().setAttribute(Boolean.TRUE); 407 } 408 409 } 410 411 /** Internal class for holding vertices while they are processed to build a BSP tree. */ 412 private static class Vertex { 413 414 /** Vertex location. */ 415 private final Vector2D location; 416 417 /** Incoming edge. */ 418 private Edge incoming; 419 420 /** Outgoing edge. */ 421 private Edge outgoing; 422 423 /** Lines bound with this vertex. */ 424 private final List<Line> lines; 425 426 /** Build a non-processed vertex not owned by any node yet. 427 * @param location vertex location 428 */ Vertex(final Vector2D location)429 Vertex(final Vector2D location) { 430 this.location = location; 431 this.incoming = null; 432 this.outgoing = null; 433 this.lines = new ArrayList<Line>(); 434 } 435 436 /** Get Vertex location. 437 * @return vertex location 438 */ getLocation()439 public Vector2D getLocation() { 440 return location; 441 } 442 443 /** Bind a line considered to contain this vertex. 444 * @param line line to bind with this vertex 445 */ bindWith(final Line line)446 public void bindWith(final Line line) { 447 lines.add(line); 448 } 449 450 /** Get the common line bound with both the instance and another vertex, if any. 451 * <p> 452 * When two vertices are both bound to the same line, this means they are 453 * already handled by node associated with this line, so there is no need 454 * to create a cut hyperplane for them. 455 * </p> 456 * @param vertex other vertex to check instance against 457 * @return line bound with both the instance and another vertex, or null if the 458 * two vertices do not share a line yet 459 */ sharedLineWith(final Vertex vertex)460 public Line sharedLineWith(final Vertex vertex) { 461 for (final Line line1 : lines) { 462 for (final Line line2 : vertex.lines) { 463 if (line1 == line2) { 464 return line1; 465 } 466 } 467 } 468 return null; 469 } 470 471 /** Set incoming edge. 472 * <p> 473 * The line supporting the incoming edge is automatically bound 474 * with the instance. 475 * </p> 476 * @param incoming incoming edge 477 */ setIncoming(final Edge incoming)478 public void setIncoming(final Edge incoming) { 479 this.incoming = incoming; 480 bindWith(incoming.getLine()); 481 } 482 483 /** Get incoming edge. 484 * @return incoming edge 485 */ getIncoming()486 public Edge getIncoming() { 487 return incoming; 488 } 489 490 /** Set outgoing edge. 491 * <p> 492 * The line supporting the outgoing edge is automatically bound 493 * with the instance. 494 * </p> 495 * @param outgoing outgoing edge 496 */ setOutgoing(final Edge outgoing)497 public void setOutgoing(final Edge outgoing) { 498 this.outgoing = outgoing; 499 bindWith(outgoing.getLine()); 500 } 501 502 /** Get outgoing edge. 503 * @return outgoing edge 504 */ getOutgoing()505 public Edge getOutgoing() { 506 return outgoing; 507 } 508 509 } 510 511 /** Internal class for holding edges while they are processed to build a BSP tree. */ 512 private static class Edge { 513 514 /** Start vertex. */ 515 private final Vertex start; 516 517 /** End vertex. */ 518 private final Vertex end; 519 520 /** Line supporting the edge. */ 521 private final Line line; 522 523 /** Node whose cut hyperplane contains this edge. */ 524 private BSPTree<Euclidean2D> node; 525 526 /** Build an edge not contained in any node yet. 527 * @param start start vertex 528 * @param end end vertex 529 * @param line line supporting the edge 530 */ Edge(final Vertex start, final Vertex end, final Line line)531 Edge(final Vertex start, final Vertex end, final Line line) { 532 533 this.start = start; 534 this.end = end; 535 this.line = line; 536 this.node = null; 537 538 // connect the vertices back to the edge 539 start.setOutgoing(this); 540 end.setIncoming(this); 541 542 } 543 544 /** Get start vertex. 545 * @return start vertex 546 */ getStart()547 public Vertex getStart() { 548 return start; 549 } 550 551 /** Get end vertex. 552 * @return end vertex 553 */ getEnd()554 public Vertex getEnd() { 555 return end; 556 } 557 558 /** Get the line supporting this edge. 559 * @return line supporting this edge 560 */ getLine()561 public Line getLine() { 562 return line; 563 } 564 565 /** Set the node whose cut hyperplane contains this edge. 566 * @param node node whose cut hyperplane contains this edge 567 */ setNode(final BSPTree<Euclidean2D> node)568 public void setNode(final BSPTree<Euclidean2D> node) { 569 this.node = node; 570 } 571 572 /** Get the node whose cut hyperplane contains this edge. 573 * @return node whose cut hyperplane contains this edge 574 * (null if edge has not yet been inserted into the BSP tree) 575 */ getNode()576 public BSPTree<Euclidean2D> getNode() { 577 return node; 578 } 579 580 /** Split the edge. 581 * <p> 582 * Once split, this edge is not referenced anymore by the vertices, 583 * it is replaced by the two half-edges and an intermediate splitting 584 * vertex is introduced to connect these two halves. 585 * </p> 586 * @param splitLine line splitting the edge in two halves 587 * @return split vertex (its incoming and outgoing edges are the two halves) 588 */ split(final Line splitLine)589 public Vertex split(final Line splitLine) { 590 final Vertex splitVertex = new Vertex(line.intersection(splitLine)); 591 splitVertex.bindWith(splitLine); 592 final Edge startHalf = new Edge(start, splitVertex, line); 593 final Edge endHalf = new Edge(splitVertex, end, line); 594 startHalf.node = node; 595 endHalf.node = node; 596 return splitVertex; 597 } 598 599 } 600 601 /** {@inheritDoc} */ 602 @Override buildNew(final BSPTree<Euclidean2D> tree)603 public PolygonsSet buildNew(final BSPTree<Euclidean2D> tree) { 604 return new PolygonsSet(tree, getTolerance()); 605 } 606 607 /** {@inheritDoc} */ 608 @Override computeGeometricalProperties()609 protected void computeGeometricalProperties() { 610 611 final Vector2D[][] v = getVertices(); 612 613 if (v.length == 0) { 614 final BSPTree<Euclidean2D> tree = getTree(false); 615 if (tree.getCut() == null && (Boolean) tree.getAttribute()) { 616 // the instance covers the whole space 617 setSize(Double.POSITIVE_INFINITY); 618 setBarycenter((Point<Euclidean2D>) Vector2D.NaN); 619 } else { 620 setSize(0); 621 setBarycenter((Point<Euclidean2D>) new Vector2D(0, 0)); 622 } 623 } else if (v[0][0] == null) { 624 // there is at least one open-loop: the polygon is infinite 625 setSize(Double.POSITIVE_INFINITY); 626 setBarycenter((Point<Euclidean2D>) Vector2D.NaN); 627 } else { 628 // all loops are closed, we compute some integrals around the shape 629 630 double sum = 0; 631 double sumX = 0; 632 double sumY = 0; 633 634 for (Vector2D[] loop : v) { 635 double x1 = loop[loop.length - 1].getX(); 636 double y1 = loop[loop.length - 1].getY(); 637 for (final Vector2D point : loop) { 638 final double x0 = x1; 639 final double y0 = y1; 640 x1 = point.getX(); 641 y1 = point.getY(); 642 final double factor = x0 * y1 - y0 * x1; 643 sum += factor; 644 sumX += factor * (x0 + x1); 645 sumY += factor * (y0 + y1); 646 } 647 } 648 649 if (sum < 0) { 650 // the polygon as a finite outside surrounded by an infinite inside 651 setSize(Double.POSITIVE_INFINITY); 652 setBarycenter((Point<Euclidean2D>) Vector2D.NaN); 653 } else { 654 setSize(sum / 2); 655 setBarycenter((Point<Euclidean2D>) new Vector2D(sumX / (3 * sum), sumY / (3 * sum))); 656 } 657 658 } 659 660 } 661 662 /** Get the vertices of the polygon. 663 * <p>The polygon boundary can be represented as an array of loops, 664 * each loop being itself an array of vertices.</p> 665 * <p>In order to identify open loops which start and end by 666 * infinite edges, the open loops arrays start with a null point. In 667 * this case, the first non null point and the last point of the 668 * array do not represent real vertices, they are dummy points 669 * intended only to get the direction of the first and last edge. An 670 * open loop consisting of a single infinite line will therefore be 671 * represented by a three elements array with one null point 672 * followed by two dummy points. The open loops are always the first 673 * ones in the loops array.</p> 674 * <p>If the polygon has no boundary at all, a zero length loop 675 * array will be returned.</p> 676 * <p>All line segments in the various loops have the inside of the 677 * region on their left side and the outside on their right side 678 * when moving in the underlying line direction. This means that 679 * closed loops surrounding finite areas obey the direct 680 * trigonometric orientation.</p> 681 * @return vertices of the polygon, organized as oriented boundary 682 * loops with the open loops first (the returned value is guaranteed 683 * to be non-null) 684 */ getVertices()685 public Vector2D[][] getVertices() { 686 if (vertices == null) { 687 if (getTree(false).getCut() == null) { 688 vertices = new Vector2D[0][]; 689 } else { 690 691 // build the unconnected segments 692 final SegmentsBuilder visitor = new SegmentsBuilder(getTolerance()); 693 getTree(true).visit(visitor); 694 final List<ConnectableSegment> segments = visitor.getSegments(); 695 696 // connect all segments, using topological criteria first 697 // and using Euclidean distance only as a last resort 698 int pending = segments.size(); 699 pending -= naturalFollowerConnections(segments); 700 if (pending > 0) { 701 pending -= splitEdgeConnections(segments); 702 } 703 if (pending > 0) { 704 pending -= closeVerticesConnections(segments); 705 } 706 707 // create the segment loops 708 final ArrayList<List<Segment>> loops = new ArrayList<List<Segment>>(); 709 for (ConnectableSegment s = getUnprocessed(segments); s != null; s = getUnprocessed(segments)) { 710 final List<Segment> loop = followLoop(s); 711 if (loop != null) { 712 if (loop.get(0).getStart() == null) { 713 // this is an open loop, we put it on the front 714 loops.add(0, loop); 715 } else { 716 // this is a closed loop, we put it on the back 717 loops.add(loop); 718 } 719 } 720 } 721 722 // transform the loops in an array of arrays of points 723 vertices = new Vector2D[loops.size()][]; 724 int i = 0; 725 726 for (final List<Segment> loop : loops) { 727 if (loop.size() < 2 || 728 (loop.size() == 2 && loop.get(0).getStart() == null && loop.get(1).getEnd() == null)) { 729 // single infinite line 730 final Line line = loop.get(0).getLine(); 731 vertices[i++] = new Vector2D[] { 732 null, 733 line.toSpace((Point<Euclidean1D>) new Vector1D(-Float.MAX_VALUE)), 734 line.toSpace((Point<Euclidean1D>) new Vector1D(+Float.MAX_VALUE)) 735 }; 736 } else if (loop.get(0).getStart() == null) { 737 // open loop with at least one real point 738 final Vector2D[] array = new Vector2D[loop.size() + 2]; 739 int j = 0; 740 for (Segment segment : loop) { 741 742 if (j == 0) { 743 // null point and first dummy point 744 double x = segment.getLine().toSubSpace((Point<Euclidean2D>) segment.getEnd()).getX(); 745 x -= FastMath.max(1.0, FastMath.abs(x / 2)); 746 array[j++] = null; 747 array[j++] = segment.getLine().toSpace((Point<Euclidean1D>) new Vector1D(x)); 748 } 749 750 if (j < (array.length - 1)) { 751 // current point 752 array[j++] = segment.getEnd(); 753 } 754 755 if (j == (array.length - 1)) { 756 // last dummy point 757 double x = segment.getLine().toSubSpace((Point<Euclidean2D>) segment.getStart()).getX(); 758 x += FastMath.max(1.0, FastMath.abs(x / 2)); 759 array[j++] = segment.getLine().toSpace((Point<Euclidean1D>) new Vector1D(x)); 760 } 761 762 } 763 vertices[i++] = array; 764 } else { 765 final Vector2D[] array = new Vector2D[loop.size()]; 766 int j = 0; 767 for (Segment segment : loop) { 768 array[j++] = segment.getStart(); 769 } 770 vertices[i++] = array; 771 } 772 } 773 774 } 775 } 776 777 return vertices.clone(); 778 779 } 780 781 /** Connect the segments using only natural follower information. 782 * @param segments segments complete segments list 783 * @return number of connections performed 784 */ naturalFollowerConnections(final List<ConnectableSegment> segments)785 private int naturalFollowerConnections(final List<ConnectableSegment> segments) { 786 int connected = 0; 787 for (final ConnectableSegment segment : segments) { 788 if (segment.getNext() == null) { 789 final BSPTree<Euclidean2D> node = segment.getNode(); 790 final BSPTree<Euclidean2D> end = segment.getEndNode(); 791 for (final ConnectableSegment candidateNext : segments) { 792 if (candidateNext.getPrevious() == null && 793 candidateNext.getNode() == end && 794 candidateNext.getStartNode() == node) { 795 // connect the two segments 796 segment.setNext(candidateNext); 797 candidateNext.setPrevious(segment); 798 ++connected; 799 break; 800 } 801 } 802 } 803 } 804 return connected; 805 } 806 807 /** Connect the segments resulting from a line splitting a straight edge. 808 * @param segments segments complete segments list 809 * @return number of connections performed 810 */ splitEdgeConnections(final List<ConnectableSegment> segments)811 private int splitEdgeConnections(final List<ConnectableSegment> segments) { 812 int connected = 0; 813 for (final ConnectableSegment segment : segments) { 814 if (segment.getNext() == null) { 815 final Hyperplane<Euclidean2D> hyperplane = segment.getNode().getCut().getHyperplane(); 816 final BSPTree<Euclidean2D> end = segment.getEndNode(); 817 for (final ConnectableSegment candidateNext : segments) { 818 if (candidateNext.getPrevious() == null && 819 candidateNext.getNode().getCut().getHyperplane() == hyperplane && 820 candidateNext.getStartNode() == end) { 821 // connect the two segments 822 segment.setNext(candidateNext); 823 candidateNext.setPrevious(segment); 824 ++connected; 825 break; 826 } 827 } 828 } 829 } 830 return connected; 831 } 832 833 /** Connect the segments using Euclidean distance. 834 * <p> 835 * This connection heuristic should be used last, as it relies 836 * only on a fuzzy distance criterion. 837 * </p> 838 * @param segments segments complete segments list 839 * @return number of connections performed 840 */ closeVerticesConnections(final List<ConnectableSegment> segments)841 private int closeVerticesConnections(final List<ConnectableSegment> segments) { 842 int connected = 0; 843 for (final ConnectableSegment segment : segments) { 844 if (segment.getNext() == null && segment.getEnd() != null) { 845 final Vector2D end = segment.getEnd(); 846 ConnectableSegment selectedNext = null; 847 double min = Double.POSITIVE_INFINITY; 848 for (final ConnectableSegment candidateNext : segments) { 849 if (candidateNext.getPrevious() == null && candidateNext.getStart() != null) { 850 final double distance = Vector2D.distance(end, candidateNext.getStart()); 851 if (distance < min) { 852 selectedNext = candidateNext; 853 min = distance; 854 } 855 } 856 } 857 if (min <= getTolerance()) { 858 // connect the two segments 859 segment.setNext(selectedNext); 860 selectedNext.setPrevious(segment); 861 ++connected; 862 } 863 } 864 } 865 return connected; 866 } 867 868 /** Get first unprocessed segment from a list. 869 * @param segments segments list 870 * @return first segment that has not been processed yet 871 * or null if all segments have been processed 872 */ getUnprocessed(final List<ConnectableSegment> segments)873 private ConnectableSegment getUnprocessed(final List<ConnectableSegment> segments) { 874 for (final ConnectableSegment segment : segments) { 875 if (!segment.isProcessed()) { 876 return segment; 877 } 878 } 879 return null; 880 } 881 882 /** Build the loop containing a segment. 883 * <p> 884 * The segment put in the loop will be marked as processed. 885 * </p> 886 * @param defining segment used to define the loop 887 * @return loop containing the segment (may be null if the loop is a 888 * degenerated infinitely thin 2 points loop 889 */ followLoop(final ConnectableSegment defining)890 private List<Segment> followLoop(final ConnectableSegment defining) { 891 892 final List<Segment> loop = new ArrayList<Segment>(); 893 loop.add(defining); 894 defining.setProcessed(true); 895 896 // add segments in connection order 897 ConnectableSegment next = defining.getNext(); 898 while (next != defining && next != null) { 899 loop.add(next); 900 next.setProcessed(true); 901 next = next.getNext(); 902 } 903 904 if (next == null) { 905 // the loop is open and we have found its end, 906 // we need to find its start too 907 ConnectableSegment previous = defining.getPrevious(); 908 while (previous != null) { 909 loop.add(0, previous); 910 previous.setProcessed(true); 911 previous = previous.getPrevious(); 912 } 913 } 914 915 // filter out spurious vertices 916 filterSpuriousVertices(loop); 917 918 if (loop.size() == 2 && loop.get(0).getStart() != null) { 919 // this is a degenerated infinitely thin closed loop, we simply ignore it 920 return null; 921 } else { 922 return loop; 923 } 924 925 } 926 927 /** Filter out spurious vertices on straight lines (at machine precision). 928 * @param loop segments loop to filter (will be modified in-place) 929 */ filterSpuriousVertices(final List<Segment> loop)930 private void filterSpuriousVertices(final List<Segment> loop) { 931 for (int i = 0; i < loop.size(); ++i) { 932 final Segment previous = loop.get(i); 933 int j = (i + 1) % loop.size(); 934 final Segment next = loop.get(j); 935 if (next != null && 936 Precision.equals(previous.getLine().getAngle(), next.getLine().getAngle(), Precision.EPSILON)) { 937 // the vertex between the two edges is a spurious one 938 // replace the two segments by a single one 939 loop.set(j, new Segment(previous.getStart(), next.getEnd(), previous.getLine())); 940 loop.remove(i--); 941 } 942 } 943 } 944 945 /** Private extension of Segment allowing connection. */ 946 private static class ConnectableSegment extends Segment { 947 948 /** Node containing segment. */ 949 private final BSPTree<Euclidean2D> node; 950 951 /** Node whose intersection with current node defines start point. */ 952 private final BSPTree<Euclidean2D> startNode; 953 954 /** Node whose intersection with current node defines end point. */ 955 private final BSPTree<Euclidean2D> endNode; 956 957 /** Previous segment. */ 958 private ConnectableSegment previous; 959 960 /** Next segment. */ 961 private ConnectableSegment next; 962 963 /** Indicator for completely processed segments. */ 964 private boolean processed; 965 966 /** Build a segment. 967 * @param start start point of the segment 968 * @param end end point of the segment 969 * @param line line containing the segment 970 * @param node node containing the segment 971 * @param startNode node whose intersection with current node defines start point 972 * @param endNode node whose intersection with current node defines end point 973 */ ConnectableSegment(final Vector2D start, final Vector2D end, final Line line, final BSPTree<Euclidean2D> node, final BSPTree<Euclidean2D> startNode, final BSPTree<Euclidean2D> endNode)974 ConnectableSegment(final Vector2D start, final Vector2D end, final Line line, 975 final BSPTree<Euclidean2D> node, 976 final BSPTree<Euclidean2D> startNode, 977 final BSPTree<Euclidean2D> endNode) { 978 super(start, end, line); 979 this.node = node; 980 this.startNode = startNode; 981 this.endNode = endNode; 982 this.previous = null; 983 this.next = null; 984 this.processed = false; 985 } 986 987 /** Get the node containing segment. 988 * @return node containing segment 989 */ getNode()990 public BSPTree<Euclidean2D> getNode() { 991 return node; 992 } 993 994 /** Get the node whose intersection with current node defines start point. 995 * @return node whose intersection with current node defines start point 996 */ getStartNode()997 public BSPTree<Euclidean2D> getStartNode() { 998 return startNode; 999 } 1000 1001 /** Get the node whose intersection with current node defines end point. 1002 * @return node whose intersection with current node defines end point 1003 */ getEndNode()1004 public BSPTree<Euclidean2D> getEndNode() { 1005 return endNode; 1006 } 1007 1008 /** Get the previous segment. 1009 * @return previous segment 1010 */ getPrevious()1011 public ConnectableSegment getPrevious() { 1012 return previous; 1013 } 1014 1015 /** Set the previous segment. 1016 * @param previous previous segment 1017 */ setPrevious(final ConnectableSegment previous)1018 public void setPrevious(final ConnectableSegment previous) { 1019 this.previous = previous; 1020 } 1021 1022 /** Get the next segment. 1023 * @return next segment 1024 */ getNext()1025 public ConnectableSegment getNext() { 1026 return next; 1027 } 1028 1029 /** Set the next segment. 1030 * @param next previous segment 1031 */ setNext(final ConnectableSegment next)1032 public void setNext(final ConnectableSegment next) { 1033 this.next = next; 1034 } 1035 1036 /** Set the processed flag. 1037 * @param processed processed flag to set 1038 */ setProcessed(final boolean processed)1039 public void setProcessed(final boolean processed) { 1040 this.processed = processed; 1041 } 1042 1043 /** Check if the segment has been processed. 1044 * @return true if the segment has been processed 1045 */ isProcessed()1046 public boolean isProcessed() { 1047 return processed; 1048 } 1049 1050 } 1051 1052 /** Visitor building segments. */ 1053 private static class SegmentsBuilder implements BSPTreeVisitor<Euclidean2D> { 1054 1055 /** Tolerance for close nodes connection. */ 1056 private final double tolerance; 1057 1058 /** Built segments. */ 1059 private final List<ConnectableSegment> segments; 1060 1061 /** Simple constructor. 1062 * @param tolerance tolerance for close nodes connection 1063 */ SegmentsBuilder(final double tolerance)1064 SegmentsBuilder(final double tolerance) { 1065 this.tolerance = tolerance; 1066 this.segments = new ArrayList<ConnectableSegment>(); 1067 } 1068 1069 /** {@inheritDoc} */ visitOrder(final BSPTree<Euclidean2D> node)1070 public Order visitOrder(final BSPTree<Euclidean2D> node) { 1071 return Order.MINUS_SUB_PLUS; 1072 } 1073 1074 /** {@inheritDoc} */ visitInternalNode(final BSPTree<Euclidean2D> node)1075 public void visitInternalNode(final BSPTree<Euclidean2D> node) { 1076 @SuppressWarnings("unchecked") 1077 final BoundaryAttribute<Euclidean2D> attribute = (BoundaryAttribute<Euclidean2D>) node.getAttribute(); 1078 final Iterable<BSPTree<Euclidean2D>> splitters = attribute.getSplitters(); 1079 if (attribute.getPlusOutside() != null) { 1080 addContribution(attribute.getPlusOutside(), node, splitters, false); 1081 } 1082 if (attribute.getPlusInside() != null) { 1083 addContribution(attribute.getPlusInside(), node, splitters, true); 1084 } 1085 } 1086 1087 /** {@inheritDoc} */ visitLeafNode(final BSPTree<Euclidean2D> node)1088 public void visitLeafNode(final BSPTree<Euclidean2D> node) { 1089 } 1090 1091 /** Add the contribution of a boundary facet. 1092 * @param sub boundary facet 1093 * @param node node containing segment 1094 * @param splitters splitters for the boundary facet 1095 * @param reversed if true, the facet has the inside on its plus side 1096 */ addContribution(final SubHyperplane<Euclidean2D> sub, final BSPTree<Euclidean2D> node, final Iterable<BSPTree<Euclidean2D>> splitters, final boolean reversed)1097 private void addContribution(final SubHyperplane<Euclidean2D> sub, 1098 final BSPTree<Euclidean2D> node, 1099 final Iterable<BSPTree<Euclidean2D>> splitters, 1100 final boolean reversed) { 1101 @SuppressWarnings("unchecked") 1102 final AbstractSubHyperplane<Euclidean2D, Euclidean1D> absSub = 1103 (AbstractSubHyperplane<Euclidean2D, Euclidean1D>) sub; 1104 final Line line = (Line) sub.getHyperplane(); 1105 final List<Interval> intervals = ((IntervalsSet) absSub.getRemainingRegion()).asList(); 1106 for (final Interval i : intervals) { 1107 1108 // find the 2D points 1109 final Vector2D startV = Double.isInfinite(i.getInf()) ? 1110 null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getInf())); 1111 final Vector2D endV = Double.isInfinite(i.getSup()) ? 1112 null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getSup())); 1113 1114 // recover the connectivity information 1115 final BSPTree<Euclidean2D> startN = selectClosest(startV, splitters); 1116 final BSPTree<Euclidean2D> endN = selectClosest(endV, splitters); 1117 1118 if (reversed) { 1119 segments.add(new ConnectableSegment(endV, startV, line.getReverse(), 1120 node, endN, startN)); 1121 } else { 1122 segments.add(new ConnectableSegment(startV, endV, line, 1123 node, startN, endN)); 1124 } 1125 1126 } 1127 } 1128 1129 /** Select the node whose cut sub-hyperplane is closest to specified point. 1130 * @param point reference point 1131 * @param candidates candidate nodes 1132 * @return node closest to point, or null if no node is closer than tolerance 1133 */ selectClosest(final Vector2D point, final Iterable<BSPTree<Euclidean2D>> candidates)1134 private BSPTree<Euclidean2D> selectClosest(final Vector2D point, final Iterable<BSPTree<Euclidean2D>> candidates) { 1135 1136 BSPTree<Euclidean2D> selected = null; 1137 double min = Double.POSITIVE_INFINITY; 1138 1139 for (final BSPTree<Euclidean2D> node : candidates) { 1140 final double distance = FastMath.abs(node.getCut().getHyperplane().getOffset(point)); 1141 if (distance < min) { 1142 selected = node; 1143 min = distance; 1144 } 1145 } 1146 1147 return min <= tolerance ? selected : null; 1148 1149 } 1150 1151 /** Get the segments. 1152 * @return built segments 1153 */ getSegments()1154 public List<ConnectableSegment> getSegments() { 1155 return segments; 1156 } 1157 1158 } 1159 1160 } 1161