edu.mines.jtk.interp.BlendedGridder2 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 edu.mines.jtk.dsp.*;
import edu.mines.jtk.util.Check;
import static edu.mines.jtk.util.ArrayMath.*;
/**
* Tensor-guided blended neighbor gridding in 2D.
* Gridding is interpolation of a set of known sample values on a
* uniformly sampled grid. Here the interpolation is performed by
* a two-step blended neighbor process described by Hale (2009).
*
* The first step is to compute for all samples the distance to the
* nearest known sample and the value of that known sample. This first
* step produces a distance map and a nearest-neighbor interpolant.
*
* The second step is to blend (smooth) the nearest-neighbor interpolant,
* where the extent of smoothing varies spatially and is proportional to
* distances in the distance map.
*
* In tensor-guided gridding, we replace distance with time. Time is a
* simple term for non-Euclidean distance measured in a metric-tensor
* field. So "nearest" now means nearest in time. In the first step we
* compute a time map by solving an eikonal equation with coefficients
* that may be both anisotropic and spatially varying. In the second
* step, we blend the nearest-neighbor interpolant with an anisotropic
* and spatially varying smoothing filter.
*
* The default tensor field is homogeneous and isotropic. In this
* special case, time is equivalent to distance, and tensor-guided
* gridding is similar to gridding with Sibson's natural neighbor
* interpolant.
*
* Reference:
*
* Hale, D., 2009, Image-guided blended neighbor interpolation, CWP-634
* @author Dave Hale, Colorado School of Mines
* @version 2009.07.21
*/
public class BlendedGridder2 implements Gridder2 {
/**
* Constructs a gridder for default tensors.
*/
public BlendedGridder2() {
this(null);
}
/**
* Constructs a gridder for default tensors and specified samples.
* The specified arrays are referenced; not copied.
* @param f array of sample values f(x1,x2).
* @param x1 array of sample x1 coordinates.
* @param x2 array of sample x2 coordinates.
*/
public BlendedGridder2(float[] f, float[] x1, float[] x2) {
this(null);
setScattered(f,x1,x2);
}
/**
* Constructs a gridder for the specified tensors.
* @param tensors the tensors.
*/
public BlendedGridder2(Tensors2 tensors) {
setTensors(tensors);
}
/**
* Constructs a gridder for the specified tensors and samples.
* The specified arrays are referenced; not copied.
* @param tensors the tensors.
* @param f array of sample values f(x1,x2).
* @param x1 array of sample x1 coordinates.
* @param x2 array of sample x2 coordinates.
*/
public BlendedGridder2(
Tensors2 tensors,
float[] f, float[] x1, float[] x2)
{
setTensors(tensors);
setScattered(f,x1,x2);
}
/**
* Sets the tensor field used by this gridder.
* The default is a homogeneous and isotropic tensor field.
* @param tensors the tensors; null for default tensors.
*/
public void setTensors(Tensors2 tensors) {
_tensors = tensors;
if (_tensors==null) {
_tensors = new Tensors2() {
public void getTensor(int i1, int i2, float[] d) {
d[0] = 1.0f;
d[1] = 0.0f;
d[2] = 1.0f;
}
};
}
}
/**
* Enables or disables blending in {@link #grid(Sampling,Sampling)}.
* If true (the default), that method will perform both of the two
* steps described; that is, it will blend (smooth) after computing
* the nearest neighbor interpolant. If false, that method will perform
* only the first step and return the nearest neighbor interpolant. The
* default is that blending is enabled.
* @param blending true, for blending; false, otherwise.
*/
public void setBlending(boolean blending) {
_blending = blending;
}
/**
* Sets the local diffusion kernel used to perform blending.
* The default kernel uses a 2x2 stencil to compute derivatives.
* @param ldk the local diffusion kernel.
*/
public void setBlendingKernel(LocalDiffusionKernel ldk) {
_ldk = ldk;
}
/**
* Sets the smoothness of the interpolation of gridded values.
* The default is 0.5, which yields an interpolant with linear precision.
* Larger values yield smoother interpolants with plateaus at known samples.
* @param smoothness the smoothness.
*/
public void setSmoothness(double smoothness) {
_c = 0.25f/(float)smoothness;
}
/**
* Sets the maximum time computed by this gridder. The gridder has
* linear precision where times are less than the maximum time.
* @param tmax the maximum time.
*/
public void setTimeMax(double tmax) {
_tmax = (float)tmax;
}
/**
* Experimental use only.
* @param tmx time marker x.
*/
public void setTimeMarkerX(boolean tmx) {
_tmx = tmx;
}
/**
* Experimental use only.
* @return time marker s.
*/
public double getTimeMarkerS() {
return _tms;
}
/**
* Computes gridded values using nearest neighbors.
* Gridded values in the array p are computed for only unknown
* samples with value equal to the specified null value. Any
* known (non-null) sample values in the array p are not changed.
*
* This method also computes and returns an array of times to
* nearest-neighbor samples. Times are zero for known samples
* and are positive for gridded samples.
* @param pnull the null value representing unknown samples.
* @param p array of sample values to be gridded.
* @return array of times to nearest known samples.
*/
public float[][] gridNearest(float pnull, float[][] p) {
int n1 = p[0].length;
int n2 = p.length;
// Make array of times, zero for samples with non-null values.
// Also, count the number of marks that we will need.
int nmark = 0;
float[][] t = new float[n2][n1];
float tnull = -FLT_MAX;
for (int i2=0; i20; --i2) {
for (int i1=n1-1; i1>0; --i1) {
s[i2][i1] = 0.25f*(s[i2 ][i1 ] +
s[i2 ][i1-1] +
s[i2-1][i1 ] +
s[i2-1][i1-1]);
}
}
}
// Construct and apply a local smoothing filter.
LocalSmoothingFilter lsf = new LocalSmoothingFilter(0.01,10000,_ldk);
lsf.setPreconditioner(true);
float pavg = sum(p)/n1/n2;
float[][] r = sub(p,pavg);
// Smoothing attenuates finite-difference errors near Nyquist.
// The problem with this smoothing is that it makes q != p when
// known samples are adjacent, as when interpolating well logs.
// This smoothing should be unnecessary for Stencil.D21.
if (_ldk.getStencil()!=LocalDiffusionKernel.Stencil.D21)
lsf.applySmoothS(r,r);
lsf.apply(_tensors,_c,s,r,q);
add(q,pavg,q);
// Restore the known sample values. Due to errors in finite-difference
// approximations, these values may have changed during smoothing.
// Even with time adjustments, this restoration is still necessary
// if we used applySmoothS above. Best to just do this in any case.
for (int i2=0; i2_tmax)
t[i2][i1] = _tmax;
}
}
}
// Adjusts times to be nearly zero in neighborhood of known samples.
private void adjustTimes(int nmark, int[][] m, float[][] t) {
int n1 = m[0].length;
int n2 = m.length;
int n1m = n1-1;
int n2m = n2-1;
float[] s = new float[nmark];
for (int i2=0; i20 )?i2-1:i2;
int i2p = (i20 )?i1-1:i1;
int i1p = (i10.0f)
t[i2][i1] = max(FLT_MIN,t[i2][i1]-s[m[i2][i1]]);
}
}
}
}