diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java index 7d13dfc48..784aa0b70 100644 --- a/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java @@ -16,6 +16,9 @@ */ package org.apache.commons.math3.geometry.partitioning; +import java.util.ArrayList; +import java.util.List; + import org.apache.commons.math3.exception.MathInternalError; import org.apache.commons.math3.geometry.Point; import org.apache.commons.math3.geometry.Space; @@ -328,6 +331,48 @@ public class BSPTree { } + /** Get the cells whose cut sub-hyperplanes are close to the point. + * @param point point to check + * @param maxOffset offset below which a cut sub-hyperplane is considered + * close to the point (in absolute value) + * @return close cells (may be empty if all cut sub-hyperplanes are farther + * than maxOffset from the point) + */ + public List> getCloseCuts(final Point point, final double maxOffset) { + final List> close = new ArrayList>(); + recurseCloseCuts(point, maxOffset, close); + return close; + } + + /** Get the cells whose cut sub-hyperplanes are close to the point. + * @param point point to check + * @param maxOffset offset below which a cut sub-hyperplane is considered + * close to the point (in absolute value) + * @param close list to fill + */ + private void recurseCloseCuts(final Point point, final double maxOffset, + final List> close) { + if (cut != null) { + + // position of the point with respect to the cut hyperplane + final double offset = cut.getHyperplane().getOffset(point); + + if (offset < -maxOffset) { + // point is on the minus side of the cut hyperplane + minus.recurseCloseCuts(point, maxOffset, close); + } else if (offset > maxOffset) { + // point is on the plus side of the cut hyperplane + plus.recurseCloseCuts(point, maxOffset, close); + } else { + // point is close to the cut hyperplane + close.add(this); + minus.recurseCloseCuts(point, maxOffset, close); + plus.recurseCloseCuts(point, maxOffset, close); + } + + } + } + /** Perform condensation on a tree. *

The condensation operation is not recursive, it must be called * explicitly from leaves to root.

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 index 3de14e626..e26b3f37d 100644 --- 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 @@ -33,10 +33,14 @@ import org.apache.commons.math3.util.Precision; */ public class ArcsSet extends AbstractRegion { + /** Tolerance below which close sub-arcs are merged together. */ + private final double tolerance; + /** Build an arcs set representing the whole circle. + * @param tolerance tolerance below which close sub-arcs are merged together */ - public ArcsSet() { - super(); + public ArcsSet(final double tolerance) { + this.tolerance = tolerance; } /** Build an arcs set corresponding to a single arc. @@ -51,9 +55,11 @@ public class ArcsSet extends AbstractRegion { *

* @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 */ - public ArcsSet(final double lower, final double upper) { - super(buildTree(lower, upper)); + public ArcsSet(final double lower, final double upper, final double tolerance) { + super(buildTree(lower, upper, tolerance)); + this.tolerance = tolerance; } /** Build an arcs set from an inside/outside BSP tree. @@ -64,9 +70,11 @@ public class ArcsSet extends AbstractRegion { * recommended to use the predefined constants * {@code Boolean.TRUE} and {@code Boolean.FALSE}

* @param tree inside/outside BSP tree representing the arcs set + * @param tolerance tolerance below which close sub-arcs are merged together */ - public ArcsSet(final BSPTree tree) { + public ArcsSet(final BSPTree tree, final double tolerance) { super(tree); + this.tolerance = tolerance; } /** Build an arcs set from a Boundary REPresentation (B-rep). @@ -87,9 +95,11 @@ public class ArcsSet extends AbstractRegion { *

If the boundary is empty, the region will represent the whole * space.

* @param boundary collection of boundary elements + * @param tolerance tolerance below which close sub-arcs are merged together */ - public ArcsSet(final Collection> boundary) { + public ArcsSet(final Collection> boundary, final double tolerance) { super(boundary); + this.tolerance = tolerance; } /** Build an inside/outside tree representing a single arc. @@ -104,9 +114,10 @@ public class ArcsSet extends AbstractRegion { *

* @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 */ - private static BSPTree buildTree(final double lower, final double upper) { + private static BSPTree buildTree(final double lower, final double upper, final double tolerance) { if (Precision.equals(lower, upper, 0)) { // the tree must cover the whole real line @@ -115,7 +126,7 @@ public class ArcsSet extends AbstractRegion { // the two boundary angles define only one cutting chord final double normalizedUpper = MathUtils.normalizeAngle(upper, lower + FastMath.PI); - return new BSPTree(new Chord(lower, normalizedUpper).wholeHyperplane(), + return new BSPTree(new Chord(lower, normalizedUpper, tolerance).wholeHyperplane(), new BSPTree(Boolean.FALSE), new BSPTree(Boolean.TRUE), null); @@ -125,7 +136,7 @@ public class ArcsSet extends AbstractRegion { /** {@inheritDoc} */ @Override public ArcsSet buildNew(final BSPTree tree) { - return new ArcsSet(tree); + return new ArcsSet(tree, tolerance); } /** {@inheritDoc} */ @@ -181,21 +192,29 @@ public class ArcsSet extends AbstractRegion { list.add(new Arc(lower, upper)); } } else { - final Chord chord = (Chord) node.getCut().getHyperplane(); - final S1Point loc = chord.getMiddle(); - double alpha = loc.getAlpha(); + final SubChord cut = (SubChord) node.getCut(); + final List cutArcs = cut.getSubArcs(); + final double cutStart = cutArcs.get(0).getInf(); + final double cutEnd = cutArcs.get(cutArcs.size() - 1).getSup(); - // make sure we explore the tree in increasing order - final BSPTree low = chord.isDirect() ? node.getMinus() : node.getPlus(); - final BSPTree high = chord.isDirect() ? node.getPlus() : node.getMinus(); - - recurseList(low, list, lower, alpha); - if ((checkPoint(low, loc) == Location.INSIDE) && - (checkPoint(high, loc) == Location.INSIDE)) { - // merge the last arc added and the first one of the high sub-tree - alpha = list.remove(list.size() - 1).getInf(); + recurseList(node.getMinus(), list, lower, cutStart); + if (list.get(list.size() - 1).checkPoint(cutStart, tolerance) == Location.INSIDE) { + // merge the last arc before the sub-chord and the sub-chord + final Arc merged = new Arc(list.remove(list.size() - 1).getInf(), + cutArcs.remove(0).getSup()); + cutArcs.add(0, merged); } - recurseList(high, list, alpha, upper); + + final List highList = new ArrayList(); + recurseList(node.getPlus(), highList, cutEnd, upper); + if (highList.get(0).checkPoint(cutEnd, tolerance) == Location.INSIDE) { + // merge the first arc after the sub-chord and the sub-chord + final Arc merged = new Arc(cutArcs.remove(cutArcs.size() - 1).getInf(), + highList.remove(0).getSup()); + highList.add(0, merged); + } + + list.addAll(highList); } diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java index 675bf9288..dd92ab14a 100644 --- a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java @@ -19,6 +19,7 @@ package org.apache.commons.math3.geometry.spherical.oned; import org.apache.commons.math3.geometry.Point; import org.apache.commons.math3.geometry.partitioning.Hyperplane; import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; /** This class represents a 1D oriented hyperplane on the circle. *

An hyperplane on the 1-sphere is a chord that splits @@ -41,15 +42,20 @@ public class Chord implements Hyperplane { /** Middle point of the chord. */ private final S1Point middle; + /** Tolerance below which close sub-arcs are merged together. */ + private final double tolerance; + /** Simple constructor. * @param start start angle of the chord * @param end end angle of the chord + * @param tolerance tolerance below which close sub-arcs are merged together */ - public Chord(final double start, final double end) { - this.start = start; - this.end = end; - this.middle = new S1Point(0.5 * (start + end)); - this.cos = FastMath.cos(0.5 * (end - start)); + public Chord(final double start, final double end, final double tolerance) { + this.start = start; + this.end = end; + this.middle = new S1Point(0.5 * (start + end)); + this.cos = FastMath.cos(0.5 * (end - start)); + this.tolerance = tolerance; } /** Copy the instance. @@ -66,6 +72,15 @@ public class Chord implements Hyperplane { return cos - middle.getVector().dotProduct(((S1Point) point).getVector()); } + /** Get the reverse of the instance. + *

Get a chord with reversed orientation with respect to the + * instance. A new object is built, the instance is untouched.

+ * @return a new chord, with orientation opposite to the instance orientation + */ + public Chord getReverse() { + return new Chord(end, MathUtils.normalizeAngle(start, end + FastMath.PI), tolerance); + } + /** Build a region covering the whole hyperplane. *

Since this class represent zero dimension spaces which does * not have lower dimension sub-spaces, this method returns a dummy @@ -86,7 +101,7 @@ public class Chord implements Hyperplane { * ArcsSet IntervalsSet} instance) */ public ArcsSet wholeSpace() { - return new ArcsSet(); + return new ArcsSet(tolerance); } /** {@inheritDoc} */ @@ -108,4 +123,11 @@ public class Chord implements Hyperplane { return end; } + /** Get the tolerance below which close sub-arcs are merged together. + * @return tolerance below which close sub-arcs are merged together + */ + public double getTolerance() { + return tolerance; + } + } diff --git a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java index 2922d972e..a0060729b 100644 --- a/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java @@ -22,6 +22,8 @@ import java.util.List; import org.apache.commons.math3.geometry.partitioning.Hyperplane; 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 sub-hyperplane for {@link Chord}. @@ -49,13 +51,24 @@ public class SubChord implements SubHyperplane { /** Simple constructor. * @param chord underlying hyperplane - * @param limits limit angles + * @param limits limit angles (the list will be copied into a new independent list) */ private SubChord(final Chord chord, final List limits) { this.chord = chord; this.limits = new ArrayList(limits); } + /** Get the sub-arcs. + * @return a newly created list with sub-arcs + */ + public List getSubArcs() { + final List subArcs = new ArrayList(limits.size() / 2); + for (int i = 0; i < limits.size(); i += 2) { + subArcs.add(new Arc(limits.get(i), limits.get(i + 1))); + } + return subArcs; + } + /** {@inheritDoc} */ public double getSize() { double sum = 0; @@ -67,23 +80,69 @@ public class SubChord implements SubHyperplane { /** {@inheritDoc} */ public Side side(final Hyperplane hyperplane) { - // TODO Auto-generated method stub - return null; + + final Chord testChord = (Chord) hyperplane; + final double reference = FastMath.PI + testChord.getStart(); + final double otherEnd = testChord.getEnd(); + + boolean inMinus = false; + boolean inPlus = false; + for (int i = 0; i < limits.size(); i += 2) { + inMinus = inMinus || (MathUtils.normalizeAngle(limits.get(i), reference) < otherEnd); + inPlus = inPlus || (MathUtils.normalizeAngle(limits.get(i + 1), reference) > otherEnd); + } + + if (inMinus) { + if (inPlus) { + return Side.BOTH; + } else { + return Side.MINUS; + } + } else { + if (inPlus) { + return Side.PLUS; + } else { + return Side.HYPER; + } + } + } /** {@inheritDoc} */ public SplitSubHyperplane split(final Hyperplane hyperplane) { - // TODO Auto-generated method stub - return null; + + final Chord testChord = (Chord) hyperplane; + final double reference = FastMath.PI + testChord.getStart(); + final double otherEnd = testChord.getEnd(); + + final List minus = new ArrayList(limits.size()); + final List plus = new ArrayList(limits.size()); + + for (int i = 0; i < limits.size(); i += 2) { + final double subStart = MathUtils.normalizeAngle(limits.get(i), reference); + final double subEnd = MathUtils.normalizeAngle(limits.get(i + 1), reference); + if (subStart < otherEnd) { + minus.add(subStart); + minus.add(FastMath.min(subEnd, otherEnd)); + } + if (subEnd > otherEnd) { + plus.add(FastMath.max(subStart, otherEnd)); + plus.add(subEnd); + } + } + + return new SplitSubHyperplane(new SubChord(chord, plus), + new SubChord(chord, minus)); + } /** {@inheritDoc} */ - public SubHyperplane copySelf() { + public SubChord copySelf() { return new SubChord(chord, limits); } /** {@inheritDoc} */ - public Hyperplane getHyperplane() { + public Chord getHyperplane() { return chord; } @@ -93,15 +152,101 @@ public class SubChord implements SubHyperplane { } /** {@inheritDoc} */ - public SubHyperplane reunite(SubHyperplane other) { - final List otherLimits = ((SubChord) other).limits; - final List merged = new ArrayList(limits.size() + otherLimits.size()); + public SubChord reunite(SubHyperplane other) { + + final List otherLimits = ((SubChord) other).limits; + + final List merged; + if (otherLimits.isEmpty()) { + merged = limits; + } else if (limits.isEmpty()) { + merged = otherLimits; + } else { + + merged = new ArrayList(limits.size() + otherLimits.size()); + final double reference = limits.get(0) + FastMath.PI; + + // initialize loop on first limits list + int i = 0; + int iEnd = limits.size() - 1; + boolean enteringI = true; + + // initialize loop on second limits list + int j = otherLimits.size() - 1; + double angleAfter = Double.POSITIVE_INFINITY; + for (int jSearch = 0; jSearch < otherLimits.size(); ++jSearch) { + // look for the first angle in the second list that lies just after first limits start + final double angleJ = MathUtils.normalizeAngle(otherLimits.get(jSearch), reference); + if (angleJ < angleAfter) { + j = jSearch; + angleAfter = angleJ; + } + } + int jEnd = (j + otherLimits.size() - 1) % otherLimits.size(); + boolean enteringJ = j % 2 == 0; + + // perform merging loop + boolean inMerged = !enteringJ; + double angleI = MathUtils.normalizeAngle(limits.get(i), reference); + double angleJ = MathUtils.normalizeAngle(otherLimits.get(j), reference); + while (i >= 0 || j >= 0) { + + if (i >= 0 && (j < 0 || angleI <= angleJ)) { + if (inMerged && (!enteringI) && enteringJ) { + // we were in a merged arc and exit from it + merged.add(angleI); + inMerged = false; + } else if (!inMerged && enteringI) { + // we were outside and enter into a merged arc + merged.add(angleI); + inMerged = true; + } + if (i == iEnd) { + i = -1; + } else { + ++i; + } + angleI = MathUtils.normalizeAngle(limits.get(i), reference); + enteringI = !enteringI; + } else { + if (inMerged && enteringI && (!enteringJ)) { + // we were in a merged arc and exit from it + merged.add(angleJ); + inMerged = false; + } else if (!inMerged && enteringJ) { + // we were outside and enter into a merged arc + merged.add(angleJ); + inMerged = true; + } + if (j == jEnd) { + j = -1; + } else { + j = (j + 1) % otherLimits.size(); + } + angleJ = MathUtils.normalizeAngle(otherLimits.get(j), reference); + enteringJ = !enteringJ; + } + + } + + if (inMerged) { + // we end the merging loop inside a merged arc, + // we have to put its start at the front of the limits list + if (merged.isEmpty()) { + // the merged arc covers all the circle + merged.add(0.0); + merged.add(2 * FastMath.PI); + } else { + double previousAngle = merged.get(merged.size() - 1) - 2 * FastMath.PI; + for (int k = 0; k < merged.size() - 1; ++k) { + final double tmp = merged.get(k); + merged.set(k, previousAngle); + previousAngle = tmp; + } + merged.set(merged.size() - 1, previousAngle); + } + } - int i = 0; - int j = 0; - boolean inside = false; - while (i < limits.size() || j < otherLimits.size()) { - // TODO: implement merging loop } return new SubChord(chord, merged); 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 index 2e7723459..57ac3cc05 100644 --- 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 @@ -23,18 +23,16 @@ 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.Chord; import org.apache.commons.math3.geometry.spherical.oned.ArcsSet; +import org.apache.commons.math3.geometry.spherical.oned.Chord; 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 circle on the 2-sphere. +/** This class represents an oriented great circle on the 2-sphere. - *

An oriented circle can be defined by a center point and an - * angular radius. The circle is the the set of points that are exactly - * at the specified angular radius from the center (which does not - * belong to the circle it defines except if angular radius is 0).

+ *

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.

*

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 @@ -56,32 +54,28 @@ public class Circle implements Hyperplane, EmbeddingThe circle is oriented in the trigonometric direction around center.

- * @param center circle enter - * @param radius cirle radius + /** Build a great circle from its pole. + *

The circle is oriented in the trigonometric direction around pole.

+ * @param pole circle pole + * @param tolerance tolerance below which close sub-arcs are merged together */ - public Circle(final Vector3D center, final double radius) { - reset(center, radius); + public Circle(final Vector3D pole, final double tolerance) { + reset(pole); + this.tolerance = tolerance; } - /** Build a great circle from a center only, radius being forced to \( \pi/2 \). - *

This constructor is recommended to build great circles as it does - * ensure the exact values of \( \cos(\pi/2) = 0 and \sin(\pi) = 1 \) are used.

- *

The circle is oriented in the trigonometric direction around center.

- * @param center circle enter + /** Build a great circle from two non-aligned points. + *

The circle is oriented from first to second point using the path smaller than \( \pi \).

+ * @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 Vector3D center) { - reset(center); + 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. @@ -89,18 +83,14 @@ public class Circle implements Hyperplane, Embedding, Embedding, EmbeddingThe circle is oriented in the trigonometric direction around center.

- * @param newCenter circle enter - * @param newRadius cirle radius + /** Reset the instance as if built from a pole. + *

The circle is oriented in the trigonometric direction around pole.

+ * @param newPole circle pole */ - public void reset(final Vector3D newCenter, final double newRadius) { - reset(newCenter, newRadius, FastMath.cos(radius), FastMath.sin(radius)); - } - - /** Reset the instance as if built from a center, radius being forced to \( \pi/2 \). - *

This constructor is recommended to build great circles as it does - * ensure the exact values of \( \cos(\pi/2) = 0 and \sin(\pi) = 1 \) are used.

- *

The circle is oriented in the trigonometric direction around center.

- * @param newCenter circle enter - */ - public void reset(final Vector3D newCenter) { - reset(newCenter, 0.5 * FastMath.PI, 0.0, 1.0); - } - - /** Reset the instance. - * @param newCenter circle enter - * @param newRadius cirle radius - * @param newCos cosine of the radius - * @param newSin sine of the radius - */ - private void reset(final Vector3D newCenter, final double newRadius, - final double newCos, final double newSin) { - this.pole = newCenter.normalize(); - this.x = newCenter.orthogonal(); - this.y = Vector3D.crossProduct(newCenter, x).normalize(); - this.radius = newRadius; - this.cos = newCos; - this.sin = newSin; + 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(); - radius = FastMath.PI - radius; - cos = -cos; + y = y.negate(); + pole = pole.negate(); } /** Get the reverse of the instance. @@ -168,7 +131,14 @@ public class Circle implements Hyperplane, Embedding, Embedding + * 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. + *

+ * @return an arbitrary x axis on the circle + * @see #getPointAt(double) + * @see #getYAxis() + * @see #getPole() */ - public S2Point[] intersection(final Circle other) { + public Vector3D getXAxis() { + return x; + } - // we look for a vector as v = a pole + b other.pole +/- c pole ^ other.pole - // and such that v angular separation with both centers is consistent - // with each circle radius, and v is also on the 2-sphere + /** Get the Y axis of the circle. + *

+ * 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. + *

+ * @return an arbitrary y axis point on the circle + * @see #getPointAt(double) + * @see #getXAxis() + * @see #getPole() + */ + public Vector3D getYAxis() { + return y; + } - final double dot = Vector3D.dotProduct(pole, other.pole); - final double f = 1.0 / (1.0 - dot * dot); - final double a = f * (cos - dot * other.cos); - final double b = f * (other.cos - dot * cos); - final Vector3D inPlane = new Vector3D(a, pole, b, other.pole); - final double omN2 = 1.0 - inPlane.getNormSq(); - if (omN2 <= 0) { - // no intersections (we include the just tangent case too) - return null; - } - - final double c = FastMath.sqrt(f * omN2); - final Vector3D outOfPlane = new Vector3D(c, Vector3D.crossProduct(pole, other.pole)); - - return new S2Point[] { - new S2Point(inPlane.add(outOfPlane)), - new S2Point(inPlane.subtract(outOfPlane)) - }; + /** Get the pole of the circle. + *

+ * As the circle is a great circle, the pole does not + * belong to it. + *

+ * @return pole of the circle + * @see #getXAxis() + * @see #getYAxis() + */ + public Vector3D getPole() { + return pole; + } + /** Get the chord of the instance that lies inside the other circle. + * @param other other circle + * @return chord of the instance that lies inside the other circle + * (guaranteed to always have a length of \( \pi \)) + */ + public Chord getChord(final Circle other) { + final double alpha = getPhase(other.pole); + final double halfPi = 0.5 * FastMath.PI; + return new Chord(alpha - halfPi, alpha + halfPi, tolerance); } /** {@inheritDoc} */ public SubCircle wholeHyperplane() { - return new SubCircle(this, new ArcsSet()); + return new SubCircle(this, new ArcsSet(tolerance)); } /** Build a region covering the whole space. @@ -252,7 +244,7 @@ public class Circle implements Hyperplane, Embedding, Embedding, Embedding { 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()); + } + /** Compute the distance (angular separation) between two points. * @param p1 first vector * @param p2 second vector 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 000000000..ca5fdcdf5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java @@ -0,0 +1,836 @@ +/* + * 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.List; + +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.AbstractRegion; +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.partitioning.SubHyperplane; +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; +import org.apache.commons.math3.util.MathUtils; + +/** This class represents a region on the 2-sphere: a set of spherical polygons. + * @version $Id$ + * @since 3.3 + */ +public class SphericalPolygonsSet extends AbstractRegion { + + /** Tolerance below which points are consider to be identical. */ + private final double tolerance; + + /** Boundary defined as an array of closed loops start vertices. */ + private List loops; + + /** Build a polygons set representing the whole real circle. + * @param tolerance below which points are consider to be identical + */ + public SphericalPolygonsSet(final double tolerance) { + this.tolerance = tolerance; + } + + /** Build a polygons set from a BSP tree. + *

The leaf nodes of the BSP tree must 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}

+ * @param tree inside/outside BSP tree representing the region + * @param tolerance below which points are consider to be identical + */ + public SphericalPolygonsSet(final BSPTree tree, final double tolerance) { + super(tree); + this.tolerance = tolerance; + } + + /** Build a polygons set from a Boundary REPresentation (B-rep). + *

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.

+ *

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.

+ *

If the boundary is empty, the region will represent the whole + * space.

+ * @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> boundary, final double tolerance) { + super(boundary); + this.tolerance = tolerance; + } + + /** Build a polygon from a simple list of vertices. + *

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.

+ *

This constructor does not handle polygons with a boundary + * forming several disconnected paths (such as polygons with holes).

+ *

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}.

+ *

If the list is empty, the region will represent the whole + * space.

+ *

+ * 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. + *

+ * @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)); + this.tolerance = hyperplaneThickness; + } + + /** Build the BSP tree of a polygons set from a simple list of vertices. + *

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.

+ *

This constructor does not handle polygons with a boundary + * forming several disconnected paths (such as polygons with holes).

+ *

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.

+ *

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}.

+ * @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 verticesToTree(final double hyperplaneThickness, + final S2Point ... vertices) { + + final int n = vertices.length; + if (n == 0) { + // the tree represents the whole space + return new BSPTree(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 edges = new ArrayList(n); + for (int i = 0; i < n; ++i) { + + // get the endpoints of the edge + final Vertex start = vArray[i]; + final Vertex end = vArray[(i + 1) % n]; + + // 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 tree = new BSPTree(); + insertEdges(hyperplaneThickness, tree, edges); + + return tree; + + } + + /** Recursively build a tree by inserting cut sub-hyperplanes. + * @param hyperplaneThickness tolerance below which points are consider 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 node, + final List 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 (inserted.getNode() == null) { + if (node.insertCut(inserted.getCircle())) { + inserted.setNode(node); + } else { + inserted = null; + } + } else { + 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 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 outsideList = new ArrayList(); + final List insideList = new ArrayList(); + 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.getMinus().setAttribute(Boolean.FALSE); + } + if (!insideList.isEmpty()) { + insertEdges(hyperplaneThickness, node.getMinus(), insideList); + } else { + node.getPlus().setAttribute(Boolean.TRUE); + } + + } + + /** Spherical polygons vertex. + * @see Edge + */ + public static 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 circles; + + /** Build a non-processed vertex not owned by any node yet. + * @param location vertex location + */ + private Vertex(final S2Point location) { + this.location = location; + this.incoming = null; + this.outgoing = null; + this.circles = new ArrayList(); + } + + /** 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 + */ + private void bindWith(final Circle circle) { + circles.add(circle); + } + + /** Get the common circle bound with both the instance and another vertex, if any. + *

+ * 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. + *

+ * @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 + */ + private 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. + *

+ * The circle supporting the incoming edge is automatically bound + * with the instance. + *

+ * @param incoming incoming edge + */ + private 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. + *

+ * The circle supporting the outgoing edge is automatically bound + * with the instance. + *

+ * @param outgoing outgoing edge + */ + private void setOutgoing(final Edge outgoing) { + this.outgoing = outgoing; + bindWith(outgoing.getCircle()); + } + + /** Get outgoing edge. + * @return outgoing edge + */ + public Edge getOutgoing() { + return outgoing; + } + + } + + /** Spherical polygons edge. + * @see Vertex + */ + public static 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; + + /** Node whose cut hyperplane contains this edge. */ + private BSPTree node; + + /** 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 + */ + private 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; + this.node = null; + + // 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. + *

+ * 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. + *

+ * @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())); + } + + /** Set the node whose cut hyperplane contains this edge. + * @param node node whose cut hyperplane contains this edge + */ + private void setNode(final BSPTree node) { + this.node = node; + } + + /** Get the node whose cut hyperplane contains this edge. + * @return node whose cut hyperplane contains this edge + */ + private BSPTree getNode() { + return node; + } + + /** Connect the instance with a following edge. + * @param next edge following the instance + */ + private void setNextEdge(final Edge next) { + end = next.getStart(); + end.setIncoming(this); + end.bindWith(getCircle()); + } + + /** Split the edge. + *

+ * 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. + *

+ * @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 + */ + private void split(final Circle splitCircle, + final List outsideList, final List insideList) { + + // get the inside chord, synchronizing its phase with the edge itself + final double edgeStart = circle.getPhase(start.getLocation().getVector()); + final double chordRelativeStart = MathUtils.normalizeAngle(circle.getChord(splitCircle).getStart(), + edgeStart + FastMath.PI) - edgeStart; + final double chordRelativeEnd = chordRelativeStart + FastMath.PI; + final double unwrappedEnd = chordRelativeStart - FastMath.PI; + + // build the sub-edges + if (chordRelativeStart < length) { + if (unwrappedEnd > 0) { + // the edge starts inside the circle, then goes outside, then comes back inside + + // create intermediate vertices + final Vertex vExit = new Vertex(new S2Point(circle.getPointAt(edgeStart + chordRelativeEnd))); + final Vertex vEnter = new Vertex(new S2Point(circle.getPointAt(edgeStart + chordRelativeStart))); + vExit.bindWith(splitCircle); + vEnter.bindWith(splitCircle); + + // create sub-edges + final Edge eStartIn = new Edge(start, vExit, unwrappedEnd, circle); + final Edge eMiddleOut = new Edge(vExit, vEnter, chordRelativeStart - unwrappedEnd, circle); + final Edge eEndIn = new Edge(vEnter, end, length - chordRelativeStart, circle); + eStartIn.setNode(node); + eMiddleOut.setNode(node); + eEndIn.setNode(node); + + // distribute the sub-edges in the appropriate lists + insideList.add(eStartIn); + insideList.add(eEndIn); + outsideList.add(eMiddleOut); + + } else { + // the edge starts outside of the circle, then comes inside + + // create intermediate vertices + final Vertex vEnter = new Vertex(new S2Point(circle.getPointAt(edgeStart + chordRelativeStart))); + vEnter.bindWith(splitCircle); + + // create sub-edges + final Edge eStartOut = new Edge(start, vEnter, chordRelativeStart, circle); + final Edge eEndIn = new Edge(vEnter, end, length - chordRelativeStart, circle); + eStartOut.setNode(node); + eEndIn.setNode(node); + + // distribute the sub-edges in the appropriate lists + outsideList.add(eStartOut); + insideList.add(eEndIn); + + } + } else { + if (unwrappedEnd > 0) { + if (unwrappedEnd > length) { + // the edge is entirely contained inside the circle + // we don't split anything + insideList.add(this); + } else { + // the edge starts inside the circle, then goes outside + + // create intermediate vertices + final Vertex vExit = new Vertex(new S2Point(circle.getPointAt(edgeStart + chordRelativeEnd))); + vExit.bindWith(splitCircle); + + // create sub-edges + final Edge eStartIn = new Edge(start, vExit, chordRelativeEnd, circle); + final Edge eEndOut = new Edge(vExit, end, length - chordRelativeEnd, circle); + eStartIn.setNode(node); + eEndOut.setNode(node); + + // distribute the sub-edges in the appropriate lists + insideList.add(eStartIn); + outsideList.add(eEndOut); + + } + } else { + // the edge is entirely outside of the circle + // we don't split anything + outsideList.add(this); + } + } + + } + + } + + /** {@inheritDoc} */ + @Override + public SphericalPolygonsSet buildNew(final BSPTree tree) { + return new SphericalPolygonsSet(tree, tolerance); + } + + /** {@inheritDoc} + * @exception MathIllegalStateException if the tolerance setting does not allow to build + * a clean non-ambiguous boundary + */ + @Override + protected void computeGeometricalProperties() throws MathIllegalStateException { + + final List boundary = getBoundaryLoops(); + + if (boundary.isEmpty()) { + final BSPTree tree = getTree(false); + if (tree.getCut() == null && (Boolean) tree.getAttribute()) { + // the instance covers the whole space + setSize(4 * FastMath.PI); + } else { + setSize(0); + } + setBarycenter(new S2Point(0, 0)); + } else { + + // compute some integrals around the shape + double sumArea = 0; + Vector3D sumB = Vector3D.ZERO; + + for (final Vertex startVertex : boundary) { + + int n = 0; + double sum = 0; + Vector3D sumP = Vector3D.ZERO; + for (Edge edge = startVertex.getOutgoing(); + edge.getEnd() != startVertex; + edge = edge.getEnd().getOutgoing()) { + final Vector3D middle = edge.getPointAt(0.5 * edge.getLength()); + sumP = new Vector3D(1, sumP, edge.getLength(), middle); + sum += Vector3D.angle(edge.getCircle().getPole(), + edge.getEnd().getOutgoing().getCircle().getPole()); + 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 + final double area = sum - (n - 2) * FastMath.PI; + sumArea += area; + + sumB = new Vector3D(1, sumB, area, sumP); + + } + + if (sumArea < 0) { + sumArea = 4 * FastMath.PI - sumArea; + sumB = sumB.negate(); + } + + setSize(sumArea); + final double norm = sumB.getNorm(); + if (norm == 0.0) { + setBarycenter(S2Point.NaN); + } else { + setBarycenter(new S2Point(new Vector3D(1.0 / norm, sumB))); + } + + } + + } + + /** Get the boundary loops of the polygon. + *

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.

+ *

If the polygon has no boundary at all, a zero length loop + * array will be returned.

+ *

If the polygon is a simple one-piece polygon, then the returned + * array will contain a single vertex. + *

+ *

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.

+ * @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 + */ + public List 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 root = getTree(true); + final EdgesBuilder visitor = new EdgesBuilder(root, tolerance); + root.visit(visitor); + final List edges = visitor.getEdges(); + + // convert the list of all edges into a list of start vertices + loops = new ArrayList(); + for (final Edge edge : edges) { + if (edge.getNode() != null) { + + // this is an edge belonging to a new loop, store it + final Vertex startVertex = edge.getStart(); + loops.add(startVertex); + + // disable all remaining edges in the same loop + for (Vertex vertex = edge.getEnd(); + vertex != startVertex; + vertex = vertex.getOutgoing().getEnd()) { + vertex.getIncoming().setNode(null); + } + startVertex.getIncoming().setNode(null); + + } + } + + } + } + + return Collections.unmodifiableList(loops); + + } + + /** Visitor building edges. */ + private static class EdgesBuilder implements BSPTreeVisitor { + + /** Root of the tree. */ + private final BSPTree root; + + /** Tolerance below which points are consider to be identical. */ + private final double tolerance; + + /** Boundary edges. */ + private final List edges; + + /** Simple constructor. + * @param root tree root + * @param tolerance below which points are consider to be identical + */ + public EdgesBuilder(final BSPTree root, final double tolerance) { + this.root = root; + this.tolerance = tolerance; + this.edges = new ArrayList(); + } + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree node) { + return Order.MINUS_SUB_PLUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree node) { + @SuppressWarnings("unchecked") + final BoundaryAttribute attribute = (BoundaryAttribute) 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 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 node) { + final Circle circle = (Circle) sub.getHyperplane(); + final List 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); + } + edge.setNode(node); + edges.add(edge); + } + } + + /** Get the edge that should naturally follow another one. + * @param previous edge to be continued + * @return other edge, starting where the other 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> candidates = root.getCloseCuts(point, tolerance); + + // filter the single other node we are interested in + BSPTree selected = null; + for (final BSPTree node : candidates) { + if (node != previous.getNode()) { + if (selected == null) { + selected = node; + } else { + throw new MathIllegalStateException(LocalizedFormats.OUTLINE_BOUNDARY_LOOP_OPEN); + } + } + } + if (selected == null) { + throw new MathIllegalStateException(LocalizedFormats.OUTLINE_BOUNDARY_LOOP_OPEN); + } + + // find the edge associated with the selected node + for (final Edge edge : edges) { + if (edge.getNode() == selected) { + if (Vector3D.angle(point.getVector(), edge.getStart().getLocation().getVector()) <= tolerance) { + return edge; + } + } + } + + // we should never reach this point + throw new MathIllegalStateException(LocalizedFormats.OUTLINE_BOUNDARY_LOOP_OPEN); + + } + + /** Get the boundary edges. + * @return boundary edges + * @exception MathIllegalStateException if there is not a single other edge + */ + public List getEdges() throws MathIllegalStateException { + + // connect the edges + for (final Edge previous : edges) { + previous.setNextEdge(getFollowingEdge(previous)); + } + + return edges; + + } + + } + +} 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 index 662cfa520..3945c497a 100644 --- 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 @@ -22,11 +22,9 @@ import org.apache.commons.math3.geometry.partitioning.Hyperplane; import org.apache.commons.math3.geometry.partitioning.Region; import org.apache.commons.math3.geometry.partitioning.Side; import org.apache.commons.math3.geometry.partitioning.SubHyperplane; -import org.apache.commons.math3.geometry.spherical.oned.Chord; 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.Chord; 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}. * @version $Id$ @@ -54,20 +52,18 @@ public class SubCircle extends AbstractSubHyperplane { @Override public Side side(final Hyperplane hyperplane) { - final Circle thisCircle = (Circle) getHyperplane(); - final Circle otherCircle = (Circle) hyperplane; - final S2Point[] crossings = thisCircle.intersection(otherCircle); + final Circle thisCircle = (Circle) getHyperplane(); + final Circle otherCircle = (Circle) hyperplane; + final Chord chord = thisCircle.getChord(otherCircle); - if (crossings == null) { + if (chord == null) { // the circles are disjoint - final double global = otherCircle.getOffset(thisCircle.getPointAt(0.0)); + final double global = otherCircle.getOffset(thisCircle.getXAxis()); return (global < -1.0e-10) ? Side.MINUS : ((global > 1.0e-10) ? Side.PLUS : Side.HYPER); } - // the circles do intersect - final boolean direct = FastMath.sin(thisCircle.getAngle() - otherCircle.getAngle()) < 0; - final S1Point x = thisCircle.toSubSpace(crossings); - return getRemainingRegion().side(new Chord(x, direct)); + // the circles do intersect each other + return getRemainingRegion().side(chord); } @@ -75,36 +71,33 @@ public class SubCircle extends AbstractSubHyperplane { @Override public SplitSubHyperplane split(final Hyperplane hyperplane) { - final Circle thisCircle = (Circle) getHyperplane(); - final Circle otherCircle = (Circle) hyperplane; - final S2Point[] crossings = thisCircle.intersection(otherCircle); + final Circle thisCircle = (Circle) getHyperplane(); + final Circle otherCircle = (Circle) hyperplane; + final Chord chord = thisCircle.getChord(otherCircle); - if (crossings == null) { + if (chord == null) { // the circles are disjoint - final double global = otherCircle.getOffset(thisCircle.getPointAt(0.0)); + final double global = otherCircle.getOffset(thisCircle.getXAxis()); return (global < -1.0e-10) ? new SplitSubHyperplane(null, this) : new SplitSubHyperplane(this, null); } // the circles do intersect - final boolean direct = FastMath.sin(thisCircle.getAngle() - otherCircle.getAngle()) < 0; - final S1Point x = thisCircle.toSubSpace(crossings); - final SubHyperplane subPlus = new Chord(x, !direct).wholeHyperplane(); - final SubHyperplane subMinus = new Chord(x, direct).wholeHyperplane(); - + final SubHyperplane subMinus = chord.wholeHyperplane(); + final SubHyperplane subPlus = chord.getReverse().wholeHyperplane(); final BSPTree splitTree = getRemainingRegion().getTree(false).split(subMinus); final BSPTree plusTree = getRemainingRegion().isEmpty(splitTree.getPlus()) ? new BSPTree(Boolean.FALSE) : new BSPTree(subPlus, new BSPTree(Boolean.FALSE), - splitTree.getPlus(), null); + splitTree.getPlus(), null); final BSPTree minusTree = getRemainingRegion().isEmpty(splitTree.getMinus()) ? new BSPTree(Boolean.FALSE) : new BSPTree(subMinus, new BSPTree(Boolean.FALSE), - splitTree.getMinus(), null); + splitTree.getMinus(), null); - return new SplitSubHyperplane(new SubCircle(thisCircle.copySelf(), new ArcsSet(plusTree)), - new SubCircle(thisCircle.copySelf(), new ArcsSet(minusTree))); + return new SplitSubHyperplane(new SubCircle(thisCircle.copySelf(), new ArcsSet(plusTree, thisCircle.getTolerance())), + new SubCircle(thisCircle.copySelf(), new ArcsSet(minusTree, thisCircle.getTolerance()))); }