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: Approximate geodesic path by Bézier curve: https://issues.apache.org/jira/browse/SIS-454 Current implémentation builds a single linear, quadratic or cubic Bézier curve. It should be improved in a future version by building a chain of Bézier curves.
Date Tue, 14 May 2019 16:04:58 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 c16edc2  Approximate geodesic path by Bézier curve: https://issues.apache.org/jira/browse/SIS-454
Current implémentation builds a single linear, quadratic or cubic Bézier curve. It should
be improved in a future version by building a chain of Bézier curves.
c16edc2 is described below

commit c16edc2e061543f490abe819bc1a4e84a91f1cba
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Tue May 14 17:59:19 2019 +0200

    Approximate geodesic path by Bézier curve: https://issues.apache.org/jira/browse/SIS-454
    Current implémentation builds a single linear, quadratic or cubic Bézier curve.
    It should be improved in a future version by building a chain of Bézier curves.
---
 .../internal/referencing/PositionTransformer.java  |  35 +++-
 .../internal/referencing/j2d/ShapeUtilities.java   |  37 ++++
 .../apache/sis/referencing/GeodeticCalculator.java | 186 ++++++++++++++-------
 .../sis/referencing/GeodeticCalculatorTest.java    |  29 ++++
 4 files changed, 224 insertions(+), 63 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
index 6252f80..707c8bd 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
@@ -23,6 +23,7 @@ import org.opengis.referencing.crs.CoordinateReferenceSystem;
 import org.opengis.referencing.operation.CoordinateOperation;
 import org.opengis.referencing.operation.CoordinateOperationFactory;
 import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
 import org.apache.sis.internal.system.DefaultFactories;
