/** * Applies this filter for specified constant scale factor. * Local smoothing for 1D arrays is a special case that requires no tensors. * All tensors are implicitly scalar values equal to one, so that filtering * is determined entirely by the specified constant scale factor. * @param c constant scale factor. * @param x input array. * @param y output array. */ public void apply(float c, float[] x, float[] y) { apply(c,null,x,y); }
private static void smoothS(float[][] x, float[][] y) { int n1 = x[0].length; int n2 = x.length; int n1m = n1-1; int n2m = n2-1; float[][] t = new float[3][n1]; scopy(x[0],t[0]); scopy(x[0],t[1]); for (int i2=0; i2<n2; ++i2) { int i2m = (i2>0)?i2-1:0; int i2p = (i2<n2m)?i2+1:n2m; int j2m = i2m%3; int j2 = i2%3; int j2p = i2p%3; scopy(x[i2p],t[j2p]); float[] x2m = t[j2m]; float[] x2p = t[j2p]; float[] x20 = t[j2]; float[] y2 = y[i2]; for (int i1=0; i1<n1; ++i1) { int i1m = (i1>0)?i1-1:0; int i1p = (i1<n1m)?i1+1:n1m; y2[i1] = 0.2500f*(x20[i1 ]) + 0.1250f*(x20[i1m]+x20[i1p]+x2m[i1 ]+x2p[i1 ]) + 0.0625f*(x2m[i1m]+x2m[i1p]+x2p[i1m]+x2p[i1p]); } } }
/** * Applies an isotropic low-pass smoothing filter L. * Input and output arrays x and y may be the same array. * @param kmax maximum wavenumber not attenuated, in cycles/sample. * @param x input array. * @param y output array. */ public void applySmoothL(double kmax, float[][] x, float[][] y) { smoothL(kmax,x,y); }
LocalSmoothingFilter lsf = new LocalSmoothingFilter(0.01,10000,_ldk); lsf.setPreconditioner(true); float pavg = sum(p)/n1/n2; float[][] r = sub(p,pavg); lsf.applySmoothS(r,r); lsf.apply(_tensors,_c,s,r,q); add(q,pavg,q);
public void testSpd2() { int n1 = 5; int n2 = 6; LocalSmoothingFilter lsf = new LocalSmoothingFilter(1.0e-6,1000); lsf.setPreconditioner(true); for (int iter=0; iter<10; ++iter) { float[][] s = randfloat(n1,n2); float[][] x = sub(randfloat(n1,n2),0.5f); float[][] y = sub(randfloat(n1,n2),0.5f); float[][] dx = zerofloat(n1,n2); float[][] dy = zerofloat(n1,n2); Tensors2 d = new RandomTensors2(n1,n2); float c = 10.0f*s[n2/2][n1/2]; lsf.apply(d,c,s,x,dx); lsf.apply(d,c,s,y,dy); float xdx = dot(x,dx); float ydy = dot(y,dy); float ydx = dot(y,dx); float xdy = dot(x,dy); assertTrue(xdx>=0.0f); assertTrue(ydy>=0.0f); assertEquals(xdy,ydx,0.0001); } }
float[][] q = new float[n2][n1]; float[][] r = new float[n2][n1]; scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); // r = b-Ax scopy(r,d); // d = r float delta = sdot(r,r); // delta = r'r float bnorm = sqrt(sdot(b,b)); float rnorm = sqrt(delta); float rnormBegin = rnorm; log.finer(" iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin); float dq = sdot(d,q); // d'q = d'Ad float alpha = delta/dq; // alpha = r'r/d'Ad saxpy( alpha,d,x); // x = x+alpha*d saxpy(-alpha,q,r); // r = r-alpha*q float deltaOld = delta; delta = sdot(r,r); // delta = r'r float beta = delta/deltaOld; sxpay(beta,r,d); // d = r+beta*d rnorm = sqrt(delta);
public void apply( Direction2 d, EigenTensors2 t, float[][] f, float[][] g) { if (_scale==0.0f) { copy(f,g); } else { int n1 = f[0].length; int n2 = f.length; float[][] au = new float[n2][n1]; float[][] av = new float[n2][n1]; float[][] sf = new float[n2][n1]; t.getEigenvalues(au,av); setEigenvalues(d,t); _lsf.applySmoothL(_kmax,f,sf); //_lsf.applySmoothS(f,sf); _lsf.apply(t,_scale,sf,g); t.setEigenvalues(au,av); } } public void apply(
/** * Applies this filter for specified tensors and scale factors. * @param d tensors. * @param c constant scale factor for tensors. * @param s array of scale factors for tensors. * @param x input array. * @param y output array. */ public void apply( Tensors2 d, float c, float[][] s, float[][] x, float[][] y) { Operator2 a = new A2(_ldk,d,c,s); scopy(x,y); if (_pc) { Operator2 m = new M2(d,c,s,x); solve(a,m,x,y); } else { solve(a,x,y); } }
private void smoothL(double kmax, float[][][] x, float[][][] y) { ensureLowpassFilter(kmax); _lpf.apply(x,y); } private void ensureLowpassFilter(double kmax) {
LocalSmoothingFilter lsf = new LocalSmoothingFilter(0.01,10000,_ldk); lsf.setPreconditioner(true); float pavg = sum(p)/n1/n2/n3; float[][][] r = sub(p,pavg); lsf.applySmoothS(r,r); lsf.apply(_tensors,_c,s,r,q); add(q,pavg,q);
float[][][] q = new float[n3][n2][n1]; float[][][] r = new float[n3][n2][n1]; scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); // r = b-Ax scopy(r,d); float delta = sdot(r,r); float bnorm = sqrt(sdot(b,b)); float rnorm = sqrt(delta); float rnormBegin = rnorm; log.finer(" iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin); a.apply(d,q); float dq = sdot(d,q); float alpha = delta/dq; saxpy( alpha,d,x); if (iter%100<99) { saxpy(-alpha,q,r); } else { scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); delta = sdot(r,r); float beta = delta/deltaOld; sxpay(beta,r,d); rnorm = sqrt(delta);
public void apply( Direction3 d, EigenTensors3 t, float[][][] f, float[][][] g) { if (_scale==0.0f) { copy(f,g); } else { int n1 = f[0][0].length; int n2 = f[0].length; int n3 = f.length; float[][][] au = new float[n3][n2][n1]; float[][][] av = new float[n3][n2][n1]; float[][][] aw = new float[n3][n2][n1]; float[][][] sf = new float[n3][n2][n1]; t.getEigenvalues(au,av,aw); setEigenvalues(d,t); _lsf.applySmoothL(_kmax,f,sf); //_lsf.applySmoothS(f,sf); _lsf.apply(t,_scale,sf,g); t.setEigenvalues(au,av,aw); } } private float _scale;
/** * Applies this filter for specified tensors and scale factors. * @param d tensors. * @param c constant scale factor for tensors. * @param s array of scale factors for tensors. * @param x input array. * @param y output array. */ public void apply( Tensors3 d, float c, float[][][] s, float[][][] x, float[][][] y) { Operator3 a = new A3(_ldk,d,c,s); scopy(x,y); if (_pc) { Operator3 m = new M3(d,c,s,x); solve(a,m,x,y); } else { solve(a,x,y); } }
private void smoothL(double kmax, float[][] x, float[][] y) { ensureLowpassFilter(kmax); _lpf.apply(x,y); } private void smoothL(double kmax, float[][][] x, float[][][] y) {
float[][] r = new float[n2][n1]; float[][] s = new float[n2][n1]; scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); // r = b-Ax float bnorm = sqrt(sdot(b,b)); float rnorm = sqrt(sdot(r,r)); float rnormBegin = rnorm; float rnormSmall = bnorm*_small; scopy(s,d); // d = s float delta = sdot(r,s); // r's = r'Mr int iter; log.fine("msolve: bnorm="+bnorm+" rnorm="+rnorm); log.finer(" iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin); float alpha = delta/sdot(d,q); // alpha = r'Mr/d'Ad saxpy( alpha,d,x); // x = x+alpha*d saxpy(-alpha,q,r); // r = r-alpha*q delta = sdot(r,s); // delta = r's = r'Mr float beta = delta/deltaOld; sxpay(beta,s,d); // d = s+beta*d rnorm = sqrt(sdot(r,r));
/** * Applies this filter for specified tensors. * @param d tensors. * @param x input array. * @param y output array. */ public void apply(Tensors2 d, float[][] x, float[][] y) { apply(d,1.0f,null,x,y); }
/** * Applies an isotropic low-pass smoothing filter L. * Input and output arrays x and y may be the same array. * @param kmax maximum wavenumber not attenuated, in cycles/sample. * @param x input array. * @param y output array. */ public void applySmoothL(double kmax, float[][][] x, float[][][] y) { smoothL(kmax,x,y); }
float[][][] r = new float[n3][n2][n1]; float[][][] s = new float[n3][n2][n1]; scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); // r = b-Ax float bnorm = sqrt(sdot(b,b)); float rnorm = sqrt(sdot(r,r)); float rnormBegin = rnorm; float rnormSmall = bnorm*_small; scopy(s,d); // d = s float delta = sdot(r,s); // r's = r'Mr int iter; log.fine("msolve: bnorm="+bnorm+" rnorm="+rnorm); log.finer(" iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin); float alpha = delta/sdot(d,q); // alpha = r'Mr/d'Ad saxpy( alpha,d,x); // x = x+alpha*d if (iter%100<99) { saxpy(-alpha,q,r); // r = r-alpha*q } else { scopy(b,r); a.apply(x,q); saxpy(-1.0f,q,r); // r = b-Ax delta = sdot(r,s); // delta = r's = r'Mr float beta = delta/deltaOld; sxpay(beta,s,d); // d = s+beta*d rnorm = sqrt(sdot(r,r));
/** * Applies this filter for identity tensors and specified scale factors. * @param c constant scale factor. * @param s array of scale factors. * @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); }