sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: First attempt to apply the fix for NADCON grids crossing anti-meridian to NetCDF localization grids crossing the anti-meridian. Current fix may be okay for grid without rotations or with a 90° rotation, but may not be okay in the general case. A better fix may require an API change in DatumShiftGrid.replaceOutsideGridCoordinate(…).
Date Wed, 19 Feb 2020 19:01:52 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new f8b3f5e  First attempt to apply the fix for NADCON grids crossing anti-meridian to
NetCDF localization grids crossing the anti-meridian. Current fix may be okay for grid without
rotations or with a 90° rotation, but may not be okay in the general case. A better fix may
require an API change in DatumShiftGrid.replaceOutsideGridCoordinate(…).
f8b3f5e is described below

commit f8b3f5ea6ea815702c9ba59a8e3f1e80c6845d76
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Feb 19 19:58:35 2020 +0100

    First attempt to apply the fix for NADCON grids crossing anti-meridian to NetCDF localization
grids crossing the anti-meridian.
    Current fix may be okay for grid without rotations or with a 90° rotation, but may not
be okay in the general case.
    A better fix may require an API change in DatumShiftGrid.replaceOutsideGridCoordinate(…).
---
 .../referencing/provider/DatumShiftGridFile.java   | 16 ++---
 .../sis/referencing/datum/DatumShiftGrid.java      |  6 +-
 .../operation/builder/LinearTransformBuilder.java  |  5 ++
 .../operation/builder/LocalizationGridBuilder.java | 40 +++++++++----
 .../operation/builder/ResidualGrid.java            | 68 +++++++++++++++++++++-
 .../operation/builder/package-info.java            |  2 +-
 .../sis/referencing/operation/builder/readme.html  |  6 +-
 .../operation/transform/LinearTransform.java       |  3 +-
 .../operation/transform/ProjectiveTransform.java   | 44 +++++++++-----
 .../operation/builder/ResidualGridTest.java        |  8 ++-
 10 files changed, 155 insertions(+), 43 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
index 7ebb2d7..6de48b8 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
@@ -124,7 +124,7 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
      *
      * @see #replaceOutsideGridCoordinate(int, double)
      */
