Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Added BSP for spherical geometry. 

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

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1554650 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2014-01-01 17:26:54 +00:00
parent 37115c5788
commit a3e4cbe0bb
8 changed files with 1236 additions and 177 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java
View File

@ -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>

63
src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java
View File

@ -33,10 +33,14 @@ import org.apache.commons.math3.util.Precision;
*/
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 AbstractRegion<Sphere1D, Sphere1D> {
* </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 AbstractRegion<Sphere1D, Sphere1D> {
* 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 AbstractRegion<Sphere1D, Sphere1D> {
* <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 AbstractRegion<Sphere1D, Sphere1D> {
* </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 AbstractRegion<Sphere1D, Sphere1D> {
// 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 AbstractRegion<Sphere1D, Sphere1D> {
/** {@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 AbstractRegion<Sphere1D, Sphere1D> {
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<Arc> 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<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();
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<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);
}
list.addAll(highList);
}

34
src/main/java/org/apache/commons/math3/geometry/spherical/oned/Chord.java
View File

@ -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.
* <p>An hyperplane on the 1-sphere is a chord that splits
@ -41,15 +42,20 @@ public class Chord implements Hyperplane<Sphere1D> {
/** 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<Sphere1D> {
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<Sphere1D> {
* ArcsSet IntervalsSet} instance)
*/
public ArcsSet wholeSpace() {
return new ArcsSet();
return new ArcsSet(tolerance);
}
/** {@inheritDoc} */
@ -108,4 +123,11 @@ public class Chord implements Hyperplane<Sphere1D> {
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;
}
}

175
src/main/java/org/apache/commons/math3/geometry/spherical/oned/SubChord.java
View File

@ -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<Sphere1D> {
/** 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 SubHyperplane<Sphere1D> {
/** {@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 SubHyperplane<Sphere1D> {
}
/** {@inheritDoc} */
public SubHyperplane<Sphere1D> reunite(SubHyperplane<Sphere1D> other) {
final List<Double> otherLimits = ((SubChord) other).limits;
final List<Double> merged = new ArrayList<Double>(limits.size() + otherLimits.size());
public SubChord reunite(SubHyperplane<Sphere1D> other) {
final List<Double> otherLimits = ((SubChord) other).limits;
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);
}
}
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);

208
src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java
View File

@ -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.
* <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 Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1
/** Second axis in the equator plane, in quadrature with respect to x. */
private Vector3D y;
/** Angular radius. */
private double radius;
/** Tolerance below which close sub-arcs are merged together. */
private final double tolerance;
/** 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
/** 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, 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 \).
* <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 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 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<Sphere2D>, Embedding<Sphere2D, Sphere1
* @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 Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1
* @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 Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1
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
/** 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 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;
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<Sphere2D>, Embedding<Sphere2D, Sphere1
* @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 Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1
* @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) {
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.
* <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;
}
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.
* <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 Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1
* SphericalPolygonsSet SphericalPolygonsSet} instance)
*/
public SphericalPolygonsSet wholeSpace() {
return new SphericalPolygonsSet();
return new SphericalPolygonsSet(tolerance);
}
/** {@inheritDoc}
@ -272,7 +264,7 @@ public class Circle implements Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1
* @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 Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1
return new Circle(rotation.applyTo(circle.pole),
rotation.applyTo(circle.x),
rotation.applyTo(circle.y),
circle.radius, circle.cos, circle.sin);
circle.tolerance);
}
/** {@inheritDoc} */

7
src/main/java/org/apache/commons/math3/geometry/spherical/twod/S2Point.java
View File

@ -162,6 +162,13 @@ public class S2Point implements Point<Sphere2D> {
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

836
src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java Normal file
View File

@ -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;
}
}
}

45
src/main/java/org/apache/commons/math3/geometry/spherical/twod/SubCircle.java
View File

@ -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<Sphere2D, Sphere1D> {
@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 AbstractSubHyperplane<Sphere2D, Sphere1D> {
@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())));
}