sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From rmarec...@apache.org
Subject svn commit: r1701977 - in /sis/branches/JDK8/core/sis-referencing/src: main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
Date Wed, 09 Sep 2015 11:00:58 GMT
Author: rmarechal
Date: Wed Sep  9 11:00:58 2015
New Revision: 1701977

URL: http://svn.apache.org/r1701977
Log:
First draft of a derivative function, not yet simplified.

Modified:
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
    sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java?rev=1701977&r1=1701976&r2=1701977&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
[UTF-8] Wed Sep  9 11:00:58 2015
@@ -20,6 +20,7 @@ import java.util.EnumMap;
 import org.opengis.parameter.ParameterDescriptor;
 import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.OperationMethod;
+import org.apache.sis.referencing.operation.matrix.Matrix2;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
 import org.apache.sis.internal.referencing.provider.TransverseMercatorSouth;
 import org.apache.sis.internal.util.DoubleDouble;
@@ -51,6 +52,7 @@ import static org.apache.sis.math.MathFu
  * all zones and a false northing of 10000000 metres is used for zones in the southern hemisphere.
  *
  * @author  Martin Desruisseaux (Geomatys)
+ * @author  Rémi Maréchal (Geomatys)
  * @since   0.6
  * @version 0.6
  * @module
@@ -213,7 +215,7 @@ public class TransverseMercator extends
     {
         final double λ  = srcPts[srcOff];
         final double φ  = srcPts[srcOff + 1];
-        final double Q  = asinh(tan(φ)) - atanh(sin(φ)*excentricity)*excentricity;
+        final double Q  = asinh(tan(φ)) - atanh(sin(φ) * excentricity) * excentricity;
         final double β  = atan(sinh(Q));
 
         // TODO: sin(atan(x)) = x / sqrt(1+x²)
@@ -221,14 +223,20 @@ public class TransverseMercator extends
         final double η0 = atanh(cos(β) * sin(λ));
         final double ξ0 = asin(sin(β) * cosh(η0));
 
-        // TODO: use trigonometric identities.
-        // See ConformalProjection for example.
+       /*
+        * Assuming that (λ, φ) ↦ Proj((λ, φ))
+        * where Proj is defined by : Proj((λ, φ)) : (η(λ, φ), ξ(λ, φ)).
+        *
+        * => (λ, φ) ↦ (η(λ, φ), ξ(λ, φ)).
+        */
+        //-- ξ(λ, φ)
         final double ξ = h4 * sin(8*ξ0) * cosh(8*η0)
                        + h3 * sin(6*ξ0) * cosh(6*η0)
                        + h2 * sin(4*ξ0) * cosh(4*η0)
                        + h1 * sin(2*ξ0) * cosh(2*η0)
                        + ξ0;
 
+        //-- η(λ, φ)
         final double η = h4 * cos(8*ξ0) * sinh(8*η0)
                        + h3 * cos(6*ξ0) * sinh(6*η0)
                        + h2 * cos(4*ξ0) * sinh(4*η0)
@@ -239,16 +247,77 @@ public class TransverseMercator extends
             dstPts[dstOff    ] = η;
             dstPts[dstOff + 1] = ξ;
         }
+        
         if (!derivate) {
             return null;
         }