@@ -109,6 +110,7 @@ public final class PositionTransformer extends GeneralDirectPosition {
      * Those fields are left to {@code null} value if the transform is identity.
      *
      * @see #setSourceCRS(CoordinateReferenceSystem)
+     * @see #inverse()
      */
     private transient MathTransform forward, inverse;
 
@@ -211,6 +213,19 @@ public final class PositionTransformer extends GeneralDirectPosition
{
     }
 
     /**
+     * Returns the inverse transform, computed when first needed.
+     */
+    private MathTransform inverse() throws TransformException {
+        if (inverse == null) {
+            if (!Utilities.equalsIgnoreMetadata(lastCRS, defaultCRS)) {
+                setSourceCRS(defaultCRS);
+            }
+            inverse = (forward != null) ? forward.inverse() : MathTransforms.identity(getDimension());
+        }
+        return inverse;
+    }
+
+    /**
      * Returns a new point with the same coordinates than this one, but transformed to the
default CRS.
      * This method never returns {@code this}, so the returned point does not need to be
cloned.
      *
@@ -218,12 +233,18 @@ public final class PositionTransformer extends GeneralDirectPosition
{
      * @throws TransformException if a coordinate transformation was required and failed.
      */
     public DirectPosition inverseTransform() throws TransformException {
-        if (!Utilities.equalsIgnoreMetadata(lastCRS, defaultCRS)) {
-            setSourceCRS(defaultCRS);
-        }
-        if (inverse == null) {
-            inverse = (forward != null) ? forward.inverse() : MathTransforms.identity(getDimension());
-        }
-        return inverse.transform(this, new GeneralDirectPosition(getCoordinateReferenceSystem()));
+        return inverse().transform(this, new GeneralDirectPosition(getCoordinateReferenceSystem()));
+    }
+
+    /**
+     * Transforms this point to the default CRS and stores the result in the given array,
and returns the derivative.
+     * The {@code target} array length should be {@code ReferencingUtilities.getDimension(defaultCRS)}.
+     *
+     * @param  target  where to store the transformed coordinates.
+     * @return the derivative (Jacobian matrix) at the location of this point.
+     * @throws TransformException if a coordinate transformation was required and failed.
+     */
+    public Matrix inverseTransform(final double[] target) throws TransformException {
+        return MathTransforms.derivativeAndTransform(inverse(), ordinates, 0, target, 0);
     }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/ShapeUtilities.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/ShapeUtilities.java
index e0785fe..3a240c7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/ShapeUtilities.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/ShapeUtilities.java
@@ -491,6 +491,43 @@ public final class ShapeUtilities extends Static {
     }
 
     /**
+     * Returns a point on the given linear, quadratic or cubic Bézier curve.
+     *
+     * @param  bezier  a {@link Line2D}, {@link QuadCurve2D} or {@link CubicCurve2D}.
+     * @param  t       a parameter from 0 to 1 inclusive.
+     * @return a point on the curve for the given <var>t</var> parameter.
+     *
+     * @see <a href="https://en.wikipedia.org/wiki/B%C3%A9zier_curve">Bézier curve
on Wikipedia</a>
+     */
+    public static Point2D.Double pointOnBezier(final Shape bezier, final double t) {
+        final double x, y;
+        final double mt = 1 - t;
+        if (bezier instanceof Line2D) {
+            final Line2D z = (Line2D) bezier;
+            x = mt * z.getX1()  +  t * z.getX2();
+            y = mt * z.getY1()  +  t * z.getY2();
+        } else if (bezier instanceof QuadCurve2D) {
+            final QuadCurve2D z = (QuadCurve2D) bezier;
+            final double a = mt * mt;
+            final double b = mt * t * 2;
+            final double c =  t * t;
+            x = a * z.getX1()  +  b * z.getCtrlX()  +  c * z.getX2();
+            y = a * z.getY1()  +  b * z.getCtrlY()  +  c * z.getY2();
+        } else if (bezier instanceof CubicCurve2D) {
+            final CubicCurve2D z = (CubicCurve2D) bezier;
+            final double a = mt * mt * mt;
+            final double b = mt * mt * t  * 3;
+            final double c = mt * (t * t) * 3;
+            final double d =  t *  t * t;
+            x = a * z.getX1()  +  b * z.getCtrlX1()  +  c * z.getCtrlX2()  +  d * z.getX2();
+            y = a * z.getY1()  +  b * z.getCtrlY1()  +  c * z.getCtrlY2()  +  d * z.getY2();
+        } else {
+            throw new IllegalArgumentException();
+        }
+        return new Point2D.Double(x, y);
+    }
+
+    /**
      * Returns the center of a circle passing by the 3 given points. The distance between
the returned
      * point and any of the given points will be constant; it is the circle radius.
      *
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 b671b0d..bd24a79 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
@@ -16,6 +16,7 @@
  */
 package org.apache.sis.referencing;
 
+import java.awt.Shape;
 import java.util.Locale;
 import java.io.IOException;
 import java.io.UncheckedIOException;
@@ -23,6 +24,7 @@ import javax.measure.Unit;
 import javax.measure.quantity.Length;
 
 import org.opengis.referencing.datum.Ellipsoid;
+import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.TransformException;
 import org.opengis.referencing.crs.GeographicCRS;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
@@ -35,6 +37,7 @@ import org.apache.sis.measure.Latitude;
 import org.apache.sis.geometry.CoordinateFormat;
 import org.apache.sis.internal.referencing.PositionTransformer;
 import org.apache.sis.internal.referencing.ReferencingUtilities;
+import org.apache.sis.internal.referencing.j2d.ShapeUtilities;
 import org.apache.sis.internal.referencing.Resources;
 import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.util.resources.Vocabulary;
@@ -91,7 +94,7 @@ public class GeodeticCalculator {
      * The ellipsoid on which geodetic computations are performed.
      * This ellipsoid is inferred from the coordinate reference system specified at construction
time.
      */
-    protected final Ellipsoid ellipsoid;
+    final Ellipsoid ellipsoid;
 
     /**
      * The radius of a hypothetical sphere having the same surface than the {@linkplain #ellipsoid}.
@@ -156,9 +159,12 @@ public class GeodeticCalculator {
      * Users should invoke {@link #create(CoordinateReferenceSystem)} instead, which will
choose a subtype
      * based on the given coordinate reference system.
      *
+     * <p>This class is currently not designed for sub-classing outside this package.
If in a future version we want to
+     * relax this restriction, we should revisit the package-private API in order to commit
to a safer protected API.</p>
+     *
      * @param  crs  the reference system for the {@link Position} arguments and return values.
      */
-    protected GeodeticCalculator(final CoordinateReferenceSystem crs) {
+    GeodeticCalculator(final CoordinateReferenceSystem crs) {
         ArgumentChecks.ensureNonNull("crs", crs);
         final GeographicCRS geographic = ReferencingUtilities.toNormalizedGeographicCRS(crs,
true, true);
         if (geographic == null) {
@@ -170,18 +176,6 @@ public class GeodeticCalculator {
     }
 
     /**
-     * Creates a new geodetic calculator associated with default ellipsoid.
-     * The ellipsoid and the geodetic datum of {@linkplain #getPositionCRS() position CRS}
-     * are the ones associated to {@link CommonCRS#defaultGeographic()}.
-     * The axes are (<var>latitude</var>, <var>longitude</var>) in
degrees.
-     *
-     * @return a new geodetic calculator using the default geodetic datum.
-     */
-    public static GeodeticCalculator create() {
-        return new GeodeticCalculator(CommonCRS.DEFAULT.geographic());
-    }
-
-    /**
      * Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
      * All {@code GeodeticCalculator} methods having a {@link Position} argument
      * or return value will use that specified CRS.
@@ -376,7 +370,7 @@ public class GeodeticCalculator {
      */
     public double getStartingAzimuth() {
         if ((validity & STARTING_AZIMUTH) == 0) {
-            computeTrack();
+            computeDistance();
         }
         return toDegrees(α1);
     }
@@ -395,7 +389,7 @@ public class GeodeticCalculator {
             if ((validity & END_POINT) == 0) {
                 computeEndPoint();                      // Compute also ending azimuth from
start point and distance.
             } else {
-                computeTrack();                         // Compute also ending azimuth from
start point and end point.
+                computeDistance();                         // Compute also ending azimuth
from start point and end point.
             }
         }
         return toDegrees(α2);
@@ -429,7 +423,7 @@ public class GeodeticCalculator {
      */
     public double getGeodesicDistance() {
         if ((validity & GEODESIC_DISTANCE) == 0) {
-            computeTrack();
+            computeDistance();
         }
         return geodesicDistance;
     }
@@ -447,44 +441,6 @@ public class GeodeticCalculator {
     }
 
     /**
-     * Computes the end point from the start point, the azimuth and the geodesic distance.
-     * 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
-     * spherical formulas. Subclasses should override if they can provide ellipsoidal formulas.</p>
-     *
-     * @throws IllegalStateException if the azimuth and the distance have not been set.
-     */
-    void computeEndPoint() {
-        if ((validity & (START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE))
-                     != (START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE))
-        {
-            throw new IllegalStateException((validity & START_POINT) == 0
-                    ? Resources.format(Resources.Keys.StartOrEndPointNotSet_1, 0)
-                    : Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
-        }
-        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Δσ_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Δσ);
-        validity |= END_POINT;
-    }
-
-    /**
      * Computes the geodetic distance and azimuths from the start point and end point.
      * This method should be invoked if the distance or an azimuth is requested while
      * {@link #STARTING_AZIMUTH}, {@link #ENDING_AZIMUTH} or {@link #GEODESIC_DISTANCE}
@@ -498,7 +454,7 @@ public class GeodeticCalculator {
      *
      * @throws IllegalStateException if the distance or azimuth has not been set.
      */
-    private void computeTrack() {
+    private void computeDistance() {
         if ((validity & (START_POINT | END_POINT))
                      != (START_POINT | END_POINT))
         {
@@ -531,6 +487,124 @@ public class GeodeticCalculator {
     }
 
     /**
+     * Computes the end point from the start point, the azimuth and the geodesic distance.
+     * 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
+     * spherical formulas. Subclasses should override if they can provide ellipsoidal formulas.</p>
+     *
+     * @throws IllegalStateException if the azimuth and the distance have not been set.
+     */
+    void computeEndPoint() {
+        if ((validity & (START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE))
+                     != (START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE))
+        {
+            throw new IllegalStateException((validity & START_POINT) == 0
+                    ? Resources.format(Resources.Keys.StartOrEndPointNotSet_1, 0)
+                    : Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
+        }
+        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Δσ_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Δσ);
+        validity |= END_POINT;
+    }
+
+    /**
+     * Creates an approximation of the geodesic track from start point to end point as a
Java2D object.
+     * The coordinates are expressed in the coordinate reference system specified at creation
time.
+     * The approximation uses linear, quadratic or cubic Bézier curves.
+     * The returned path has the following characteristics:
+     *
+     * <ol>
+     *   <li>The first point is {@link #getStartPoint()}.</li>
+     *   <li>The beginning of the curve (more specifically, the tangent at starting
point) is oriented toward the direction given
+     *       by {@linkplain #getStartingAzimuth()}, adjusted for the map projection (if any)
deformation at that location.</li>
+     *   <li>The curve passes at least by the midway point.</li>
+     *   <li>The end of the curve (more specifically, the tangent at ending point)
is oriented toward the direction given by
+     *       {@linkplain #getEndingAzimuth()}, adjusted for the map projection (if any) deformation
at that location.</li>
+     *   <li>The last point is {@link #getEndPoint()}.</li>
+     * </ol>
+     *
+     * <b>Limitations:</b>
+     * current implementation builds a single linear, quadratic or cubic Bézier curve. It
does not yet create a chain
+     * of Bézier curves. Consequently the errors may be larger than the given {@code tolerance}
threshold. Another
+     * limitation is that this method depends on the presence of {@code java.desktop} module.
Those limitations may be
+     * addressed in a future version (see <a href="https://issues.apache.org/jira/browse/SIS-453">SIS-453</a>).
+     *
+     * @param  tolerance  maximal error between the approximated curve and actual geodesic
track
+     *                    in the units of measurement given by {@link #getDistanceUnit()}.
+     *                    <em>See limitations above</em>.
+     * @return an approximation of geodesic track as Bézier curves in a Java2D object.
+     * @throws TransformException if the coordinates can not be transformed to {@linkplain
#getPositionCRS() position CRS}.
+     * @throws IllegalStateException if some required properties have not been specified.
+     */
+    public Shape toGeodesicPath2D(double tolerance) throws TransformException {
+        if ((validity & (START_POINT | STARTING_AZIMUTH | END_POINT | ENDING_AZIMUTH
| GEODESIC_DISTANCE))
+                     != (START_POINT | STARTING_AZIMUTH | END_POINT | ENDING_AZIMUTH | GEODESIC_DISTANCE))
+        {
+            if ((validity & END_POINT) == 0) {
+                computeEndPoint();
+            } else {
+                computeDistance();
+            }
+        }
+        tolerance *= (180/PI) / radius;                                     // Angular tolerance
in degrees.
+        final double d1, x1, y1, d2, x2, y2;                                // Parameters
for the Bezier curve.
+        final double[] transformed = new double[ReferencingUtilities.getDimension(userToGeodetic.defaultCRS)];
+        d1 = slope(α1, geographic(φ1, λ1).inverseTransform(transformed)); x1 = transformed[0];
y1 = transformed[1];
+        d2 = slope(α2, geographic(φ2, λ2).inverseTransform(transformed)); x2 = transformed[0];
y2 = transformed[1];
+        final double sφ2 = φ2;                                              // Save setting
before modification.
+        final double sλ2 = λ2;
+        final double sα2 = α2;
+        final double sd  = geodesicDistance;
+        try {
+            geodesicDistance /= 2;
+            computeEndPoint();
+            final Matrix d = geographic(φ2, λ2).inverseTransform(transformed);      //
Coordinates of midway point.
+            double εx;                                                              // Tolerance
for φ (first coordinate).
+            double εy = tolerance / cos(φ2);                                        //
Tolerance for λ (second coordinate).
+            εx = d.getElement(0,0)*tolerance + d.getElement(0,1)*εy;                //
Tolerance for x in user CRS.
+            εy = d.getElement(1,0)*tolerance + d.getElement(1,1)*εy;                //
Tolerance for y in user CRS.
+            return ShapeUtilities.bezier(x1, y1, transformed[0], transformed[1], x2, y2,
d1, d2, εx, εy);
+        } finally {
+            φ2 = sφ2;                                                               //
Restore the setting previously saved.
+            λ2 = sλ2;
+            α2 = sα2;
+            geodesicDistance = sd;
+        }
+    }
+
+    /**
+     * Returns the tangent of the given angle converted to the user CRS space.
+     *
+     * @param  α  azimuth angle in radians, with 0 pointing toward north and values increasing
clockwise.
+     * @param  d  Jacobian matrix from (φ,λ) to the user coordinate reference system.
+     * @return ∂y/∂x.
+     */
+    private double slope(final double α, final Matrix d) {
+        final double dx = cos(α);     // sin(π/2 - α) = -sin(α - π/2) = cos(α)
+        final double dy = sin(α);     // cos(π/2 - α) = +cos(α - π/2) = sin(α)
+        final double tx = d.getElement(0,0)*dx + d.getElement(0,1)*dy;
+        final double ty = d.getElement(1,0)*dx + d.getElement(1,1)*dy;
+        return ty / tx;
+    }
+
+    /**
      * Returns a string representation of start point, end point, azimuths and distance.
      *
      * @return a string representation of this calculator state.
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 512a00b..25f0d54 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
@@ -16,9 +16,12 @@
  */
 package org.apache.sis.referencing;
 
+import java.awt.Shape;
+import java.awt.geom.Point2D;
 import java.util.Random;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.referencing.j2d.ShapeUtilities;
 import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.geometry.DirectPosition2D;
 import org.apache.sis.measure.Units;
@@ -169,4 +172,30 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
         c.setGeodesicDistance(18743000);
         assertPositionEquals(121.8, 31.4, c.getEndPoint(), 0.01);
     }
+
+    /**
+     * Tests {@link GeodeticCalculator#toGeodesicPath2D(double)}. This method uses a CRS
+     * that swap axis order as a way to verify that user-specified CRS is taken in account.
+     *
+     * @throws TransformException if an error occurred while transforming coordinates.
+     */
+    @Test
+    @DependsOnMethod("testUsingTransform")
+    public void testToGeodesicPath2D() throws TransformException {
+        final GeodeticCalculator c = new GeodeticCalculator(CommonCRS.SPHERE.normalizedGeographic());
+        c.setStartPoint(-33.0, -71.6);                  // Valparaíso
+        c.setEndPoint  ( 31.4, 121.8);                  // Shanghai
+        final Shape path = c.toGeodesicPath2D(1000);
+        assertPointEquals( -71.6, -33.0, ShapeUtilities.pointOnBezier(path, 0));
+        assertPointEquals( 121.8,  31.4, ShapeUtilities.pointOnBezier(path, 1));
+        assertPointEquals(-159.2,  -6.8, ShapeUtilities.pointOnBezier(path, 0.5));
+    }
+
+    /**
+     * Asserts that a Java2D point is equal to the expected value.
+     */
+    private static void assertPointEquals(final double x, final double y, final Point2D actual)
{
+        assertEquals("x", x, actual.getX(), 0.05);
+        assertEquals("y", y, actual.getY(), 0.05);
+    }
 }


Mime
View raw message