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
114
115
116
117
118
|
// Copyright (C) 1999-2019
// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
// For conditions of distribution and use, see copyright notice in "copyright"
#include <math.h>
#include <iostream>
#include <sstream>
#include <iomanip>
using namespace std;
#include "convolve.h"
void boxcar(double* kernel, int k)
{
int kk = 2*k+1;
int kk2 = kk*kk;
for (int yy=-k; yy<=k; yy++)
for (int xx=-k; xx<=k; xx++)
kernel[(yy+k)*kk+(xx+k)] = 1;
// normalize kernel
for (int ii=0; ii<kk2; ii++)
kernel[ii] /= kk2;
// dumpKernel(kernel, k);
}
void tophat(double* kernel, int k)
{
int k2 = k*k;
int kk = 2*k+1;
int kk2 = kk*kk;
int cnt =0;
for (int yy=-k; yy<=k; yy++)
for (int xx=-k; xx<=k; xx++)
if ((xx*xx + yy*yy)/k2 <= 1) {
kernel[(yy+k)*kk+(xx+k)] = 1;
cnt++;
}
// normalize kernel
if (cnt)
for (int ii=0; ii<kk2; ii++)
kernel[ii] /= cnt;
// dumpKernel(kernel, k);
}
void gaussian(double* kernel, int k, double ss)
{
int kk = 2*k+1;
int k2 = k*k;
int kk2 = kk*kk;
double s2 = ss*ss;
double tt =0;
for (int yy=-k; yy<=k; yy++)
for (int xx=-k; xx<=k; xx++) {
if ((xx*xx + yy*yy) <= k2) {
double vv = exp(-.5*((xx*xx + yy*yy)/s2));
kernel[(yy+k)*kk+(xx+k)] = vv;
tt += vv;
}
}
// normalize kernel
if (tt)
for (int ii=0; ii<kk2; ii++)
kernel[ii] /= tt;
// dumpKernel(kernel, k);
}
void elliptic(double* kernel, 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 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);
double tt =0;
for (int yy=-k; yy<=k; yy++) {
for (int xx=-k; xx<=k; xx++) {
double dd = xx*cos(aa)+yy*sin(aa);
double ee = xx*sin(aa)-yy*cos(aa);
if ((dd*dd)/(k*k) + (ee*ee)/(rm*rm) <= 1) {
double vv = exp(-(a*xx*xx + 2*b*xx*yy + c*yy*yy));
kernel[(yy+k)*kk+(xx+k)] = vv;
tt += vv;
}
}
}
// normalize kernel
if (tt)
for (int ii=0; ii<kk2; ii++)
kernel[ii] /= tt;
// dumpKernel(kernel, k);
}
void dumpKernel(double* kernel, int k)
{
int kk = 2*k+1;
for (int yy=-k; yy<=k; yy++)
for (int xx=-k; xx<=k; xx++)
cerr << '(' << xx << ',' << yy << ")="
<< kernel[(yy+k)*kk+(xx+k)] << endl;
}
|