-    private final double cycle;
+    private final double periodX;
 
     /**
      * The best translation accuracy that we can expect from this file.
@@ -186,9 +186,9 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
         this.files      = files;
         this.nx         = nx;
         if (Units.isAngular(coordinateUnit)) {
-            cycle = Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE) / Math.abs(Δx));
+            periodX = Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE) / Math.abs(Δx));
         } else {
-            cycle = 0;
+            periodX = 0;
             /*
              * Note: non-angular source coordinates are currently never used in this package.
              * If it continue to be like that in the future, we should remove the check for
@@ -210,7 +210,7 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
         nx         = other.nx;
         accuracy   = other.accuracy;
         subgrids   = other.subgrids;
-        cycle      = other.cycle;
+        periodX    = other.periodX;
     }
 
     /**
@@ -233,8 +233,8 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
         descriptor = other.descriptor;
         files      = other.files;
         this.nx    = nx;
-        cycle      = (other.cycle == 0) ? 0 :
-                Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE) / AffineTransforms2D.getScaleX0(gridToCRS));
+        periodX    = (other.periodX == 0) ? 0 : Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE)
+                                              / AffineTransforms2D.getScaleX0(gridToCRS));
         // Accuracy to be set by caller. Initial value needs to be zero.
     }
 
@@ -396,8 +396,8 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
      */
     @Override
     protected double replaceOutsideGridCoordinate(final int dimension, final double gridCoordinate)
{
-        if (dimension == 0 && cycle != 0) {
-            return Math.IEEEremainder(gridCoordinate, cycle);
+        if (dimension == 0 && periodX != 0) {
+            return Math.IEEEremainder(gridCoordinate, periodX);
         }
         return super.replaceOutsideGridCoordinate(dimension, gridCoordinate);
     }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
index 158540a..87f074f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
@@ -767,8 +767,8 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T
extends Quantity<T
      *
      * <p>The default implementation returns {@code gridCoordinate} unchanged. Subclasses
need to override this
      * method if they want to handle longitude cycles. Note that the coordinate value is
a grid index, not a
-     * longitude value. So the cycle to add or remove is the number of cells that the grid
would have if it was
-     * spanning 360° of longitude.</p>
+     * longitude value. So the period to add or remove is the number of cells that the grid
would have if it
+     * was spanning 360° of longitude.</p>
      *
      * <div class="note"><b>Example:</b>
      * this method may be implemented as below:
@@ -777,7 +777,7 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T
extends Quantity<T
      *     &#64;Override
      *     protected double replaceOutsideGridCoordinate(int dimension, double gridCoordinate)
{
      *         if (dimension == 0) {
-     *             return Math.IEEEremainder(gridCoordinate, cycle);
+     *             return Math.IEEEremainder(gridCoordinate, periodX);
      *         } else {
      *             return super.replaceOutsideGridCoordinate(dimension, gridCoordinate);
      *         }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
index 412f7d7..d5849ab 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
@@ -1324,6 +1324,11 @@ search:         for (int j=domain(); --j >= 0;) {
      * Computes the matrix of the linear approximation. This is the implementation of {@link
#create(MathTransformFactory)}
      * without the step creating the {@link LinearTransform} from a matrix. The {@link #correlations}
field is set as a side
      * effect of this method call.
+     *
+     * <p>In current implementation, the transform represented by the returned matrix
is always affine
+     * (i.e. the last row is fixed to [0 0 … 0 1]). If this is no longer the case in a
future version,
+     * some codes may not work anymore. Search for {@code isAffine()} statements for locating
codes
+     * that depend on affine transform assumption.</p>
      */
     @SuppressWarnings("serial")
     private MatrixSIS fit() throws FactoryException {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
index 8952f30..c8046ab 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
@@ -84,7 +84,7 @@ import static org.apache.sis.referencing.operation.builder.ResidualGrid.SOURCE_D
  * See the <cite>Linearizers</cite> section in {@link LinearTransformBuilder}
for more discussion.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @see InterpolatedTransform
  * @see LinearTransform
@@ -139,6 +139,17 @@ public class LocalizationGridBuilder extends TransformBuilder {
     private MathTransform transform;
 
     /**
+     * If the coordinates in some dimensions are cyclic, their periods. Otherwise {@code
null}.
+     * Values are in units of the target CRS. For longitude wraparounds, the period is typically
360°.
+     * Array length shall be {@code linear.getTargetDimensions()} and non-cyclic dimensions
shall have
+     * a period of zero (not {@link Double#NaN}, because we will use this array as a displacement
vector).
+     *
+     * @see #resolveWraparoundAxis(int, int, double)
+     * @see ResidualGrid#cellPeriods
+     */
+    private double[] periods;
+
+    /**
      * Creates a new, initially empty, builder for a localization grid of the given size.
      *
      * @param width   the number of columns in the grid of target positions.
@@ -549,7 +560,8 @@ public class LocalizationGridBuilder extends TransformBuilder {
      *                    The recommended direction is the direction of most stable values,
typically 1 (rows) for longitudes.
      * @param  period     that wraparound range (typically 360° for longitudes). Must be
strictly positive.
      * @return the range of coordinate values in the specified dimension after correction
for wraparound values.
-     * @throws IllegalStateException if {@link #create(MathTransformFactory) create(…)}
has already been invoked.
+     * @throws IllegalStateException if this method has already been invoked for the same
dimension,
+     *         or if {@link #create(MathTransformFactory) create(…)} has already been invoked.
      *
      * @since 1.0
      */
@@ -558,6 +570,14 @@ public class LocalizationGridBuilder extends TransformBuilder {
         ArgumentChecks.ensureBetween("dimension", 0, linear.getTargetDimensions() - 1, dimension);
         ArgumentChecks.ensureBetween("direction", 0, linear.getSourceDimensions() - 1, direction);
         ArgumentChecks.ensureStrictlyPositive("period", period);
+        if (periods == null) {
+            periods = new double[linear.getTargetDimensions()];
+        }
+        if (periods[dimension] != 0) {
+            throw new IllegalStateException(Errors.format(
+                    Errors.Keys.ValueAlreadyDefined_1, Strings.bracket("periods", dimension)));
+        }
+        periods[dimension] = period;
         return linear.resolveWraparoundAxis(dimension, direction, period);
     }
 
@@ -670,16 +690,16 @@ public class LocalizationGridBuilder extends TransformBuilder {
                             residual[k++] = (float) dy;
                         }
                     }
+                    if (isLinear) {
+                        step = MathTransforms.concatenate(sourceToGrid, gridToCoord);
+                    } else {
+                        step = InterpolatedTransform.createGeodeticTransformation(nonNull(factory),
+                                new ResidualGrid(sourceToGrid, gridToCoord, width, height,
residual,
+                                (gridPrecision > 0) ? gridPrecision : DEFAULT_PRECISION,
periods));
+                    }
                 } catch (TransformException e) {
                     throw new FactoryException(e);                                      
   // Should never happen.
                 }
-                if (isLinear) {
-                    step = MathTransforms.concatenate(sourceToGrid, gridToCoord);
-                } else {
-                    step = InterpolatedTransform.createGeodeticTransformation(nonNull(factory),
-                            new ResidualGrid(sourceToGrid, gridToCoord, width, height, residual,
-                            (gridPrecision > 0) ? gridPrecision : DEFAULT_PRECISION));
-                }
             }
             /*
              * At this point we finished to compute the transformation to target coordinates.
@@ -693,7 +713,7 @@ public class LocalizationGridBuilder extends TransformBuilder {
                 throw new InvalidGeodeticParameterException(Resources.format(
                         Resources.Keys.NonInvertibleOperation_1, linear.linearizerID()),
e);
             }
-            transform = step;
+            transform = step;                               // Set only after everything
succeeded.
         }
         return transform;
     }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ResidualGrid.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ResidualGrid.java
index b8e4304..85ffe2a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ResidualGrid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ResidualGrid.java
@@ -22,6 +22,7 @@ import javax.measure.quantity.Dimensionless;
 import org.opengis.parameter.ParameterDescriptor;
 import org.opengis.parameter.ParameterDescriptorGroup;
 import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.parameter.ParameterBuilder;
 import org.apache.sis.referencing.datum.DatumShiftGrid;
@@ -42,7 +43,7 @@ import org.apache.sis.measure.Units;
  * The residuals after an affine approximation has been created for a set of matching control
point pairs.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.8
  * @module
  */
@@ -133,21 +134,68 @@ final class ResidualGrid extends DatumShiftGrid<Dimensionless,Dimensionless>
{
     private final double accuracy;
 
     /**
+     * If grid coordinates in some target dimensions are cyclic, their periods in number
of cells.
+     * Otherwise {@code null}. If non-null, non-cyclic dimensions shall have a period of
zero.
+     *
+     * <p>Note that the unit of measurement is different than {@link LocalizationGridBuilder#periods}.
+     * This fields contain values in number of cells, which must have been converted from
values in
+     * degrees (for example). Those numbers may be non-integers.</p>
+     *
+     * @see LocalizationGridBuilder#periods
+     */
+    private final double[] cellPeriods;
+
+    /**
      * Creates a new residual grid.
      *
      * @param sourceToGrid  conversion from the "real world" source coordinates to grid indices
including fractional parts.
      * @param gridToTarget  conversion from grid coordinates to the final "real world" coordinates.
      * @param residuals     the residual data, as translations to apply on the result of
affine transform.
      * @param precision     desired precision of inverse transformations in unit of grid
cells.
+     * @param periods       if grid coordinates in some dimensions are cyclic, their periods
in units of target CRS.
      */
     ResidualGrid(final LinearTransform sourceToGrid, final LinearTransform gridToTarget,
-            final int nx, final int ny, final float[] residuals, final double precision)
+            final int nx, final int ny, final float[] residuals, final double precision,
+            final double[] periods) throws TransformException
     {
         super(Units.UNITY, sourceToGrid, new int[] {nx, ny}, true, Units.UNITY);
         this.gridToTarget = gridToTarget;
         this.offsets      = residuals;
         this.accuracy     = precision;
         this.nx           = nx;
+        double[] cellPeriods = null;
+        if (periods != null && gridToTarget.isAffine()) {
+            /*
+             * We require the transform to be affine because it makes the Jacobian independent
of
+             * coordinate values. It allows us to replace a period in target units by periods
in
+             * grid units without having to take the coordinate values in account.
+             */
+            final MatrixSIS m = MatrixSIS.castOrCopy(gridToTarget.inverse().derivative(null));
+            cellPeriods = m.multiply(periods);
+            /*
+             * Ideally we would be done. But if we have more than one wraparound dimension
+             * (i.e. if the `periods` array contains more than 1 non-zero value) and if the
+             * `gridToTarget` transform combines the coordinates in many dimensions, we may
+             * have the contributions of many dimensions mixed in a way that we will not
be
+             * able to use. This loop verify that each `cellPeriods` value depends on only
+             * one wraparound dimension of the target CRS. If this is not the case, we have
+             * ambiguity; value will be reset to 0 for disabling wraparound on that dimension.
+             */
+            for (int j=0; j<cellPeriods.length; j++) {
+                int contributions = 0;
+                final double cp = Math.abs(cellPeriods[j]);
+                if (cp > 0 && cp >= getGridSize()[j]) {
+                    final double tolerance = Math.ulp(cp);
+                    for (int i=0; i<periods.length; i++) {
+                        if (!(Math.abs(periods[i] * m.getElement(j,i)) <= tolerance))
{     // ! for counting NaNs.
+                            if (++contributions >= 2) break;                         
      // No need to continue.
+                        }
+                    }
+                }
+                cellPeriods[j] = (contributions == 1) ? cp : 0;
+            }
+        }
+        this.cellPeriods = cellPeriods;
     }
 
     /**
@@ -196,6 +244,22 @@ final class ResidualGrid extends DatumShiftGrid<Dimensionless,Dimensionless>
{
     }
 
     /**
+     * Invoked when a {@code gridX} or {@code gridY} coordinate is outside the range of valid
grid coordinates.
+     * If the coordinate outside the range is a longitude value and if we handle those values
as cyclic, brings
+     * that coordinate inside the range.
+     */
+    @Override
+    protected double replaceOutsideGridCoordinate(final int dimension, final double gridCoordinate)
{
+        if (cellPeriods != null) {
+            final double p = cellPeriods[dimension];
+            if (p != 0) {
+                return Math.IEEEremainder(gridCoordinate, p);
+            }
+        }
+        return super.replaceOutsideGridCoordinate(dimension, gridCoordinate);
+    }
+
+    /**
      * View over one target dimension of the localization grid. Used for populating the {@link
ParameterDescriptorGroup}
      * that describes the {@code MathTransform}. Those parameters are themselves used for
formatting Well Known Text.
      * Current implementation can be used only when the number of grid dimensions is {@value
#INTERPOLATED_DIMENSIONS}.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
index 5c0774b..2ffb543 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
@@ -30,7 +30,7 @@
  * convenience.</p>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.5
  * @module
  */
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/readme.html
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/readme.html
index bef92a0..6cd3ec1 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/readme.html
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/readme.html
@@ -21,15 +21,17 @@
       typical linear regression method assumes that <var>x</var> values are exact
and all errors are in <var>y</var> values.
       Applied to the <code>LocalizationGridBuilder</code> context, it means that
linear regression method
       assumes that grid indices are exact and all errors are in geospatial coordinates.
-      This is a reasonable assumption if the linear regression is used directly.
+    </p><p>
+      The assumption that all errors are in geospatial coordinates is reasonable if the linear
regression is used directly.
       But in <code>LocalizationGridBuilder</code> context, having the smallest
errors on geospatial coordinates
       may not be so important because those errors are corrected by the residual grids during
<em>forward</em> transformations.
       However during <em>inverse</em> transformations, it may be useful that
grid indices estimated by the linear regression
       are as close as possible to the real grid indices in order to allow iterations to converge
faster
       (such iterations exist only in inverse operations, not in forward operations).
       For that reason, <code>LocalizationGridBuilder</code> may want to minimize
errors on grid indices instead than geospatial coordinates.
-      We can achieve this result by computing linear regression coefficients for an equation
of the form
+      We could achieve this result by computing linear regression coefficients for an equation
of the form
       <var>x</var> = <var>C₀</var> + <var>C₁</var> × <var>y</var>,
then inverting that equation.
+    </p><p>
       This approach was attempted in a previous version and gave encouraging results, but
we nevertheless reverted that change.
       One reason is that even if inverting the linear regression calculation allowed iterations
to converge more often with curved
       localization grids, it was not sufficient anyway: we still had too many cases where
inverse transformations did not converge.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform.java
index 9cac821..3776331 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform.java
@@ -75,7 +75,8 @@ import org.opengis.referencing.operation.NoninvertibleTransformException;
 public interface LinearTransform extends MathTransform {
     /**
      * Returns {@code true} if this transform is affine.
-     * An affine transform preserves parallelism.
+     * An affine transform preserves parallelism and has the same
+     * {@linkplain AbstractMathTransform#derivative derivative} at every locations.
      *
      * @return {@code true} if this transform is affine.
      *
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
index a4efca3..7516f1f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
@@ -34,7 +34,7 @@ import org.apache.sis.util.ArgumentChecks;
  * lines in the source is preserved in the output.</p>
  *
  * @author  Martin Desruisseaux (IRD, Geomatys)
- * @version 0.7
+ * @version 1.1
  *
  * @see java.awt.geom.AffineTransform
  *
@@ -263,7 +263,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
         while (--numPts >= 0) {
             int mix = 0;
             for (int j=0; j<numRow; j++) {
-                double sum = elt[mix + srcDim];
+                double sum = elt[mix + srcDim];         // Initialize to translation term.
                 for (int i=0; i<srcDim; i++) {
                     final double e = elt[mix++];
                     if (e != 0) {
@@ -282,7 +282,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
             }
             final double w = buffer[dstDim];
             for (int j=0; j<dstDim; j++) {
-                // 'w' is equal to 1 if the transform is affine.
+                // `w` is equal to 1 if the transform is affine.
                 dstPts[dstOff + j] = buffer[j] / w;
             }
             srcOff += srcInc;
@@ -309,8 +309,8 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
     @Override
     public final void transform(float[] srcPts, int srcOff, final float[] dstPts, int dstOff,
int numPts) {
         final int srcDim, dstDim;
-        int srcInc = srcDim = numCol-1;
-        int dstInc = dstDim = numRow-1;
+        int srcInc = srcDim = numCol - 1;
+        int dstInc = dstDim = numRow - 1;
         if (srcPts == dstPts) {
             switch (IterationStrategy.suggest(srcOff, srcDim, dstOff, dstDim, numPts)) {
                 case ASCENDING: {
@@ -337,7 +337,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
                 double sum = elt[mix + srcDim];
                 for (int i=0; i<srcDim; i++) {
                     final double e = elt[mix++];
-                    if (e != 0) { // See comment in transform(double[], ...)
+                    if (e != 0) {                                   // See comment in transform(double[],
...)
                         sum += srcPts[srcOff + i] * e;
                     }
                 }
@@ -364,8 +364,8 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
      */
     @Override
     public final void transform(final double[] srcPts, int srcOff, final float[] dstPts,
int dstOff, int numPts) {
-        final int srcDim = numCol-1;
-        final int dstDim = numRow-1;
+        final int srcDim = numCol - 1;
+        final int dstDim = numRow - 1;
         final double[] buffer = new double[numRow];
         while (--numPts >= 0) {
             int mix = 0;
@@ -424,20 +424,38 @@ class ProjectiveTransform extends AbstractLinearTransform implements
ExtendedPre
     }
 
     /**
-     * Gets the derivative of this transform at a point.
-     * For a matrix transform, the derivative is the same everywhere.
+     * Gets the derivative of this transform at a point. For an affine transform,
+     * the derivative is the same everywhere and the given point is ignored.
+     * For a perspective transform, the given point is used.
      *
-     * @param  point  ignored (can be {@code null}).
+     * @param  point  the coordinate point where to evaluate the derivative.
      */
     @Override
     public final Matrix derivative(final DirectPosition point) {
         final int srcDim = numCol - 1;
         final int dstDim = numRow - 1;
+        /*
+         * In the `transform(…)` method, all coordinate values are divided by a `w` coefficient
+         * which depends on the position. We need to reproduce that division here. Note that
`w`
+         * coefficient is different than 1 only if the transform is non-affine.
+         */
+        int mix = dstDim * numCol;
+        double w = elt[mix + srcDim];                   // `w` is equal to 1 if the transform
is affine.
+        for (int i=0; i<srcDim; i++) {
+            final double e = elt[mix++];
+            if (e != 0) {                               // For avoiding NullPointerException
if affine.
+                w += point.getOrdinate(i) * e;
+            }
+        }
+        /*
+         * In the usual affine case (w=1), the derivative is a copy of this matrix
+         * with last row and last column omitted.
+         */
+        mix = 0;
         final MatrixSIS matrix = Matrices.createZero(dstDim, srcDim);
-        int mix = 0;
         for (int j=0; j<dstDim; j++) {
             for (int i=0; i<srcDim; i++) {
-                matrix.setElement(j, i, elt[mix++]);
+                matrix.setElement(j, i, elt[mix++] / w);
             }
             mix++;                                      // Skip translation column.
         }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/ResidualGridTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/ResidualGridTest.java
index 189348d..5cda046 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/ResidualGridTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/ResidualGridTest.java
@@ -31,7 +31,7 @@ import static org.opengis.test.Assert.*;
  * Tests {@link ResidualGrid}.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   1.0
  * @module
  */
@@ -45,13 +45,15 @@ public final strictfp class ResidualGridTest extends TestCase {
      * Creates a new test case with a 3×4 grid with 2 values in each cells.
      * Those two values are typically the horizontal components of translation vectors.
      * The grid has no "source to grid" or "grid to CRS" transformations.
+     *
+     * @throws TransformException if an error occurred while handling a wraparound axis.
      */
-    public ResidualGridTest() {
+    public ResidualGridTest() throws TransformException {
         grid = new ResidualGrid(MathTransforms.identity(2), MathTransforms.identity(2), 3,
4, new float[] {
                 0,2  ,  1,2  ,  2,1,
                 1,3  ,  2,2  ,  1,1,
                 0,4  ,  2,3  ,  3,2,
-                1,4  ,  3,3  ,  3,2}, 0.1);
+                1,4  ,  3,3  ,  3,2}, 0.1, null);
     }
 
     /**


Mime
View raw message