edu.mines.jtk.interp.SplinesGridder2 Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2010, 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 java.util.logging.Logger;
import edu.mines.jtk.dsp.*;
import edu.mines.jtk.util.Check;
import static edu.mines.jtk.util.ArrayMath.*;
/**
* Tensor-guided 2D gridding with bi-harmonic and harmonic splines.
* At locations where gridded values are not constrained by specified
* (known) sample values, the gridded values q satisfy the equation
* (G'DGG'DG+tG'DG)q = 0, where G is a finite-difference approximation
* of the gradient operator, G' is its transpose, D is a tensor field,
* and t is a scalar constant that controls the tension, the weight of
* the harmonic G'DG operator relative to the bi-harmonic G'DGG'DG
* operator.
*
* In the isotropic and homogeneous case, where D = I is the identity
* matrix, this gridder is an implementation of Smith and Wessel's
* (1990) gridding with continuous curvature splines in tension. The
* finite-difference approximations to derivatives in G are different,
* however, because of the tensors D between G' and G in G'DG. These
* tensors D must be symmetric and positive-semidefinite (SPSD) so that
* the finite-difference operators G'DG and G'DGG'DG are SPSD as well.
*
* Another difference between this implementation and that of Smith and
* Wessel (1990) is that a multigrid method is not used here to find the
* gridded values q. Multigrid methods are ineffective for tensors fields
* D that are anisotropic and inhomogeneous. Instead, this implementation
* uses conjugate-gradient (CG) iterations.
*
* The gridded values q must be obtained by solving iteratively the large
* sparse system of equations (G'DGG'DG+tG'DG)q = 0. To facilitate a CG
* solver, these equations are rewritten as (K+MAM)q = (K-MAK)q, where
* A = G'DGG'DG+tG'DG, M is a diagonal matrix operator with ones where
* sample values are missing and zeros where values are known, and
* K = I-M is a diagonal matrix with ones where sample values are known
* and zero where they are missing. The matrix K+MAM on the left-hand side
* is symmetric and positive-definite (SPD), so that CG iterations can be
* used to solve for the unknown values in q. Only known values in q appear
* in the right-hand side column vector (K-MAK)q.
*
* See Smith, W.H.F., and P. Wessel, 1990, Gridding with continuous
* curvature splines in tension: Geophysics, 55, 293--305.
*
* @author Dave Hale, Colorado School of Mines
* @version 2010.01.19
*/
public class SplinesGridder2 implements Gridder2 {
/**
* Constructs a gridder for default tensors.
*/
public SplinesGridder2() {
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 SplinesGridder2(float[] f, float[] x1, float[] x2) {
this(null);
setScattered(f,x1,x2);
}
/**
* Constructs a gridder for the specified tensors.
* @param tensors the tensors.
*/
public SplinesGridder2(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 SplinesGridder2(
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;
}
};
}
}
/**
* Sets the tension, the weight for the harmonic spline.
* The default tension is 0.0, for a purely bi-harmonic spline.
*
* This tension parameter is not equivalent to the parameter t in
* the class documentation above. The latter would be infinite if
* the tension specified here could be set to 1.0.
* @param tension the tension; must be in the range [0:1).
*/
public void setTension(double tension) {
Check.argument(0<=tension,"0<=tension");
Check.argument(tension<1,"tension<1");
_tension = (float)tension;
}
/**
* Sets the maximum number of conjugate-gradient iterations.
* The default maximum number of iterations is 10,000.
* @param niter the maximum number of iterations.
*/
public void setMaxIterations(int niter) {
_niter = niter;
}
/**
* Returns the number of conjugate-gradient iterations required.
* The number returned corresponds to the last use of this gridder.
* @return the number of iterations.
*/
public int getIterationCount() {
return _residuals.size()-1;
}
/**
* Gets the initial residual and one residual for each iteration.
* The residuals returned correspond to the last use of this gridder.
*
* Residuals are normalized root-mean-square differences between
* the left and right sides of the system of equations that are
* solved iteratively when computing gridded values. The returned
* residuals are normalized, so that the zeroth residual (before any
* conjugate-gradient iterations are performed) is one.
* @return array of residuals.
*/
public float[] getResiduals() {
int n = _residuals.size();
float[] r = new float[n];
for (int i=0; i _residuals = new ArrayList();
private LocalDiffusionKernel _ldk =
new LocalDiffusionKernel(LocalDiffusionKernel.Stencil.D22);
private static Logger log =
Logger.getLogger(SplinesGridder2.class.getName());
private static interface Operator2 {
public void apply(float[][] x, float[][] y);
}
// Smoothing operator SS used for preconditioning. This operator
// attenuates frequencies near the Nyquist limit for which
// finite-difference approximations in G'DG are poor.
private static class SmoothOperator2 implements Operator2 {
public void apply(float[][] x, float[][] y) {
smoothS(x,y);
smoothS(y,y);
}
}
// The left-hand-side operator K + M(G'DGG'DG + tG'DG)M.
// Can also apply the right-hand-side operator K - M(...)K.
private static class LaplaceOperator2 implements Operator2 {
LaplaceOperator2(
LocalDiffusionKernel ldk, Tensors2 d, float t, boolean[][] m)
{
_ldk = ldk;
_d = d;
_t = t;
_m = m;
_z = new float[m.length][m[0].length];
}
public void apply(float[][] x, float[][] y) {
int n1 = x[0].length;
int n2 = x.length;
// z = Mx
for (int i2=0; i2© 2015 - 2025 Weber Informatics LLC | Privacy Policy