From fca9d92d00a9df9b4bb0facfe13e2e8a79ed0073 Mon Sep 17 00:00:00 2001
From: William Joye <wjoye@cfa.harvard.edu>
Date: Fri, 15 Mar 2019 14:25:22 -0400
Subject: thread contour

---
 tksao/frame/fitsanalysis.C |  14 ++----
 tksao/frame/fvcontour.C    | 112 +++++++--------------------------------------
 tksao/frame/fvcontour.h    |   2 -
 tksao/util/convolve.C      |  51 +++++++++------------
 tksao/util/convolve.h      |   8 ++--
 5 files changed, 46 insertions(+), 141 deletions(-)

diff --git a/tksao/frame/fitsanalysis.C b/tksao/frame/fitsanalysis.C
index 45e14ba..00b2f85 100644
--- a/tksao/frame/fitsanalysis.C
+++ b/tksao/frame/fitsanalysis.C
@@ -83,23 +83,19 @@ void FitsImage::smooth(pthread_t* thread, t_smooth_arg* targ)
   double* dest = (double*)analysis_->data();
 
   // kernel
-  // create kernel
-  int rr = 2*r+1;
-  double* kernel = new double[rr*rr];
-  memset(kernel, 0, rr*rr*sizeof(double));
-
+  double* kernel =NULL;
   switch (context_->smoothFunction()) {
   case Context::BOXCAR:
-    boxcar(kernel,r);
+    kernel = boxcar(r);
     break;
   case Context::TOPHAT:
-    tophat(kernel,r);
+    kernel = tophat(r);
     break;
   case Context::GAUSSIAN:
-    gaussian(kernel,r,ss);
+    kernel = gaussian(r,ss);
     break;
   case Context::ELLIPTIC:
-    elliptic(kernel,r,mm,ss,sm,aa);
+    kernel = elliptic(r,mm,ss,sm,aa);
     break;
   }
 
diff --git a/tksao/frame/fvcontour.C b/tksao/frame/fvcontour.C
index f9d8a6c..1f37ce9 100644
--- a/tksao/frame/fvcontour.C
+++ b/tksao/frame/fvcontour.C
@@ -228,35 +228,21 @@ void FVContour::nobin(FitsImage* fits)
 
   // generate kernel
   int r = smooth_-1;
-  double* kernel = NULL;
-
-  if (0)
-    kernel = gaussian(r);
-  else {
-    int rr = 2*r+1;
-    double sigma = r/2.;
-    int ksz = rr*rr;
-    kernel = new double[ksz];
-    memset(kernel, 0, ksz*sizeof(double));
-    ::gaussian(kernel, r, sigma);
-  }
+  double* kernel = ::gaussian(r, r/2.);
 
   // convolve
-  if (0)
-    convolve(fits,kernel,img,r);
-  else {
-    double* src = new double[size];
-    if (!img) {
-      internalError("FVContour could not allocate enough memory");
-      return;
-    }
-    for (long ii=0; ii<size; ii++)
-      src[ii] = FLT_MIN;
+  double* src = new double[size];
+  if (!img) {
+    internalError("FVContour could not allocate enough memory");
+    return;
+  }
+  for (long ii=0; ii<size; ii++)
+    src[ii] = FLT_MIN;
 
-    FitsBound* params = 
-      fits->getDataParams(((Base*)parent_)->currentContext->secMode());
+  FitsBound* params = 
+    fits->getDataParams(((Base*)parent_)->currentContext->secMode());
 
-    SETSIGBUS
+  SETSIGBUS
     for(long jj=params->ymin; jj<params->ymax; jj++) {
       for(long ii=params->xmin; ii<params->xmax; ii++) {
 	long kk = jj*width + ii;
@@ -265,13 +251,12 @@ void FVContour::nobin(FitsImage* fits)
 	  src[kk] = vv;
       }
     }
-    CLEARSIGBUS
+  CLEARSIGBUS
 
-    ::convolve(kernel, src, img,
-	       params->xmin, params->ymin, params->xmax, params->ymax,
-	       width, r);
-    delete [] src;
-  }
+  ::convolve(kernel, src, img,
+	     params->xmin, params->ymin, params->xmax, params->ymax,
+	     width, r);
+  delete [] src;
   
   // now, do contours
   build(width, height, img, fits->dataToRef);
