001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: *
005: * (C) 2003-2006, Geotools Project Managment Committee (PMC)
006: * (C) 2002, Institut de Recherche pour le Développement
007: *
008: * This library is free software; you can redistribute it and/or
009: * modify it under the terms of the GNU Lesser General Public
010: * License as published by the Free Software Foundation; either
011: * version 2.1 of the License, or (at your option) any later version.
012: *
013: * This library is distributed in the hope that it will be useful,
014: * but WITHOUT ANY WARRANTY; without even the implied warranty of
015: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
016: * Lesser General Public License for more details.
017: */
018: package org.geotools.referencing.operation.builder;
019:
020: // J2SE dependencies
021: import java.awt.Dimension;
022: import java.awt.Point;
023: import java.awt.Rectangle;
024: import java.awt.geom.AffineTransform;
025: import java.awt.geom.Point2D;
026: import java.awt.image.WritableRaster; // For javadoc
027: import java.util.Arrays;
028:
029: // JAI dependencies
030: import javax.media.jai.Warp; // For javadoc
031: import javax.media.jai.WarpGrid; // For javadoc
032: import javax.media.jai.WarpPolynomial; // For javadoc
033: import javax.media.jai.RasterFactory; // For javadoc
034:
035: // OpenGIS dependencies
036: import org.opengis.coverage.grid.GridGeometry; // For javadoc
037: import org.opengis.referencing.operation.MathTransform2D;
038:
039: // Geotools dependencies
040: import org.geotools.referencing.crs.DefaultDerivedCRS; // For javadoc
041: import org.geotools.referencing.crs.DefaultGeographicCRS; // For javadoc
042: import org.geotools.referencing.cs.DefaultCartesianCS; // For javadoc
043: import org.geotools.referencing.datum.DefaultGeodeticDatum; // For javadoc
044: import org.geotools.referencing.operation.DefaultOperationMethod; // For javadoc
045: import org.geotools.referencing.operation.transform.WarpTransform2D;
046: import org.geotools.referencing.operation.transform.ProjectiveTransform;
047:
048: /**
049: * A factory for {@link MathTransform2D} backed by a <cite>grid of localization</cite>.
050: * A grid of localization is a two-dimensional array of coordinate points. The grid size
051: * is {@code width} × {@code height}. Input coordinates are
052: * (<var>i</var>,<var>j</var>) index in the grid, where <var>i</var> must be in the range
053: * {@code [0..width-1]} and <var>j</var> in the range {@code [0..height-1]} inclusive.
054: * Output coordinates are the values stored in the grid of localization at the specified index.
055: * <p>
056: * The {@code LocalizationGrid} class is usefull when the
057: * "{@linkplain GridGeometry#getGridToCoordinateSystem grid to coordinate system}"
058: * transform for a coverage is not some kind of global mathematical relationship like an
059: * {@linkplain AffineTransform affine transform}. Instead, the "real world" coordinates
060: * are explicitly specified for each pixels. If the real world coordinates are know only for some
061: * pixels at a fixed interval, then a transformation can be constructed by the concatenation of
062: * an affine transform with a grid of localization.
063: * <p>
064: * After a {@code LocalizationGrid} object has been fully constructed (i.e. real world coordinates
065: * have been specified for all grid cells), a transformation from grid coordinates to "real world"
066: * coordinates can be obtained with the {@link #getMathTransform} method. If this transformation is
067: * close enough to an affine transform, then an instance of {@link AffineTransform} is returned.
068: * Otherwise, a transform backed by the localization grid is returned.
069: * <p>
070: * The example below goes through the steps of constructing a coordinate reference system for a grid
071: * coverage from its grid of localization. This example assumes that the "real world" coordinates
072: * are longitudes and latitudes on the {@linkplain DefaultGeodeticDatum#WGS84 WGS84} ellipsoid.
073: *
074: * <blockquote><table border='2' cellpadding='6'><tr><td><pre>
075: * <FONT color='#008000'>//
076: * // Constructs a localization grid of size 10×10.
077: * //</FONT>
078: * LocalizationGrid grid = new LocalizationGrid(10,10);
079: * for (int j=0; j<10; j++) {
080: * for (int i=0; i<10; i++) {
081: * double x = ...; <FONT color='#008000'>// Set longitude here</FONT>
082: * double y = ...; <FONT color='#008000'>// Set latitude here</FONT>
083: * grid.{@linkplain #setLocalizationPoint(int,int,double,double) setLocalizationPoint}(i,j,x,y);
084: * }
085: * }
086: * <FONT color='#008000'>//
087: * // Constructs the grid coordinate reference system. <var>degree</var> is the polynomial
088: * // degree (e.g. 2) for a math transform that approximately map the grid of localization.
089: * // For a more accurate (but not always better) math transform backed by the whole grid,
090: * // invokes {@linkplain #getMathTransform()} instead, or use the special value of 0 for the degree
091: * // argument.
092: * //</FONT>
093: * MathTransform2D realToGrid = grid.{@linkplain #getPolynomialTransform(int) getPolynomialTransform}(degree).inverse();
094: * CoordinateReferenceSystem realCRS = DefaultGeographicCRS.WGS84;
095: * CoordinateReferenceSystem gridCRS = new {@linkplain DefaultDerivedCRS}("The grid CRS",
096: * new {@linkplain DefaultOperationMethod#DefaultOperationMethod(MathTransform) DefaultOperationMethod}(realToGrid),
097: * realCRS, <FONT color='#008000'>// The target ("real world") CRS</FONT>
098: * realToGrid, <FONT color='#008000'>// How the grid CRS relates to the "real world" CRS</FONT>
099: * {@linkplain DefaultCartesianCS#GRID});
100: *
101: * <FONT color='#008000'>//
102: * // Constructs the grid coverage using the grid coordinate system (not the "real world"
103: * // one). It is usefull to display the coverage in its native CRS before we resample it.
104: * // Note that if the grid of localization does not define the geographic location for
105: * // all pixels, then we need to specify some affine transform in place of the call to
106: * // IdentityTransform. For example if the grid of localization defines the location of
107: * // 1 pixel, then skip 3, then defines the location of 1 pixel, etc., then the affine
108: * // transform should be AffineTransform.getScaleInstance(0.25, 0.25).
109: * //</FONT>
110: * {@linkplain WritableRaster} raster = {@linkplain RasterFactory}.createBandedRaster(DataBuffer.TYPE_FLOAT,
111: * width, height, 1, null);
112: * for (int y=0; y<height; y++) {
113: * for (int x=0; x<width; x++) {
114: * raster.setSample(x, y, 0, <cite>some_value</cite>);
115: * }
116: * }
117: * GridCoverageFactory factory = FactoryFinder.getGridCoverageFactory(null);
118: * GridCoverage coverage = factory.create("My grayscale coverage", raster, gridCRS,
119: * IdentityTransform.create(2), null, null, null, null, null);
120: * coverage.show();
121: * <FONT color='#008000'>//
122: * // Projects the coverage from its current 'gridCS' to the 'realCS'. If the grid of
123: * // localization was built from the orbit of some satellite, then the projected
124: * // coverage will tpypically have a curved aspect.
125: * //</FONT>
126: * coverage = (Coverage2D) Operations.DEFAULT.resample(coverage, realCRS);
127: * coverage.show();
128: * </pre></td></tr></table></blockquote>
129: *
130: * @since 2.4
131: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/builder/LocalizationGrid.java $
132: * @version $Id: LocalizationGrid.java 29101 2008-02-06 12:11:19Z desruisseaux $
133: * @author Remi Eve
134: * @author Martin Desruisseaux
135: * @author Alessio Fabiani
136: *
137: * @see org.opengis.referencing.crs.DerivedCRS
138: */
139: public class LocalizationGrid {
140: /**
141: * <var>x</var> (usually longitude) offset relative to an entry.
142: * Points are stored in {@link #grid} as {@code (x,y)} pairs.
143: */
144: private static final int X_OFFSET = LocalizationGridTransform2D.X_OFFSET;
145:
146: /**
147: * <var>y</var> (usually latitude) offset relative to an entry.
148: * Points are stored in {@link #grid} as {@code (x,y)} pairs.
149: */
150: private static final int Y_OFFSET = LocalizationGridTransform2D.Y_OFFSET;
151:
152: /**
153: * Length of an entry in the {@link #grid} array. This lenght
154: * is equals to the dimension of output coordinate points.
155: */
156: private static final int CP_LENGTH = LocalizationGridTransform2D.CP_LENGTH;
157:
158: /**
159: * Number of grid's columns.
160: */
161: private final int width;
162:
163: /**
164: * Number of grid's rows.
165: */
166: private final int height;
167:
168: /**
169: * Grid of coordinate points.
170: * Points are stored as {@code (x,y)} pairs.
171: */
172: private double[] grid;
173:
174: /**
175: * A global affine transform for the whole grid. This affine transform
176: * will be computed when first requested using a "least squares" fitting.
177: */
178: private transient AffineTransform global;
179:
180: /**
181: * Math transforms from grid to "real world" data for various degrees. By convention,
182: * {@code transforms[0]} is the transform backed by the whole grid. Other index are fittings
183: * using different polynomial degrees ({@code transforms[1]} for affine, {@code transforms[2]}
184: * for quadratic, <cite>etc.</cite>). Will be computed only when first needed.
185: */
186: private transient MathTransform2D[] transforms;
187:
188: /**
189: * Constructs an initially empty localization grid. All "real worlds"
190: * coordinates are initially set to {@code (NaN,NaN)}.
191: *
192: * @param width Number of grid's columns.
193: * @param height Number of grid's rows.
194: */
195: public LocalizationGrid(final int width, final int height) {
196: if (width < 2) {
197: throw new IllegalArgumentException(String.valueOf(width));
198: }
199: if (height < 2) {
200: throw new IllegalArgumentException(String.valueOf(height));
201: }
202: this .width = width;
203: this .height = height;
204: this .grid = new double[width * height * CP_LENGTH];
205: Arrays.fill(grid, Float.NaN);
206: }
207:
208: /**
209: * Calcule l'indice d'un enregistrement dans la grille.
210: *
211: * @param row Coordonnee x du point.
212: * @param col Coordonnee y du point.
213: * @return l'indice de l'enregistrement ou du point dans la matrice.
214: */
215: private int computeOffset(final int col, final int row) {
216: if (col < 0 || col >= width) {
217: throw new IndexOutOfBoundsException(String.valueOf(col));
218: }
219: if (row < 0 || row >= height) {
220: throw new IndexOutOfBoundsException(String.valueOf(row));
221: }
222: return (col + row * width) * CP_LENGTH;
223: }
224:
225: /**
226: * Returns the grid size. Grid coordinates are always in the range
227: * <code>x<sub>input</sub> = [0..width-1]</code> and
228: * <code>y<sub>input</sub> = [0..height-1]</code> inclusive.
229: */
230: public Dimension getSize() {
231: return new Dimension(width, height);
232: }
233:
234: /**
235: * Returns the "real world" coordinates for the specified grid coordinates.
236: * Grid coordinates must be integers inside this grid's range. For general
237: * transformations involving non-integer grid coordinates and/or coordinates
238: * outside this grid's range, use {@link #getMathTransform} instead.
239: *
240: * @param source The point in grid coordinates.
241: * @return target The corresponding point in "real world" coordinates.
242: * @throws IndexOutOfBoundsException If the source point is not in this grid's range.
243: */
244: public synchronized Point2D getLocalizationPoint(final Point source) {
245: final int offset = computeOffset(source.x, source.y);
246: return new Point2D.Double(grid[offset + X_OFFSET], grid[offset
247: + Y_OFFSET]);
248: }
249:
250: /**
251: * Set a point in this localization grid.
252: *
253: * @param source The point in grid coordinates.
254: * @param target The corresponding point in "real world" coordinates.
255: * @throws IndexOutOfBoundsException If the source point is not in this grid's range.
256: */
257: public void setLocalizationPoint(final Point source,
258: final Point2D target) {
259: setLocalizationPoint(source.x, source.y, target.getX(), target
260: .getY());
261: }
262:
263: /**
264: * Set a point in this localization grid.
265: *
266: * @param sourceX <var>x</var> coordinates in grid coordinates,
267: * in the range {@code [0..width-1]} inclusive.
268: * @param sourceY <var>y</var> coordinates in grid coordinates.
269: * in the range {@code [0..height-1]} inclusive.
270: * @param targetX <var>x</var> coordinates in "real world" coordinates.
271: * @param targetY <var>y</var> coordinates in "real world" coordinates.
272: * @throws IndexOutOfBoundsException If the source coordinates is not in this grid's range.
273: */
274: public synchronized void setLocalizationPoint(int sourceX,
275: int sourceY, double targetX, double targetY) {
276: final int offset = computeOffset(sourceX, sourceY);
277: notifyChange();
278: global = null;
279: grid[offset + X_OFFSET] = targetX;
280: grid[offset + Y_OFFSET] = targetY;
281: }
282:
283: /**
284: * Apply a transformation to every "real world" coordinate points in a sub-region
285: * of this grid.
286: *
287: * @param transform The transform to apply.
288: * @param region The bounding rectangle (in grid coordinate) for region where to
289: * apply the transform, or {@code null} to transform the whole grid.
290: */
291: public synchronized void transform(final AffineTransform transform,
292: final Rectangle region) {
293: assert X_OFFSET == 0 : X_OFFSET;
294: assert Y_OFFSET == 1 : Y_OFFSET;
295: assert CP_LENGTH == 2 : CP_LENGTH;
296: if (region == null) {
297: transform.transform(grid, 0, grid, 0, width * height);
298: return;
299: }
300: computeOffset(region.x, region.y); // Range check.
301: int j = region.x + region.width;
302: if (j > width) {
303: throw new IndexOutOfBoundsException(String.valueOf(j));
304: }
305: j = region.y + region.height; // Range check performed in the loop.
306: while (--j >= region.y) {
307: final int offset = computeOffset(region.x, j);
308: notifyChange();
309: transform.transform(grid, offset, grid, offset,
310: region.width);
311: }
312: global = null;
313: }
314:
315: /**
316: * Returns {@code true} if this localization grid
317: * contains at least one {@code NaN} value.
318: */
319: public synchronized boolean isNaN() {
320: for (int i = grid.length; --i >= 0;) {
321: if (Double.isNaN(grid[i])) {
322: return true;
323: }
324: }
325: return false;
326: }
327:
328: /**
329: * Returns {@code true} if all coordinates in this grid are increasing or decreasing.
330: * More specifically, returns {@code true} if the following conditions are meets:
331: * <ul>
332: * <li>Coordinates in a row must be increasing or decreasing. If {@code strict} is
333: * {@code true}, then coordinates must be strictly increasing or decreasing (i.e.
334: * equals value are not accepted). {@code NaN} values are always ignored.</li>
335: * <li>Coordinates in all rows must be increasing, or coordinates in all rows must be
336: * decreasing.</li>
337: * <li>Idem for columns (Coordinates in a columns must be increasing or decreasing,
338: * etc.).</li>
339: * </ul>
340: *
341: * <var>x</var> and <var>y</var> coordinates are tested independently.
342: *
343: * @param strict {@code true} to require strictly increasing or decreasing order,
344: * or {@code false} to accept values that are equals.
345: * @return {@code true} if coordinates are increasing or decreasing in the same
346: * direction for all rows and columns.
347: */
348: public synchronized boolean isMonotonic(final boolean strict) {
349: int orderX = INCREASING | DECREASING;
350: int orderY = INCREASING | DECREASING;
351: if (!strict) {
352: orderX |= EQUALS;
353: orderY |= EQUALS;
354: }
355: for (int i = 0; i < width; i++) {
356: final int offset = computeOffset(i, 0);
357: final int s = CP_LENGTH * width;
358: if ((orderX = testOrder(grid, offset + X_OFFSET, height, s,
359: orderX)) == 0)
360: return false;
361: if ((orderY = testOrder(grid, offset + Y_OFFSET, height, s,
362: orderY)) == 0)
363: return false;
364: }
365: orderX = INCREASING | DECREASING;
366: orderY = INCREASING | DECREASING;
367: if (!strict) {
368: orderX |= EQUALS;
369: orderY |= EQUALS;
370: }
371: for (int j = 0; j < height; j++) {
372: final int offset = computeOffset(0, j);
373: final int s = CP_LENGTH;
374: if ((orderX = testOrder(grid, offset + X_OFFSET, width, s,
375: orderX)) == 0)
376: return false;
377: if ((orderY = testOrder(grid, offset + Y_OFFSET, width, s,
378: orderY)) == 0)
379: return false;
380: }
381: return true;
382: }
383:
384: /** Constant for {@link #testOrder}. */
385: private static final int INCREASING = 1;
386: /** Constant for {@link #testOrder}. */
387: private static final int DECREASING = 2;
388: /** Constant for {@link #testOrder}. */
389: private static final int EQUALS = 4;
390:
391: /**
392: * Checks the ordering of elements in a sub-array. {@link Float#NaN} values are ignored.
393: *
394: * @param grid The {link #grid} array.
395: * @param offset The first element to test.
396: * @param num The number of elements to test.
397: * @param step The amount to increment {@code offset} in order to reach the next element.
398: * @param flags A combinaison of {@link #INCREASING}, {@link #DECREASING} and {@link #EQUALS}
399: * that specify which ordering are accepted.
400: * @return 0 if the array is unordered. Otherwise, returns {@code flags} with maybe
401: * one of {@link #INCREASING} or {@link #DECREASING} flags cleared.
402: */
403: private static int testOrder(final double[] grid, int offset,
404: int num, final int step, int flags) {
405: // We will check (num-1) combinaisons of coordinates.
406: for (--num; --num >= 0; offset += step) {
407: final double v1 = grid[offset];
408: if (Double.isNaN(v1))
409: continue;
410: while (true) {
411: final double v2 = grid[offset + step];
412: final int required, clear;
413: if (v1 == v2) {
414: required = EQUALS; // "equals" must be accepted.
415: clear = ~0; // Do not clear anything.
416: } else if (v2 > v1) {
417: required = INCREASING; // "increasing" must be accepted.
418: clear = ~DECREASING; // do not accepts "decreasing" anymore.
419: } else if (v2 < v1) {
420: required = DECREASING; // "decreasing" must be accepted.
421: clear = ~INCREASING; // do not accepts "increasing" anymore.
422: } else {
423: // 'v2' is NaN. Search for the next element.
424: if (--num < 0) {
425: return flags;
426: }
427: offset += step;
428: continue; // Mimic the "goto" statement.
429: }
430: if ((flags & required) == 0) {
431: return 0;
432: }
433: flags &= clear;
434: break;
435: }
436: }
437: return flags;
438: }
439:
440: /**
441: * Makes sure that the grid doesn't contains identical consecutive ordinates. If many
442: * consecutives ordinates are found to be identical in a row or in a column, then
443: * the first one is left inchanged and the other ones are linearly interpolated.
444: */
445: public void removeSingularities() {
446: removeSingularities(X_OFFSET, false);
447: removeSingularities(X_OFFSET, true);
448: removeSingularities(Y_OFFSET, false);
449: removeSingularities(Y_OFFSET, true);
450: }
451:
452: /**
453: * Applies a linear interpolation on consecutive identical ordinates.
454: *
455: * @param index The offset of the ordinate to test.
456: * Should be {@link #X_OFFSET} or {@link #Y_OFFSET}.
457: * @param vertical {@code true} to scan the grid vertically, or
458: * {@code false} to scan the grid horizontally.
459: */
460: private void removeSingularities(final int index,
461: final boolean vertical) {
462: final int step, val1, val2;
463: if (vertical) {
464: step = CP_LENGTH * width;
465: val1 = width;
466: val2 = height;
467: } else {
468: step = CP_LENGTH;
469: val1 = height;
470: val2 = width;
471: }
472: for (int i = 0; i < val1; i++) {
473: final int offset;
474: if (vertical) {
475: offset = computeOffset(i, 0) + index;
476: } else {
477: offset = computeOffset(0, i) + index;
478: }
479: int singularityOffset = -1;
480: for (int j = 1; j < val2; j++) {
481: final int previousOffset = offset + step * (j - 1);
482: final int currentOffset = previousOffset + step;
483: if (grid[previousOffset] == grid[currentOffset]) {
484: if (singularityOffset == -1) {
485: singularityOffset = (previousOffset == offset) ? previousOffset
486: : previousOffset - step;
487: }
488: } else if (singularityOffset != -1) {
489: final int num = (currentOffset - singularityOffset)
490: / step + 1;
491: replaceSingularity(grid, singularityOffset, num,
492: step);
493: singularityOffset = -1;
494: }
495: }
496: if (singularityOffset != -1) {
497: final int currentOffset = offset + step * (val2 - 1);
498: final int num = (currentOffset - singularityOffset)
499: / step + 1;
500: replaceSingularity(grid, singularityOffset, num, step);
501: }
502: }
503: }
504:
505: /**
506: * Replace consecutive singularity by linear values in sub-array.
507: *
508: * Example (we consider a grid of five element with singularity) :
509: *
510: * before
511: * *--*--*--*--*--*
512: * |07|08|08|08|11|
513: * *--*--*--*--*--*
514: *
515: * Params are : offset = 0, num = 5, step = 1
516: *
517: * after
518: * *--*--*--*--*--*
519: * |07|08|09|10|11|
520: * *--*--*--*--*--*
521: * | |
522: * | |
523: * linear values are
524: * computed with these
525: * values
526: *
527: * @param grid The {link #grid} array.
528: * @param offset The first element.
529: * @param num The number of element.
530: * @param step The amount to increment {@code offset} in order to reach the next element.
531: */
532: private static void replaceSingularity(final double[] grid,
533: int offset, int num, final int step) {
534: final double increment = (grid[offset + (num - 1) * step] - grid[offset])
535: / ((double) (num - 1));
536: final double value = grid[offset];
537: offset += step;
538: for (int i = 0; i < (num - 2); i++, offset += step) {
539: grid[offset] = value + (increment * (i + 1));
540: }
541: }
542:
543: /**
544: * Returns an affine transform for the whole grid. This transform is only an approximation
545: * for this localization grid. It is fitted (like "curve fitting") to grid data using the
546: * "least squares" method.
547: *
548: * @return A global affine transform as an approximation for the whole localization grid.
549: */
550: public synchronized AffineTransform getAffineTransform() {
551: if (global == null) {
552: final double[] matrix = new double[6];
553: fitPlane(X_OFFSET, matrix);
554: assert X_OFFSET == 0 : X_OFFSET;
555: fitPlane(Y_OFFSET, matrix);
556: assert Y_OFFSET == 1 : Y_OFFSET;
557: global = new AffineTransform(matrix);
558: }
559: return (AffineTransform) global.clone();
560: }
561:
562: /**
563: * Fit a plane through the longitude or latitude values. More specifically, find
564: * coefficients <var>c</var>, <var>cx</var> and <var>cy</var> for the following
565: * equation:
566: *
567: * <pre>[longitude or latitude] = c + cx*x + cy*y</pre>.
568: *
569: * where <var>x</var> and <var>cx</var> are grid coordinates.
570: * Coefficients are computed using the least-squares method.
571: *
572: * @param offset {@link X_OFFSET} for fitting longitude values, or {@link X_OFFSET}
573: * for fitting latitude values (assuming tha "real world" coordinates
574: * are longitude and latitude values).
575: * @param coeff An array of length 6 in which to store plane's coefficients.
576: * Coefficients will be store in the following order:
577: *
578: * {@code coeff[0 + offset] = cx;}
579: * {@code coeff[2 + offset] = cy;}
580: * {@code coeff[4 + offset] = c;}
581: */
582: private void fitPlane(final int offset, final double[] coeff) {
583: /*
584: * 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
585: * values. When possible, we will avoid to compute the sum inside the loop and use the
586: * following identities instead:
587: *
588: * 1 + 2 + 3 ... + n = n*(n+1)/2 (arithmetic series)
589: * 1˛ + 2˛ + 3˛ ... + n˛ = n*(n+0.5)*(n+1)/3
590: */
591: double x, y, z, xx, yy, xy, zx, zy;
592: z = zx = zy = 0; // To be computed in the loop.
593: int n = offset;
594: for (int yi = 0; yi < height; yi++) {
595: for (int xi = 0; xi < width; xi++) {
596: assert computeOffset(xi, yi) + offset == n : n;
597: final double zi = grid[n];
598: z += zi;
599: zx += zi * xi;
600: zy += zi * yi;
601: n += CP_LENGTH;
602: }
603: }
604: n = (n - offset) / CP_LENGTH;
605: assert n == width * height : n;
606: x = (n * (double) (width - 1)) / 2;
607: y = (n * (double) (height - 1)) / 2;
608: xx = (n * (width - 0.5) * (width - 1)) / 3;
609: yy = (n * (height - 0.5) * (height - 1)) / 3;
610: xy = (n * (double) ((height - 1) * (width - 1))) / 4;
611: /*
612: * Solve the following equations for cx and cy:
613: *
614: * ( zx - z*x ) = cx*(xx - x*x) + cy*(xy - x*y)
615: * ( zy - z*y ) = cx*(xy - x*y) + cy*(yy - y*y)
616: */
617: zx -= z * x / n;
618: zy -= z * y / n;
619: xx -= x * x / n;
620: xy -= x * y / n;
621: yy -= y * y / n;
622: final double den = (xy * xy - xx * yy);
623: final double cy = (zx * xy - zy * xx) / den;
624: final double cx = (zy * xy - zx * yy) / den;
625: final double c = (z - (cx * x + cy * y)) / n;
626: coeff[0 + offset] = cx;
627: coeff[2 + offset] = cy;
628: coeff[4 + offset] = c;
629: }
630:
631: /**
632: * Returns this localization grid and its inverse as warp objects. This method tries to fit a
633: * {@linkplain WarpPolynomial polynomial warp} to the gridded coordinates. The resulting Warp
634: * is wrapped into a {@link WarpTransform2D}.
635: */
636: private MathTransform2D fitWarps(final int degree) {
637: final float[] srcCoords = new float[width * height * 2];
638: final float[] dstCoords = new float[srcCoords.length];
639: int gridOffset = 0;
640: int warpOffset = 0;
641: for (int yi = 0; yi < height; yi++) {
642: for (int xi = 0; xi < width; xi++) {
643: assert gridOffset == computeOffset(xi, yi);
644: final float x = (float) grid[gridOffset + X_OFFSET];
645: final float y = (float) grid[gridOffset + Y_OFFSET];
646: if (!Float.isNaN(x) && !Float.isNaN(y)) {
647: srcCoords[warpOffset] = xi;
648: srcCoords[warpOffset + 1] = yi;
649: dstCoords[warpOffset] = x;
650: dstCoords[warpOffset + 1] = y;
651: warpOffset += 2;
652: }
653: gridOffset += CP_LENGTH;
654: }
655: }
656: return new WarpTransform2D(null, srcCoords, 0, null, dstCoords,
657: 0, warpOffset / 2, degree);
658: }
659:
660: /**
661: * Returns a math transform from grid to "real world" coordinates using a polynomial fitting
662: * of the specified degree. By convention, a {@code degree} of 0 will returns the
663: * {@linkplain #getMathTransform() math transform backed by the whole grid}. Greater values
664: * will use a fitted polynomial ({@linkplain #getAffineTransform affine transform} for
665: * degree 1, quadratic transform for degree 2, cubic transform for degree 3, etc.).
666: *
667: * @param degree The polynomial degree for the fitting, or 0 for a transform backed by the
668: * whole grid.
669: */
670: public synchronized MathTransform2D getPolynomialTransform(
671: final int degree) {
672: if (degree < 0 || degree >= WarpTransform2D.MAX_DEGREE + 1) {
673: // TODO: provides a localized error message.
674: throw new IllegalArgumentException();
675: }
676: if (transforms == null) {
677: transforms = new MathTransform2D[WarpTransform2D.MAX_DEGREE + 1];
678: }
679: if (transforms[degree] == null) {
680: final MathTransform2D tr;
681: switch (degree) {
682: case 0: {
683: // NOTE: 'grid' is not cloned. This GridLocalization's grid will need to be
684: // cloned if a "set" method is invoked after the math transform creation.
685: tr = new LocalizationGridTransform2D(width, height,
686: grid, getAffineTransform());
687: break;
688: }
689: case 1: {
690: tr = (MathTransform2D) ProjectiveTransform
691: .create(getAffineTransform());
692: break;
693: }
694: default: {
695: tr = fitWarps(degree);
696: break;
697: }
698: }
699: transforms[degree] = tr;
700: }
701: return transforms[degree];
702: }
703:
704: /**
705: * Returns a math transform from grid to "real world" coordinates. The math transform is
706: * backed by the full grid of localization. In terms of JAI's {@linkplain Warp image warp}
707: * operations, this math transform is backed by a {@link WarpGrid} while the previous methods
708: * return math transforms backed by {@link WarpPolynomial}.
709: */
710: public final MathTransform2D getMathTransform() {
711: return getPolynomialTransform(0);
712: }
713:
714: /**
715: * Notify this localization grid that a coordinate is about to be changed. This method
716: * invalidate any transforms previously created.
717: */
718: private void notifyChange() {
719: if (transforms != null) {
720: if (transforms[0] != null) {
721: // Clones is required only for the grid-backed transform.
722: grid = (double[]) grid.clone();
723: }
724: // Signal that all transforms need to be recomputed.
725: transforms = null;
726: }
727: }
728: }
|