sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1713002 - in /sis/branches/JDK8/core/sis-referencing/src: main/java/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java test/java/org/apache/sis/referencing/operation/transform/MolodenskyTransformTest.java
Date Fri, 06 Nov 2015 18:00:10 GMT
Author: desruisseaux
Date: Fri Nov  6 18:00:10 2015
New Revision: 1713002

URL: http://svn.apache.org/viewvc?rev=1713002&view=rev
Log:
Initial port of the Molodensky derivative formulas.

Modified:
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java
    sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MolodenskyTransformTest.java

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java?rev=1713002&r1=1713001&r2=1713002&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MolodenskyTransform.java
[UTF-8] Fri Nov  6 18:00:10 2015
@@ -30,6 +30,7 @@ import org.opengis.referencing.operation
 import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.MathTransformFactory;
 import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.datum.DefaultEllipsoid;
 import org.apache.sis.internal.referencing.provider.Molodensky;
 import org.apache.sis.internal.referencing.provider.AbridgedMolodensky;
@@ -516,32 +517,94 @@ public class MolodenskyTransform extends
         final double sinφ  = sin(φ);
         final double cosφ  = cos(φ);
         final double sin2φ = sinφ * sinφ;
-              double ρden  = 1 - excentricitySquared * sin2φ;               // Denominator
of ρ (completed later)
-        final double νden  = sqrt(ρden);                                    // Denominator
of ν
-        double ρ = semiMajor * (1 - excentricitySquared) / (ρden *= νden);  // (also complete
calculation of ρden)
-        double ν = semiMajor / νden;
+        final double ν2den = 1 - excentricitySquared*sin2φ;                 // Square of
the denominator of ν
+        final double νden  = sqrt(ν2den);                                   // Denominator
of ν
+        final double ρden  = ν2den * νden;                                  // Denominator
of ρ
+        double ρ = semiMajor * (1 - excentricitySquared) / ρden;            // Other notation:
Rm = ρ
+        double ν = semiMajor / νden;                                        // Other notation:
Rn = ν
         double t = Δfmod * 2;                                               // A term in
the calculation of Δφ
         if (!isAbridged) {
             ρ += h;
             ν += h;
-            t = t*(0.5/νden + 0.5/ρden) + Δa*excentricitySquared/νden;
+            t = t*(0.5/νden + 0.5/ρden)                 // = Δf⋅[ν⋅(b/a) + ρ⋅(a/b)]
    (without the +h in ν and ρ)
+                    + Δa*excentricitySquared/νden;      // = Δa⋅[ℯ²⋅ν/a]
         }
