sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Try harder to avoid unreasonably long chain of Bézier curves. First implementation of createCircularRegion2D(…), not yet working because of infinite derivative values.
Date Wed, 22 May 2019 18:49:15 GMT
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 <martin.desruisseaux@geomatys.com>
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
 build.xml.stylesheet.CRC32=28e38971@1.53.1.46
-nbproject/build-impl.xml.data.CRC32=790389f0
+nbproject/build-impl.xml.data.CRC32=15917a2f
 nbproject/build-impl.xml.script.CRC32=a7689f96
 nbproject/build-impl.xml.stylesheet.CRC32=3a2fa800@1.89.1.48
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>


Mime
View raw message