edu.mines.jtk.dsp.LocalDiffusionKernel 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.dsp;
import edu.mines.jtk.util.*;
import static edu.mines.jtk.util.ArrayMath.*;
/**
* A local diffusion kernel for use in anisotropic diffusion filtering.
*
* This kernel is a filter that computes y += G'DGx where G is a
* gradient operator, G' is its adjoint, and D is a local diffusion
* tensor field that determines for each image sample the filter
* coefficients.
*
* A local diffusion kernel is typically used in combinations with others.
* For example, the filter implied by (I+G'DG)y = G'DGx acts as a notch
* filter. It attenuates features for which G'DG is zero while preserving
* other features. The diffusion tensors in D control the width, orientation,
* and anisotropy of the spectral notch. Note that application of this filter
* requires solution of a sparse symmetric positive-definite system of
* equations.
*
* An even simpler example is the filter implied by (I+G'DG)y = x. This
* filter smooths features in the directions implied by the tensors D.
* Again, application of this filter requires solving a sparse symmetric
* positive-definite system of equations.
*
* The accumulation of the kernel output in y = y+G'DGx is useful when
* constructing such filter combinations. Given y = 0, this kernel
* computes y = G'DGx. Given y = x, it computes y = (I+G'DG)x.
*
* @author Dave Hale, Colorado School of Mines
* @version 2009.12.31
*/
public class LocalDiffusionKernel {
/**
* The stencil used in finite-difference approximation of derivatives.
* In each stencil name, the first digit equals the number of samples
* used in the direction of the derivative, and the second digit equals
* the number of samples in the orthogonal direction. Names correspond
* to 2D stencils, but each has a natural 3D extension.
*
* Note that the stencil implied by G'DG is larger than that used to
* approximate the derivatives in G. For example, a 2x2 derivative
* approximation implies a 3x3 stencil for G'DG.
*/
public enum Stencil {
/**
* A 2x1 stencil.
* The 2D version has 3 non-zero coefficients.
* The 3D version has 4 non-zero coefficients.
* This stencil should be specified for isotropic diffusion only.
* When using this stencil, any specified tensors D are ignored.
*/
D21,
/**
* A 2x2 stencil.
* The 2D version has 4 non-zero coefficients.
* The 3D version has 8 non-zero coefficients.
* This stencil is the default.
*/
D22,
/**
* A 2x4 stencil.
* The 2D version has 8 non-zero coefficients.
* The 3D version has 24 non-zero coefficients.
* The 3D version is not yet implemented.
*/
D24,
/**
* A 3x3 stencil.
* The 2D version has 6 non-zero coefficients.
* The 3D version has 18 non-zero coefficients.
*/
D33,
/** A 7x1 stencil.
* Both 2D and 3D versions have 6 non-zero coefficients.
*/
D71,
/** A 9x1 stencil.
* Both 2D and 3D versions have 8 non-zero coefficients.
* The 3D version is not yet implemented.
*/
D91,
}
/**
* Constructs a local diffusion kernel with default 2x2 stencil.
*/
public LocalDiffusionKernel() {
this(Stencil.D22);
}
/**
* Constructs a local diffusion kernel with specified stencil.
* @param s the stencil used to approximate a derivative.
*/
public LocalDiffusionKernel(Stencil s) {
_stencil = s;
}
/**
* Gets the stencil used by this local diffusion kernel.
* @return the stencil.
*/
public Stencil getStencil() {
return _stencil;
}
/**
* Sets the number of kernel passes in each apply of this filter.
* For example, if npass = 2, then the output is computed in two
* passes: (1) y += G'DGx, (2) y += G'DGy.
* The default is one pass.
* @param npass the number of passes.
*/
public void setNumberOfPasses(int npass) {
_npass = npass;
}
/**
* Applies this filter for constant isotropic identity tensor.
* @param x input array.
* @param y output array.
*/
public void apply(float[][] x, float[][] y)
{
apply(null,1.0f,x,y);
}
/**
* Applies this filter for specified tensor coefficients.
* @param d tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(Tensors2 d, float[][] x, float[][] y)
{
apply(d,1.0f,x,y);
}
/**
* Applies this filter for specified scale factor.
* Uses a constant isotropic identity tensor.
* @param c constant scale factor for tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(float c, float[][] x, float[][] y) {
apply(null,c,null,x,y);
}
/**
* Applies this filter for specified tensor coefficients and scale factor.
* @param d tensor coefficients.
* @param c constant scale factor for tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(Tensors2 d, float c, float[][] x, float[][] y) {
apply(d,c,null,x,y);
}
/**
* Applies this filter for specified scale factors.
* Uses a constant isotropic identity tensor.
* @param c constant scale factor for tensor coefficients.
* @param s array of scale factors for tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(float c, float[][] s, float[][] x, float[][] y) {
apply(null,c,s,x,y);
}
/**
* Applies this filter for specified tensor coefficients and scale factors.
* @param d tensor coefficients.
* @param c constant scale factor for tensor coefficients.
* @param s array of scale factors for tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(
Tensors2 d, float c, float[][] s, float[][] x, float[][] y)
{
for (int ipass=0; ipass<_npass; ++ipass) {
if (ipass>0)
x = copy(y);
if (d==null)
d = IDENTITY_TENSORS2;
if (_stencil==Stencil.D21) {
apply21(c,s,x,y);
} else if (_stencil==Stencil.D22) {
apply22(d,c,s,x,y);
} else if (_stencil==Stencil.D24) {
apply24(d,c,s,x,y);
} else if (_stencil==Stencil.D33) {
apply33(d,c,s,x,y);
} else if (_stencil==Stencil.D71) {
apply71(d,c,s,x,y);
} else if (_stencil==Stencil.D91) {
apply91(d,c,s,x,y);
}
}
}
/**
* Applies this filter for a constant isotropic identity tensor.
* @param x input array.
* @param y output array.
*/
public void apply(float[][][] x, float[][][] y) {
apply(null,1.0f,x,y);
}
/**
* Applies this filter for specified tensor coefficients.
* @param d tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(Tensors3 d, float[][][] x, float[][][] y) {
apply(d,1.0f,x,y);
}
/**
* Applies this filter for specified scale factor.
* Uses a constant isotropic identity tensor.
* @param c constant scale factor for tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(float c, float[][][] x, float[][][] y) {
apply(null,c,null,x,y);
}
/**
* Applies this filter for specified tensor coefficients and scale factor.
* @param d tensor coefficients.
* @param c constant scale factor for tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(Tensors3 d, float c, float[][][] x, float[][][] y) {
apply(d,c,null,x,y);
}
/**
* Applies this filter for specified scale factors.
* Uses a constant isotropic identity tensor.
* @param c constant scale factor for tensor coefficients.
* @param s array of scale factors for tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(float c, float[][][] s, float[][][] x, float[][][] y) {
apply(null,c,s,x,y);
}
/**
* Applies this filter for specified tensor coefficients and scale factors.
* @param d tensor coefficients.
* @param c constant scale factor for tensor coefficients.
* @param s array of scale factors for tensor coefficients.
* @param x input array.
* @param y output array.
*/
public void apply(
Tensors3 d, float c, float[][][] s, float[][][] x, float[][][] y)
{
int n3 = x.length;
int i3start = 0;
int i3step = 1;
int i3stop = n3;
for (int ipass=0; ipass<_npass; ++ipass) {
if (ipass>0)
x = copy(y);
if (d==null)
d = IDENTITY_TENSORS3;
if (_stencil==Stencil.D21) {
i3start = 0; i3step = 2; i3stop = n3;
} else if (_stencil==Stencil.D22) {
i3start = 1; i3step = 2; i3stop = n3;
} else if (_stencil==Stencil.D24) {
i3start = 1; i3step = 4; i3stop = n3;
} else if (_stencil==Stencil.D33) {
i3start = 1; i3step = 3; i3stop = n3-1;
} else if (_stencil==Stencil.D71) {
i3start = 0; i3step = 7; i3stop = n3;
}
if (_parallel) {
applyParallel(i3start,i3step,i3stop,d,c,s,x,y);
} else {
applySerial(i3start,1,i3stop,d,c,s,x,y);
}
}
}
///////////////////////////////////////////////////////////////////////////
// private
private static Tensors2 IDENTITY_TENSORS2 = new Tensors2() {
public void getTensor(int i1, int i2, float[] d) {
d[0] = 1.0f;
d[1] = 0.0f;
d[2] = 1.0f;
}
};
private static Tensors3 IDENTITY_TENSORS3 = new Tensors3() {
public void getTensor(int i1, int i2, int i3, float[] d) {
d[0] = 1.0f;
d[1] = 0.0f;
d[2] = 0.0f;
d[3] = 1.0f;
d[4] = 0.0f;
d[5] = 1.0f;
}
};
private Stencil _stencil;
private int _npass = 1;
private boolean _parallel = true;
private static void trace(String s) {
System.out.println(s);
}
private void apply(
int i3, Tensors3 d, float c, float[][][] s, float[][][] x, float[][][] y)
{
if (_stencil==Stencil.D21) {
apply21(i3,c,s,x,y);
} else if (_stencil==Stencil.D22) {
apply22(i3,d,c,s,x,y);
} else if (_stencil==Stencil.D24) {
//apply24(i3,d,c,s,x,y);
throw new UnsupportedOperationException(
"Stencil.D24 not supported for 3D arrays");
} else if (_stencil==Stencil.D33) {
apply33(i3,d,c,s,x,y);
} else if (_stencil==Stencil.D71) {
apply71(i3,d,c,s,x,y);
} else if (_stencil==Stencil.D91) {
//apply91(i3,d,c,s,x,y);
throw new UnsupportedOperationException(
"Stencil.D91 not supported for 3D arrays");
}
}
private void applySerial(
int i3start, int i3step, int i3stop,
Tensors3 d, float c, float[][][] s, float[][][] x, float[][][] y)
{
for (int i3=i3start; i30)?i2-1:0;
for (int i1=0; i10)?i1-1:0;
float cs1 = c;
float cs2 = c;
if (s!=null) {
cs1 *= 0.5f*(s[i2][i1]+s[i2][m1]);
cs2 *= 0.5f*(s[i2][i1]+s[m2][i1]);
}
float x1 = x[i2][i1]-x[i2][m1];
float x2 = x[i2][i1]-x[m2][i1];
float y1 = cs1*x1;
float y2 = cs2*x2;
y[i2][i1] += y1+y2;
y[i2][m1] -= y1;
y[m2][i1] -= y2;
}
}
}
private void apply21(
int i3, float c, float[][][] s, float[][][] x, float[][][] y)
{
int n1 = x[0][0].length;
int n2 = x[0].length;
int m3 = (i3>0)?i3-1:0;
for (int i2=0; i20)?i2-1:0;
for (int i1=0; i10)?i1-1:0;
float cs1 = c;
float cs2 = c;
float cs3 = c;
if (s!=null) {
cs1 *= 0.5f*(s[i3][i2][i1]+s[i3][i2][m1]);
cs2 *= 0.5f*(s[i3][i2][i1]+s[i3][m2][i1]);
cs3 *= 0.5f*(s[i3][i2][i1]+s[m3][i2][i1]);
}
float x1 = x[i3][i2][i1]-x[i3][i2][m1];
float x2 = x[i3][i2][i1]-x[i3][m2][i1];
float x3 = x[i3][i2][i1]-x[m3][i2][i1];
float y1 = cs1*x1;
float y2 = cs2*x2;
float y3 = cs3*x3;
y[i3][i2][i1] += y1+y2+y3;
y[i3][i2][m1] -= y1;
y[i3][m2][i1] -= y2;
y[m3][i2][i1] -= y3;
}
}
}
///////////////////////////////////////////////////////////////////////////
// D22
private void apply22(
Tensors2 d, float c, float[][] s, float[][] x, float[][] y)
{
c *= 0.25f;
int n1 = x[0].length;
int n2 = x.length;
float[] di = new float[3];
for (int i2=1; i2=n1) i2p1 = n1-1;
float[] xm2=x[i2m2], xm1=x[i2m1], xp0=x[i2p0], xp1=x[i2p1];
float[] ym2=y[i2m2], ym1=y[i2m1], yp0=y[i2p0], yp1=y[i2p1];
int m2, m1 = 0, p0 = 0, p1 = 1;
for (int i1=1; i1=n1) p1 = n1-1;
d.getTensor(i1,i2,di);
float csi = (s!=null)?c*s[i2][i1]:c;
float d11 = di[0]*csi;
float d12 = di[1]*csi;
float d22 = di[2]*csi;
float xa = xp0[p0]-xm1[m1];
float xb = xm1[p0]-xp0[m1];
float x1 = xa+xb+b*(xp1[p0]+xm2[p0]-xp1[m1]-xm2[m1]);
float x2 = xa-xb+b*(xp0[p1]+xp0[m2]-xm1[p1]-xm1[m2]);
float y1 = d11*x1+d12*x2;
float y2 = d12*x1+d22*x2;
float ya = y1+y2;
float yb = y1-y2;
float yc = b*y1;
float yd = b*y2;
yp0[p0] += ya; ym1[m1] -= ya;
ym1[p0] += yb; yp0[m1] -= yb;
yp1[p0] += yc; ym2[m1] -= yc;
ym2[p0] += yc; yp1[m1] -= yc;
yp0[p1] += yd; ym1[m2] -= yd;
yp0[m2] += yd; ym1[p1] -= yd;
}
}
}
///////////////////////////////////////////////////////////////////////////
// D33
private void apply33(
Tensors2 d, float c, float[][] s, float[][] x, float[][] y)
{
float p = 0.182962f; // Scharr, best for high anisotropy
float a = 0.5f-p; // ~ 10/32
float b = 0.5f*p; // ~ 3/32
b /= a;
c *= a*a;
int n1 = x[0].length;
int n2 = x.length;
float[] di = new float[3];
for (int i2=1; i2=n2) i2p1 = n2-1;
if (i2p2>=n2) i2p2 = n2-1;
if (i2p3>=n2) i2p3 = n2-1;
float[] xm3 = x[i2m3], xm2 = x[i2m2], xm1 = x[i2m1];
float[] xp3 = x[i2p3], xp2 = x[i2p2], xp1 = x[i2p1];
float[] xp0 = x[i2p0];
float[] ym3 = y[i2m3], ym2 = y[i2m2], ym1 = y[i2m1];
float[] yp3 = y[i2p3], yp2 = y[i2p2], yp1 = y[i2p1];
float[] yp0 = y[i2p0];
int m3,m2=0,m1=0,p0=0,p1=0,p2=1,p3=2;
for (int i1=0; i1=n1) p1 = n1-1;
if (p2>=n1) p2 = n1-1;
if (p3>=n1) p3 = n1-1;
d.getTensor(i1,i2,di);
float csi = (s!=null)?c*s[i2][i1]:c;
float d11 = di[0]*csi;
float d12 = di[1]*csi;
float d22 = di[2]*csi;
float x1 = c1*(xp0[p1]-xp0[m1]) +
c2*(xp0[p2]-xp0[m2]) +
c3*(xp0[p3]-xp0[m3]);
float x2 = c1*(xp1[p0]-xm1[p0]) +
c2*(xp2[p0]-xm2[p0]) +
c3*(xp3[p0]-xm3[p0]);
float y1 = d11*x1+d12*x2;
float y2 = d12*x1+d22*x2;
float c1y1 = c1*y1; yp0[p1] += c1y1; yp0[m1] -= c1y1;
float c2y1 = c2*y1; yp0[p2] += c2y1; yp0[m2] -= c2y1;
float c3y1 = c3*y1; yp0[p3] += c3y1; yp0[m3] -= c3y1;
float c1y2 = c1*y2; yp1[p0] += c1y2; ym1[p0] -= c1y2;
float c2y2 = c2*y2; yp2[p0] += c2y2; ym2[p0] -= c2y2;
float c3y2 = c3*y2; yp3[p0] += c3y2; ym3[p0] -= c3y2;
}
}
}
private void apply71X(
int i3, Tensors3 d, float c, float[][][] s, float[][][] x, float[][][] y)
{
final float c1 = C71[1], c2 = C71[2], c3 = C71[3];
int n1 = x[0][0].length;
int n2 = x[0].length;
int n3 = x.length;
float[] di = new float[6];
int i3m3 = max(0,i3-3), i3p3 = min(n3-1,i3+3);
int i3m2 = max(0,i3-2), i3p2 = min(n3-1,i3+2);
int i3m1 = max(0,i3-1), i3p1 = min(n3-1,i3+1);
float[][] g1 = new float[n2][n1];
float[][] g2 = new float[n2][n1];
gf(C71,x[i3],g1,g2);
for (int i2=0; i2=n3) i3p1 = n3-1;
int i3p2 = i3+2; if (i3p2>=n3) i3p2 = n3-1;
int i3p3 = i3+3; if (i3p3>=n3) i3p3 = n3-1;
int i2m3,i2m2=0,i2m1=0,i2p0=0,i2p1=0,i2p2=1,i2p3=2;
for (int i2=0; i2=n2) i2p1 = n2-1;
if (i2p2>=n2) i2p2 = n2-1;
if (i2p3>=n2) i2p3 = n2-1;
float[] xp0p0 = x[i3p0][i2p0], yp0p0 = y[i3p0][i2p0];
float[] xp0m3 = x[i3p0][i2m3], yp0m3 = y[i3p0][i2m3];
float[] xp0m2 = x[i3p0][i2m2], yp0m2 = y[i3p0][i2m2];
float[] xp0m1 = x[i3p0][i2m1], yp0m1 = y[i3p0][i2m1];
float[] xp0p1 = x[i3p0][i2p1], yp0p1 = y[i3p0][i2p1];
float[] xp0p2 = x[i3p0][i2p2], yp0p2 = y[i3p0][i2p2];
float[] xp0p3 = x[i3p0][i2p3], yp0p3 = y[i3p0][i2p3];
float[] xm3p0 = x[i3m3][i2p0], ym3p0 = y[i3m3][i2p0];
float[] xm2p0 = x[i3m2][i2p0], ym2p0 = y[i3m2][i2p0];
float[] xm1p0 = x[i3m1][i2p0], ym1p0 = y[i3m1][i2p0];
float[] xp1p0 = x[i3p1][i2p0], yp1p0 = y[i3p1][i2p0];
float[] xp2p0 = x[i3p2][i2p0], yp2p0 = y[i3p2][i2p0];
float[] xp3p0 = x[i3p3][i2p0], yp3p0 = y[i3p3][i2p0];
int m3,m2=0,m1=0,p0=0,p1=0,p2=1,p3=2;
for (int i1=0; i1=n1) p1 = n1-1;
if (p2>=n1) p2 = n1-1;
if (p3>=n1) p3 = n1-1;
d.getTensor(i1,i2,i3,di);
float csi = (s!=null)?c*s[i3][i2][i1]:c;
float d11 = di[0]*csi;
float d12 = di[1]*csi;
float d13 = di[2]*csi;
float d22 = di[3]*csi;
float d23 = di[4]*csi;
float d33 = di[5]*csi;
float x1 = c1*(xp0p0[p1]-xp0p0[m1]) +
c2*(xp0p0[p2]-xp0p0[m2]) +
c3*(xp0p0[p3]-xp0p0[m3]);
float x2 = c1*(xp0p1[p0]-xp0m1[p0]) +
c2*(xp0p2[p0]-xp0m2[p0]) +
c3*(xp0p3[p0]-xp0m3[p0]);
float x3 = c1*(xp1p0[p0]-xm1p0[p0]) +
c2*(xp2p0[p0]-xm2p0[p0]) +
c3*(xp3p0[p0]-xm3p0[p0]);
float y1 = d11*x1+d12*x2+d13*x3;
float y2 = d12*x1+d22*x2+d23*x3;
float y3 = d13*x1+d23*x2+d33*x3;
float c1y1 = c1*y1; yp0p0[p1] += c1y1; yp0p0[m1] -= c1y1;
float c2y1 = c2*y1; yp0p0[p2] += c2y1; yp0p0[m2] -= c2y1;
float c3y1 = c3*y1; yp0p0[p3] += c3y1; yp0p0[m3] -= c3y1;
float c1y2 = c1*y2; yp0p1[p0] += c1y2; yp0m1[p0] -= c1y2;
float c2y2 = c2*y2; yp0p2[p0] += c2y2; yp0m2[p0] -= c2y2;
float c3y2 = c3*y2; yp0p3[p0] += c3y2; yp0m3[p0] -= c3y2;
float c1y3 = c1*y3; yp1p0[p0] += c1y3; ym1p0[p0] -= c1y3;
float c2y3 = c2*y3; yp2p0[p0] += c2y3; ym2p0[p0] -= c2y3;
float c3y3 = c3*y3; yp3p0[p0] += c3y3; ym3p0[p0] -= c3y3;
}
}
}
///////////////////////////////////////////////////////////////////////////
// D91
private static final float[] C91 = {
0.0f, 0.8947167f, -0.3153471f, 0.1096895f, -0.0259358f
};
private void apply91(
Tensors2 d, float c, float[][] s, float[][] x, float[][] y)
{
float c1 = C91[1], c2 = C91[2], c3 = C91[3], c4 = C91[4];
int n1 = x[0].length;
int n2 = x.length;
float[] di = new float[3];
int i2m4,i2m3=0,i2m2=0,i2m1=0,i2p0=0,i2p1=0,i2p2=1,i2p3=2,i2p4=3;
for (int i2=0; i2=n2) i2p1 = n2-1;
if (i2p2>=n2) i2p2 = n2-1;
if (i2p3>=n2) i2p3 = n2-1;
if (i2p4>=n2) i2p4 = n2-1;
float[] xm4 = x[i2m4], xm3 = x[i2m3], xm2 = x[i2m2], xm1 = x[i2m1];
float[] xp4 = x[i2p4], xp3 = x[i2p3], xp2 = x[i2p2], xp1 = x[i2p1];
float[] xp0 = x[i2p0];
float[] ym4 = y[i2m4], ym3 = y[i2m3], ym2 = y[i2m2], ym1 = y[i2m1];
float[] yp4 = y[i2p4], yp3 = y[i2p3], yp2 = y[i2p2], yp1 = y[i2p1];
float[] yp0 = y[i2p0];
int m4,m3=0,m2=0,m1=0,p0=0,p1=0,p2=1,p3=2,p4=3;
for (int i1=0; i1=n1) p1 = n1-1;
if (p2>=n1) p2 = n1-1;
if (p3>=n1) p3 = n1-1;
if (p4>=n1) p4 = n1-1;
d.getTensor(i1,i2,di);
float csi = (s!=null)?c*s[i2][i1]:c;
float d11 = di[0]*csi;
float d12 = di[1]*csi;
float d22 = di[2]*csi;
float x1 = c1*(xp0[p1]-xp0[m1]) +
c2*(xp0[p2]-xp0[m2]) +
c3*(xp0[p3]-xp0[m3]) +
c4*(xp0[p4]-xp0[m4]);
float x2 = c1*(xp1[p0]-xm1[p0]) +
c2*(xp2[p0]-xm2[p0]) +
c3*(xp3[p0]-xm3[p0]) +
c4*(xp4[p0]-xm4[p0]);
float y1 = d11*x1+d12*x2;
float y2 = d12*x1+d22*x2;
float c1y1 = c1*y1; yp0[p1] += c1y1; yp0[m1] -= c1y1;
float c2y1 = c2*y1; yp0[p2] += c2y1; yp0[m2] -= c2y1;
float c3y1 = c3*y1; yp0[p3] += c3y1; yp0[m3] -= c3y1;
float c4y1 = c4*y1; yp0[p4] += c4y1; yp0[m4] -= c4y1;
float c1y2 = c1*y2; yp1[p0] += c1y2; ym1[p0] -= c1y2;
float c2y2 = c2*y2; yp2[p0] += c2y2; ym2[p0] -= c2y2;
float c3y2 = c3*y2; yp3[p0] += c3y2; ym3[p0] -= c3y2;
float c4y2 = c4*y2; yp4[p0] += c4y2; ym4[p0] -= c4y2;
}
}
}
//////////////////////////////////////////////////////////////////////////
// gradients (forward and transpose) for stencils with mx1 coefficients
private static void gf(float[] c, float[] x, float[] y) {
int nc = c.length-1;
int n1 = x.length;
int n1m1 = n1-1, n1nc = n1-nc;
for (int i1=0; i1n1m1) ip = n1m1;
yi += ci*(x[ip]-x[im]);
}
y[i1] = yi;
}
if (nc==3 && n1>6) { // middle part optimized for nc = 3
final float c1 = c[1], c2 = c[2], c3 = c[3];
float xm3, xm2 = x[0], xm1 = x[1], xp0 = x[2],
xp1 = x[3], xp2 = x[4], xp3 = x[5];
for (int i1=3; i1n1m1) ip = n1m1;
yi += ci*(x[ip]-x[im]);
}
y[i1] = yi;
}
}
private static void gt(float[] c, float[] x, float[] y) {
int nc = c.length-1;
int n1 = x.length;
int n1m1 = n1-1, n1nc = n1-nc;
for (int i1=0; i1n1m1) ip = n1m1;
if (im6) { // middle part optimized for nc = 3
final float c1 = c[1], c2 = c[2], c3 = c[3];
float xm3, xm2 = x[0], xm1 = x[1], xp0 = x[2],
xp1 = x[3], xp2 = x[4], xp3 = x[5];
for (int i1=3; i1n1m1) ip = n1m1;
if (im>=n1nc) y[im] -= ci*xi;
if (ip>=n1nc) y[ip] += ci*xi;
}
}
}
private static void gf1(float[] c, float[][] x, float[][] g1) {
int n2 = x.length;
for (int i2=0; i2=3)?x[i2-3]:x[0];
float[] xm2 = (i2>=2)?x[i2-2]:x[0];
float[] xm1 = (i2>=1)?x[i2-1]:x[0];
float[] xp1 = (i2=ic)?x[i2-ic]:x[0];
float[] xp = (i2n2m1) ip = n2m1;
if (im6) { // middle part optimized for nc = 3
final float c1 = c[1], c2 = c[2], c3 = c[3];
for (int i2=3; i2n2m1) ip = n2m1;
if (im>=n2nc) {
float[] x2m = x[im];
for (int i1=0; i1=n2nc) {
float[] x2p = x[ip];
for (int i1=0; i1© 2015 - 2025 Weber Informatics LLC | Privacy Policy