@@ -281,71 +266,6 @@ void FVContour::nobin(FitsImage* fits)
   delete [] img;
 }
 
-void FVContour::convolve(FitsImage* fits, double* kernel, double* dest, int r)
-{
-  FitsBound* params = 
-    fits->getDataParams(((Base*)parent_)->currentContext->secMode());
-  long width = fits->width();
-  int rr = 2*r+1;
-
-  SETSIGBUS
-  for (long jj=params->ymin; jj<params->ymax; jj++) {
-    for (long ii=params->xmin; ii<params->xmax; ii++) {
-      long ir  = ii-r;
-      long irr = ii+r;
-      long jr = jj-r;
-      long jrr = jj+r;
-
-      for (long n=jr, nn=0; n<=jrr; n++, nn++) {
-	if (n>=params->ymin && n<params->ymax) {
-	  for (long m=ir, mm=0; m<=irr; m++, mm++) {
-	    if (m>=params->xmin && m<params->xmax) {
-	      double vv = fits->getValueDouble(n*width+m);
-	      if (isfinite(vv)) {
-		double kk = kernel[nn*rr+mm];
-		double* ptr = dest+(jj*width+ii);
-		if (*ptr == FLT_MIN)
-		  *ptr  = vv*kk;
-		else
-		  *ptr += vv*kk;
-	      }
-	    }
-	  }
-	}
-      }
-    }
-  }
-  CLEARSIGBUS
-}
-
-double* FVContour::gaussian(int r)
-{
-  int rr = 2*r+1;
-  int ksz = rr*rr;
-  double sigma = r/2.;
-  double* kernel = new double[ksz];
-  memset(kernel, 0, ksz*sizeof(double));
-  
-  double kt = 0;
-  double aa = 1./(sigma*sigma);
-  double cc = 1./(sigma*sigma);
-  for (int yy=-r; yy<=r; yy++) {
-    for (int xx=-r; xx<=r; xx++) { 
-      if ((xx*xx + yy*yy) <= r*r) {
-	double vv = exp(-.5*(aa*xx*xx + cc*yy*yy));
-	kernel[(yy+r)*rr+(xx+r)] = vv;
-	kt += vv;
-      }
-    }
-  }
-
-  // normalize kernel
-  for (int aa=0; aa<ksz; aa++)
-    kernel[aa] /= kt;
-
-  return kernel;
-}
-
 void FVContour::bin(FitsImage* fits)
 {
   FitsBound* params = 
diff --git a/tksao/frame/fvcontour.h b/tksao/frame/fvcontour.h
index c8068b2..79f4384 100644
--- a/tksao/frame/fvcontour.h
+++ b/tksao/frame/fvcontour.h
@@ -38,8 +38,6 @@ class FVContour {
   void unity(FitsImage*);
   void bin(FitsImage*);
   void nobin(FitsImage*);
-  void convolve(FitsImage*, double*, double*, int);
-  double* gaussian(int);
   void build(long xdim, long ydim, double *image, Matrix&);
   void trace(long xdim, long ydim, double cntr,
 	     long xCell, long yCell, int side,
diff --git a/tksao/util/convolve.C b/tksao/util/convolve.C
index 3cec687..d09a69e 100644
--- a/tksao/util/convolve.C
+++ b/tksao/util/convolve.C
@@ -11,11 +11,14 @@ using namespace std;
 
 #include "convolve.h"
 
-void boxcar(double* kernel, int k)
+double* boxcar(int k)
 {
   int kk = 2*k+1;
   int kk2 = kk*kk;
 
+  double* kernel = new double[kk*kk];
+  memset(kernel, 0, kk*kk*sizeof(double));
+
   for (int yy=-k; yy<=k; yy++)
     for (int xx=-k; xx<=k; xx++)
 	kernel[(yy+k)*kk+(xx+k)] = 1;
@@ -24,15 +27,18 @@ void boxcar(double* kernel, int k)
   for (int ii=0; ii<kk2; ii++)
     kernel[ii] /= kk2;
 
-  //  dumpKernel(kernel, k);
+  return kernel;
 }
 
-void tophat(double* kernel, int k)
+double* tophat(int k)
 {
-  int k2 = k*k;
   int kk = 2*k+1;
+  int k2 = k*k;
   int kk2 = kk*kk;
 
+  double* kernel = new double[kk*kk];
+  memset(kernel, 0, kk*kk*sizeof(double));
+
   int cnt =0;
   for (int yy=-k; yy<=k; yy++)
     for (int xx=-k; xx<=k; xx++)
@@ -46,16 +52,19 @@ void tophat(double* kernel, int k)
     for (int ii=0; ii<kk2; ii++)
       kernel[ii] /= cnt;
 
-  //  dumpKernel(kernel, k);
+  return kernel;
 }
 
-void gaussian(double* kernel, int k, double ss)
+double* gaussian(int k, double ss)
 {
   int kk = 2*k+1;
   int k2 = k*k;
   int kk2 = kk*kk;
   double s2 = ss*ss;
 
+  double* kernel = new double[kk*kk];
+  memset(kernel, 0, kk*kk*sizeof(double));
+
   double tt =0;
   for (int yy=-k; yy<=k; yy++)
     for (int xx=-k; xx<=k; xx++) {
@@ -71,16 +80,19 @@ void gaussian(double* kernel, int k, double ss)
     for (int ii=0; ii<kk2; ii++)
       kernel[ii] /= tt;
 
-  // dumpKernel(kernel, k);
+  return kernel;
 }
 
-void elliptic(double* kernel, int k, int rm, double ss, double sm, double aa)
+double* elliptic(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* kernel = new double[kk*kk];
+  memset(kernel, 0, kk*kk*sizeof(double));
+
   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);
@@ -103,7 +115,7 @@ void elliptic(double* kernel, int k, int rm, double ss, double sm, double aa)
     for (int ii=0; ii<kk2; ii++)
       kernel[ii] /= tt;
 
-  // dumpKernel(kernel, k);
+  return kernel;
 }
 
 void dumpKernel(double* kernel, int k)
@@ -121,26 +133,6 @@ void* convolve(double* kernel, double* src, double* dest,
 {
   int kk = 2*k+1;
 
-  if (0) {
-  double* dptr = dest;
-  int height = ymax;
-  for (int jj=0; jj<height; jj++) {
-    for (int ii=0; ii<width; ii++, dptr++) {
-
-      for (int nn=jj-k, qq=0; nn<=jj+k; nn++, qq++) {
-	if (nn>=0 && nn<height) {
-	  register int nd = nn*width;
-	  register int qd = qq*kk;
-	  for (int mm=ii-k, pp=0; mm<=ii+k; mm++, pp++) {
-	    if (mm>=0 && mm<width)
-	      *dptr += src[nd+mm]*kernel[qd+pp];
-	  }
-	}
-      }
-    }
-  }
-  }
-  else {
   for (int jj=ymin; jj<ymax; jj++) {
     for (int ii=xmin; ii<xmax; ii++) {
       register int dd = jj*width+ii;
@@ -156,7 +148,6 @@ void* convolve(double* kernel, double* src, double* dest,
       }
     }
   }
-  }
 
   return NULL;
 }
diff --git a/tksao/util/convolve.h b/tksao/util/convolve.h
index 7465833..637cf97 100644
--- a/tksao/util/convolve.h
+++ b/tksao/util/convolve.h
@@ -5,10 +5,10 @@
 #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);
+double* boxcar(int r);
+double* tophat(int r);
+double* gaussian(int r, double sigma);
+double* elliptic(int r, int m, double ss, double sm, double aa);
 
 void* convolve(double* kernel, double* src, double* dest,
 	       int xmin, int ymin, int xmax, int ymax, int width, int k);
-- 
cgit v0.12