edu.mines.jtk.interp.SibsonInterpolator2 Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2009, 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 java.util.ArrayList;
import edu.mines.jtk.dsp.Sampling;
import edu.mines.jtk.la.DMatrix;
import edu.mines.jtk.la.DMatrixQrd;
import edu.mines.jtk.mesh.Geometry;
import edu.mines.jtk.mesh.TriMesh;
import edu.mines.jtk.util.Check;
/**
* Sibson interpolation of scattered samples of 2D functions f(x1,x2).
* Sibson's (1981) interpolant at any point (x1,x2) is a weighted sum of
* values for a nearby subset of samples, the so-called natural neighbors
* of that point. Sibson interpolation is often called "natural neighbor
* interpolation."
*
* The interpolation weights, also called "Sibson coordinates", are areas
* of overlapping Voronoi polygons, normalized so that they sum to one for
* any interpolation point (x1,x2). Various implementations of Sibson
* interpolation differ primarily in how those areas are computed.
*
* The basic Sibson interpolant is C1 (that is, it's gradient is continuous)
* at all points (x1,x2) except at the sample points, where it is C0.
* Sibson (1981) also described an extension of his interpolant that is
* everywhere C1 and therefore smoother. This smoother interpolant requires
* gradients at the sample points, and those gradients can be estimated or
* specified explicitly.
*
* The use of gradients is controlled by a gradient power. If this power
* is zero (the default), then gradients are not used. Sibson's (1981)
* smoother C1 interpolant corresponds to a power of 1.0. Larger powers
* cause the interpolant to more rapidly approach the linear functions
* defined by the values and gradients at the sample points.
*
* Sibson's interpolant is undefined at points on or outside the convex
* hull of sample points. In this sense, Sibson interpolation does not
* extrapolate; the interpolant is implicitly bounded by the convex hull,
* and null values are returned when attempting to interpolate outside
* those bounds.
*
* To extend the interpolant outside the convex hull, this class enables
* bounds to be set explicitly. When bounds are set, extra ghost samples
* are added outside the convex hull. These ghost samples have no values,
* but they create a larger convex hull so that Sibson coordinates can be
* computed anywhere within the specified bounds. While often useful, this
* extrapolation feature should be used with care, because the added ghost
* samples may alter the Sibson interpolant at points inside but near the
* original convex hull.
*
* References:
*
-
* Braun, J. and M. Sambridge, 1995, A numerical method for solving partial
* differential equations on highly irregular evolving grids: Nature, 376,
* 655--660.
*
-
* Lasserre J.B., 1983, An analytical expression and an algorithm for the
* volume of a convex polyhedron in R^n: Journal of Optimization Theory
* and Applications, 39, 363--377.
*
-
* Sambridge, M., J. Braun, and H. McQueen, 1995, Geophysical
* parameterization and interpolation of irregular data using
* natural neighbors.
*
-
* Sibson, R., 1981, A brief description of natural neighbor interpolation,
* in V. Barnett, ed., Interpreting Multivariate Data, John Wiley and Sons,
* 21--36.
*
-
* Watson, D.F., 1992, Contouring: a guide to the analysis and display
* of spatial data: Pergamon, Oxford.
*
*
* @author Dave Hale, Colorado School of Mines
* @version 2009.06.14
*/
public class SibsonInterpolator2 {
/**
* The implementation method. Methods differ in the algorithm by which
* Sibson coordinates (polygon areas) are computed for natural neighbors.
*/
public enum Method {
/**
* The Hale-Liang algorithm.
* Developed by Dave Hale and Luming Liang at the Colorado School of
* Mines. Accurate, robust and fast, this method is the default.
*/
HALE_LIANG,
/**
* The Braun-Sambridge algorithm. Developed by Braun and Sambridge (1995),
* this method uses Lasserre's (1983) algorithm for computing the areas
* of polygons without computing their vertices. This method is much
* slower than the other methods provided here. It may also suffer from
* numerical instability due to divisions in Lasserre's algorithm.
*/
BRAUN_SAMBRIDGE,
/**
* The Watson-Sambridge algorithm. Developed by Watson (1992) and
* described further by Sambridge et al. (1995). Though simplest to
* implement and fast, this method is inaccurate (sometimes wildly so)
* at interpolation points (x1,x2) that lie on or near edges of a
* Delaunay triangulation of the scattered sample points.
*/
WATSON_SAMBRIDGE,
}
/**
* Sample index and corresponding interpolation weight (Sibson coordinate).
*/
public static class IndexWeight {
/** Index for one sample. */
public int index;
/** Interpolation weight for one sample; in the range [0,1].*/
public float weight;
IndexWeight(int index, float weight) {
this.index = index;
this.weight = weight;
}
}
/**
* Constructs an interpolator with specified sample coordinates.
* Function values f(x1,x2) are not set and are assumed to be zero.
* Uses the most accurate and efficient implementation.
* Coordinates for each sample must be unique.
* @param x1 array of sample x1 coordinates.
* @param x2 array of sample x2 coordinates.
*/
public SibsonInterpolator2(float[] x1, float[] x2) {
this(Method.HALE_LIANG,null,x1,x2);
}
/**
* Constructs an interpolator with specified samples.
* Uses the most accurate and efficient implementation.
* Coordinates for each sample must be unique.
* @param f array of sample values f(x1,x2).
* @param x1 array of sample x1 coordinates.
* @param x2 array of sample x2 coordinates.
*/
public SibsonInterpolator2(float[] f, float[] x1, float[] x2) {
this(Method.HALE_LIANG,f,x1,x2);
}
/**
* Constructs an interpolator with specified method and sample coordinates.
* Function values f(x1,x2) are not set and are assumed to be zero.
* Coordinates for each sample must be unique.
*
* This constructor is provided primarily for testing.
* The default Hale-Liang method is accurate and fast.
* @param method the implementation method.
* @param x1 array of sample x1 coordinates.
* @param x2 array of sample x2 coordinates.
*/
public SibsonInterpolator2(Method method, float[] x1, float[] x2) {
this(method,null,x1,x2);
}
/**
* Constructs an interpolator with specified method and samples.
* Coordinates for each sample must be unique.
*
* This constructor is provided primarily for testing.
* The default Hale-Liang method is accurate and fast.
* @param method the implementation method.
* @param f array of sample values f(x1,x2).
* @param x1 array of sample x1 coordinates.
* @param x2 array of sample x2 coordinates.
*/
public SibsonInterpolator2(
Method method, float[] f, float[] x1, float[] x2)
{
makeMesh(f,x1,x2);
_nodeList = new TriMesh.NodeList();
_triList = new TriMesh.TriList();
if (method==Method.WATSON_SAMBRIDGE) {
_va = new WatsonSambridge();
} else if (method==Method.BRAUN_SAMBRIDGE) {
_va = new BraunSambridge();
} else if (method==Method.HALE_LIANG) {
_va = new HaleLiang();
}
}
/**
* Sets the samples to be interpolated.
* Any sample coordinates, values or gradients set previously are forgotten.
* @param f array of sample values f(x1,x2).
* @param x1 array of sample x1 coordinates.
* @param x2 array of sample x2 coordinates.
*/
public void setSamples(float[] f, float[] x1, float[] x2) {
makeMesh(f,x1,x2);
_haveGradients = false;
if (_gradientPower>0.0)
estimateGradients();
}
/**
* Sets the null value for this interpolator.
* This null value is returned when interpolation is attempted at a
* point that lies outside the bounds of this interpolator. Those
* bounds are by default the convex hull of the sample points, but
* may also be set explicitly. The default null value is zero.
* @param fnull the null value.
*/
public void setNullValue(float fnull) {
_fnull = fnull;
}
/**
* Sets gradients for all samples. If the gradient power is currently
* zero, then this method also sets the gradient power to one. To later
* ignore gradients that have been set, the gradient power can be reset
* to zero.
* @param g1 array of 1st components of gradients.
* @param g2 array of 2nd components of gradients.
*/
public void setGradients(float[] g1, float[] g2) {
int n = g1.length;
for (int i=0; i
* If the gradient power is set to a non-zero value, and if gradients
* have not been set explicitly, then this method will also estimate
* gradients for all samples, as described by Sibson (1981).
*
* If bounds are set explicitly, gradient estimates can be improved
* by setting the bounds before calling this method.
* @param gradientPower the gradient power.
*/
public void setGradientPower(double gradientPower) {
if (!_haveGradients && gradientPower>0.0)
estimateGradients();
_gradientPower = gradientPower;
}
/**
* Sets a bounding box for this interpolator.
* Sibson interpolation is undefined for points (x1,x2) outside the
* convex hull of sample points, so the default bounds are that convex
* hull. This method enables extrapolation for points outside the convex
* hull, while restricting interpolation to points inside the box.
*
* If gradients are to be computed (not specified explicitly), it is best
* to set bounds by calling this method before computing gradients.
* @param x1min lower bound on x1.
* @param x1max upper bound on x1.
* @param x2min lower bound on x2.
* @param x2max upper bound on x2.
*/
public void setBounds(float x1min, float x1max, float x2min, float x2max) {
Check.argument(x1min
* If gradients are to be computed (not specified explicitly), it is best
* to set bounds by calling this method before computing gradients.
* @param s1 sampling of x1.
* @param s2 sampling of x2.
*/
public void setBounds(Sampling s1, Sampling s2) {
setBounds((float)s1.getFirst(),(float)s1.getLast(),
(float)s2.getFirst(),(float)s2.getLast());
}
/**
* If bounds have been set explicitly, this method unsets them,
* so that the convex hull of sample points will be used instead.
*/
public void useConvexHullBounds() {
_useBoundingBox = false;
removeGhostNodes();
}
/**
* Returns a value interpolated at the specified point.
* @param x1 the x1 coordinate of the point.
* @param x2 the x2 coordinate of the point.
* @return the interpolated value.
*/
public float interpolate(float x1, float x2) {
if (!inBounds(x1,x2))
return _fnull;
double asum = computeAreas(x1,x2);
if (asum<=0.0)
return _fnull;
if (usingGradients()) {
return interpolate1(asum,x1,x2);
} else {
return interpolate0(asum);
}
}
/**
* Returns an array of interpolated values sampled on a grid.
* @param s1 the sampling of n1 x1 coordinates.
* @param s2 the sampling of n2 x2 coordinates.
* @return array[n2][n1] of interpolated values.
*/
public float[][] interpolate(Sampling s1, Sampling s2) {
int n1 = s1.getCount();
int n2 = s2.getCount();
float[][] f = new float[n2][n1];
for (int i2=0; i2
* Indices and weights are especially useful in applications where they
* can be reused, say, to interpolate multiple function values at a
* single point.
* @param x1 the x1 coordinate of the point.
* @param x2 the x2 coordinate of the point.
* @return array of sample indices and weights; null if none.
*/
public IndexWeight[] getIndexWeights(float x1, float x2) {
if (!inBounds(x1,x2))
return null;
float wsum = (float)computeAreas(x1,x2);
if (wsum==0.0f)
return null;
float wscl = 1.0f/wsum;
int nnode = _nodeList.nnode();
TriMesh.Node[] nodes = _nodeList.nodes();
IndexWeight[] iw = new IndexWeight[nnode];
for (int inode=0; inode
* If bounds have not been set explicitly, then this method will return
* a null value if the validated sample is on the convex hull of samples.
*
* This method does not recompute gradients that may have been estimated
* using the sample to be validated. Therefore, validation should be
* performed without using gradients.
* @param i the index of the sample to validate.
* @return the value interpolated at the validated sample point.
*/
public float validate(int i) {
return validate(new int[]{i})[0];
}
/**
* Interpolates at specified sample points without using those samples.
* This method implements a form of cross-validation. Differences
* between the values of the specified samples and the returned
* interpolated values are measures of errors for those samples.
*
* If bounds have not been set explicitly, then this method will return
* null values if the validated sample is on the convex hull of samples.
*
* This method does not recompute gradients that may have been estimated
* using the samples to be validated. Therefore, validation should be
* performed without using gradients.
* @param i array of indices of samples to validate.
* @return array of values interpolated at validated sample points.
*/
public float[] validate(int[] i) {
int nv = i.length;
for (int iv=0; iv_x1max) _x1max = x1[i];
if (x2[i]<_x2min) _x2min = x2[i];
if (x2[i]>_x2max) _x2max = x2[i];
}
// ArrayMath of nodes, one for each specified sample.
_nodes = new TriMesh.Node[n];
// Construct the mesh with nodes at sample points.
_mesh = new TriMesh();
for (int i=0; i0.0;
}
// Adds ghost nodes to the mesh without any values or gradients.
// The arrays x[ng] and y[ng] contain coordinates of ng ghost nodes.
private void addGhostNodes(float[] x, float[] y) {
int ng = x.length;
for (int ig=0; ig gnodes = new ArrayList(4);
TriMesh.NodeIterator ni = _mesh.getNodes();
while (ni.hasNext()) {
TriMesh.Node n = ni.next();
if (ghost(n))
gnodes.add(n);
}
// Then remove the ghost nodes from the mesh.
for (TriMesh.Node gnode:gnodes)
_mesh.removeNode(gnode);
}
// Computes Sibson areas for the specified point (x,y).
// Returns true, if successful; false, otherwise.
private double computeAreas(float x, float y) {
if (!getNaturalNabors(x,y))
return 0.0;
return _va.accumulateAreas(x,y,_mesh,_nodeList,_triList);
}
// Returns true if not using bounding box or if point is inside the box.
private boolean inBounds(float x1, float x2) {
return !_useBoundingBox ||
_x1bmn<=x1 && x1<=_x1bmx &&
_x2bmn<=x2 && x2<=_x2bmx;
}
// Gets lists of natural neighbor nodes and tris of point (x,y).
// Before building the lists, node and tri marks are cleared. Then,
// as nodes and tris are added to the lists, they are marked, and
// node areas are initialized to zero.
// Returns true, if the lists are not empty; false, otherwise.
private boolean getNaturalNabors(float x, float y) {
_mesh.clearNodeMarks();
_mesh.clearTriMarks();
_nodeList.clear();
_triList.clear();
TriMesh.PointLocation pl = _mesh.locatePoint(x,y);
if (pl.isOutside())
return false;
addTri(x,y,pl.tri());
return true;
}
private void addTri(double xp, double yp, TriMesh.Tri tri) {
_mesh.mark(tri);
_triList.add(tri);
addNode(tri.nodeA());
addNode(tri.nodeB());
addNode(tri.nodeC());
TriMesh.Tri ta = tri.triA();
TriMesh.Tri tb = tri.triB();
TriMesh.Tri tc = tri.triC();
if (needTri(xp,yp,ta)) addTri(xp,yp,ta);
if (needTri(xp,yp,tb)) addTri(xp,yp,tb);
if (needTri(xp,yp,tc)) addTri(xp,yp,tc);
}
private void addNode(TriMesh.Node node) {
if (_mesh.isMarked(node))
return;
_mesh.mark(node);
_nodeList.add(node);
NodeData data = data(node);
data.area = 0.0;
}
private boolean needTri(double xp, double yp, TriMesh.Tri tri) {
if (tri==null || _mesh.isMarked(tri))
return false;
TriMesh.Node na = tri.nodeA();
TriMesh.Node nb = tri.nodeB();
TriMesh.Node nc = tri.nodeC();
double xa = na.xp(), ya = na.yp();
double xb = nb.xp(), yb = nb.yp();
double xc = nc.xp(), yc = nc.yp();
return Geometry.inCircle(xa,ya,xb,yb,xc,yc,xp,yp)>0.0;
}
// C0 interpolation; does not use gradients.
private float interpolate0(double asum) {
double afsum = 0.0;
int nnode = _nodeList.nnode();
TriMesh.Node[] nodes = _nodeList.nodes();
for (int inode=0; inode0.0) {
int nm = _nodeList.nnode();
TriMesh.Node[] ms = _nodeList.nodes();
double hxx = 0.0, hxy = 0.0, hyy = 0.0;
double px = 0.0, py = 0.0;
double nr = 0; // number of real (not ghost) natural neighbor nodes
for (int im=0; im1) { // need at least two real natural neighbor nodes
double lxx = Math.sqrt(hxx); // Cholesky decomposition
double lxy = hxy/lxx;
double dyy = hyy-lxy*lxy;
double lyy = Math.sqrt(dyy);
double qx = px/lxx; // forward elimination
double qy = (py-lxy*qx)/lyy;
double gy = qy/lyy; // back substitution
double gx = (qx-lxy*gy)/lxx;
data.gx = (float)gx;
data.gy = (float)gy;
}
}
}
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
// Given a point (xp,yp) at which to interpolate, an implementation of
// natural neighbor interpolation must accumulate areas for all natural
// neighbor nodes in the the specified node list. This abstract base
// class maintains the total area accumulated for all nodes.
private static abstract class AreaAccumulator {
public abstract double accumulateAreas(
double xp, double yp,
TriMesh mesh, TriMesh.NodeList nodeList, TriMesh.TriList triList);
protected void clear() {
_sum = 0.0;
}
protected double sum() {
return _sum;
}
protected void accumulate(TriMesh.Node node, double area) {
if (ghost(node)) return; // ignore ghost nodes!
NodeData data = data(node);
data.area += area;
_sum += area;
}
private double _sum;
}
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
private static class WatsonSambridge extends AreaAccumulator {
public double accumulateAreas(
double xp, double yp,
TriMesh mesh, TriMesh.NodeList nodeList, TriMesh.TriList triList)
{
clear();
int ntri = triList.ntri();
TriMesh.Tri[] tris = triList.tris();
for (int itri=0; itri edges = _edgeList.edges();
for (int iedge=0; iedge edges = _edgeList.edges();
for (int iedge=0; iedge _edges = new ArrayList(12);
int nedge() {
return _nedge;
}
ArrayList edges() {
return _edges;
}
void clear() {
_nedge = 0;
}
void add(
TriMesh.Node na, TriMesh.Node nb,
double xf, double yf,
double xr, double yr)
{
if (_nedge==_edges.size()) // rarely must we construct a
_edges.add(new Edge()); // new edge like this, because we
Edge edge = _edges.get(_nedge); // can reuse an existing edge
Edge eb = null;
int nfound = 0;
for (int iedge=0; iedge<_nedge && nfound<2; ++iedge) {
Edge ei = _edges.get(iedge);
if (nb==ei.na) {
eb = ei;
++nfound;
}
if (na==ei.nb) {
ei.eb = edge;
++nfound;
}
}
edge.na = na; edge.nb = nb;
edge.xf = xf; edge.yf = yf;
edge.xr = xr; edge.yr = yr;
edge.eb = eb;
++_nedge;
}
}
}
}