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: Azimuths computed by GeodeticCalculator can be used for angles in a shape only if the projection is conformal. We use the ∂y/∂φ factor of Mercator projection for applying a correction. Also force the circular shape to use cubic Bézier (disallowing quadratic Bézier in that particular case), which improve quality a lot.
Date Mon, 03 Jun 2019 13:22:29 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 f9009a1  Azimuths computed by GeodeticCalculator can be used for angles in a shape
only if the projection is conformal. We use the ∂y/∂φ factor of Mercator projection for
applying a correction. Also force the circular shape to use cubic Bézier (disallowing quadratic
Bézier in that particular case), which improve quality a lot.
f9009a1 is described below

commit f9009a105e33ce37a9847936f7096c373a2afedd
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Jun 3 15:17:13 2019 +0200

    Azimuths computed by GeodeticCalculator can be used for angles in a shape only if the
projection is conformal.
    We use the ∂y/∂φ factor of Mercator projection for applying a correction.
    Also force the circular shape to use cubic Bézier (disallowing quadratic Bézier in that
particular case), which improve quality a lot.
---
 .../sis/internal/referencing/j2d/Bezier.java       | 21 ++++++----
 .../apache/sis/referencing/GeodeticCalculator.java | 46 ++++++++++++++++------
 .../operation/projection/ConformalProjection.java  |  2 +-
 .../sis/referencing/GeodeticCalculatorTest.java    |  4 +-
 4 files changed, 52 insertions(+), 21 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 dcae453..5ecfc5b 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
