sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/02: Replace α₁ and α₂ angles in GeodeticCalculator by dφ₁/dλ₁ and dφ₂/dλ₂ vector components. This is the complement of previous commit that applied similar change in Bezier class. It resolves issues with infinities, which fix GeodeticCalculator.createCircularRegion2D(…).
Date Fri, 24 May 2019 08:07:30 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 1ead312c9570ce83e3fe7c79c1491337c43069ec
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Fri May 24 09:46:39 2019 +0200

    Replace α₁ and α₂ angles in GeodeticCalculator by dφ₁/dλ₁ and dφ₂/dλ₂
vector components.
    This is the complement of previous commit that applied similar change in Bezier class.
    It resolves issues with infinities, which fix GeodeticCalculator.createCircularRegion2D(…).
---
 .../org/apache/sis/distance/LatLonPointRadius.java |   4 +-
 .../apache/sis/referencing/GeodeticCalculator.java | 117 +++++++++++----------
 .../sis/referencing/GeodeticCalculatorTest.java    |   7 +-
 3 files changed, 70 insertions(+), 58 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/distance/LatLonPointRadius.java
b/core/sis-referencing/src/main/java/org/apache/sis/distance/LatLonPointRadius.java
index 0dc59dc..267346f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/distance/LatLonPointRadius.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/distance/LatLonPointRadius.java
@@ -36,9 +36,9 @@ import org.apache.sis.referencing.GeodeticCalculator;
  * Represents a 2D point associated with a radius to enable great circle
  * estimation on earth surface.
  *
