diff options
Diffstat (limited to 'tksao/util/smooth.C')
-rw-r--r-- | tksao/util/smooth.C | 97 |
1 files changed, 60 insertions, 37 deletions
diff --git a/tksao/util/smooth.C b/tksao/util/smooth.C index 4b7e49e..f87d9f2 100644 --- a/tksao/util/smooth.C +++ b/tksao/util/smooth.C @@ -11,76 +11,57 @@ using namespace std; #include "smooth.h" -void boxcar(double* kernel, int k, int r) -{ - if (r>k) - r=k; +static void dumpKernel(double* kernel, int k); +void boxcar(double* kernel, int k) +{ int kk = 2*k+1; - int ksz = kk*kk; + int kk2 = kk*kk; - int cnt =0; for (int yy=-k; yy<=k; yy++) for (int xx=-k; xx<=k; xx++) - if (abs(yy) <= r && abs(xx) <= r) { kernel[(yy+k)*kk+(xx+k)] = 1; - cnt++; - } // normalize kernel - if (cnt) - for (int ii=0; ii<ksz; ii++) - kernel[ii] /= cnt; + for (int ii=0; ii<kk2; ii++) + kernel[ii] /= kk2; - for (int yy=-k; yy<=k; yy++) - for (int xx=-k; xx<=k; xx++) - cerr << '(' << xx << ',' << yy << ")=" - << kernel[(yy+k)*kk+(xx+k)] << endl; + // dumpKernel(kernel, k); } -void tophat(double* kernel, int k, int r) +void tophat(double* kernel, int k) { - if (r>k) - r=k; - + int k2 = k*k; int kk = 2*k+1; - int ksz = kk*kk; - int rr = r*r; + int kk2 = kk*kk; int cnt =0; for (int yy=-k; yy<=k; yy++) for (int xx=-k; xx<=k; xx++) - if ((xx*xx + yy*yy) <= rr) { + if ((xx*xx + yy*yy)/k2 <= 1) { kernel[(yy+k)*kk+(xx+k)] = 1; cnt++; } // normalize kernel if (cnt) - for (int ii=0; ii<ksz; ii++) + for (int ii=0; ii<kk2; ii++) kernel[ii] /= cnt; - for (int yy=-k; yy<=k; yy++) - for (int xx=-k; xx<=k; xx++) - cerr << '(' << xx << ',' << yy << ")=" - << kernel[(yy+k)*kk+(xx+k)] << endl; + // dumpKernel(kernel, k); } -void gaussian(double* kernel, int k, int r) +void gaussian(double* kernel, int k, double ss) { - if (r>k) - r=k; - int k2 = k*k; int kk = 2*k+1; - int ksz = kk*kk; - double sigma = r/2.; - double s2 = sigma*sigma; + int kk2 = kk*kk; + double s2 = ss*ss; double total =0; for (int yy=-k; yy<=k; yy++) for (int xx=-k; xx<=k; xx++) - if ((xx*xx + yy*yy) <= k2) { + if ((xx*xx + yy*yy)/k2 <= 1) { double vv = exp(-.5*((xx*xx + yy*yy)/s2)); kernel[(yy+k)*kk+(xx+k)] = vv; total += vv; @@ -88,7 +69,49 @@ void gaussian(double* kernel, int k, int r) // normalize kernel if (total) - for (int ii=0; ii<ksz; ii++) + for (int ii=0; ii<kk2; ii++) kernel[ii] /= total; + + // dumpKernel(kernel, k); +} + +void elliptic(double* kernel, int k, int rm, double ss, double sm, double aa) +{ + int kk = 2*k+1; + int kk2 = kk*kk; + double s2 = ss*ss; + double sm2 = sm*sm; + + double a = cos(aa)*cos(aa)/(2*s2) + sin(aa)*sin(aa)/(2*sm2); + double b = -sin(2*aa)/(4*s2) + sin(2*aa)/(4*sm2); + double c = sin(aa)*sin(aa)/(2*s2) + cos(aa)*cos(aa)/(2*sm2); + + double tt =0; + for (int yy=-k; yy<=k; yy++) { + for (int xx=-k; xx<=k; xx++) { + double dd = xx*cos(aa)+yy*sin(aa); + double ee = xx*sin(aa)-yy*cos(aa); + if (dd*dd/(k*k) + ee*ee/(rm*rm)<= 1) { + double vv = exp(-(a*xx*xx + 2*b*xx*yy + c*yy*yy)); + kernel[(yy+k)*kk+(xx+k)] = vv; + tt += vv; + } + } + } + // normalize kernel + if (tt) + for (int ii=0; ii<kk2; ii++) + kernel[ii] /= tt; + + // dumpKernel(kernel, k); } +void dumpKernel(double* kernel, int k) +{ + int kk = 2*k+1; + + for (int yy=-k; yy<=k; yy++) + for (int xx=-k; xx<=k; xx++) + cerr << '(' << xx << ',' << yy << ")=" + << kernel[(yy+k)*kk+(xx+k)] << endl; +} |