/*
* Geotools 2 - OpenSource mapping toolkit
* (C) 2003, Geotools Project Management Committee (PMC)
* (C) 2002, Institut de Recherche pour le Développement
*
* 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; either
* version 2.1 of the License, or (at your option) any later version.
*
* 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.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
package org.geotools.gc;
// J2SE dependencies
import java.awt.Dimension;
import java.awt.Point;
import java.awt.Rectangle;
import java.awt.geom.AffineTransform;
import java.awt.geom.Point2D;
import java.util.Arrays;
import org.geotools.cs.FittedCoordinateSystem;
import org.geotools.cs.GeographicCoordinateSystem;
import org.geotools.ct.MathTransform2D;
/**
* A factory for {@link MathTransform2D} backed by a grid of localization.
* A grid of localization is a two-dimensional array (or a matrix) of coordinate points.
* The grid size is width
× height
. Input coordinates
* are (i,j) index in the grid, where i must be in the range
* [0..width-1]
and j in the range [0..height-1]
inclusive.
* Output coordinates are the values stored in the grid of localization at the specified index.
*
* The LocalizationGrid
class is usefull when the
* "{@linkplain GridGeometry#getGridToCoordinateSystem grid to coordinate system}"
* transform for a coverage is not some kind of global mathematical relationship like an
* {@linkplain AffineTransform affine transform}. Instead, the "real world" coordinates
* are explicitly specified for each pixels. If the real world coordinates are know only for some
* pixels at a fixed interval, then a transformation can be constructed by the concatenation of
* an affine transform with a grid of localization. After a LocalizationGrid
object
* has been fully constructed (i.e. real world coordinates have been specified for all grid cells),
* a transformation from grid coordinates to "real world" coordinates can be obtained with
* the {@link #getMathTransform} method. If this transformation is close enough to an affine
* transform, then an instance of {@link AffineTransform} is returned. Otherwise, a transform
* backed by the localization grid is returned.
*
* The example below goes through the steps of constructing a coordinate system for a grid coverage
* from its grid of localization. This example assumes that the "real world" coordinates are
* longitudes and latitudes on the {@linkplain GeographicCoordinateSystem#WGS84 WGS84} ellipsoid.
*
*
* * @version $Id$ * @author Remi Eve * @author Martin Desruisseaux * * @see FittedCoordinateSystem * * @deprecated Replaced by {@link org.geotools.referencing.operation.transform.LocalizationGrid}. */ public class LocalizationGrid { /** * x (usually longitude) offset relative to an entry. * Points are stored in {@link #grid} as* // * // Constructs a localization grid of size 10×10. * // * LocalizationGrid grid = new LocalizationGrid(10,10); * for (int j=0; j<10; j++) { * for (int i=0; i<10; i++) { * double x = ...; // Set longitude here * double y = ...; // Set latitude here * grid.{@link #setLocalizationPoint(int,int,double,double) setLocalizationPoint}(i,j,x,y); * } * } * // * // Constructs the grid coordinate system. * // * CoordinateSystem realCS = GeographicCoordinateSystem.WGS84; * CoordinateSystem gridCS = new {@link FittedCoordinateSystem}("The grid CS", realCS, grid.{@link #getMathTransform()}, * new AxisInfo[] { * AxisInfo.X, * AxisInfo.Y}); * // * // Constructs the grid coverage using the grid coordinate system (not the "real world" * // one). It is usefull to display the coverage in its native CS before we resample it. * // Note that if the grid of localization does not define the geographic location for * // all pixels, then we need to specify some affine transform in place of the IDENTITY * // argument. For example if the grid of localization defines the location of 1 pixel, * // then skip 3, then defines the location of 1 pixel, etc., then the affine transform * // should be AffineTransform.getScaleInstance(0.25, 0.25). * // * GridCoverage coverage; * coverage = new GridCoverage("The grid coverage", theRaster, gridCS, * MathTransform2D.IDENTITY, ...); * FrameFactory.show(coverage); * // * // Project the coverage from its current 'gridCS' to the 'realCS'. If the grid of * // localization was built from the orbit of some satellite, then the projected * // coverage will tpypically have a curved aspect. * // * GridCoverageProcessor processor = GridCoverageProcessor.getDefault(); * coverage = processor.doOperation("Resample", coverage, * "CoordinateSystem", realCS); * FrameFactory.show(coverage); *
(x,y)
pairs.
*/
private static final int X_OFFSET = LocalizationGridTransform2D.X_OFFSET;
/**
* y (usually latitude) offset relative to an entry.
* Points are stored in {@link #grid} as (x,y)
pairs.
*/
private static final int Y_OFFSET = LocalizationGridTransform2D.Y_OFFSET;
/**
* Length of an entry in the {@link #grid} array. This lenght
* is equals to the dimension of output coordinate points.
*/
private static final int CP_LENGTH = LocalizationGridTransform2D.CP_LENGTH;
/**
* 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 (x,y)
pairs.
*/
private double[] grid;
/**
* A global affine transform for the whole grid. This affine transform
* will be computed when first requested using a "least squares" fitting.
*/
private transient AffineTransform global;
/**
* A math transform from grid to "real world" data.
* Will be computed only when first needed.
*/
private transient MathTransform2D transform;
/**
* Construct an initially empty localization grid. All "real worlds"
* coordinates are initially set to (NaN,NaN)
.
*
* @param width Number of grid's columns.
* @param height Number of grid's rows.
*/
public LocalizationGrid(final int width, final int height) {
if (width < 2) {
throw new IllegalArgumentException(String.valueOf(width));
}
if (height < 2) {
throw new IllegalArgumentException(String.valueOf(height));
}
this.width = width;
this.height = height;
this.grid = new double[width * height * CP_LENGTH];
Arrays.fill(grid, Float.NaN);
}
/**
* Calcule l'indice d'un enregistrement dans la grille.
*
* @param row Coordonnee x du point.
* @param col Coordonnee y du point.
* @return l'indice de l'enregistrement ou du point dans la matrice.
*/
private int computeOffset(final int col, final int row) {
if (col<0 || col>=width) {
throw new IndexOutOfBoundsException(String.valueOf(col));
}
if (row<0 || row>=height) {
throw new IndexOutOfBoundsException(String.valueOf(row));
}
return (col + row * width) * CP_LENGTH;
}
/**
* Returns the grid size. Grid coordinates are always in the range
* xinput = [0..width-1]
and
* yinput = [0..height-1]
inclusive.
*/
public Dimension getSize() {
return new Dimension(width, height);
}
/**
* Returns the "real world" coordinates for the specified grid coordinates.
* Grid coordinates must be integers inside this grid's range. For general
* transformations involving non-integer grid coordinates and/or coordinates
* outside this grid's range, use {@link #getMathTransform} instead.
*
* @param source The point in grid coordinates.
* @return target The corresponding point in "real world" coordinates.
* @throws IndexOutOfBoundsException If the source point is not in this grid's range.
*/
public synchronized Point2D getLocalizationPoint(final Point source) {
final int offset = computeOffset(source.x, source.y);
return new Point2D.Double(grid[offset + X_OFFSET],
grid[offset + Y_OFFSET]);
}
/**
* Set a point in this localization grid.
*
* @param source The point in grid coordinates.
* @param target The corresponding point in "real world" coordinates.
* @throws IndexOutOfBoundsException If the source point is not in this grid's range.
*/
public void setLocalizationPoint(final Point source, final Point2D target) {
setLocalizationPoint(source.x, source.y, target.getX(), target.getY());
}
/**
* Set a point in this localization grid.
*
* @param sourceX x coordinates in grid coordinates,
* in the range [0..width-1]
inclusive.
* @param sourceY y coordinates in grid coordinates.
* in the range [0..height-1]
inclusive.
* @param targetX x coordinates in "real world" coordinates.
* @param targetY y coordinates in "real world" coordinates.
* @throws IndexOutOfBoundsException If the source coordinates is not in this grid's range.
*/
public synchronized void setLocalizationPoint(int sourceX, int sourceY,
double targetX, double targetY)
{
final int offset = computeOffset(sourceX, sourceY);
if (transform != null) {
transform = null;
grid = (double[]) grid.clone();
}
global = null;
grid[offset + X_OFFSET] = targetX;
grid[offset + Y_OFFSET] = targetY;
}
/**
* Apply a transformation to every "real world" coordinate points in a sub-region
* of this grid.
*
* @param transform The transform to apply.
* @param region The bounding rectangle (in grid coordinate) for region where to
* apply the transform, or null
to transform the whole grid.
*/
public synchronized void transform(final AffineTransform transform, final Rectangle region) {
assert X_OFFSET == 0 : X_OFFSET;
assert Y_OFFSET == 1 : Y_OFFSET;
assert CP_LENGTH == 2 : CP_LENGTH;
if (region == null) {
transform.transform(grid, 0, grid, 0, width*height);
return;
}
computeOffset(region.x, region.y); // Range check.
int j = region.x + region.width;
if (j > width) {
throw new IndexOutOfBoundsException(String.valueOf(j));
}
j = region.y + region.height; // Range check performed in the loop.
while (--j >= region.y) {
final int offset = computeOffset(region.x, j);
if (this.transform != null) {
this.transform = null;
grid = (double[]) grid.clone();
}
transform.transform(grid, offset, grid, offset, region.width);
}
global = null;
}
/**
* Returns true
if this localization grid
* contains at least one NaN
value.
*/
public synchronized boolean isNaN() {
for (int i=grid.length; --i>=0;) {
if (Double.isNaN(grid[i])) {
return true;
}
}
return false;
}
/**
* Returns true
if all coordinates in this grid are increasing or decreasing.
* More specifically, returns true
if the following conditions are meets:
* strict
is
* true
, then coordinates must be strictly increasing or decreasing (i.e.
* equals value are not accepted). NaN
values are always ignored.true
to require strictly increasing or decreasing order,
* or false
to accept values that are equals.
* @return true
if coordinates are increasing or decreasing in the same
* direction for all rows and columns.
*/
public synchronized boolean isMonotonic(final boolean strict) {
int orderX = INCREASING|DECREASING;
int orderY = INCREASING|DECREASING;
if (!strict) {
orderX |= EQUALS;
orderY |= EQUALS;
}
for (int i=0; iflags
with maybe
* one of {@link #INCREASING} or {@link #DECREASING} flags cleared.
*/
private static int testOrder(final double[] grid, int offset, int num, final int step, int flags)
{
// We will check (num-1) combinaisons of coordinates.
for (--num; --num>=0; offset += step) {
final double v1 = grid[offset];
if (Double.isNaN(v1)) continue;
while (true) {
final double v2 = grid[offset + step];
final int required, clear;
if (v1 == v2) {
required = EQUALS; // "equals" must be accepted.
clear = ~0; // Do not clear anything.
} else if (v2 > v1) {
required = INCREASING; // "increasing" must be accepted.
clear = ~DECREASING; // do not accepts "decreasing" anymore.
} else if (v2 < v1) {
required = DECREASING; // "decreasing" must be accepted.
clear = ~INCREASING; // do not accepts "increasing" anymore.
} else {
// 'v2' is NaN. Search for the next element.
if (--num < 0) {
return flags;
}
offset += step;
continue; // Mimic the "goto" statement.
}
if ((flags & required) == 0) {
return 0;
}
flags &= clear;
break;
}
}
return flags;
}
/**
* Make sure that the grid doesn't contains identical consecutive ordinates. If many
* consecutives ordinates are found to be identical in a row or in a column, then
* the first one is left inchanged and the other ones are linearly interpolated.
*/
public void removeSingularities() {
removeSingularities(X_OFFSET, false);
removeSingularities(X_OFFSET, true );
removeSingularities(Y_OFFSET, false);
removeSingularities(Y_OFFSET, true );
}
/**
* Apply a linear interpolation on consecutive identical ordinates.
*
* @param index The offset of the ordinate to test.
* Should be {@link #X_OFFSET} or {@link #Y_OFFSET}.
* @param vertical true
to scan the grid vertically, or
* false
to scan the grid horizontally.
*/
private void removeSingularities(final int index, final boolean vertical) {
final int step, val1, val2;
if (vertical) {
step = CP_LENGTH*width;
val1 = width;
val2 = height;
} else {
step = CP_LENGTH;
val1 = height;
val2 = width;
}
for (int i=0; i[longitude or latitude] = c + cx*x + cy*y. * * where x and cx are grid coordinates. * Coefficients are computed using the least-squares method. * * @param offset {@link X_OFFSET} for fitting longitude values, or {@link X_OFFSET} * for fitting latitude values (assuming tha "real world" coordinates * are longitude and latitude values). * @param coeff An array of length 6 in which to store plane's coefficients. * Coefficients will be store in the following order: * *
coeff[0 + offset] = cx;
* coeff[2 + offset] = cy;
* coeff[4 + offset] = c;
*/
private void fitPlane(final int offset, final double[] coeff) {
/*
* Compute the sum of x, y and z values. Compute also the sum of x*x, y*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:
*
* 1 + 2 + 3 ... + n = n*(n+1)/2 (arithmetic series)
* 1˛ + 2˛ + 3˛ ... + n˛ = n*(n+0.5)*(n+1)/3
*/
double x,y,z, xx,yy, xy, zx,zy;
z = zx = zy = 0; // To be computed in the loop.
int n=offset;
for (int yi=0; yi