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->width =0;
targ->height =0;
targ->k =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->width, tt->height, tt->k);
return NULL;
}
void FitsImage::smooth(pthread_t* thread, t_smooth_arg* targ)
{
int rr = context_->smoothRadius();
int mm = context_->smoothRadiusMinor();
double ss = context_->smoothSigma();
double sm = context_->smoothSigmaMinor();
double aa = context_->smoothAngle();
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 kk = 2*rr+1;
double* kernel = new double[kk*kk];
memset(kernel, 0, kk*kk*sizeof(double));
switch (context_->smoothFunction()) {
case Context::BOXCAR:
boxcar(kernel,rr);
break;
case Context::TOPHAT:
tophat(kernel,rr);
break;
case Context::GAUSSIAN:
gaussian(kernel,rr,ss);
break;
case Context::ELLIPTIC:
elliptic(kernel,rr,mm,ss,sm,aa);
break;
}
// convolve
targ->kernel = kernel;
targ->src = src;
targ->dest = dest;
targ->width = ww;
targ->height = hh;
targ->k = rr;
int result = pthread_create(thread, NULL, convolveThread, targ);
if (result)
internalError("Unable to Create Thread");
}
|