/* * 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. * *
 * //
 * // 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);
 * 
* * @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 (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: * * * x and y coordinates are tested independently. * * @param strict 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; ioffset in order to reach the next element. * @param flags A combinaison of {@link #INCREASING}, {@link #DECREASING} and {@link #EQUALS} * that specify which ordering are accepted. * @return 0 if the array is unordered. Otherwise, returns flags 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; ioffset in order to reach the next element. */ private static void replaceSingularity(final double[] grid, int offset, int num, final int step) { final double increment = (grid[offset+(num-1)*step] - grid[offset])/((double)(num-1)); final double value = grid[offset]; offset+= step; for (int i=0; i<(num-2); i++, offset += step) { grid[offset] = value + (increment * (i+1)); } } /** * Returns an affine transform for the whole grid. This transform is only an approximation * for this localization grid. It is fitted (like "curve fitting") to grid data using the * "least squares" method. * * @return A global affine transform as an approximation for the whole localization grid. */ public synchronized AffineTransform getAffineTransform() { if (global == null) { final double[] matrix = new double[6]; fitPlane(X_OFFSET, matrix); assert X_OFFSET==0 : X_OFFSET; fitPlane(Y_OFFSET, matrix); assert Y_OFFSET==1 : Y_OFFSET; global = new AffineTransform(matrix); } return (AffineTransform) global.clone(); } /** * Fit a plane through the longitude or latitude values. More specifically, find * coefficients c, cx and cy for the following * equation: * *
[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