Initial commit

1 files changed, 575 insertions, 0 deletions

diff --git a/tksao/frame/fvcontour.C b/tksao/frame/fvcontour.C
new file mode 100644
index 0000000..8b049dc
--- /dev/null
+++ b/tksao/frame/fvcontour.C
@@ -0,0 +1,575 @@
+// Copyright (C) 1999-2016
+// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+// For conditions of distribution and use, see copyright notice in "copyright"
+
+#include "context.h"
+#include "base.h"
+#include "fitsimage.h"
+#include "fvcontour.h"
+
+#include "sigbus.h"
+
+#ifdef INTERP
+#undef INTERP
+#endif
+#define INTERP (((Base*)parent_)->interp)
+
+static const char* methodName_[] = { 
+    "smooth",
+    "block"
+  };
+
+// It is a modified version of contour code found in Fv 2.4
+// Fv may be obtained from the HEASARC (High Energy Astrophysics Science
+// Archive Research Center) FTOOLS Web site at:
+// http://heasarc.gsfc.nasa.gov/ftools/fv.html
+
+// The original author is unknown.
+
+FVContour::FVContour()
+{
+  parent_ =NULL;
+
+  colorName_ = dupstr("green");
+  lineWidth_ =1;
+  dash_ =0;
+  dlist_[0] =8;
+  dlist_[1] =3;
+
+  method_ = BLOCK;
+  smooth_ =4;
+  numLevel_ =5;
+
+  colorScaleType_ = FrScale::LINEARSCALE;
+  clipMode_ = (float)FrScale::MINMAX;
+  expo_ =1000;
+  limits_ = Vector(0,100);
+
+  level_ =NULL;
+  scale_ =NULL;
+}
+
+FVContour::~FVContour()
+{
+  if (colorName_)
+    delete [] colorName_;
+
+  if (level_)
+    delete [] level_;
+
+  if (scale_)
+    delete scale_;
+}
+
+void FVContour::create(Base* pp, FitsImage* fits, FrScale* fr,
+		       const char* cc, int ww, int dd, 
+		       Method mm, int nn, int rr, 
+		       const char* ll,
+		       FrScale::ColorScaleType sc, float exp, 
+		       float cm, Vector lim)
+{
+  lcontourlevel_.deleteAll();
+
+  parent_ = pp;
+
+  colorName_ = dupstr(cc);
+  lineWidth_ = ww;
+  dash_ = dd;
+
+  method_ = mm;
+  smooth_ = rr;
+  numLevel_ = nn;
+
+  colorScaleType_ = sc;
+  clipMode_ = cm;
+  expo_ = exp;
+  limits_ = lim;
+
+  level_ = dupstr(ll);
+
+  if (level_ && strlen(level_)>0) {
+    int cnt = 0;
+    double levels[100];
+    string x(level_);
+    istringstream str(x);
+
+    while ((cnt<100) && (str >> levels[cnt]))
+      cnt++;
+
+    scale_ = new InverseScale(cnt, levels);
+  }
+  else
+    buildScale(fits, fr);
+
+  append(fits);
+}
+
+void FVContour::buildScale(FitsImage* fits, FrScale* fr)
+{
+  switch (colorScaleType_) {
+  case FrScale::LINEARSCALE:
+    scale_ = new LinearInverseScale(numLevel_, limits_[0], limits_[1]);
+    break;
+  case FrScale::LOGSCALE:
+    scale_ =  new LogInverseScale(numLevel_, limits_[0], limits_[1], expo_);
+    break;
+  case FrScale::POWSCALE:
+    scale_ =  new PowInverseScale(numLevel_, limits_[0], limits_[1], expo_);
+    break;
+  case FrScale::SQRTSCALE:
+    scale_ =  new SqrtInverseScale(numLevel_, limits_[0], limits_[1]);
+    break;
+  case FrScale::SQUAREDSCALE:
+    scale_ =  new SquaredInverseScale(numLevel_, limits_[0], limits_[1]);
+    break;
+  case FrScale::ASINHSCALE:
+    scale_ =  new AsinhInverseScale(numLevel_, limits_[0], limits_[1]);
+    break;
+  case FrScale::SINHSCALE:
+    scale_ =  new SinhInverseScale(numLevel_, limits_[0], limits_[1]);
+    break;
+  case FrScale::HISTEQUSCALE:
+    scale_ =  new HistEquInverseScale(numLevel_, limits_[0], limits_[1], 
+				      fr->histequ(fits), HISTEQUSIZE);
+    break;
+  case FrScale::IISSCALE:
+    scale_ = new IISInverseScale(numLevel_, limits_[0], limits_[1],
+				 fits->iisz());
+    break;
+  }
+}
+
+void FVContour::update(FitsImage* fits)
+{
+  if (lcontourlevel_.isEmpty()) 
+    return;
+
+  lcontourlevel_.deleteAll();
+  append(fits);
+}
+
+void FVContour::update(FitsImage* fits, FrScale* fr)
+{
+  if (lcontourlevel_.isEmpty())
+    return;
+
+  lcontourlevel_.deleteAll();
+
+  if (scale_)
+    delete scale_;
+  limits_ = Vector(fr->low(),fr->high());
+  expo_ = fr->expo();
+
+  buildScale(fits, fr);
+
+  if (level_)
+    delete [] level_;
+  ostringstream str;
+  str << *scale_ << ends;
+  level_ = dupstr(str.str().c_str());
+
+  append(fits);
+}
+
+const char* FVContour::methodName()
+{
+  return methodName_[method_];
+}
+
+void FVContour::append(FitsImage* fits)
+{
+  if (smooth_ == 1)
+    unity(fits);
+  else
+    switch (method_) {
+    case SMOOTH:
+      nobin(fits);
+      break;
+    case BLOCK:
+      bin(fits);
+      break;
+    }
+}
+
+void FVContour::unity(FitsImage* fits)
+{
+  FitsBound* params = 
+    fits->getDataParams(((Base*)parent_)->currentContext->secMode());
+  long width = fits->width();
+  long height = fits->height();
+
+  // blank img
+  long size = width*height;
+  double* img = new double[size];
+  if (!img) {
+    internalError("FVContour could not allocate enough memory");
+    return;
+  }
+  for (long ii=0; ii<size; ii++)
+    img[ii] = FLT_MIN;
+
+  // fill img
+  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))
+	img[kk] = vv;
+    }
+  }
+  CLEARSIGBUS
+
+  // do contours
+  build(width, height, img, fits->dataToRef);
+
+  // clean up
+  delete [] img;
+}
+
+void FVContour::nobin(FitsImage* fits)
+{
+  long width = fits->width();
+  long height = fits->height();
+
+  // blank img
+  long size = width*height;
+  double* img = new double[size];
+  if (!img) {
+    internalError("FVContour could not allocate enough memory");
+    return;
+  }
+  for (long ii=0; ii<size; ii++)
+    img[ii] = FLT_MIN;
+
+  // generate kernal
+  int r = smooth_-1;
+  double* kernal = gaussian(r);
+
+  // convolve
+  convolve(fits,kernal,img,r);
+  
+  // now, do contours
+  build(width, height, img, fits->dataToRef);
+
+  // cleanup
+  delete kernal;
+  delete [] img;
+}
+
+void FVContour::convolve(FitsImage* fits, double* kernal, 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 = kernal[nn*rr+mm];
+		double* ptr = dest+(jj*width+ii);
+		if (*ptr == FLT_MIN)
+		  *ptr  = vv*kk;
+		else
+		  *ptr += vv*kk;
+	      }
+	    }
+	  }
+	}
+      }
+    }
+  }
+  CLEARSIGBUS
+}
+
+double* FVContour::tophat(int r)
+{
+  int rr = 2*r+1;
+  int ksz = rr*rr;
+  double* kernal = new double[ksz];
+  memset(kernal, 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) {
+	kernal[(yy+r)*rr+(xx+r)] = 1;
+	kt++;
+      }
+    }
+  }
+
+  // normalize kernal
+  for (int aa=0; aa<ksz; aa++)
+    kernal[aa] /= kt;
+
+  return kernal;
+}
+
+double* FVContour::gaussian(int r)
+{
+  int rr = 2*r+1;
+  int ksz = rr*rr;
+  double sigma = r/2.;
+  double* kernal = new double[ksz];
+  memset(kernal, 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));
+	kernal[(yy+r)*rr+(xx+r)] = vv;
+	kt += vv;
+      }
+    }
+  }
+
+  // normalize kernal
+  for (int aa=0; aa<ksz; aa++)
+    kernal[aa] /= kt;
+
+  return kernal;
+}
+
+void FVContour::bin(FitsImage* fits)
+{
+  FitsBound* params = 
+    fits->getDataParams(((Base*)parent_)->currentContext->secMode());
+  long width = fits->width();
+  long height = fits->height();
+
+  int rr = smooth_;
+
+  long w2 = (long)(width/rr);
+  long h2 = (long)(height/rr);
+
+  Matrix m = 
+    Translate((Vector(-width,-height)/2).floor()) * 
+    Scale(1./rr) * 
+    Translate((Vector(w2,h2)/2).floor());
+  Matrix n = m.invert();
+  double* mm = m.mm();
+
+  double* img = new double[w2 * h2];
+  {
+    for (long jj=0; jj<h2; jj++)
+      for (long ii=0; ii<w2; ii++)
+	img[jj*w2 + ii] = FLT_MIN;
+  }
+
+  short* count = new short[w2 * h2];
+  memset(count, 0, w2*h2*sizeof(short));
+
+  SETSIGBUS
+  for (long jj=params->ymin; jj<params->ymax; jj++) {
+    for (long ii=params->xmin; ii<params->xmax; ii++) {
+      double xx = ii*mm[0] + jj*mm[3] + mm[6];
+      double yy = ii*mm[1] + jj*mm[4] + mm[7];
+
+      if (xx >= 0 && xx < w2 && yy >= 0 && yy < h2) {
+	long kk = (long(yy)*w2 + long(xx));
+	double v = fits->getValueDouble(jj*width + ii);
+	if (isfinite(v)) {
+	  if (count[kk])
+	    img[kk] += v;
+	  else
+	    img[kk] = v;
+
+	  count[kk]++;
+	}
+      }
+    }
+  }
+
+  for (long kk=0; kk<w2*h2; kk++)
+    if (count[kk])
+      img[kk] /= count[kk];
+
+  CLEARSIGBUS
+
+  delete [] count;
+
+  Matrix w = n * fits->dataToRef;
+  build(w2, h2, img, w);
+
+  delete [] img;
+}
+
+void FVContour::build(long xdim, long ydim, double *image, Matrix& mx)
+{
+  long nelem = xdim*ydim;
+  char* usedGrid = new char[nelem];
+  double** rows = new double*[ydim];
+
+  for (long jj=0; jj<ydim; jj++)
+    rows[jj] = image + jj*xdim;
+
+  for (long c=0; c<scale_->size(); c++) {
+    double cntour = scale_->level(c);
+
+    ContourLevel* cl =new ContourLevel(parent_, cntour, colorName_, lineWidth_, 
+				       dash_, dlist_);
+    memset(usedGrid,0,nelem);
+
+    //  Search outer edge
+    long ii,jj;
+
+    //  Search top
+    for (jj=0, ii=0; ii<xdim-1; ii++)
+      if (rows[jj][ii]<cntour && cntour<=rows[jj][ii+1])
+	trace(xdim, ydim, cntour, ii, jj, top, rows, usedGrid, mx, cl);
+
+    //  Search right
+    for (jj=0; jj<ydim-1; jj++)
+      if (rows[jj][ii]<cntour && cntour<=rows[jj+1][ii])
+	trace(xdim, ydim, cntour, ii-1, jj, right, rows, usedGrid, mx, cl);
+
+    //  Search Bottom
+    for (ii--; ii>=0; ii--)
+      if (rows[jj][ii+1]<cntour && cntour<=rows[jj][ii])
+	trace(xdim, ydim, cntour, ii, jj-1, bottom, rows, usedGrid, mx, cl);
+
+    //  Search Left
+    for (ii=0, jj--; jj>=0; jj--)
+      if (rows[jj+1][ii]<cntour && cntour<=rows[jj][ii])
+	trace(xdim, ydim, cntour, ii, jj, left, rows, usedGrid, mx, cl);
+
+    //  Search each row of the image
+    for (jj=1; jj<ydim-1; jj++)
+      for (ii=0; ii<xdim-1; ii++)
+	if (!usedGrid[jj*xdim + ii] && 
+	    rows[jj][ii]<cntour && 
+	    cntour<=rows[jj][ii+1])
+	  trace(xdim, ydim, cntour, ii, jj, top, rows, usedGrid, mx, cl);
+
+    if (!cl->lcontour().isEmpty())
+      lcontourlevel_.append(cl);
+  }
+
+  delete [] usedGrid;
+  delete [] rows;
+}
+
+void FVContour::trace(long xdim, long ydim, double cntr,
+		      long xCell, long yCell, int side, 
+		      double** rows, char* usedGrid, 
+		      Matrix& mx, ContourLevel* cl)
+{
+  long ii = xCell;
+  long jj = yCell;
+  int origSide = side;
+
+  int init = 1;
+  int done = (ii<0 || ii>=xdim-1 || (jj<0 && jj>=ydim-1));
+
+  Contour* cc = new Contour(cl);
+
+  while (!done) {
+    int flag = 0;
+    double a = rows[jj][ii];
+    double b = rows[jj][ii+1];
+    double c = rows[jj+1][ii+1];
+    double d = rows[jj+1][ii];
+
+    double X, Y;
+    if (init) {
+      init = 0;
+      switch (side) {
+      case top:
+	X = (cntr-a) / (b-a) + ii;
+	Y = jj;
+	break;
+      case right:
+	X = ii+1;
+	Y = (cntr-b) / (c-b) + jj;
+	break;
+      case bottom:
+	X = (cntr-c) / (d-c) + ii;
+	Y = jj+1;
+	break;
+      case left:
+	X = ii;
+	Y = (cntr-a) / (d-a) + jj;
+	break;
+      }
+
+    }
+    else {
+      if (side==top)
+	usedGrid[jj*xdim + ii] = 1;
+
+      do {
+	if (++side == none)
+	  side = top;
+
+	switch (side) {
+	case top:
+	  if (a>=cntr && cntr>b) {
+	    flag = 1;
+	    X = (cntr-a) / (b-a) + ii;
+	    Y = jj;
+	    jj--;
+	  }
+	  break;
+	case right:
+	  if( b>=cntr && cntr>c ) {
+	    flag = 1;
+	    X = ii+1;
+	    Y = (cntr-b) / (c-b) + jj;
+	    ii++;
+	  }
+	  break;
+	case bottom:
+	  if( c>=cntr && cntr>d ) {
+	    flag = 1;
+	    X = (cntr-d) / (c-d) + ii;
+	    Y = jj+1;
+	    jj++;
+	  }
+	  break;
+	case left:
+	  if( d>=cntr && cntr>a ) {
+	    flag = 1;
+	    X = ii;
+	    Y = (cntr-a) / (d-a) + jj;
+	    ii--;
+	  }
+	  break;
+	}
+      } while (!flag);
+
+      if (++side == none)
+	side = top;
+      if (++side == none)
+	side = top;
+      if (ii==xCell && jj==yCell && side==origSide)
+	done = 1;
+      if (ii<0 || ii>=xdim-1 || jj<0 || jj>=ydim-1)
+	done = 1;
+    }
+
+    cc->lvertex().append(new Vertex(Vector(X+.5,Y+.5)*mx));
+  }
+
+  if (!cc->lvertex().isEmpty())
+    cl->lcontour().append(cc);
+  else
+    delete cc;
+}
+