thread contour

2 files changed, 131 insertions, 0 deletions

diff --git a/tksao/util/convolve.C b/tksao/util/convolve.C
new file mode 100644
index 0000000..bdc6739
--- /dev/null
+++ b/tksao/util/convolve.C
@@ -0,0 +1,118 @@
+// Copyright (C) 1999-2019
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#include <math.h>
+
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+using namespace std;
+
+#include "convolve.h"
+
+void boxcar(double* kernel, int k)
+{
+  int kk = 2*k+1;
+  int kk2 = kk*kk;
+
+  for (int yy=-k; yy<=k; yy++)
+    for (int xx=-k; xx<=k; xx++)
+	kernel[(yy+k)*kk+(xx+k)] = 1;
+
+  // normalize kernel
+  for (int ii=0; ii<kk2; ii++)
+    kernel[ii] /= kk2;
+
+  //  dumpKernel(kernel, k);
+}
+
+void tophat(double* kernel, int k)
+{
+  int k2 = k*k;
+  int kk = 2*k+1;
+  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)/k2 <= 1) {
+	kernel[(yy+k)*kk+(xx+k)] = 1;
+	cnt++;
+      }
+
+  // normalize kernel
+  if (cnt) 
+    for (int ii=0; ii<kk2; ii++)
+      kernel[ii] /= cnt;
+
+  //  dumpKernel(kernel, k);
+}
+
+void gaussian(double* kernel, int k, double ss)
+{
+  int kk = 2*k+1;
+  int k2 = k*k;
+  int kk2 = kk*kk;
+  double s2 = ss*ss;
+
+  double tt =0;
+  for (int yy=-k; yy<=k; yy++)
+    for (int xx=-k; xx<=k; xx++) {
+      if ((xx*xx + yy*yy) <= k2) {
+	double vv = exp(-.5*((xx*xx + yy*yy)/s2));
+	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 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/convolve.h b/tksao/util/convolve.h
new file mode 100644
index 0000000..90b8f74
--- /dev/null
+++ b/tksao/util/convolve.h
@@ -0,0 +1,13 @@
+// Copyright (C) 1999-2019
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#ifndef __convolve_h__
+#define __convolve_h__
+
+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