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

---
 tksao/frame/fitsanalysis.C |  3 ++-
 tksao/frame/fvcontour.C    | 55 +++++++++++++++++++++++++---------------------
 tksao/frame/fvcontour.h    |  1 -
 tksao/util/convolve.C      | 22 ++++++++++++++++++-
 tksao/util/convolve.h      |  2 +-
 5 files changed, 54 insertions(+), 29 deletions(-)

diff --git a/tksao/frame/fitsanalysis.C b/tksao/frame/fitsanalysis.C
index cad6cb7..679c4c0 100644
--- a/tksao/frame/fitsanalysis.C
+++ b/tksao/frame/fitsanalysis.C
@@ -52,7 +52,8 @@ void FitsImage::analysis(int which, pthread_t* thread, t_smooth_arg* targ)
 void* convolveThread(void* vv)
 {
   t_smooth_arg* tt = (t_smooth_arg*)vv;
-  convolve(tt->kernel, tt->src, tt->dest, tt->width, tt->height, tt->k);
+  convolve(tt->kernel, tt->src, tt->dest, 0, 0, tt->width, tt->height,
+	   tt->width, tt->k);
   return NULL;
 }
 
diff --git a/tksao/frame/fvcontour.C b/tksao/frame/fvcontour.C
index 5f39176..f9d8a6c 100644
--- a/tksao/frame/fvcontour.C
+++ b/tksao/frame/fvcontour.C
@@ -242,7 +242,36 @@ void FVContour::nobin(FitsImage* fits)
   }
 
   // convolve
-  convolve(fits,kernel,img,r);
+  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;
+
+    FitsBound* params = 
+      fits->getDataParams(((Base*)parent_)->currentContext->secMode());
+
+    SETSIGBUS
+    for(long jj=params->ymin; jj<params->ymax; jj++) {
+      for(long ii=params->xmin; ii<params->xmax; ii++) {
+	long kk = jj*width + ii;
+	double vv = fits->getValueDouble(kk);
+	if (isfinite(vv))
+	  src[kk] = vv;
+      }
+    }
+    CLEARSIGBUS
+
+    ::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);
@@ -289,30 +318,6 @@ void FVContour::convolve(FitsImage* fits, double* kernel, double* dest, int r)
   CLEARSIGBUS
 }
 
-double* FVContour::tophat(int r)
-{
-  int rr = 2*r+1;
-  int ksz = rr*rr;
-  double* kernel = new double[ksz];
-  memset(kernel, 0, ksz*sizeof(double));
-  
-  double kt = 0;
-  for (int yy=-r; yy<=r; yy++) {
-    for (int xx=-r; xx<=r; xx++) { 
-      if ((xx*xx + yy*yy) <= r*r) {
-	kernel[(yy+r)*rr+(xx+r)] = 1;
-	kt++;
-      }
-    }
-  }
-
-  // normalize kernel
-  for (int aa=0; aa<ksz; aa++)
-    kernel[aa] /= kt;
-
-  return kernel;
-}
-
 double* FVContour::gaussian(int r)
 {
   int rr = 2*r+1;
diff --git a/tksao/frame/fvcontour.h b/tksao/frame/fvcontour.h
index b547eb0..c8068b2 100644
--- a/tksao/frame/fvcontour.h
+++ b/tksao/frame/fvcontour.h
@@ -39,7 +39,6 @@ class FVContour {
   void bin(FitsImage*);
   void nobin(FitsImage*);
   void convolve(FitsImage*, double*, double*, int);
-  double* tophat(int);
   double* gaussian(int);
   void build(long xdim, long ydim, double *image, Matrix&);
   void trace(long xdim, long ydim, double cntr,
diff --git a/tksao/util/convolve.C b/tksao/util/convolve.C
index 51b5f15..3cec687 100644
--- a/tksao/util/convolve.C
+++ b/tksao/util/convolve.C
@@ -117,11 +117,13 @@ void dumpKernel(double* kernel, int k)
 }
 
 void* convolve(double* kernel, double* src, double* dest,
-	       int width, int height, int k)
+	       int xmin, int ymin, int xmax, int ymax, int width, int k)
 {
   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++) {
 
@@ -137,6 +139,24 @@ void* convolve(double* kernel, double* src, double* dest,
       }
     }
   }
+  }
+  else {
+  for (int jj=ymin; jj<ymax; jj++) {
+    for (int ii=xmin; ii<xmax; ii++) {
+      register int dd = jj*width+ii;
+      for (int nn=jj-k, qq=0; nn<=jj+k; nn++, qq++) {
+	if (nn>=ymin && nn<ymax) {
+	  register int nd = nn*width;
+	  register int qd = qq*kk;
+	  for (int mm=ii-k, pp=0; mm<=ii+k; mm++, pp++) {
+	    if (mm>=xmin && mm<xmax)
+	      dest[dd] += src[nd+mm]*kernel[qd+pp];
+	  }
+	}
+      }
+    }
+  }
+  }
 
   return NULL;
 }
diff --git a/tksao/util/convolve.h b/tksao/util/convolve.h
index 5aeb551..7465833 100644
--- a/tksao/util/convolve.h
+++ b/tksao/util/convolve.h
@@ -11,6 +11,6 @@ void gaussian(double* kernel, int r, double sigma);
 void elliptic(double* kernel, int r, int m, double ss, double sm, double aa);
 
 void* convolve(double* kernel, double* src, double* dest,
-	       int width, int height, int k);
+	       int xmin, int ymin, int xmax, int ymax, int width, int k);
 
 #endif
-- 
cgit v0.12