diff options
Diffstat (limited to 'tksao/frame/fitsanalysis.C')
| -rw-r--r-- | tksao/frame/fitsanalysis.C | 231 |
1 files changed, 231 insertions, 0 deletions
diff --git a/tksao/frame/fitsanalysis.C b/tksao/frame/fitsanalysis.C new file mode 100644 index 0000000..d3bc830 --- /dev/null +++ b/tksao/frame/fitsanalysis.C @@ -0,0 +1,231 @@ +// Copyright (C) 1999-2016 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#ifndef __WIN32 +#include <pthread.h> +#endif + +#include "fitsimage.h" +#include "analysis.h" +#include "smooth.h" +#include "context.h" + +void* convolve(void* tt); + +#ifdef __WIN32 + +void FitsImage::analysis(int which) +{ + if (DebugPerf) + cerr << "FitsImage::analysis()" << endl; + + if (manageAnalysis_) { + if (analysis_) + delete analysis_; + if (analysisdata_) + delete analysisdata_; + } + manageAnalysis_ =0; + analysis_ = block_; + analysisdata_ = blockdata_; + + if (which) { + analysis_ = new FitsAnalysis(block_); + if (analysis_->isValid()) { + analysisdata_ = new FitsDatam<double>(analysis_, interp_); + + smooth(); + manageAnalysis_ =1; + } + else { + delete analysis_; + analysis_ = block_; + } + } + + image_ = analysis_; + data_ = analysisdata_; +} + +void FitsImage::smooth() +{ + // radius + int r = context_->smoothRadius(); + + int ww = analysis_->head()->naxis(0); + int hh = analysis_->head()->naxis(1); + + // src + double* src = new double[ww*hh]; + double* ptr = src; + for (long jj=0; jj<hh; jj++) + for (long ii=0; ii<ww; ii++, ptr++) + *ptr = blockdata_->getValueDouble(jj*ww+ii); + + // dest + 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)); + + switch (context_->smoothFunction()) { + case Context::BOXCAR: + boxcar(kernel,r); + break; + case Context::TOPHAT: + tophat(kernel,r); + break; + case Context::GAUSSIAN: + gaussian(kernel,r); + break; + } + + // convolve + t_smooth_arg* targ = new t_smooth_arg; + targ->kernel = kernel; + targ->src = src; + targ->dest = dest; + targ->width = ww; + targ->height = hh; + targ->radius = r; + + convolve(targ); + + // clean up + if (targ->kernel) + delete [] targ->kernel; + if (targ->src) + delete [] targ->src; + if (targ) + delete targ; +} + +#else + +void FitsImage::analysis(int which, pthread_t* thread, t_smooth_arg* targ) +{ + if (DebugPerf) + cerr << "FitsImage::analysis()" << endl; + + targ->kernel =NULL; + targ->src =NULL; + targ->dest =NULL; + targ->width =0; + targ->height =0; + targ->radius =0; + + if (manageAnalysis_) { + if (analysis_) + delete analysis_; + if (analysisdata_) + delete analysisdata_; + } + manageAnalysis_ =0; + analysis_ = block_; + analysisdata_ = blockdata_; + + if (which) { + analysis_ = new FitsAnalysis(block_); + if (analysis_->isValid()) { + analysisdata_ = new FitsDatam<double>(analysis_, interp_); + + smooth(thread, targ); + manageAnalysis_ =1; + } + else { + delete analysis_; + analysis_ = block_; + } + } + + image_ = analysis_; + data_ = analysisdata_; +} + +void FitsImage::smooth(pthread_t* thread, t_smooth_arg* targ) +{ + // radius + int r = context_->smoothRadius(); + + int ww = analysis_->head()->naxis(0); + int hh = analysis_->head()->naxis(1); + + // src + double* src = new double[ww*hh]; + double* ptr = src; + for (long jj=0; jj<hh; jj++) + for (long ii=0; ii<ww; ii++, ptr++) + *ptr = blockdata_->getValueDouble(jj*ww+ii); + + // dest + 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)); + + switch (context_->smoothFunction()) { + case Context::BOXCAR: + boxcar(kernel,r); + break; + case Context::TOPHAT: + tophat(kernel,r); + break; + case Context::GAUSSIAN: + gaussian(kernel,r); + break; + } + + // convolve + targ->kernel = kernel; + targ->src = src; + targ->dest = dest; + targ->width = ww; + targ->height = hh; + targ->radius = r; + + int result = pthread_create(thread, NULL, convolve, targ); + if (result) + internalError("Unable to Create Thread"); +} + +#endif + +void* convolve(void* tt) +{ + t_smooth_arg* targ = (t_smooth_arg*)tt; + double* kernel = targ->kernel; + double* src = targ->src; + double* dest = targ->dest; + int width = targ->width; + int height = targ->height; + int r = targ->radius; + + int rr = 2*r+1; + + double* dptr = dest; + for (int jj=0; jj<height; jj++) { + for (int ii=0; ii<width; ii++, dptr++) { + + for (int nn=jj-r, qq=0; nn<=jj+r; nn++, qq++) { + if (nn>=0 && nn<height) { + register int nd = nn*width; + register int qd = qq*rr; + for (int mm=ii-r, pp=0; mm<=ii+r; mm++, pp++) { + if (mm>=0 && mm<width) + *dptr += src[nd+mm]*kernel[qd+pp]; + } + } + } + } + } + + return NULL; +} + |