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: Replace α₀ and α₄ ratios by (dx₀,dy₀) and (dx₄,dy₄) vector components in order to avoid divison by zero.
Date Thu, 23 May 2019 12:31:31 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 0f797b8  Replace α₀ and α₄ ratios by (dx₀,dy₀) and (dx₄,dy₄) vector
components in order to avoid divison by zero.
0f797b8 is described below

commit 0f797b886acfc441dfa6ad6d2af6962771239a8d
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Thu May 23 14:26:29 2019 +0200

    Replace α₀ and α₄ ratios by (dx₀,dy₀) and (dx₄,dy₄) vector components in
order to avoid divison by zero.
---
 .../sis/internal/referencing/j2d/Bezier.java       | 222 ++++++++++++---------
 .../apache/sis/referencing/GeodeticCalculator.java |  32 ++-
 .../referencing/j2d/ShapeUtilitiesExt.java         |   5 +-
 3 files changed, 142 insertions(+), 117 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 1c9d580..78694c5 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
@@ -23,7 +23,6 @@ import java.awt.geom.CubicCurve2D;
 import org.opengis.referencing.operation.TransformException;
 
 import static java.lang.Math.abs;
-import static java.lang.Double.isInfinite;
 
 
 /**
@@ -61,17 +60,30 @@ public abstract class Bezier {
     private final Path2D path;
 
     /**
-     * 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.
+     * Maximal distance (approximate) on <var>x</var> and <var>y</var>
axis between Bézier curves and desired curve.
+     * Also used for deciding if a cubic Bézier curve can be simplified to quadratic curve
or straight line. 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
negative value for forcing
+     * unconditional divisions of Bézier curves in two sub-curves.
      */
     protected double εx, εy;
 
     /**
-     * (<var>x</var> and <var>y</var>) coordinates and (∂y/∂x)
derivative of starting point at <var>t</var>=0.
+     * (<var>x</var> and <var>y</var>) coordinates at starting point
<var>t</var>=0.
+     * This is initialized by {@link #build()} and updated by {@link #sequence} when moving
along the path.
      */
-    private double x0, y0, α0;
+    private double x0, y0;
+
+    /**
+     * Components of α=(∂y/∂x) derivative at starting point <var>t</var>=0.
+     * Initialized and updated in same time than {@link #x0} and {@link #y0}.
+     */
+    private double dx0, dy0;
+
+    /**
+     * Components of α=(∂y/∂x) derivative at the point evaluated by {@link #evaluateAt(double)}.
+     */
+    protected double dx, dy;
 
     /**
      * A buffer used by subclasses for storing results of {@link #evaluateAt(double)}.
@@ -104,19 +116,18 @@ public abstract class Bezier {
      * This method will also be invoked at least for values ¼ and ¾, and potentially for
other values too.
      *
      * <p>This method shall store the point coordinates in the {@link #point} array
with <var>x</var> coordinate in the
-     * first element and <var>y</var> coordinate in the second element. The returned
value is the derivative (∂y/∂x) at
-     * that location. If this method can not compute a coordinate, it can store {@link Double#NaN}
values except for
-     * <var>t</var>=0, ½ and 1 where finite coordinates are mandatory.</p>
+     * first element and <var>y</var> coordinate in the second element. This
method shall also store derivative (∂y/∂x)
+     * at that location in the {@link #dx} and {@link #dy} fields. If this method can not
compute a coordinate, it can
+     * store {@link Double#NaN} values except for <var>t</var>=0, ½ and 1 where
finite coordinates are mandatory.</p>
      *
      * <p>Subclasses can optionally update the {@link #εx} and {@link #εy} values
if the tolerance thresholds change as a
      * result of this method call, for example because we come closer to a pole. The tolerance
values used for each Bézier
      * curve are the ones computed at <var>t</var>=¼ and <var>t</var>=¾
of that curve.</p>
      *
      * @param  t  desired point on the curve, from 0 (start point) to 1 (end point) inclusive.
-     * @return derivative (∂y/∂x) at the point.
      * @throws TransformException if the point coordinates can not be computed.
      */
