001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math3.geometry.spherical.twod; 018 019import java.util.ArrayList; 020import java.util.Collection; 021import java.util.Collections; 022import java.util.Iterator; 023import java.util.List; 024 025import org.apache.commons.math3.exception.MathIllegalStateException; 026import org.apache.commons.math3.geometry.enclosing.EnclosingBall; 027import org.apache.commons.math3.geometry.enclosing.WelzlEncloser; 028import org.apache.commons.math3.geometry.euclidean.threed.Euclidean3D; 029import org.apache.commons.math3.geometry.euclidean.threed.Rotation; 030import org.apache.commons.math3.geometry.euclidean.threed.RotationConvention; 031import org.apache.commons.math3.geometry.euclidean.threed.SphereGenerator; 032import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; 033import org.apache.commons.math3.geometry.partitioning.AbstractRegion; 034import org.apache.commons.math3.geometry.partitioning.BSPTree; 035import org.apache.commons.math3.geometry.partitioning.BoundaryProjection; 036import org.apache.commons.math3.geometry.partitioning.RegionFactory; 037import org.apache.commons.math3.geometry.partitioning.SubHyperplane; 038import org.apache.commons.math3.geometry.spherical.oned.Sphere1D; 039import org.apache.commons.math3.util.FastMath; 040import org.apache.commons.math3.util.MathUtils; 041 042/** This class represents a region on the 2-sphere: a set of spherical polygons. 043 * @since 3.3 044 */ 045public class SphericalPolygonsSet extends AbstractRegion<Sphere2D, Sphere1D> { 046 047 /** Boundary defined as an array of closed loops start vertices. */ 048 private List<Vertex> loops; 049 050 /** Build a polygons set representing the whole real 2-sphere. 051 * @param tolerance below which points are consider to be identical 052 */ 053 public SphericalPolygonsSet(final double tolerance) { 054 super(tolerance); 055 } 056 057 /** Build a polygons set representing a hemisphere. 058 * @param pole pole of the hemisphere (the pole is in the inside half) 059 * @param tolerance below which points are consider to be identical 060 */ 061 public SphericalPolygonsSet(final Vector3D pole, final double tolerance) { 062 super(new BSPTree<Sphere2D>(new Circle(pole, tolerance).wholeHyperplane(), 063 new BSPTree<Sphere2D>(Boolean.FALSE), 064 new BSPTree<Sphere2D>(Boolean.TRUE), 065 null), 066 tolerance); 067 } 068 069 /** Build a polygons set representing a regular polygon. 070 * @param center center of the polygon (the center is in the inside half) 071 * @param meridian point defining the reference meridian for first polygon vertex 072 * @param outsideRadius distance of the vertices to the center 073 * @param n number of sides of the polygon 074 * @param tolerance below which points are consider to be identical 075 */ 076 public SphericalPolygonsSet(final Vector3D center, final Vector3D meridian, 077 final double outsideRadius, final int n, 078 final double tolerance) { 079 this(tolerance, createRegularPolygonVertices(center, meridian, outsideRadius, n)); 080 } 081 082 /** Build a polygons set from a BSP tree. 083 * <p>The leaf nodes of the BSP tree <em>must</em> have a 084 * {@code Boolean} attribute representing the inside status of 085 * the corresponding cell (true for inside cells, false for outside 086 * cells). In order to avoid building too many small objects, it is 087 * recommended to use the predefined constants 088 * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p> 089 * @param tree inside/outside BSP tree representing the region 090 * @param tolerance below which points are consider to be identical 091 */ 092 public SphericalPolygonsSet(final BSPTree<Sphere2D> tree, final double tolerance) { 093 super(tree, tolerance); 094 } 095 096 /** Build a polygons set from a Boundary REPresentation (B-rep). 097 * <p>The boundary is provided as a collection of {@link 098 * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the 099 * 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 */ 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 */ 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 */ 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 */ 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 */ 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 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 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 */ 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 */ 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 */ 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 */ 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}