edu.mines.jtk.interp.CubicInterpolator Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2006, Colorado School of Mines and others.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
****************************************************************************/
package edu.mines.jtk.interp;
import edu.mines.jtk.util.Check;
import static edu.mines.jtk.util.ArrayMath.*;
/**
* Piecewise cubic interpolation of a function y(x).
*
* Piecewise cubic interpolators differ in the method they use to compute
* slopes y'(x) at specified x (knots). The classic cubic spline computes the
* slopes to obtain a continuous second derivative at the knots. These splines
* often yield unacceptable wiggliness (overshoot) between the knots. A linear
* spline yields no overshoot, but has discontinuous first (and higher)
* derivatives. A monotonic spline has continuous first derivatives and yields
* monotonic interpolation (with no overshoot) where function values at the
* knots are monotonically increasing or decreasing.
*
* For x outside the range of values specified when an interpolator was
* constructed, the interpolator extrapolates using the cubic
* polynomial corresponding to the knot nearest to x.
* @author Dave Hale, Colorado School of Mines
* @version 2012.12.27
*/
public class CubicInterpolator {
/**
* The method used to compute 1st derivatives y'(x).
*/
public enum Method {
/**
* The interpolated y(x) is continuous, but has discontinuous 1st and
* higher derivatives. This method is equivalent to (though less efficient
* than) simple piecewise linear interpolation.
*/
LINEAR,
/**
* The interpolated y(x) is continuous with continuous 1st derivative, but
* may have discontinuous 2nd and higher-order derivatives. This method
* preserves monotonicity. In intervals where specified y(x) are
* monotonic, the interpolated values y(x) are also monotonic.
*/
MONOTONIC,
/**
* The interpolated y(x) is continuous with continuous 1st and 2nd
* derivatives, but may have discontinuous 3rd and higher-order
* derivatives.
*/
SPLINE
}
/**
* Constructs an interpolator with specified 1st derivatives y'(x).
* @param x array of values at which y(x) are specified.
* These values must be monotonically increasing or decreasing,
* with no equal values. (In other words, the array must be
* monotonic-definite.)
* @param y array of function values y(x).
* @param y1 array of 1st derivatives y'(x).
*/
public CubicInterpolator(float[] x, float[] y, float[] y1) {
Check.argument(isMonotonic(x), "array x is monotonic");
int n = x.length;
_xd = new float[n];
_yd = new float[n][4];
for (int i=0; i
* The Fritsch-Carlson method yields continuous 1st derivatives, but 2nd
* and 3rd derivatives are discontinuous. The method will yield a
* monotonic interpolant for monotonic data. 1st derivatives are set to
* zero wherever first divided differences change sign.
*
* For more information, see Fritsch, F. N., and Carlson, R. E., 1980,
* Monotone piecewise cubic interpolation: SIAM J. Numer. Anal., v. 17,
* n. 2, p. 238-246.
*
* Also, see the book by Kahaner, D., Moler, C., and Nash, S., 1989,
* Numerical Methods and Software, Prentice Hall. This function was
* derived from SUBROUTINE PCHEZ contained on the diskette that comes
* with the book.
*/
private static void initMonotonic(int n, float[] x, float[][] y) {
// If n=1, then use constant interpolation.
if (n==1) {
y[0][1] = 0.0f;
y[0][2] = 0.0f;
y[0][3] = 0.0f;
return;
// Else, if n=2, then use linear interpolation.
} else if (n==2) {
y[0][1] = y[1][1] = (y[1][0]-y[0][0])/(x[1]-x[0]);
y[0][2] = y[1][2] = 0.0f;
y[0][3] = y[1][3] = 0.0f;
return;
}
// Set left end derivative via shape-preserving 3-point formula.
float h1 = x[1]-x[0];
float h2 = x[2]-x[1];
float hsum = h1+h2;
float del1 = (y[1][0]-y[0][0])/h1;
float del2 = (y[2][0]-y[1][0])/h2;
float w1 = (h1+hsum)/hsum;
float w2 = -h1/hsum;
y[0][1] = w1*del1+w2*del2;
if (y[0][1]*del1<=0.0f)
y[0][1] = 0.0f;
else if (del1*del2<0.0f) {
float dmax = 3.0f*del1;
if (abs(y[0][1])>abs(dmax)) y[0][1] = dmax;
}
// Loop over interior points.
for (int i=1; iabs(dmax)) y[n-1][1] = dmax;
}
compute2ndAnd3rdDerivatives(x,y);
}
/**
* Computes cubic spline interpolation coefficients for interpolation
* with continuous second derivatives.
*/
private static void initSpline(int n, float[] x, float[][] y) {
// If n=1, then use constant interpolation.
if (n==1) {
y[0][1] = 0.0f;
y[0][2] = 0.0f;
y[0][3] = 0.0f;
return;
// Else, if n=2, then use linear interpolation.
} else if (n==2) {
y[0][1] = y[1][1] = (y[1][0]-y[0][0])/(x[1]-x[0]);
y[0][2] = y[1][2] = 0.0f;
y[0][3] = y[1][3] = 0.0f;
return;
}
// Set left end derivative via shape-preserving 3-point formula.
float h1 = x[1]-x[0];
float h2 = x[2]-x[1];
float hsum = h1+h2;
float del1 = (y[1][0]-y[0][0])/h1;
float del2 = (y[2][0]-y[1][0])/h2;
float w1 = (h1+hsum)/hsum;
float w2 = -h1/hsum;
float sleft = w1*del1+w2*del2;
if (sleft*del1<=0.0f)
sleft = 0.0f;
else if (del1*del2<0.0f) {
float dmax = 3.0f*del1;
if (abs(sleft)>abs(dmax)) sleft = dmax;
}
// Set right end derivative via shape-preserving 3-point formula.
h1 = x[n-2]-x[n-3];
h2 = x[n-1]-x[n-2];
hsum = h1+h2;
del1 = (y[n-2][0]-y[n-3][0])/h1;
del2 = (y[n-1][0]-y[n-2][0])/h2;
w1 = -h2/hsum;
w2 = (h2+hsum)/hsum;
float sright = w1*del1+w2*del2;
if (sright*del2<=0.0f)
sright = 0.0f;
else if (del1*del2<0.0f) {
float dmax = 3.0f*del2;
if (abs(sright)>abs(dmax)) sright = dmax;
}
// Compute tridiagonal system coefficients and right-hand-side.
float[] work = new float[n];
work[0] = 1.0f;
y[0][2] = 2.0f*sleft;
for (int i=1; i=0; --i)
y[i][1] -= y[i+1][3]*y[i+1][1];
compute2ndAnd3rdDerivatives(x,y);
}
private static void compute2ndAnd3rdDerivatives(float[] x, float[][] y) {
int n = x.length;
for (int i=0; i