sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From rmarec...@apache.org
Subject svn commit: r1740348 - in /sis/branches/ImageDatum/core/sis-referencing/src: main/java/org/apache/sis/referencing/operation/builder/ test/java/org/apache/sis/referencing/operation/builder/
Date Thu, 21 Apr 2016 16:08:38 GMT
Author: rmarechal
Date: Thu Apr 21 16:08:37 2016
New Revision: 1740348

URL: http://svn.apache.org/viewvc?rev=1740348&view=rev
Log:
Add Person correlation and relative tests

Added:
    sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/DoubleDatumShiftGrid.java
Modified:
    sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
    sis/branches/ImageDatum/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilderTest.java

Added: sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/DoubleDatumShiftGrid.java
URL: http://svn.apache.org/viewvc/sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/DoubleDatumShiftGrid.java?rev=1740348&view=auto
==============================================================================
--- sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/DoubleDatumShiftGrid.java
(added)
+++ sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/DoubleDatumShiftGrid.java
Thu Apr 21 16:08:37 2016
@@ -0,0 +1,79 @@
+/*
+ * Copyright 2016 rmarechal.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.referencing.operation.builder;
+
+import javax.measure.quantity.Quantity;
+import javax.measure.unit.Unit;
+import org.apache.sis.referencing.datum.DatumShiftGrid;
+import org.apache.sis.referencing.operation.transform.LinearTransform;
+import org.apache.sis.util.ArgumentChecks;
+
+/**
+     * An implementation of {@link DatumShiftGridFile} which stores the offset values in
{@code double[]} arrays.
+     *
+     * @author  Martin Desruisseaux (Geomatys)
+     * @author  Remi Marechal (Geomatys)
+     * @since   0.7
+     * @version 0.7
+     * @module
+     */
+strictfp class DoubleDatumShiftGrid<C extends Quantity, T extends Quantity> extends
DatumShiftGrid<C, T> {
+
+    /**
+     *
+     */
+    private final double[][] targets;
+
+    /**
+     *
+     */
+    private final double accuracy;
+
+    /**
+     *
+     * @param coordinateUnit
+     * @param sourcePointCoordinateToGrid
+     * @param gridSize
+     * @param translationUnit
+     * @param targets
+     * @param cellPrecision
+     */
+    public DoubleDatumShiftGrid(final Unit<C> coordinateUnit, final LinearTransform
sourcePointCoordinateToGrid,
+            int[] gridSize, /*final boolean isCellValueRatio,*/ final Unit<T> translationUnit,
+            final double[][] targets, double cellPrecision) {
+        super(coordinateUnit, sourcePointCoordinateToGrid, gridSize, false, translationUnit);
+        ArgumentChecks.ensureNonNull("Targets points", targets);
+        ArgumentChecks.ensurePositive("cellPrecision", cellPrecision);
+        this.targets  = targets;
+        this.accuracy = cellPrecision;
+    }
+
+    @Override
+    public int getTranslationDimensions() {
+        return targets.length;
+    }
+
+    @Override
+    public double getCellValue(int dim, int gridX, int gridY) {
+        return targets[dim][gridY * getGridSize()[0] + gridX];
+    }
+
+    @Override
+    public double getCellPrecision() {
+        return accuracy;
+    }
+
+}

Modified: sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
URL: http://svn.apache.org/viewvc/sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java?rev=1740348&r1=1740347&r2=1740348&view=diff
==============================================================================
--- sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
[UTF-8] (original)
+++ sis/branches/ImageDatum/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
[UTF-8] Thu Apr 21 16:08:37 2016
@@ -18,11 +18,16 @@ package org.apache.sis.referencing.opera
 
 import java.io.IOException;
 import java.util.Arrays;
+import javax.measure.quantity.Quantity;
+import javax.measure.unit.Unit;
+
 import org.opengis.geometry.DirectPosition;
 import org.opengis.geometry.MismatchedDimensionException;
+
 import org.apache.sis.io.TableAppender;
 import org.apache.sis.math.Line;
 import org.apache.sis.math.Plane;
+import org.apache.sis.referencing.datum.DatumShiftGrid;
 import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
 import org.apache.sis.referencing.operation.transform.LinearTransform;
@@ -33,7 +38,9 @@ import org.apache.sis.util.resources.Err
 import org.apache.sis.util.Classes;
 import org.apache.sis.util.Debug;
 import org.apache.sis.util.NullArgumentException;
+import org.opengis.referencing.operation.TransformException;
 
