001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2006, Geotools Project Managment Committee (PMC)
005: * (C) 2001, Institut de Recherche pour le Développement
006: *
007: * This library is free software; you can redistribute it and/or
008: * modify it under the terms of the GNU Lesser General Public
009: * License as published by the Free Software Foundation; either
010: * version 2.1 of the License, or (at your option) any later version.
011: *
012: * This library is distributed in the hope that it will be useful,
013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
015: * Lesser General Public License for more details.
016: */
017: package org.geotools.coverage.grid;
018:
019: // J2SE dependencies
020: import java.awt.Rectangle;
021: import java.awt.geom.Point2D;
022: import java.awt.image.RenderedImage;
023: import java.lang.reflect.Array;
024: import java.util.ArrayList;
025: import java.util.List;
026:
027: // JAI dependencies
028: import javax.media.jai.Interpolation;
029: import javax.media.jai.InterpolationNearest;
030: import javax.media.jai.iterator.RectIter;
031: import javax.media.jai.iterator.RectIterFactory;
032:
033: // OpenGIS dependencies
034: import org.opengis.coverage.CannotEvaluateException;
035: import org.opengis.coverage.PointOutsideCoverageException;
036: import org.opengis.referencing.operation.MathTransform2D;
037: import org.opengis.referencing.operation.NoninvertibleTransformException;
038: import org.opengis.referencing.operation.TransformException;
039:
040: /**
041: * A grid coverage using an {@linkplain Interpolation interpolation} for evaluating points. This
042: * interpolator is not used for {@linkplain InterpolationNearest nearest-neighbor interpolation}
043: * (use the plain {@link GridCoverage2D} class for that). It should work for other kinds of
044: * interpolation however.
045: *
046: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/coverage/src/main/java/org/geotools/coverage/grid/Interpolator2D.java $
047: * @version $Id: Interpolator2D.java 22710 2006-11-12 18:04:54Z desruisseaux $
048: * @author Martin Desruisseaux
049: *
050: * @since 2.2
051: */
052: public final class Interpolator2D extends GridCoverage2D {
053: /**
054: * The greatest value smaller than 1 representable as a {@code float} number.
055: * This value can be obtained with {@code org.geotools.resources.XMath.previous(1f)}.
056: */
057: private static final float ONE_EPSILON = 0.99999994f;
058:
059: /**
060: * Default interpolations, in preference order. Will be constructed only when first needed.
061: */
062: private static Interpolation[] DEFAULTS;
063:
064: /**
065: * Transform from "real world" coordinates to grid coordinates.
066: * This transform maps coordinates to pixel <em>centers</em>.
067: */
068: private final MathTransform2D toGrid;
069:
070: /**
071: * The interpolation method.
072: */
073: private final Interpolation interpolation;
074:
075: /**
076: * Second interpolation method to use if this one failed. May be {@code null} if there
077: * is no fallback. By convention, {@code this} means that interpolation should fallback
078: * on {@code super.evaluate(...)} (i.e. nearest neighbor).
079: */
080: private final Interpolator2D fallback;
081:
082: /**
083: * Image bounds. Bounds have been reduced by {@link Interpolation}'s padding.
084: */
085: private final int xmin, ymin, xmax, ymax;
086:
087: /**
088: * Interpolation padding.
089: */
090: private final int top, left;
091:
092: /**
093: * The interpolation bounds. Interpolation will use pixel inside this rectangle.
094: * This rectangle is passed as an argument to {@link RectIterFactory}.
095: */
096: private final Rectangle bounds;
097:
098: /**
099: * Arrays to use for passing arguments to interpolation.
100: * This array will be constructed only when first needed.
101: */
102: private transient double[][] doubles;
103:
104: /**
105: * Arrays to use for passing arguments to interpolation.
106: * This array will be constructed only when first needed.
107: */
108: private transient float[][] floats;
109:
110: /**
111: * Arrays to use for passing arguments to interpolation.
112: * This array will be constructed only when first needed.
113: */
114: private transient int[][] ints;
115:
116: /**
117: * Constructs a new interpolator using default interpolations.
118: *
119: * @param coverage The coverage to interpolate.
120: */
121: public static GridCoverage2D create(final GridCoverage2D coverage) {
122: // No need to synchronize: not a big deal if two arrays are created.
123: if (DEFAULTS == null) {
124: DEFAULTS = new Interpolation[] {
125: Interpolation
126: .getInstance(Interpolation.INTERP_BICUBIC),
127: Interpolation
128: .getInstance(Interpolation.INTERP_BILINEAR),
129: Interpolation
130: .getInstance(Interpolation.INTERP_NEAREST) };
131: }
132: return create(coverage, DEFAULTS);
133: }
134:
135: /**
136: * Constructs a new interpolator for a single interpolation.
137: *
138: * @param coverage The coverage to interpolate.
139: * @param interpolation The interpolation to use.
140: */
141: public static GridCoverage2D create(final GridCoverage2D coverage,
142: final Interpolation interpolation) {
143: return create(coverage, new Interpolation[] { interpolation });
144: }
145:
146: /**
147: * Constructs a new interpolator for an interpolation and its fallbacks. The fallbacks
148: * are used if the primary interpolation failed because of {@linkplain Float#NaN NaN}
149: * values in the interpolated point neighbor.
150: *
151: * @param coverage The coverage to interpolate.
152: * @param interpolations The interpolation to use and its fallback (if any).
153: */
154: public static GridCoverage2D create(GridCoverage2D coverage,
155: final Interpolation[] interpolations) {
156: while (coverage instanceof Interpolator2D) {
157: coverage = ((Interpolator2D) coverage).getSource();
158: }
159: if (interpolations.length == 0
160: || (interpolations[0] instanceof InterpolationNearest)) {
161: return coverage;
162: }
163: return new Interpolator2D(coverage, interpolations, 0);
164: }
165:
166: /**
167: * Constructs a new interpolator for the specified interpolation.
168: *
169: * @param coverage The coverage to interpolate.
170: * @param interpolations The interpolations to use and its fallback
171: * (if any). This array must have at least 1 element.
172: * @param index The index of interpolation to use in the {@code interpolations} array.
173: */
174: private Interpolator2D(final GridCoverage2D coverage,
175: final Interpolation[] interpolations, final int index) {
176: super (null, coverage);
177: this .interpolation = interpolations[index];
178: if (index + 1 < interpolations.length) {
179: if (interpolations[index + 1] instanceof InterpolationNearest) {
180: // By convention, 'fallback==this' is for 'super.evaluate(...)'
181: // (i.e. "NearestNeighbor").
182: this .fallback = this ;
183: } else {
184: this .fallback = new Interpolator2D(coverage,
185: interpolations, index + 1);
186: }
187: } else {
188: this .fallback = null;
189: }
190: /*
191: * Computes the affine transform from "real world" coordinates to grid coordinates.
192: * This transform maps coordinates to pixel <em>centers</em>. If this transform has
193: * already be created during fallback construction, reuse the fallback's instance
194: * instead of creating a new identical one.
195: */
196: if (fallback != null && fallback != this ) {
197: this .toGrid = fallback.toGrid;
198: } else
199: try {
200: final MathTransform2D transform = gridGeometry
201: .getGridToCRS2D();
202: // Note: If we want nearest-neighbor interpolation, we need to add the
203: // following line (assuming the transform is an 'AffineTransform'):
204: //
205: // transform.translate(-0.5, -0.5);
206: //
207: // This is because we need to cancel the last 'translate(0.5, 0.5)' that appears
208: // in GridGeometry's constructor (we must remember that OpenGIS's transform maps
209: // pixel CENTER, while JAI transforms maps pixel UPPER LEFT corner). For exemple
210: // the (12.4, 18.9) coordinates still lies on the [12,9] pixel. Since the JAI's
211: // nearest-neighbor interpolation use 'Math.floor' operation instead of
212: // 'Math.round', we must follow this convention.
213: //
214: // For other kinds of interpolation, we want to maps pixel values to pixel center.
215: // For example, coordinate (12.5, 18.5) (in floating-point coordinates) lies at
216: // the center of pixel [12,18] (in integer coordinates); the evaluated value
217: // should be the exact pixel's value. On the other hand, coordinate (12.5, 19)
218: // (in floating-point coordinates) lies exactly at the edge between pixels
219: // [12,19] and [12,20]; the evaluated value should be a mid-value between those
220: // two pixels. If we want center of mass located at pixel centers, we must keep
221: // the (0.5, 0.5) translation provided by 'GridGeometry' for interpolation other
222: // than nearest-neighbor.
223: toGrid = (MathTransform2D) transform.inverse();
224: } catch (NoninvertibleTransformException exception) {
225: final IllegalArgumentException e = new IllegalArgumentException();
226: e.initCause(exception);
227: throw e;
228: }
229:
230: final int left = interpolation.getLeftPadding();
231: final int right = interpolation.getRightPadding();
232: final int top = interpolation.getTopPadding();
233: final int bottom = interpolation.getBottomPadding();
234:
235: this .top = top;
236: this .left = left;
237:
238: final int x = image.getMinX();
239: final int y = image.getMinY();
240:
241: this .xmin = x + left;
242: this .ymin = y + top;
243: this .xmax = x + image.getWidth() - right;
244: this .ymax = y + image.getHeight() - bottom;
245:
246: bounds = new Rectangle(0, 0, interpolation.getWidth(),
247: interpolation.getHeight());
248: }
249:
250: /**
251: * Returns the source grid coverage.
252: */
253: private GridCoverage2D getSource() {
254: final List sources = getSources();
255: assert sources.size() == 1 : sources;
256: return (GridCoverage2D) sources.get(0);
257: }
258:
259: /**
260: * Invoked by {@link #geophysics(boolean)} when the packed or geophysics companion of this
261: * grid coverage need to be created. This method apply to the new grid coverage the same
262: * interpolation than this grid coverage.
263: *
264: * @param geo {@code true} to get a grid coverage with sample values equals to
265: * geophysics values, or {@code false} to get the packed version.
266: * @return The newly created grid coverage.
267: */
268: protected GridCoverage2D createGeophysics(final boolean geo) {
269: return create(getSource().geophysics(geo), getInterpolations());
270: }
271:
272: /**
273: * Returns interpolations. The first array's element is the interpolation for
274: * this grid coverage. Other elements (if any) are fallbacks.
275: */
276: public Interpolation[] getInterpolations() {
277: final List interp = new ArrayList();
278: Interpolator2D scan = this ;
279: do {
280: interp.add(interpolation);
281: if (scan.fallback == scan) {
282: interp.add(Interpolation
283: .getInstance(Interpolation.INTERP_NEAREST));
284: break;
285: }
286: scan = scan.fallback;
287: } while (scan != null);
288: return (Interpolation[]) interp
289: .toArray(new Interpolation[interp.size()]);
290: }
291:
292: /**
293: * Returns the primary interpolation used by this {@code Interpolator2D}.
294: */
295: public Interpolation getInterpolation() {
296: return interpolation;
297: }
298:
299: /**
300: * Returns a sequence of integer values for a given two-dimensional point in the coverage.
301: *
302: * @param coord The coordinate point where to evaluate.
303: * @param dest An array in which to store values, or {@code null}.
304: * @return An array containing values.
305: * @throws CannotEvaluateException if the values can't be computed at the specified coordinate.
306: * More specifically, {@link PointOutsideCoverageException} is thrown if the evaluation
307: * failed because the input point has invalid coordinates.
308: */
309: public int[] evaluate(final Point2D coord, int[] dest)
310: throws CannotEvaluateException {
311: if (fallback != null) {
312: dest = super .evaluate(coord, dest);
313: }
314: try {
315: final Point2D pixel = toGrid.transform(coord, null);
316: final double x = pixel.getX();
317: final double y = pixel.getY();
318: if (!Double.isNaN(x) && !Double.isNaN(y)) {
319: dest = interpolate(x, y, dest, 0, image.getNumBands());
320: if (dest != null) {
321: return dest;
322: }
323: }
324: } catch (TransformException exception) {
325: throw new CannotEvaluateException(/*coord*/null, exception); // TODO
326: }
327: throw new PointOutsideCoverageException(
328: pointOutsideCoverage(coord));
329: }
330:
331: /**
332: * Returns a sequence of float values for a given two-dimensional point in the coverage.
333: *
334: * @param coord The coordinate point where to evaluate.
335: * @param dest An array in which to store values, or {@code null}.
336: * @return An array containing values.
337: * @throws CannotEvaluateException if the values can't be computed at the specified coordinate.
338: * More specifically, {@link PointOutsideCoverageException} is thrown if the evaluation
339: * failed because the input point has invalid coordinates.
340: */
341: public float[] evaluate(final Point2D coord, float[] dest)
342: throws CannotEvaluateException {
343: if (fallback != null) {
344: dest = super .evaluate(coord, dest);
345: }
346: try {
347: final Point2D pixel = toGrid.transform(coord, null);
348: final double x = pixel.getX();
349: final double y = pixel.getY();
350: if (!Double.isNaN(x) && !Double.isNaN(y)) {
351: dest = interpolate(x, y, dest, 0, image.getNumBands());
352: if (dest != null) {
353: return dest;
354: }
355: }
356: } catch (TransformException exception) {
357: throw new CannotEvaluateException(/*coord*/null, exception); // TODO
358: }
359: throw new PointOutsideCoverageException(
360: pointOutsideCoverage(coord));
361: }
362:
363: /**
364: * Returns a sequence of double values for a given two-dimensional point in the coverage.
365: *
366: * @param coord The coordinate point where to evaluate.
367: * @param dest An array in which to store values, or {@code null}.
368: * @return An array containing values.
369: * @throws CannotEvaluateException if the values can't be computed at the specified coordinate.
370: * More specifically, {@link PointOutsideCoverageException} is thrown if the evaluation
371: * failed because the input point has invalid coordinates.
372: */
373: public double[] evaluate(final Point2D coord, double[] dest)
374: throws CannotEvaluateException {
375: if (fallback != null) {
376: dest = super .evaluate(coord, dest);
377: }
378: try {
379: final Point2D pixel = toGrid.transform(coord, null);
380: final double x = pixel.getX();
381: final double y = pixel.getY();
382: if (!Double.isNaN(x) && !Double.isNaN(y)) {
383: dest = interpolate(x, y, dest, 0, image.getNumBands());
384: if (dest != null) {
385: return dest;
386: }
387: }
388: } catch (TransformException exception) {
389: throw new CannotEvaluateException(/*coord*/null, exception); // TODO
390: }
391: throw new PointOutsideCoverageException(
392: pointOutsideCoverage(coord));
393: }
394:
395: /**
396: * Interpolates at the specified position. If {@code fallback!=null},
397: * then {@code dest} <strong>must</strong> have been initialized with
398: * {@code super.evaluate(...)} prior to invoking this method.
399: *
400: * @param x The x position in pixel's coordinates.
401: * @param y The y position in pixel's coordinates.
402: * @param dest The destination array, or null.
403: * @param band The first band's index to interpolate.
404: * @param bandUp The last band's index+1 to interpolate.
405: * @return {@code null} if point is outside grid coverage.
406: */
407: private synchronized int[] interpolate(final double x,
408: final double y, int[] dest, int band, final int bandUp) {
409: final double x0 = Math.floor(x);
410: final double y0 = Math.floor(y);
411: final int ix = (int) x0;
412: final int iy = (int) y0;
413: if (!(ix >= xmin && ix < xmax && iy >= ymin && iy < ymax)) {
414: if (fallback == null)
415: return null;
416: if (fallback == this )
417: return dest; // super.evaluate(...) succeed prior to this call.
418: return fallback.interpolate(x, y, dest, band, bandUp);
419: }
420: /*
421: * Creates buffers, if not already created.
422: */
423: int[][] samples = ints;
424: if (samples == null) {
425: final int rowCount = interpolation.getHeight();
426: final int colCount = interpolation.getWidth();
427: ints = samples = new int[rowCount][];
428: for (int i = 0; i < rowCount; i++) {
429: samples[i] = new int[colCount];
430: }
431: }
432: if (dest == null) {
433: dest = new int[bandUp];
434: }
435: /*
436: * Builds up a RectIter and use it for interpolating all bands.
437: * There is very few points, so the cost of creating a RectIter
438: * may be important. But it seems to still lower than query tiles
439: * many time (which may involve more computation than necessary).
440: */
441: bounds.x = ix - left;
442: bounds.y = iy - top;
443: final RectIter iter = RectIterFactory.create(image, bounds);
444: for (; band < bandUp; band++) {
445: iter.startLines();
446: int j = 0;
447: do {
448: iter.startPixels();
449: final int[] row = samples[j++];
450: int i = 0;
451: do {
452: row[i++] = iter.getSample(band);
453: } while (!iter.nextPixelDone());
454: assert i == row.length;
455: } while (!iter.nextLineDone());
456: assert j == samples.length;
457: final int xfrac = (int) ((x - x0) * (1 << interpolation
458: .getSubsampleBitsH()));
459: final int yfrac = (int) ((y - y0) * (1 << interpolation
460: .getSubsampleBitsV()));
461: dest[band] = interpolation.interpolate(samples, xfrac,
462: yfrac);
463: }
464: return dest;
465: }
466:
467: /**
468: * Interpolates at the specified position. If {@code fallback!=null},
469: * then {@code dest} <strong>must</strong> have been initialized with
470: * {@code super.evaluate(...)} prior to invoking this method.
471: *
472: * @param x The x position in pixel's coordinates.
473: * @param y The y position in pixel's coordinates.
474: * @param dest The destination array, or null.
475: * @param band The first band's index to interpolate.
476: * @param bandUp The last band's index+1 to interpolate.
477: * @return {@code null} if point is outside grid coverage.
478: */
479: private synchronized float[] interpolate(final double x,
480: final double y, float[] dest, int band, final int bandUp) {
481: final double x0 = Math.floor(x);
482: final double y0 = Math.floor(y);
483: final int ix = (int) x0;
484: final int iy = (int) y0;
485: if (!(ix >= xmin && ix < xmax && iy >= ymin && iy < ymax)) {
486: if (fallback == null)
487: return null;
488: if (fallback == this )
489: return dest; // super.evaluate(...) succeed prior to this call.
490: return fallback.interpolate(x, y, dest, band, bandUp);
491: }
492: /*
493: * Create buffers, if not already created.
494: */
495: float[][] samples = floats;
496: if (samples == null) {
497: final int rowCount = interpolation.getHeight();
498: final int colCount = interpolation.getWidth();
499: floats = samples = new float[rowCount][];
500: for (int i = 0; i < rowCount; i++) {
501: samples[i] = new float[colCount];
502: }
503: }
504: if (dest == null) {
505: dest = new float[bandUp];
506: }
507: /*
508: * Builds up a RectIter and use it for interpolating all bands.
509: * There is very few points, so the cost of creating a RectIter
510: * may be important. But it seems to still lower than query tiles
511: * many time (which may involve more computation than necessary).
512: */
513: bounds.x = ix - left;
514: bounds.y = iy - top;
515: final RectIter iter = RectIterFactory.create(image, bounds);
516: for (; band < bandUp; band++) {
517: iter.startLines();
518: int j = 0;
519: do {
520: iter.startPixels();
521: final float[] row = samples[j++];
522: int i = 0;
523: do {
524: row[i++] = iter.getSampleFloat(band);
525: } while (!iter.nextPixelDone());
526: assert i == row.length;
527: } while (!iter.nextLineDone());
528: assert j == samples.length;
529: float dx = (float) (x - x0);
530: if (dx == 1)
531: dx = ONE_EPSILON;
532: float dy = (float) (y - y0);
533: if (dy == 1)
534: dy = ONE_EPSILON;
535: final float value = interpolation.interpolate(samples, dx,
536: dy);
537: if (Float.isNaN(value)) {
538: if (fallback == this )
539: continue; // 'dest' was set by 'super.evaluate(...)'.
540: if (fallback != null) {
541: fallback.interpolate(x, y, dest, band, band + 1);
542: continue;
543: }
544: // If no fallback was specified, then 'dest' is not required to
545: // have been initialized. It may contains random value. Set it
546: // to the NaN value...
547: }
548: dest[band] = value;
549: }
550: return dest;
551: }
552:
553: /**
554: * Interpolate at the specified position. If {@code fallback!=null},
555: * then {@code dest} <strong>must</strong> have been initialized with
556: * {@code super.evaluate(...)} prior to invoking this method.
557: *
558: * @param x The x position in pixel's coordinates.
559: * @param y The y position in pixel's coordinates.
560: * @param dest The destination array, or null.
561: * @param band The first band's index to interpolate.
562: * @param bandUp The last band's index+1 to interpolate.
563: * @return {@code null} if point is outside grid coverage.
564: */
565: private synchronized double[] interpolate(final double x,
566: final double y, double[] dest, int band, final int bandUp) {
567: final double x0 = Math.floor(x);
568: final double y0 = Math.floor(y);
569: final int ix = (int) x0;
570: final int iy = (int) y0;
571: if (!(ix >= xmin && ix < xmax && iy >= ymin && iy < ymax)) {
572: if (fallback == null)
573: return null;
574: if (fallback == this )
575: return dest; // super.evaluate(...) succeed prior to this call.
576: return fallback.interpolate(x, y, dest, band, bandUp);
577: }
578: /*
579: * Creates buffers, if not already created.
580: */
581: double[][] samples = doubles;
582: if (samples == null) {
583: final int rowCount = interpolation.getHeight();
584: final int colCount = interpolation.getWidth();
585: doubles = samples = new double[rowCount][];
586: for (int i = 0; i < rowCount; i++) {
587: samples[i] = new double[colCount];
588: }
589: }
590: if (dest == null) {
591: dest = new double[bandUp];
592: }
593: /*
594: * Builds up a RectIter and use it for interpolating all bands.
595: * There is very few points, so the cost of creating a RectIter
596: * may be important. But it seems to still lower than query tiles
597: * many time (which may involve more computation than necessary).
598: */
599: bounds.x = ix - left;
600: bounds.y = iy - top;
601: final RectIter iter = RectIterFactory.create(image, bounds);
602: for (; band < bandUp; band++) {
603: iter.startLines();
604: int j = 0;
605: do {
606: iter.startPixels();
607: final double[] row = samples[j++];
608: int i = 0;
609: do {
610: row[i++] = iter.getSampleDouble(band);
611: } while (!iter.nextPixelDone());
612: assert i == row.length;
613: } while (!iter.nextLineDone());
614: assert j == samples.length;
615: float dx = (float) (x - x0);
616: if (dx == 1)
617: dx = ONE_EPSILON;
618: float dy = (float) (y - y0);
619: if (dy == 1)
620: dy = ONE_EPSILON;
621: final double value = interpolation.interpolate(samples, dx,
622: dy);
623: if (Double.isNaN(value)) {
624: if (fallback == this )
625: continue; // 'dest' was set by 'super.evaluate(...)'.
626: if (fallback != null) {
627: fallback.interpolate(x, y, dest, band, band + 1);
628: continue;
629: }
630: // If no fallback was specified, then 'dest' is not required to
631: // have been initialized. It may contains random value. Set it
632: // to the NaN value...
633: }
634: dest[band] = value;
635: }
636: return dest;
637: }
638: }
|