commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1554650 - in /commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry: partitioning/ spherical/oned/ spherical/twod/
Date Wed, 01 Jan 2014 17:26:55 GMT
Author: luc
Date: Wed Jan  1 17:26:54 2014
New Revision: 1554650

URL: http://svn.apache.org/r1554650
Log:
Added BSP for spherical geometry.

This is work in progress! It is not tested yet.

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java   (with props)
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java?rev=1554650&r1=1554649&r2=1554650&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java Wed Jan  1 17:26:54 2014
@@ -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<S extends Space> {
 
     }
 
+    /** 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<BSPTree<S>> getCloseCuts(final Point<S> point, final double maxOffset) {
+        final List<BSPTree<S>> close = new ArrayList<BSPTree<S>>();
+        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<S> point, final double maxOffset,
+                                  final List<BSPTree<S>> 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.
      * <p>The condensation operation is not recursive, it must be called
      * explicitly from leaves to root.</p>

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java?rev=1554650&r1=1554649&r2=1554650&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java Wed Jan  1 17:26:54 2014
@@ -33,10 +33,14 @@ import org.apache.commons.math3.util.Pre
  */
 public class ArcsSet extends AbstractRegion<Sphere1D, Sphere1D> {
 
+    /** 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 AbstractReg
      * </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
      */
-    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 AbstractReg
      * 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
      */
-    public ArcsSet(final BSPTree<Sphere1D> tree) {
+    public ArcsSet(final BSPTree<Sphere1D> 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 AbstractReg
      * <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
      */
-    public ArcsSet(final Collection<SubHyperplane<Sphere1D>> boundary) {
+    public ArcsSet(final Collection<SubHyperplane<Sphere1D>> 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 AbstractReg
      * </p>
      * @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<Sphere1D> buildTree(final double lower, final double upper) {
+    private static BSPTree<Sphere1D> 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 AbstractReg
 
         // the two boundary angles define only one cutting chord
         final double normalizedUpper = MathUtils.normalizeAngle(upper, lower + FastMath.PI);
-        return new BSPTree<Sphere1D>(new Chord(lower, normalizedUpper).wholeHyperplane(),
+        return new BSPTree<Sphere1D>(new Chord(lower, normalizedUpper, tolerance).wholeHyperplane(),
                                      new BSPTree<Sphere1D>(Boolean.FALSE),
                                      new BSPTree<Sphere1D>(Boolean.TRUE),
                                      null);
@@ -125,7 +136,7 @@ public class ArcsSet extends AbstractReg
     /** {@inheritDoc} */
     @Override
     public ArcsSet buildNew(final BSPTree<Sphere1D> tree) {
-        return new ArcsSet(tree);
+        return new ArcsSet(tree, tolerance);
     }
 
     /** {@inheritDoc} */
@@ -181,21 +192,29 @@ public class ArcsSet extends AbstractReg
                 list.add(new Arc(lower, upper));
             }
         } else {
-            final Chord   chord = (Chord) node.getCut().getHyperplane();
-            final S1Point loc   = chord.getMiddle();
-            double        alpha = loc.getAlpha();
-
-            // make sure we explore the tree in increasing order
-            final BSPTree<Sphere1D> low  = chord.isDirect() ? node.getMinus() : node.getPlus();
-            final BSPTree<Sphere1D> 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();
+            final SubChord  cut      = (SubChord) node.getCut();
+            final List<Arc> cutArcs  = cut.getSubArcs();
+            final double    cutStart = cutArcs.get(0).getInf();
+            final double    cutEnd   = cutArcs.get(cutArcs.size() - 1).getSup();
+
+            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);
+            }
+
+            final List<Arc> highList = new ArrayList<Arc>();
+            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);
             }
-            recurseList(high, list, alpha, upper);
+
+            list.addAll(highList);
 
         }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java?rev=1554650&r1=1554649&r2=1554650&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java Wed Jan  1 17:26:54 2014
@@ -19,6 +19,7 @@ package org.apache.commons.math3.geometr
 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.
  * <p>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.
+     * <p>Get a chord with reversed orientation with respect to the
+     * instance. A new object is built, the instance is untouched.</p>
+     * @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.
      * <p>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;
+    }
+
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java?rev=1554650&r1=1554649&r2=1554650&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java Wed Jan  1 17:26:54 2014
@@ -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 SubHype
 
     /** 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<Double> limits) {
         this.chord = chord;
         this.limits = new ArrayList<Double>(limits);
     }
 
+    /** Get the sub-arcs.
+     * @return a newly created list with sub-arcs
+     */
+    public List<Arc> getSubArcs() {
+        final List<Arc> subArcs = new ArrayList<Arc>(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 SubHype
 
     /** {@inheritDoc} */
     public Side side(final Hyperplane<Sphere1D> 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<Sphere1D> split(final Hyperplane<Sphere1D> 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<Double> minus = new ArrayList<Double>(limits.size());
+        final List<Double> plus  = new ArrayList<Double>(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<Sphere1D>(new SubChord(chord, plus),
+                                                new SubChord(chord, minus));
+
     }
 
     /** {@inheritDoc} */
-    public SubHyperplane<Sphere1D> copySelf() {
+    public SubChord copySelf() {
         return new SubChord(chord, limits);
     }
 
     /** {@inheritDoc} */
-    public Hyperplane<Sphere1D> getHyperplane() {
+    public Chord getHyperplane() {
         return chord;
     }
 
@@ -93,15 +152,101 @@ public class SubChord implements SubHype
     }
 
     /** {@inheritDoc} */
-    public SubHyperplane<Sphere1D> reunite(SubHyperplane<Sphere1D> other) {
+    public SubChord reunite(SubHyperplane<Sphere1D> other) {
+
         final List<Double> otherLimits = ((SubChord) other).limits;
-        final List<Double> merged = new ArrayList<Double>(limits.size() + otherLimits.size());
 
-        int i = 0;
-        int j = 0;
-        boolean inside = false;
-        while (i < limits.size() || j < otherLimits.size()) {
-            // TODO: implement merging loop
+        final List<Double> merged;
+        if (otherLimits.isEmpty()) {
+            merged = limits;
+        } else if (limits.isEmpty()) {
+            merged = otherLimits;
+        } else {
+
+            merged = new ArrayList<Double>(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);
+                }
+            }
+
         }
 
         return new SubChord(chord, merged);

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java?rev=1554650&r1=1554649&r2=1554650&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java Wed Jan  1 17:26:54 2014
@@ -23,18 +23,16 @@ import org.apache.commons.math3.geometry
 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.
 
- * <p>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).</p>
+ * <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
@@ -56,32 +54,28 @@ public class Circle implements Hyperplan
     /** Second axis in the equator plane, in quadrature with respect to x. */
     private Vector3D y;
 
-    /** Angular radius. */
-    private double radius;
-
-    /** Cosine of the radius. */
-    private double cos;
-
-    /** Sine of the radius. */
-    private double sin;
-
-    /** Build a circle from a center and a radius.
-     * <p>The circle is oriented in the trigonometric direction around center.</p>
-     * @param center circle enter
-     * @param radius cirle radius
-     */
-    public Circle(final Vector3D center, final double radius) {
-        reset(center, radius);
-    }
+    /** Tolerance below which close sub-arcs are merged together. */
+    private final double tolerance;
 
-    /** Build a great circle from a center only, radius being forced to \( \pi/2 \).
-     * <p>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.</p>
-     * <p>The circle is oriented in the trigonometric direction around center.</p>
-     * @param center circle enter
+    /** 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 center) {
-        reset(center);
+    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.
@@ -89,18 +83,14 @@ public class Circle implements Hyperplan
      * @param pole circle pole
      * @param x first axis in the equator plane
      * @param y second axis in the equator plane
-     * @param radius cirle radius
-     * @param cos cosine of the radius
-     * @param sin sine of the radius
+     * @param tolerance tolerance below which close sub-arcs are merged together
      */
     private Circle(final Vector3D pole, final Vector3D x, final Vector3D y,
-                   final double radius, final double cos, final double sin) {
-        this.pole   = pole;
-        this.x      = x;
-        this.y      = y;
-        this.radius = radius;
-        this.cos    = cos;
-        this.sin    = sin;
+                   final double tolerance) {
+        this.pole      = pole;
+        this.x         = x;
+        this.y         = y;
+        this.tolerance = tolerance;
     }
 
     /** Copy constructor.
@@ -109,7 +99,7 @@ public class Circle implements Hyperplan
      * @param circle circle to copy
      */
     public Circle(final Circle circle) {
-        this(circle.pole, circle.x, circle.y, circle.radius, circle.cos, circle.sin);
+        this(circle.pole, circle.x, circle.y, circle.tolerance);
     }
 
     /** {@inheritDoc} */
@@ -117,49 +107,22 @@ public class Circle implements Hyperplan
         return new Circle(this);
     }
 
-    /** Reset the instance as if built from a center and a radius.
-     * <p>The circle is oriented in the trigonometric direction around center.</p>
-     * @param newCenter circle enter
-     * @param newRadius cirle radius
-     */
-    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 \).
-     * <p>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.</p>
-     * <p>The circle is oriented in the trigonometric direction around center.</p>
-     * @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;
+    /** 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();
-        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 Hyperplan
      * @return a new circle, with orientation opposite to the instance orientation
      */
     public Circle getReverse() {
-        return new Circle(pole.negate(), x, y.negate(), FastMath.PI - radius, -cos, sin);
+        return new Circle(pole.negate(), x, y.negate(), tolerance);
+    }
+
+    /** 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;
     }
 
     /** {@inheritDoc}
@@ -203,48 +173,70 @@ public class Circle implements Hyperplan
      * @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) {
-        final double cosAlpha = FastMath.cos(alpha);
-        final double sinAlpha = FastMath.sin(alpha);
-        return new Vector3D(cosAlpha * sin, x, sinAlpha * sin, y, cos, pole);
+        return new Vector3D(FastMath.cos(alpha), x, FastMath.sin(alpha), y);
     }
 
-    /** Get the intersection points of the instance and another circle.
-     * @param other other circle
-     * @return intersection points of the instance and the other circle
-     * or null if there are no intersection points
+    /** 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 S2Point[] intersection(final Circle other) {
-
-        // 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
-
-        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;
-        }
+    public Vector3D getXAxis() {
+        return x;
+    }
 
-        final double c = FastMath.sqrt(f * omN2);
-        final Vector3D outOfPlane = new Vector3D(c, Vector3D.crossProduct(pole, other.pole));
+    /** 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;
+    }
 
-        return new S2Point[] {
-            new S2Point(inPlane.add(outOfPlane)),
-            new S2Point(inPlane.subtract(outOfPlane))
-        };
+    /** 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 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 Hyperplan
      * SphericalPolygonsSet SphericalPolygonsSet} instance)
      */
     public SphericalPolygonsSet wholeSpace() {
-        return new SphericalPolygonsSet();
+        return new SphericalPolygonsSet(tolerance);
     }
 
     /** {@inheritDoc}
@@ -272,7 +264,7 @@ public class Circle implements Hyperplan
      * @see #getOffset(Point)
      */
     public double getOffset(final Vector3D direction) {
-        return Vector3D.angle(pole, direction) - radius;
+        return Vector3D.angle(pole, direction) - 0.5 * FastMath.PI;
     }
 
     /** {@inheritDoc} */
@@ -325,7 +317,7 @@ public class Circle implements Hyperplan
             return new Circle(rotation.applyTo(circle.pole),
                               rotation.applyTo(circle.x),
                               rotation.applyTo(circle.y),
-                              circle.radius, circle.cos, circle.sin);
+                              circle.tolerance);
         }
 
         /** {@inheritDoc} */

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java?rev=1554650&r1=1554649&r2=1554650&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java Wed Jan  1 17:26:54 2014
@@ -162,6 +162,13 @@ public class S2Point implements Point<Sp
         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

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java?rev=1554650&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java Wed Jan  1 17:26:54 2014
@@ -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<Sphere2D, Sphere1D> {
+
+    /** 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<Vertex> 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.
+     * <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);
+        this.tolerance = 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);
+        this.tolerance = 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 #PolygonsSet(Collection) 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
+     * &pi; 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));
+        this.tolerance = hyperplaneThickness;
+    }
+
+    /** 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);
+        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<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 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<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 (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<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.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<Circle> 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<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
+         */
+        private 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
+         */
+        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.
+         * <p>
+         * The circle supporting the incoming edge is automatically bound
+         * with the instance.
+         * </p>
+         * @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.
+         * <p>
+         * The circle supporting the outgoing edge is automatically bound
+         * with the instance.
+         * </p>
+         * @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<Sphere2D> 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.
+         * <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()));
+        }
+
+        /** Set the node whose cut hyperplane contains this edge.
+         * @param node node whose cut hyperplane contains this edge
+         */
+        private void setNode(final BSPTree<Sphere2D> node) {
+            this.node = node;
+        }
+
+        /** Get the node whose cut hyperplane contains this edge.
+         * @return node whose cut hyperplane contains this edge
+         */
+        private BSPTree<Sphere2D> 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.
+         * <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
+         */
+        private void split(final Circle splitCircle,
+                           final List<Edge> outsideList, final List<Edge> 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<Sphere2D> 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<Vertex> boundary = getBoundaryLoops();
+
+        if (boundary.isEmpty()) {
+            final BSPTree<Sphere2D> 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.
+     * <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
+     */
+    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, tolerance);
+                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>();
+                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<Sphere2D> {
+
+        /** Root of the tree. */
+        private final BSPTree<Sphere2D> root;
+
+        /** Tolerance below which points are consider to be identical. */
+        private final double tolerance;
+
+        /** Boundary edges. */
+        private final List<Edge> edges;
+
+        /** Simple constructor.
+         * @param root tree root
+         * @param tolerance below which points are consider to be identical
+         */
+        public EdgesBuilder(final BSPTree<Sphere2D> root, final double tolerance) {
+            this.root      = root;
+            this.tolerance = tolerance;
+            this.edges     = new ArrayList<Edge>();
+        }
+
+        /** {@inheritDoc} */
+        public Order visitOrder(final BSPTree<Sphere2D> node) {
+            return Order.MINUS_SUB_PLUS;
+        }
+
+        /** {@inheritDoc} */
+        public void visitInternalNode(final BSPTree<Sphere2D> node) {
+            @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);
+                }
+                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<BSPTree<Sphere2D>> candidates = root.getCloseCuts(point, tolerance);
+
+            // filter the single other node we are interested in
+            BSPTree<Sphere2D> selected = null;
+            for (final BSPTree<Sphere2D> 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<Edge> getEdges() throws MathIllegalStateException {
+
+            // connect the edges
+            for (final Edge previous : edges) {
+                previous.setNextEdge(getFollowingEdge(previous));
+            }
+
+            return edges;
+
+        }
+
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java?rev=1554650&r1=1554649&r2=1554650&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java Wed Jan  1 17:26:54 2014
@@ -22,11 +22,9 @@ import org.apache.commons.math3.geometry
 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 AbstractS
     @Override
     public Side side(final Hyperplane<Sphere2D> 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 AbstractS
     @Override
     public SplitSubHyperplane<Sphere2D> split(final Hyperplane<Sphere2D> 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<Sphere2D>(null, this) :
                    new SplitSubHyperplane<Sphere2D>(this, null);
         }
 
         // the circles do intersect
-        final boolean direct = FastMath.sin(thisCircle.getAngle() - otherCircle.getAngle()) < 0;
-        final S1Point x      = thisCircle.toSubSpace(crossings);
-        final SubHyperplane<Sphere1D> subPlus  = new Chord(x, !direct).wholeHyperplane();
-        final SubHyperplane<Sphere1D> subMinus = new Chord(x,  direct).wholeHyperplane();
-
+        final SubHyperplane<Sphere1D> subMinus = chord.wholeHyperplane();
+        final SubHyperplane<Sphere1D> subPlus  = chord.getReverse().wholeHyperplane();
         final BSPTree<Sphere1D> splitTree = getRemainingRegion().getTree(false).split(subMinus);
         final BSPTree<Sphere1D> plusTree  = getRemainingRegion().isEmpty(splitTree.getPlus()) ?
                                                new BSPTree<Sphere1D>(Boolean.FALSE) :
                                                new BSPTree<Sphere1D>(subPlus, new BSPTree<Sphere1D>(Boolean.FALSE),
-                                                                        splitTree.getPlus(), null);
+                                                                     splitTree.getPlus(), null);
         final BSPTree<Sphere1D> minusTree = getRemainingRegion().isEmpty(splitTree.getMinus()) ?
                                                new BSPTree<Sphere1D>(Boolean.FALSE) :
                                                new BSPTree<Sphere1D>(subMinus, new BSPTree<Sphere1D>(Boolean.FALSE),
-                                                                        splitTree.getMinus(), null);
+                                                                     splitTree.getMinus(), null);
 
-        return new SplitSubHyperplane<Sphere2D>(new SubCircle(thisCircle.copySelf(), new ArcsSet(plusTree)),
-                                                new SubCircle(thisCircle.copySelf(), new ArcsSet(minusTree)));
+        return new SplitSubHyperplane<Sphere2D>(new SubCircle(thisCircle.copySelf(), new ArcsSet(plusTree, thisCircle.getTolerance())),
+                                                new SubCircle(thisCircle.copySelf(), new ArcsSet(minusTree, thisCircle.getTolerance())));
 
     }
 



Mime
View raw message