+import static java.lang.Math.sqrt;
 
 /**
  * Creates a linear (usually affine) transform which will map approximatively the given source
points to
@@ -79,7 +86,7 @@ public class LinearTransformBuilder {
 
     /**
      * Define the current index of inserted source and target points.
-     * @see #addNoRegularPoints(double[], double[]) 
+     * @see #addNoRegularPoints(double[], double[])
      */
     private int noneRegularPointPosition = 0;
 
@@ -403,7 +410,7 @@ public class LinearTransformBuilder {
      *               {@code coeff[2 + offset] = cy;}
      *               {@code coeff[4 + offset] = c;}
      */
-    private void fitPlane(final double[] grid, final int offset, final double[] coeff) {
+    private double fitPlane(final double[] grid, final int offset, final double[] coeff)
{
         final int width  = gridSize[0];
         final int height = gridSize[1];
         /*
@@ -414,8 +421,8 @@ public class LinearTransformBuilder {
          *           1 + 2 + 3 ... + n    =    n*(n+1)/2              (arithmetic series)
          *        1² + 2² + 3² ... + n²   =    n*(n+0.5)*(n+1)/3
          */
-        double x,y,z, xx,yy, xy, zx,zy;
-        z = zx = zy = 0; // To be computed in the loop.
+        double x,y,z, z2, xx,yy, xy, zx,zy;
+        z = zx = zy = z2 = 0; // To be computed in the loop.
         int n = 0;
         for (int yi=0; yi<height; yi++) {
             for (int xi=0; xi<width; xi++) {
@@ -424,6 +431,7 @@ public class LinearTransformBuilder {
                 if (Double.isNaN(zi))
                     throw new IllegalStateException("The point at coordinate : ("+xi+", "+yi+")
is not referenced.");
                 z  += zi;
+                z2 += zi * zi; //-- prepare correlation computing
                 zx += zi*xi;
                 zy += zi*yi;
                 n++;
@@ -442,18 +450,38 @@ public class LinearTransformBuilder {
          *    ( zx - z*x )  =  cx*(xx - x*x) + cy*(xy - x*y)
          *    ( zy - z*y )  =  cx*(xy - x*y) + cy*(yy - y*y)
          */
-        zx -= z*x/n;
-        zy -= z*y/n;
-        xx -= x*x/n;
-        xy -= x*y/n;
-        yy -= y*y/n;
-        final double den= (xy*xy - xx*yy);
-        final double cy = (zx*xy - zy*xx) / den;
-        final double cx = (zy*xy - zx*yy) / den;
-        final double c  = (z - (cx*x + cy*y)) / n;
+        //-- plan coefficients computing
+        final double pzx = zx - z*x/n;
+        final double pzy = zy - z*y/n;
+        final double pxx = xx - x*x/n;
+        final double pxy = xy - x*y/n;
+        final double pyy = yy - y*y/n;
+        final double den = (pxy * pxy - pxx * pyy);
+        final double cy  = (pzx * pxy - pzy * pxx) / den;
+        final double cx  = (pzy * pxy - pzx * pyy) / den;
+        final double c   = (z - (cx * x + cy * y)) / n;
         coeff[0 + offset] = cx;
         coeff[1 + offset] = cy;
         coeff[2 + offset] = c;
+
+        /*
+         * At this point, the model is computed. Now computes an estimation of the Pearson
+         * correlation coefficient. Note that both the z array and the z computed from the
+         * model have the same average, called sum_z below (the name is not true anymore).
+         *
+         * We do not use double-double arithmetic here since the Pearson coefficient is
+         * for information purpose (quality estimation).
+         */
+        final double mean_x = x / n;
+        final double mean_y = y / n;
+        final double mean_z = z / n;
+
+        final double sum_ds2 = cx*cx*xx + 2*cx*cy*xy + cy*cy*yy
+                                + (cx*mean_x + cy*mean_y) * (cx*(n*mean_x - 2*x) + cy*(n*mean_y
- 2*y));
+        final double sum_dz2 = z2 + mean_z * (n * mean_z - 2 * z);
+        final double sum_dsz = cx*(zx - mean_z*x) + cy*(zy - mean_z*y) - (cy*mean_y + cx*mean_x)*(n*mean_z
- z);
+
+        return Math.min(sum_dsz / sqrt(sum_ds2 * sum_dz2), 1);
     }
 
 
@@ -820,8 +848,8 @@ public class LinearTransformBuilder {
                         //-- means regular grid
                         //-- correlation ??
                         double[] elements = new double[(targetDim + 1)*(sourceDim + 1)];
-                        for (int j=0; j<targets.length; j++) {
-                            fitPlane(targets[j], j * (sourceDim + 1), elements);
+                        for (int j = 0; j < targets.length; j++) {
+                            correlation[j] = fitPlane(targets[j], j * (sourceDim + 1), elements);
                         }
                         matrix.setElements(elements);
                     } else {
@@ -844,6 +872,27 @@ public class LinearTransformBuilder {
     }
 
     /**
+     *
+     * @param translationUnit
+     * @return
+     * @throws TransformException
+     */
+    public DatumShiftGrid getResidus(final Unit<? extends Quantity> translationUnit)
+            throws TransformException {
+        if (gridSize == null)
+            throw new IllegalStateException("Impossible to compute Datum grid from none regular
grid, the grid should be regular.");
+
+        if (transform == null)
+            transform = create();
+        assert isValid();
+
+        final DoubleDatumShiftGrid ddsg = new DoubleDatumShiftGrid(Unit.ONE, MathTransforms.identity(gridSize.length),
+                                                                    gridSize, translationUnit,
+                                                                    targets, COORDS_TOLERANCE);
+        return ddsg;
+    }
+
+    /**
      * Returns the correlation coefficients of the last transform created by {@link #create()},
      * or {@code null} if none. If non-null, the array length is equals to the number of
target
      * dimensions.

Modified: sis/branches/ImageDatum/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilderTest.java
URL: http://svn.apache.org/viewvc/sis/branches/ImageDatum/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilderTest.java?rev=1740348&r1=1740347&r2=1740348&view=diff
==============================================================================
--- sis/branches/ImageDatum/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilderTest.java
[UTF-8] (original)
+++ sis/branches/ImageDatum/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilderTest.java
[UTF-8] Thu Apr 21 16:08:37 2016
@@ -19,12 +19,14 @@ package org.apache.sis.referencing.opera
 import java.util.Random;
 import java.awt.geom.AffineTransform;
 import java.util.Arrays;
+
 import org.opengis.referencing.operation.Matrix;
 import org.apache.sis.geometry.DirectPosition1D;
 import org.apache.sis.geometry.DirectPosition2D;
 import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.TestUtilities;
 import org.apache.sis.test.TestCase;
+
 import org.junit.Assert;
 import org.junit.Test;
 
@@ -124,7 +126,7 @@ public final strictfp class LinearTransf
         assertEquals("m₁₁",  3, m.getElement(1, 1), STRICT);
         assertEquals("m₁₂", -1, m.getElement(1, 2), STRICT);
 
-//        assertArrayEquals("correlation", new double[] {1, 1}, builder.correlation(), STRICT);
+        assertArrayEquals("correlation", new double[] {1, 1}, builder.correlation(), STRICT);
     }
 
     /**
@@ -162,7 +164,7 @@ public final strictfp class LinearTransf
         assertEquals("m₁₁",  3, m.getElement(1, 1), STRICT);
         assertEquals("m₁₂", -1, m.getElement(1, 2), STRICT);
 
-//        assertArrayEquals("correlation", new double[] {1, 1}, builder.correlation(), STRICT);
+        assertArrayEquals("correlation", new double[] {1, 1}, builder.correlation(), STRICT);
 
                                 //--------------//
 
@@ -198,7 +200,7 @@ public final strictfp class LinearTransf
     }
 
     /**
-     * Test method {@link LinearTransformBuilder#setModelTiePoints(int, int, int, int, int,
double[]) .
+     * Test method {@link LinearTransformBuilder#setModelTiePoints(int, int, int, int, int,
double[]) }.
      */
     @Test
     public void setTiePointTest() {
@@ -229,7 +231,7 @@ public final strictfp class LinearTransf
         assertEquals("m₁₁",  3, m.getElement(1, 1), STRICT);
         assertEquals("m₁₂", -1, m.getElement(1, 2), STRICT);
 
-//        assertArrayEquals("correlation", new double[] {1, 1}, builder.correlation(), STRICT);
+        assertArrayEquals("correlation", new double[] {1, 1}, builder.correlation(), STRICT);
 
         //-- regular
         builder  = new LinearTransformBuilder(2, 2);
@@ -249,7 +251,7 @@ public final strictfp class LinearTransf
         assertEquals("m₁₁",  3, m.getElement(1, 1), STRICT);
         assertEquals("m₁₂", -1, m.getElement(1, 2), STRICT);
 
-//        assertArrayEquals("correlation", new double[] {1, 1}, builder.correlation(), STRICT);
+        assertArrayEquals("correlation", new double[] {1, 1}, builder.correlation(), STRICT);
     }
 
 



Mime
View raw message