@@ -107,6 +107,13 @@ public abstract class Bezier {
     protected int depth;
 
     /**
+     * Whether to force the creation of {@link QuadCurve2D}. The default value is {@code
false}, which allows simplification
+     * to {@link CubicCurve2D} or {@link Line2D}. This flag can be set to {@code true} when
building circular shapes, in which
+     * case simplifications to {@link CubicCurve2D} produce bad results.
+     */
+    protected boolean forceCubic;
+
+    /**
      * Creates a new builder.
      *
      * @param  dimension  length of the {@link #point} array. Must be at least 2.
@@ -201,8 +208,8 @@ public abstract class Bezier {
         if (tx < εx) εx = tx;                                       // Take smallest
tolerance values.
         if (ty < εy) εy = ty;
         if (curve(x1, y1, x2, y2, x3, y3, x4, y4, dx4, dy4)) {
-            x0  = x4;
-            y0  = y4;
+            x0  =  x4;
+            y0  =  y4;
             dx0 = dx4;
             dy0 = dy4;
         } else {
@@ -268,8 +275,8 @@ public abstract class Bezier {
         y1 -= y0;  y2 -= y0;  y3 -= y0;
         final double Δx = x4 - x0;
         final double Δy = y4 - y0;
-        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))
+        if (!forceCubic && (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))
         {
             /*
              * Verify that all points are on the line joining P₀ to P₄. If any coordinate
tested below is NaN,
@@ -417,12 +424,12 @@ public abstract class Bezier {
          *
          *     Δ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
+         * We multiply tolerance factor by 2 because 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) <= 2*εx &&      // P₀ is zero.
-            abs(Δqy = (3*(by - ay) - Δy)/2) <= 2*εy)
+        if (!forceCubic && 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 db51946..8457227 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
@@ -33,7 +33,6 @@ import org.opengis.referencing.crs.CoordinateReferenceSystem;
 import org.opengis.geometry.coordinate.Position;
 import org.opengis.geometry.DirectPosition;
 
-import org.apache.sis.io.TableAppender;
 import org.apache.sis.measure.AngleFormat;
 import org.apache.sis.measure.Latitude;
 import org.apache.sis.measure.Units;
@@ -48,6 +47,7 @@ import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.util.resources.Vocabulary;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.ArgumentChecks;
+import org.apache.sis.io.TableAppender;
 
 import static java.lang.Math.*;
 
@@ -494,6 +494,21 @@ public class GeodeticCalculator {
     }
 
     /**
+     * Computes (∂y/∂φ)⁻¹ where (∂y/∂φ) is the partial derivative of Northing
values in a Mercator projection
+     * at the given latitude on an ellipsoid with semi-major axis length of 1. There is no
method for partial
+     * derivative of Easting values since it is 1 everywhere. This derivative is cos(φ)
on a sphere and close
+     * but slightly different on an ellipsoid.
+     *
+     * @param  φ  the latitude in radians.
+     * @return the northing derivative of a Mercator projection at the given latitude on
an ellipsoid with a=1.
+     *
+     * @see org.apache.sis.referencing.operation.projection.ConformalProjection#dy_dφ
+     */
+    private double dφ_dy(final double φ) {
+        return cos(φ);
+    }
+
+    /**
      * Computes the length of rhumb line from start point to end point.
      *
      * @see <a href="https://en.wikipedia.org/wiki/Rhumb_line">Rhumb line on Wikipedia</a>
@@ -780,20 +795,28 @@ public class GeodeticCalculator {
             if ((λ2 - λ1) * dλ1 < 0) {            // Reminder: Δλ or dλ₁ may be
zero.
                 λ2 += 2*PI * signum(dλ1);         // We need λ₁ < λ₂ if heading
east, or λ₁ > λ₂ if heading west.
             }
-            final Matrix d = geographic(φ2, λ2).inverseTransform(point);    // Coordinates
and Jacobian of point.
+            final Matrix d = geographic(φ2, λ2).inverseTransform(point);    // `point`
coordinates in user-specified CRS.
+            final double dφ_dy = dφ_dy(φ2);                                 // ∂φ/∂y
= cos(φ) for Mercator on a sphere of radius 1.
             final double m00 = d.getElement(0,0);
             final double m01 = d.getElement(0,1);
             final double m10 = d.getElement(1,0);
             final double m11 = d.getElement(1,1);
-            εy = tolerance / cos(abs(φ2));                                  // Tolerance
for λ (second coordinate).
-            εx = m00*tolerance + m01*εy;                                    // Tolerance
for x in user CRS.
-            εy = m10*tolerance + m11*εy;                                    // Tolerance
for y in user CRS.
+            double t = tolerance / dφ_dy;                                   // Tolerance
for λ (second coordinate).
+            εx = m00*tolerance + m01*t;                                     // Tolerance
for x in user CRS.
+            εy = m10*tolerance + m11*t;                                     // Tolerance
for y in user CRS.
             /*
-             * Returns the tangent of the ending azimuth converted to the user CRS space.
-             * d is the Jacobian matrix from (φ,λ) to the user coordinate reference system.
+             * Return the tangent of the ending azimuth converted to the user CRS space.
Note that if we draw the shape on
+             * screen with (longitude, latitude) axes, the angles seen on screen are not
the real angles measured on Earth.
+             * In order to see the "real" angles, we need to draw the shape on a conformal
projection such as Mercator.
+             * Said otherwise, the angle value computed from the (dx,dy) vector is "real"
only in a conformal projection.
+             * Consequently if the output CRS is a Mercator projection, then the angle computed
from the (dx,dy) vector
+             * at the end of this method should be the ending azimuth angle unchanged. We
achieve this equivalence by
+             * multiplying dφ by a factor which will cancel the ∂y/∂φ factor of Mercator
projection at that latitude.
+             * Note that there is no need to modify dλ since ∂x/∂λ = 1 everywhere on
Mercator projection with a=1.
              */
-            dx = m00*dφ2 + m01*dλ2;                                         // Reminder:
coordinates in (φ,λ) order.
-            dy = m10*dφ2 + m11*dλ2;
+            t  = dφ2 * dφ_dy;
+            dx = m00*t + m01*dλ2;                                           // Reminder:
coordinates in (φ,λ) order.
+            dy = m10*t + m11*dλ2;
         }
 
         /**
@@ -818,12 +841,13 @@ public class GeodeticCalculator {
         private final double dφi, dλi;
 
         /**
-         * Creates a builder for the given tolerance at equator in degrees.
+         * Creates a builder for the given tolerance in degrees at equator.
          */
         CircularPath(final double εx) {
             super(εx);
             dφi = dφ1;
             dλi = dλ1;
+            forceCubic = true;
         }
 
         /**
@@ -836,7 +860,7 @@ public class GeodeticCalculator {
          */
         @Override
         protected void evaluateAt(final double t) throws TransformException {
-            final double α1 = t * (2*PI) + PI/4;        // Add 45° for hiding little corners
in multiple of 90°.
+            final double α1 = t * (2*PI);
             dλ1 = sin(α1);
             dφ1 = cos(α1);
             validity |= STARTING_AZIMUTH;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
index 2fa54ff..1e119ed 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
@@ -317,7 +317,7 @@ abstract class ConformalProjection extends NormalizedProjection {
      * other projections, since Mercator can be considered as a special case of Lambert Conic
Conformal for instance.
      *
      * <p>In order to get the derivative of the {@link #expΨ(double, double)} function,
caller can multiply
-     * the returned value by by {@code expΨ}.</p>
+     * the returned value by {@code expΨ}.</p>
      *
      * @param  sinφ  the sine of latitude.
      * @param  cosφ  the cosine of latitude.
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 c846cea..b969d95 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
@@ -232,8 +232,8 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
         final Rectangle2D bounds = region.getBounds2D();
         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.
+        assertEquals("width",  2.1, bounds.getWidth(),   0.1);
+        assertEquals("height", 1.8, bounds.getHeight(),  0.1);
     }
 
     /**


Mime
View raw message