+        
+        //-- dQ_dλ = 0;
+        final double dQ_dφ = sqrt(1 + tan(φ) * tan(φ)) - excentricitySquared * cos(φ)
/ (1 - excentricitySquared * sin(φ) * sin(φ));
+        
+        //-- dβ_dλ = 0;
+        final double dβ_dφ = dQ_dφ * cosh(Q) / (1 + sinh(Q) * sinh(Q));
+
+        // TODO: sin(β) = sin(atan(sinh(Q))) = sinh(Q) / sqrt(1 + sinh²(Q))
+        //       cos(β) = cos(atan(sinh(Q))) = 1       / sqrt(1 + sinh²(Q))
+        final double a     =   cos(β) * sin(λ);
+        final double da_dλ =   cos(β) * cos(λ);
+        final double da_dφ = - sin(λ) * dβ_dφ * sin(β);
+        
+        final double dη0_dλ = da_dλ / (1 - a * a);
+        final double dη0_dφ = da_dφ / (1 - a * a);
+        
+        final double b     = sin(β) * cosh(η0);
+        final double db_dλ = sin(β) * dη0_dλ * sinh(η0);
+        final double db_dφ = dβ_dφ  * cos(β) * cosh(η0) + dη0_dφ * sinh(η0) * sin(β);
+        
+        final double dξ0_dλ = db_dλ / sqrt(1 - b * b);
+        final double dξ0_dφ = db_dφ / sqrt(1 - b * b);
+
+       /*
+        * Assuming that Jac(Proj((λ, φ))) is the Jacobian matrix of Proj((λ, φ)) function.
+        *
+        * So derivative Proj((λ, φ)) is defined by :
+        *                     _                            _
+        *                    | dη(λ, φ) / dλ, dη(λ, φ) / dφ |
+        * Jac              = |                              |
+        *    (Proj(λ, φ))    | dξ(λ, φ) / dλ, dξ(λ, φ) / dφ |
+        *                    |_                            _|
+        */
+
+        //-- dξ(λ, φ) / dλ
+        final double dξ_dλ = h4 * (8 * dξ0_dλ * cos(8*ξ0) * cosh(8*η0) + 8 * dη0_dλ
* sinh(8*η0) * sin(8*ξ0))
+                           + h3 * (6 * dξ0_dλ * cos(6*ξ0) * cosh(6*η0) + 6 * dη0_dλ
* sinh(6*η0) * sin(6*ξ0))
+                           + h2 * (4 * dξ0_dλ * cos(4*ξ0) * cosh(4*η0) + 4 * dη0_dλ
* sinh(4*η0) * sin(4*ξ0))
+                           + h1 * (2 * dξ0_dλ * cos(2*ξ0) * cosh(2*η0) + 2 * dη0_dλ
* sinh(2*η0) * sin(2*ξ0))
+                           + dξ0_dλ;
+        //-- dξ(λ, φ) / dφ
+        final double dξ_dφ = h4 * (8 * dξ0_dφ * cos(8*ξ0) * cosh(8*η0) + 8 * dη0_dφ
* sinh(8*η0) * sin(8*ξ0))
+                           + h3 * (6 * dξ0_dφ * cos(6*ξ0) * cosh(6*η0) + 6 * dη0_dφ
* sinh(6*η0) * sin(6*ξ0))
+                           + h2 * (4 * dξ0_dφ * cos(4*ξ0) * cosh(4*η0) + 4 * dη0_dφ
* sinh(4*η0) * sin(4*ξ0))
+                           + h1 * (2 * dξ0_dφ * cos(2*ξ0) * cosh(2*η0) + 2 * dη0_dφ
* sinh(2*η0) * sin(2*ξ0))
+                           + dξ0_dφ;
+
+        //-- dη(λ, φ) / dλ
+        final double dη_dλ = h4 * (-8 * dξ0_dλ * sin(8 * ξ0) * sinh(8*η0) + 8 * dη0_dλ
* cosh(8*η0) * cos(8*ξ0))
+                           + h3 * (-6 * dξ0_dλ * sin(6 * ξ0) * sinh(6*η0) + 6 * dη0_dλ
* cosh(6*η0) * cos(6*ξ0))
+                           + h2 * (-4 * dξ0_dλ * sin(4 * ξ0) * sinh(4*η0) + 4 * dη0_dλ
* cosh(4*η0) * cos(4*ξ0))
+                           + h1 * (-2 * dξ0_dλ * sin(2 * ξ0) * sinh(2*η0) + 2 * dη0_dλ
* cosh(2*η0) * cos(2*ξ0))
+                           + dη0_dλ;
+
+        //-- dη(λ, φ) / dφ
+        final double dη_dφ = h4 * (-8 * dξ0_dφ * sin(8 * ξ0) * sinh(8*η0) + 8 * dη0_dφ
* cosh(8*η0) * cos(8*ξ0))
+                           + h3 * (-6 * dξ0_dφ * sin(6 * ξ0) * sinh(6*η0) + 6 * dη0_dφ
* cosh(6*η0) * cos(6*ξ0))
+                           + h2 * (-4 * dξ0_dφ * sin(4 * ξ0) * sinh(4*η0) + 4 * dη0_dφ
* cosh(4*η0) * cos(4*ξ0))
+                           + h1 * (-2 * dξ0_dφ * sin(2 * ξ0) * sinh(2*η0) + 2 * dη0_dφ
* cosh(2*η0) * cos(2*ξ0))
+                           + dη0_dφ;
 
-        // TODO: compute projection derivative.
-        return null;
+        return new Matrix2(dη_dλ, dη_dφ,
+                           dξ_dλ, dξ_dφ);
     }
 
     /**
-     * Transforms the specified (η,ξ) coordinates and stores the result in {@code dstPts}
(angles in radians).
+     * Transforms the specified (η, ξ) coordinates and stores the result in {@code dstPts}
(angles in radians).
      *
      * @throws ProjectionException if the point can not be converted.
      */

Modified: sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java?rev=1701977&r1=1701976&r2=1701977&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
[UTF-8] Wed Sep  9 11:00:58 2015
@@ -36,6 +36,7 @@ import static java.lang.StrictMath.toRad
  */
 @DependsOn(NormalizedProjectionTest.class)
 public final strictfp class TransverseMercatorTest extends MapProjectionTestCase {
+    
     /**
      * Creates a new instance of {@link TransverseMercator}.
      *
@@ -61,7 +62,6 @@ public final strictfp class TransverseMe
      * @see org.opengis.test.referencing.ParameterizedTransformTest#testTransverseMercator()
      */
     @Test
-    @org.junit.Ignore("Missing implementation of the projection derivative.")
     public void testTransverseMercator() throws FactoryException, TransformException {
         createGeoApiTest(new org.apache.sis.internal.referencing.provider.TransverseMercator()).testTransverseMercator();
     }
@@ -76,7 +76,6 @@ public final strictfp class TransverseMe
      * @see org.opengis.test.referencing.ParameterizedTransformTest#testTransverseMercatorSouthOrientated()
      */
     @Test
-    @org.junit.Ignore("Missing implementation of the projection derivative.")
     public void testTransverseMercatorSouthOrientated() throws FactoryException, TransformException
{
         createGeoApiTest(new TransverseMercatorSouth()).testTransverseMercatorSouthOrientated();
     }
@@ -87,7 +86,6 @@ public final strictfp class TransverseMe
      * @throws TransformException Should never happen.
      */
     @Test
-    @org.junit.Ignore("Missing implementation of the projection derivative.")
     public void testSphericalDerivative() throws TransformException {
         createNormalizedProjection(false, 0);
         tolerance = 1E-9;
@@ -105,7 +103,6 @@ public final strictfp class TransverseMe
      * @throws TransformException Should never happen.
      */
     @Test
-    @org.junit.Ignore("Missing implementation of the projection derivative.")
     public void testEllipsoidalDerivative() throws TransformException {
         createNormalizedProjection(true, 0);
         tolerance = 1E-9;



Mime
View raw message