From f11ec933aa77ef7b53e584b4878438cd30a73bc4 Mon Sep 17 00:00:00 2001
From: William Joye <wjoye@cfa.harvard.edu>
Date: Thu, 14 Mar 2019 16:06:50 -0400
Subject: thread contour

---
 tksao/configure            |   1 +
 tksao/configure.ac         |   1 +
 tksao/frame/fitsanalysis.C | 112 +-----------------------------------------
 tksao/util/convolve.C      | 118 +++++++++++++++++++++++++++++++++++++++++++++
 tksao/util/convolve.h      |  13 +++++
 5 files changed, 135 insertions(+), 110 deletions(-)
 create mode 100644 tksao/util/convolve.C
 create mode 100644 tksao/util/convolve.h

diff --git a/tksao/configure b/tksao/configure
index d9d5e50..59f3c04 100755
--- a/tksao/configure
+++ b/tksao/configure
@@ -5349,6 +5349,7 @@ panner/panner.C
 panner/pannertrue.C
 panner/parser.C
 util/attribute.C
+util/convolve.C
 util/grf.C
 util/grf3d.C
 util/grid25dbase.C
diff --git a/tksao/configure.ac b/tksao/configure.ac
index a4735e4..507f7dd 100644
--- a/tksao/configure.ac
+++ b/tksao/configure.ac
@@ -252,6 +252,7 @@ panner/panner.C
 panner/pannertrue.C
 panner/parser.C
 util/attribute.C
+util/convolve.C
 util/grf.C
 util/grf3d.C
 util/grid25dbase.C
diff --git a/tksao/frame/fitsanalysis.C b/tksao/frame/fitsanalysis.C
index 99e1e58..eb823b4 100644
--- a/tksao/frame/fitsanalysis.C
+++ b/tksao/frame/fitsanalysis.C
@@ -4,14 +4,11 @@
 
 #include <pthread.h>
 
-#include "fitsimage.h"
 #include "analysis.h"
+#include "fitsimage.h"
 #include "context.h"
+#include "convolve.h"
 
-static void boxcar(double* kernel, int r);
-static void tophat(double* kernel, int r);
-static void gaussian(double* kernel, int r, double sigma);
-static void elliptic(double* kernel, int r, int m, double ss, double sm, double aa);
 static void* convolve(void* tt);
 
 void FitsImage::analysis(int which, pthread_t* thread, t_smooth_arg* targ)
@@ -109,111 +106,6 @@ void FitsImage::smooth(pthread_t* thread, t_smooth_arg* targ)
     internalError("Unable to Create Thread");
 }
 
-static 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 k2 = k*k;
-  int kk = 2*k+1;
-  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;
-}
-
 void* convolve(void* tt)
 {
   t_smooth_arg* targ = (t_smooth_arg*)tt;
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
-- 
cgit v0.12