diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/geometry/spherical')
17 files changed, 3571 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Arc.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Arc.java new file mode 100644 index 0000000..af0388e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Arc.java @@ -0,0 +1,132 @@ +/* + * 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.oned; + +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.geometry.partitioning.Region.Location; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; +import org.apache.commons.math3.util.Precision; + + +/** This class represents an arc on a circle. + * @see ArcsSet + * @since 3.3 + */ +public class Arc { + + /** The lower angular bound of the arc. */ + private final double lower; + + /** The upper angular bound of the arc. */ + private final double upper; + + /** Middle point of the arc. */ + private final double middle; + + /** Tolerance below which angles are considered identical. */ + private final double tolerance; + + /** Simple constructor. + * <p> + * If either {@code lower} is equals to {@code upper} or + * the interval exceeds \( 2 \pi \), the arc is considered + * to be the full circle and its initial defining boundaries + * will be forgotten. {@code lower} is not allowed to be + * greater than {@code upper} (an exception is thrown in this case). + * {@code lower} will be canonicalized between 0 and \( 2 \pi \), and + * upper shifted accordingly, so the {@link #getInf()} and {@link #getSup()} + * may not return the value used at instance construction. + * </p> + * @param lower lower angular bound of the arc + * @param upper upper angular bound of the arc + * @param tolerance tolerance below which angles are considered identical + * @exception NumberIsTooLargeException if lower is greater than upper + */ + public Arc(final double lower, final double upper, final double tolerance) + throws NumberIsTooLargeException { + this.tolerance = tolerance; + if (Precision.equals(lower, upper, 0) || (upper - lower) >= MathUtils.TWO_PI) { + // the arc must cover the whole circle + this.lower = 0; + this.upper = MathUtils.TWO_PI; + this.middle = FastMath.PI; + } else if (lower <= upper) { + this.lower = MathUtils.normalizeAngle(lower, FastMath.PI); + this.upper = this.lower + (upper - lower); + this.middle = 0.5 * (this.lower + this.upper); + } else { + throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL, + lower, upper, true); + } + } + + /** Get the lower angular bound of the arc. + * @return lower angular bound of the arc, + * always between 0 and \( 2 \pi \) + */ + public double getInf() { + return lower; + } + + /** Get the upper angular bound of the arc. + * @return upper angular bound of the arc, + * always between {@link #getInf()} and {@link #getInf()} \( + 2 \pi \) + */ + public double getSup() { + return upper; + } + + /** Get the angular size of the arc. + * @return angular size of the arc + */ + public double getSize() { + return upper - lower; + } + + /** Get the barycenter of the arc. + * @return barycenter of the arc + */ + public double getBarycenter() { + return middle; + } + + /** Get the tolerance below which angles are considered identical. + * @return tolerance below which angles are considered identical + */ + public double getTolerance() { + return tolerance; + } + + /** Check a point with respect to the arc. + * @param point point to check + * @return a code representing the point status: either {@link + * Location#INSIDE}, {@link Location#OUTSIDE} or {@link Location#BOUNDARY} + */ + public Location checkPoint(final double point) { + final double normalizedPoint = MathUtils.normalizeAngle(point, middle); + if (normalizedPoint < lower - tolerance || normalizedPoint > upper + tolerance) { + return Location.OUTSIDE; + } else if (normalizedPoint > lower + tolerance && normalizedPoint < upper - tolerance) { + return Location.INSIDE; + } else { + return (getSize() >= MathUtils.TWO_PI - tolerance) ? Location.INSIDE : Location.BOUNDARY; + } + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java new file mode 100644 index 0000000..0a00aa7 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java @@ -0,0 +1,955 @@ +/* + * 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.oned; + +import java.util.ArrayList; +import java.util.Collection; +import java.util.Iterator; +import java.util.List; +import java.util.NoSuchElementException; + +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.geometry.Point; +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.Side; +import org.apache.commons.math3.geometry.partitioning.SubHyperplane; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; +import org.apache.commons.math3.util.Precision; + +/** This class represents a region of a circle: a set of arcs. + * <p> + * Note that due to the wrapping around \(2 \pi\), barycenter is + * ill-defined here. It was defined only in order to fulfill + * the requirements of the {@link + * org.apache.commons.math3.geometry.partitioning.Region Region} + * interface, but its use is discouraged. + * </p> + * @since 3.3 + */ +public class ArcsSet extends AbstractRegion<Sphere1D, Sphere1D> implements Iterable<double[]> { + + /** Build an arcs set representing the whole circle. + * @param tolerance tolerance below which close sub-arcs are merged together + */ + public ArcsSet(final double tolerance) { + super(tolerance); + } + + /** Build an arcs set corresponding to a single arc. + * <p> + * If either {@code lower} is equals to {@code upper} or + * the interval exceeds \( 2 \pi \), the arc is considered + * to be the full circle and its initial defining boundaries + * will be forgotten. {@code lower} is not allowed to be greater + * than {@code upper} (an exception is thrown in this case). + * </p> + * @param lower lower bound of the arc + * @param upper upper bound of the arc + * @param tolerance tolerance below which close sub-arcs are merged together + * @exception NumberIsTooLargeException if lower is greater than upper + */ + public ArcsSet(final double lower, final double upper, final double tolerance) + throws NumberIsTooLargeException { + super(buildTree(lower, upper, tolerance), tolerance); + } + + /** Build an arcs set from an inside/outside 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 arcs set + * @param tolerance tolerance below which close sub-arcs are merged together + * @exception InconsistentStateAt2PiWrapping if the tree leaf nodes are not + * consistent across the \( 0, 2 \pi \) crossing + */ + public ArcsSet(final BSPTree<Sphere1D> tree, final double tolerance) + throws InconsistentStateAt2PiWrapping { + super(tree, tolerance); + check2PiConsistency(); + } + + /** Build an arcs 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 disjoints polyhedrons 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 + * @param tolerance tolerance below which close sub-arcs are merged together + * @exception InconsistentStateAt2PiWrapping if the tree leaf nodes are not + * consistent across the \( 0, 2 \pi \) crossing + */ + public ArcsSet(final Collection<SubHyperplane<Sphere1D>> boundary, final double tolerance) + throws InconsistentStateAt2PiWrapping { + super(boundary, tolerance); + check2PiConsistency(); + } + + /** Build an inside/outside tree representing a single arc. + * @param lower lower angular bound of the arc + * @param upper upper angular bound of the arc + * @param tolerance tolerance below which close sub-arcs are merged together + * @return the built tree + * @exception NumberIsTooLargeException if lower is greater than upper + */ + private static BSPTree<Sphere1D> buildTree(final double lower, final double upper, + final double tolerance) + throws NumberIsTooLargeException { + + if (Precision.equals(lower, upper, 0) || (upper - lower) >= MathUtils.TWO_PI) { + // the tree must cover the whole circle + return new BSPTree<Sphere1D>(Boolean.TRUE); + } else if (lower > upper) { + throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL, + lower, upper, true); + } + + // this is a regular arc, covering only part of the circle + final double normalizedLower = MathUtils.normalizeAngle(lower, FastMath.PI); + final double normalizedUpper = normalizedLower + (upper - lower); + final SubHyperplane<Sphere1D> lowerCut = + new LimitAngle(new S1Point(normalizedLower), false, tolerance).wholeHyperplane(); + + if (normalizedUpper <= MathUtils.TWO_PI) { + // simple arc starting after 0 and ending before 2 \pi + final SubHyperplane<Sphere1D> upperCut = + new LimitAngle(new S1Point(normalizedUpper), true, tolerance).wholeHyperplane(); + return new BSPTree<Sphere1D>(lowerCut, + new BSPTree<Sphere1D>(Boolean.FALSE), + new BSPTree<Sphere1D>(upperCut, + new BSPTree<Sphere1D>(Boolean.FALSE), + new BSPTree<Sphere1D>(Boolean.TRUE), + null), + null); + } else { + // arc wrapping around 2 \pi + final SubHyperplane<Sphere1D> upperCut = + new LimitAngle(new S1Point(normalizedUpper - MathUtils.TWO_PI), true, tolerance).wholeHyperplane(); + return new BSPTree<Sphere1D>(lowerCut, + new BSPTree<Sphere1D>(upperCut, + new BSPTree<Sphere1D>(Boolean.FALSE), + new BSPTree<Sphere1D>(Boolean.TRUE), + null), + new BSPTree<Sphere1D>(Boolean.TRUE), + null); + } + + } + + /** Check consistency. + * @exception InconsistentStateAt2PiWrapping if the tree leaf nodes are not + * consistent across the \( 0, 2 \pi \) crossing + */ + private void check2PiConsistency() throws InconsistentStateAt2PiWrapping { + + // start search at the tree root + BSPTree<Sphere1D> root = getTree(false); + if (root.getCut() == null) { + return; + } + + // find the inside/outside state before the smallest internal node + final Boolean stateBefore = (Boolean) getFirstLeaf(root).getAttribute(); + + // find the inside/outside state after the largest internal node + final Boolean stateAfter = (Boolean) getLastLeaf(root).getAttribute(); + + if (stateBefore ^ stateAfter) { + throw new InconsistentStateAt2PiWrapping(); + } + + } + + /** Get the first leaf node of a tree. + * @param root tree root + * @return first leaf node (i.e. node corresponding to the region just after 0.0 radians) + */ + private BSPTree<Sphere1D> getFirstLeaf(final BSPTree<Sphere1D> root) { + + if (root.getCut() == null) { + return root; + } + + // find the smallest internal node + BSPTree<Sphere1D> smallest = null; + for (BSPTree<Sphere1D> n = root; n != null; n = previousInternalNode(n)) { + smallest = n; + } + + return leafBefore(smallest); + + } + + /** Get the last leaf node of a tree. + * @param root tree root + * @return last leaf node (i.e. node corresponding to the region just before \( 2 \pi \) radians) + */ + private BSPTree<Sphere1D> getLastLeaf(final BSPTree<Sphere1D> root) { + + if (root.getCut() == null) { + return root; + } + + // find the largest internal node + BSPTree<Sphere1D> largest = null; + for (BSPTree<Sphere1D> n = root; n != null; n = nextInternalNode(n)) { + largest = n; + } + + return leafAfter(largest); + + } + + /** Get the node corresponding to the first arc start. + * @return smallest internal node (i.e. first after 0.0 radians, in trigonometric direction), + * or null if there are no internal nodes (i.e. the set is either empty or covers the full circle) + */ + private BSPTree<Sphere1D> getFirstArcStart() { + + // start search at the tree root + BSPTree<Sphere1D> node = getTree(false); + if (node.getCut() == null) { + return null; + } + + // walk tree until we find the smallest internal node + node = getFirstLeaf(node).getParent(); + + // walk tree until we find an arc start + while (node != null && !isArcStart(node)) { + node = nextInternalNode(node); + } + + return node; + + } + + /** Check if an internal node corresponds to the start angle of an arc. + * @param node internal node to check + * @return true if the node corresponds to the start angle of an arc + */ + private boolean isArcStart(final BSPTree<Sphere1D> node) { + + if ((Boolean) leafBefore(node).getAttribute()) { + // it has an inside cell before it, it may end an arc but not start it + return false; + } + + if (!(Boolean) leafAfter(node).getAttribute()) { + // it has an outside cell after it, it is a dummy cut away from real arcs + return false; + } + + // the cell has an outside before and an inside after it + // it is the start of an arc + return true; + + } + + /** Check if an internal node corresponds to the end angle of an arc. + * @param node internal node to check + * @return true if the node corresponds to the end angle of an arc + */ + private boolean isArcEnd(final BSPTree<Sphere1D> node) { + + if (!(Boolean) leafBefore(node).getAttribute()) { + // it has an outside cell before it, it may start an arc but not end it + return false; + } + + if ((Boolean) leafAfter(node).getAttribute()) { + // it has an inside cell after it, it is a dummy cut in the middle of an arc + return false; + } + + // the cell has an inside before and an outside after it + // it is the end of an arc + return true; + + } + + /** Get the next internal node. + * @param node current internal node + * @return next internal node in trigonometric order, or null + * if this is the last internal node + */ + private BSPTree<Sphere1D> nextInternalNode(BSPTree<Sphere1D> node) { + + if (childAfter(node).getCut() != null) { + // the next node is in the sub-tree + return leafAfter(node).getParent(); + } + + // there is nothing left deeper in the tree, we backtrack + while (isAfterParent(node)) { + node = node.getParent(); + } + return node.getParent(); + + } + + /** Get the previous internal node. + * @param node current internal node + * @return previous internal node in trigonometric order, or null + * if this is the first internal node + */ + private BSPTree<Sphere1D> previousInternalNode(BSPTree<Sphere1D> node) { + + if (childBefore(node).getCut() != null) { + // the next node is in the sub-tree + return leafBefore(node).getParent(); + } + + // there is nothing left deeper in the tree, we backtrack + while (isBeforeParent(node)) { + node = node.getParent(); + } + return node.getParent(); + + } + + /** Find the leaf node just before an internal node. + * @param node internal node at which the sub-tree starts + * @return leaf node just before the internal node + */ + private BSPTree<Sphere1D> leafBefore(BSPTree<Sphere1D> node) { + + node = childBefore(node); + while (node.getCut() != null) { + node = childAfter(node); + } + + return node; + + } + + /** Find the leaf node just after an internal node. + * @param node internal node at which the sub-tree starts + * @return leaf node just after the internal node + */ + private BSPTree<Sphere1D> leafAfter(BSPTree<Sphere1D> node) { + + node = childAfter(node); + while (node.getCut() != null) { + node = childBefore(node); + } + + return node; + + } + + /** Check if a node is the child before its parent in trigonometric order. + * @param node child node considered + * @return true is the node has a parent end is before it in trigonometric order + */ + private boolean isBeforeParent(final BSPTree<Sphere1D> node) { + final BSPTree<Sphere1D> parent = node.getParent(); + if (parent == null) { + return false; + } else { + return node == childBefore(parent); + } + } + + /** Check if a node is the child after its parent in trigonometric order. + * @param node child node considered + * @return true is the node has a parent end is after it in trigonometric order + */ + private boolean isAfterParent(final BSPTree<Sphere1D> node) { + final BSPTree<Sphere1D> parent = node.getParent(); + if (parent == null) { + return false; + } else { + return node == childAfter(parent); + } + } + + /** Find the child node just before an internal node. + * @param node internal node at which the sub-tree starts + * @return child node just before the internal node + */ + private BSPTree<Sphere1D> childBefore(BSPTree<Sphere1D> node) { + if (isDirect(node)) { + // smaller angles are on minus side, larger angles are on plus side + return node.getMinus(); + } else { + // smaller angles are on plus side, larger angles are on minus side + return node.getPlus(); + } + } + + /** Find the child node just after an internal node. + * @param node internal node at which the sub-tree starts + * @return child node just after the internal node + */ + private BSPTree<Sphere1D> childAfter(BSPTree<Sphere1D> node) { + if (isDirect(node)) { + // smaller angles are on minus side, larger angles are on plus side + return node.getPlus(); + } else { + // smaller angles are on plus side, larger angles are on minus side + return node.getMinus(); + } + } + + /** Check if an internal node has a direct limit angle. + * @param node internal node to check + * @return true if the limit angle is direct + */ + private boolean isDirect(final BSPTree<Sphere1D> node) { + return ((LimitAngle) node.getCut().getHyperplane()).isDirect(); + } + + /** Get the limit angle of an internal node. + * @param node internal node to check + * @return limit angle + */ + private double getAngle(final BSPTree<Sphere1D> node) { + return ((LimitAngle) node.getCut().getHyperplane()).getLocation().getAlpha(); + } + + /** {@inheritDoc} */ + @Override + public ArcsSet buildNew(final BSPTree<Sphere1D> tree) { + return new ArcsSet(tree, getTolerance()); + } + + /** {@inheritDoc} */ + @Override + protected void computeGeometricalProperties() { + if (getTree(false).getCut() == null) { + setBarycenter(S1Point.NaN); + setSize(((Boolean) getTree(false).getAttribute()) ? MathUtils.TWO_PI : 0); + } else { + double size = 0.0; + double sum = 0.0; + for (final double[] a : this) { + final double length = a[1] - a[0]; + size += length; + sum += length * (a[0] + a[1]); + } + setSize(size); + if (Precision.equals(size, MathUtils.TWO_PI, 0)) { + setBarycenter(S1Point.NaN); + } else if (size >= Precision.SAFE_MIN) { + setBarycenter(new S1Point(sum / (2 * size))); + } else { + final LimitAngle limit = (LimitAngle) getTree(false).getCut().getHyperplane(); + setBarycenter(limit.getLocation()); + } + } + } + + /** {@inheritDoc} + * @since 3.3 + */ + @Override + public BoundaryProjection<Sphere1D> projectToBoundary(final Point<Sphere1D> point) { + + // get position of test point + final double alpha = ((S1Point) point).getAlpha(); + + boolean wrapFirst = false; + double first = Double.NaN; + double previous = Double.NaN; + for (final double[] a : this) { + + if (Double.isNaN(first)) { + // remember the first angle in case we need it later + first = a[0]; + } + + if (!wrapFirst) { + if (alpha < a[0]) { + // the test point lies between the previous and the current arcs + // offset will be positive + if (Double.isNaN(previous)) { + // we need to wrap around the circle + wrapFirst = true; + } else { + final double previousOffset = alpha - previous; + final double currentOffset = a[0] - alpha; + if (previousOffset < currentOffset) { + return new BoundaryProjection<Sphere1D>(point, new S1Point(previous), previousOffset); + } else { + return new BoundaryProjection<Sphere1D>(point, new S1Point(a[0]), currentOffset); + } + } + } else if (alpha <= a[1]) { + // the test point lies within the current arc + // offset will be negative + final double offset0 = a[0] - alpha; + final double offset1 = alpha - a[1]; + if (offset0 < offset1) { + return new BoundaryProjection<Sphere1D>(point, new S1Point(a[1]), offset1); + } else { + return new BoundaryProjection<Sphere1D>(point, new S1Point(a[0]), offset0); + } + } + } + previous = a[1]; + } + + if (Double.isNaN(previous)) { + + // there are no points at all in the arcs set + return new BoundaryProjection<Sphere1D>(point, null, MathUtils.TWO_PI); + + } else { + + // the test point if before first arc and after last arc, + // somewhere around the 0/2 \pi crossing + if (wrapFirst) { + // the test point is between 0 and first + final double previousOffset = alpha - (previous - MathUtils.TWO_PI); + final double currentOffset = first - alpha; + if (previousOffset < currentOffset) { + return new BoundaryProjection<Sphere1D>(point, new S1Point(previous), previousOffset); + } else { + return new BoundaryProjection<Sphere1D>(point, new S1Point(first), currentOffset); + } + } else { + // the test point is between last and 2\pi + final double previousOffset = alpha - previous; + final double currentOffset = first + MathUtils.TWO_PI - alpha; + if (previousOffset < currentOffset) { + return new BoundaryProjection<Sphere1D>(point, new S1Point(previous), previousOffset); + } else { + return new BoundaryProjection<Sphere1D>(point, new S1Point(first), currentOffset); + } + } + + } + + } + + /** Build an ordered list of arcs representing the instance. + * <p>This method builds this arcs set as an ordered list of + * {@link Arc Arc} elements. An empty tree will build an empty list + * while a tree representing the whole circle will build a one + * element list with bounds set to \( 0 and 2 \pi \).</p> + * @return a new ordered list containing {@link Arc Arc} elements + */ + public List<Arc> asList() { + final List<Arc> list = new ArrayList<Arc>(); + for (final double[] a : this) { + list.add(new Arc(a[0], a[1], getTolerance())); + } + return list; + } + + /** {@inheritDoc} + * <p> + * The iterator returns the limit angles pairs of sub-arcs in trigonometric order. + * </p> + * <p> + * The iterator does <em>not</em> support the optional {@code remove} operation. + * </p> + */ + public Iterator<double[]> iterator() { + return new SubArcsIterator(); + } + + /** Local iterator for sub-arcs. */ + private class SubArcsIterator implements Iterator<double[]> { + + /** Start of the first arc. */ + private final BSPTree<Sphere1D> firstStart; + + /** Current node. */ + private BSPTree<Sphere1D> current; + + /** Sub-arc no yet returned. */ + private double[] pending; + + /** Simple constructor. + */ + SubArcsIterator() { + + firstStart = getFirstArcStart(); + current = firstStart; + + if (firstStart == null) { + // all the leaf tree nodes share the same inside/outside status + if ((Boolean) getFirstLeaf(getTree(false)).getAttribute()) { + // it is an inside node, it represents the full circle + pending = new double[] { + 0, MathUtils.TWO_PI + }; + } else { + pending = null; + } + } else { + selectPending(); + } + } + + /** Walk the tree to select the pending sub-arc. + */ + private void selectPending() { + + // look for the start of the arc + BSPTree<Sphere1D> start = current; + while (start != null && !isArcStart(start)) { + start = nextInternalNode(start); + } + + if (start == null) { + // we have exhausted the iterator + current = null; + pending = null; + return; + } + + // look for the end of the arc + BSPTree<Sphere1D> end = start; + while (end != null && !isArcEnd(end)) { + end = nextInternalNode(end); + } + + if (end != null) { + + // we have identified the arc + pending = new double[] { + getAngle(start), getAngle(end) + }; + + // prepare search for next arc + current = end; + + } else { + + // the final arc wraps around 2\pi, its end is before the first start + end = firstStart; + while (end != null && !isArcEnd(end)) { + end = previousInternalNode(end); + } + if (end == null) { + // this should never happen + throw new MathInternalError(); + } + + // we have identified the last arc + pending = new double[] { + getAngle(start), getAngle(end) + MathUtils.TWO_PI + }; + + // there won't be any other arcs + current = null; + + } + + } + + /** {@inheritDoc} */ + public boolean hasNext() { + return pending != null; + } + + /** {@inheritDoc} */ + public double[] next() { + if (pending == null) { + throw new NoSuchElementException(); + } + final double[] next = pending; + selectPending(); + return next; + } + + /** {@inheritDoc} */ + public void remove() { + throw new UnsupportedOperationException(); + } + + } + + /** Compute the relative position of the instance with respect + * to an arc. + * <p> + * The {@link Side#MINUS} side of the arc is the one covered by the arc. + * </p> + * @param arc arc to check instance against + * @return one of {@link Side#PLUS}, {@link Side#MINUS}, {@link Side#BOTH} + * or {@link Side#HYPER} + * @deprecated as of 3.6, replaced with {@link #split(Arc)}.{@link Split#getSide()} + */ + @Deprecated + public Side side(final Arc arc) { + return split(arc).getSide(); + } + + /** Split the instance in two parts by an arc. + * @param arc splitting arc + * @return an object containing both the part of the instance + * on the plus side of the arc and the part of the + * instance on the minus side of the arc + */ + public Split split(final Arc arc) { + + final List<Double> minus = new ArrayList<Double>(); + final List<Double> plus = new ArrayList<Double>(); + + final double reference = FastMath.PI + arc.getInf(); + final double arcLength = arc.getSup() - arc.getInf(); + + for (final double[] a : this) { + final double syncedStart = MathUtils.normalizeAngle(a[0], reference) - arc.getInf(); + final double arcOffset = a[0] - syncedStart; + final double syncedEnd = a[1] - arcOffset; + if (syncedStart < arcLength) { + // the start point a[0] is in the minus part of the arc + minus.add(a[0]); + if (syncedEnd > arcLength) { + // the end point a[1] is past the end of the arc + // so we leave the minus part and enter the plus part + final double minusToPlus = arcLength + arcOffset; + minus.add(minusToPlus); + plus.add(minusToPlus); + if (syncedEnd > MathUtils.TWO_PI) { + // in fact the end point a[1] goes far enough that we + // leave the plus part of the arc and enter the minus part again + final double plusToMinus = MathUtils.TWO_PI + arcOffset; + plus.add(plusToMinus); + minus.add(plusToMinus); + minus.add(a[1]); + } else { + // the end point a[1] is in the plus part of the arc + plus.add(a[1]); + } + } else { + // the end point a[1] is in the minus part of the arc + minus.add(a[1]); + } + } else { + // the start point a[0] is in the plus part of the arc + plus.add(a[0]); + if (syncedEnd > MathUtils.TWO_PI) { + // the end point a[1] wraps around to the start of the arc + // so we leave the plus part and enter the minus part + final double plusToMinus = MathUtils.TWO_PI + arcOffset; + plus.add(plusToMinus); + minus.add(plusToMinus); + if (syncedEnd > MathUtils.TWO_PI + arcLength) { + // in fact the end point a[1] goes far enough that we + // leave the minus part of the arc and enter the plus part again + final double minusToPlus = MathUtils.TWO_PI + arcLength + arcOffset; + minus.add(minusToPlus); + plus.add(minusToPlus); + plus.add(a[1]); + } else { + // the end point a[1] is in the minus part of the arc + minus.add(a[1]); + } + } else { + // the end point a[1] is in the plus part of the arc + plus.add(a[1]); + } + } + } + + return new Split(createSplitPart(plus), createSplitPart(minus)); + + } + + /** Add an arc limit to a BSP tree under construction. + * @param tree BSP tree under construction + * @param alpha arc limit + * @param isStart if true, the limit is the start of an arc + */ + private void addArcLimit(final BSPTree<Sphere1D> tree, final double alpha, final boolean isStart) { + + final LimitAngle limit = new LimitAngle(new S1Point(alpha), !isStart, getTolerance()); + final BSPTree<Sphere1D> node = tree.getCell(limit.getLocation(), getTolerance()); + if (node.getCut() != null) { + // this should never happen + throw new MathInternalError(); + } + + node.insertCut(limit); + node.setAttribute(null); + node.getPlus().setAttribute(Boolean.FALSE); + node.getMinus().setAttribute(Boolean.TRUE); + + } + + /** Create a split part. + * <p> + * As per construction, the list of limit angles is known to have + * an even number of entries, with start angles at even indices and + * end angles at odd indices. + * </p> + * @param limits limit angles of the split part + * @return split part (may be null) + */ + private ArcsSet createSplitPart(final List<Double> limits) { + if (limits.isEmpty()) { + return null; + } else { + + // collapse close limit angles + for (int i = 0; i < limits.size(); ++i) { + final int j = (i + 1) % limits.size(); + final double lA = limits.get(i); + final double lB = MathUtils.normalizeAngle(limits.get(j), lA); + if (FastMath.abs(lB - lA) <= getTolerance()) { + // the two limits are too close to each other, we remove both of them + if (j > 0) { + // regular case, the two entries are consecutive ones + limits.remove(j); + limits.remove(i); + i = i - 1; + } else { + // special case, i the the last entry and j is the first entry + // we have wrapped around list end + final double lEnd = limits.remove(limits.size() - 1); + final double lStart = limits.remove(0); + if (limits.isEmpty()) { + // the ends were the only limits, is it a full circle or an empty circle? + if (lEnd - lStart > FastMath.PI) { + // it was full circle + return new ArcsSet(new BSPTree<Sphere1D>(Boolean.TRUE), getTolerance()); + } else { + // it was an empty circle + return null; + } + } else { + // we have removed the first interval start, so our list + // currently starts with an interval end, which is wrong + // we need to move this interval end to the end of the list + limits.add(limits.remove(0) + MathUtils.TWO_PI); + } + } + } + } + + // build the tree by adding all angular sectors + BSPTree<Sphere1D> tree = new BSPTree<Sphere1D>(Boolean.FALSE); + for (int i = 0; i < limits.size() - 1; i += 2) { + addArcLimit(tree, limits.get(i), true); + addArcLimit(tree, limits.get(i + 1), false); + } + + if (tree.getCut() == null) { + // we did not insert anything + return null; + } + + return new ArcsSet(tree, getTolerance()); + + } + } + + /** Class holding the results of the {@link #split split} method. + */ + public static class Split { + + /** Part of the arcs set on the plus side of the splitting arc. */ + private final ArcsSet plus; + + /** Part of the arcs set on the minus side of the splitting arc. */ + private final ArcsSet minus; + + /** Build a Split from its parts. + * @param plus part of the arcs set on the plus side of the + * splitting arc + * @param minus part of the arcs set on the minus side of the + * splitting arc + */ + private Split(final ArcsSet plus, final ArcsSet minus) { + this.plus = plus; + this.minus = minus; + } + + /** Get the part of the arcs set on the plus side of the splitting arc. + * @return part of the arcs set on the plus side of the splitting arc + */ + public ArcsSet getPlus() { + return plus; + } + + /** Get the part of the arcs set on the minus side of the splitting arc. + * @return part of the arcs set on the minus side of the splitting arc + */ + public ArcsSet getMinus() { + return minus; + } + + /** Get the side of the split arc with respect to its splitter. + * @return {@link Side#PLUS} if only {@link #getPlus()} returns non-null, + * {@link Side#MINUS} if only {@link #getMinus()} returns non-null, + * {@link Side#BOTH} if both {@link #getPlus()} and {@link #getMinus()} + * return non-null or {@link Side#HYPER} if both {@link #getPlus()} and + * {@link #getMinus()} return null + * @since 3.6 + */ + public Side getSide() { + if (plus != null) { + if (minus != null) { + return Side.BOTH; + } else { + return Side.PLUS; + } + } else if (minus != null) { + return Side.MINUS; + } else { + return Side.HYPER; + } + } + + } + + /** Specialized exception for inconsistent BSP tree state inconsistency. + * <p> + * This exception is thrown at {@link ArcsSet} construction time when the + * {@link org.apache.commons.math3.geometry.partitioning.Region.Location inside/outside} + * state is not consistent at the 0, \(2 \pi \) crossing. + * </p> + */ + public static class InconsistentStateAt2PiWrapping extends MathIllegalArgumentException { + + /** Serializable UID. */ + private static final long serialVersionUID = 20140107L; + + /** Simple constructor. + */ + public InconsistentStateAt2PiWrapping() { + super(LocalizedFormats.INCONSISTENT_STATE_AT_2_PI_WRAPPING); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/LimitAngle.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/LimitAngle.java new file mode 100644 index 0000000..748a142 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/LimitAngle.java @@ -0,0 +1,127 @@ +/* + * 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.oned; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.partitioning.Hyperplane; + +/** This class represents a 1D oriented hyperplane on the circle. + * <p>An hyperplane on the 1-sphere is an angle with an orientation.</p> + * <p>Instances of this class are guaranteed to be immutable.</p> + * @since 3.3 + */ +public class LimitAngle implements Hyperplane<Sphere1D> { + + /** Angle location. */ + private S1Point location; + + /** Orientation. */ + private boolean direct; + + /** Tolerance below which angles are considered identical. */ + private final double tolerance; + + /** Simple constructor. + * @param location location of the hyperplane + * @param direct if true, the plus side of the hyperplane is towards + * angles greater than {@code location} + * @param tolerance tolerance below which angles are considered identical + */ + public LimitAngle(final S1Point location, final boolean direct, final double tolerance) { + this.location = location; + this.direct = direct; + this.tolerance = tolerance; + } + + /** Copy the instance. + * <p>Since instances are immutable, this method directly returns + * the instance.</p> + * @return the instance itself + */ + public LimitAngle copySelf() { + return this; + } + + /** {@inheritDoc} */ + public double getOffset(final Point<Sphere1D> point) { + final double delta = ((S1Point) point).getAlpha() - location.getAlpha(); + return direct ? delta : -delta; + } + + /** Check if the hyperplane orientation is direct. + * @return true if the plus side of the hyperplane is towards + * angles greater than hyperplane location + */ + public boolean isDirect() { + return direct; + } + + /** Get the reverse of the instance. + * <p>Get a limit angle with reversed orientation with respect to the + * instance. A new object is built, the instance is untouched.</p> + * @return a new limit angle, with orientation opposite to the instance orientation + */ + public LimitAngle getReverse() { + return new LimitAngle(location, !direct, tolerance); + } + + /** Build a region covering the whole hyperplane. + * <p>Since this class represent zero dimension spaces which does + * not have lower dimension sub-spaces, this method returns a dummy + * implementation of a {@link + * org.apache.commons.math3.geometry.partitioning.SubHyperplane SubHyperplane}. + * This implementation is only used to allow the {@link + * org.apache.commons.math3.geometry.partitioning.SubHyperplane + * SubHyperplane} class implementation to work properly, it should + * <em>not</em> be used otherwise.</p> + * @return a dummy sub hyperplane + */ + public SubLimitAngle wholeHyperplane() { + return new SubLimitAngle(this, null); + } + + /** Build a region covering the whole space. + * @return a region containing the instance (really an {@link + * ArcsSet IntervalsSet} instance) + */ + public ArcsSet wholeSpace() { + return new ArcsSet(tolerance); + } + + /** {@inheritDoc} */ + public boolean sameOrientationAs(final Hyperplane<Sphere1D> other) { + return !(direct ^ ((LimitAngle) other).direct); + } + + /** Get the hyperplane location on the circle. + * @return the hyperplane location + */ + public S1Point getLocation() { + return location; + } + + /** {@inheritDoc} */ + public Point<Sphere1D> project(Point<Sphere1D> point) { + return location; + } + + /** {@inheritDoc} */ + public double getTolerance() { + return tolerance; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/S1Point.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/S1Point.java new file mode 100644 index 0000000..263a559 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/S1Point.java @@ -0,0 +1,157 @@ +/* + * 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.oned; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; + +/** This class represents a point on the 1-sphere. + * <p>Instances of this class are guaranteed to be immutable.</p> + * @since 3.3 + */ +public class S1Point implements Point<Sphere1D> { + + // CHECKSTYLE: stop ConstantName + /** A vector with all coordinates set to NaN. */ + public static final S1Point NaN = new S1Point(Double.NaN, Vector2D.NaN); + // CHECKSTYLE: resume ConstantName + + /** Serializable UID. */ + private static final long serialVersionUID = 20131218L; + + /** Azimuthal angle \( \alpha \). */ + private final double alpha; + + /** Corresponding 2D normalized vector. */ + private final Vector2D vector; + + /** Simple constructor. + * Build a vector from its coordinates + * @param alpha azimuthal angle \( \alpha \) + * @see #getAlpha() + */ + public S1Point(final double alpha) { + this(MathUtils.normalizeAngle(alpha, FastMath.PI), + new Vector2D(FastMath.cos(alpha), FastMath.sin(alpha))); + } + + /** Build a point from its internal components. + * @param alpha azimuthal angle \( \alpha \) + * @param vector corresponding vector + */ + private S1Point(final double alpha, final Vector2D vector) { + this.alpha = alpha; + this.vector = vector; + } + + /** Get the azimuthal angle \( \alpha \). + * @return azimuthal angle \( \alpha \) + * @see #S1Point(double) + */ + public double getAlpha() { + return alpha; + } + + /** Get the corresponding normalized vector in the 2D euclidean space. + * @return normalized vector + */ + public Vector2D getVector() { + return vector; + } + + /** {@inheritDoc} */ + public Space getSpace() { + return Sphere1D.getInstance(); + } + + /** {@inheritDoc} */ + public boolean isNaN() { + return Double.isNaN(alpha); + } + + /** {@inheritDoc} */ + public double distance(final Point<Sphere1D> point) { + return distance(this, (S1Point) 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(S1Point p1, S1Point p2) { + return Vector2D.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 S1Point) { + final S1Point rhs = (S1Point) other; + if (rhs.isNaN()) { + return this.isNaN(); + } + + return alpha == rhs.alpha; + } + + 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 1759 * MathUtils.hash(alpha); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Sphere1D.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Sphere1D.java new file mode 100644 index 0000000..ce5c7cd --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Sphere1D.java @@ -0,0 +1,106 @@ +/* + * 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.oned; + +import java.io.Serializable; + +import org.apache.commons.math3.exception.MathUnsupportedOperationException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.geometry.Space; + +/** + * This class implements a one-dimensional sphere (i.e. a circle). + * <p> + * We use here the topologists definition of the 1-sphere (see + * <a href="http://mathworld.wolfram.com/Sphere.html">Sphere</a> on + * MathWorld), i.e. the 1-sphere is the one-dimensional closed curve + * defined in 2D as x<sup>2</sup>+y<sup>2</sup>=1. + * </p> + * @since 3.3 + */ +public class Sphere1D implements Serializable, Space { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 20131218L; + + /** Private constructor for the singleton. + */ + private Sphere1D() { + } + + /** Get the unique instance. + * @return the unique instance + */ + public static Sphere1D getInstance() { + return LazyHolder.INSTANCE; + } + + /** {@inheritDoc} */ + public int getDimension() { + return 1; + } + + /** {@inheritDoc} + * <p> + * As the 1-dimension sphere does not have proper sub-spaces, + * this method always throws a {@link NoSubSpaceException} + * </p> + * @return nothing + * @throws NoSubSpaceException in all cases + */ + public Space getSubSpace() throws NoSubSpaceException { + throw new NoSubSpaceException(); + } + + // 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 Sphere1D INSTANCE = new Sphere1D(); + } + // CHECKSTYLE: resume HideUtilityClassConstructor + + /** Handle deserialization of the singleton. + * @return the singleton instance + */ + private Object readResolve() { + // return the singleton instance + return LazyHolder.INSTANCE; + } + + /** Specialized exception for inexistent sub-space. + * <p> + * This exception is thrown when attempting to get the sub-space of a one-dimensional space + * </p> + */ + public static class NoSubSpaceException extends MathUnsupportedOperationException { + + /** Serializable UID. */ + private static final long serialVersionUID = 20140225L; + + /** Simple constructor. + */ + public NoSubSpaceException() { + super(LocalizedFormats.NOT_SUPPORTED_IN_DIMENSION_N, 1); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubLimitAngle.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubLimitAngle.java new file mode 100644 index 0000000..ebd3627 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubLimitAngle.java @@ -0,0 +1,66 @@ +/* + * 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.oned; + +import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane; +import org.apache.commons.math3.geometry.partitioning.Hyperplane; +import org.apache.commons.math3.geometry.partitioning.Region; + +/** This class represents sub-hyperplane for {@link LimitAngle}. + * <p>Instances of this class are guaranteed to be immutable.</p> + * @since 3.3 + */ +public class SubLimitAngle extends AbstractSubHyperplane<Sphere1D, Sphere1D> { + + /** Simple constructor. + * @param hyperplane underlying hyperplane + * @param remainingRegion remaining region of the hyperplane + */ + public SubLimitAngle(final Hyperplane<Sphere1D> hyperplane, + final Region<Sphere1D> remainingRegion) { + super(hyperplane, remainingRegion); + } + + /** {@inheritDoc} */ + @Override + public double getSize() { + return 0; + } + + /** {@inheritDoc} */ + @Override + public boolean isEmpty() { + return false; + } + + /** {@inheritDoc} */ + @Override + protected AbstractSubHyperplane<Sphere1D, Sphere1D> buildNew(final Hyperplane<Sphere1D> hyperplane, + final Region<Sphere1D> remainingRegion) { + return new SubLimitAngle(hyperplane, remainingRegion); + } + + /** {@inheritDoc} */ + @Override + public SplitSubHyperplane<Sphere1D> split(final Hyperplane<Sphere1D> hyperplane) { + final double global = hyperplane.getOffset(((LimitAngle) getHyperplane()).getLocation()); + return (global < -1.0e-10) ? + new SplitSubHyperplane<Sphere1D>(null, this) : + new SplitSubHyperplane<Sphere1D>(this, null); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/package-info.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/package-info.java new file mode 100644 index 0000000..d54bc0b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/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 1-sphere. + * </p> + * <p> + * We use here the topologists definition of the 1-sphere (see + * <a href="http://mathworld.wolfram.com/Sphere.html">Sphere</a> on + * MathWorld), i.e. the 1-sphere is the one-dimensional closed curve + * defined in 2D as x<sup>2</sup>+y<sup>2</sup>=1. + * </p> + * + */ +package org.apache.commons.math3.geometry.spherical.oned; 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; |