This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
new 09fa9bf Try harder to avoid unreasonably long chain of Bézier curves.
First implementation of createCircularRegion2D(…), not yet working because of
infinite derivative values.
09fa9bf is described below
commit 09fa9bf7e972ac665e9a07b3ca63fcb4a62e2b2d
Author: Martin Desruisseaux <[email protected]>
AuthorDate: Wed May 22 20:48:12 2019 +0200
Try harder to avoid unreasonably long chain of Bézier curves.
First implementation of createCircularRegion2D(…), not yet working because
of infinite derivative values.
---
.../internal/referencing/PositionTransformer.java | 18 ++
.../sis/internal/referencing/j2d/Bezier.java | 218 ++++++++++++++-------
.../apache/sis/referencing/GeodeticCalculator.java | 173 +++++++++++++---
.../sis/referencing/GeodeticCalculatorTest.java | 88 ++++++---
ide-project/NetBeans/nbproject/genfiles.properties | 2 +-
ide-project/NetBeans/nbproject/project.xml | 1 +
6 files changed, 375 insertions(+), 125 deletions(-)
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
index 87c2661..6ad5976 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
@@ -181,6 +181,24 @@ public final class PositionTransformer extends
GeneralDirectPosition {
}
/**
+ * Transforms the given position from the CRS of this position to the
default CRS.
+ * The result is stored in the given array.
+ *
+ * @param point the coordinates of the point to transform in-place.
+ * @throws TransformException if a coordinate transformation was required
and failed.
+ */
+ public void transform(final double[] point) throws TransformException {
+ if (point != null) {
+ if (lastCRS != defaultCRS) {
+ setSourceCRS(defaultCRS);
+ }
+ if (forward != null) {
+ forward.transform(point, 0, point, 0, 1);
+ }
+ }
+ }
+
+ /**
* Transforms a given position from its CRS to the CRS of this {@code
PositionTransformer}.
* If the CRS associated to the given position is {@code null}, then that
CRS is assumed to
* be the default CRS specified at construction time. Otherwise if that
CRS is not equal to
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/Bezier.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/Bezier.java
index 9ea2294..1c9d580 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/Bezier.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/Bezier.java
@@ -17,6 +17,9 @@
package org.apache.sis.internal.referencing.j2d;
import java.awt.geom.Path2D;
+import java.awt.geom.Line2D;
+import java.awt.geom.QuadCurve2D;
+import java.awt.geom.CubicCurve2D;
import org.opengis.referencing.operation.TransformException;
import static java.lang.Math.abs;
@@ -24,7 +27,17 @@ import static java.lang.Double.isInfinite;
/**
- * Helper class for appending Bézier curves to a path.
+ * Helper class for appending Bézier curves to a path. This class computes
cubic Bézier from 3 points and two slopes.
+ * The three points are the start point (t=0), middle point (t=½) and end
point (t=1). The slopes are at start point
+ * and end point. The Bézier curve will obey exactly to those conditions (up
to rounding errors).
+ *
+ * <p>After creating each Bézier curve, this class performs a quality control
using 2 additional points located at
+ * one quarter (t≈¼) and three quarters (t≈¾) of the curve. If the distance
between given points and a close point
+ * on the curve is greater than {@linkplain #εx} and {@linkplain #εy}
thresholds, then the curve is divided in two
+ * smaller curves and the process is repeated until curves meet the quality
controls.</p>
+ *
+ * <p>If a quadratic curve degenerates to a cubic curve or a straight line,
ignoring errors up to {@linkplain #εx}
+ * and {@linkplain #εy}, then {@link CubicCurve2D} are replaced by {@link
QuadCurve2D} or {@link Line2D}.</p>
*
* @author Martin Desruisseaux (Geomatys)
* @version 1.0
@@ -36,14 +49,22 @@ import static java.lang.Double.isInfinite;
*/
public abstract class Bezier {
/**
+ * Limit the maximal depth for avoiding shapes of unreasonable size.
+ *
+ * @see #depth
+ */
+ private static final int DEPTH_LIMIT = 10;
+
+ /**
* The path where to append Bézier curves.
*/
private final Path2D path;
/**
- * Maximal distance on <var>x</var> and <var>y</var> axis between the
cubic Bézier curve and quadratic
- * or linear simplifications. Initial value is zero. Subclasses should set
the values either in their
- * constructor or at {@link #evaluateAt(double)} method call.
+ * Maximal distance (approximate) on <var>x</var> and <var>y</var> axis
between the cubic Bézier curve and quadratic
+ * or linear simplifications. Initial value is zero. Subclasses should set
the values either in their constructor or
+ * at {@link #evaluateAt(double)} method call. Can be set to infinity or
NaN for disabling the quality checks, or to
+ * a negative value for forcing an unconditional division of a Bézier
curve in two more curves.
*/
protected double εx, εy;
@@ -60,6 +81,13 @@ public abstract class Bezier {
protected final double[] point;
/**
+ * The number of times a curve has been divided in smaller curves.
+ *
+ * @see #DEPTH_LIMIT
+ */
+ protected int depth;
+
+ /**
* Creates a new builder.
*
* @param dimension length of the {@link #point} array. Must be at least
2.
@@ -91,6 +119,26 @@ public abstract class Bezier {
protected abstract double evaluateAt(double t) throws TransformException;
/**
+ * Returns whether we should accept the given coordinates. This method is
invoked when this {@code Bezier} class thinks
+ * that the given point is not valid, but can not said for sure because
computing an exact answer would be expansive
+ * (because of the difficulty to map Bézier parameter <var>t</var>=¼ and ¾
to a distance on the curve).
+ *
+ * <p>The default implementation returns always {@code false}, since this
method is invoked only when this class
+ * thinks that the point is not valid.</p>
+ *
+ * @param x first coordinate value of a point on the curve to evaluate
for fitness.
+ * @param y second coordinate value of a point on the curve to evaluate
for fitness.
+ * @return whether the given point is close enough to a valid point.
+ * @throws TransformException if a coordinate operations was required and
failed.
+ *
+ * @todo This method is a hack. We should replace it by a computation of
the Bézier <var>t</var> parameter
+ * of the point closest to the given (x₁,y₁) and (x₃,y₃) points.
+ */
+ protected boolean isValid(double x, double y) throws TransformException {
+ return false;
+ }
+
+ /**
* Creates a sequence of Bézier curves from the position given by {@code
evaluateAt(0)} to the position given
* by {@code evaluateAt(0)}. This method determines the number of
intermediate points required for achieving
* the precision requested by the {@code εx} and {@code εy} parameters
given at construction time.
@@ -140,20 +188,24 @@ public abstract class Bezier {
final double tx = εx;
final double ty = εy;
final double α1 = evaluateAt(0.75*tmin + 0.25*tmax); //
`point` become (x₁,y₁) with this method call.
+ final double xi = point[0]; //
Must be after the call to `evaluateAt(t₁)`.
+ final double yi = point[1];
if (tx < εx) εx = tx; //
Take smallest tolerance values.
if (ty < εy) εy = ty;
- if (curve(point[0], point[1], x2, y2, x3, y3, x4, y4, α4)) { //
Must be after the call to `evaluateAt(t₁)`.
+ if (curve(xi, yi, x2, y2, x3, y3, x4, y4, α4)) {
x0 = x4;
y0 = y4;
α0 = α4;
} else {
+ depth++;
final double εxo = εx;
final double εyo = εy;
final double th = 0.5 * (tmin + tmax);
- sequence(tmin, th, point[0], point[1], α1, x2, y2, α2);
- sequence(th, tmax, x3, y3, α3, x4, y4, α4);
+ sequence(tmin, th, xi, yi, α1, x2, y2, α2);
+ sequence(th, tmax, x3, y3, α3, x4, y4, α4);
εx = εxo;
εy = εyo;
+ depth--;
}
}
@@ -175,9 +227,9 @@ public abstract class Bezier {
* If the curve is not accurate enough, then this method does nothing and
return {@code false}. In that case, the caller
* should split the curve in two smaller parts and invoke this method
again.</p>
*
- * <p>If the full equation is required for representing the curve, then
this method appends a {@link java.awt.geom.CubicCurve2D}.
- * If the same curve can be represented by a quadratic curve, then this
method appends a {@link java.awt.geom.QuadCurve2D}.
- * If the curve is actually a straight line, then this method appends a
{@link java.awt.geom.Line2D}.</p>
+ * <p>If the full equation is required for representing the curve, then
this method appends a {@link CubicCurve2D}.
+ * If the same curve can be represented by a quadratic curve, then this
method appends a {@link QuadCurve2D}.
+ * If the curve is actually a straight line, then this method appends a
{@link Line2D}.</p>
*
* @param x1 <var>x</var> coordinate at <var>t</var>≈¼ (quality control,
may be NaN).
* @param y1 <var>y</var> coordinate at <var>t</var>≈¼ (quality control,
may be NaN).
@@ -193,7 +245,8 @@ public abstract class Bezier {
private boolean curve(double x1, double y1,
double x2, double y2,
double x3, double y3,
- final double x4, final double y4, final double α4)
+ final double x4, final double y4,
+ final double α4) throws TransformException
{
/*
* Equations in this method are simplified as if (x₀,y₀) coordinates
are (0,0).
@@ -224,11 +277,13 @@ public abstract class Bezier {
* they are not used in Line2D computation. This allows
`evaluateAt(t)` method to return NaN if it can
* not provide a point for a given `t` value.
*/
- final double slope = Δy / Δx;
- if ((!isVertical && (abs(x2*slope - y2) > εy || abs(x1*slope
- y1) > εy || abs(x3*slope - y3) > εy)) ||
- (!isHorizontal && (abs(y2/slope - x2) > εx || abs(y1/slope
- x1) > εx || abs(y3/slope - x3) > εx)))
- {
- return false;
+ if (depth < DEPTH_LIMIT) {
+ final double slope = Δy / Δx;
+ if ((!isVertical && (abs(x2*slope - y2) > εy ||
abs(x1*slope - y1) > εy || abs(x3*slope - y3) > εy)) ||
+ (!isHorizontal && (abs(y2/slope - x2) > εx ||
abs(y1/slope - x1) > εx || abs(y3/slope - x3) > εx)))
+ {
+ return false;
+ }
}
path.lineTo(x4, y4);
return true;
@@ -281,66 +336,74 @@ public abstract class Bezier {
* If any of (x₁,y₁) and (x₃,y₃) coordinates is NaN, then this method
accepts the curve as valid.
* This allows `evaluateAt(t)` to return NaN if it can not provide a
point for a given `t` value.
*/
- double xi = 27./64*ax + 9./64*bx + 1./64*Δx; // "xi" is for "x
interpolated (on curve)".
- double yi = 27./64*ay + 9./64*by + 1./64*Δy;
- if (abs(xi - x1) > εx || abs(yi - y1) > εy) {
+ if (depth < DEPTH_LIMIT) {
+ double xi = 27./64*ax + 9./64*bx + 1./64*Δx; // "xi" is for
"x interpolated (on curve)".
+ double yi = 27./64*ay + 9./64*by + 1./64*Δy;
+ if (abs(xi - x1) > εx || abs(yi - y1) > εy) {
+ /*
+ * Above code tested (x,y) coordinates at t=¼ exactly (we will
test t=¾ later). However this t value does not
+ * necessarily correspond to one quarter of the distance,
because the speed at which t varies is not the same
+ * than the speed at which Bézier curve length increases.
Unfortunately computing the t values at a given arc
+ * length is complicated. We tested an approach based on
computing the y value on the curve for a given x value
+ * by starting from the Bézier curve equation:
+ *
+ * x(t) = x₀(1-t)³ + 3aₓ(1-t)²t + 3bₓ(1-t)t² + x₄t³
+ *
+ * rearranged as:
+ *
+ * (-x₀ + 3aₓ - 3bₓ + x₄)t³ + (3x₀ - 6aₓ + 3bₓ)t² + (-3x₀
+ 3aₓ)t + x₀ - x = 0
+ *
+ * and finding the roots with the CubicCurve2D.solveCubic(…)
method. However the results were worst than using
+ * fixed t values. If we want to improve on that in a future
version, we would need a function for computing
+ * arc length (for example based on
https://pomax.github.io/bezierinfo/#arclength), then use iterative method
+ * like
https://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf
(retrieved May 2019).
+ *
+ * Instead we perform another test using the tangent of the
curve at point P₁ (and later P₃).
+ *
+ * x′(t) = 3(1-t)²(aₓ - x₀) + 6t(1-t)(bₓ - aₓ) + 3t²(x₄ -
bₓ) and same for y′(t).
+ *
+ * The derivatives give us a straight line (the tangent) as an
approximation of the curve around P₁.
+ * We can then compute the point on that line which is nearest
to P₁. It should be close to current
+ * (xi,yi) coordinates, but may vary a little bit.
+ */
+ double slope = (27./16*ay + 18./16*(by-ay) + 3./16*(y4-by))
+ / (27./16*ax + 18./16*(bx-ax) + 3./16*(x4-bx));
// ∂y/∂x at t=¼.
+ double offset = (yi - slope*xi);
// Value of y at x=0.
+ xi = ((y1 - offset) * slope + x1) / (slope*slope + 1);
// NaN if slope is infinite.
+ yi = offset + xi*slope;
// Closer (xi,yi) coordinates.
+ if (abs(xi - x1) > εx || abs(yi - y1) > εy) {
+ if (!isValid(xi + x0, yi + y0)) return false;
+ }
+ /*
+ * At this point we consider (x₁,y₁) close enough even if the
initial test considered it too far.
+ * This decision is based on the assumption that the straight
line is an approximation good enough
+ * in the vicinity of P₁. We did not verified that assumption.
If we want to improve on that in a
+ * future version, we could use the second derivative:
+ *
+ * x″(t) = 6(1-t)(x₀ - 2aₓ + bₓ) + 6t(aₓ - 2bₓ + x₄)
and same for y″(t).
+ *
+ * Applying chain rule: ∂²y/∂x² = y″/x′² + y′/x″
+ *
+ * We could then estimate the change of slope at the new
(xi,yi) compared to the initial (xi,yi)
+ * and verify that this change is below some threshold. We do
not perform that check for now on
+ * the assumption that the Bézier curve is smooth enough in
the context of map projections.
+ */
+ }
/*
- * Above code tested (x,y) coordinates at t=¼ exactly (we will
test t=¾ later). However this t value does not
- * necessarily correspond to one quarter of the distance, because
the speed at which t varies is not the same
- * than the speed at which Bézier curve length increases.
Unfortunately computing the t values at a given arc
- * length is complicated. We tested an approach based on computing
the y value on the curve for a given x value
- * by starting from the Bézier curve equation:
- *
- * x(t) = x₀(1-t)³ + 3a(1-t)²t + 3b(1-t)t² + x₄t³
- *
- * rearranged as:
- *
- * (-x₀ + 3a - 3b + x₄)t³ + (3x₀ - 6a + 3b)t² + (-3x₀ + 3a)t +
x₀ - x = 0
- *
- * and finding the roots with the CubicCurve2D.solveCubic(…)
method. However the results were worst than using
- * fixed t values. If we want to improve on that in a future
version, we would need a function for computing
- * arc length (for example based on
https://pomax.github.io/bezierinfo/#arclength), then use iterative method
- * like
https://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf
(retrieved May 2019).
- *
- * Instead we perform another test using the tangent of the curve
at point P₁ (and later P₃).
- *
- * x′(t) = 3(1-t)²(a-x₀) + 6(1-t)t(b-a) + 3t²(x₄ - b)
and same for y′(t).
- *
- * The derivatives give us a straight line (the tangent) as an
approximation of the curve around P₁.
- * We can then compute the point on that line which is nearest to
P₁. It should be close to current
- * (xi,yi) coordinates, but may vary a little bit.
+ * Same than above, but with point P₁ replaced by P₃ and t=¼
replaced by t=¾.
+ * The change of t value changes the coefficients in formulas
below.
*/
- double slope = (27./16*ay + 18./16*(by-ay) + 3./16*(y4-by))
- / (27./16*ax + 18./16*(bx-ax) + 3./16*(x4-bx));
// ∂y/∂x at t=¼.
- double offset = (yi - slope*xi);
// Value of y at x=0.
- xi = ((y1 - offset) * slope + x1) / (slope*slope + 1);
// Closer (xi,yi) coordinates.
- yi = offset + xi*slope;
- xi = abs(xi - x1);
// At this point (xi,yi) are distances.
- yi = abs(yi - y1);
- if ((xi > εx && xi != Double.POSITIVE_INFINITY) ||
- (yi > εy && yi != Double.POSITIVE_INFINITY))
- {
- return false;
- }
- }
- /*
- * Same than above, but with point P₁ replaced by P₃ and t=¼ replaced
by t=¾.
- * The change of t value changes the coefficients in formulas below.
- */
- xi = 9./64*ax + 27./64*bx + 27./64*Δx;
- yi = 9./64*ay + 27./64*by + 27./64*Δy;
- if (abs(xi - x3) > εx || abs(yi - y3) > εy) {
- double slope = (3./16*ay + 18./16*(by-ay) + 27./16*(y4-by))
- / (3./16*ax + 18./16*(bx-ax) + 27./16*(x4-bx));
- double offset = (yi - slope*xi);
- xi = ((y3 - offset) * slope + x3) / (slope*slope + 1);
- yi = offset + xi*slope;
- xi = abs(xi - x3);
- yi = abs(yi - y3);
- if ((xi > εx && xi != Double.POSITIVE_INFINITY) ||
- (yi > εy && yi != Double.POSITIVE_INFINITY))
- {
- return false;
+ xi = 9./64*ax + 27./64*bx + 27./64*Δx;
+ yi = 9./64*ay + 27./64*by + 27./64*Δy;
+ if (abs(xi - x3) > εx || abs(yi - y3) > εy) {
+ double slope = (3./16*ay + 18./16*(by-ay) + 27./16*(y4-by))
+ / (3./16*ax + 18./16*(bx-ax) + 27./16*(x4-bx));
+ double offset = (yi - slope*xi);
+ xi = ((y3 - offset) * slope + x3) / (slope*slope + 1);
+ yi = offset + xi*slope;
+ if (abs(xi - x3) > εx || abs(yi - y3) > εy) {
+ if (!isValid(xi + x0, yi + y0)) return false;
+ }
}
}
/*
@@ -358,10 +421,13 @@ public abstract class Bezier {
* We compute C both ways and check if they are close enough to each
other:
*
* ΔC = (3⋅(B - A) - (P₄ - P₀))/2
+ *
+ * We multiply tolerance factor by 2 because of moving the quadratic
curve control point by 1 can move the closest
+ * point on the curve by at most ½.
*/
final double Δqx, Δqy;
- if (abs(Δqx = (3*(bx - ax) - Δx)/2) <= εx && // P₀ is zero.
- abs(Δqy = (3*(by - ay) - Δy)/2) <= εy)
+ if (abs(Δqx = (3*(bx - ax) - Δx)/2) <= 2*εx && // P₀ is zero.
+ abs(Δqy = (3*(by - ay) - Δy)/2) <= 2*εy)
{
final double qx = (3*ax + Δqx)/2; // Take average of
2 control points.
final double qy = (3*ay + Δqy)/2;
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
index 9bce0d0..2ab39b2 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
@@ -17,6 +17,7 @@
package org.apache.sis.referencing;
import java.awt.Shape;
+import java.awt.geom.Path2D;
import java.util.Locale;
import java.io.IOException;
import java.io.UncheckedIOException;
@@ -628,17 +629,22 @@ public class GeodeticCalculator {
* <li>The last point is {@link #getEndPoint()}, potentially with 360°
added or subtracted to the longitude.</li>
* </ol>
*
+ * This method tries to stay within the given tolerance threshold of the
geodesic track.
+ * The {@code tolerance} parameter should not be too small for avoiding
creation of unreasonably long chain of Bézier curves.
+ * For example a value of 1/10 of geodesic length may be sufficient.
+ *
* <b>Limitations:</b>
* This method depends on the presence of {@code java.desktop} module.
This limitations may be addressed
* in a future Apache SIS version (see <a
href="https://issues.apache.org/jira/browse/SIS-453";>SIS-453</a>).
*
* @param tolerance maximal error between the approximated curve and
actual geodesic track
* in the units of measurement given by {@link
#getDistanceUnit()}.
+ * This is approximate; the actual errors may vary
around that value.
* @return an approximation of geodesic track as Bézier curves in a Java2D
object.
* @throws TransformException if the coordinates can not be transformed to
{@linkplain #getPositionCRS() position CRS}.
* @throws IllegalStateException if some required properties have not been
specified.
*/
- public Shape toGeodesicPath2D(final double tolerance) throws
TransformException {
+ public Shape createGeodesicPath2D(final double tolerance) throws
TransformException {
ArgumentChecks.ensureStrictlyPositive("tolerance", tolerance);
if (isInvalid(START_POINT | STARTING_AZIMUTH | END_POINT |
ENDING_AZIMUTH | GEODESIC_DISTANCE)) {
if (isInvalid(END_POINT)) {
@@ -647,8 +653,8 @@ public class GeodeticCalculator {
computeDistance();
}
}
- final PathBuilder bezier = new PathBuilder(toDegrees(tolerance /
radius));
- final Shape path;
+ final PathBuilder bezier = new PathBuilder(tolerance);
+ final Path2D path;
try {
path = bezier.build();
} finally {
@@ -658,17 +664,53 @@ public class GeodeticCalculator {
}
/**
+ * Creates an approximation of the region at a constant geodesic distance
around the start point.
+ * The returned shape is circlelike with the {@linkplain #getStartPoint()
start point} in its center.
+ * The coordinates are expressed in the coordinate reference system
specified at creation time.
+ * The approximation uses cubic Bézier curves.
+ *
+ * <p>This method tries to stay within the given tolerance threshold of
the geodesic track.
+ * The {@code tolerance} parameter should not be too small for avoiding
creation of unreasonably long chain of Bézier curves.
+ * For example a value of 1/10 of geodesic length may be sufficient.</p>
+ *
+ * <b>Limitations:</b>
+ * This method depends on the presence of {@code java.desktop} module.
This limitations may be addressed
+ * in a future Apache SIS version (see <a
href="https://issues.apache.org/jira/browse/SIS-453";>SIS-453</a>).
+ *
+ * @param tolerance maximal error in the units of measurement given by
{@link #getDistanceUnit()}.
+ * This is approximate; the actual errors may vary
around that value.
+ * @return an approximation of circular region as a Java2D object.
+ * @throws TransformException if the coordinates can not be transformed to
{@linkplain #getPositionCRS() position CRS}.
+ * @throws IllegalStateException if some required properties have not been
specified.
+ */
+ public Shape createCircularRegion2D(final double tolerance) throws
TransformException {
+ ArgumentChecks.ensureStrictlyPositive("tolerance", tolerance);
+ if (isInvalid(START_POINT | GEODESIC_DISTANCE)) {
+ computeDistance();
+ }
+ final CircularPath bezier = new CircularPath(tolerance);
+ final Path2D path;
+ try {
+ path = bezier.build();
+ } finally {
+ bezier.reset();
+ }
+ path.closePath();
+ return path;
+ }
+
+ /**
* Builds a geodesic path as a sequence of Bézier curves. The start point
and end points are the points
* in enclosing {@link GeodeticCalculator} at the time this class is
instantiated. The start coordinates
* given by {@link #φ1} and {@link #λ1} shall never change for this whole
builder lifetime. However the
* end coordinates ({@link #φ2}, {@link #λ2}) will vary at each step.
*/
- private final class PathBuilder extends Bezier {
+ private class PathBuilder extends Bezier {
/**
- * End point coordinates and derivative, together with geodesic and
loxodromic distances.
+ * The initial (i) and final (f) coordinates and derivatives, together
with geodesic and loxodromic distances.
* Saved for later restoration by {@link #reset()}.
*/
- private final double φf, λf, αf, distance, length;
+ private final double αi, αf, φf, λf, distance, length;
/**
* {@link #validity} flags at the time {@code PathBuilder} is
instantiated.
@@ -682,17 +724,19 @@ public class GeodeticCalculator {
private final double tolerance;
/**
- * Creates a builder for the given tolerance at equator in degrees.
+ * Creates a builder for the given tolerance at equator in metres.
*/
- PathBuilder(final double tolerance) {
+ PathBuilder(final double εx) {
super(ReferencingUtilities.getDimension(userToGeodetic.defaultCRS));
- this.tolerance = tolerance;
- φf = φ2;
- λf = λ2;
- αf = α2;
- distance = geodesicDistance;
- length = rhumblineLength;
- flags = validity;
+ αi = α1;
+ αf = α2;
+ φf = φ2;
+ λf = λ2;
+ α1 = IEEEremainder(α1, 2*PI); // For reliable
signum(α1) result.
+ tolerance = toDegrees(εx / radius);
+ distance = geodesicDistance;
+ length = rhumblineLength;
+ flags = validity;
}
/**
@@ -717,11 +761,20 @@ public class GeodeticCalculator {
α2 = αf;
} else {
geodesicDistance = distance * t;
+ validity |= GEODESIC_DISTANCE;
computeEndPoint();
}
- final double sign = signum(α1);
- if (sign == signum(λ1 - λ2)) {
- λ2 += 2*PI * sign; // We need λ₁ < λ₂ if heading
east, or λ₁ > λ₂ if heading west.
+ return evaluateAtEndPoint();
+ }
+
+ /**
+ * Implementation of {@link #evaluateAt(double)} using the current φ₂,
λ₂ and α₂ values.
+ * This method stores the projected coordinates in the {@link #point}
array and returns
+ * the derivative ∂y/∂x.
+ */
+ final double evaluateAtEndPoint() throws TransformException {
+ if ((λ2 - λ1) * α1 < 0) {
+ λ2 += 2*PI * signum(α1); // We need λ₁ < λ₂ if
heading east, or λ₁ > λ₂ if heading west.
}
final Matrix d = geographic(φ2, λ2).inverseTransform(point); //
Coordinates and Jacobian of point.
final double m00 = d.getElement(0,0);
@@ -736,20 +789,54 @@ public class GeodeticCalculator {
* α2 is azimuth angle in radians, with 0 pointing toward north
and values increasing clockwise.
* d is the Jacobian matrix from (φ,λ) to the user coordinate
reference system.
*/
- final double dx = cos(α2); // sin(π/2 - α) = -sin(α -
π/2) = cos(α)
- final double dy = sin(α2); // cos(π/2 - α) = +cos(α -
π/2) = sin(α)
+ final double dx = cos(α2); // sin(π/2 - α) = -sin(α -
π/2) = cos(α)
+ final double dy = sin(α2); // cos(π/2 - α) = +cos(α -
π/2) = sin(α)
final double tx = m00*dx + m01*dy;
final double ty = m10*dx + m11*dy;
return ty / tx;
}
/**
+ * Returns whether the point at given (<var>x</var>, <var>y</var>)
coordinates is close to the geodesic path.
+ * This method is invoked when the {@link Bezier} helper class thinks
that the point is not on the path, but
+ * could be wrong because of the difficulty to evaluate the Bézier
<var>t</var> parameter of closest point.
+ */
+ @Override
+ protected boolean isValid(final double x, final double y) throws
TransformException {
+ /*
+ * Following code is equivalent to `setEndPoint(new
DirectPosition2D(x, y))` but without checks
+ * for argument validity and with less temporary objects creation
(we recycle the `point` array).
+ */
+ point[0] = x;
+ point[1] = y;
+ for (int i=2; i<point.length; i++) {
+ point[i] = 0;
+ }
+ userToGeodetic.transform(point);
+ φ2 = toRadians(point[0]);
+ λ2 = toRadians(point[1]);
+ validity |= END_POINT;
+ /*
+ * Computes the azimuth to the given point. This azimuth may be
different than the azimuth of the
+ * geodesic path we are building. Compute the point that we would
have if the azimuth was correct
+ * and check the distance between those two points.
+ */
+ computeDistance();
+ α1 = αi;
+ computeEndPoint();
+ final DirectPosition p = geographic(φ2, λ2).inverseTransform();
+ return abs(p.getOrdinate(0) - x) <= εx &&
+ abs(p.getOrdinate(1) - y) <= εy;
+ }
+
+ /**
* Restores the enclosing {@link GeodeticCalculator} to the state that
it has at {@code PathBuilder} instantiation time.
*/
- void reset() {
+ final void reset() {
+ α1 = αi;
+ α2 = αf;
φ2 = φf;
λ2 = λf;
- α2 = αf;
geodesicDistance = distance;
rhumblineLength = length;
validity = flags;
@@ -757,6 +844,48 @@ public class GeodeticCalculator {
}
/**
+ * Builds a circular region around the start point. The shape is created
as a sequence of Bézier curves.
+ */
+ private final class CircularPath extends PathBuilder {
+ /**
+ * Creates a builder for the given tolerance at equator in degrees.
+ */
+ CircularPath(final double εx) {
+ super(εx);
+ }
+
+ /**
+ * Invoked for computing a new point on the circular path. This method
is invoked with a <var>t</var> value varying from
+ * 0 to 1 inclusive. The <var>t</var> value is multiplied by 2π for
getting an angle. This method stores the coordinates
+ * in the {@link #point} array and returns the derivative (∂y/∂x) at
that location.
+ *
+ * @param t angle fraction from 0 to 1 inclusive.
+ * @return derivative (∂y/∂x) at the point.
+ * @throws TransformException if the point coordinates can not be
computed.
+ */
+ @Override
+ protected double evaluateAt(final double t) throws TransformException {
+ α1 = IEEEremainder((t - 0.5) * (2*PI), 2*PI);
+ validity |= STARTING_AZIMUTH;
+ computeEndPoint();
+ final double d = evaluateAtEndPoint();
+ if (depth <= 1) {
+ // Force division of the curve in two smaller curves. We want
at least 4 Bézier curves in an ellipse.
+ εx = εy = -1;
+ }
+ return d;
+ }
+
+ /**
+ * No additional test (compared to {@link Bezier} base class) for
determining if the point is close enough.
+ */
+ @Override
+ protected boolean isValid(final double x, final double y) throws
TransformException {
+ return false;
+ }
+ }
+
+ /**
* Returns a string representation of start point, end point, azimuths and
distance.
* The text representation is implementation-specific and may change in
any future version.
* Current implementation is like below:
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
index 70751ec..c904168 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
@@ -31,6 +31,8 @@ import
org.apache.sis.internal.referencing.j2d.ShapeUtilitiesExt;
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.geometry.DirectPosition2D;
import org.apache.sis.util.CharSequences;
+import org.apache.sis.math.StatisticsFormat;
+import org.apache.sis.math.Statistics;
import org.apache.sis.measure.Units;
import org.apache.sis.test.widget.VisualCheck;
import org.apache.sis.test.OptionalTestData;
@@ -213,9 +215,27 @@ public final strictfp class GeodeticCalculatorTest extends
TestCase {
}
/**
- * Tests {@link GeodeticCalculator#toGeodesicPath2D(double)}. This method
uses a CRS that swap axis order
- * as a way to verify that user-specified CRS is taken in account. The
start point and end point are the
- * same than in {@link #testGeodesicDistanceAndAzimuths()}. Note that this
path crosses the anti-meridian,
+ * Tests {@link GeodeticCalculator#createCircularRegion2D(double)}.
+ *
+ * @throws TransformException if an error occurred while transforming
coordinates.
+ */
+ @Test
+ @DependsOnMethod("testUsingTransform")
+ public void testCircularRegion2D() throws TransformException {
+ final GeodeticCalculator c = create(true);
+ c.setStartPoint(-33.0, -71.6); // Valparaíso
+ c.setGeodesicDistance(100000); // 100 km
+ Shape region = c.createCircularRegion2D(10000);
+ if (VisualCheck.SHOW_WIDGET) {
+ VisualCheck.show(region);
+ }
+ // TODO: test bounding box.
+ }
+
+ /**
+ * Tests {@link GeodeticCalculator#createGeodesicPath2D(double)}. This
method uses a CRS that swap axis order
+ * as a way to verify that user-specified CRS is taken in account. The
start point and end point are the same
+ * than in {@link #testGeodesicDistanceAndAzimuths()}. Note that this path
crosses the anti-meridian,
* so the end point needs to be shifted by 360°.
*
* @throws TransformException if an error occurred while transforming
coordinates.
@@ -225,10 +245,10 @@ public final strictfp class GeodeticCalculatorTest
extends TestCase {
public void testGeodesicPath2D() throws TransformException {
final GeodeticCalculator c = create(true);
final double tolerance = 0.05;
- c.setStartPoint(-33.0, -71.6);
// Valparaíso
- c.setEndPoint ( 31.4, 121.8);
// Shanghai
- final Shape singleCurve = c.toGeodesicPath2D(Double.POSITIVE_INFINITY);
- final Shape multiCurves = c.toGeodesicPath2D(10000);
// 10 km tolerance.
+ c.setStartPoint(-33.0, -71.6);
// Valparaíso
+ c.setEndPoint ( 31.4, 121.8);
// Shanghai
+ final Shape singleCurve =
c.createGeodesicPath2D(Double.POSITIVE_INFINITY);
+ final Shape multiCurves = c.createGeodesicPath2D(10000);
// 10 km tolerance.
/*
* The approximation done by a single curve is not very good, but is
easier to test.
*/
@@ -257,7 +277,7 @@ public final strictfp class GeodeticCalculatorTest extends
TestCase {
c.setEndPoint (0, 12);
assertEquals(-90, c.getStartingAzimuth(), tolerance);
assertEquals(-90, c.getEndingAzimuth(), tolerance);
- final Shape geodeticCurve = c.toGeodesicPath2D(1);
+ final Shape geodeticCurve = c.createGeodesicPath2D(1);
final double[] coords = new double[2];
for (final PathIterator it = geodeticCurve.getPathIterator(null, 1);
!it.isDone(); it.next()) {
it.currentSegment(coords);
@@ -277,6 +297,7 @@ public final strictfp class GeodeticCalculatorTest extends
TestCase {
final Random random = TestUtilities.createRandomNumberGenerator();
final GeodeticCalculator c = create(false);
final Geodesic reference = new
Geodesic(c.ellipsoid.getSemiMajorAxis(), 1/c.ellipsoid.getInverseFlattening());
+ final StatisticsFormat sf = VERBOSE ? StatisticsFormat.getInstance() :
null;
for (int i=0; i<100; i++) {
final double φ1 = random.nextDouble() * 180 - 90;
final double λ1 = random.nextDouble() * 360 - 180;
@@ -292,38 +313,53 @@ public final strictfp class GeodeticCalculatorTest
extends TestCase {
assertEquals("Starting azimuth", expected.azi1,
c.getStartingAzimuth(), Formulas.ANGULAR_TOLERANCE);
assertEquals("Ending azimuth", expected.azi2,
c.getEndingAzimuth(), Formulas.ANGULAR_TOLERANCE);
assertTrue ("Rhumb ≧ geodesic", rhumbLine >= geodesic);
- if (false) {
- // Disabled because currently too inaccurate - see
https://issues.apache.org/jira/browse/SIS-453
- assertEquals("Distance measured along geodesic path",
geodesic, length(c), Formulas.ANGULAR_TOLERANCE);
+ if (sf != null) {
+ // Checks the geodesic path on only 10% of test data, because
this computation is expensive.
+ if ((i % 10) == 0) {
+ out.println(c);
+ out.println(sf.format(geodesicPathFitness(c, 1000)));
+ }
}
}
}
/**
- * Measures an estimation of the length of the path returned by {@link
GeodeticCalculator#toGeodesicPath2D(double)}.
- * This method iterates over line segments and use the given calculator
for computing the geodesic distance of each
- * segment. The state of the given calculator is modified by this method.
+ * Estimates the differences between the points on the Bézier curves and
the points computed by geodetic calculator.
+ * This method estimates the length of the path returned by {@link
GeodeticCalculator#createGeodesicPath2D(double)}
+ * and compares with the expected distance and azimuth at each point, by
iterating over line segments and computing
+ * the geodesic distance of each segment. The state of the given
calculator is modified by this method.
+ *
+ * @param resolution tolerance threshold for the curve approximation, in
metres.
+ * @return statistics about errors relative to the resolution.
*/
- private static double length(final GeodeticCalculator c) throws
TransformException {
- final PathIterator iterator =
c.toGeodesicPath2D(10).getPathIterator(null, 100);
- final double[] buffer = new double[2];
- double length=0;
+ private static Statistics[] geodesicPathFitness(final GeodeticCalculator
c, final double resolution) throws TransformException {
+ final PathIterator iterator =
c.createGeodesicPath2D(resolution).getPathIterator(null,
Formulas.ANGULAR_TOLERANCE);
+ final Statistics xError = new Statistics("Δx/r");
+ final Statistics yError = new Statistics("Δy/r");
+ final Statistics aErrors = new Statistics("Δα (°)");
+ final double azimuth = c.getStartingAzimuth();
+ final double toMetres = (PI/180) *
Formulas.getAuthalicRadius(c.ellipsoid);
+ final double[] buffer = new double[2];
while (!iterator.isDone()) {
switch (iterator.currentSegment(buffer)) {
- default: fail("Unexpected path"); break;
- case PathIterator.SEG_MOVETO: {
- c.setStartPoint(buffer[0], buffer[1]);
- break;
- }
+ default: fail("Unexpected segment"); break;
+ case PathIterator.SEG_MOVETO: break;
case PathIterator.SEG_LINETO: {
c.setEndPoint(buffer[0], buffer[1]);
- length += c.getGeodesicDistance();
- c.moveToEndPoint();
+ aErrors.accept(abs(c.getStartingAzimuth() - azimuth));
+ c.setStartingAzimuth(azimuth);
+ DirectPosition endPoint = c.getEndPoint();
+ final double φ = endPoint.getOrdinate(0);
+ final double λ = endPoint.getOrdinate(1);
+ double dy = (buffer[0] - φ) * toMetres;
+ double dx = IEEEremainder(buffer[1] - λ, 360) * toMetres *
cos(toRadians(φ));
+ yError.accept(abs(dy) / resolution);
+ xError.accept(abs(dx) / resolution);
}
}
iterator.next();
}
- return length;
+ return new Statistics[] {xError, yError, aErrors};
}
/**
diff --git a/ide-project/NetBeans/nbproject/genfiles.properties
b/ide-project/NetBeans/nbproject/genfiles.properties
index b11d782..89d141d 100644
--- a/ide-project/NetBeans/nbproject/genfiles.properties
+++ b/ide-project/NetBeans/nbproject/genfiles.properties
@@ -3,6 +3,6 @@
build.xml.data.CRC32=58e6b21c
build.xml.script.CRC32=462eaba0
[email protected]
-nbproject/build-impl.xml.data.CRC32=790389f0
+nbproject/build-impl.xml.data.CRC32=15917a2f
nbproject/build-impl.xml.script.CRC32=a7689f96
nbproject/[email protected]
diff --git a/ide-project/NetBeans/nbproject/project.xml
b/ide-project/NetBeans/nbproject/project.xml
index 233ef9b..c22b81f 100644
--- a/ide-project/NetBeans/nbproject/project.xml
+++ b/ide-project/NetBeans/nbproject/project.xml
@@ -86,6 +86,7 @@
<word>boolean</word>
<word>Bézier</word>
<word>centimetre</word>
+ <word>circlelike</word>
<word>classname</word>
<word>classnames</word>
<word>classpath</word>
