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     * &pi; 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}