diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/geometry/spherical/twod')
10 files changed, 1998 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java new file mode 100644 index 0000000..a34db6d --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java @@ -0,0 +1,326 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.spherical.twod; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.euclidean.threed.Rotation; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.geometry.partitioning.Embedding; +import org.apache.commons.math3.geometry.partitioning.Hyperplane; +import org.apache.commons.math3.geometry.partitioning.SubHyperplane; +import org.apache.commons.math3.geometry.partitioning.Transform; +import org.apache.commons.math3.geometry.spherical.oned.Arc; +import org.apache.commons.math3.geometry.spherical.oned.ArcsSet; +import org.apache.commons.math3.geometry.spherical.oned.S1Point; +import org.apache.commons.math3.geometry.spherical.oned.Sphere1D; +import org.apache.commons.math3.util.FastMath; + +/** This class represents an oriented great circle on the 2-sphere. + + * <p>An oriented circle can be defined by a center point. The circle + * is the the set of points that are in the normal plan the center.</p> + + * <p>Since it is oriented the two spherical caps at its two sides are + * unambiguously identified as a left cap and a right cap. This can be + * used to identify the interior and the exterior in a simple way by + * local properties only when part of a line is used to define part of + * a spherical polygon boundary.</p> + + * @since 3.3 + */ +public class Circle implements Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1D> { + + /** Pole or circle center. */ + private Vector3D pole; + + /** First axis in the equator plane, origin of the phase angles. */ + private Vector3D x; + + /** Second axis in the equator plane, in quadrature with respect to x. */ + private Vector3D y; + + /** Tolerance below which close sub-arcs are merged together. */ + private final double tolerance; + + /** Build a great circle from its pole. + * <p>The circle is oriented in the trigonometric direction around pole.</p> + * @param pole circle pole + * @param tolerance tolerance below which close sub-arcs are merged together + */ + public Circle(final Vector3D pole, final double tolerance) { + reset(pole); + this.tolerance = tolerance; + } + + /** Build a great circle from two non-aligned points. + * <p>The circle is oriented from first to second point using the path smaller than \( \pi \).</p> + * @param first first point contained in the great circle + * @param second second point contained in the great circle + * @param tolerance tolerance below which close sub-arcs are merged together + */ + public Circle(final S2Point first, final S2Point second, final double tolerance) { + reset(first.getVector().crossProduct(second.getVector())); + this.tolerance = tolerance; + } + + /** Build a circle from its internal components. + * <p>The circle is oriented in the trigonometric direction around center.</p> + * @param pole circle pole + * @param x first axis in the equator plane + * @param y second axis in the equator plane + * @param tolerance tolerance below which close sub-arcs are merged together + */ + private Circle(final Vector3D pole, final Vector3D x, final Vector3D y, + final double tolerance) { + this.pole = pole; + this.x = x; + this.y = y; + this.tolerance = tolerance; + } + + /** Copy constructor. + * <p>The created instance is completely independent from the + * original instance, it is a deep copy.</p> + * @param circle circle to copy + */ + public Circle(final Circle circle) { + this(circle.pole, circle.x, circle.y, circle.tolerance); + } + + /** {@inheritDoc} */ + public Circle copySelf() { + return new Circle(this); + } + + /** Reset the instance as if built from a pole. + * <p>The circle is oriented in the trigonometric direction around pole.</p> + * @param newPole circle pole + */ + public void reset(final Vector3D newPole) { + this.pole = newPole.normalize(); + this.x = newPole.orthogonal(); + this.y = Vector3D.crossProduct(newPole, x).normalize(); + } + + /** Revert the instance. + */ + public void revertSelf() { + // x remains the same + y = y.negate(); + pole = pole.negate(); + } + + /** Get the reverse of the instance. + * <p>Get a circle with reversed orientation with respect to the + * instance. A new object is built, the instance is untouched.</p> + * @return a new circle, with orientation opposite to the instance orientation + */ + public Circle getReverse() { + return new Circle(pole.negate(), x, y.negate(), tolerance); + } + + /** {@inheritDoc} */ + public Point<Sphere2D> project(Point<Sphere2D> point) { + return toSpace(toSubSpace(point)); + } + + /** {@inheritDoc} */ + public double getTolerance() { + return tolerance; + } + + /** {@inheritDoc} + * @see #getPhase(Vector3D) + */ + public S1Point toSubSpace(final Point<Sphere2D> point) { + return new S1Point(getPhase(((S2Point) point).getVector())); + } + + /** Get the phase angle of a direction. + * <p> + * The direction may not belong to the circle as the + * phase is computed for the meridian plane between the circle + * pole and the direction. + * </p> + * @param direction direction for which phase is requested + * @return phase angle of the direction around the circle + * @see #toSubSpace(Point) + */ + public double getPhase(final Vector3D direction) { + return FastMath.PI + FastMath.atan2(-direction.dotProduct(y), -direction.dotProduct(x)); + } + + /** {@inheritDoc} + * @see #getPointAt(double) + */ + public S2Point toSpace(final Point<Sphere1D> point) { + return new S2Point(getPointAt(((S1Point) point).getAlpha())); + } + + /** Get a circle point from its phase around the circle. + * @param alpha phase around the circle + * @return circle point on the sphere + * @see #toSpace(Point) + * @see #getXAxis() + * @see #getYAxis() + */ + public Vector3D getPointAt(final double alpha) { + return new Vector3D(FastMath.cos(alpha), x, FastMath.sin(alpha), y); + } + + /** Get the X axis of the circle. + * <p> + * This method returns the same value as {@link #getPointAt(double) + * getPointAt(0.0)} but it does not do any computation and always + * return the same instance. + * </p> + * @return an arbitrary x axis on the circle + * @see #getPointAt(double) + * @see #getYAxis() + * @see #getPole() + */ + public Vector3D getXAxis() { + return x; + } + + /** Get the Y axis of the circle. + * <p> + * This method returns the same value as {@link #getPointAt(double) + * getPointAt(0.5 * FastMath.PI)} but it does not do any computation and always + * return the same instance. + * </p> + * @return an arbitrary y axis point on the circle + * @see #getPointAt(double) + * @see #getXAxis() + * @see #getPole() + */ + public Vector3D getYAxis() { + return y; + } + + /** Get the pole of the circle. + * <p> + * As the circle is a great circle, the pole does <em>not</em> + * belong to it. + * </p> + * @return pole of the circle + * @see #getXAxis() + * @see #getYAxis() + */ + public Vector3D getPole() { + return pole; + } + + /** Get the arc of the instance that lies inside the other circle. + * @param other other circle + * @return arc of the instance that lies inside the other circle + */ + public Arc getInsideArc(final Circle other) { + final double alpha = getPhase(other.pole); + final double halfPi = 0.5 * FastMath.PI; + return new Arc(alpha - halfPi, alpha + halfPi, tolerance); + } + + /** {@inheritDoc} */ + public SubCircle wholeHyperplane() { + return new SubCircle(this, new ArcsSet(tolerance)); + } + + /** Build a region covering the whole space. + * @return a region containing the instance (really a {@link + * SphericalPolygonsSet SphericalPolygonsSet} instance) + */ + public SphericalPolygonsSet wholeSpace() { + return new SphericalPolygonsSet(tolerance); + } + + /** {@inheritDoc} + * @see #getOffset(Vector3D) + */ + public double getOffset(final Point<Sphere2D> point) { + return getOffset(((S2Point) point).getVector()); + } + + /** Get the offset (oriented distance) of a direction. + * <p>The offset is defined as the angular distance between the + * circle center and the direction minus the circle radius. It + * is therefore 0 on the circle, positive for directions outside of + * the cone delimited by the circle, and negative inside the cone.</p> + * @param direction direction to check + * @return offset of the direction + * @see #getOffset(Point) + */ + public double getOffset(final Vector3D direction) { + return Vector3D.angle(pole, direction) - 0.5 * FastMath.PI; + } + + /** {@inheritDoc} */ + public boolean sameOrientationAs(final Hyperplane<Sphere2D> other) { + final Circle otherC = (Circle) other; + return Vector3D.dotProduct(pole, otherC.pole) >= 0.0; + } + + /** Get a {@link org.apache.commons.math3.geometry.partitioning.Transform + * Transform} embedding a 3D rotation. + * @param rotation rotation to use + * @return a new transform that can be applied to either {@link + * Point Point}, {@link Circle Line} or {@link + * org.apache.commons.math3.geometry.partitioning.SubHyperplane + * SubHyperplane} instances + */ + public static Transform<Sphere2D, Sphere1D> getTransform(final Rotation rotation) { + return new CircleTransform(rotation); + } + + /** Class embedding a 3D rotation. */ + private static class CircleTransform implements Transform<Sphere2D, Sphere1D> { + + /** Underlying rotation. */ + private final Rotation rotation; + + /** Build a transform from a {@code Rotation}. + * @param rotation rotation to use + */ + CircleTransform(final Rotation rotation) { + this.rotation = rotation; + } + + /** {@inheritDoc} */ + public S2Point apply(final Point<Sphere2D> point) { + return new S2Point(rotation.applyTo(((S2Point) point).getVector())); + } + + /** {@inheritDoc} */ + public Circle apply(final Hyperplane<Sphere2D> hyperplane) { + final Circle circle = (Circle) hyperplane; + return new Circle(rotation.applyTo(circle.pole), + rotation.applyTo(circle.x), + rotation.applyTo(circle.y), + circle.tolerance); + } + + /** {@inheritDoc} */ + public SubHyperplane<Sphere1D> apply(final SubHyperplane<Sphere1D> sub, + final Hyperplane<Sphere2D> original, + final Hyperplane<Sphere2D> transformed) { + // as the circle is rotated, the limit angles are rotated too + return sub; + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Edge.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Edge.java new file mode 100644 index 0000000..a9ccb08 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Edge.java @@ -0,0 +1,222 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.spherical.twod; + +import java.util.List; + +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.geometry.spherical.oned.Arc; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; + +/** Spherical polygons boundary edge. + * @see SphericalPolygonsSet#getBoundaryLoops() + * @see Vertex + * @since 3.3 + */ +public class Edge { + + /** Start vertex. */ + private final Vertex start; + + /** End vertex. */ + private Vertex end; + + /** Length of the arc. */ + private final double length; + + /** Circle supporting the edge. */ + private final Circle circle; + + /** Build an edge not contained in any node yet. + * @param start start vertex + * @param end end vertex + * @param length length of the arc (it can be greater than \( \pi \)) + * @param circle circle supporting the edge + */ + Edge(final Vertex start, final Vertex end, final double length, final Circle circle) { + + this.start = start; + this.end = end; + this.length = length; + this.circle = circle; + + // connect the vertices back to the edge + start.setOutgoing(this); + end.setIncoming(this); + + } + + /** Get start vertex. + * @return start vertex + */ + public Vertex getStart() { + return start; + } + + /** Get end vertex. + * @return end vertex + */ + public Vertex getEnd() { + return end; + } + + /** Get the length of the arc. + * @return length of the arc (can be greater than \( \pi \)) + */ + public double getLength() { + return length; + } + + /** Get the circle supporting this edge. + * @return circle supporting this edge + */ + public Circle getCircle() { + return circle; + } + + /** Get an intermediate point. + * <p> + * The angle along the edge should normally be between 0 and {@link #getLength()} + * in order to remain within edge limits. However, there are no checks on the + * value of the angle, so user can rebuild the full circle on which an edge is + * defined if they want. + * </p> + * @param alpha angle along the edge, counted from {@link #getStart()} + * @return an intermediate point + */ + public Vector3D getPointAt(final double alpha) { + return circle.getPointAt(alpha + circle.getPhase(start.getLocation().getVector())); + } + + /** Connect the instance with a following edge. + * @param next edge following the instance + */ + void setNextEdge(final Edge next) { + end = next.getStart(); + end.setIncoming(this); + end.bindWith(getCircle()); + } + + /** Split the edge. + * <p> + * Once split, this edge is not referenced anymore by the vertices, + * it is replaced by the two or three sub-edges and intermediate splitting + * vertices are introduced to connect these sub-edges together. + * </p> + * @param splitCircle circle splitting the edge in several parts + * @param outsideList list where to put parts that are outside of the split circle + * @param insideList list where to put parts that are inside the split circle + */ + void split(final Circle splitCircle, + final List<Edge> outsideList, final List<Edge> insideList) { + + // get the inside arc, synchronizing its phase with the edge itself + final double edgeStart = circle.getPhase(start.getLocation().getVector()); + final Arc arc = circle.getInsideArc(splitCircle); + final double arcRelativeStart = MathUtils.normalizeAngle(arc.getInf(), edgeStart + FastMath.PI) - edgeStart; + final double arcRelativeEnd = arcRelativeStart + arc.getSize(); + final double unwrappedEnd = arcRelativeEnd - MathUtils.TWO_PI; + + // build the sub-edges + final double tolerance = circle.getTolerance(); + Vertex previousVertex = start; + if (unwrappedEnd >= length - tolerance) { + + // the edge is entirely contained inside the circle + // we don't split anything + insideList.add(this); + + } else { + + // there are at least some parts of the edge that should be outside + // (even is they are later be filtered out as being too small) + double alreadyManagedLength = 0; + if (unwrappedEnd >= 0) { + // the start of the edge is inside the circle + previousVertex = addSubEdge(previousVertex, + new Vertex(new S2Point(circle.getPointAt(edgeStart + unwrappedEnd))), + unwrappedEnd, insideList, splitCircle); + alreadyManagedLength = unwrappedEnd; + } + + if (arcRelativeStart >= length - tolerance) { + // the edge ends while still outside of the circle + if (unwrappedEnd >= 0) { + previousVertex = addSubEdge(previousVertex, end, + length - alreadyManagedLength, outsideList, splitCircle); + } else { + // the edge is entirely outside of the circle + // we don't split anything + outsideList.add(this); + } + } else { + // the edge is long enough to enter inside the circle + previousVertex = addSubEdge(previousVertex, + new Vertex(new S2Point(circle.getPointAt(edgeStart + arcRelativeStart))), + arcRelativeStart - alreadyManagedLength, outsideList, splitCircle); + alreadyManagedLength = arcRelativeStart; + + if (arcRelativeEnd >= length - tolerance) { + // the edge ends while still inside of the circle + previousVertex = addSubEdge(previousVertex, end, + length - alreadyManagedLength, insideList, splitCircle); + } else { + // the edge is long enough to exit outside of the circle + previousVertex = addSubEdge(previousVertex, + new Vertex(new S2Point(circle.getPointAt(edgeStart + arcRelativeStart))), + arcRelativeStart - alreadyManagedLength, insideList, splitCircle); + alreadyManagedLength = arcRelativeStart; + previousVertex = addSubEdge(previousVertex, end, + length - alreadyManagedLength, outsideList, splitCircle); + } + } + + } + + } + + /** Add a sub-edge to a list if long enough. + * <p> + * If the length of the sub-edge to add is smaller than the {@link Circle#getTolerance()} + * tolerance of the support circle, it will be ignored. + * </p> + * @param subStart start of the sub-edge + * @param subEnd end of the sub-edge + * @param subLength length of the sub-edge + * @param splitCircle circle splitting the edge in several parts + * @param list list where to put the sub-edge + * @return end vertex of the edge ({@code subEnd} if the edge was long enough and really + * added, {@code subStart} if the edge was too small and therefore ignored) + */ + private Vertex addSubEdge(final Vertex subStart, final Vertex subEnd, final double subLength, + final List<Edge> list, final Circle splitCircle) { + + if (subLength <= circle.getTolerance()) { + // the edge is too short, we ignore it + return subStart; + } + + // really add the edge + subEnd.bindWith(splitCircle); + final Edge edge = new Edge(subStart, subEnd, subLength, circle); + list.add(edge); + return subEnd; + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/EdgesBuilder.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/EdgesBuilder.java new file mode 100644 index 0000000..844cfb1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/EdgesBuilder.java @@ -0,0 +1,169 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.spherical.twod; + +import java.util.ArrayList; +import java.util.IdentityHashMap; +import java.util.List; +import java.util.Map; + +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.geometry.partitioning.BSPTree; +import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor; +import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute; +import org.apache.commons.math3.geometry.spherical.oned.Arc; +import org.apache.commons.math3.geometry.spherical.oned.ArcsSet; +import org.apache.commons.math3.geometry.spherical.oned.S1Point; + +/** Visitor building edges. + * @since 3.3 + */ +class EdgesBuilder implements BSPTreeVisitor<Sphere2D> { + + /** Root of the tree. */ + private final BSPTree<Sphere2D> root; + + /** Tolerance below which points are consider to be identical. */ + private final double tolerance; + + /** Built edges and their associated nodes. */ + private final Map<Edge, BSPTree<Sphere2D>> edgeToNode; + + /** Reversed map. */ + private final Map<BSPTree<Sphere2D>, List<Edge>> nodeToEdgesList; + + /** Simple constructor. + * @param root tree root + * @param tolerance below which points are consider to be identical + */ + EdgesBuilder(final BSPTree<Sphere2D> root, final double tolerance) { + this.root = root; + this.tolerance = tolerance; + this.edgeToNode = new IdentityHashMap<Edge, BSPTree<Sphere2D>>(); + this.nodeToEdgesList = new IdentityHashMap<BSPTree<Sphere2D>, List<Edge>>(); + } + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree<Sphere2D> node) { + return Order.MINUS_SUB_PLUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree<Sphere2D> node) { + nodeToEdgesList.put(node, new ArrayList<Edge>()); + @SuppressWarnings("unchecked") + final BoundaryAttribute<Sphere2D> attribute = (BoundaryAttribute<Sphere2D>) node.getAttribute(); + if (attribute.getPlusOutside() != null) { + addContribution((SubCircle) attribute.getPlusOutside(), false, node); + } + if (attribute.getPlusInside() != null) { + addContribution((SubCircle) attribute.getPlusInside(), true, node); + } + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree<Sphere2D> node) { + } + + /** Add the contribution of a boundary edge. + * @param sub boundary facet + * @param reversed if true, the facet has the inside on its plus side + * @param node node to which the edge belongs + */ + private void addContribution(final SubCircle sub, final boolean reversed, + final BSPTree<Sphere2D> node) { + final Circle circle = (Circle) sub.getHyperplane(); + final List<Arc> arcs = ((ArcsSet) sub.getRemainingRegion()).asList(); + for (final Arc a : arcs) { + final Vertex start = new Vertex((S2Point) circle.toSpace(new S1Point(a.getInf()))); + final Vertex end = new Vertex((S2Point) circle.toSpace(new S1Point(a.getSup()))); + start.bindWith(circle); + end.bindWith(circle); + final Edge edge; + if (reversed) { + edge = new Edge(end, start, a.getSize(), circle.getReverse()); + } else { + edge = new Edge(start, end, a.getSize(), circle); + } + edgeToNode.put(edge, node); + nodeToEdgesList.get(node).add(edge); + } + } + + /** Get the edge that should naturally follow another one. + * @param previous edge to be continued + * @return other edge, starting where the previous one ends (they + * have not been connected yet) + * @exception MathIllegalStateException if there is not a single other edge + */ + private Edge getFollowingEdge(final Edge previous) + throws MathIllegalStateException { + + // get the candidate nodes + final S2Point point = previous.getEnd().getLocation(); + final List<BSPTree<Sphere2D>> candidates = root.getCloseCuts(point, tolerance); + + // the following edge we are looking for must start from one of the candidates nodes + double closest = tolerance; + Edge following = null; + for (final BSPTree<Sphere2D> node : candidates) { + for (final Edge edge : nodeToEdgesList.get(node)) { + if (edge != previous && edge.getStart().getIncoming() == null) { + final Vector3D edgeStart = edge.getStart().getLocation().getVector(); + final double gap = Vector3D.angle(point.getVector(), edgeStart); + if (gap <= closest) { + closest = gap; + following = edge; + } + } + } + } + + if (following == null) { + final Vector3D previousStart = previous.getStart().getLocation().getVector(); + if (Vector3D.angle(point.getVector(), previousStart) <= tolerance) { + // the edge connects back to itself + return previous; + } + + // this should never happen + throw new MathIllegalStateException(LocalizedFormats.OUTLINE_BOUNDARY_LOOP_OPEN); + + } + + return following; + + } + + /** Get the boundary edges. + * @return boundary edges + * @exception MathIllegalStateException if there is not a single other edge + */ + public List<Edge> getEdges() throws MathIllegalStateException { + + // connect the edges + for (final Edge previous : edgeToNode.keySet()) { + previous.setNextEdge(getFollowingEdge(previous)); + } + + return new ArrayList<Edge>(edgeToNode.keySet()); + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/PropertiesComputer.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/PropertiesComputer.java new file mode 100644 index 0000000..593180f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/PropertiesComputer.java @@ -0,0 +1,173 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.spherical.twod; + +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.geometry.partitioning.BSPTree; +import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; + +/** Visitor computing geometrical properties. + * @since 3.3 + */ +class PropertiesComputer implements BSPTreeVisitor<Sphere2D> { + + /** Tolerance below which points are consider to be identical. */ + private final double tolerance; + + /** Summed area. */ + private double summedArea; + + /** Summed barycenter. */ + private Vector3D summedBarycenter; + + /** List of points strictly inside convex cells. */ + private final List<Vector3D> convexCellsInsidePoints; + + /** Simple constructor. + * @param tolerance below which points are consider to be identical + */ + PropertiesComputer(final double tolerance) { + this.tolerance = tolerance; + this.summedArea = 0; + this.summedBarycenter = Vector3D.ZERO; + this.convexCellsInsidePoints = new ArrayList<Vector3D>(); + } + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree<Sphere2D> node) { + return Order.MINUS_SUB_PLUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree<Sphere2D> node) { + // nothing to do here + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree<Sphere2D> node) { + if ((Boolean) node.getAttribute()) { + + // transform this inside leaf cell into a simple convex polygon + final SphericalPolygonsSet convex = + new SphericalPolygonsSet(node.pruneAroundConvexCell(Boolean.TRUE, + Boolean.FALSE, + null), + tolerance); + + // extract the start of the single loop boundary of the convex cell + final List<Vertex> boundary = convex.getBoundaryLoops(); + if (boundary.size() != 1) { + // this should never happen + throw new MathInternalError(); + } + + // compute the geometrical properties of the convex cell + final double area = convexCellArea(boundary.get(0)); + final Vector3D barycenter = convexCellBarycenter(boundary.get(0)); + convexCellsInsidePoints.add(barycenter); + + // add the cell contribution to the global properties + summedArea += area; + summedBarycenter = new Vector3D(1, summedBarycenter, area, barycenter); + + } + } + + /** Compute convex cell area. + * @param start start vertex of the convex cell boundary + * @return area + */ + private double convexCellArea(final Vertex start) { + + int n = 0; + double sum = 0; + + // loop around the cell + for (Edge e = start.getOutgoing(); n == 0 || e.getStart() != start; e = e.getEnd().getOutgoing()) { + + // find path interior angle at vertex + final Vector3D previousPole = e.getCircle().getPole(); + final Vector3D nextPole = e.getEnd().getOutgoing().getCircle().getPole(); + final Vector3D point = e.getEnd().getLocation().getVector(); + double alpha = FastMath.atan2(Vector3D.dotProduct(nextPole, Vector3D.crossProduct(point, previousPole)), + -Vector3D.dotProduct(nextPole, previousPole)); + if (alpha < 0) { + alpha += MathUtils.TWO_PI; + } + sum += alpha; + n++; + } + + // compute area using extended Girard theorem + // see Spherical Trigonometry: For the Use of Colleges and Schools by I. Todhunter + // article 99 in chapter VIII Area Of a Spherical Triangle. Spherical Excess. + // book available from project Gutenberg at http://www.gutenberg.org/ebooks/19770 + return sum - (n - 2) * FastMath.PI; + + } + + /** Compute convex cell barycenter. + * @param start start vertex of the convex cell boundary + * @return barycenter + */ + private Vector3D convexCellBarycenter(final Vertex start) { + + int n = 0; + Vector3D sumB = Vector3D.ZERO; + + // loop around the cell + for (Edge e = start.getOutgoing(); n == 0 || e.getStart() != start; e = e.getEnd().getOutgoing()) { + sumB = new Vector3D(1, sumB, e.getLength(), e.getCircle().getPole()); + n++; + } + + return sumB.normalize(); + + } + + /** Get the area. + * @return area + */ + public double getArea() { + return summedArea; + } + + /** Get the barycenter. + * @return barycenter + */ + public S2Point getBarycenter() { + if (summedBarycenter.getNormSq() == 0) { + return S2Point.NaN; + } else { + return new S2Point(summedBarycenter); + } + } + + /** Get the points strictly inside convex cells. + * @return points strictly inside convex cells + */ + public List<Vector3D> getConvexCellsInsidePoints() { + return convexCellsInsidePoints; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java new file mode 100644 index 0000000..677e830 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java @@ -0,0 +1,237 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.spherical.twod; + +import org.apache.commons.math3.exception.MathArithmeticException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; + +/** This class represents a point on the 2-sphere. + * <p> + * We use the mathematical convention to use the azimuthal angle \( \theta \) + * in the x-y plane as the first coordinate, and the polar angle \( \varphi \) + * as the second coordinate (see <a + * href="http://mathworld.wolfram.com/SphericalCoordinates.html">Spherical + * Coordinates</a> in MathWorld). + * </p> + * <p>Instances of this class are guaranteed to be immutable.</p> + * @since 3.3 + */ +public class S2Point implements Point<Sphere2D> { + + /** +I (coordinates: \( \theta = 0, \varphi = \pi/2 \)). */ + public static final S2Point PLUS_I = new S2Point(0, 0.5 * FastMath.PI, Vector3D.PLUS_I); + + /** +J (coordinates: \( \theta = \pi/2, \varphi = \pi/2 \))). */ + public static final S2Point PLUS_J = new S2Point(0.5 * FastMath.PI, 0.5 * FastMath.PI, Vector3D.PLUS_J); + + /** +K (coordinates: \( \theta = any angle, \varphi = 0 \)). */ + public static final S2Point PLUS_K = new S2Point(0, 0, Vector3D.PLUS_K); + + /** -I (coordinates: \( \theta = \pi, \varphi = \pi/2 \)). */ + public static final S2Point MINUS_I = new S2Point(FastMath.PI, 0.5 * FastMath.PI, Vector3D.MINUS_I); + + /** -J (coordinates: \( \theta = 3\pi/2, \varphi = \pi/2 \)). */ + public static final S2Point MINUS_J = new S2Point(1.5 * FastMath.PI, 0.5 * FastMath.PI, Vector3D.MINUS_J); + + /** -K (coordinates: \( \theta = any angle, \varphi = \pi \)). */ + public static final S2Point MINUS_K = new S2Point(0, FastMath.PI, Vector3D.MINUS_K); + + // CHECKSTYLE: stop ConstantName + /** A vector with all coordinates set to NaN. */ + public static final S2Point NaN = new S2Point(Double.NaN, Double.NaN, Vector3D.NaN); + // CHECKSTYLE: resume ConstantName + + /** Serializable UID. */ + private static final long serialVersionUID = 20131218L; + + /** Azimuthal angle \( \theta \) in the x-y plane. */ + private final double theta; + + /** Polar angle \( \varphi \). */ + private final double phi; + + /** Corresponding 3D normalized vector. */ + private final Vector3D vector; + + /** Simple constructor. + * Build a vector from its spherical coordinates + * @param theta azimuthal angle \( \theta \) in the x-y plane + * @param phi polar angle \( \varphi \) + * @see #getTheta() + * @see #getPhi() + * @exception OutOfRangeException if \( \varphi \) is not in the [\( 0; \pi \)] range + */ + public S2Point(final double theta, final double phi) + throws OutOfRangeException { + this(theta, phi, vector(theta, phi)); + } + + /** Simple constructor. + * Build a vector from its underlying 3D vector + * @param vector 3D vector + * @exception MathArithmeticException if vector norm is zero + */ + public S2Point(final Vector3D vector) throws MathArithmeticException { + this(FastMath.atan2(vector.getY(), vector.getX()), Vector3D.angle(Vector3D.PLUS_K, vector), + vector.normalize()); + } + + /** Build a point from its internal components. + * @param theta azimuthal angle \( \theta \) in the x-y plane + * @param phi polar angle \( \varphi \) + * @param vector corresponding vector + */ + private S2Point(final double theta, final double phi, final Vector3D vector) { + this.theta = theta; + this.phi = phi; + this.vector = vector; + } + + /** Build the normalized vector corresponding to spherical coordinates. + * @param theta azimuthal angle \( \theta \) in the x-y plane + * @param phi polar angle \( \varphi \) + * @return normalized vector + * @exception OutOfRangeException if \( \varphi \) is not in the [\( 0; \pi \)] range + */ + private static Vector3D vector(final double theta, final double phi) + throws OutOfRangeException { + + if (phi < 0 || phi > FastMath.PI) { + throw new OutOfRangeException(phi, 0, FastMath.PI); + } + + final double cosTheta = FastMath.cos(theta); + final double sinTheta = FastMath.sin(theta); + final double cosPhi = FastMath.cos(phi); + final double sinPhi = FastMath.sin(phi); + + return new Vector3D(cosTheta * sinPhi, sinTheta * sinPhi, cosPhi); + + } + + /** Get the azimuthal angle \( \theta \) in the x-y plane. + * @return azimuthal angle \( \theta \) in the x-y plane + * @see #S2Point(double, double) + */ + public double getTheta() { + return theta; + } + + /** Get the polar angle \( \varphi \). + * @return polar angle \( \varphi \) + * @see #S2Point(double, double) + */ + public double getPhi() { + return phi; + } + + /** Get the corresponding normalized vector in the 3D euclidean space. + * @return normalized vector + */ + public Vector3D getVector() { + return vector; + } + + /** {@inheritDoc} */ + public Space getSpace() { + return Sphere2D.getInstance(); + } + + /** {@inheritDoc} */ + public boolean isNaN() { + return Double.isNaN(theta) || Double.isNaN(phi); + } + + /** Get the opposite of the instance. + * @return a new vector which is opposite to the instance + */ + public S2Point negate() { + return new S2Point(-theta, FastMath.PI - phi, vector.negate()); + } + + /** {@inheritDoc} */ + public double distance(final Point<Sphere2D> point) { + return distance(this, (S2Point) point); + } + + /** Compute the distance (angular separation) between two points. + * @param p1 first vector + * @param p2 second vector + * @return the angular separation between p1 and p2 + */ + public static double distance(S2Point p1, S2Point p2) { + return Vector3D.angle(p1.vector, p2.vector); + } + + /** + * Test for the equality of two points on the 2-sphere. + * <p> + * If all coordinates of two points are exactly the same, and none are + * <code>Double.NaN</code>, the two points are considered to be equal. + * </p> + * <p> + * <code>NaN</code> coordinates are considered to affect globally the vector + * and be equals to each other - i.e, if either (or all) coordinates of the + * 2D vector are equal to <code>Double.NaN</code>, the 2D vector is equal to + * {@link #NaN}. + * </p> + * + * @param other Object to test for equality to this + * @return true if two points on the 2-sphere objects are equal, false if + * object is null, not an instance of S2Point, or + * not equal to this S2Point instance + * + */ + @Override + public boolean equals(Object other) { + + if (this == other) { + return true; + } + + if (other instanceof S2Point) { + final S2Point rhs = (S2Point) other; + if (rhs.isNaN()) { + return this.isNaN(); + } + + return (theta == rhs.theta) && (phi == rhs.phi); + } + return false; + } + + /** + * Get a hashCode for the 2D vector. + * <p> + * All NaN values have the same hash code.</p> + * + * @return a hash code value for this object + */ + @Override + public int hashCode() { + if (isNaN()) { + return 542; + } + return 134 * (37 * MathUtils.hash(theta) + MathUtils.hash(phi)); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Sphere2D.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Sphere2D.java new file mode 100644 index 0000000..93ff04e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Sphere2D.java @@ -0,0 +1,80 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.geometry.spherical.twod; + +import java.io.Serializable; + +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.spherical.oned.Sphere1D; + +/** + * This class implements a two-dimensional sphere (i.e. the regular sphere). + * <p> + * We use here the topologists definition of the 2-sphere (see + * <a href="http://mathworld.wolfram.com/Sphere.html">Sphere</a> on + * MathWorld), i.e. the 2-sphere is the two-dimensional surface + * defined in 3D as x<sup>2</sup>+y<sup>2</sup>+z<sup>2</sup>=1. + * </p> + * @since 3.3 + */ +public class Sphere2D implements Serializable, Space { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 20131218L; + + /** Private constructor for the singleton. + */ + private Sphere2D() { + } + + /** Get the unique instance. + * @return the unique instance + */ + public static Sphere2D getInstance() { + return LazyHolder.INSTANCE; + } + + /** {@inheritDoc} */ + public int getDimension() { + return 2; + } + + /** {@inheritDoc} */ + public Sphere1D getSubSpace() { + return Sphere1D.getInstance(); + } + + // CHECKSTYLE: stop HideUtilityClassConstructor + /** Holder for the instance. + * <p>We use here the Initialization On Demand Holder Idiom.</p> + */ + private static class LazyHolder { + /** Cached field instance. */ + private static final Sphere2D INSTANCE = new Sphere2D(); + } + // CHECKSTYLE: resume HideUtilityClassConstructor + + /** Handle deserialization of the singleton. + * @return the singleton instance + */ + private Object readResolve() { + // return the singleton instance + return LazyHolder.INSTANCE; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java new file mode 100644 index 0000000..8e41659 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java @@ -0,0 +1,565 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.spherical.twod; + +import java.util.ArrayList; +import java.util.Collection; +import java.util.Collections; +import java.util.Iterator; +import java.util.List; + +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.geometry.enclosing.EnclosingBall; +import org.apache.commons.math3.geometry.enclosing.WelzlEncloser; +import org.apache.commons.math3.geometry.euclidean.threed.Euclidean3D; +import org.apache.commons.math3.geometry.euclidean.threed.Rotation; +import org.apache.commons.math3.geometry.euclidean.threed.RotationConvention; +import org.apache.commons.math3.geometry.euclidean.threed.SphereGenerator; +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.geometry.partitioning.AbstractRegion; +import org.apache.commons.math3.geometry.partitioning.BSPTree; +import org.apache.commons.math3.geometry.partitioning.BoundaryProjection; +import org.apache.commons.math3.geometry.partitioning.RegionFactory; +import org.apache.commons.math3.geometry.partitioning.SubHyperplane; +import org.apache.commons.math3.geometry.spherical.oned.Sphere1D; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; + +/** This class represents a region on the 2-sphere: a set of spherical polygons. + * @since 3.3 + */ +public class SphericalPolygonsSet extends AbstractRegion<Sphere2D, Sphere1D> { + + /** Boundary defined as an array of closed loops start vertices. */ + private List<Vertex> loops; + + /** Build a polygons set representing the whole real 2-sphere. + * @param tolerance below which points are consider to be identical + */ + public SphericalPolygonsSet(final double tolerance) { + super(tolerance); + } + + /** Build a polygons set representing a hemisphere. + * @param pole pole of the hemisphere (the pole is in the inside half) + * @param tolerance below which points are consider to be identical + */ + public SphericalPolygonsSet(final Vector3D pole, final double tolerance) { + super(new BSPTree<Sphere2D>(new Circle(pole, tolerance).wholeHyperplane(), + new BSPTree<Sphere2D>(Boolean.FALSE), + new BSPTree<Sphere2D>(Boolean.TRUE), + null), + tolerance); + } + + /** Build a polygons set representing a regular polygon. + * @param center center of the polygon (the center is in the inside half) + * @param meridian point defining the reference meridian for first polygon vertex + * @param outsideRadius distance of the vertices to the center + * @param n number of sides of the polygon + * @param tolerance below which points are consider to be identical + */ + public SphericalPolygonsSet(final Vector3D center, final Vector3D meridian, + final double outsideRadius, final int n, + final double tolerance) { + this(tolerance, createRegularPolygonVertices(center, meridian, outsideRadius, n)); + } + + /** Build a polygons set from a BSP tree. + * <p>The leaf nodes of the BSP tree <em>must</em> have a + * {@code Boolean} attribute representing the inside status of + * the corresponding cell (true for inside cells, false for outside + * cells). In order to avoid building too many small objects, it is + * recommended to use the predefined constants + * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p> + * @param tree inside/outside BSP tree representing the region + * @param tolerance below which points are consider to be identical + */ + public SphericalPolygonsSet(final BSPTree<Sphere2D> tree, final double tolerance) { + super(tree, tolerance); + } + + /** Build a polygons set from a Boundary REPresentation (B-rep). + * <p>The boundary is provided as a collection of {@link + * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the + * interior part of the region on its minus side and the exterior on + * its plus side.</p> + * <p>The boundary elements can be in any order, and can form + * several non-connected sets (like for example polygons with holes + * or a set of disjoint polygons considered as a whole). In + * fact, the elements do not even need to be connected together + * (their topological connections are not used here). However, if the + * boundary does not really separate an inside open from an outside + * open (open having here its topological meaning), then subsequent + * calls to the {@link + * org.apache.commons.math3.geometry.partitioning.Region#checkPoint(org.apache.commons.math3.geometry.Point) + * checkPoint} method will not be meaningful anymore.</p> + * <p>If the boundary is empty, the region will represent the whole + * space.</p> + * @param boundary collection of boundary elements, as a + * collection of {@link SubHyperplane SubHyperplane} objects + * @param tolerance below which points are consider to be identical + */ + public SphericalPolygonsSet(final Collection<SubHyperplane<Sphere2D>> boundary, final double tolerance) { + super(boundary, tolerance); + } + + /** Build a polygon from a simple list of vertices. + * <p>The boundary is provided as a list of points considering to + * represent the vertices of a simple loop. The interior part of the + * region is on the left side of this path and the exterior is on its + * right side.</p> + * <p>This constructor does not handle polygons with a boundary + * forming several disconnected paths (such as polygons with holes).</p> + * <p>For cases where this simple constructor applies, it is expected to + * be numerically more robust than the {@link #SphericalPolygonsSet(Collection, + * double) general constructor} using {@link SubHyperplane subhyperplanes}.</p> + * <p>If the list is empty, the region will represent the whole + * space.</p> + * <p> + * Polygons with thin pikes or dents are inherently difficult to handle because + * they involve circles with almost opposite directions at some vertices. Polygons + * whose vertices come from some physical measurement with noise are also + * difficult because an edge that should be straight may be broken in lots of + * different pieces with almost equal directions. In both cases, computing the + * circles intersections is not numerically robust due to the almost 0 or almost + * π angle. Such cases need to carefully adjust the {@code hyperplaneThickness} + * parameter. A too small value would often lead to completely wrong polygons + * with large area wrongly identified as inside or outside. Large values are + * often much safer. As a rule of thumb, a value slightly below the size of the + * most accurate detail needed is a good value for the {@code hyperplaneThickness} + * parameter. + * </p> + * @param hyperplaneThickness tolerance below which points are considered to + * belong to the hyperplane (which is therefore more a slab) + * @param vertices vertices of the simple loop boundary + */ + public SphericalPolygonsSet(final double hyperplaneThickness, final S2Point ... vertices) { + super(verticesToTree(hyperplaneThickness, vertices), hyperplaneThickness); + } + + /** Build the vertices representing a regular polygon. + * @param center center of the polygon (the center is in the inside half) + * @param meridian point defining the reference meridian for first polygon vertex + * @param outsideRadius distance of the vertices to the center + * @param n number of sides of the polygon + * @return vertices array + */ + private static S2Point[] createRegularPolygonVertices(final Vector3D center, final Vector3D meridian, + final double outsideRadius, final int n) { + final S2Point[] array = new S2Point[n]; + final Rotation r0 = new Rotation(Vector3D.crossProduct(center, meridian), + outsideRadius, RotationConvention.VECTOR_OPERATOR); + array[0] = new S2Point(r0.applyTo(center)); + + final Rotation r = new Rotation(center, MathUtils.TWO_PI / n, RotationConvention.VECTOR_OPERATOR); + for (int i = 1; i < n; ++i) { + array[i] = new S2Point(r.applyTo(array[i - 1].getVector())); + } + + return array; + } + + /** Build the BSP tree of a polygons set from a simple list of vertices. + * <p>The boundary is provided as a list of points considering to + * represent the vertices of a simple loop. The interior part of the + * region is on the left side of this path and the exterior is on its + * right side.</p> + * <p>This constructor does not handle polygons with a boundary + * forming several disconnected paths (such as polygons with holes).</p> + * <p>This constructor handles only polygons with edges strictly shorter + * than \( \pi \). If longer edges are needed, they need to be broken up + * in smaller sub-edges so this constraint holds.</p> + * <p>For cases where this simple constructor applies, it is expected to + * be numerically more robust than the {@link #PolygonsSet(Collection) general + * constructor} using {@link SubHyperplane subhyperplanes}.</p> + * @param hyperplaneThickness tolerance below which points are consider to + * belong to the hyperplane (which is therefore more a slab) + * @param vertices vertices of the simple loop boundary + * @return the BSP tree of the input vertices + */ + private static BSPTree<Sphere2D> verticesToTree(final double hyperplaneThickness, + final S2Point ... vertices) { + + final int n = vertices.length; + if (n == 0) { + // the tree represents the whole space + return new BSPTree<Sphere2D>(Boolean.TRUE); + } + + // build the vertices + final Vertex[] vArray = new Vertex[n]; + for (int i = 0; i < n; ++i) { + vArray[i] = new Vertex(vertices[i]); + } + + // build the edges + List<Edge> edges = new ArrayList<Edge>(n); + Vertex end = vArray[n - 1]; + for (int i = 0; i < n; ++i) { + + // get the endpoints of the edge + final Vertex start = end; + end = vArray[i]; + + // get the circle supporting the edge, taking care not to recreate it + // if it was already created earlier due to another edge being aligned + // with the current one + Circle circle = start.sharedCircleWith(end); + if (circle == null) { + circle = new Circle(start.getLocation(), end.getLocation(), hyperplaneThickness); + } + + // create the edge and store it + edges.add(new Edge(start, end, + Vector3D.angle(start.getLocation().getVector(), + end.getLocation().getVector()), + circle)); + + // check if another vertex also happens to be on this circle + for (final Vertex vertex : vArray) { + if (vertex != start && vertex != end && + FastMath.abs(circle.getOffset(vertex.getLocation())) <= hyperplaneThickness) { + vertex.bindWith(circle); + } + } + + } + + // build the tree top-down + final BSPTree<Sphere2D> tree = new BSPTree<Sphere2D>(); + insertEdges(hyperplaneThickness, tree, edges); + + return tree; + + } + + /** Recursively build a tree by inserting cut sub-hyperplanes. + * @param hyperplaneThickness tolerance below which points are considered to + * belong to the hyperplane (which is therefore more a slab) + * @param node current tree node (it is a leaf node at the beginning + * of the call) + * @param edges list of edges to insert in the cell defined by this node + * (excluding edges not belonging to the cell defined by this node) + */ + private static void insertEdges(final double hyperplaneThickness, + final BSPTree<Sphere2D> node, + final List<Edge> edges) { + + // find an edge with an hyperplane that can be inserted in the node + int index = 0; + Edge inserted = null; + while (inserted == null && index < edges.size()) { + inserted = edges.get(index++); + if (!node.insertCut(inserted.getCircle())) { + inserted = null; + } + } + + if (inserted == null) { + // no suitable edge was found, the node remains a leaf node + // we need to set its inside/outside boolean indicator + final BSPTree<Sphere2D> parent = node.getParent(); + if (parent == null || node == parent.getMinus()) { + node.setAttribute(Boolean.TRUE); + } else { + node.setAttribute(Boolean.FALSE); + } + return; + } + + // we have split the node by inserting an edge as a cut sub-hyperplane + // distribute the remaining edges in the two sub-trees + final List<Edge> outsideList = new ArrayList<Edge>(); + final List<Edge> insideList = new ArrayList<Edge>(); + for (final Edge edge : edges) { + if (edge != inserted) { + edge.split(inserted.getCircle(), outsideList, insideList); + } + } + + // recurse through lower levels + if (!outsideList.isEmpty()) { + insertEdges(hyperplaneThickness, node.getPlus(), outsideList); + } else { + node.getPlus().setAttribute(Boolean.FALSE); + } + if (!insideList.isEmpty()) { + insertEdges(hyperplaneThickness, node.getMinus(), insideList); + } else { + node.getMinus().setAttribute(Boolean.TRUE); + } + + } + + /** {@inheritDoc} */ + @Override + public SphericalPolygonsSet buildNew(final BSPTree<Sphere2D> tree) { + return new SphericalPolygonsSet(tree, getTolerance()); + } + + /** {@inheritDoc} + * @exception MathIllegalStateException if the tolerance setting does not allow to build + * a clean non-ambiguous boundary + */ + @Override + protected void computeGeometricalProperties() throws MathIllegalStateException { + + final BSPTree<Sphere2D> tree = getTree(true); + + if (tree.getCut() == null) { + + // the instance has a single cell without any boundaries + + if (tree.getCut() == null && (Boolean) tree.getAttribute()) { + // the instance covers the whole space + setSize(4 * FastMath.PI); + setBarycenter(new S2Point(0, 0)); + } else { + setSize(0); + setBarycenter(S2Point.NaN); + } + + } else { + + // the instance has a boundary + final PropertiesComputer pc = new PropertiesComputer(getTolerance()); + tree.visit(pc); + setSize(pc.getArea()); + setBarycenter(pc.getBarycenter()); + + } + + } + + /** Get the boundary loops of the polygon. + * <p>The polygon boundary can be represented as a list of closed loops, + * each loop being given by exactly one of its vertices. From each loop + * start vertex, one can follow the loop by finding the outgoing edge, + * then the end vertex, then the next outgoing edge ... until the start + * vertex of the loop (exactly the same instance) is found again once + * the full loop has been visited.</p> + * <p>If the polygon has no boundary at all, a zero length loop + * array will be returned.</p> + * <p>If the polygon is a simple one-piece polygon, then the returned + * array will contain a single vertex. + * </p> + * <p>All edges in the various loops have the inside of the region on + * their left side (i.e. toward their pole) and the outside on their + * right side (i.e. away from their pole) when moving in the underlying + * circle direction. This means that the closed loops obey the direct + * trigonometric orientation.</p> + * @return boundary of the polygon, organized as an unmodifiable list of loops start vertices. + * @exception MathIllegalStateException if the tolerance setting does not allow to build + * a clean non-ambiguous boundary + * @see Vertex + * @see Edge + */ + public List<Vertex> getBoundaryLoops() throws MathIllegalStateException { + + if (loops == null) { + if (getTree(false).getCut() == null) { + loops = Collections.emptyList(); + } else { + + // sort the arcs according to their start point + final BSPTree<Sphere2D> root = getTree(true); + final EdgesBuilder visitor = new EdgesBuilder(root, getTolerance()); + root.visit(visitor); + final List<Edge> edges = visitor.getEdges(); + + + // convert the list of all edges into a list of start vertices + loops = new ArrayList<Vertex>(); + while (!edges.isEmpty()) { + + // this is an edge belonging to a new loop, store it + Edge edge = edges.get(0); + final Vertex startVertex = edge.getStart(); + loops.add(startVertex); + + // remove all remaining edges in the same loop + do { + + // remove one edge + for (final Iterator<Edge> iterator = edges.iterator(); iterator.hasNext();) { + if (iterator.next() == edge) { + iterator.remove(); + break; + } + } + + // go to next edge following the boundary loop + edge = edge.getEnd().getOutgoing(); + + } while (edge.getStart() != startVertex); + + } + + } + } + + return Collections.unmodifiableList(loops); + + } + + /** Get a spherical cap enclosing the polygon. + * <p> + * This method is intended as a first test to quickly identify points + * that are guaranteed to be outside of the region, hence performing a full + * {@link #checkPoint(org.apache.commons.math3.geometry.Vector) checkPoint} + * only if the point status remains undecided after the quick check. It is + * is therefore mostly useful to speed up computation for small polygons with + * complex shapes (say a country boundary on Earth), as the spherical cap will + * be small and hence will reliably identify a large part of the sphere as outside, + * whereas the full check can be more computing intensive. A typical use case is + * therefore: + * </p> + * <pre> + * // compute region, plus an enclosing spherical cap + * SphericalPolygonsSet complexShape = ...; + * EnclosingBall<Sphere2D, S2Point> cap = complexShape.getEnclosingCap(); + * + * // check lots of points + * for (Vector3D p : points) { + * + * final Location l; + * if (cap.contains(p)) { + * // we cannot be sure where the point is + * // we need to perform the full computation + * l = complexShape.checkPoint(v); + * } else { + * // no need to do further computation, + * // we already know the point is outside + * l = Location.OUTSIDE; + * } + * + * // use l ... + * + * } + * </pre> + * <p> + * In the special cases of empty or whole sphere polygons, special + * spherical caps are returned, with angular radius set to negative + * or positive infinity so the {@link + * EnclosingBall#contains(org.apache.commons.math3.geometry.Point) ball.contains(point)} + * method return always false or true. + * </p> + * <p> + * This method is <em>not</em> guaranteed to return the smallest enclosing cap. + * </p> + * @return a spherical cap enclosing the polygon + */ + public EnclosingBall<Sphere2D, S2Point> getEnclosingCap() { + + // handle special cases first + if (isEmpty()) { + return new EnclosingBall<Sphere2D, S2Point>(S2Point.PLUS_K, Double.NEGATIVE_INFINITY); + } + if (isFull()) { + return new EnclosingBall<Sphere2D, S2Point>(S2Point.PLUS_K, Double.POSITIVE_INFINITY); + } + + // as the polygons is neither empty nor full, it has some boundaries and cut hyperplanes + final BSPTree<Sphere2D> root = getTree(false); + if (isEmpty(root.getMinus()) && isFull(root.getPlus())) { + // the polygon covers an hemisphere, and its boundary is one 2π long edge + final Circle circle = (Circle) root.getCut().getHyperplane(); + return new EnclosingBall<Sphere2D, S2Point>(new S2Point(circle.getPole()).negate(), + 0.5 * FastMath.PI); + } + if (isFull(root.getMinus()) && isEmpty(root.getPlus())) { + // the polygon covers an hemisphere, and its boundary is one 2π long edge + final Circle circle = (Circle) root.getCut().getHyperplane(); + return new EnclosingBall<Sphere2D, S2Point>(new S2Point(circle.getPole()), + 0.5 * FastMath.PI); + } + + // gather some inside points, to be used by the encloser + final List<Vector3D> points = getInsidePoints(); + + // extract points from the boundary loops, to be used by the encloser as well + final List<Vertex> boundary = getBoundaryLoops(); + for (final Vertex loopStart : boundary) { + int count = 0; + for (Vertex v = loopStart; count == 0 || v != loopStart; v = v.getOutgoing().getEnd()) { + ++count; + points.add(v.getLocation().getVector()); + } + } + + // find the smallest enclosing 3D sphere + final SphereGenerator generator = new SphereGenerator(); + final WelzlEncloser<Euclidean3D, Vector3D> encloser = + new WelzlEncloser<Euclidean3D, Vector3D>(getTolerance(), generator); + EnclosingBall<Euclidean3D, Vector3D> enclosing3D = encloser.enclose(points); + final Vector3D[] support3D = enclosing3D.getSupport(); + + // convert to 3D sphere to spherical cap + final double r = enclosing3D.getRadius(); + final double h = enclosing3D.getCenter().getNorm(); + if (h < getTolerance()) { + // the 3D sphere is centered on the unit sphere and covers it + // fall back to a crude approximation, based only on outside convex cells + EnclosingBall<Sphere2D, S2Point> enclosingS2 = + new EnclosingBall<Sphere2D, S2Point>(S2Point.PLUS_K, Double.POSITIVE_INFINITY); + for (Vector3D outsidePoint : getOutsidePoints()) { + final S2Point outsideS2 = new S2Point(outsidePoint); + final BoundaryProjection<Sphere2D> projection = projectToBoundary(outsideS2); + if (FastMath.PI - projection.getOffset() < enclosingS2.getRadius()) { + enclosingS2 = new EnclosingBall<Sphere2D, S2Point>(outsideS2.negate(), + FastMath.PI - projection.getOffset(), + (S2Point) projection.getProjected()); + } + } + return enclosingS2; + } + final S2Point[] support = new S2Point[support3D.length]; + for (int i = 0; i < support3D.length; ++i) { + support[i] = new S2Point(support3D[i]); + } + + final EnclosingBall<Sphere2D, S2Point> enclosingS2 = + new EnclosingBall<Sphere2D, S2Point>(new S2Point(enclosing3D.getCenter()), + FastMath.acos((1 + h * h - r * r) / (2 * h)), + support); + + return enclosingS2; + + } + + /** Gather some inside points. + * @return list of points known to be strictly in all inside convex cells + */ + private List<Vector3D> getInsidePoints() { + final PropertiesComputer pc = new PropertiesComputer(getTolerance()); + getTree(true).visit(pc); + return pc.getConvexCellsInsidePoints(); + } + + /** Gather some outside points. + * @return list of points known to be strictly in all outside convex cells + */ + private List<Vector3D> getOutsidePoints() { + final SphericalPolygonsSet complement = + (SphericalPolygonsSet) new RegionFactory<Sphere2D>().getComplement(this); + final PropertiesComputer pc = new PropertiesComputer(getTolerance()); + complement.getTree(true).visit(pc); + return pc.getConvexCellsInsidePoints(); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java new file mode 100644 index 0000000..97164cc --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java @@ -0,0 +1,72 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.spherical.twod; + +import org.apache.commons.math3.geometry.euclidean.threed.Vector3D; +import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane; +import org.apache.commons.math3.geometry.partitioning.Hyperplane; +import org.apache.commons.math3.geometry.partitioning.Region; +import org.apache.commons.math3.geometry.spherical.oned.Arc; +import org.apache.commons.math3.geometry.spherical.oned.ArcsSet; +import org.apache.commons.math3.geometry.spherical.oned.Sphere1D; +import org.apache.commons.math3.util.FastMath; + +/** This class represents a sub-hyperplane for {@link Circle}. + * @since 3.3 + */ +public class SubCircle extends AbstractSubHyperplane<Sphere2D, Sphere1D> { + + /** Simple constructor. + * @param hyperplane underlying hyperplane + * @param remainingRegion remaining region of the hyperplane + */ + public SubCircle(final Hyperplane<Sphere2D> hyperplane, + final Region<Sphere1D> remainingRegion) { + super(hyperplane, remainingRegion); + } + + /** {@inheritDoc} */ + @Override + protected AbstractSubHyperplane<Sphere2D, Sphere1D> buildNew(final Hyperplane<Sphere2D> hyperplane, + final Region<Sphere1D> remainingRegion) { + return new SubCircle(hyperplane, remainingRegion); + } + + /** {@inheritDoc} */ + @Override + public SplitSubHyperplane<Sphere2D> split(final Hyperplane<Sphere2D> hyperplane) { + + final Circle thisCircle = (Circle) getHyperplane(); + final Circle otherCircle = (Circle) hyperplane; + final double angle = Vector3D.angle(thisCircle.getPole(), otherCircle.getPole()); + + if (angle < thisCircle.getTolerance() || angle > FastMath.PI - thisCircle.getTolerance()) { + // the two circles are aligned or opposite + return new SplitSubHyperplane<Sphere2D>(null, null); + } else { + // the two circles intersect each other + final Arc arc = thisCircle.getInsideArc(otherCircle); + final ArcsSet.Split split = ((ArcsSet) getRemainingRegion()).split(arc); + final ArcsSet plus = split.getPlus(); + final ArcsSet minus = split.getMinus(); + return new SplitSubHyperplane<Sphere2D>(plus == null ? null : new SubCircle(thisCircle.copySelf(), plus), + minus == null ? null : new SubCircle(thisCircle.copySelf(), minus)); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Vertex.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Vertex.java new file mode 100644 index 0000000..3003da8 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Vertex.java @@ -0,0 +1,124 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.spherical.twod; + +import java.util.ArrayList; +import java.util.List; + +/** Spherical polygons boundary vertex. + * @see SphericalPolygonsSet#getBoundaryLoops() + * @see Edge + * @since 3.3 + */ +public class Vertex { + + /** Vertex location. */ + private final S2Point location; + + /** Incoming edge. */ + private Edge incoming; + + /** Outgoing edge. */ + private Edge outgoing; + + /** Circles bound with this vertex. */ + private final List<Circle> circles; + + /** Build a non-processed vertex not owned by any node yet. + * @param location vertex location + */ + Vertex(final S2Point location) { + this.location = location; + this.incoming = null; + this.outgoing = null; + this.circles = new ArrayList<Circle>(); + } + + /** Get Vertex location. + * @return vertex location + */ + public S2Point getLocation() { + return location; + } + + /** Bind a circle considered to contain this vertex. + * @param circle circle to bind with this vertex + */ + void bindWith(final Circle circle) { + circles.add(circle); + } + + /** Get the common circle bound with both the instance and another vertex, if any. + * <p> + * When two vertices are both bound to the same circle, this means they are + * already handled by node associated with this circle, so there is no need + * to create a cut hyperplane for them. + * </p> + * @param vertex other vertex to check instance against + * @return circle bound with both the instance and another vertex, or null if the + * two vertices do not share a circle yet + */ + Circle sharedCircleWith(final Vertex vertex) { + for (final Circle circle1 : circles) { + for (final Circle circle2 : vertex.circles) { + if (circle1 == circle2) { + return circle1; + } + } + } + return null; + } + + /** Set incoming edge. + * <p> + * The circle supporting the incoming edge is automatically bound + * with the instance. + * </p> + * @param incoming incoming edge + */ + void setIncoming(final Edge incoming) { + this.incoming = incoming; + bindWith(incoming.getCircle()); + } + + /** Get incoming edge. + * @return incoming edge + */ + public Edge getIncoming() { + return incoming; + } + + /** Set outgoing edge. + * <p> + * The circle supporting the outgoing edge is automatically bound + * with the instance. + * </p> + * @param outgoing outgoing edge + */ + void setOutgoing(final Edge outgoing) { + this.outgoing = outgoing; + bindWith(outgoing.getCircle()); + } + + /** Get outgoing edge. + * @return outgoing edge + */ + public Edge getOutgoing() { + return outgoing; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/twod/package-info.java b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/package-info.java new file mode 100644 index 0000000..3f3c5b0 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/package-info.java @@ -0,0 +1,30 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** + * + * <p> + * This package provides basic geometry components on the 2-sphere. + * </p> + * <p> + * We use here the topologists definition of the 2-sphere (see + * <a href="http://mathworld.wolfram.com/Sphere.html">Sphere</a> on + * MathWorld), i.e. the 2-sphere is the two-dimensional surface + * defined in 3D as x<sup>2</sup>+y<sup>2</sup>+z<sup>2</sup>=1. + * </p> + * + */ +package org.apache.commons.math3.geometry.spherical.twod; |