- * <div class="warning"><b>Warning:</b> This class may be refactored as
a geometric object in a future SIS version.
- * Current implementation does not verify the CRS of circle center or the datum.</div>
+ * @deprecated Replaced by {@link org.apache.sis.referencing.GeodeticCalculator#createCircularRegion2D(double)}.
  */
+@Deprecated
 public class LatLonPointRadius {
 
   private final DirectPosition center;
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 eed90ff..5123b3b 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
@@ -127,13 +127,22 @@ public class GeodeticCalculator {
     private double φ2, λ2;
 
     /**
-     * The azimuth at start point and end point, in radians between -π and +π.
-     * 0° point toward North and values are increasing clockwise.
+     * The azimuth at start point and end point as vector components.
+     * Angles can be obtained as below (reminder: tan(π/2 − α) = 1/tan(α)):
+     *
+     * <ul>
+     *   <li>Arithmetic angle: {@code atan2(dφ, dλ)} (radians increasing anticlockwise)</li>
+     *   <li>Geographic angle: {@code atan2(dλ, dφ)} gives the azimuth in radians
between -π and +π
+     *       with 0° pointing toward North and values increasing clockwise.</li>
+     * </ul>
+     *
+     * Those vectors may not be normalized to unitary vectors. This is often not needed because
most formulas are
+     * written in a way that cancel the magnitude. If nevertheless needed, normalization
will be applied in formulas.
      *
      * @see #STARTING_AZIMUTH
      * @see #ENDING_AZIMUTH
      */
-    private double α1, α2;
+    private double dφ1, dλ1, dφ2, dλ2;
 
     /**
      * The shortest distance from the starting point ({@link #φ1},{@link #λ1}) to the end
point ({@link #φ2},{@link #λ2}).
@@ -155,9 +164,9 @@ public class GeodeticCalculator {
 
     /**
      * A bitmask specifying which information are valid. For example if the {@link #END_POINT}
bit is not set,
-     * then {@link #φ2} and {@link #λ2} need to be computed, which implies the computation
of {@link #α1} and
-     * {@link #α2} as well. If the {@link #GEODESIC_DISTANCE} bit is not set, then {@link
#geodesicDistance}
-     * needs to be computed, which implies recomputation of {@link #α1} and {@link #α2}
as well.
+     * then {@link #φ2} and {@link #λ2} need to be computed, which implies the computation
of ∂φ/∂λ as well.
+     * If the {@link #GEODESIC_DISTANCE} bit is not set, then {@link #geodesicDistance} needs
to be computed,
+     * which implies recomputation of ∂φ/∂λ as well.
      */
     private int validity;
 
@@ -196,6 +205,10 @@ public class GeodeticCalculator {
      * or return value will use that specified CRS.
      * That CRS is the value returned by {@link #getPositionCRS()}.
      *
+     * <p><b>Limitation:</b>
+     * current implementation uses only spherical formulas.
+     * Implementation using ellipsoidal formulas will be provided in a future Apache SIS
release.</p>
+     *
      * @param  crs  the reference system for the {@link Position} objects.
      * @return a new geodetic calculator using the specified CRS.
      */
@@ -377,7 +390,7 @@ public class GeodeticCalculator {
         if (isInvalid(STARTING_AZIMUTH)) {
             computeDistance();
         }
-        return toDegrees(α1);
+        return toDegrees(atan2(dλ1, dφ1));                  // tan(π/2 − θ)  =  1/tan(θ)
     }
 
     /**
@@ -391,9 +404,11 @@ public class GeodeticCalculator {
      *
      * @see #setGeodesicDistance(double)
      */
-    public void setStartingAzimuth(final double azimuth) {
+    public void setStartingAzimuth(double azimuth) {
         ArgumentChecks.ensureFinite("azimuth", azimuth);
-        α1 = toRadians(IEEEremainder(azimuth, 360));
+        azimuth = toRadians(azimuth);
+        dφ1 = cos(azimuth);                                 // sin(π/2 − θ)  =  cos(θ)
+        dλ1 = sin(azimuth);                                 // cos(π/2 − θ)  =  sin(θ)
         validity |= STARTING_AZIMUTH;
         validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH);
     }
@@ -415,7 +430,7 @@ public class GeodeticCalculator {
                 computeDistance();                         // Compute also ending azimuth
from start point and end point.
             }
         }
-        return toDegrees(α2);
+        return toDegrees(atan2(dλ2, dφ2));
     }
 
     /**
@@ -542,17 +557,16 @@ public class GeodeticCalculator {
 
         final double cosφ1_sinφ2 = cosφ1 * sinφ2;
         final double cosφ2_sinφ1 = cosφ2 * sinφ1;
-        final double α1y = cosφ2 * sinΔλ;
-        final double α1x = cosφ1_sinφ2 - cosφ2_sinφ1*cosΔλ;
-
-        α1 = atan2(α1y, α1x);
-        α2 = atan2(cosφ1*sinΔλ, cosφ1_sinφ2*cosΔλ - cosφ2_sinφ1);
+        dλ1 = cosφ2*sinΔλ;
+        dλ2 = cosφ1*sinΔλ;
+        dφ1 = cosφ1_sinφ2 - cosφ2_sinφ1*cosΔλ;
+        dφ2 = cosφ1_sinφ2*cosΔλ - cosφ2_sinφ1;
         /*
          * Δσ = acos(sinφ₁⋅sinφ₂ + cosφ₁⋅cosφ₂⋅cosΔλ) is a first estimation
inaccurate for small distances.
          * Δσ = atan2(…) computes the same value but with better accuracy.
          */
         double Δσ = sinφ1*sinφ2 + cosφ1*cosφ2*cosΔλ;        // Actually Δσ = acos(…).
-        Δσ = atan2(hypot(α1x, α1y), Δσ);
+        Δσ = atan2(hypot(dλ1, dφ1), Δσ);
         geodesicDistance = radius * Δσ;
         validity |= (STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
     }
@@ -562,7 +576,7 @@ public class GeodeticCalculator {
      * This method should be invoked if the end point or ending azimuth is requested while
      * {@link #END_POINT} validity flag is not set.
      *
-     * <p>The default implementation computes {@link #φ2}, {@link #λ2} and {@link
#α2} using
+     * <p>The default implementation computes {@link #φ2}, {@link #λ2} and ∂φ/∂λ
derivatives using
      * spherical formulas. Subclasses should override if they can provide ellipsoidal formulas.</p>
      *
      * @throws IllegalStateException if the azimuth and the distance have not been set.
@@ -573,23 +587,24 @@ public class GeodeticCalculator {
                     ? Resources.format(Resources.Keys.StartOrEndPointNotSet_1, 0)
                     : Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
         }
+        final double vm    = hypot(dφ1, dλ1);
         final double Δσ    = geodesicDistance / radius;
         final double sinΔσ = sin(Δσ);
         final double cosΔσ = cos(Δσ);
         final double sinφ1 = sin(φ1);
         final double cosφ1 = cos(φ1);
-        final double sinα1 = sin(α1);
-        final double cosα1 = cos(α1);
-
+        final double sinα1 = dλ1 / vm;          // α₁ is the azimuth at starting point.
+        final double cosα1 = dφ1 / vm;
         final double sinΔσ_cosα1 = sinΔσ * cosα1;
         final double Δλy = sinΔσ * sinα1;
         final double Δλx = cosΔσ*cosφ1 - sinφ1*sinΔσ_cosα1;
         final double Δλ  = atan2(Δλy, Δλx);
 
         final double sinφ2 = sinφ1*cosΔσ + cosφ1*sinΔσ_cosα1;       // First estimation
of φ2.
-        φ2 = atan(sinφ2 / hypot(Δλx, Δλy));                         // Improve accuracy
close to poles.
-        λ2 = IEEEremainder(λ1 + Δλ, 2*PI);
-        α2 = atan2(sinα1, cosΔσ*cosα1 - sinφ1/cosφ1 * sinΔσ);
+        φ2  = atan(sinφ2 / hypot(Δλx, Δλy));                        // Improve accuracy
close to poles.
+        λ2  = IEEEremainder(λ1 + Δλ, 2*PI);
+        dφ2 = cosΔσ*cosα1 - sinφ1/cosφ1 * sinΔσ;
+        dλ2 = sinα1;
         validity |= END_POINT;
     }
 
@@ -605,9 +620,8 @@ public class GeodeticCalculator {
         if (isInvalid(END_POINT)) {
             computeEndPoint();
         }
-        φ1 = φ2;
-        λ1 = λ2;
-        α1 = α2;
+        φ1  = φ2;  dφ1 = dφ2;
+        λ1  = λ2;  dλ1 = dλ2;
         validity = START_POINT | STARTING_AZIMUTH;
     }
 
@@ -708,7 +722,7 @@ public class GeodeticCalculator {
          * 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 αi, αf, φf, λf, distance, length;
+        private final double dφi, dλi, dφf, dλf, φf, λf, distance, length;
 
         /**
          * {@link #validity} flags at the time {@code PathBuilder} is instantiated.
@@ -726,11 +740,8 @@ public class GeodeticCalculator {
          */
         PathBuilder(final double εx) {
             super(ReferencingUtilities.getDimension(userToGeodetic.defaultCRS));
-            αi        = α1;
-            αf        = α2;
-            φf        = φ2;
-            λf        = λ2;
-            α1        = IEEEremainder(α1, 2*PI);            // For reliable signum(α1)
result.
+            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;
             length    = rhumblineLength;
@@ -748,13 +759,11 @@ public class GeodeticCalculator {
         @Override
         protected void evaluateAt(final double t) throws TransformException {
             if (t == 0) {
-                φ2 = φ1;                        // Start point requested.
-                λ2 = λ1;
-                α2 = α1;
+                φ2 = φ1;  dφ2 = dφ1;            // Start point requested.
+                λ2 = λ1;  dλ2 = dλ1;
             } else if (t == 1) {
-                φ2 = φf;                        // End point requested.
-                λ2 = λf;
-                α2 = αf;
+                φ2 = φf;  dφ2 = dφf;            // End point requested.
+                λ2 = λf;  dλ2 = dλf;
             } else {
                 geodesicDistance = distance * t;
                 validity |= GEODESIC_DISTANCE;
@@ -764,13 +773,13 @@ public class GeodeticCalculator {
         }
 
         /**
-         * 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.
+         * Implementation of {@link #evaluateAt(double)} using the current φ₂, λ₂ and
∂φ₂/∂λ₂ values.
+         * This method stores the projected coordinates in the {@link #point} array and stores
the
+         * derivative ∂y/∂x in the {@link #dx}, {@link #dy} fields.
          */
         final void evaluateAtEndPoint() throws TransformException {
-            if ((λ2 - λ1) * α1 < 0) {
-                λ2 += 2*PI * signum(α1);          // We need λ₁ < λ₂ if heading
east, or λ₁ > λ₂ if heading west.
+            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 double m00 = d.getElement(0,0);
@@ -782,13 +791,10 @@ public class GeodeticCalculator {
             εy = m10*tolerance + m11*εy;                                    // Tolerance
for y in user CRS.
             /*
              * Returns the tangent of the ending azimuth converted to the user CRS space.
-             * α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.
+             * d is the Jacobian matrix from (φ,λ) to the user coordinate reference system.
              */
-            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;
+            dx = m00*dφ2 + m01*dλ2;                                         // Reminder:
coordinates in (φ,λ) order.
+            dy = m10*dφ2 + m11*dλ2;
         }
 
         /**
@@ -811,7 +817,8 @@ public class GeodeticCalculator {
              * and check the distance between those two points.
              */
             computeDistance();
-            α1 = αi;
+            dφ1 = dφi;
+            dλ1 = dλi;
             computeEndPoint();
             final DirectPosition p = geographic(φ2, λ2).inverseTransform();
             return abs(p.getOrdinate(0) - x) <= εx &&
@@ -822,10 +829,8 @@ public class GeodeticCalculator {
          * Restores the enclosing {@link GeodeticCalculator} to the state that it has at
{@code PathBuilder} instantiation time.
          */
         final void reset() {
-            α1 = αi;
-            α2 = αf;
-            φ2 = φf;
-            λ2 = λf;
+            dφ1 = dφi;  dφ2 = dφf;  φ2 = φf;
+            dλ1 = dλi;  dλ2 = dλf;  λ2 = λf;
             geodesicDistance = distance;
             rhumblineLength  = length;
             validity         = flags;
@@ -853,7 +858,9 @@ public class GeodeticCalculator {
          */
         @Override
         protected void evaluateAt(final double t) throws TransformException {
-            α1 = IEEEremainder((t - 0.5) * (2*PI), 2*PI);
+            double α1 = IEEEremainder((t - 0.5) * (2*PI), 2*PI);
+            dλ1 = sin(α1);
+            dφ1 = cos(α1);
             validity |= STARTING_AZIMUTH;
             computeEndPoint();
             evaluateAtEndPoint();
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 c904168..b628331 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
@@ -19,6 +19,7 @@ package org.apache.sis.referencing;
 import java.awt.Shape;
 import java.awt.geom.Path2D;
 import java.awt.geom.Point2D;
+import java.awt.geom.Rectangle2D;
 import java.awt.geom.PathIterator;
 import java.util.Arrays;
 import java.util.Random;
@@ -229,7 +230,11 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
         if (VisualCheck.SHOW_WIDGET) {
             VisualCheck.show(region);
         }
-        // TODO: test bounding box.
+        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);
     }
 
     /**


Mime
View raw message