edu.mines.jtk.interp.TrilinearInterpolator3 Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2013, 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.dsp.Sampling;
import edu.mines.jtk.util.Check;
import static edu.mines.jtk.util.ArrayMath.*;
/**
* Piecewise trilinear polynomial interpolation of a function y(x1,x2,x3).
*
* The function y(x1,x2,x3) is specified by samples on a regular grid, which
* need not be uniform. The regular grid is specified by one-dimensional
* arrays of monotonically increasing coordinates x1, x2 and x3, such that
* gridded x1 are identical for all gridded x2 and x3, gridded x2 are
* identical for all gridded x1 and x2, and gridded x3 are identical for all
* gridded x1 and x2.
*
* Extrapolation (that is, interpolation outside the specified grid of
* (x1,x2,x3) coordinates), is performed using linear polynomials for the
* nearest grid samples. Extrapolated values can be well outside the [min,max]
* range of interpolated values, and should typically be avoided.
*
* @author Dave Hale, Colorado School of Mines
* @version 2013.01.22
*/
public class TrilinearInterpolator3 {
/**
* Constructs an interpolator for specified values y(x1,x2).
* @param x1 array[n1] of x1 coordinates; monotonically increasing.
* @param x2 array[n2] of x2 coordinates; monotonically increasing.
* @param x3 array[n3] of x3 coordinates; monotonically increasing.
* @param y array[n3][n2][n1] of sampled values y(x1,x2,x3).
*/
public TrilinearInterpolator3(
float[] x1, float[] x2, float[] x3, float[][][] y)
{
this(x1.length,x2.length,x3.length,x1,x2,x3,y);
}
/**
* Constructs an interpolator for specified values.
* @param n1 number of x1 coordinates.
* @param n2 number of x2 coordinates.
* @param n3 number of x3 coordinates.
* @param x1 array[n1] of x1 coordinates; monotonically increasing.
* @param x2 array[n2] of x2 coordinates; monotonically increasing.
* @param x3 array[n3] of x3 coordinates; monotonically increasing.
* @param y array[n3][n2][n1] of sampled values y(x1,x2,x3).
*/
public TrilinearInterpolator3(
int n1, int n2, int n3,
float[] x1, float[] x2, float[] x3, float[][][] y)
{
Check.argument(isMonotonic(x1), "array x1 is monotonic");
Check.argument(isMonotonic(x2), "array x2 is monotonic");
Check.argument(isMonotonic(x3), "array x3 is monotonic");
_x1 = copy(n1,x1);
_x2 = copy(n2,x2);
_x3 = copy(n3,x3);
makeCoefficients(x1,x2,x3,y);
}
/**
* Returns interpolated value y.
* Same as {@link #interpolate000(float,float,float)}.
* @param x1 coordinate x1.
* @param x2 coordinate x2.
* @param x3 coordinate x3.
* @return interpolated y(x1,x2,x3).
*/
public float interpolate(float x1, float x2, float x3) {
return interpolate000(x1,x2,x3);
}
/**
* Returns interpolated value y.
* @param x1 coordinate x1.
* @param x2 coordinate x2.
* @param x3 coordinate x3.
* @return interpolated y(x1,x2,x3).
*/
public float interpolate000(float x1, float x2, float x3) {
return interpolate000(x1,x2,x3,_ks);
}
/**
* Returns interpolated partial derivative dy/d1.
* @param x1 coordinate x1.
* @param x2 coordinate x2.
* @param x3 coordinate x3.
* @return interpolated dy/d1.
*/
public float interpolate100(float x1, float x2, float x3) {
return interpolate100(x1,x2,x3,_ks);
}
/**
* Returns interpolated partial derivative dy/d2.
* @param x1 coordinate x1.
* @param x2 coordinate x2.
* @param x3 coordinate x3.
* @return interpolated dy/d2.
*/
public float interpolate010(float x1, float x2, float x3) {
return interpolate010(x1,x2,x3,_ks);
}
/**
* Returns interpolated partial derivative dy/d3.
* @param x1 coordinate x1.
* @param x2 coordinate x2.
* @param x3 coordinate x3.
* @return interpolated dy/d3.
*/
public float interpolate001(float x1, float x2, float x3) {
return interpolate001(x1,x2,x3,_ks);
}
/**
* Returns an array of interpolated values y.
* Same as {@link #interpolate000(Sampling,Sampling,Sampling)}.
* @param s1 sampling of coordinate x1.
* @param s2 sampling of coordinate x2.
* @param s3 sampling of coordinate x3.
* @return interpolated y(x1,x2,x3).
*/
public float[][][] interpolate(Sampling s1, Sampling s2, Sampling s3) {
return interpolate000(s1,s2,s3);
}
/**
* Returns an array of interpolated values y.
* @param s1 sampling of coordinate x1.
* @param s2 sampling of coordinate x2.
* @param s3 sampling of coordinate x3.
* @return interpolated y(x1,x2,x3).
*/
public float[][][] interpolate000(Sampling s1, Sampling s2, Sampling s3) {
int n1 = s1.getCount();
int n2 = s2.getCount();
int n3 = s3.getCount();
float[][][] y = new float[n3][n2][n1];
interpolate000(s1,s2,s3,y);
return y;
}
/**
* Computes an array of interpolated values y.
* Same as {@link #interpolate000(Sampling,Sampling,Sampling,float[][][])}.
* @param s1 sampling of coordinate x1.
* @param s2 sampling of coordinate x2.
* @param s3 sampling of coordinate x3.
* @param y output array of interpolated y(x1,x2,x3).
*/
public void interpolate(
Sampling s1, Sampling s2, Sampling s3, float[][][] y)
{
interpolate000(s1,s2,s3,y);
}
/**
* Computes an array of interpolated values y.
* @param s1 sampling of coordinate x1.
* @param s2 sampling of coordinate x2.
* @param s3 sampling of coordinate x3.
* @param y output array of interpolated y(x1,x2,x3).
*/
public void interpolate000(
Sampling s1, Sampling s2, Sampling s3, float[][][] y)
{
int n1 = s1.getCount();
int n2 = s2.getCount();
int n3 = s3.getCount();
int[] k1 = makeIndices(s1,_x1);
int[] k2 = makeIndices(s2,_x2);
int[] k3 = makeIndices(s3,_x3);
for (int i3=0; i3=xs.length-1)
i = xs.length-2;
return i;
}
private void updateIndices(float x1, float x2, float x3, int[] ks) {
ks[0] = index(x1,_x1,ks[0]);
ks[1] = index(x2,_x2,ks[1]);
ks[2] = index(x3,_x3,ks[2]);
}
private static int[] makeIndices(float[] xi, float[] xs) {
int n = xi.length;
int[] ki = new int[n];
ki[0] = index(xi[0],xs,0);
for (int i=1; i