/*
* GeoTools - The Open Source Java GIS Toolkit
* http://geotools.org
*
* (C) 2002-2008, Open Source Geospatial Foundation (OSGeo)
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation;
* version 2.1 of the License.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*/
package org.geotools.referencing.operation.builder;
import java.awt.Point;
import java.awt.geom.AffineTransform;
import java.awt.geom.NoninvertibleTransformException;
import java.awt.geom.Point2D;
import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.Serializable;
import java.util.Arrays;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.parameter.ParameterDescriptor;
import org.opengis.parameter.ParameterDescriptorGroup;
import org.opengis.parameter.ParameterNotFoundException;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.MathTransform2D;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.Transformation;
import org.opengis.referencing.operation.TransformException;
import org.geotools.metadata.iso.citation.Citations;
import org.geotools.referencing.NamedIdentifier;
import org.geotools.referencing.operation.MathTransformProvider;
import org.geotools.referencing.operation.matrix.Matrix2;
import org.geotools.referencing.operation.transform.AbstractMathTransform;
import org.geotools.util.logging.Logging;
import org.geotools.util.Utilities;
import org.geotools.resources.i18n.Errors;
import org.geotools.resources.i18n.ErrorKeys;
/**
* Transform a set of coordinate points using a grid of localization.
* Input coordinates are index in this two-dimensional array.
* Those input coordinates (or index) should be in the range
*
* xinput = [0..width-1]
and
* yinput = [0..height-1]
inclusive,
*
* where {@code width} and {@code height} are the number of columns and
* rows in the grid of localization. Output coordinates are the values stored in
* the grid of localization at the specified index. If input coordinates (index)
* are non-integer values, then output coordinates are interpolated using a bilinear
* interpolation. If input coordinates are outside the grid range, then output
* coordinates are extrapolated.
*
* @since 2.0
* @source $URL$
* @version $Id$
* @author Remi Eve
* @author Martin Desruisseaux (IRD)
*
* @todo This class should extends {@link WarpTransform2D} and constructs a
* {@link javax.media.jai.WarpGrid} the first time the {@link WarpTransform2D#getWarp()}
* method is invoked (GEOT-522).
*/
final class LocalizationGridTransform2D extends AbstractMathTransform
implements MathTransform2D, Serializable
{
/**
* Serial number for interoperability with different versions.
*/
private static final long serialVersionUID = 1067560328828441295L;
/**
* Maximal number of iterations to try before to fail
* during an inverse transformation.
*/
private static final int MAX_ITER = 40;
/**
* Set to {@code true} for a conservative (and maybe slower) algorithm
* in {@link #inverseTransform}.
*/
private static final boolean CONSERVATIVE = true;
/**
* Set to {@code true} for forcing {@link #inverseTransform} to returns
* a value instead of throwing an exception if the transform do not converge.
* This is a temporary flag until we find why the inverse transform fails to
* converge in some case.
*/
private static final boolean MASK_NON_CONVERGENCE;
static {
String property;
try {
property = System.getProperty("org.geotools.referencing.forceConvergence", "false");
} catch (SecurityException exception) {
// We are running in an applet.
property = "false";
}
MASK_NON_CONVERGENCE = property.equalsIgnoreCase("true");
}
/**
* x (usually longitude) offset relative to an entry.
* Points are stored in {@link #grid} as {@code (x,y)} pairs.
*/
static final int X_OFFSET = 0;
/**
* y (usually latitude) offset relative to an entry.
* Points are stored in {@link #grid} as {@code (x,y)} pairs.
*/
static final int Y_OFFSET = 1;
/**
* Length of an entry in the {@link #grid} array. This lenght
* is equals to the dimension of output coordinate points.
*/
static final int CP_LENGTH = 2;
/**
* Number of grid's columns.
*/
private final int width;
/**
* Number of grid's rows.
*/
private final int height;
/**
* Grid of coordinate points.
* Points are stored as {@code (x,y)} pairs.
*/
private final double[] grid;
/**
* A global affine transform for the whole grid.
*/
private final AffineTransform global;
/**
* The inverse math transform. Will be constructed only when first requested.
*/
private transient MathTransform2D inverse;
/**
* Constructs a localization grid using the specified data.
*
* @param width Number of grid's columns.
* @param height Number of grid's rows.
* @param grid The localization grid as an array of {@code (x,y)} coordinates.
* This array is not cloned; this is the caller's responsability to ensure
* that it will not be modified as long as this transformation is strongly
* reachable.
* @param global A global affine transform for the whole grid.
*/
protected LocalizationGridTransform2D(final int width, final int height, final double[] grid,
final AffineTransform global)
{
this.width = width;
this.height = height;
this.grid = grid;
this.global = global;
}
/**
* Returns the parameter descriptors for this math transform.
*/
@Override
public ParameterDescriptorGroup getParameterDescriptors() {
return Provider.PARAMETERS;
}
/**
* Returns the dimension of input points.
*/
public int getSourceDimensions() {
return 2;
}
/**
* Returns the dimension of output points.
*/
public int getTargetDimensions() {
return 2;
}
/**
* Tests if this transform is the identity transform.
*/
@Override
public boolean isIdentity() {
return false;
}
/**
* Gets the derivative of this transform at a point.
*/
@Override
public Matrix derivative(final Point2D point) {
final AffineTransform tr = new AffineTransform();
getAffineTransform(point.getX(), point.getY(), tr);
return new Matrix2(tr.getScaleX(), tr.getShearX(),
tr.getShearY(), tr.getScaleY());
}
/**
* Transforme des coordonnées sources (généralement des index de pixels) en coordonnées
* destinations (généralement des degrés de longitude et latitude). Les transformations
* feront intervenir des interpolations linéaires si les coordonnées sources ne sont pas
* entières.
*
* @param srcPts Points d'entrée.
* @param srcOff Index du premier point d'entrée à transformer.
* @param dstPts Points de sortie.
* @param dstOff Index du premier point de sortie.
* @param numPts Nombre de points à transformer.
*/
@Override
public void transform(final float[] srcPts, int srcOff,
final float[] dstPts, int dstOff, int numPts)
{
transform(srcPts, null, srcOff, dstPts, null, dstOff, numPts);
}
/**
* Transforme des coordonnées sources (généralement des index de pixels) en coordonnées
* destinations (généralement des degrés de longitude et latitude). Les transformations
* feront intervenir des interpolations linéaires si les coordonnées sources ne sont pas
* entières.
*
* @param srcPts Points d'entrée.
* @param srcOff Index du premier point d'entrée à transformer.
* @param dstPts Points de sortie.
* @param dstOff Index du premier point de sortie.
* @param numPts Nombre de points à transformer.
*/
public void transform(final double[] srcPts, int srcOff,
final double[] dstPts, int dstOff, int numPts)
{
transform(null, srcPts, srcOff, null, dstPts, dstOff, numPts);
}
/**
* Implementation of direct transformation.
*/
private void transform(final float[] srcPts1, final double[] srcPts2, int srcOff,
final float[] dstPts1, final double[] dstPts2, int dstOff, int numPts)
{
final int minCol = 0;
final int minRow = 0;
final int maxCol = width - 2;
final int maxRow = height - 2;
int postIncrement = 0;
if (srcOff < dstOff) {
if ((srcPts2!=null) ? srcPts2==dstPts2 : srcPts1==dstPts1) {
srcOff += (numPts-1)*2;
dstOff += (numPts-1)*2;
postIncrement = -4;
}
}
while (--numPts >= 0) {
final double xi, yi;
if (srcPts2 != null) {
xi = srcPts2[srcOff++];
yi = srcPts2[srcOff++];
} else {
xi = srcPts1[srcOff++];
yi = srcPts1[srcOff++];
}
final int col = Math.max(Math.min((int)xi, maxCol), minCol);
final int row = Math.max(Math.min((int)yi, maxRow), minRow);
final int offset00 = (col + row*width)*CP_LENGTH;
final int offset01 = offset00 + CP_LENGTH*width; // Une ligne plus bas
final int offset10 = offset00 + CP_LENGTH; // Une colonne à droite
final int offset11 = offset01 + CP_LENGTH; // Une colonne à droite, une ligne plus bas
/*
* Interpole les coordonnées de destination [00]--.(x0,y0)----[10]
* sur la ligne courante (x0,y0) ainsi que | |
* sur la ligne suivante (x1,y1). Exemple | .(x,y) |
* ci-contre: les coordonnées sources sont | |
* entre crochets, et les coordonnées de la [01]--.(x1,y1)----[11]
* sortie (à calculer) sont entre parenthèses.
*/
final double x0 = linearInterpolation(col+0, grid[offset00 + X_OFFSET],
col+1, grid[offset10 + X_OFFSET], xi);
final double y0 = linearInterpolation(col+0, grid[offset00 + Y_OFFSET],
col+1, grid[offset10 + Y_OFFSET], xi);
final double x1 = linearInterpolation(col+0, grid[offset01 + X_OFFSET],
col+1, grid[offset11 + X_OFFSET], xi);
final double y1 = linearInterpolation(col+0, grid[offset01 + Y_OFFSET],
col+1, grid[offset11 + Y_OFFSET], xi);
/*
* Interpole maintenant les coordonnées (x,y) entre les deux lignes.
*/
final double xf = linearInterpolation(row, x0, row+1, x1, yi);
final double yf = linearInterpolation(row, y0, row+1, y1, yi);
if (dstPts2 != null) {
dstPts2[dstOff++] = xf;
dstPts2[dstOff++] = yf;
} else {
dstPts1[dstOff++] = (float) xf;
dstPts1[dstOff++] = (float) yf;
}
srcOff += postIncrement;
dstOff += postIncrement;
if (false) {
final java.io.PrintStream out = System.out;
out.print("TD ==> xi : "); out.print (xi);
out.print( " / yi : "); out.print (yi);
out.print(" ---> xo : "); out.print (xf);
out.print( " / yo : "); out.println(yf);
}
}
}
/**
* Interpole/extrapole entre deux points.
*
* @param x1 Coordonnee x du premier point.
* @param y1 Coordonnee y du premier point.
* @param x2 Coordonnee x du second point.
* @param y2 Coordonnee y du second point.
* @param x Position x à laquelle calculer la valeur de y.
* @return La valeur y interpolée entre les deux points.
*/
private static double linearInterpolation(final double x1, final double y1,
final double x2, final double y2, final double x)
{
return y1 + (y2-y1)/(x2-x1) * (x-x1);
}
/**
* Retourne une approximation de la transformation affine à la position indiquée.
*
* @param col Coordonnee x du point.
* @param row Coordonnee y du point.
* @param dest Matrice dans laquelle écrire la transformation affine.
*/
private void getAffineTransform(double x, double y, final AffineTransform dest) {
int col = (int) x;
int row = (int) y;
if (col > width -2) col = width -2;
if (row > height-2) row = height-2;
if (col < 0) col = 0;
if (row < 0) row = 0;
final int sgnCol;
final int sgnRow;
if (x-col > 0.5) {
sgnCol = -1;
col++;
} else sgnCol = +1;
if (y-row > 0.5) {
sgnRow = -1;
row++;
} else sgnRow = +1;
/*
* Le calcul de la transformation affine comprend 6 P00------P10
* inconnues. Sa solution recquiert donc 6 équations. | |
* Nous les obtenons en utilisant 3 points, chaque | |
* points ayant 2 coordonnées. Voir exemple ci-contre: P01----(ignoré)
*/
final int offset00 = (col + row*width)*CP_LENGTH;
final int offset01 = offset00 + sgnRow*CP_LENGTH*width;
final int offset10 = offset00 + sgnCol*CP_LENGTH;
x = grid[offset00 + X_OFFSET];
y = grid[offset00 + Y_OFFSET];
final double dxCol = (grid[offset10 + X_OFFSET] - x) * sgnCol;
final double dyCol = (grid[offset10 + Y_OFFSET] - y) * sgnCol;
final double dxRow = (grid[offset01 + X_OFFSET] - x) * sgnRow;
final double dyRow = (grid[offset01 + Y_OFFSET] - y) * sgnRow;
dest.setTransform(dxCol, dyCol, dxRow, dyRow,
x - dxCol*col - dxRow*row,
y - dyCol*col - dyRow*row);
/*
* Si l'on transforme les 3 points qui ont servit à déterminer la transformation
* affine, on devrait obtenir un résultat identique (aux erreurs d'arrondissement
* près) peu importe que l'on utilise la transformation affine ou la grille de
* localisation.
*/
assert distance(new Point(col, row ), dest) < 1E-5;
assert distance(new Point(col+sgnCol, row ), dest) < 1E-5;
assert distance(new Point(col, row+sgnRow), dest) < 1E-5;
}
/**
* Transform a point using the localization grid, transform it back using the inverse
* of the specified affine transform, and returns the distance between the source and
* the resulting point. This is used for assertions only.
*
* @param index The source point to test.
* @param tr The affine transform to test.
* @return The distance in grid coordinate. Should be close to 0.
*/
private double distance(final Point2D index, final AffineTransform tr) {
try {
Point2D geoCoord = transform(index, null);
geoCoord = tr.inverseTransform(geoCoord, geoCoord);
return geoCoord.distance(index);
} catch (TransformException exception) {
// Should not happen
throw new AssertionError(exception);
} catch (NoninvertibleTransformException exception) {
// Not impossible. What should we do? Open question...
throw new AssertionError(exception);
}
}
/**
* Apply the inverse transform to a set of points. More specifically, this method transform
* "real world" coordinates to grid coordinates. This method use an iterative algorithm for
* that purpose. A {@link TransformException} is thrown in the computation do not converge.
* The algorithm applied by this method and its callers is:
*
*