-    protected abstract double evaluateAt(double t) throws TransformException;
+    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
@@ -140,22 +151,26 @@ public abstract class Bezier {
 
     /**
      * 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
+     * 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.
      *
      * @return the sequence of Bézier curves.
      * @throws TransformException if the coordinates of a point can not be computed.
      */
     public final Path2D build()  throws TransformException {
-        α0 = evaluateAt(0);
-        x0 = point[0];
-        y0 = point[1];
+        evaluateAt(0);
+        x0  = point[0];
+        y0  = point[1];
+        dx0 = dx;
+        dy0 = dy;
         path.moveTo(x0, y0);
-        final double α2 = evaluateAt(0.5);
-        final double x2 = point[0];
-        final double y2 = point[1];
-        final double α4 = evaluateAt(1);                                // `point` become
(x₄,y₄) with this method call.
-        sequence(0, 1, x2, y2, α2, point[0], point[1], α4);             // Must be after
the call to `evaluateAt(t₄)`.
+        evaluateAt(0.5);
+        final double  x2 = point[0];
+        final double  y2 = point[1];
+        final double dx2 = dx;
+        final double dy2 = dy;
+        evaluateAt(1);                                                  // `point` become
(x₄,y₄) with this method call.
+        sequence(0, 1, x2, y2, dx2, dy2, point[0], point[1], dx, dy);   // Must be after
the call to `evaluateAt(t₄)`.
         return path;
     }
 
@@ -164,45 +179,52 @@ public abstract class Bezier {
      * {@code evaluateAt(tmax)}. This method invokes itself recursively as long as more points
are needed for achieving
      * the precision requested by the {@code εx} and {@code εy} parameters given at construction
time.
      *
-     * <p>The {@link #x0}, {@link #y0} and {@link #α0} fields must be set to the start
point before this method is invoked.
-     * On method completion those fields will be updated to the coordinates and derivative
of end point, thus enabling
-     * chaining with the next {@code sequence(…)} call.</p>
+     * <p>The {@link #x0}, {@link #y0}, {@link #dx0} and {@link #dy0} fields must be
set to the start point before this
+     * method is invoked. On method completion those fields will be updated to the coordinates
and derivative of end point,
+     * thus enabling chaining with the next {@code sequence(…)} call.</p>
      *
      * @param  tmin  <var>t</var> value of the start point (initially 0).
      * @param  tmax  <var>t</var> value of the end point   (initially 1).
      * @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  α2    the derivative (∂y/∂x) at mid-point.
+     * @param  dx2   <var>x</var> component of the the derivative (∂y/∂x)
at mid-point.
+     * @param  dy2   <var>y</var> component of the the derivative (∂y/∂x)
at mid-point.
      * @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 end point.
+     * @param  dx4   <var>x</var> component of the derivative (∂y/∂x) at
end point.
+     * @param  dy4   <var>y</var> component of the derivative (∂y/∂x) at
end point.
      * @throws TransformException if the coordinates of a point can not be computed.
      */
     private void sequence(final double tmin, final double tmax,
-            final double x2, final double y2, final double α2,
-            final double x4, final double y4, final double α4) throws TransformException
+            final double x2, final double y2, final double dx2, final double dy2,
+            final double x4, final double y4, final double dx4, final double dy4) throws
TransformException
     {
-        final double α3 = evaluateAt(0.25*tmin + 0.75*tmax);
-        final double x3 = point[0];
-        final double y3 = point[1];
-        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.
+        evaluateAt(0.25*tmin + 0.75*tmax);
+        final double  x3 = point[0];
+        final double  y3 = point[1];
+        final double dx3 = dx;
+        final double dy3 = dy;
+        final double  tx = εx;
+        final double  ty = εy;
+        evaluateAt(0.75*tmin + 0.25*tmax);                          // `point` become (x₁,y₁)
with this method call.
+        final double  x1 = point[0];                                // Must be after the
call to `evaluateAt(t₁)`.
+        final double  y1 = point[1];
+        final double dx1 = dx;
+        final double dy1 = dy;
+        if (tx < εx) εx = tx;                                       // Take smallest
tolerance values.
         if (ty < εy) εy = ty;
-        if (curve(xi, yi, x2, y2, x3, y3, x4, y4, α4)) {
-            x0 = x4;
-            y0 = y4;
-            α0 = α4;
+        if (curve(x1, y1, x2, y2, x3, y3, x4, y4, dx4, dy4)) {
+            x0  = x4;
+            y0  = y4;
+            dx0 = dx4;
+            dy0 = dy4;
         } else {
             depth++;
             final double εxo = εx;
             final double εyo = εy;
             final double th = 0.5 * (tmin + tmax);
-            sequence(tmin, th, xi, yi, α1, x2, y2, α2);
-            sequence(th, tmax, x3, y3, α3, x4, y4, α4);
+            sequence(tmin, th, x1, y1, dx1, dy1, x2, y2, dx2, dy2);
+            sequence(th, tmax, x3, y3, dx3, dy3, x4, y4, dx4, dy4);
             εx = εxo;
             εy = εyo;
             depth--;
@@ -213,7 +235,7 @@ public abstract class Bezier {
      * Computes a Bézier curve passing by the given points and with the given derivatives
at end points.
      * The cubic curve equation is:
      *
-     * <blockquote>B(t) = (1-t)³⋅P₀ + 3(1-t)²t⋅A + 3(1-t)t²⋅B + t³⋅P₄</blockquote>
+     * <blockquote>B(t) = (1−t)³⋅P₀ + 3(1−t)²t⋅A + 3(1−t)t²⋅B + t³⋅P₄</blockquote>
      *
      * where t ∈ [0…1], P₀ and P₄ are start and end points of the curve and A and
B are control points generally not on the curve.
      * For any point <var>P<sub>i</sub></var>, the index <var>i</var>
gives the value of the parameter <var>t</var> where the point
@@ -231,22 +253,23 @@ public abstract class Bezier {
      * 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).
-     * @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  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.
+     * @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  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  dx4  <var>x</var> component of the derivative α₄=∂y/∂x
at ending point.
+     * @param  dy4  <var>y</var> component of 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.
      */
     private boolean curve(double x1, double y1,
                           double x2, double y2,
                           double x3, double y3,
-                          final double x4, final double y4,
-                          final double α4) throws TransformException
+                          final double  x4, final double  y4,
+                          final double dx4, final double dy4) throws TransformException
     {
         /*
          * Equations in this method are simplified as if (x₀,y₀) coordinates are (0,0).
@@ -258,41 +281,34 @@ public abstract class Bezier {
         y1 -= y0;  y2 -= y0;  y3 -= y0;
         final double Δx = x4 - x0;
         final double Δy = y4 - y0;
-        final boolean isVertical = abs(Δx) <= εx;
-        if (isVertical || ((isInfinite(α0) || abs(Δx*α0 - Δy) <= εy)
-                       &&  (isInfinite(α4) || abs(Δx*α4 - Δy) <= εy)))
+        if ((abs(dx0) >= abs(dy0) ? abs(Δx*(dy0/dx0) - Δy) <= εy : abs(Δy*(dx0/dy0)
- Δx) <= εx) &&
+            (abs(dx4) >= abs(dy4) ? abs(Δx*(dy4/dx4) - Δy) <= εy : abs(Δy*(dx4/dy4)
- Δx) <= εx))
         {
             /*
-             * Following tests are partially redundant with above tests, but may detect a
larger
-             * error than above tests did. They are also necessary if a derivative was infinite,
-             * or if `isVertical` was true.
+             * Verify that all points are on the line joining P₀ to P₄. If any coordinate
tested below is NaN,
+             * ignore it since 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 boolean isHorizontal = abs(Δy) <= εy;
-            if (isHorizontal || ((α0 == 0 || abs(Δy/α0 - Δx) <= εx)
-                             &&  (α4 == 0 || abs(Δy/α4 - Δx) <= εx)))
-            {
-                /*
-                 * Verify that all points are on the line joining P₀ to P₄. We test P₂
first since it is the most likely
-                 * to be away if the curve is not almost a line. If any coordinate tested
below is NaN, ignore it since
-                 * 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.
-                 */
-                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)))
-                    {
+            if (depth < DEPTH_LIMIT) {
+                if (abs(Δx) >= abs(Δy)) {
+                    final double m = Δy/Δx;
+                    if (abs(m*x2 - y2) > εy || abs(m*x1 - y1) > εy || abs(m*x3 -
y3) > εy) {
+                        return false;
+                    }
+                } else {
+                    final double m = Δx/Δy;
+                    if (abs(m*y2 - x2) > εx || abs(m*y1 - x1) > εx || abs(m*y3 -
x3) > εx) {
                         return false;
                     }
                 }
-                path.lineTo(x4, y4);
-                return true;
             }
+            path.lineTo(x4, y4);
+            return true;
         }
         /*
          * Bezier curve equation for starting point P₀, ending point P₄ and control points
A and B
          *
-         * t ∈ [0…1]:  B(t)  =  (1-t)³⋅P₀ + 3(1-t)²t⋅A + 3(1-t)t²⋅B + t³⋅P₄
+         * t ∈ [0…1]:  B(t)  =  (1−t)³⋅P₀ + 3(1−t)²t⋅A + 3(1−t)t²⋅B
+ t³⋅P₄
          * Midpoint:   B(½)  =  ⅛⋅(P₀ + P₄) + ⅜⋅(A + B)
          *
          * Notation:   (x₀, y₀)   are coordinates of P₀ (same rule for P₄).
@@ -303,28 +319,38 @@ public abstract class Bezier {
          * Some relationships:
          *
          *     x₂ = ⅛⋅(x₀ + x₄) + ⅜⋅(Ax + Bx)
-         *     (Ay - y₀) / (Ax - x₀) = α₀
-         *     (y₄ - By) / (x₄ - Bx) = α₄
+         *     (Ay − y₀) / (Ax − x₀) = α₀
+         *     (y₄ − By) / (x₄ − Bx) = α₄
          *
          * Setting (x₀,y₀) = (0,0) for simplicity and rearranging above equations:
          *
-         *     Ax = (8⋅x₂ - x₄)/3 - Bx              where    Bx = x₄ - (y₄ - By)/α₄
-         *        = (8⋅x₂ - 4⋅x₄)/3 + (y₄ - By)/α₄
+         *     Ax = (8⋅x₂ − x₄)/3 − Bx              where    Bx = x₄ − (y₄
− By)/α₄
+         *        = (8⋅x₂ − 4⋅x₄)/3 + (y₄ − By)/α₄
          *
          * Doing similar rearrangement for y:
          *
-         *     By = (8⋅y₂ - y₄)/3 - Ay    where    Ay = Ax⋅α₀
-         *        = (8⋅y₂ - y₄)/3 - Ax⋅α₀
+         *     By = (8⋅y₂ − y₄)/3 − Ay    where    Ay = Ax⋅α₀
+         *        = (8⋅y₂ − y₄)/3 − Ax⋅α₀
          *
          * Putting together and isolating Ax:
          *
-         *      Ax = (8⋅x₂ - 4⋅x₄)/3 + (Ax⋅α₀ - (8⋅y₂ - 4⋅y₄)/3)/α₄
-         *         = (8⋅x₂ - 4⋅x₄ - (8⋅y₂ - 4⋅y₄)/α₄) / 3(1 - α₀/α₄)
+         *      Ax = (8⋅x₂ − 4⋅x₄)/3 + (Ax⋅α₀ − (8⋅y₂ − 4⋅y₄)/3)/α₄
+         *         = (8⋅x₂ − 4⋅x₄ − (8⋅y₂ − 4⋅y₄)/α₄) / 3(1 −
α₀/α₄)
+         *
+         * We substitute α₀=dy₀/dx₀ and α₄=dy₄/dx₄ and rearrange the dx and
dy terms in order
+         * to avoid division by zero. The results are the equations in the code below.
          */
-        final double ax = ((8*x2 - 4*Δx)*α4 - (8*y2 - 4*Δy)) / (3*(α4 - α0));
-        final double ay = ax * α0;
-        final double bx = (8*x2 - Δx)/3 - ax;
-        final double by = Δy - (Δx - bx)*α4;
+        final double m  = ((8*x2 - 4*Δx)*dy4 - (8*y2 - 4*Δy)*dx4) / (3*(dx0*dy4 - dy0*dx4));
+        final double ax = m * dx0;
+        final double ay = m * dy0;
+        final double bx, by;
+        if (abs(dx4) >= abs(dy4)) {
+            bx = (8*x2 - Δx)/3 - ax;
+            by = Δy - (Δx - bx)*(dy4/dx4);
+        } else {
+            by = (8*y2 - Δy)/3 - ay;
+            bx = Δx - (Δy - by)*(dx4/dy4);
+        }
         /*
          * At this point we got the control points A and B. Verify if the curve is a good
approximation.
          * We compute the points on the curve at t=¼ and t=¾ as below (simplified with
P₀ = (0,0)) and
@@ -347,11 +373,11 @@ public abstract class Bezier {
                  * 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³
+                 *     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
+                 *     (−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
@@ -360,7 +386,7 @@ public abstract class Bezier {
                  *
                  * 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).
+                 *     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
@@ -380,7 +406,7 @@ public abstract class Bezier {
                  * 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).
+                 *     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″
                  *
@@ -415,12 +441,12 @@ public abstract class Bezier {
          *
          * We want C instead, which can be computed in two ways:
          *
-         *     C   =  (3A - P₀)/2
-         *     C   =  (3B - P₄)/2
+         *     C   =  (3A − P₀)/2
+         *     C   =  (3B − P₄)/2
          *
          * We compute C both ways and check if they are close enough to each other:
          *
-         *     ΔC  =  (3⋅(B - A) - (P₄ - P₀))/2
+         *     Δ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 ½.
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 2ab39b2..e7744ab 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
@@ -741,16 +741,14 @@ public class GeodeticCalculator {
 
         /**
          * Invoked for computing a new point on the Bézier curve. This method is invoked
with a <var>t</var> value varying from
-         * 0 (start point) to 1 (end point) inclusive. This method stores the point coordinates
in the {@link #point} array with
-         * <var>x</var> coordinate in the first element and <var>y</var>
coordinate in the second element. The returned value is
-         * the derivative (∂y/∂x) at that location.
+         * 0 (start point) to 1 (end point) inclusive. This method stores the coordinates
in the {@link #point} array and the
+         * derivative (∂y/∂x) in the {@linkplain #dx dx} and {@linkplain #dy dy} fields.
          *
          * @param  t  desired point on the curve, from 0 (start point) to 1 (end point) 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 {
+        protected void evaluateAt(final double t) throws TransformException {
             if (t == 0) {
                 φ2 = φ1;                        // Start point requested.
                 λ2 = λ1;
@@ -764,7 +762,7 @@ public class GeodeticCalculator {
                 validity |= GEODESIC_DISTANCE;
                 computeEndPoint();
             }
-            return evaluateAtEndPoint();
+            evaluateAtEndPoint();
         }
 
         /**
@@ -772,7 +770,7 @@ public class GeodeticCalculator {
          * This method stores the projected coordinates in the {@link #point} array and returns
          * the derivative ∂y/∂x.
          */
-        final double evaluateAtEndPoint() throws TransformException {
+        final void evaluateAtEndPoint() throws TransformException {
             if ((λ2 - λ1) * α1 < 0) {
                 λ2 += 2*PI * signum(α1);          // We need λ₁ < λ₂ if heading
east, or λ₁ > λ₂ if heading west.
             }
@@ -789,11 +787,10 @@ 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 tx = m00*dx + m01*dy;
-            final double ty = m10*dx + m11*dy;
-            return ty / tx;
+            final double αx = cos(α2);              // sin(π/2 - α) = -sin(α - π/2)
= cos(α)
+            final double αy = sin(α2);              // cos(π/2 - α) = +cos(α - π/2)
= sin(α)
+            dx = m00*αx + m01*αy;
+            dy = m10*αx + m11*αy;
         }
 
         /**
@@ -857,23 +854,24 @@ public class GeodeticCalculator {
         /**
          * 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.
+         * in the {@link #point} array and the derivative (∂y/∂x) in the {@linkplain
#dx dx} and {@linkplain #dy dy} fields.
          *
          * @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 {
+        protected void evaluateAt(final double t) throws TransformException {
             α1 = IEEEremainder((t - 0.5) * (2*PI), 2*PI);
             validity |= STARTING_AZIMUTH;
             computeEndPoint();
-            final double d = evaluateAtEndPoint();
+            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;
+            final double d = dx;
+            dx = dy;
+            dy = -d;        // Perpendicular direction.
         }
 
         /**
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/j2d/ShapeUtilitiesExt.java
b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/j2d/ShapeUtilitiesExt.java
index b33a665..eb578dc 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/j2d/ShapeUtilitiesExt.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/j2d/ShapeUtilitiesExt.java
@@ -75,7 +75,7 @@ public final class ShapeUtilitiesExt {
                                final double εx, final double εy)
     {
         final Bezier bezier = new Bezier(2) {
-            @Override protected double evaluateAt(final double t) {
+            @Override protected void evaluateAt(final double t) {
                 final double x, y, α;
                 if (t == 0) {
                     x = x1;
@@ -96,7 +96,8 @@ public final class ShapeUtilitiesExt {
                 }
                 point[0] = x;
                 point[1] = y;
-                return α;
+                dx = 2.5;           // Artificial factor for testing purpose. Otherwise could
be 1.
+                dy = dx * α;
             }
         };
         bezier.εx = εx;


Mime
View raw message