summaryrefslogtreecommitdiffstats
path: root/tksao/frame/fitsanalysis.C
blob: 283b18575635771781431706e63f40e5914580a8 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
// Copyright (C) 1999-2018
// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
// For conditions of distribution and use, see copyright notice in "copyright"

#include <pthread.h>

#include "analysis.h"
#include "fitsimage.h"
#include "context.h"
#include "convolve.h"

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->xmin =0;
  targ->xmax =0;
  targ->ymin =0;
  targ->ymax =0;
  targ->width =0;
  targ->r =0;

  if (manageAnalysis_) {
    if (analysis_)
      delete analysis_;
    if (analysisdata_)
      delete analysisdata_;
  }
  manageAnalysis_ =0;
  analysis_ = block_;
  analysisdata_ = blockdata_;

  if (which) {
    analysis_ = new FitsAnalysis(block_, -64);
    if (analysis_->isValid()) {
      analysisdata_ = new FitsDatam<double>(analysis_, interp_);

      smooth(thread, targ);
      manageAnalysis_ =1;
    }
    else {
      delete analysis_;
      analysis_ = block_;
    }
  }

  image_ = analysis_;
  data_ = analysisdata_;
}

void* convolveThread(void* vv)
{
  t_smooth_arg* tt = (t_smooth_arg*)vv;
  convolve(tt->kernel, tt->src, tt->dest,
 	   tt->xmin, tt->ymin, tt->xmax, tt->ymax, tt->width, tt->r);
  return NULL;
}

void FitsImage::smooth(pthread_t* thread, t_smooth_arg* targ)
{
  FitsBound* params = getDataParams(context_->secMode());
  int width = analysis_->head()->naxis(0);
  int height = analysis_->head()->naxis(1);

  // src
  double* src = new double[width*height];
  double* ptr = src;
  for (long jj=0; jj<height; jj++)
    for (long ii=0; ii<width; ii++, ptr++)
      *ptr = blockdata_->getValueDouble(jj*width+ii);

  // dest
  double* dest = (double*)analysis_->data();

  // kernel
  double* kernel =NULL;
  switch (context_->smoothFunction()) {
  case Context::BOXCAR:
    kernel = boxcar(context_->smoothRadius());
    break;
  case Context::TOPHAT:
    kernel = tophat(context_->smoothRadius());
    break;
  case Context::GAUSSIAN:
    kernel = gaussian(context_->smoothRadius(), context_->smoothSigma());
    break;
  case Context::ELLIPTIC:
    kernel = elliptic(context_->smoothRadius(), context_->smoothRadiusMinor(),
		      context_->smoothSigma(), context_->smoothSigmaMinor(),
		      context_->smoothAngle());
    break;
  }

  // convolve
  targ->kernel = kernel;
  targ->src = src;
  targ->dest = dest;
  targ->xmin = params->xmin;
  targ->xmax = params->xmax;
  targ->ymin = params->ymin;
  targ->ymax = params->ymax;
  targ->width = width;
  targ->r = context_->smoothRadius();

  int result = pthread_create(thread, NULL, convolveThread, targ);
  if (result)
    internalError("Unable to Create Thread");
}