sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1787227 - in /sis/branches/JDK8/core: sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ sis-utility/src/main/java/org/apache/sis/math/...
Date Thu, 16 Mar 2017 18:32:29 GMT
Author: desruisseaux
Date: Thu Mar 16 18:32:29 2017
New Revision: 1787227

URL: http://svn.apache.org/viewvc?rev=1787227&view=rev
Log:
LinearTransformBuilder should be able to take advantage of the knownledge that source positions are distributed on a grid (when this is the case).

Modified:
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java
    sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/ArrayVector.java
    sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/CompoundDirectPositions.java
    sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Line.java
    sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java
    sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
    sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
    sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties
    sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/math/PlaneTest.java

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -16,12 +16,16 @@
  */
 package org.apache.sis.referencing.operation.builder;
 
+import java.util.Map;
+import java.util.Arrays;
 import java.io.IOException;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.geometry.MismatchedDimensionException;
+import org.opengis.geometry.coordinate.Position;
 import org.apache.sis.io.TableAppender;
 import org.apache.sis.math.Line;
 import org.apache.sis.math.Plane;
+import org.apache.sis.math.Vector;
 import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
 import org.apache.sis.referencing.operation.transform.LinearTransform;
@@ -34,21 +38,26 @@ import org.apache.sis.util.Debug;
 
 
 /**
- * Creates a linear (usually affine) transform which will map approximatively the given source points to
- * the given target points. The transform coefficients are determined using a <cite>least squares</cite>
- * estimation method, with the assumption that source points are precise and all uncertainty is in the
- * target points.
+ * Creates an affine transform which will map approximatively the given source positions to the given target positions.
+ * In many cases, the <em>source</em> positions are grid indices and the <em>target</em> positions are geographic or
+ * projected coordinates, but this is not mandatory. If the source positions are known to be grid indices,
+ * then a builder created by the {@link #LinearTransformBuilder(int...)} constructor will be more efficient.
+ * Otherwise a builder created by the {@link #LinearTransformBuilder()} constructor will be able to handle
+ * randomly distributed coordinates.
+ *
+ * <p>The transform coefficients are determined using a <cite>least squares</cite> estimation method,
+ * with the assumption that source positions are exact and all the uncertainty is in the target positions.</p>
  *
  * <div class="note"><b>Implementation note:</b>
  * The quantity that current implementation tries to minimize is not strictly the squared Euclidian distance.
- * The current implementation rather processes each target dimension independently, which may not give the same
- * result than if we tried to minimize the squared Euclidian distances by taking all dimensions in account together.
- * This algorithm may change in future SIS versions.
+ * The current implementation rather processes each <em>target</em> dimension independently, which may not give
+ * the same result than if we tried to minimize the squared Euclidian distances by taking all target dimensions
+ * in account together. This algorithm may change in future SIS versions.
  * </div>
  *
  * @author  Martin Desruisseaux (Geomatys)
  * @since   0.5
- * @version 0.5
+ * @version 0.8
  * @module
  *
  * @see LinearTransform
@@ -57,141 +66,365 @@ import org.apache.sis.util.Debug;
  */
 public class LinearTransformBuilder {
     /**
-     * The arrays of source ordinate values, for example (x[], y[]).
-     * This is {@code null} if not yet specified.
+     * Number of grid columns and rows, or {@code null} if the coordinates are not distributed on a regular grid.
+     * If the grid size is known, then the {@link #sources} coordinates do not need to be specified.
+     */
+    private final int[] gridSize;
+
+    /**
+     * The arrays of source ordinate values. Accessed with indices in that order: {@code sources[dimension][point]}.
+     * This layout allows to create only a few (typically two) large arrays instead of a multitude of small arrays.
+     * Example: {x[], y[]}.
+     *
+     * <p>In the special case where {@link #gridSize} is non-null, then this array does not need to be specified
+     * and can be {@code null}. In such case, the source coordinates are implicitly:</p>
+     *
+     * <blockquote>
+     * (0,0), (1,0), (2,0), (3,0) … ({@link #gridSize}[0]-1, 0),<br>
+     * (0,1), (1,1), (2,1), (3,1) … ({@link #gridSize}[0]-1, 1),<br>
+     * (0,2), (1,2), (2,2), (3,2) … ({@link #gridSize}[0]-1, 2),<br>
+     * (0,{@link #gridSize}[1]-1) … ({@link #gridSize}[0]-1, {@link #gridSize}[1]-1).
+     * </blockquote>
      */
     private double[][] sources;
 
     /**
-     * The arrays of target ordinate values, for example (x[], y[], z[]).
+     * The arrays of target ordinate values. Accessed with indices in that order: {@code targets[dimension][point]}.
+     * This layout allows to create only a few (typically two) large arrays instead of a multitude of small arrays.
+     * Example: {x[], y[], z[]}.
      * This is {@code null} if not yet specified.
      */
     private double[][] targets;
 
     /**
+     * The product of all {@link #gridSize} values, or 0 if none if {@link #gridSize} is null.
+     * If non-zero, then this is the length of {@link #targets} arrays to create.
+     */
+    private final int gridLength;
+
+    /**
+     * Number of valid positions in the {@link #sources} or {@link #targets} arrays.
+     * Note that the "valid" positions may contain {@link Double#NaN} ordinate values.
+     */
+    private int numPoints;
+
+    /**
      * The transform created by the last call to {@link #create()}.
+     * This is reset to {@code null} when coordinates are modified.
      */
-    private LinearTransform transform;
+    private transient LinearTransform transform;
 
     /**
      * An estimation of the Pearson correlation coefficient for each target dimension.
-     * This is {@code null} if not yet specified.
+     * This is {@code null} if not yet computed.
      */
-    private double[] correlation;
+    private transient double[] correlation;
 
     /**
-     * Creates a new linear transform builder.
+     * Creates a new linear transform builder for randomly distributed positions.
+     *
+     * <div class="note"><b>Tip:</b>
+     * if the source coordinates are grid indices, then
+     * the {@link #LinearTransformBuilder(int, int)} constructor will create a more efficient builder.
+     * </div>
      */
     public LinearTransformBuilder() {
+        gridSize = null;
+        gridLength = 0;
+    }
+
+    /**
+     * Creates a new linear transform builder for source positions distributed on a regular grid.
+     * This constructor notifies {@code LinearTransformBuilder} that ordinate values of all source positions will
+     * be integers in the [0 … {@code gridSize[0]}-1] range for the first dimension (typically column indices),
+     * in the [0 … {@code gridSize[1]}-1] range for the second dimension (typically row indices), <i>etc.</i>
+     * The dimension of all source positions is the length of the given {@code gridSize} array.
+     *
+     * <p>An empty array is equivalent to invoking the no-argument constructor,
+     * i.e. no restriction is put on the source coordinates.</p>
+     *
+     * @param  gridSize  the number of integer ordinate values in each grid dimension.
+     * @throws IllegalArgumentException if a grid size is not strictly positive, or if the product
+     *         of all values (∏{@code gridSize}) is greater than {@link Integer#MAX_VALUE}.
+     *
+     * @since 0.8
+     */
+    public LinearTransformBuilder(int... gridSize) {
+        ArgumentChecks.ensureNonNull("gridSize", gridSize);
+        if (gridSize.length == 0) {
+            this.gridSize = null;
+            this.gridLength = 0;
+        } else {
+            gridSize = gridSize.clone();
+            long length = 1;
+            for (int s : gridSize) {
+                ArgumentChecks.ensureStrictlyPositive("gridSize", s);
+                length = Math.multiplyExact(length, s);
+            }
+            if (length > Integer.MAX_VALUE) {
+                throw new IllegalArgumentException(Errors.format(Errors.Keys.ValueOutOfRange_4,
+                        "∏gridSize", 1, Integer.MAX_VALUE, length));
+            }
+            this.gridSize = gridSize;
+            gridLength = (int) length;
+        }
     }
 
     /**
-     * Extracts the ordinate values of the given points into separated arrays, one for each dimension.
+     * Returns the offset where to store a target position for the given source position.
+     * This method should be invoked only when this {@code LinearTransformBuilder} has
+     * been constructed for a grid of known size.
      *
-     * @param  points     the points from which to extract the ordinate values.
-     * @param  dimension  the expected number of dimensions.
+     * @throws IllegalArgumentException if an ordinate value is illegal.
      */
-    private static double[][] toArrays(final DirectPosition[] points, final int dimension) {
-        final int length = points.length;
-        final double[][] ordinates = new double[dimension][length];
-        for (int j=0; j<length; j++) {
-            final DirectPosition p = points[j];
-            final int d = p.getDimension();
-            if (d != dimension) {
-                throw new MismatchedDimensionException(Errors.format(
-                        Errors.Keys.MismatchedDimension_3, "points[" + j + ']', dimension, d));
+    private int offset(final DirectPosition src) {
+        int offset = 1;
+        for (int j = gridSize.length; j != 0;) {
+            final int s = gridSize[--j];
+            final double ordinate = src.getOrdinate(j);
+            final int i = (int) ordinate;
+            if (i != ordinate) {
+                throw new IllegalArgumentException(Errors.format(Errors.Keys.NotAnInteger_1, ordinate));
             }
-            for (int i=0; i<dimension; i++) {
-                ordinates[i][j] = p.getOrdinate(i);
+            if (i < 0 || i >= s) {
+                throw new IllegalArgumentException(Errors.format(Errors.Keys.ValueOutOfRange_4, "source", 0, s-1, i));
             }
+            offset = offset * s + i;
         }
-        return ordinates;
+        return offset;
+    }
+
+    /**
+     * Builds the exception message for an unexpected position dimension. This method assumes
+     * that positions are stored in this builder as they are read from user-provided collection.
+     */
+    private String mismatchedDimension(final String name, final int expected, final int actual) {
+        return Errors.format(Errors.Keys.MismatchedDimension_3, name + '[' + numPoints + ']', expected, actual);
     }
 
     /**
-     * Sets the source points. The number of points shall be the same than the number of target points.
+     * Returns the direct position of the given position, or {@code null} if none.
+     */
+    private static DirectPosition position(final Position p) {
+        return (p != null) ? p.getDirectPosition() : null;
+    }
+
+    /**
+     * Sets the source and target positions, overwriting any previous setting. The source positions are the keys in
+     * the given map, and the target positions are the associated values in the map. The map should not contain two
+     * entries with the same source position. Null positions are silently ignored.
+     *
+     * <p>All source positions shall have the same number of dimensions (the <cite>source dimension</cite>),
+     * and all target positions shall have the same number of dimensions (the <cite>target dimension</cite>).
+     * However the source dimension does not need to be the same the target dimension.
+     * Apache SIS currently supports only one- or two-dimensional source positions,
+     * together with arbitrary target dimension.</p>
+     *
+     * <p>If this builder has been created with the {@link #LinearTransformBuilder(int...)} constructor,
+     * then the ordinate values of all source positions shall be integers in the [0 … {@code gridSize[0]}-1]
+     * range for the first dimension (typically column indices), in the [0 … {@code gridSize[1]}-1] range for
+     * the second dimension (typically row indices), <i>etc</i>. This constraint does not apply for builders
+     * created with the {@link #LinearTransformBuilder()} constructor.</p>
      *
-     * <p><b>Limitation:</b> in current implementation, the source points must be one or two-dimensional.
-     * This restriction may be removed in a future SIS version.</p>
+     * @param  sourceToTarget  a map of source positions to target positions.
+     *         Source positions are assumed precise and target positions are assumed uncertain.
+     * @throws MismatchedDimensionException if some positions do not have the expected number of dimensions.
+     *
+     * @since 0.8
+     */
+    public void setPositions(final Map<? extends Position, ? extends Position> sourceToTarget)
+            throws MismatchedDimensionException
+    {
+        ArgumentChecks.ensureNonNull("sourceToTarget", sourceToTarget);
+        pendingSources = null;
+        pendingTargets = null;
+        transform   = null;
+        correlation = null;
+        sources     = null;
+        targets     = null;
+        numPoints   = 0;
+        int srcDim  = 0;
+        int tgtDim  = 0;
+        for (final Map.Entry<? extends Position, ? extends Position> entry : sourceToTarget.entrySet()) {
+            final DirectPosition src = position(entry.getKey());   if (src == null) continue;
+            final DirectPosition tgt = position(entry.getValue()); if (tgt == null) continue;
+            /*
+             * The first time that we get a non-null source and target coordinate, allocate the arrays.
+             * The sources arrays are allocated only if the source coordiantes are randomly distributed.
+             */
+            if (targets == null) {
+                tgtDim = tgt.getDimension();
+                if (tgtDim <= 0) {
+                    throw new MismatchedDimensionException(mismatchedDimension("target", 2, tgtDim));
+                }
+                if (gridSize == null) {
+                    srcDim = src.getDimension();
+                    if (srcDim <= 0) {
+                        throw new MismatchedDimensionException(mismatchedDimension("source", 2, srcDim));
+                    }
+                    final int length = sourceToTarget.size();
+                    sources = new double[srcDim][length];
+                    targets = new double[tgtDim][length];
+                } else {
+                    srcDim  = gridSize.length;
+                    targets = new double[tgtDim][gridLength];
+                    for (final double[] r : targets) {
+                        Arrays.fill(r, Double.NaN);
+                    }
+                }
+            }
+            /*
+             * Verify that the source and target coordinates have the expected number of dimensions before to store
+             * the coordinates. If the grid size is known, we do not need to store the source coordinates. Instead,
+             * we compute its index in the fixed-size target arrays.
+             */
+            int d;
+            if ((d = src.getDimension()) != srcDim) throw new MismatchedDimensionException(mismatchedDimension("source", srcDim, d));
+            if ((d = tgt.getDimension()) != tgtDim) throw new MismatchedDimensionException(mismatchedDimension("target", tgtDim, d));
+            int index;
+            if (gridSize != null) {
+                index = offset(src);
+            } else {
+                index = numPoints;
+                for (int i=0; i<srcDim; i++) {
+                    sources[i][index] = src.getOrdinate(i);
+                }
+            }
+            for (int i=0; i<tgtDim; i++) {
+                targets[i][index] = tgt.getOrdinate(i);
+            }
+            numPoints++;
+        }
+    }
+
+    /**
+     * Sets the source points, overwriting any previous setting. The number of source points will need to be the same
+     * than the number of {@linkplain #setTargetPoints target points} when the {@link #create()} method will be invoked.
+     * In current Apache SIS implementation, the source points must be one or two-dimensional.
+     *
+     * <p>If this builder has been created with the {@link #LinearTransformBuilder(int, int)} constructor,
+     * then all given points must be two-dimensional and all ordinate values must be integers in the
+     * [0 … <var>width</var>-1] or [0 … <var>height</var>-1] range for the first and second dimension
+     * respectively. This constraint does not apply if this builder has been created with the
+     * {@link #LinearTransformBuilder()} constructor.</p>
+     *
+     * <p>It is caller's responsibility to ensure that no source point is duplicated.
+     * If the same source point is repeated twice, then {@code LinearTransformBuilder} behavior is undefined.</p>
      *
      * @param  points  the source points, assumed precise.
      * @throws MismatchedDimensionException if at least one point does not have the expected number of dimensions.
+     *
+     * @deprecated Replaced by {@link #setPositions(Map)}.
      */
+    @Deprecated
     public void setSourcePoints(final DirectPosition... points) throws MismatchedDimensionException {
         ArgumentChecks.ensureNonNull("points", points);
-        if (points.length != 0) {
-            sources = toArrays(points, points[0].getDimension() == 1 ? 1 : 2);
-        } else {
-            sources = null;
-        }
         transform   = null;
         correlation = null;
+        sources     = null;
+        targets     = null;
+        numPoints   = 0;
+        pendingSources = points.clone();
     }
 
     /**
-     * Sets the target points. The number of points shall be the same than the number of source points.
+     * Sets the target points, overwriting any previous setting. The number of target points will need to be the same
+     * than the number of {@linkplain #setSourcePoints source points} when the {@link #create()} method will be invoked.
      * Target points can have any number of dimensions (not necessarily 2), but all points shall have
      * the same number of dimensions.
      *
      * @param  points  the target points, assumed uncertain.
      * @throws MismatchedDimensionException if not all points have the same number of dimensions.
+     *
+     * @deprecated Replaced by {@link #setPositions(Map)}.
      */
+    @Deprecated
     public void setTargetPoints(final DirectPosition... points) throws MismatchedDimensionException {
         ArgumentChecks.ensureNonNull("points", points);
-        if (points.length != 0) {
-            targets = toArrays(points, points[0].getDimension());
-        } else {
-            targets = null;
-        }
         transform   = null;
         correlation = null;
+        sources     = null;
+        targets     = null;
+        numPoints   = 0;
+        pendingTargets = points.clone();
     }
 
-    /*
-     * No getters yet because we did not determined what they should return.
-     * Array? Collection? Map<source,target>?
-     */
+    @Deprecated
+    private transient DirectPosition[] pendingSources, pendingTargets;
+
+    @Deprecated
+    private void processPendings() {
+        if (pendingSources != null || pendingTargets != null) {
+            if (pendingSources == null || pendingTargets == null) {
+                throw new IllegalStateException(Errors.format(
+                        Errors.Keys.MissingValueForProperty_1, (pendingSources == null) ? "sources" : "targets"));
+            }
+            final int length = pendingSources.length;
+            if (pendingTargets.length != length) {
+                throw new IllegalStateException(Errors.format(Errors.Keys.MismatchedArrayLengths));
+            }
+            final Map<DirectPosition,DirectPosition> sourceToTarget = new java.util.HashMap<>(length);
+            for (int i=0; i<length; i++) {
+                sourceToTarget.put(pendingSources[i], pendingTargets[i]);
+            }
+            setPositions(sourceToTarget);
+        }
+    }
 
     /**
-     * Creates a linear transform approximation from the source points to the target points.
-     * This method assumes that source points are precise and all uncertainty is in the target points.
+     * Creates a linear transform approximation from the source positions to the target positions.
+     * This method assumes that source positions are precise and that all uncertainties are in the
+     * target positions.
      *
      * @return the fitted linear transform.
+     * @throws IllegalStateException if the source or target points have not be specified
+     *         or if those two sets do not have the same number of points.
      */
     public LinearTransform create() {
         if (transform == null) {
+            processPendings();
             final double[][] sources = this.sources;                    // Protect from changes.
             final double[][] targets = this.targets;
-            if (sources == null || targets == null) {
-                throw new IllegalStateException(Errors.format(
-                        Errors.Keys.MissingValueForProperty_1, (sources == null) ? "sources" : "targets"));
+            if (targets == null) {
+                throw new IllegalStateException(Errors.format(Errors.Keys.MissingValueForProperty_1, "sourceToTarget"));
             }
-            final int sourceDim = sources.length;
+            final int sourceDim = (sources != null) ? sources.length : gridSize.length;
             final int targetDim = targets.length;
             correlation = new double[targetDim];
             final MatrixSIS matrix = Matrices.createZero(targetDim + 1, sourceDim + 1);
             matrix.setElement(targetDim, sourceDim, 1);
-            switch (sourceDim) {
-                case 1: {
-                    final Line line = new Line();
-                    for (int j=0; j < targets.length; j++) {
-                        correlation[j] = line.fit(sources[0], targets[j]);
+            for (int j=0; j < targetDim; j++) {
+                final double c;
+                switch (sourceDim) {
+                    case 1: {
+                        final Line line = new Line();
+                        if (sources != null) {
+                            c = line.fit(sources[0], targets[j]);
+                        } else {
+                            c = line.fit(Vector.createSequence(0, 1, gridSize[0]),
+                                         Vector.create(targets[j], false));
+                        }
                         matrix.setElement(j, 0, line.slope());
                         matrix.setElement(j, 1, line.y0());
+                        break;
                     }
-                    break;
-                }
-                case 2: {
-                    final Plane plan = new Plane();
-                    for (int j=0; j < targets.length; j++) {
-                        correlation[j] = plan.fit(sources[0], sources[1], targets[j]);
+                    case 2: {
+                        final Plane plan = new Plane();
+                        if (sources != null) {
+                            c = plan.fit(sources[0], sources[1], targets[j]);
+                        } else {
+                            c = plan.fit(gridSize[0], gridSize[1], Vector.create(targets[j], false));
+                        }
                         matrix.setElement(j, 0, plan.slopeX());
                         matrix.setElement(j, 1, plan.slopeY());
                         matrix.setElement(j, 2, plan.z0());
+                        break;
+                    }
+                    default: {
+                        throw new UnsupportedOperationException(Errors.format(Errors.Keys.ExcessiveNumberOfDimensions_1, sourceDim));
                     }
-                    break;
                 }
-                default: throw new AssertionError(sourceDim);   // Should have been verified by setSourcePoints(…) method.
+                correlation[j] = c;
             }
             transform = MathTransforms.linear(matrix);
         }
@@ -227,13 +460,9 @@ public class LinearTransformBuilder {
             buffer.append(':').append(lineSeparator);
             final TableAppender table = new TableAppender(buffer, " ");
             table.setMultiLinesCells(true);
-            table.append(Matrices.toString(transform.getMatrix()));
-            table.nextColumn();
-            table.append(lineSeparator);
-            table.append("  ");
-            table.append(Vocabulary.format(Vocabulary.Keys.Correlation));
-            table.append(" =");
-            table.nextColumn();
+            table.append(Matrices.toString(transform.getMatrix())).nextColumn();
+            table.append(lineSeparator).append("  ")
+                 .append(Vocabulary.format(Vocabulary.Keys.Correlation)).append(" =").nextColumn();
             table.append(Matrices.create(correlation.length, 1, correlation).toString());
             try {
                 table.flush();

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -31,7 +31,7 @@
  *
  * @author  Martin Desruisseaux (Geomatys)
  * @version 0.5
- * @since   0.5
+ * @since   0.8
  * @module
  */
 package org.apache.sis.referencing.operation.builder;

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -55,6 +55,9 @@ import org.opengis.parameter.ParameterVa
  * current implementation can only be backed by the Transverse Mercator projection,
  * but future versions could apply to some other projections if needed.</div>
  *
+ * <p>Examples of CRS using this projection are <cite>WGS 84 / UTM grid system</cite>
+ * EPSG:32600 (northern hemisphere) and EPSG:32700 (southern hemisphere).</p>
+ *
  * @author  Martin Desruisseaux (Geomatys)
  * @since   0.8
  * @version 0.8

Modified: sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/ArrayVector.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/ArrayVector.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/ArrayVector.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/ArrayVector.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -210,7 +210,7 @@ abstract class ArrayVector<E extends Num
     /**
      * A vector backed by an array of type {@code double[]}.
      */
-    private static final class Doubles extends ArrayVector<Double> {
+    static final class Doubles extends ArrayVector<Double> {
         /** For cross-version compatibility. */
         private static final long serialVersionUID = -2900375382498345812L;
 

Modified: sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/CompoundDirectPositions.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/CompoundDirectPositions.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/CompoundDirectPositions.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/CompoundDirectPositions.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -35,14 +35,14 @@ import org.apache.sis.util.resources.Err
  *
  * @author  Martin Desruisseaux (Geomatys)
  * @since   0.5
- * @version 0.5
+ * @version 0.8
  * @module
  */
 final class CompoundDirectPositions implements DirectPosition, Iterable<DirectPosition>, Iterator<DirectPosition> {
     /**
      * The arrays of ordinate values, for example (x[], y[], z[]).
      */
-    private final double[][] ordinates;
+    private final Vector[] ordinates;
 
     /**
      * Length of all ordinate values minus one.
@@ -57,15 +57,16 @@ final class CompoundDirectPositions impl
     /**
      * Wraps the given array of ordinate values.
      */
-    CompoundDirectPositions(final double[]... ordinates) {
+    CompoundDirectPositions(final Vector... ordinates) {
         this.ordinates = ordinates;
-        final int length = ordinates[0].length;
+        final int length = ordinates[0].size();
         for (int i=1; i<ordinates.length; i++) {
-            if (ordinates[i].length != length) {
+            if (ordinates[i].size() != length) {
                 throw new IllegalArgumentException(Errors.format(Errors.Keys.MismatchedArrayLengths));
             }
         }
-        last = length - 1;
+        last  = length - 1;
+        index = length;
     }
 
     /**
@@ -75,6 +76,9 @@ final class CompoundDirectPositions impl
      */
     @Override
     public Iterator<DirectPosition> iterator() {
+        if (hasNext()) {
+            throw new UnsupportedOperationException(Errors.format(Errors.Keys.CanIterateOnlyOnce));
+        }
         index = -1;
         return this;
     }
@@ -131,7 +135,7 @@ final class CompoundDirectPositions impl
      */
     @Override
     public double getOrdinate(final int dimension) {
-        return ordinates[dimension][index];
+        return ordinates[dimension].doubleValue(index);
     }
 
     /**

Modified: sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Line.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Line.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Line.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Line.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -22,6 +22,7 @@ import org.opengis.geometry.MismatchedDi
 import org.apache.sis.internal.util.DoubleDouble;
 import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.util.resources.Errors;
+import org.apache.sis.util.ArgumentChecks;
 
 import static java.lang.Double.*;
 
@@ -42,7 +43,7 @@ import static java.lang.Double.*;
  *
  * @author  Martin Desruisseaux (MPO, IRD)
  * @since   0.5
- * @version 0.5
+ * @version 0.8
  * @module
  *
  * @see Plane
@@ -222,6 +223,8 @@ public class Line implements Cloneable,
      * least-squares senses. This method assume that the <var>x</var> values are precise and all uncertainty
      * is in <var>y</var>.
      *
+     * <p>The default implementation delegates to {@link #fit(Vector, Vector)}.</p>
+     *
      * @param  x  vector of <var>x</var> values (independent variable).
      * @param  y  vector of <var>y</var> values (dependent variable).
      * @return estimation of the correlation coefficient. The closer this coefficient is to +1 or -1, the better the fit.
@@ -229,8 +232,31 @@ public class Line implements Cloneable,
      * @throws IllegalArgumentException if <var>x</var> and <var>y</var> do not have the same length.
      */
     public double fit(final double[] x, final double[] y) {
-        // Do not invoke an overrideable method with our tricky iterable.
-        return doFit(new CompoundDirectPositions(x, y));
+        ArgumentChecks.ensureNonNull("x", x);
+        ArgumentChecks.ensureNonNull("y", y);
+        return fit(new ArrayVector.Doubles(x), new ArrayVector.Doubles(y));
+    }
+
+    /**
+     * Given a set of data points <var>x</var>[0 … <var>n</var>-1], <var>y</var>[0 … <var>n</var>-1],
+     * fits them to a straight line <var>y</var> = <var>slope</var>⋅<var>x</var> + <var>y₀</var> in a
+     * least-squares senses. This method assume that the <var>x</var> values are precise and all uncertainty
+     * is in <var>y</var>.
+     *
+     * <p>The default implementation delegates to {@link #fit(Iterable)}.</p>
+     *
+     * @param  x  vector of <var>x</var> values (independent variable).
+     * @param  y  vector of <var>y</var> values (dependent variable).
+     * @return estimation of the correlation coefficient. The closer this coefficient is to +1 or -1, the better the fit.
+     *
+     * @throws IllegalArgumentException if <var>x</var> and <var>y</var> do not have the same length.
+     *
+     * @since 0.8
+     */
+    public double fit(final Vector x, final Vector y) {
+        ArgumentChecks.ensureNonNull("x", x);
+        ArgumentChecks.ensureNonNull("y", y);
+        return fit(new CompoundDirectPositions(x, y));
     }
 
     /**
@@ -246,13 +272,7 @@ public class Line implements Cloneable,
      * @throws MismatchedDimensionException if a point is not two-dimensional.
      */
     public double fit(final Iterable<? extends DirectPosition> points) {
-        return doFit(points);
-    }
-
-    /**
-     * Implementation of public {@code fit(…)} methods.
-     */
-    private double doFit(final Iterable<? extends DirectPosition> points) {
+        ArgumentChecks.ensureNonNull("points", points);
         int i = 0, n = 0;
         final DoubleDouble mean_x = new DoubleDouble();
         final DoubleDouble mean_y = new DoubleDouble();

Modified: sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -16,12 +16,14 @@
  */
 package org.apache.sis.math;
 
+import java.util.Iterator;
 import java.io.Serializable;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.geometry.MismatchedDimensionException;
 import org.apache.sis.internal.util.DoubleDouble;
 import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.util.resources.Errors;
+import org.apache.sis.util.ArgumentChecks;
 
 import static java.lang.Math.abs;
 import static java.lang.Math.sqrt;
@@ -44,7 +46,7 @@ import static java.lang.Math.ulp;
  * @author  Martin Desruisseaux (MPO, IRD)
  * @author  Howard Freeland (MPO, for algorithmic inspiration)
  * @since   0.5
- * @version 0.5
+ * @version 0.8
  * @module
  *
  * @see Line
@@ -212,12 +214,11 @@ public class Plane implements Cloneable,
 
     /**
      * Computes the plane's coefficients from the given ordinate values.
-     * This method uses a linear regression in the least-square sense,
-     * with the assumption that the (<var>x</var>,<var>y</var>) values are precise
-     * and all uncertainty is in <var>z</var>.
+     * This method uses a linear regression in the least-square sense, with the assumption that
+     * the (<var>x</var>,<var>y</var>) values are precise and all uncertainty is in <var>z</var>.
+     * {@link Double#NaN} values are ignored. The result is undetermined if all points are colinear.
      *
-     * <p>{@link Double#NaN} values are ignored.
-     * The result is undetermined if all points are colinear.</p>
+     * <p>The default implementation delegates to {@link #fit(Vector, Vector, Vector)}.</p>
      *
      * @param  x  vector of <var>x</var> coordinates.
      * @param  y  vector of <var>y</var> coordinates.
@@ -226,154 +227,308 @@ public class Plane implements Cloneable,
      * @throws IllegalArgumentException if <var>x</var>, <var>y</var> and <var>z</var> do not have the same length.
      */
     public double fit(final double[] x, final double[] y, final double[] z) {
-        // Do not invoke an overrideable method with our tricky iterable.
-        return fit(new CompoundDirectPositions(x, y, z), true, true, true);
+        ArgumentChecks.ensureNonNull("x", x);
+        ArgumentChecks.ensureNonNull("y", y);
+        ArgumentChecks.ensureNonNull("z", z);
+        return fit(new ArrayVector.Doubles(x), new ArrayVector.Doubles(y), new ArrayVector.Doubles(z));
+    }
+
+    /**
+     * Computes the plane's coefficients from the given ordinate values.
+     * This method uses a linear regression in the least-square sense, with the assumption that
+     * the (<var>x</var>,<var>y</var>) values are precise and all uncertainty is in <var>z</var>.
+     * {@link Double#NaN} values are ignored. The result is undetermined if all points are colinear.
+     *
+     * <p>The default implementation delegates to {@link #fit(Iterable)}.</p>
+     *
+     * @param  x  vector of <var>x</var> coordinates.
+     * @param  y  vector of <var>y</var> coordinates.
+     * @param  z  vector of <var>z</var> values.
+     * @return an estimation of the Pearson correlation coefficient.
+     * @throws IllegalArgumentException if <var>x</var>, <var>y</var> and <var>z</var> do not have the same length.
+     *
+     * @since 0.8
+     */
+    public double fit(final Vector x, final Vector y, final Vector z) {
+        ArgumentChecks.ensureNonNull("x", x);
+        ArgumentChecks.ensureNonNull("y", y);
+        ArgumentChecks.ensureNonNull("z", z);
+        return fit(new CompoundDirectPositions(x, y, z));
+    }
+
+    /**
+     * Computes the plane's coefficients from values distributed on a regular grid. Invoking this method
+     * is equivalent to invoking {@link #fit(Vector, Vector, Vector)} where all vectors have a length of
+     * {@code nx} × {@code ny} and the <var>x</var> and <var>y</var> vectors have the following content:
+     *
+     * <table class="compact" summary="x and y vector content">
+     *   <tr>
+     *     <th><var>x</var> vector</th>
+     *     <th><var>y</var> vector</th>
+     *   </tr><tr>
+     *     <td>{@preformat math
+     *          0 1 2 3 4 5 … n<sub>x</sub>-1
+     *          0 1 2 3 4 5 … n<sub>x</sub>-1
+     *          0 1 2 3 4 5 … n<sub>x</sub>-1
+     *          …
+     *          0 1 2 3 4 5 … n<sub>x</sub>-1
+     *     }</td>
+     *     <td>{@preformat math
+     *          0 0 0 0 0 0 … 0
+     *          1 1 1 1 1 1 … 1
+     *          2 2 2 2 2 2 … 2
+     *          n<sub>y</sub>-1 n<sub>y</sub>-1 n<sub>y</sub>-1 … n<sub>y</sub>-1
+     *     }</td>
+     *   </tr>
+     * </table>
+     *
+     * This method uses a linear regression in the least-square sense, with the assumption that
+     * the (<var>x</var>,<var>y</var>) values are precise and all uncertainty is in <var>z</var>.
+     * {@link Double#NaN} values are ignored. The result is undetermined if all points are colinear.
+     *
+     * @param  nx  number of columns.
+     * @param  ny  number of rows.
+     * @param  z   values of a matrix of {@code nx} columns by {@code ny} rows organized in a row-major fashion.
+     * @return an estimation of the Pearson correlation coefficient.
+     * @throws IllegalArgumentException if <var>z</var> does not have the expected length.
+     *
+     * @since 0.8
+     */
+    public double fit(final int nx, final int ny, final Vector z) {
+        ArgumentChecks.ensureStrictlyPositive("nx", nx);
+        ArgumentChecks.ensureStrictlyPositive("ny", ny);
+        ArgumentChecks.ensureNonNull("z", z);
+        final int length = Math.multiplyExact(nx, ny);
+        if (z.size() != length) {
+            throw new IllegalArgumentException(Errors.format(Errors.Keys.UnexpectedArrayLength_2, length, z.size()));
+        }
+        final Fit r = new Fit(nx, ny, z);
+        r.resolve();
+        final double p = r.correlation(nx, length, z, null);
+        setEquation(r.sx.value, r.sy.value, r.z0.value);
+        return p;
     }
 
     /**
      * Computes the plane's coefficients from the given sequence of points.
-     * This method uses a linear regression in the least-square sense,
-     * with the assumption that the (<var>x</var>,<var>y</var>) values are precise
-     * and all uncertainty is in <var>z</var>.
-     *
-     * <p>Points shall be three dimensional with ordinate values in the (<var>x</var>,<var>y</var>,<var>z</var>) order.
-     * {@link Double#NaN} ordinate values are ignored.
-     * The result is undetermined if all points are colinear.</p>
+     * This method uses a linear regression in the least-square sense, with the assumption that
+     * the (<var>x</var>,<var>y</var>) values are precise and all uncertainty is in <var>z</var>.
+     * Points shall be three dimensional with ordinate values in the (<var>x</var>,<var>y</var>,<var>z</var>) order.
+     * {@link Double#NaN} values are ignored. The result is undetermined if all points are colinear.
      *
      * @param  points  the three-dimensional points.
      * @return an estimation of the Pearson correlation coefficient.
      * @throws MismatchedDimensionException if a point is not three-dimensional.
      */
     public double fit(final Iterable<? extends DirectPosition> points) {
-        return fit(points, true, true, true);
+        ArgumentChecks.ensureNonNull("points", points);
+        final Fit r = new Fit(points);
+        r.resolve();
+        final double p = r.correlation(0, 0, null, points.iterator());
+        /*
+         * Store the result only when we are done, so we have a "all or nothing" behavior.
+         * We invoke the setEquation(sx, sy, z₀) method in case the user override it.
+         */
+        setEquation(r.sx.value, r.sy.value, r.z0.value);
+        return p;
     }
 
     /**
-     * Implementation of public {@code fit(…)} methods.
-     * This method needs to iterate over the points two times:
+     * Computes the plane coefficients. This class needs to iterate over the points two times:
      * one for computing the coefficients, and one for the computing the Pearson coefficient.
      * The second pass can also opportunistically checks if some small coefficients can be replaced by zero.
      */
-    private double fit(final Iterable<? extends DirectPosition> points,
-            boolean detectZeroSx, boolean detectZeroSy, boolean detectZeroZ0)
-    {
-        int i = 0, n = 0;
-        final DoubleDouble sum_x  = new DoubleDouble();
-        final DoubleDouble sum_y  = new DoubleDouble();
-        final DoubleDouble sum_z  = new DoubleDouble();
-        final DoubleDouble sum_xx = new DoubleDouble();
-        final DoubleDouble sum_yy = new DoubleDouble();
-        final DoubleDouble sum_xy = new DoubleDouble();
-        final DoubleDouble sum_zx = new DoubleDouble();
-        final DoubleDouble sum_zy = new DoubleDouble();
-        final DoubleDouble xx     = new DoubleDouble();
-        final DoubleDouble yy     = new DoubleDouble();
-        final DoubleDouble xy     = new DoubleDouble();
-        final DoubleDouble zx     = new DoubleDouble();
-        final DoubleDouble zy     = new DoubleDouble();
-        for (final DirectPosition p : points) {
-            final int dimension = p.getDimension();
-            if (dimension != DIMENSION) {
-                throw new MismatchedDimensionException(Errors.format(
-                        Errors.Keys.MismatchedDimension_3, "points[" + i + ']', DIMENSION, dimension));
+    private static final class Fit {
+        private final DoubleDouble sum_x  = new DoubleDouble();
+        private final DoubleDouble sum_y  = new DoubleDouble();
+        private final DoubleDouble sum_z  = new DoubleDouble();
+        private final DoubleDouble sum_xx = new DoubleDouble();
+        private final DoubleDouble sum_yy = new DoubleDouble();
+        private final DoubleDouble sum_xy = new DoubleDouble();
+        private final DoubleDouble sum_zx = new DoubleDouble();
+        private final DoubleDouble sum_zy = new DoubleDouble();
+        private final DoubleDouble xx     = new DoubleDouble();
+        private final DoubleDouble yy     = new DoubleDouble();
+        private final DoubleDouble xy     = new DoubleDouble();
+        private final DoubleDouble zx     = new DoubleDouble();
+        private final DoubleDouble zy     = new DoubleDouble();
+        private final int n;
+
+        /** Solution of the plane equation. */ DoubleDouble sx, sy, z0;
+
+        /**
+         * Computes the values of all {@code sum_*} fields from randomly distributed points.
+         * Value of all other fields are undetermined..
+         */
+        Fit(final Iterable<? extends DirectPosition> points) {
+            int i = 0, n = 0;
+            for (final DirectPosition p : points) {
+                final int dimension = p.getDimension();
+                if (dimension != DIMENSION) {
+                    throw new MismatchedDimensionException(Errors.format(
+                            Errors.Keys.MismatchedDimension_3, "points[" + i + ']', DIMENSION, dimension));
+                }
+                i++;
+                final double x = p.getOrdinate(0); if (Double.isNaN(x)) continue;
+                final double y = p.getOrdinate(1); if (Double.isNaN(y)) continue;
+                final double z = p.getOrdinate(2); if (Double.isNaN(z)) continue;
+                xx.setToProduct(x, x);
+                yy.setToProduct(y, y);
+                xy.setToProduct(x, y);
+                zx.setToProduct(z, x);
+                zy.setToProduct(z, y);
+                sum_x .add(x );
+                sum_y .add(y );
+                sum_z .add(z );
+                sum_xx.add(xx);
+                sum_yy.add(yy);
+                sum_xy.add(xy);
+                sum_zx.add(zx);
+                sum_zy.add(zy);
+                n++;
             }
-            i++;
-            final double x = p.getOrdinate(0); if (Double.isNaN(x)) continue;
-            final double y = p.getOrdinate(1); if (Double.isNaN(y)) continue;
-            final double z = p.getOrdinate(2); if (Double.isNaN(z)) continue;
-            xx.setToProduct(x, x);
-            yy.setToProduct(y, y);
-            xy.setToProduct(x, y);
-            zx.setToProduct(z, x);
-            zy.setToProduct(z, y);
-            sum_x.add(x);
-            sum_y.add(y);
-            sum_z.add(z);
-            sum_xx.add(xx);
-            sum_yy.add(yy);
-            sum_xy.add(xy);
-            sum_zx.add(zx);
-            sum_zy.add(zy);
-            n++;
+            this.n = n;
         }
-        /*
-         *    ( sum_zx - sum_z*sum_x )  =  sx*(sum_xx - sum_x*sum_x) + sy*(sum_xy - sum_x*sum_y)
-         *    ( sum_zy - sum_z*sum_y )  =  sx*(sum_xy - sum_x*sum_y) + sy*(sum_yy - sum_y*sum_y)
-         */
-        zx.setFrom(sum_x); zx.divide(-n, 0); zx.multiply(sum_z); zx.add(sum_zx);    // zx = sum_zx - sum_z*sum_x/n
-        zy.setFrom(sum_y); zy.divide(-n, 0); zy.multiply(sum_z); zy.add(sum_zy);    // zy = sum_zy - sum_z*sum_y/n
-        xx.setFrom(sum_x); xx.divide(-n, 0); xx.multiply(sum_x); xx.add(sum_xx);    // xx = sum_xx - sum_x*sum_x/n
-        xy.setFrom(sum_y); xy.divide(-n, 0); xy.multiply(sum_x); xy.add(sum_xy);    // xy = sum_xy - sum_x*sum_y/n
-        yy.setFrom(sum_y); yy.divide(-n, 0); yy.multiply(sum_y); yy.add(sum_yy);    // yy = sum_yy - sum_y*sum_y/n
-        /*
-         * den = (xy*xy - xx*yy)
-         */
-        final DoubleDouble tmp = new DoubleDouble(xx); tmp.multiply(yy);
-        final DoubleDouble den = new DoubleDouble(xy); den.multiply(xy);
-        den.subtract(tmp);
-        /*
-         * sx = (zy*xy - zx*yy) / den
-         * sy = (zx*xy - zy*xx) / den
-         * z₀ = (sum_z - (sx*sum_x + sy*sum_y)) / n
-         */
-        final DoubleDouble sx = new DoubleDouble(zy); sx.multiply(xy); tmp.setFrom(zx); tmp.multiply(yy); sx.subtract(tmp); sx.divide(den);
-        final DoubleDouble sy = new DoubleDouble(zx); sy.multiply(xy); tmp.setFrom(zy); tmp.multiply(xx); sy.subtract(tmp); sy.divide(den);
-        final DoubleDouble z0 = new DoubleDouble(sy);
-        z0.multiply(sum_y);
-        tmp.setFrom(sx);
-        tmp.multiply(sum_x);
-        tmp.add(z0);
-        z0.setFrom(sum_z);
-        z0.subtract(tmp);
-        z0.divide(n, 0);
-        /*
-         * 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).
+
+        /**
+         * Computes the values of all {@code sum_*} fields from the <var>z</var> values on a regular grid.
+         * Value of all other fields are undetermined..
          */
-        final double mean_x = sum_x.value / n;
-        final double mean_y = sum_y.value / n;
-        final double mean_z = sum_z.value / n;
-        final double offset = abs((sx.value * mean_x + sy.value * mean_y) + z0.value); // Offsetted z₀ - see comment before usage.
-        double sum_ds2 = 0, sum_dz2 = 0, sum_dsz = 0;
-        for (final DirectPosition p : points) {
-            final double x = (p.getOrdinate(0) - mean_x) * sx.value;
-            final double y = (p.getOrdinate(1) - mean_y) * sy.value;
-            final double z = (p.getOrdinate(2) - mean_z);
-            final double s = x + y;
-            if (!Double.isNaN(s) && !Double.isNaN(z)) {
-                sum_ds2 += s * s;
-                sum_dz2 += z * z;
-                sum_dsz += s * z;
-            }
+        Fit(final int nx, final int ny, final Vector vz) {
             /*
-             * Algorithm for detecting if a coefficient should be zero:
-             * If for every points given by the user, adding (sx⋅x) in (sx⋅x + sy⋅y + z₀) does not make any difference
-             * because (sx⋅x) is smaller than 1 ULP of (sy⋅y + z₀), then it is not worth adding it and  sx  can be set
-             * to zero. The same rational applies to (sy⋅y) and z₀.
-             *
-             * Since we work with differences from the means, the  z = sx⋅x + sy⋅y + z₀  equation can be rewritten as:
+             * Computes the sum of x, y and z values. Computes also the sum of x², y², x⋅y, z⋅x and z⋅y values.
+             * When possible, we will avoid to compute the sum inside the loop and use the following identities
+             * instead:
              *
-             *     Δz = sx⋅Δx + sy⋅Δy + (sx⋅mx + sy⋅my + z₀ - mz)    where the term between (…) is close to zero.
+             *           1 + 2 + 3 ... + n    =    n⋅(n+1)/2              (arithmetic series)
+             *        1² + 2² + 3² ... + n²   =    n⋅(n+0.5)⋅(n+1)/3
              *
-             * The check for (sx⋅Δx) and (sy⋅Δy) below ignore the (…) term since it is close to zero.
-             * The check for  z₀  is derived from an equation without the  -mz  term.
+             * Note that for exclusive upper bound, we need to replace n by n-1 in above formulas.
              */
-            if (detectZeroSx && abs(x) >= ulp(y * ZERO_THRESHOLD)) detectZeroSx = false;
-            if (detectZeroSy && abs(y) >= ulp(x * ZERO_THRESHOLD)) detectZeroSy = false;
-            if (detectZeroZ0 && offset >= ulp(s * ZERO_THRESHOLD)) detectZeroZ0 = false;
+            int i = 0, n = 0;
+            for (int y=0; y<ny; y++) {
+                for (int x=0; x<nx; x++) {
+                    final double z = vz.doubleValue(i++);
+                    if (!Double.isNaN(z)) {
+                        zx.setToProduct(z, x);
+                        zy.setToProduct(z, y);
+                        sum_z .add(z );
+                        sum_zx.add(zx);
+                        sum_zy.add(zy);
+                        n++;
+                    }
+                }
+            }
+            sum_x .value = n/2d;  sum_x .multiply(nx-1,   0);                     // Division by 2 is exact.
+            sum_y .value = n/2d;  sum_y .multiply(ny-1,   0);
+            sum_xx.value = n;     sum_xx.multiply(nx-0.5, 0); sum_xx.multiply(nx-1, 0); sum_xx.divide(3, 0);
+            sum_yy.value = n;     sum_yy.multiply(ny-0.5, 0); sum_yy.multiply(ny-1, 0); sum_yy.divide(3, 0);
+            sum_xy.value = n/4d;  sum_xy.multiply(ny-1,   0); sum_xy.multiply(nx-1, 0);
+            this.n = n;
         }
-        /*
-         * Store the result only when we are done, so we have a "all or nothing" behavior.
-         * We invoke the setEquation(sx, sy, z₀) method in case the user override it.
+
+        /**
+         * Computes the {@link #sx}, {@link #sy} and {@link #z0} values using the sums computed by the constructor.
          */
-        setEquation(detectZeroSx ? 0 : sx.value,
-                    detectZeroSy ? 0 : sy.value,
-                    detectZeroZ0 ? 0 : z0.value);
-        return Math.min(sum_dsz / sqrt(sum_ds2 * sum_dz2), 1);
+        private void resolve() {
+            /*
+             *    ( sum_zx - sum_z⋅sum_x )  =  sx⋅(sum_xx - sum_x⋅sum_x) + sy⋅(sum_xy - sum_x⋅sum_y)
+             *    ( sum_zy - sum_z⋅sum_y )  =  sx⋅(sum_xy - sum_x⋅sum_y) + sy⋅(sum_yy - sum_y⋅sum_y)
+             */
+            zx.setFrom(sum_x); zx.divide(-n, 0); zx.multiply(sum_z); zx.add(sum_zx);    // zx = sum_zx - sum_z⋅sum_x/n
+            zy.setFrom(sum_y); zy.divide(-n, 0); zy.multiply(sum_z); zy.add(sum_zy);    // zy = sum_zy - sum_z⋅sum_y/n
+            xx.setFrom(sum_x); xx.divide(-n, 0); xx.multiply(sum_x); xx.add(sum_xx);    // xx = sum_xx - sum_x⋅sum_x/n
+            xy.setFrom(sum_y); xy.divide(-n, 0); xy.multiply(sum_x); xy.add(sum_xy);    // xy = sum_xy - sum_x⋅sum_y/n
+            yy.setFrom(sum_y); yy.divide(-n, 0); yy.multiply(sum_y); yy.add(sum_yy);    // yy = sum_yy - sum_y⋅sum_y/n
+            /*
+             * den = (xy⋅xy - xx⋅yy)
+             */
+            final DoubleDouble tmp = new DoubleDouble(xx); tmp.multiply(yy);
+            final DoubleDouble den = new DoubleDouble(xy); den.multiply(xy);
+            den.subtract(tmp);
+            /*
+             * sx = (zy⋅xy - zx⋅yy) / den
+             * sy = (zx⋅xy - zy⋅xx) / den
+             * z₀ = (sum_z - (sx⋅sum_x + sy⋅sum_y)) / n
+             */
+            sx = new DoubleDouble(zy); sx.multiply(xy); tmp.setFrom(zx); tmp.multiply(yy); sx.subtract(tmp); sx.divide(den);
+            sy = new DoubleDouble(zx); sy.multiply(xy); tmp.setFrom(zy); tmp.multiply(xx); sy.subtract(tmp); sy.divide(den);
+            z0 = new DoubleDouble(sy);
+            z0.multiply(sum_y);
+            tmp.setFrom(sx);
+            tmp.multiply(sum_x);
+            tmp.add(z0);
+            z0.setFrom(sum_z);
+            z0.subtract(tmp);
+            z0.divide(n, 0);
+        }
+
+        /**
+         * Computes an estimation of the Pearson correlation coefficient. We do not use double-double arithmetic
+         * here since the Pearson coefficient is for information purpose (quality estimation).
+         *
+         * <p>Only one of ({@code nx}, {@code length}, {@code z}) tuple or {@code points} argument should be non-null.</p>
+         */
+        double correlation(final int nx, final int length, final Vector vz,
+                           final Iterator<? extends DirectPosition> points)
+        {
+            boolean detectZeroSx = true;
+            boolean detectZeroSy = true;
+            boolean detectZeroZ0 = true;
+            final double sx     = this.sx.value;
+            final double sy     = this.sy.value;
+            final double z0     = this.z0.value;
+            final double mean_x = sum_x.value / n;
+            final double mean_y = sum_y.value / n;
+            final double mean_z = sum_z.value / n;
+            final double offset = abs((sx * mean_x + sy * mean_y) + z0);    // Offsetted z₀ - see comment before usage.
+            int index = 0;
+            double sum_ds2 = 0, sum_dz2 = 0, sum_dsz = 0;
+            for (;;) {
+                double x, y, z;
+                if (vz != null) {
+                    if (index >= length) break;
+                    x = index % nx;
+                    y = index / nx;
+                    z = vz.doubleValue(index++);
+                } else {
+                    if (!points.hasNext()) break;
+                    final DirectPosition p = points.next();
+                    x = p.getOrdinate(0);
+                    y = p.getOrdinate(1);
+                    z = p.getOrdinate(2);
+                }
+                x = (x - mean_x) * sx;
+                y = (y - mean_y) * sy;
+                z = (z - mean_z);
+                final double s = x + y;
+                if (!Double.isNaN(s) && !Double.isNaN(z)) {
+                    sum_ds2 += s * s;
+                    sum_dz2 += z * z;
+                    sum_dsz += s * z;
+                }
+                /*
+                 * Algorithm for detecting if a coefficient should be zero:
+                 * If for every points given by the user, adding (sx⋅x) in (sx⋅x + sy⋅y + z₀) does not make any difference
+                 * because (sx⋅x) is smaller than 1 ULP of (sy⋅y + z₀), then it is not worth adding it and  sx  can be set
+                 * to zero. The same rational applies to (sy⋅y) and z₀.
+                 *
+                 * Since we work with differences from the means, the  z = sx⋅x + sy⋅y + z₀  equation can be rewritten as:
+                 *
+                 *     Δz = sx⋅Δx + sy⋅Δy + (sx⋅mx + sy⋅my + z₀ - mz)    where the term between (…) is close to zero.
+                 *
+                 * The check for (sx⋅Δx) and (sy⋅Δy) below ignore the (…) term since it is close to zero.
+                 * The check for  z₀  is derived from an equation without the  -mz  term.
+                 */
+                if (detectZeroSx && abs(x) >= ulp(y * ZERO_THRESHOLD)) detectZeroSx = false;
+                if (detectZeroSy && abs(y) >= ulp(x * ZERO_THRESHOLD)) detectZeroSy = false;
+                if (detectZeroZ0 && offset >= ulp(s * ZERO_THRESHOLD)) detectZeroZ0 = false;
+            }
+            if (detectZeroSx) this.sx.clear();
+            if (detectZeroSy) this.sy.clear();
+            if (detectZeroZ0) this.z0.clear();
+            return Math.min(sum_dsz / sqrt(sum_ds2 * sum_dz2), 1);
+        }
     }
 
     /**

Modified: sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -70,6 +70,11 @@ public final class Errors extends Indexe
         public static final short AmbiguousName_3 = 1;
 
         /**
+         * This object can iterate only once.
+         */
+        public static final short CanIterateOnlyOnce = 172;
+
+        /**
          * No element can be added to this set because properties ‘{0}’ and ‘{1}’ are mutually
          * exclusive.
          */
@@ -640,6 +645,11 @@ public final class Errors extends Indexe
         public static final short NotAUnicodeIdentifier_1 = 112;
 
         /**
+         * {0} is not an integer value.
+         */
+        public static final short NotAnInteger_1 = 171;
+
+        /**
          * Argument ‘{0}’ shall not be null.
          */
         public static final short NullArgument_1 = 113;

Modified: sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties [ISO-8859-1] (original)
+++ sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties [ISO-8859-1] Thu Mar 16 18:32:29 2017
@@ -25,6 +25,7 @@
 # can reorder the parameters as they want.
 #
 AmbiguousName_3                   = Name \u201c{2}\u201d is ambiguous because it can be understood as either \u201c{0}\u201d or \u201c{1}\u201d.
+CanIterateOnlyOnce                = This object can iterate only once.
 CanNotAddToExclusiveSet_2         = No element can be added to this set because properties \u2018{0}\u2019 and \u2018{1}\u2019 are mutually exclusive.
 CanNotAssign_2                    = Can not assign \u201c{1}\u201d to \u201c{0}\u201d.
 CanNotAssignUnitToDimension_2     = Can not assign units \u201c{1}\u201d to dimension \u201c{0}\u201d.
@@ -133,6 +134,7 @@ NonTemporalUnit_1                 = \u20
 NonSystemUnit_1                   = \u201c{0}\u201d is not a fundamental or derived unit.
 NonRatioUnit_1                    = The scale of measurement for \u201c{0}\u201d unit is not a ratio scale.
 NotABackwardReference_1           = No element for the \u201c{0}\u201d identifier, or the identifier is a forward reference.
+NotAnInteger_1                    = {0} is not an integer value.
 NotAKeyValuePair_1                = \u201c{0}\u201d is not a key-value pair.
 NotANumber_1                      = Argument \u2018{0}\u2019 shall not be NaN (Not-a-Number).
 NotAPrimitiveWrapper_1            = Class \u2018{0}\u2019 is not a primitive type wrapper.

Modified: sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties [ISO-8859-1] (original)
+++ sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties [ISO-8859-1] Thu Mar 16 18:32:29 2017
@@ -22,6 +22,7 @@
 #   U+00A0 NO-BREAK SPACE         before  :
 #
 AmbiguousName_3                   = Le nom \u00ab\u202f{2}\u202f\u00bb est ambigu\u00eb car il peut \u00eatre interpr\u00e9t\u00e9 comme \u00ab\u202f{0}\u202f\u00bb ou \u00ab\u202f{1}\u202f\u00bb.
+CanIterateOnlyOnce                = Cet objet ne peut it\u00e9rer qu\u2019une seule fois.
 CanNotAddToExclusiveSet_2         = Aucun \u00e9l\u00e9ment ne peut \u00eatre ajout\u00e9 \u00e0 cet ensemble car les propri\u00e9t\u00e9s \u2018{0}\u2019 et \u2018{1}\u2019 sont mutuellement exclusives.
 CanNotAssign_2                    = Ne peut pas assigner \u00ab\u202f{1}\u202f\u00bb \u00e0 \u00ab\u202f{0}\u202f\u00bb.
 CanNotAssignUnitToDimension_2     = Ne peut pas assigner les unit\u00e9s \u00ab\u202f{1}\u202f\u00bb \u00e0 la dimension \u00ab\u202f{0}\u202f\u00bb.
@@ -130,6 +131,7 @@ NonTemporalUnit_1                 = \u00
 NonSystemUnit_1                   = \u00ab\u202f{0}\u202f\u00bb n\u2019est pas une unit\u00e9 fondamentale ou d\u00e9riv\u00e9e.
 NonRatioUnit_1                    = L\u2019\u00e9chelle de mesure de l\u2019unit\u00e9 \u00ab\u202f{0}\u202f\u00bb n\u2019est pas une \u00e9chelle de rapports.
 NotABackwardReference_1           = Il n\u2019y a pas d\u2019\u00e9l\u00e9ment pour l\u2019identifiant \u201c{0}\u201d, ou l\u2019identifiant est une r\u00e9f\u00e9rence vers l\u2019avant.
+NotAnInteger_1                    = {0} n\u2019est pas un nombre entier.
 NotAKeyValuePair_1                = \u00ab\u202f{0}\u202f\u00bb n\u2019est pas une paire cl\u00e9-valeur.
 NotANumber_1                      = L\u2019argument \u2018{0}\u2019 ne doit pas \u00eatre NaN (Not-a-Number).
 NotAPrimitiveWrapper_1            = La classe \u2018{0}\u2019 n\u2019est pas un adaptateur d\u2019un type primitif.

Modified: sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/math/PlaneTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/math/PlaneTest.java?rev=1787227&r1=1787226&r2=1787227&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/math/PlaneTest.java [UTF-8] (original)
+++ sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/math/PlaneTest.java [UTF-8] Thu Mar 16 18:32:29 2017
@@ -31,12 +31,18 @@ import static org.apache.sis.test.Assert
  *
  * @author  Martin Desruisseaux (MPO, IRD)
  * @since   0.5
- * @version 0.5
+ * @version 0.8
  * @module
  */
 @DependsOn(org.apache.sis.internal.util.DoubleDoubleTest.class)
 public final strictfp class PlaneTest extends TestCase {
     /**
+     * The Pearson coefficient computed by the last call to
+     * {@link #assertFitEquals(double, double[], double[], double[])}.
+     */
+    private double pearson;
+
+    /**
      * Invokes {@link Plane#fit(DirectPosition[])} with the given arrays,
      * and compares the fitted values against the original values.
      *
@@ -44,9 +50,9 @@ public final strictfp class PlaneTest ex
      *
      * @param  tolerance  the maximal difference allowed between the fitted and the original values.
      */
-    private static Plane assertFitEquals(final double tolerance, final double[] x, final double[] y, final double[] z) {
+    private Plane assertFitEquals(final double tolerance, final double[] x, final double[] y, final double[] z) {
         final Plane plan = new Plane();
-        final double pearson = plan.fit(x, y, z);
+        pearson = plan.fit(x, y, z);
         assertEquals("Pearson correlation coefficient", 1, pearson, 0.01);
         for (int i = 0; i < z.length; i++) {
             assertEquals("z(x,y)", z[i], plan.z(x[i], y[i]), tolerance);
@@ -59,11 +65,11 @@ public final strictfp class PlaneTest ex
      *
      * @param  rd      the random number generator to use.
      * @param  length  the desired array length.
-     * @param  offset  the minimal value allowed in the returned array.
      * @param  scale   the difference between the minimal and maximal allowed values.
+     * @param  offset  the minimal value allowed in the returned array.
      * @return the array of random values.
      */
-    private static double[] random(final Random rd, final int length, final double offset, final double scale) {
+    private static double[] random(final Random rd, final int length, final double scale, final double offset) {
         final double[] x = new double[length];
         for (int i=0; i<length; i++) {
             x[i] = offset + scale * rd.nextDouble();
@@ -77,10 +83,9 @@ public final strictfp class PlaneTest ex
      */
     @Test
     public void testFit3Points() {
-        // We fix the random seed in order to avoid colinear points.
-        final Random rd = new Random(4762038814364076443L);
+        final Random rd = new Random(4762038814364076443L);     // Fix the random seed for avoiding colinear points.
         for (int n=0; n<10; n++) {
-            assertFitEquals(1E-9, // We expect close to exact matches.
+            assertFitEquals(1E-9,                               // We expect close to exact matches.
                     random(rd, 3, 100, -25),
                     random(rd, 3,  80, -20),
                     random(rd, 3, 150, -40));
@@ -96,16 +101,67 @@ public final strictfp class PlaneTest ex
         final Random  rd = new Random(7241997054993322309L);
         final double x[] = random(rd, 4000, 100, -25);
         final double y[] = random(rd, 4000,  80, -20);
-        final double z[] = random(rd, 4000, 150, -40);
-        final Plane plan = new Plane(2, 3, -4);
+        final double z[] = random(rd, 4000,  10,  -5);
+        final Plane reference = new Plane(2, 3, -4);
         for (int i=0; i<z.length; i++) {
-            // Compute points with random displacements above or below a known plane.
-            z[i] = plan.z(x[i], y[i]) + (10 * rd.nextDouble() - 5);
+            z[i] += reference.z(x[i], y[i]);    // Compute points with random displacements above or below a known plane.
         }
+        /*
+         * For a random seed fixed to 7241997054993322309, the coefficients are:
+         *
+         *   - reference:     z(x,y) = 2.0⋅x + 3.0⋅y + -4.0
+         *   - fitted:        z(x,y) = 2.001595888896693⋅x + 3.0021028196088055⋅y + -4.105960575835259
+         */
         final Plane fitted = assertFitEquals(6, x, y, z);
-        assertEquals("sx", plan.slopeX(), fitted.slopeX(), 0.01);
-        assertEquals("sy", plan.slopeY(), fitted.slopeY(), 0.01);
-        assertEquals("z₀", plan.z0(),     fitted.z0(),     1);
+        assertEquals("sx", reference.slopeX(), fitted.slopeX(), 0.002);
+        assertEquals("sy", reference.slopeY(), fitted.slopeY(), 0.003);
+        assertEquals("z₀", reference.z0(),     fitted.z0(),     0.2);
+    }
+
+    /**
+     * Verifies that {@link Plane#fit(int, int, Vector)} produces the same result than
+     * {@link Plane#fit(double[], double[], double[])}.
+     *
+     * @since 0.8
+     */
+    @Test
+    @DependsOnMethod("testFitManyPoints")
+    public void testFitGrid() {
+        final int nx = 20;
+        final int ny = 30;
+        final Random  rd = new Random(1224444984079270867L);
+        final double z[] = random(rd, nx*ny, 12, -6);
+        final double x[] = new double[z.length];
+        final double y[] = new double[z.length];
+        final Plane reference = new Plane(-46, 17, 35);
+        for (int i=0; i<z.length; i++) {
+            x[i] = i % nx;
+            y[i] = i / nx;
+            z[i] += reference.z(x[i], y[i]);
+        }
+        /*
+         * Opportunistically verify the result of fit(double[], double[], double[]), but it is
+         * not the purpose of this test (it should have been verified by testFitManyPoints()).
+         * For a random seed fixed to 7241997054993322309, the coefficients are:
+         *
+         *   - reference:     z(x,y) = -46.0⋅x + 17.0⋅y + 35.0
+         *   - fitted:        z(x,y) = -45.96034442769177⋅x + 16.981792790278828⋅y + 34.901440258861314
+         */
+        final Plane fitted = assertFitEquals(7, x, y, z);
+        assertEquals("sx", reference.slopeX(), fitted.slopeX(), 0.05);
+        assertEquals("sy", reference.slopeY(), fitted.slopeY(), 0.02);
+        assertEquals("z₀", reference.z0(),     fitted.z0(),     0.1);
+        final double ep = pearson;
+        /*
+         * Verify that the optimized code path for a grid produces the exact same result than the
+         * generic code path.
+         */
+        final Plane gf = new Plane();
+        gf.fit(nx, ny, Vector.create(z, false));
+        assertEquals("sx", fitted.slopeX(), gf.slopeX(), STRICT);
+        assertEquals("sy", fitted.slopeY(), gf.slopeY(), STRICT);
+        assertEquals("z₀", fitted.z0(),     gf.z0(),     STRICT);
+        assertEquals("Pearson", ep,         pearson,     STRICT);
     }
 
     /**



Mime
View raw message