Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG_field_util_impl
34 : #define __H__UG_field_util_impl
35 :
36 : #include <algorithm>
37 : #include <deque>
38 : #include <queue>
39 :
40 : namespace ug{
41 :
42 : template <class T>
43 0 : void BlurField(Field<T>& field, number alpha, size_t numIterations, const T& noDataValue)
44 : {
45 : using namespace std;
46 0 : for(size_t mainIter = 0; mainIter < numIterations; ++mainIter){
47 0 : for(int iy = 0; iy < (int)field.height(); ++iy){
48 0 : for(int ix = 0; ix < (int)field.width(); ++ix){
49 0 : if(field.at(ix, iy) != noDataValue){
50 : T val = 0;
51 : number num = 0;
52 0 : for(int iny = max<int>(iy - 1, 0); iny < min<int>(iy + 2, (int)field.height()); ++iny){
53 0 : for(int inx = max<int>(ix - 1, 0); inx < min<int>(ix + 2, (int)field.width()); ++inx){
54 0 : if(!(inx == 0 && iny == 0) && (field.at(inx, iny) != noDataValue)){
55 0 : val += field.at(inx, iny);
56 0 : ++num;
57 : }
58 : }
59 : }
60 :
61 0 : if(num > 0){
62 0 : val *= alpha / num;
63 0 : field.at(ix, iy) *= (1.-alpha);
64 0 : field.at(ix, iy) += val;
65 : }
66 : }
67 : }
68 : }
69 : }
70 0 : }
71 :
72 : namespace fieldutil{
73 : template <class T>
74 : struct Cell{
75 0 : Cell(int _x, int _y, T _val) : x(_x), y(_y), value(_val) {}
76 : int x;
77 : int y;
78 : T value;
79 : };
80 : }
81 :
82 : template <class T>
83 0 : bool EliminateInvalidCells(Field<T>& field, const T& noDataValue)
84 : {
85 : using namespace std;
86 : typedef fieldutil::Cell<T> Cell;
87 :
88 : deque<Cell> cells;
89 0 : number inProgressValue = -noDataValue;
90 :
91 : const int numNbrs = 8;
92 0 : const int xadd[numNbrs] = {-1, 0, 1, -1, 1, -1, 0, 1};
93 0 : const int yadd[numNbrs] = {-1, -1, -1, 0, 0, 1, 1, 1};
94 :
95 : const int maxNumSteps = 4;
96 0 : const int minNumValidNbrsInStep[maxNumSteps] = {4, 3, 2, 1};
97 :
98 : // initially count the number of invalid cells
99 : size_t numInvalidCells = 0;
100 0 : for(int iy = 0; iy < (int)field.height(); ++iy){
101 0 : for(int ix = 0; ix < (int)field.width(); ++ix){
102 0 : if(field.at(ix, iy) == noDataValue)
103 0 : ++numInvalidCells;
104 : }
105 : }
106 :
107 : // find the initial cells which contain no-data-values and which are neighbors
108 : // of valid cells
109 : // We do this in several steps to better smear out the values and to avoid sharp features
110 : // in the smeared out regions
111 0 : for(int istep = 0; (istep < maxNumSteps) && (numInvalidCells > 0); ++istep){
112 0 : const int minNumValidNbrs = minNumValidNbrsInStep[istep];
113 :
114 0 : for(int iy = 0; iy < (int)field.height(); ++iy){
115 0 : for(int ix = 0; ix < (int)field.width(); ++ix){
116 0 : if(field.at(ix, iy) != noDataValue)
117 0 : continue;
118 :
119 : int numValidNbrs = 0;
120 0 : for(int i = 0; i < numNbrs; ++i){
121 0 : const int inx = ix + xadd[i];
122 0 : const int iny = iy + yadd[i];
123 :
124 0 : if((inx >= 0) && (inx < (int)field.width())
125 0 : && (iny >= 0) && (iny < (int)field.height())
126 0 : && (field.at(inx, iny) != noDataValue)
127 0 : && (field.at(inx, iny) != inProgressValue))
128 : {
129 0 : ++numValidNbrs;
130 : }
131 : }
132 :
133 0 : if(numValidNbrs >= minNumValidNbrs){
134 0 : cells.push_back(Cell(ix, iy, inProgressValue));
135 0 : field.at(ix, iy) = inProgressValue;
136 0 : break;
137 : }
138 : }
139 : }
140 :
141 0 : while(!cells.empty()){
142 : // iterate over all entries in the queue and calculate their correct values
143 0 : for(typename deque<Cell>::iterator cellIter = cells.begin(); cellIter != cells.end(); ++cellIter)
144 : {
145 : Cell& cell = *cellIter;
146 0 : const int ix = cell.x;
147 0 : const int iy = cell.y;
148 :
149 : // get average value of valid neighbor cells
150 : T avVal = 0;
151 : number numValidNbrs = 0;
152 0 : for(int i = 0; i < numNbrs; ++i){
153 0 : const int inx = ix + xadd[i];
154 0 : const int iny = iy + yadd[i];
155 :
156 0 : if((inx >= 0) && (inx < (int)field.width())
157 0 : && (iny >= 0) && (iny < (int)field.height()))
158 : {
159 0 : if((field.at(inx, iny) != noDataValue)
160 0 : && (field.at(inx, iny) != inProgressValue))
161 : {
162 0 : avVal += field.at(inx, iny);
163 0 : ++numValidNbrs;
164 : }
165 : }
166 : }
167 :
168 0 : UG_COND_THROW(numValidNbrs < (number)minNumValidNbrs, "Implementation error!");
169 0 : avVal *= 1. / numValidNbrs;
170 0 : cell.value = avVal;
171 0 : --numInvalidCells;
172 : }
173 :
174 : // copy values to the field and collect new candidates
175 0 : while(!(cells.empty() || (cells.front().value == inProgressValue)))
176 : {
177 0 : const int ix = cells.front().x;
178 0 : const int iy = cells.front().y;
179 0 : field.at(ix, iy) = cells.front().value;
180 : cells.pop_front();
181 :
182 0 : for(int i = 0; i < numNbrs; ++i){
183 0 : const int inx = ix + xadd[i];
184 0 : const int iny = iy + yadd[i];
185 :
186 0 : if((inx >= 0) && (inx < (int)field.width())
187 0 : && (iny >= 0) && (iny < (int)field.height()))
188 : {
189 0 : if(field.at(inx, iny) == noDataValue){
190 : // the nbr-cell is a possible new candidate.
191 : // count the number of valid nbrs-of that cell and
192 : // add it to the queue if there are enough.
193 : int numValidNbrsOfNbr = 0;
194 0 : for(int j = 0; j < numNbrs; ++j){
195 0 : const int inxNbr = inx + xadd[j];
196 0 : const int inyNbr = iny + yadd[j];
197 :
198 0 : if((inxNbr >= 0) && (inxNbr < (int)field.width())
199 0 : && (inyNbr >= 0) && (inyNbr < (int)field.height())
200 0 : && (field.at(inxNbr, inyNbr) != noDataValue)
201 0 : && (field.at(inxNbr, inyNbr) != inProgressValue))
202 : {
203 0 : ++numValidNbrsOfNbr;
204 : }
205 : }
206 0 : if(numValidNbrsOfNbr >= minNumValidNbrs){
207 0 : cells.push_back(Cell(inx, iny, inProgressValue));
208 0 : field.at(inx, iny) = inProgressValue;
209 : }
210 : }
211 : }
212 : }
213 : }
214 : }
215 : }
216 :
217 0 : return numInvalidCells == 0;
218 : }
219 :
220 :
221 : template <class T>
222 0 : void InvalidateSmallLenses(Field<T>& field, size_t thresholdCellCount,
223 : const T& noDataValue)
224 : {
225 : using namespace std;
226 :
227 : // const int numNbrs = 8;
228 : // const int xadd[numNbrs] = {-1, 0, 1, -1, 1, -1, 0, 1};
229 : // const int yadd[numNbrs] = {-1, -1, -1, 0, 0, 1, 1, 1};
230 : const size_t numNbrs = 4;
231 0 : const int xadd[numNbrs] = {0, -1, 1, 0};
232 0 : const int yadd[numNbrs] = {-1, 0, 0, 1};
233 :
234 : // this field stores whether we already visited the given cell
235 0 : Field<bool> visited(field.width(), field.height(), false);
236 : vector<pair<int, int> > cells;
237 :
238 0 : const int fwidth = (int)field.width();
239 0 : const int fheight = (int)field.height();
240 :
241 0 : for(int outerIy = 0; outerIy < fheight; ++outerIy){
242 0 : for(int outerIx = 0; outerIx < fwidth; ++outerIx){
243 0 : if(visited.at(outerIx, outerIy)
244 0 : || (field.at(outerIx, outerIy) == noDataValue))
245 : {
246 0 : continue;
247 : }
248 :
249 : cells.clear();
250 0 : cells.push_back(make_pair(outerIx, outerIy));
251 : size_t curCell = 0;
252 0 : while(curCell < cells.size()){
253 0 : int ix = cells[curCell].first;
254 0 : int iy = cells[curCell].second;
255 :
256 0 : for(size_t inbr = 0; inbr < numNbrs; ++inbr){
257 0 : int nx = ix + xadd[inbr];
258 0 : int ny = iy + yadd[inbr];
259 0 : if((nx >= 0 && nx < fwidth && ny >= 0 && ny < fheight)
260 0 : && !visited.at(nx, ny))
261 : {
262 0 : visited.at(nx, ny) = true;
263 0 : if(field.at(nx, ny) != noDataValue){
264 0 : cells.push_back(make_pair(nx, ny));
265 : }
266 : }
267 : }
268 0 : ++curCell;
269 : }
270 :
271 0 : if(cells.size() < thresholdCellCount){
272 0 : for(size_t i = 0; i < cells.size(); ++i){
273 0 : int ix = cells[i].first;
274 0 : int iy = cells[i].second;
275 0 : field.at(ix, iy) = noDataValue;
276 : }
277 : }
278 : }
279 : }
280 0 : }
281 :
282 : }// end of namespace
283 :
284 : #endif //__H__UG_field_util_impl
|
