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);
}
