Hi Jody et al.
Don't know if the code below is useful in this context, but if it is
feel free to use it (suitably hacked). It's a routine to fit a smooth
curve through a set of points using cubic spline interpolation. Works
for me but you might want to test it out.
ciao
Michael
/**
* A utility class for cubic spline interpolation in two dimensions.
*
* @author Michael Bedward <[EMAIL PROTECTED]>
*/
public class SplineInterpolator {
/**
* Cubic spline interpolation through the set of points with
* coordinates given by x and y input arrays.
* <P>
* Returns the coords of a line interpolated through np points
* (including the first and last control points) which are
* approximately evenly-spaced (fudged using linear interp).
* <P>
* Adapted from an algorithm in Press et al. 1992
* "Numeical Recipes in C". Cambridge University Press
*
* @param x x-coords of control points
* @param y y-coords of control points
* @param np number of interpolation points
* @return a 2D double array of interpolated point locations with
* x-xoords in array[0][] and y-coords in array[1][]
*
* @throws java.lang.IllegalArgumentException
*/
public static double[][] interp(double[] x, double[] y, int np)
throws IllegalArgumentException {
int N = x.length;
if (N == 0 || N != y.length) {
throw new IllegalArgumentException("length of input array
must be equal and > 0");
}
if (np < N) {
throw new IllegalArgumentException("np should be larger
than length of input arrays");
}
// x values need to be considered in ascending order
int[] o = getOrder(x);
// Interpolation (as done below) will fail if any of the input
x coords are equal
double prec = 1.0e-8;
for (int i = 1; i < N; i++) {
if (Math.abs(x[o[i]] - x[o[i - 1]]) < prec) {
throw new IllegalArgumentException("x values must be distinct");
}
}
// first derivatives
double[] dy = new double[N];
// second derivatives
double[] ddy = new double[N];
// in this implementation we set the value of the second
derivatives at the end-points
// to 0 (the so-called natural cublic spline). See Press et
al. 1992 for alternative
// approaches.
ddy[0] = ddy[N - 1] = dy[0] = 0.0;
//
========================================================================
// tridiagonal algorithm to calculate the piece-wise splines
double sig, p;
for (int i = 1; i < N - 1; i++) {
sig = (x[o[i]] - x[o[i - 1]]) / (x[o[i + 1]] - x[o[i - 1]]);
p = sig * ddy[i - 1] + 2.0;
ddy[i] = (sig - 1.0) / p;
dy[i] = (y[o[i + 1]] - y[o[i]]) / (x[o[i + 1]] - x[o[i]])
- (y[o[i]] - y[o[i - 1]]) / (x[o[i]] - x[o[i - 1]]);
dy[i] = (6.0d * dy[i] / (x[o[i + 1]] - x[o[i - 1]]) - sig
* dy[i - 1]) / p;
}
for (int i = N - 2; i >= 0; i--) {
ddy[i] = ddy[i] * ddy[i + 1] + dy[i];
}
//
========================================================================
// interpolation
// first find approximate distance along curve by linear interp
double[] dist = new double[N];
double[] xdiff = new double[N];
dist[0] = xdiff[0] = 0.0;
for (int i = 1; i < N; i++) {
xdiff[i] = x[o[i]] - x[o[i - 1]];
double ydiff = y[o[i]] - y[o[i - 1]];
dist[i] = dist[i - 1] + Math.sqrt(xdiff[i] * xdiff[i] +
ydiff * ydiff);
}
// spacing
double dIncr = dist[N - 1] / (np - 1);
double[][] xyInterp = new double[2][np];
xyInterp[0][0] = x[o[0]];
xyInterp[1][0] = y[o[0]];
xyInterp[0][np - 1] = x[o[N - 1]];
xyInterp[1][np - 1] = y[o[N - 1]];
double curDist = dIncr;
int k = 1;
for (int i = 1; i < np - 1; i++, curDist += dIncr) {
while (curDist > dist[k]) {
k++;
}
// locate the x coord of the interpolated point
double delta = (curDist - dist[k - 1]) / (dist[k] - dist[k - 1]);
xyInterp[0][i] = x[o[k - 1]] + delta * xdiff[k];
// locate the y coord by evaluating the cubic spline polynomial
double h = x[o[k]] - x[o[k - 1]];
double a = (x[o[k]] - xyInterp[0][i]) / h;
double b = (xyInterp[0][i] - x[o[k - 1]]) / h;
xyInterp[1][i] = a * y[o[k - 1]] + b * y[o[k]] +
((a * a * a - a) * ddy[k - 1] + (b * b * b - b) *
ddy[k]) * (h * h) / 6.0;
}
return xyInterp;
}
/*
* Get order of values in the input array.
*/
private static int[] getOrder(double[] x) {
if (x.length == 0) {
return null;
}
int[] order = new int[x.length];
for (int i = 0; i < x.length; i++) {
order[i] = i;
}
// using shellsort algorithm as modified from Sedgewick 1990...
// Algorithms in C
//
int[] incs = new int[]{
1391376, 463792, 198768, 86961, 33936, 13776,
4592, 1968, 861, 336, 112, 48, 21, 7, 3, 1
};
int i0, j;
double v;
for (int k = 0; k < incs.length; k++) {
for (int h = incs[k], i = h; i < x.length; i++) {
j = i;
i0 = order[j];
v = x[i0];
while (j > h - 1 && x[order[j - h]] > v) {
order[j] = order[j - h];
j -= h;
}
order[j] = i0;
}
}
return order;
}
}
-------------------------------------------------------------------------
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/
_______________________________________________
Geotools-gt2-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
