Line data Source code
1 : // This file is part of the Mamico project. For conditions of distribution
2 : // and use, please see the copyright notice in Mamico's main folder, or at
3 : // www5.in.tum.de/mamico
4 :
5 : /** Implementation of Gauss.h
6 : * @author Felix Maurer
7 : */
8 :
9 : // Member functions of Gauss.h:
10 : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
11 : // TODO: This is very much hardcoded for kernel radius = 1
12 0 : void coupling::filtering::Gauss<dim, scope...>::operator()() {
13 : #ifdef DEBUG_GAUSS
14 : std::cout << " GAUSS(Dim: " << _dim << "): Applying filter..." << std::endl;
15 : #endif
16 :
17 : using namespace coupling::indexing;
18 :
19 : ScalarIndex index;
20 : ScalarIndex indexAbove;
21 : ScalarIndex indexBelow;
22 :
23 0 : for (unsigned int raw_index = 0; raw_index < _inputCells.size(); raw_index++) {
24 :
25 : // construct CellIndex using raw index
26 0 : index = ScalarIndex{raw_index};
27 :
28 : // make use of implicit conversion between linear <-> vector indices
29 0 : indexAbove = getIndexAbove(index, _dim);
30 0 : indexBelow = getIndexBelow(index, _dim);
31 :
32 : /*
33 : std::cout << "INDEX: " << index << std::endl;
34 : std::cout << "ABOVE: " << indexAbove << std::endl;
35 : std::cout << "BELOW: " << indexBelow << std::endl;
36 : */
37 :
38 : //[0] = below, [1] = at index, [2] = above
39 0 : double weights[3] = {_kernel[0], _kernel[1], _kernel[2]};
40 :
41 : // only one of these two cases can occur at once
42 0 : if (indexBelow == index) {
43 : // set weight of out of bounds index to zero
44 0 : weights[0] = 0;
45 :
46 0 : switch (_extrapolationStrategy) {
47 : case NONE:
48 : break;
49 0 : case MIRROR:
50 0 : weights[2] *= 2;
51 0 : break;
52 0 : case REFLECT:
53 0 : weights[1] *= 2;
54 0 : break;
55 : }
56 :
57 : // normalize the other two weights
58 0 : double weightsSum = weights[1] + weights[2];
59 0 : weights[1] /= weightsSum;
60 0 : weights[2] /= weightsSum;
61 : }
62 :
63 : // Since the domain has to be larger than one unit in all directions, not
64 : // both the index above and below can be illegal at the same time.
65 0 : else if (indexAbove == index) {
66 : // set weight of out of bounds index to zero
67 0 : weights[2] = 0;
68 :
69 0 : switch (_extrapolationStrategy) {
70 : case NONE:
71 : break;
72 0 : case MIRROR:
73 0 : weights[0] *= 2;
74 0 : break;
75 0 : case REFLECT:
76 0 : weights[1] *= 2;
77 0 : break;
78 : }
79 :
80 : // normalize the other two weights
81 0 : double weightsSum = weights[0] + weights[1];
82 0 : weights[0] /= weightsSum;
83 0 : weights[1] /= weightsSum;
84 : }
85 :
86 : // apply to scalars
87 0 : for (const auto scalarProperty : _scalarAccessFunctionPairs) {
88 0 : (_outputCells[raw_index]->*scalarProperty.set)((_inputCells[indexBelow.get()]->*scalarProperty.get)() * weights[0] +
89 0 : (_inputCells[raw_index]->*scalarProperty.get)() * weights[1] +
90 0 : (_inputCells[indexAbove.get()]->*scalarProperty.get)() * weights[2]);
91 : }
92 : // apply to vectors
93 0 : for (const auto vectorProperty : _vectorAccessFunctionPairs) {
94 0 : (_outputCells[raw_index]->*vectorProperty.set)((_inputCells[indexBelow.get()]->*vectorProperty.get)() * weights[0] +
95 0 : (_inputCells[raw_index]->*vectorProperty.get)() * weights[1] +
96 0 : (_inputCells[indexAbove.get()]->*vectorProperty.get)() * weights[2]);
97 : }
98 : } // index
99 :
100 : #ifdef DEBUG_GAUSS
101 : std::cout << "done." << std::endl;
102 : #endif
103 0 : }
104 :
105 : // Private functions of Gauss.h:
106 :
107 : /**
108 : * @brief Computes the one-dimensional Gaussian kernel.
109 : * Kernel size is 1 + 2 * GAUSS_KERNEL_RADIUS, because there is one central cell a
110 : * and GAUSS_KERNEL_RADIUS neighboring cells before and after that.
111 : *
112 : * @returns an array with kernel weights, normalized to sum(weights)=1
113 : */
114 : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
115 0 : std::array<double, 1 + 2 * GAUSS_KERNEL_RADIUS> coupling::filtering::Gauss<dim, scope...>::generateKernel() {
116 : // std::cout << "gauss: kernel radius: " << GAUSS_KERNEL_RADIUS << std::endl;
117 :
118 : std::array<double, 1 + 2 * GAUSS_KERNEL_RADIUS> kernel;
119 0 : unsigned int i = 0; // index in kernel
120 0 : double sum = 0; // used for normalization
121 :
122 : // fill kernel using Gaussian error function
123 0 : for (int x = -GAUSS_KERNEL_RADIUS; x <= GAUSS_KERNEL_RADIUS; x++) {
124 0 : kernel[i] = gaussianDensityFunction(x);
125 :
126 0 : sum += kernel[i];
127 0 : i++;
128 : }
129 :
130 : // normalize kernel. by definition, sum cannot be 0.
131 0 : for (auto& k : kernel) {
132 0 : k *= 1 / sum;
133 :
134 : #ifdef DEBUG_GAUSS
135 : std::cout << " GAUSS(Dim: " << _dim << "): Add kernel entry: " << k << std::endl;
136 : #endif
137 : }
138 :
139 0 : return kernel;
140 : }
141 :
142 : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
143 0 : constexpr double coupling::filtering::Gauss<dim, scope...>::gaussianDensityFunction(int x) {
144 0 : double ans = (1 / (_sigma * sqrt(2 * M_PI))) * exp(-0.5 * (x / _sigma) * (x / _sigma));
145 0 : return ans;
146 : }
147 :
148 : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
149 : typename coupling::filtering::Gauss<dim, scope...>::VectorIndex
150 0 : coupling::filtering::Gauss<dim, scope...>::getIndexAbove(const coupling::filtering::Gauss<dim, scope...>::VectorIndex index, unsigned int d) {
151 : // check for border indices
152 0 : const coupling::indexing::BaseIndex<dim> boundary = coupling::filtering::Gauss<dim, scope...>::VectorIndex::upperBoundary;
153 0 : for (unsigned int d = 0; d < dim; d++) {
154 : // convert index to BaseIndex and compare with boundary.
155 0 : if ((coupling::indexing::BaseIndex<dim>{index}).get()[d] >= boundary.get()[d])
156 : // if index is out of bounds, return it unmodified. this gets detected in
157 : // operator()'s boundary handling.
158 0 : return index;
159 : }
160 :
161 : // get raw vector index
162 0 : tarch::la::Vector<dim, int> ans = index.get();
163 : // get upper index in dimension d
164 0 : ans[d] += 1;
165 :
166 0 : return VectorIndex{ans};
167 : }
168 :
169 : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
170 : typename coupling::filtering::Gauss<dim, scope...>::VectorIndex
171 0 : coupling::filtering::Gauss<dim, scope...>::getIndexBelow(const coupling::filtering::Gauss<dim, scope...>::VectorIndex index, unsigned int d) {
172 : // check for border indices
173 0 : const coupling::indexing::BaseIndex<dim> boundary = coupling::filtering::Gauss<dim, scope...>::VectorIndex::lowerBoundary;
174 0 : for (unsigned int d = 0; d < dim; d++) {
175 : // convert index to BaseIndex and compare with boundary.
176 0 : if ((coupling::indexing::BaseIndex<dim>{index}).get()[d] <= boundary.get()[d])
177 : // if index is out of bounds, return it unmodified. this gets detected in
178 : // operator()'s boundary handling.
179 0 : return index;
180 : }
181 :
182 : // copy input index
183 0 : tarch::la::Vector<dim, int> ans = index.get();
184 : // get lower index in dimension d
185 0 : ans[d] -= 1;
186 :
187 0 : return VectorIndex{ans};
188 : }
|