-        final double tXY = tY*sinλ + tX*cosλ;
-        dstPts[dstOff++] = λ + ANGULAR_SCALE * (tY*cosλ - tX*sinλ) / (ν*cosφ);
-        dstPts[dstOff++] = φ + ANGULAR_SCALE * ((t*cosφ - tXY)*sinφ + tZ*cosφ) / ρ;
-        if (isTarget3D) {
-            t = Δfmod * sin2φ;                                              // A term in
the calculation of Δh
-            double d = Δa;
-            if (!isAbridged) {
-                t /= νden;
-                d *= νden;
+        final double spcλ = tY*sinλ + tX*cosλ;                      // "spc" stands for
"sin plus cos"
+        final double cmsλ = tY*cosλ - tX*sinλ;                      // "cms" stands for
"cos minus sin"
+        final double cmsφ = (tZ + t*sinφ)*cosφ - spcλ*sinφ;
+        final double scaleX = ANGULAR_SCALE / (ν*cosφ);
+        final double scaleY = ANGULAR_SCALE / ρ;
+        if (dstPts != null) {
+            dstPts[dstOff++] = λ + (cmsλ * scaleX);
+            dstPts[dstOff++] = φ + (cmsφ * scaleY);
+            if (isTarget3D) {
+                double t1 = Δfmod * sin2φ;          // A term in the calculation of Δh
+                double t2 = Δa;
+                if (!isAbridged) {
+                    t1 /= νden;                     // = Δf⋅(b/a)⋅ν⋅sin²φ
+                    t2 *= νden;                     // = Δa⋅(a/ν)
+                }
+                dstPts[dstOff++] = h + spcλ*cosφ + tZ*sinφ + t1 - t2;
             }
-            dstPts[dstOff++] = h + tXY*cosφ + tZ*sinφ + t - d;
         }
         if (!derivate) {
             return null;
         }
-        return null;  // TODO
+        /*
+         * At this point the (Abridged) Molodensky transformation is finished.
+         * Code below this point is only for computing the derivative, if requested.
+         * Note: variable names do not necessarily tell all the terms that they contain.
+         */
+        final Matrix matrix   = Matrices.createDiagonal(getTargetDimensions(), getSourceDimensions());
+        final double sinφcosφ = sinφ * cosφ;
+        final double dν       = excentricitySquared*sinφcosφ / ν2den;
+        final double dν3ρ     = 3*dν * (1 - excentricitySquared) / ν2den;
+        //    double dXdλ     = spcλ;
+        final double dYdλ     = cmsλ * sinφ;
+        final double dZdλ     = cmsλ * cosφ;
+              double dXdφ     = dYdλ / cosφ;
+              double dYdφ     = -tZ*sinφ - cosφ*spcλ  +  t*(1 - 2*sin2φ);
+              double dZdφ     =  tZ*cosφ - sinφ*spcλ;
+        if (isAbridged) {
+            /*
+             *   Δfmod  =  (a⋅Δf) + (f⋅Δa)
+             *   t      =  2⋅Δfmod
+             *   dXdh   =  0  so no need to set the matrix element.
+             *   dYdh   =  0  so no need to set the matrix element.
+             */
+            dXdφ -= cmsλ * dν;
+            dYdφ -= cmsφ * dν3ρ;
+            dZdφ += t*cosφ*sinφ;
+        } else {
+            /*
+             *   Δfmod  =  b⋅Δf
+             *   t      =  Δf⋅[ν⋅(b/a) + ρ⋅(a/b)]    (real ν and ρ, without +
h)
+             *   ν         is actually ν + h
+             *   ρ         is actually ρ + h
+             */
+            final double dρ = dν3ρ * νden * (semiMajor / ρ);    // Reminder: that ρ
contains a h term.
+            dXdφ -= dν * cmsλ * semiMajor / (νden*ν);           // Reminder: that ν
contains a h term.
+            dYdφ -= dρ * dZdφ - (Δfmod*(dν*2/(1 - excentricitySquared) + (1 + 1/ν2den)*(dν
- dρ))
+                                  + Δa*(dν + 1)*excentricitySquared) * sinφcosφ / νden;
+            if (isSource3D) {
+                final double dXdh =  cmsλ / ν;
+                final double dYdh = -cmsφ / ρ;
+                matrix.setElement(0, 2, -dXdh * scaleX);
+                matrix.setElement(1, 2, +dYdh * scaleY);
+            }
+            final double t1 = Δfmod * (dν*sin2φ + 2*sinφcosφ);
+            final double t2 = Δa * dν;
+            dZdφ += t1/νden + t2*νden;
+        }
+        matrix.setElement(0, 0, 1 - spcλ * scaleX);
+        matrix.setElement(1, 1, 1 + dYdφ * scaleY);
+        matrix.setElement(0, 1,   + dXdφ * scaleX);
+        matrix.setElement(1, 0,   - dYdλ * scaleY);
+        if (isTarget3D) {
+            matrix.setElement(2, 0, dZdλ);
+            matrix.setElement(2, 1, dZdφ);
+        }
+        return matrix;
     }
 
     /**
@@ -600,19 +663,19 @@ public class MolodenskyTransform extends
             final double sinφ  = sin(φ);
             final double cosφ  = cos(φ);
             final double sin2φ = sinφ * sinφ;
-                  double ρden  = 1 - excentricitySquared * sin2φ;
-            final double νden  = sqrt(ρden);
-            double ρ = semiMajor * (1 - excentricitySquared) / (ρden *= νden);
-            double ν = semiMajor / νden;
-            double t = Δfmod * 2;
+                  double ρden  = 1 - excentricitySquared * sin2φ;               // Denominator
of ρ (completed later)
+            final double νden  = sqrt(ρden);                                    // Denominator
of ν
+            double ρ = semiMajor * (1 - excentricitySquared) / (ρden *= νden);  // (also
complete calculation of ρden)
+            double ν = semiMajor / νden;                                        // Other
notation: Rm = ρ and Rn = ν
+            double t = Δfmod * 2;                                               // A term
in the calculation of Δφ
             if (!isAbridged) {
                 ρ += h;
                 ν += h;
                 t = t*(0.5/νden + 0.5/ρden) + Δa*excentricitySquared/νden;
             }
-            final double tXY = tY*sinλ + tX*cosλ;
+            final double spcλ = tY*sinλ + tX*cosλ;
             dstPts[dstOff++] = λ + ANGULAR_SCALE * (tY*cosλ - tX*sinλ) / (ν*cosφ);
-            dstPts[dstOff++] = φ + ANGULAR_SCALE * ((t*cosφ - tXY)*sinφ + tZ*cosφ) /
ρ;
+            dstPts[dstOff++] = φ + ANGULAR_SCALE * ((t*cosφ - spcλ)*sinφ + tZ*cosφ)
/ ρ;
             if (isTarget3D) {
                 t = Δfmod * sin2φ;
                 double d = Δa;
@@ -620,7 +683,7 @@ public class MolodenskyTransform extends
                     t /= νden;
                     d *= νden;
                 }
-                dstPts[dstOff++] = h + tXY*cosφ + tZ*sinφ + t - d;
+                dstPts[dstOff++] = h + spcλ*cosφ + tZ*sinφ + t - d;
             }
             srcOff += srcInc;
             dstOff += dstInc;

Modified: sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MolodenskyTransformTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MolodenskyTransformTest.java?rev=1713002&r1=1713001&r2=1713002&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MolodenskyTransformTest.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MolodenskyTransformTest.java
[UTF-8] Fri Nov  6 18:00:10 2015
@@ -73,6 +73,36 @@ public final strictfp class MolodenskyTr
         zTolerance = GeocentricTranslationTest.precision(3);        // Required precision
for h
         zDimension = new int[] {2};                                 // Dimension of h where
to apply zTolerance
         assertFalse(transform.isIdentity());
+        validate();
+    }
+
+    /**
+     * Tests the derivative of the Abridged Molodensky transformation.
+     *
+     * @throws FactoryException if an error occurred while creating a transform step.
+     * @throws TransformException if the transformation failed.
+     */
+    @Test
+    public void testAbridgedMolodenskyDerivative() throws FactoryException, TransformException
{
+        create(true);
+        verifyDerivative( 0,  0,  0);
+        verifyDerivative(-3, 30,  7);
+        verifyDerivative(+6, 60, 20);
+    }
+
+    /**
+     * Tests the derivative of the Molodensky transformation.
+     *
+     * @throws FactoryException if an error occurred while creating a transform step.
+     * @throws TransformException if the transformation failed.
+     */
+    @Test
+    @DependsOnMethod("testAbridgedMolodenskyDerivative")
+    public void testMolodenskyDerivative() throws FactoryException, TransformException {
+        create(false);
+        verifyDerivative( 0,  0,  0);
+        verifyDerivative(-3, 30,  7);
+        verifyDerivative(+6, 60, 20);
     }
 
     /**
@@ -88,10 +118,9 @@ public final strictfp class MolodenskyTr
      * @throws TransformException if the transformation failed.
      */
     @Test
+    @DependsOnMethod("testAbridgedMolodenskyDerivative")
     public void testAbridgedMolodensky() throws FactoryException, TransformException {
-        isDerivativeSupported = false;
         create(true);
-        validate();
         final double[] sample   = GeocentricTranslationTest.samplePoint(1);
         final double[] expected = GeocentricTranslationTest.samplePoint(5);
         isInverseTransformSupported = false;
@@ -118,11 +147,9 @@ public final strictfp class MolodenskyTr
      * @throws TransformException if the transformation failed.
      */
     @Test
-    @DependsOnMethod("testAbridgedMolodensky")
+    @DependsOnMethod({"testAbridgedMolodensky", "testMolodenskyDerivative"})
     public void testMolodensky() throws FactoryException, TransformException {
-        isDerivativeSupported = false;
         create(false);
-        validate();
         final double[] sample   = GeocentricTranslationTest.samplePoint(1);
         final double[] expected = GeocentricTranslationTest.samplePoint(4);
         isInverseTransformSupported = false;
@@ -147,7 +174,6 @@ public final strictfp class MolodenskyTr
     @Test
     @DependsOnMethod("testMolodensky")
     public void testRandomPoints() throws FactoryException, TransformException {
-        isDerivativeSupported = false;
         create(false);
         tolerance  = Formulas.LINEAR_TOLERANCE * 3;     // To be converted in degrees by
ToleranceModifier.GEOGRAPHIC
         zTolerance = Formulas.LINEAR_TOLERANCE * 2;



Mime
View raw message