sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/02: Use Newton iterative method for finding if a given (xp,yp) point is near enough to the Bézier curve. This removes the need for https://issues.apache.org/jira/browse/SIS-455 and produces much better results than what we had before.
Date Sat, 25 May 2019 12:55:55 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

commit ce2eebd96610b5e561f81d616ee5cfa1640df4fc
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Sat May 25 14:54:18 2019 +0200

    Use Newton iterative method for finding if a given (xp,yp) point is near enough to the
Bézier curve.
    This removes the need for https://issues.apache.org/jira/browse/SIS-455 and produces much
better results than what we had before.
---
 .../sis/internal/referencing/j2d/Bezier.java       | 189 ++++++++++++++-------
 .../apache/sis/referencing/GeodeticCalculator.java |  58 ++-----
 .../sis/referencing/GeodeticCalculatorTest.java    |   8 +-
 .../org/apache/sis/test/widget/ShapeViewer.java    |   2 +-
 4 files changed, 148 insertions(+), 109 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 464f8a4..dcae453 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
@@ -52,7 +52,14 @@ public abstract class Bezier {
      *
      * @see #depth
      */
-    private static final int DEPTH_LIMIT = 10;
+    private static final int DEPTH_LIMIT = 12;
+
+    /**
+     * Maximum number of iterations for iterative computations. We use a small value because
failure to provide
+     * an answer with a small number of iterations probably means that the curve needs to
be separated again in
+     * two smaller curves.
+     */
+    private static final int MAXIMUM_ITERATIONS = 4;
 
     /**
      * The path where to append Bézier curves.
@@ -130,28 +137,6 @@ public abstract class Bezier {
     protected abstract void 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.
-     *
-     * @see <a href="https://issues.apache.org/jira/browse/SIS-455">SIS-455</a>
-     */
-    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(1)}. 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.
@@ -271,7 +256,7 @@ public abstract class Bezier {
                           double x2, double y2,
                           double x3, double y3,
                           final double  x4, final double  y4,
-                          final double dx4, final double dy4) throws TransformException
+                          final double dx4, final double dy4)
     {
         /*
          * Equations in this method are simplified as if (x₀,y₀) coordinates are (0,0).
@@ -365,9 +350,9 @@ public abstract class Bezier {
          * This allows `evaluateAt(t)` to return NaN if it can not provide a point for a
given `t` value.
          */
         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) {
+            double x = 27./64*ax + 9./64*bx + 1./64*Δx;
+            double y = 27./64*ay + 9./64*by + 1./64*Δy;
+            if (abs(x - x1) > εx || abs(y - 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
@@ -382,12 +367,7 @@ public abstract class Bezier {
                  *     (−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).
-                 * See https://issues.apache.org/jira/browse/SIS-455 for more information.
-                 *
-                 * Instead we perform another test using the tangent of the curve at point
P₁ (and later P₃).
+                 * fixed t values. 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).
                  *
@@ -395,43 +375,29 @@ public abstract class Bezier {
                  * 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;
+                if (!isNear(x,  y,  0.25,           // Point at t=¼
+                           ax, ay, 27./16,          // Coefficient for t=¼
+                           bx, by, 18./16,
+                           Δx, Δy,  3./16,
+                           x1, y1))                 // Point for which to evaluate the distance.
+                {
+                    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.
-                 */
             }
             /*
              * 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;
-                if (abs(xi - x3) > εx || abs(yi - y3) > εy) {
-                    if (!isValid(xi + x0, yi + y0)) return false;
+            x = 9./64*ax + 27./64*bx + 27./64*Δx;
+            y = 9./64*ay + 27./64*by + 27./64*Δy;
+            if (abs(x - x3) > εx || abs(y - y3) > εy) {
+                if (!isNear(x,  y,  0.75,
+                           ax, ay,  3./16,
+                           bx, by, 18./16,
+                           Δx, Δy, 27./16,
+                           x3, y3))
+                {
+                    return false;
                 }
             }
         }
@@ -466,4 +432,99 @@ public abstract class Bezier {
         }
         return true;
     }
+
+    /**
+     * Returns {@code true} if the point specified by the ({@code xp}, {@code yp}) coordinates
is close to the Bézier curve
+     * specified by the given control points. All given control points shall be shifted as
if (x₀,y₀) coordinates are (0,0).
+     * If this method can not determine for sure that the given point is close, it conservatively
returns {@code false}.
+     * This is okay since failure to determine whether the point is close may mean that the
curvature is strong, in which
+     * case we should probably split the curves in two smaller curves again.
+     *
+     * @param  ax  the <var>x</var> coordinate of first control point.
+     * @param  ay  the <var>y</var> coordinate of first control point.
+     * @param  bx  the <var>x</var> coordinate of mid-point.
+     * @param  by  the <var>y</var> coordinate of mid-point.
+     * @param  x4  the <var>x</var> coordinate of last control point.
+     * @param  y4  the <var>y</var> coordinate of last control point.
+     * @param  xp  the <var>x</var> coordinate of the point to test for proximity.
+     * @param  yp  the <var>y</var> coordinate of the point to test for proximity.
+     * @param   x  the <var>x</var> coordinate of a close point on the curve.
+     * @param   y  the <var>y</var> coordinate of a close point on the curve.
+     * @param   t  the <var>t</var> parameter of Bézier curve where (<var>x</var>,<var>y</var>),
and <var>m</var> are evaluated.
+     * @param  am  the 3(1−t)²      multiplication factor of {@code ax} and {@code ay}
coordinates for given <var>t</var> value.
+     * @param  bm  the 6t(1−t)      multiplication factor of {@code bx} and {@code by}
coordinates for given <var>t</var> value.
+     * @param  cm  the 3t²          multiplication factor of {@code x4} and {@code y4} coordinates
for given <var>t</var> value.
+     * @return {@code true} if the given point is close to the Bézier curve, or {@code false}
otherwise or in case of doubt.
+     *
+     * @see <a href="https://issues.apache.org/jira/browse/SIS-455">SIS-455</a>
+     */
+    private boolean isNear(double x, double y, double t,
+            final double ax, final double ay, double am,
+            final double bx, final double by, double bm,
+            final double x4, final double y4, double cm,
+            final double xp, final double yp)
+    {
+        final double Δax = bx - ax;
+        final double Δay = by - ay;
+        final double Δbx = x4 - bx;
+        final double Δby = y4 - by;
+        int n = MAXIMUM_ITERATIONS;
+        for (;;) {
+            /*
+             * Tangent of the curve:
+             *
+             *     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 t.
+             * Then we compute the point on that line which is nearest to the given point
(xp,yp). It should be
+             * close to current (x,y) coordinates, but may vary a little bit.
+             */
+            final double tx = am * ax  +  bm * Δax  +  cm * Δbx;            // ∂x/∂t
+            final double ty = am * ay  +  bm * Δay  +  cm * Δby;            // ∂y/∂t
+            double xi, yi;
+            if (abs(tx) >= abs(ty)) {
+                final double slope  = ty / tx;                              // ∂y/∂x
at t=¼.
+                final double offset = (y - slope*x);                        // Value of y
at x=0.
+                xi = ((yp - offset) * slope + xp) / (slope*slope + 1);      // Move closer
to (xp,yp).
+                t += (xi - x) / tx;                                         // Using ∆x≈(∂x/∂t)⋅∆t
+            } else {
+                final double slope  = tx / ty;
+                final double offset = (x - slope*y);
+                yi = ((xp - offset) * slope + yp) / (slope*slope + 1);
+                t += (yi - y) / ty;
+            }
+            if      (t < 0) t = 0;
+            else if (t > 1) t = 1;
+            /*
+             * At this point we got the xi or yi coordinate or the closest point on the linear
approximation,
+             * and an estimation of the t Bézier parameter for that point. Use that t value
for finding the
+             * (xi,yi) coordinates on the Bézier curve.
+             */
+            am = 1 - t;                                             // New coefficients for
the x(t) and y(t).
+            bm = 3*am;
+            cm = t*t;
+            xi = ((ax*am + t*bx)*bm + x4*cm)*t;                     // 3(1−t)²t⋅aₓ
+ 3(1−t)t²⋅bₓ + t³⋅x₄
+            yi = ((ay*am + t*by)*bm + y4*cm)*t;
+            if (abs(xi - xp) <= εx && abs(yi - yp) <= εy) {
+                return true;                                        // Point is close enough.
+            }
+            if (--n >= 0 && (abs(xi - x) > εx/2 || abs(yi - y) > εy/2))
{
+                /*
+                 * We need one more iteration. Change the coefficients for making them apply
+                 * to the x′(t) and y′(t) derivatives instead of the x(t) and y(t) functions.
+                 */
+                am *= bm;           // 3(1−t)²
+                bm *= 2*t;          // 6t(1−t)
+                cm *= 3;            // 3t²
+                x = xi;
+                y = yi;
+            } else {
+                /*
+                 * Too many iterations, or a coordinate value is NaN, or we find the nearest
point
+                 * but that point is too far from the given (xp,yp) coordinates.
+                 */
+                return false;
+            }
+        }
+    }
 }
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 8bbcbae..db51946 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
@@ -719,10 +719,10 @@ public class GeodeticCalculator {
      */
     private class PathBuilder extends Bezier {
         /**
-         * The initial (i) and final (f) coordinates and derivatives, together with geodesic
and loxodromic distances.
+         * The final (f) coordinates and derivatives, together with geodesic and loxodromic
distances.
          * Saved for later restoration by {@link #reset()}.
          */
-        private final double dφi, dλi, dφf, dλf, φf, λf, distance, length;
+        private final double dφf, dλf, φf, λf, distance, length;
 
         /**
          * {@link #validity} flags at the time {@code PathBuilder} is instantiated.
@@ -740,7 +740,6 @@ public class GeodeticCalculator {
          */
         PathBuilder(final double εx) {
             super(ReferencingUtilities.getDimension(userToGeodetic.defaultCRS));
-            dφi = dφ1;  dλi = dλ1;
             dφf = dφ2;  dλf = dλ2;  φf = φ2;  λf = λ2;
             tolerance = toDegrees(εx / radius);
             distance  = geodesicDistance;
@@ -798,41 +797,11 @@ public class GeodeticCalculator {
         }
 
         /**
-         * 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.
-         *
-         * @see <a href="https://issues.apache.org/jira/browse/SIS-455">SIS-455</a>
-         */
-        @Override
-        protected boolean isValid(final double x, final double y) throws TransformException
{
-            point[0] = x;
-            point[1] = y;
-            for (int i=2; i<point.length; i++) {
-                point[i] = 0;
-            }
-            userToGeodetic.transform(point);
-            setEndGeographicPoint(point[0], point[1]);
-            /*
-             * 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();
-            dφ1 = dφi;
-            dλ1 = dλ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.
          */
-        final void reset() {
-            dφ1 = dφi;  dφ2 = dφf;  φ2 = φf;
-            dλ1 = dλi;  dλ2 = dλf;  λ2 = λf;
+        void reset() {
+            dφ2 = dφf;  φ2 = φf;
+            dλ2 = dλf;  λ2 = λf;
             geodesicDistance = distance;
             rhumblineLength  = length;
             validity         = flags;
@@ -844,10 +813,17 @@ public class GeodeticCalculator {
      */
     private final class CircularPath extends PathBuilder {
         /**
+         * The initial (i) derivatives, saved for later restoration by {@link #reset()}.
+         */
+        private final double dφi, dλi;
+
+        /**
          * Creates a builder for the given tolerance at equator in degrees.
          */
         CircularPath(final double εx) {
             super(εx);
+            dφi = dφ1;
+            dλi = dλ1;
         }
 
         /**
@@ -860,7 +836,7 @@ public class GeodeticCalculator {
          */
         @Override
         protected void evaluateAt(final double t) throws TransformException {
-            double α1 = IEEEremainder((t - 0.5) * (2*PI), 2*PI);
+            final double α1 = t * (2*PI) + PI/4;        // Add 45° for hiding little corners
in multiple of 90°.
             dλ1 = sin(α1);
             dφ1 = cos(α1);
             validity |= STARTING_AZIMUTH;
@@ -876,11 +852,13 @@ public class GeodeticCalculator {
         }
 
         /**
-         * No additional test (compared to {@link Bezier} base class) for determining if
the point is close enough.
+         * Restores the enclosing {@link GeodeticCalculator} to the state that it has at
{@code PathBuilder} instantiation time.
          */
         @Override
-        protected boolean isValid(final double x, final double y) throws TransformException
{
-            return false;
+        void reset() {
+            dφ1 = dφi;
+            dλ1 = dλi;
+            super.reset();
         }
     }
 
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 43b871d..c846cea 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
@@ -230,10 +230,10 @@ public final strictfp class GeodeticCalculatorTest extends TestCase
{
             VisualCheck.show(region);
         }
         final Rectangle2D bounds = region.getBounds2D();
-        assertEquals("xmin", -72.67228, bounds.getMinX(), 5E-6);
-        assertEquals("ymin", -33.89932, bounds.getMinY(), 5E-6);
-        assertEquals("xmax", -70.52772, bounds.getMaxX(), 5E-6);
-        assertEquals("ymax", -32.10068, bounds.getMaxY(), 5E-6);
+        assertEquals("x",    -71.6, bounds.getCenterX(), 1E-3);
+        assertEquals("y",    -33.0, bounds.getCenterY(), 1E-3);
+        assertEquals("width",  2.8, bounds.getWidth(),   0.1);      // Would have expected
2.1 if box was tight.
+        assertEquals("height", 2.0, bounds.getHeight(),  0.1);      // Would have expected
1.8 if box was tight.
     }
 
     /**
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/widget/ShapeViewer.java
b/core/sis-referencing/src/test/java/org/apache/sis/test/widget/ShapeViewer.java
index dc471ea..e93e36d 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/test/widget/ShapeViewer.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/test/widget/ShapeViewer.java
@@ -50,7 +50,7 @@ final strictfp class ShapeViewer extends JPanel {
      * Colors to use for drawing shapes.
      */
     private final Color[] colors = {
-        Color.RED, Color.GREEN, Color.BLUE
+        Color.RED, Color.ORANGE, Color.YELLOW, Color.GREEN, Color.CYAN, Color.BLUE, Color.MAGENTA
     };
 
     /**


Mime
View raw message