85 {
86
89 double data = 0;
90
91 vector<int> vec_ix;
92 vector<int> vec_iy;
93 vector<double> vec_data;
94 vector<double> vec_dist;
95
98
99 int ix = ceil(x /
sCale);
100 int iy = ceil(y /
sCale);
101
102 vec_ix.resize(2 * nx);
103 vec_iy.resize(2 * ny);
104
105 int ii = 0;
106 for (
int i = 0;
i < 2 * nx;
i++) {
107 int idx = ix - nx +
i;
108 if (idx >= size_x / 2 || idx < -size_x / 2)
109 continue;
110 vec_ix[ii++] = idx;
111 }
112 vec_ix.resize(ii);
113 ii = 0;
114 for (
int i = 0;
i < 2 * ny;
i++) {
115 int idx = iy - ny +
i;
116 if (idx >= size_y / 2 || idx < -size_y / 2)
117 continue;
118 vec_iy[ii++] = idx;
119 }
120 vec_iy.resize(ii);
121 ii = 0;
122
123 vec_data.resize(vec_iy.size() * vec_ix.size());
124 vec_dist.resize(vec_iy.size() * vec_ix.size());
125
126 for (vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();
127 it_iy++) {
128 for (vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();
129 it_ix++) {
130 vec_data[ii] =
131 main_vector[*it_iy + size_y / 2][*it_ix + size_x / 2];
132 vec_dist[ii] = ((*it_iy *
sCale - y) * (*it_iy *
sCale - y) +
134
135 ii++;
136 }
137 }
138
139
140 vector<double> kernel;
141 kernel.resize(vec_data.size());
143 double sigma = 10;
144 double sum = 0;
145 const double m = (1 / sqrt(M_PI * 2 * sigma * sigma));
146 for (int ii = 0; ii < vec_dist.size(); ii++) {
147 kernel[
i] =
m * exp(-(vec_dist[ii]) /
148 (2 * sigma * sigma));
150 }
151 ii = 0;
152 for (vector<double>::iterator vec_itr = vec_data.begin();
153 vec_itr != vec_data.end(); vec_itr++) {
154 kernel[ii] /= sum;
155 data += (*vec_itr) * kernel[ii++];
156 }
157 return data;
158 }
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'm', 3 > m
vector< vector< double > > main_vector
