diff options
Diffstat (limited to 'tksao')
-rw-r--r-- | tksao/frame/fitsanalysis.C | 14 | ||||
-rw-r--r-- | tksao/frame/fvcontour.C | 112 | ||||
-rw-r--r-- | tksao/frame/fvcontour.h | 2 | ||||
-rw-r--r-- | tksao/util/convolve.C | 51 | ||||
-rw-r--r-- | 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); |