add elliptical gaussian smooth

2 files changed, 64 insertions, 40 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;
+}
diff --git a/tksao/util/smooth.h b/tksao/util/smooth.h
index 0caf192..53e66dc 100644
--- a/tksao/util/smooth.h
+++ b/tksao/util/smooth.h
@@ -5,8 +5,9 @@
 #ifndef __smooth_h__
 #define __smooth_h__
 
-void boxcar(double* kernel, int k, int r);
-void tophat(double* kernel, int k, int r);
-void gaussian(double* kernel, int k, int r);
+void boxcar(double* kernel, int r);
+void tophat(double* kernel, int r);
+void gaussian(double* kernel, int r, double sigma);
+void elliptic(double* kernel, int r, int m, double ss, double sm, double aa);
 
 #endif