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.spherical.twod; 18 19 import java.util.ArrayList; 20 import java.util.Collection; 21 import java.util.Collections; 22 import java.util.Iterator; 23 import java.util.List; 24 25 import org.apache.commons.math3.exception.MathIllegalStateException; 26 import org.apache.commons.math3.geometry.enclosing.EnclosingBall; 27 import org.apache.commons.math3.geometry.enclosing.WelzlEncloser; 28 import org.apache.commons.math3.geometry.euclidean.threed.Euclidean3D; 29 import org.apache.commons.math3.geometry.euclidean.threed.Rotation; 30 import org.apache.commons.math3.geometry.euclidean.threed.RotationConvention; 31 import org.apache.commons.math3.geometry.euclidean.threed.SphereGenerator; 32 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; 33 import org.apache.commons.math3.geometry.partitioning.AbstractRegion; 34 import org.apache.commons.math3.geometry.partitioning.BSPTree; 35 import org.apache.commons.math3.geometry.partitioning.BoundaryProjection; 36 import org.apache.commons.math3.geometry.partitioning.RegionFactory; 37 import org.apache.commons.math3.geometry.partitioning.SubHyperplane; 38 import org.apache.commons.math3.geometry.spherical.oned.Sphere1D; 39 import org.apache.commons.math3.util.FastMath; 40 import org.apache.commons.math3.util.MathUtils; 41 42 /** This class represents a region on the 2-sphere: a set of spherical polygons. 43 * @since 3.3 44 */ 45 public class SphericalPolygonsSet extends AbstractRegion<Sphere2D, Sphere1D> { 46 47 /** Boundary defined as an array of closed loops start vertices. */ 48 private List<Vertex> loops; 49 50 /** Build a polygons set representing the whole real 2-sphere. 51 * @param tolerance below which points are consider to be identical 52 */ SphericalPolygonsSet(final double tolerance)53 public SphericalPolygonsSet(final double tolerance) { 54 super(tolerance); 55 } 56 57 /** Build a polygons set representing a hemisphere. 58 * @param pole pole of the hemisphere (the pole is in the inside half) 59 * @param tolerance below which points are consider to be identical 60 */ SphericalPolygonsSet(final Vector3D pole, final double tolerance)61 public SphericalPolygonsSet(final Vector3D pole, final double tolerance) { 62 super(new BSPTree<Sphere2D>(new Circle(pole, tolerance).wholeHyperplane(), 63 new BSPTree<Sphere2D>(Boolean.FALSE), 64 new BSPTree<Sphere2D>(Boolean.TRUE), 65 null), 66 tolerance); 67 } 68 69 /** Build a polygons set representing a regular polygon. 70 * @param center center of the polygon (the center is in the inside half) 71 * @param meridian point defining the reference meridian for first polygon vertex 72 * @param outsideRadius distance of the vertices to the center 73 * @param n number of sides of the polygon 74 * @param tolerance below which points are consider to be identical 75 */ SphericalPolygonsSet(final Vector3D center, final Vector3D meridian, final double outsideRadius, final int n, final double tolerance)76 public SphericalPolygonsSet(final Vector3D center, final Vector3D meridian, 77 final double outsideRadius, final int n, 78 final double tolerance) { 79 this(tolerance, createRegularPolygonVertices(center, meridian, outsideRadius, n)); 80 } 81 82 /** Build a polygons set from a BSP tree. 83 * <p>The leaf nodes of the BSP tree <em>must</em> have a 84 * {@code Boolean} attribute representing the inside status of 85 * the corresponding cell (true for inside cells, false for outside 86 * cells). In order to avoid building too many small objects, it is 87 * recommended to use the predefined constants 88 * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p> 89 * @param tree inside/outside BSP tree representing the region 90 * @param tolerance below which points are consider to be identical 91 */ SphericalPolygonsSet(final BSPTree<Sphere2D> tree, final double tolerance)92 public SphericalPolygonsSet(final BSPTree<Sphere2D> tree, final double tolerance) { 93 super(tree, tolerance); 94 } 95 96 /** Build a polygons set from a Boundary REPresentation (B-rep). 97 * <p>The boundary is provided as a collection of {@link 98 * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the 99 * interior part of the region on its minus side and the exterior on 100 * its plus side.</p> 101 * <p>The boundary elements can be in any order, and can form 102 * several non-connected sets (like for example polygons with holes 103 * or a set of disjoint polygons considered as a whole). In 104 * fact, the elements do not even need to be connected together 105 * (their topological connections are not used here). However, if the 106 * boundary does not really separate an inside open from an outside 107 * open (open having here its topological meaning), then subsequent 108 * calls to the {@link 109 * org.apache.commons.math3.geometry.partitioning.Region#checkPoint(org.apache.commons.math3.geometry.Point) 110 * checkPoint} method will not be meaningful anymore.</p> 111 * <p>If the boundary is empty, the region will represent the whole 112 * space.</p> 113 * @param boundary collection of boundary elements, as a 114 * collection of {@link SubHyperplane SubHyperplane} objects 115 * @param tolerance below which points are consider to be identical 116 */ SphericalPolygonsSet(final Collection<SubHyperplane<Sphere2D>> boundary, final double tolerance)117 public SphericalPolygonsSet(final Collection<SubHyperplane<Sphere2D>> boundary, final double tolerance) { 118 super(boundary, tolerance); 119 } 120 121 /** Build a polygon from a simple list of vertices. 122 * <p>The boundary is provided as a list of points considering to 123 * represent the vertices of a simple loop. The interior part of the 124 * region is on the left side of this path and the exterior is on its 125 * right side.</p> 126 * <p>This constructor does not handle polygons with a boundary 127 * forming several disconnected paths (such as polygons with holes).</p> 128 * <p>For cases where this simple constructor applies, it is expected to 129 * be numerically more robust than the {@link #SphericalPolygonsSet(Collection, 130 * double) general constructor} using {@link SubHyperplane subhyperplanes}.</p> 131 * <p>If the list is empty, the region will represent the whole 132 * space.</p> 133 * <p> 134 * Polygons with thin pikes or dents are inherently difficult to handle because 135 * they involve circles with almost opposite directions at some vertices. Polygons 136 * whose vertices come from some physical measurement with noise are also 137 * difficult because an edge that should be straight may be broken in lots of 138 * different pieces with almost equal directions. In both cases, computing the 139 * circles intersections is not numerically robust due to the almost 0 or almost 140 * π angle. Such cases need to carefully adjust the {@code hyperplaneThickness} 141 * parameter. A too small value would often lead to completely wrong polygons 142 * with large area wrongly identified as inside or outside. Large values are 143 * often much safer. As a rule of thumb, a value slightly below the size of the 144 * most accurate detail needed is a good value for the {@code hyperplaneThickness} 145 * parameter. 146 * </p> 147 * @param hyperplaneThickness tolerance below which points are considered to 148 * belong to the hyperplane (which is therefore more a slab) 149 * @param vertices vertices of the simple loop boundary 150 */ SphericalPolygonsSet(final double hyperplaneThickness, final S2Point ... vertices)151 public SphericalPolygonsSet(final double hyperplaneThickness, final S2Point ... vertices) { 152 super(verticesToTree(hyperplaneThickness, vertices), hyperplaneThickness); 153 } 154 155 /** Build the vertices representing a regular polygon. 156 * @param center center of the polygon (the center is in the inside half) 157 * @param meridian point defining the reference meridian for first polygon vertex 158 * @param outsideRadius distance of the vertices to the center 159 * @param n number of sides of the polygon 160 * @return vertices array 161 */ createRegularPolygonVertices(final Vector3D center, final Vector3D meridian, final double outsideRadius, final int n)162 private static S2Point[] createRegularPolygonVertices(final Vector3D center, final Vector3D meridian, 163 final double outsideRadius, final int n) { 164 final S2Point[] array = new S2Point[n]; 165 final Rotation r0 = new Rotation(Vector3D.crossProduct(center, meridian), 166 outsideRadius, RotationConvention.VECTOR_OPERATOR); 167 array[0] = new S2Point(r0.applyTo(center)); 168 169 final Rotation r = new Rotation(center, MathUtils.TWO_PI / n, RotationConvention.VECTOR_OPERATOR); 170 for (int i = 1; i < n; ++i) { 171 array[i] = new S2Point(r.applyTo(array[i - 1].getVector())); 172 } 173 174 return array; 175 } 176 177 /** Build the BSP tree of a polygons set from a simple list of vertices. 178 * <p>The boundary is provided as a list of points considering to 179 * represent the vertices of a simple loop. The interior part of the 180 * region is on the left side of this path and the exterior is on its 181 * right side.</p> 182 * <p>This constructor does not handle polygons with a boundary 183 * forming several disconnected paths (such as polygons with holes).</p> 184 * <p>This constructor handles only polygons with edges strictly shorter 185 * than \( \pi \). If longer edges are needed, they need to be broken up 186 * in smaller sub-edges so this constraint holds.</p> 187 * <p>For cases where this simple constructor applies, it is expected to 188 * be numerically more robust than the {@link #PolygonsSet(Collection) general 189 * constructor} using {@link SubHyperplane subhyperplanes}.</p> 190 * @param hyperplaneThickness tolerance below which points are consider to 191 * belong to the hyperplane (which is therefore more a slab) 192 * @param vertices vertices of the simple loop boundary 193 * @return the BSP tree of the input vertices 194 */ verticesToTree(final double hyperplaneThickness, final S2Point ... vertices)195 private static BSPTree<Sphere2D> verticesToTree(final double hyperplaneThickness, 196 final S2Point ... vertices) { 197 198 final int n = vertices.length; 199 if (n == 0) { 200 // the tree represents the whole space 201 return new BSPTree<Sphere2D>(Boolean.TRUE); 202 } 203 204 // build the vertices 205 final Vertex[] vArray = new Vertex[n]; 206 for (int i = 0; i < n; ++i) { 207 vArray[i] = new Vertex(vertices[i]); 208 } 209 210 // build the edges 211 List<Edge> edges = new ArrayList<Edge>(n); 212 Vertex end = vArray[n - 1]; 213 for (int i = 0; i < n; ++i) { 214 215 // get the endpoints of the edge 216 final Vertex start = end; 217 end = vArray[i]; 218 219 // get the circle supporting the edge, taking care not to recreate it 220 // if it was already created earlier due to another edge being aligned 221 // with the current one 222 Circle circle = start.sharedCircleWith(end); 223 if (circle == null) { 224 circle = new Circle(start.getLocation(), end.getLocation(), hyperplaneThickness); 225 } 226 227 // create the edge and store it 228 edges.add(new Edge(start, end, 229 Vector3D.angle(start.getLocation().getVector(), 230 end.getLocation().getVector()), 231 circle)); 232 233 // check if another vertex also happens to be on this circle 234 for (final Vertex vertex : vArray) { 235 if (vertex != start && vertex != end && 236 FastMath.abs(circle.getOffset(vertex.getLocation())) <= hyperplaneThickness) { 237 vertex.bindWith(circle); 238 } 239 } 240 241 } 242 243 // build the tree top-down 244 final BSPTree<Sphere2D> tree = new BSPTree<Sphere2D>(); 245 insertEdges(hyperplaneThickness, tree, edges); 246 247 return tree; 248 249 } 250 251 /** Recursively build a tree by inserting cut sub-hyperplanes. 252 * @param hyperplaneThickness tolerance below which points are considered to 253 * belong to the hyperplane (which is therefore more a slab) 254 * @param node current tree node (it is a leaf node at the beginning 255 * of the call) 256 * @param edges list of edges to insert in the cell defined by this node 257 * (excluding edges not belonging to the cell defined by this node) 258 */ insertEdges(final double hyperplaneThickness, final BSPTree<Sphere2D> node, final List<Edge> edges)259 private static void insertEdges(final double hyperplaneThickness, 260 final BSPTree<Sphere2D> node, 261 final List<Edge> edges) { 262 263 // find an edge with an hyperplane that can be inserted in the node 264 int index = 0; 265 Edge inserted = null; 266 while (inserted == null && index < edges.size()) { 267 inserted = edges.get(index++); 268 if (!node.insertCut(inserted.getCircle())) { 269 inserted = null; 270 } 271 } 272 273 if (inserted == null) { 274 // no suitable edge was found, the node remains a leaf node 275 // we need to set its inside/outside boolean indicator 276 final BSPTree<Sphere2D> parent = node.getParent(); 277 if (parent == null || node == parent.getMinus()) { 278 node.setAttribute(Boolean.TRUE); 279 } else { 280 node.setAttribute(Boolean.FALSE); 281 } 282 return; 283 } 284 285 // we have split the node by inserting an edge as a cut sub-hyperplane 286 // distribute the remaining edges in the two sub-trees 287 final List<Edge> outsideList = new ArrayList<Edge>(); 288 final List<Edge> insideList = new ArrayList<Edge>(); 289 for (final Edge edge : edges) { 290 if (edge != inserted) { 291 edge.split(inserted.getCircle(), outsideList, insideList); 292 } 293 } 294 295 // recurse through lower levels 296 if (!outsideList.isEmpty()) { 297 insertEdges(hyperplaneThickness, node.getPlus(), outsideList); 298 } else { 299 node.getPlus().setAttribute(Boolean.FALSE); 300 } 301 if (!insideList.isEmpty()) { 302 insertEdges(hyperplaneThickness, node.getMinus(), insideList); 303 } else { 304 node.getMinus().setAttribute(Boolean.TRUE); 305 } 306 307 } 308 309 /** {@inheritDoc} */ 310 @Override buildNew(final BSPTree<Sphere2D> tree)311 public SphericalPolygonsSet buildNew(final BSPTree<Sphere2D> tree) { 312 return new SphericalPolygonsSet(tree, getTolerance()); 313 } 314 315 /** {@inheritDoc} 316 * @exception MathIllegalStateException if the tolerance setting does not allow to build 317 * a clean non-ambiguous boundary 318 */ 319 @Override computeGeometricalProperties()320 protected void computeGeometricalProperties() throws MathIllegalStateException { 321 322 final BSPTree<Sphere2D> tree = getTree(true); 323 324 if (tree.getCut() == null) { 325 326 // the instance has a single cell without any boundaries 327 328 if (tree.getCut() == null && (Boolean) tree.getAttribute()) { 329 // the instance covers the whole space 330 setSize(4 * FastMath.PI); 331 setBarycenter(new S2Point(0, 0)); 332 } else { 333 setSize(0); 334 setBarycenter(S2Point.NaN); 335 } 336 337 } else { 338 339 // the instance has a boundary 340 final PropertiesComputer pc = new PropertiesComputer(getTolerance()); 341 tree.visit(pc); 342 setSize(pc.getArea()); 343 setBarycenter(pc.getBarycenter()); 344 345 } 346 347 } 348 349 /** Get the boundary loops of the polygon. 350 * <p>The polygon boundary can be represented as a list of closed loops, 351 * each loop being given by exactly one of its vertices. From each loop 352 * start vertex, one can follow the loop by finding the outgoing edge, 353 * then the end vertex, then the next outgoing edge ... until the start 354 * vertex of the loop (exactly the same instance) is found again once 355 * the full loop has been visited.</p> 356 * <p>If the polygon has no boundary at all, a zero length loop 357 * array will be returned.</p> 358 * <p>If the polygon is a simple one-piece polygon, then the returned 359 * array will contain a single vertex. 360 * </p> 361 * <p>All edges in the various loops have the inside of the region on 362 * their left side (i.e. toward their pole) and the outside on their 363 * right side (i.e. away from their pole) when moving in the underlying 364 * circle direction. This means that the closed loops obey the direct 365 * trigonometric orientation.</p> 366 * @return boundary of the polygon, organized as an unmodifiable list of loops start vertices. 367 * @exception MathIllegalStateException if the tolerance setting does not allow to build 368 * a clean non-ambiguous boundary 369 * @see Vertex 370 * @see Edge 371 */ getBoundaryLoops()372 public List<Vertex> getBoundaryLoops() throws MathIllegalStateException { 373 374 if (loops == null) { 375 if (getTree(false).getCut() == null) { 376 loops = Collections.emptyList(); 377 } else { 378 379 // sort the arcs according to their start point 380 final BSPTree<Sphere2D> root = getTree(true); 381 final EdgesBuilder visitor = new EdgesBuilder(root, getTolerance()); 382 root.visit(visitor); 383 final List<Edge> edges = visitor.getEdges(); 384 385 386 // convert the list of all edges into a list of start vertices 387 loops = new ArrayList<Vertex>(); 388 while (!edges.isEmpty()) { 389 390 // this is an edge belonging to a new loop, store it 391 Edge edge = edges.get(0); 392 final Vertex startVertex = edge.getStart(); 393 loops.add(startVertex); 394 395 // remove all remaining edges in the same loop 396 do { 397 398 // remove one edge 399 for (final Iterator<Edge> iterator = edges.iterator(); iterator.hasNext();) { 400 if (iterator.next() == edge) { 401 iterator.remove(); 402 break; 403 } 404 } 405 406 // go to next edge following the boundary loop 407 edge = edge.getEnd().getOutgoing(); 408 409 } while (edge.getStart() != startVertex); 410 411 } 412 413 } 414 } 415 416 return Collections.unmodifiableList(loops); 417 418 } 419 420 /** Get a spherical cap enclosing the polygon. 421 * <p> 422 * This method is intended as a first test to quickly identify points 423 * that are guaranteed to be outside of the region, hence performing a full 424 * {@link #checkPoint(org.apache.commons.math3.geometry.Vector) checkPoint} 425 * only if the point status remains undecided after the quick check. It is 426 * is therefore mostly useful to speed up computation for small polygons with 427 * complex shapes (say a country boundary on Earth), as the spherical cap will 428 * be small and hence will reliably identify a large part of the sphere as outside, 429 * whereas the full check can be more computing intensive. A typical use case is 430 * therefore: 431 * </p> 432 * <pre> 433 * // compute region, plus an enclosing spherical cap 434 * SphericalPolygonsSet complexShape = ...; 435 * EnclosingBall<Sphere2D, S2Point> cap = complexShape.getEnclosingCap(); 436 * 437 * // check lots of points 438 * for (Vector3D p : points) { 439 * 440 * final Location l; 441 * if (cap.contains(p)) { 442 * // we cannot be sure where the point is 443 * // we need to perform the full computation 444 * l = complexShape.checkPoint(v); 445 * } else { 446 * // no need to do further computation, 447 * // we already know the point is outside 448 * l = Location.OUTSIDE; 449 * } 450 * 451 * // use l ... 452 * 453 * } 454 * </pre> 455 * <p> 456 * In the special cases of empty or whole sphere polygons, special 457 * spherical caps are returned, with angular radius set to negative 458 * or positive infinity so the {@link 459 * EnclosingBall#contains(org.apache.commons.math3.geometry.Point) ball.contains(point)} 460 * method return always false or true. 461 * </p> 462 * <p> 463 * This method is <em>not</em> guaranteed to return the smallest enclosing cap. 464 * </p> 465 * @return a spherical cap enclosing the polygon 466 */ getEnclosingCap()467 public EnclosingBall<Sphere2D, S2Point> getEnclosingCap() { 468 469 // handle special cases first 470 if (isEmpty()) { 471 return new EnclosingBall<Sphere2D, S2Point>(S2Point.PLUS_K, Double.NEGATIVE_INFINITY); 472 } 473 if (isFull()) { 474 return new EnclosingBall<Sphere2D, S2Point>(S2Point.PLUS_K, Double.POSITIVE_INFINITY); 475 } 476 477 // as the polygons is neither empty nor full, it has some boundaries and cut hyperplanes 478 final BSPTree<Sphere2D> root = getTree(false); 479 if (isEmpty(root.getMinus()) && isFull(root.getPlus())) { 480 // the polygon covers an hemisphere, and its boundary is one 2π long edge 481 final Circle circle = (Circle) root.getCut().getHyperplane(); 482 return new EnclosingBall<Sphere2D, S2Point>(new S2Point(circle.getPole()).negate(), 483 0.5 * FastMath.PI); 484 } 485 if (isFull(root.getMinus()) && isEmpty(root.getPlus())) { 486 // the polygon covers an hemisphere, and its boundary is one 2π long edge 487 final Circle circle = (Circle) root.getCut().getHyperplane(); 488 return new EnclosingBall<Sphere2D, S2Point>(new S2Point(circle.getPole()), 489 0.5 * FastMath.PI); 490 } 491 492 // gather some inside points, to be used by the encloser 493 final List<Vector3D> points = getInsidePoints(); 494 495 // extract points from the boundary loops, to be used by the encloser as well 496 final List<Vertex> boundary = getBoundaryLoops(); 497 for (final Vertex loopStart : boundary) { 498 int count = 0; 499 for (Vertex v = loopStart; count == 0 || v != loopStart; v = v.getOutgoing().getEnd()) { 500 ++count; 501 points.add(v.getLocation().getVector()); 502 } 503 } 504 505 // find the smallest enclosing 3D sphere 506 final SphereGenerator generator = new SphereGenerator(); 507 final WelzlEncloser<Euclidean3D, Vector3D> encloser = 508 new WelzlEncloser<Euclidean3D, Vector3D>(getTolerance(), generator); 509 EnclosingBall<Euclidean3D, Vector3D> enclosing3D = encloser.enclose(points); 510 final Vector3D[] support3D = enclosing3D.getSupport(); 511 512 // convert to 3D sphere to spherical cap 513 final double r = enclosing3D.getRadius(); 514 final double h = enclosing3D.getCenter().getNorm(); 515 if (h < getTolerance()) { 516 // the 3D sphere is centered on the unit sphere and covers it 517 // fall back to a crude approximation, based only on outside convex cells 518 EnclosingBall<Sphere2D, S2Point> enclosingS2 = 519 new EnclosingBall<Sphere2D, S2Point>(S2Point.PLUS_K, Double.POSITIVE_INFINITY); 520 for (Vector3D outsidePoint : getOutsidePoints()) { 521 final S2Point outsideS2 = new S2Point(outsidePoint); 522 final BoundaryProjection<Sphere2D> projection = projectToBoundary(outsideS2); 523 if (FastMath.PI - projection.getOffset() < enclosingS2.getRadius()) { 524 enclosingS2 = new EnclosingBall<Sphere2D, S2Point>(outsideS2.negate(), 525 FastMath.PI - projection.getOffset(), 526 (S2Point) projection.getProjected()); 527 } 528 } 529 return enclosingS2; 530 } 531 final S2Point[] support = new S2Point[support3D.length]; 532 for (int i = 0; i < support3D.length; ++i) { 533 support[i] = new S2Point(support3D[i]); 534 } 535 536 final EnclosingBall<Sphere2D, S2Point> enclosingS2 = 537 new EnclosingBall<Sphere2D, S2Point>(new S2Point(enclosing3D.getCenter()), 538 FastMath.acos((1 + h * h - r * r) / (2 * h)), 539 support); 540 541 return enclosingS2; 542 543 } 544 545 /** Gather some inside points. 546 * @return list of points known to be strictly in all inside convex cells 547 */ getInsidePoints()548 private List<Vector3D> getInsidePoints() { 549 final PropertiesComputer pc = new PropertiesComputer(getTolerance()); 550 getTree(true).visit(pc); 551 return pc.getConvexCellsInsidePoints(); 552 } 553 554 /** Gather some outside points. 555 * @return list of points known to be strictly in all outside convex cells 556 */ getOutsidePoints()557 private List<Vector3D> getOutsidePoints() { 558 final SphericalPolygonsSet complement = 559 (SphericalPolygonsSet) new RegionFactory<Sphere2D>().getComplement(this); 560 final PropertiesComputer pc = new PropertiesComputer(getTolerance()); 561 complement.getTree(true).visit(pc); 562 return pc.getConvexCellsInsidePoints(); 563 } 564 565 } 566