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: Better map for determining if a cubic Bézier curve is an acceptable approximation of the geodesic path. Previous implementation assumed that point at t=1/4 on the Bézier curve was located at 1/4 of geodesic path. This is not true, because the t parameter in Bézier equations does not vary at the same speed than distance. The relationship between `t` and distance is too complicated for the needs of this method (it requires numerical integration and iterative method since there is no exact [...]
Date Tue, 21 May 2019 22:44:32 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 3d2a248  Better map for determining if a cubic Bézier curve is an acceptable approximation
of the geodesic path. Previous implementation assumed that point at t=1/4 on the Bézier curve
was located at 1/4 of geodesic path. This is not true, because the t parameter in Bézier
equations does not vary at the same speed than distance. The relationship between `t` and
distance is too complicated for the needs of this method (it requires numerical integration
and iterative method since t [...]
3d2a248 is described below

commit 3d2a248d9f1bc23178d3c050e2d761847585d17d
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed May 22 00:32:36 2019 +0200

    Better map for determining if a cubic Bézier curve is an acceptable approximation of
the geodesic path.
    Previous implementation assumed that point at t=1/4 on the Bézier curve was located at
1/4 of geodesic path.
    This is not true, because the t parameter in Bézier equations does not vary at the same
speed than distance.
    The relationship between `t` and distance is too complicated for the needs of this method
    (it requires numerical integration and iterative method since there is no exact solution).
    They are close, but not close enough for ignoring the difference. This commit adds the
use of
    tangent in calculation that determine whether a point is considered close to the Bézier
curve.
    This reduces the number of Bézier curves by 1/3 in the test performed by GeodeticCalculatorTest.
---
 .../sis/internal/referencing/j2d/Bezier.java       | 85 ++++++++++++++++++----
 .../apache/sis/referencing/GeodeticCalculator.java |  1 +
 2 files changed, 72 insertions(+), 14 deletions(-)

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 90f3999..9ea2294 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
@@ -170,7 +170,7 @@ public abstract class Bezier {
      * The (x₂,y₂) arguments give the coordinates of the point at <var>t</var>=½
(2/4).
      * The (x₄,y₄) arguments give the coordinates of point P₄ at <var>t</var>=1
(4/4).
      *
-     * <p>The P₁ and P₃ points are for quality control. They shall be points at
<var>t</var>=¼ <var>t</var>=¾ respectively.
+     * <p>The P₁ and P₃ points are for quality control. They should be points at
<var>t</var>≈¼ <var>t</var>≈¾ respectively.
      * Those points are not used for computing the curve, but are used for checking if the
curve is an accurate approximation.
      * 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>
@@ -179,21 +179,21 @@ public abstract class Bezier {
      * 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>
      *
-     * @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).
+     * @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).
      * @param  x2  <var>x</var> coordinate at <var>t</var>=½ (mid-point).
      * @param  y2  <var>y</var> coordinate at <var>t</var>=½ (mid-point).
-     * @param  x3  <var>x</var> coordinate at <var>t</var>=¾ (quality
control, may be NaN).
-     * @param  y3  <var>y</var> coordinate at <var>t</var>=¾ (quality
control, may be NaN).
+     * @param  x3  <var>x</var> coordinate at <var>t</var>≈¾ (quality
control, may be NaN).
+     * @param  y3  <var>y</var> coordinate at <var>t</var>≈¾ (quality
control, may be NaN).
      * @param  x4  <var>x</var> coordinate at <var>t</var>=1 (end
point).
      * @param  y4  <var>y</var> coordinate at <var>t</var>=1 (end
point).
      * @param  α4  the derivative (∂y/∂x) at ending point.
      * @return {@code true} if the curve has been added to the path, or {@code false} if
the approximation is not sufficient.
      */
-    final boolean curve(double x1, double y1,
-                        double x2, double y2,
-                        double x3, double y3,
-                        final double x4, final double y4, final double α4)
+    private boolean curve(double x1, double y1,
+                          double x2, double y2,
+                          double x3, double y3,
+                          final double x4, final double y4, final double α4)
     {
         /*
          * Equations in this method are simplified as if (x₀,y₀) coordinates are (0,0).
@@ -279,12 +279,69 @@ public abstract class Bezier {
          *     P(¾) = (   P₀ +  9⋅A + 27⋅B + 27⋅P₄)/64
          *
          * If any of (x₁,y₁) and (x₃,y₃) coordinates is NaN, then this method accepts
the curve as valid.
-         * This allows `evaluateAt(t)` method to return NaN if it can not provide a point
for a given `t` value.
+         * This allows `evaluateAt(t)` to return NaN if it can not provide a point for a
given `t` value.
          */
-        if (abs(27./64*ax + 9./64*bx + 1./64*Δx - x1) > εx || abs(9./64*ax + 27./64*bx
+ 27./64*Δx - x3) > εx ||
-            abs(27./64*ay + 9./64*by + 1./64*Δy - y1) > εy || abs(9./64*ay + 27./64*by
+ 27./64*Δy - y3) > εy)
-        {
-            return false;
+        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₀) + 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.
+             */
+            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;
+            }
         }
         /*
          * Verify if we can simplify cubic curve to a quadratic curve. If we were elevating
the Bézier curve degree from
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 1b88c27..9bce0d0 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
@@ -639,6 +639,7 @@ public class GeodeticCalculator {
      * @throws IllegalStateException if some required properties have not been specified.
      */
     public Shape toGeodesicPath2D(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)) {
                 computeEndPoint();


Mime
View raw message