/*
* GeoTools - The Open Source Java GIS Toolkit
* http://geotools.org
*
* (C) 2011, Open Source Geospatial Foundation (OSGeo)
* (C) 2008-2011 TOPP - www.openplans.org.
*
* 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.process.vector;
/**
* Interpolates a grid to a grid of different dimensions using bilinear interpolation.
* <p>
* NO_DATA cell values are supported in the source grid.
* There are two ways of handling the boundary between NO_DATA cells and data cells:
* <ul>
* <li><b>Truncate</b> - If any source cell is NO_DATA, the dest value is NO_DATA.
* This is simple and fast, but does make the data boundaries look a bit ragged.
* <li><b>Smooth</b> - If only one source cell is NO_DATA, the value is interpolated
* using the 3 valid values, across one-half of the interpolated cells.
* This smoothes off the boundary.
* If 2 or more source cells are NO_DATA, then Truncation is used.
* </ul>
* <p>
* Reference: http://en.wikipedia.org/wiki/Bilinear_interpolation.
*
* @author Martin Davis - OpenGeo
*
*/
public class BilinearInterpolator {
private static final float NULL_NO_DATA = Float.NaN;
private final float[][] src;
private float noDataValue = NULL_NO_DATA;
/**
* Creates a new interpolator on a given source grid.
*
* @param src the source grid
*/
public BilinearInterpolator(final float[][] src)
{
this(src, NULL_NO_DATA);
}
/**
* Creates a new interpolator on a given source grid.
*
* @param src the source grid
* @param noDataValue the NO_DATA value (if none use Float.NaN)
*/
public BilinearInterpolator(final float[][] src, final float noDataValue)
{
this.src = src;
this.noDataValue = noDataValue;
}
/**
* Interpolates the source grid to a new grid of different dimensions.
*
* @param width the width of the destination grid
* @param height the height of the destination grid
* @param smoothBoundary true if boundary smoothing should be performed
* @return the interpolated grid
*/
public float[][] interpolate(final int width, final int height, boolean smoothBoundary)
{
int srcWidth = src.length;
int srcHeight = src[0].length;
float[][] dest = new float[width][height];
float xRatio = ((float) srcWidth - 1) / width ;
float yRatio = ((float) srcHeight - 1) / height ;
for (int i = 0; i < width; i++) {
for (int j = 0; j < height; j++) {
float x = i * xRatio;
float y = j * yRatio;
int ix = (int) x;
int iy = (int) y;
float xfrac = x - ix;
float yfrac = y - iy;
float val;
if (ix < srcWidth - 1 && iy < srcHeight - 1) {
// interpolate if src cell is in grid
float v00 = src[ix][iy];
float v10 = src[ix + 1][iy];
float v01 = src[ix][iy + 1];
float v11 = src[ix + 1][iy + 1];
if (v00 == noDataValue
|| v10 == noDataValue
|| v01 == noDataValue
|| v11 == noDataValue) {
// handle src cell with NO_DATA value(s)
if (smoothBoundary) {
val = interpolateBoundaryCell(xfrac, yfrac, v00, v10, v01, v11);
}
else {
val = noDataValue;
}
}
else {
// All src cell corners have values
// Compute bilinear interpolation over the src cell
val = ( v00*(1-xfrac)*(1-yfrac) + v10*(xfrac)*(1-yfrac) +
v01*(yfrac)*(1-xfrac) + v11*(xfrac*yfrac)
) ;
}
}
else {
// dest index at edge of grid
val = src[ix][iy];
}
dest[i][j] = val;
}
}
return dest;
}
/**
* Interpolates across an edge grid cell, which has 1 or more NO_DATA values.
* Grid cells with 2 or or NO_DATA values still receive the value NO_DATA.
* Otherwise, the cell is interpolated across the triangle defined by the
* 3 valid corner values.
* This produces a much smoother edge appearance, containing 45-degree lines
* instead of a jagged stairstep boundary.
*
* @param xfrac the fractional x location of the interpolation point within the cell
* @param yfrac the fractional y location of the interpolation point within the cell
* @param v00 the lower left value
* @param v10 the lower right value
* @param v01 the upper left value
* @param v11 the upper right value
* @return the interpolated value
*/
private float interpolateBoundaryCell(float xfrac, float yfrac, float v00, float v10, float v01, float v11)
{
// count noData values
int count = 0;
if (v00 == noDataValue) count++;
if (v10 == noDataValue) count++;
if (v01 == noDataValue) count++;
if (v11 == noDataValue) count++;
/**
* Cells with >= 2 noData ==> noData
*/
if (count > 1) return noDataValue;
/**
* Now only one cell has noData value.
* Compute interpolation over cell, with vertex layout normalized to put noData in NE.
* This is done by flipping the cell across the X or Y axis, or both
* (and transforming the point offsets likewise)
*/
if (v00 == noDataValue) return interpolateBoundaryCellNorm(1.0f - yfrac, 1.0f - xfrac, v11, v10, v01);
if (v11 == noDataValue) return interpolateBoundaryCellNorm(xfrac, yfrac, v00, v10, v01);
if (v10 == noDataValue) return interpolateBoundaryCellNorm(xfrac, 1.0f - yfrac, v01, v11, v00);
if (v01 == noDataValue) return interpolateBoundaryCellNorm(1.0f - xfrac, yfrac, v10, v00, v11);
// should never reach here
return noDataValue;
}
/**
* Computes an interpolated value across a grid cell which has a single
* NO_DATA value in the NE corner (<tt>v11</tt>), and valid data values
* in the three other corners.
* The interpolated value is computed using linear interpolation across
* the 3D triangle defined by the three valid corners.
*
* @param xfrac the fractional x location of the interpolation point within the cell
* @param yfrac the fractional y location of the interpolation point within the cell
* @param v00 the lower left value
* @param v10 the lower right value
* @param v01 the upper left value
* @return the interpolated value
*/
private float interpolateBoundaryCellNorm(float xfrac, float yfrac, float v00, float v10, float v01)
{
// if point is in NE triangle, value is NO_DATA
if (xfrac + yfrac > 1) return noDataValue;
// interpolate across plane defined by SW triangle and values
float dx = v10 - v00;
float dy = v01 - v00;
float v = v00 + (xfrac * dx) + (yfrac * dy);
return v;
}
}