Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
field_util_impl.h
Go to the documentation of this file.
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
40namespace ug{
41
42template <class T>
43void BlurField(Field<T>& field, number alpha, size_t numIterations, const T& noDataValue)
44{
45 using namespace std;
46 for(size_t mainIter = 0; mainIter < numIterations; ++mainIter){
47 for(int iy = 0; iy < (int)field.height(); ++iy){
48 for(int ix = 0; ix < (int)field.width(); ++ix){
49 if(field.at(ix, iy) != noDataValue){
50 T val = 0;
51 number num = 0;
52 for(int iny = max<int>(iy - 1, 0); iny < min<int>(iy + 2, (int)field.height()); ++iny){
53 for(int inx = max<int>(ix - 1, 0); inx < min<int>(ix + 2, (int)field.width()); ++inx){
54 if(!(inx == 0 && iny == 0) && (field.at(inx, iny) != noDataValue)){
55 val += field.at(inx, iny);
56 ++num;
57 }
58 }
59 }
60
61 if(num > 0){
62 val *= alpha / num;
63 field.at(ix, iy) *= (1.-alpha);
64 field.at(ix, iy) += val;
65 }
66 }
67 }
68 }
69 }
70}
71
72namespace fieldutil{
73template <class T>
74struct Cell{
75 Cell(int _x, int _y, T _val) : x(_x), y(_y), value(_val) {}
76 int x;
77 int y;
79};
80}
81
82template <class T>
83bool EliminateInvalidCells(Field<T>& field, const T& noDataValue)
84{
85 using namespace std;
86 typedef fieldutil::Cell<T> Cell;
87
88 deque<Cell> cells;
89 number inProgressValue = -noDataValue;
90
91 const int numNbrs = 8;
92 const int xadd[numNbrs] = {-1, 0, 1, -1, 1, -1, 0, 1};
93 const int yadd[numNbrs] = {-1, -1, -1, 0, 0, 1, 1, 1};
94
95 const int maxNumSteps = 4;
96 const int minNumValidNbrsInStep[maxNumSteps] = {4, 3, 2, 1};
97
98// initially count the number of invalid cells
99 size_t numInvalidCells = 0;
100 for(int iy = 0; iy < (int)field.height(); ++iy){
101 for(int ix = 0; ix < (int)field.width(); ++ix){
102 if(field.at(ix, iy) == noDataValue)
103 ++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 for(int istep = 0; (istep < maxNumSteps) && (numInvalidCells > 0); ++istep){
112 const int minNumValidNbrs = minNumValidNbrsInStep[istep];
113
114 for(int iy = 0; iy < (int)field.height(); ++iy){
115 for(int ix = 0; ix < (int)field.width(); ++ix){
116 if(field.at(ix, iy) != noDataValue)
117 continue;
118
119 int numValidNbrs = 0;
120 for(int i = 0; i < numNbrs; ++i){
121 const int inx = ix + xadd[i];
122 const int iny = iy + yadd[i];
123
124 if((inx >= 0) && (inx < (int)field.width())
125 && (iny >= 0) && (iny < (int)field.height())
126 && (field.at(inx, iny) != noDataValue)
127 && (field.at(inx, iny) != inProgressValue))
128 {
129 ++numValidNbrs;
130 }
131 }
132
133 if(numValidNbrs >= minNumValidNbrs){
134 cells.push_back(Cell(ix, iy, inProgressValue));
135 field.at(ix, iy) = inProgressValue;
136 break;
137 }
138 }
139 }
140
141 while(!cells.empty()){
142 // iterate over all entries in the queue and calculate their correct values
143 for(typename deque<Cell>::iterator cellIter = cells.begin(); cellIter != cells.end(); ++cellIter)
144 {
145 Cell& cell = *cellIter;
146 const int ix = cell.x;
147 const int iy = cell.y;
148
149 // get average value of valid neighbor cells
150 T avVal = 0;
151 number numValidNbrs = 0;
152 for(int i = 0; i < numNbrs; ++i){
153 const int inx = ix + xadd[i];
154 const int iny = iy + yadd[i];
155
156 if((inx >= 0) && (inx < (int)field.width())
157 && (iny >= 0) && (iny < (int)field.height()))
158 {
159 if((field.at(inx, iny) != noDataValue)
160 && (field.at(inx, iny) != inProgressValue))
161 {
162 avVal += field.at(inx, iny);
163 ++numValidNbrs;
164 }
165 }
166 }
167
168 UG_COND_THROW(numValidNbrs < (number)minNumValidNbrs, "Implementation error!");
169 avVal *= 1. / numValidNbrs;
170 cell.value = avVal;
171 --numInvalidCells;
172 }
173
174 // copy values to the field and collect new candidates
175 while(!(cells.empty() || (cells.front().value == inProgressValue)))
176 {
177 const int ix = cells.front().x;
178 const int iy = cells.front().y;
179 field.at(ix, iy) = cells.front().value;
180 cells.pop_front();
181
182 for(int i = 0; i < numNbrs; ++i){
183 const int inx = ix + xadd[i];
184 const int iny = iy + yadd[i];
185
186 if((inx >= 0) && (inx < (int)field.width())
187 && (iny >= 0) && (iny < (int)field.height()))
188 {
189 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 for(int j = 0; j < numNbrs; ++j){
195 const int inxNbr = inx + xadd[j];
196 const int inyNbr = iny + yadd[j];
197
198 if((inxNbr >= 0) && (inxNbr < (int)field.width())
199 && (inyNbr >= 0) && (inyNbr < (int)field.height())
200 && (field.at(inxNbr, inyNbr) != noDataValue)
201 && (field.at(inxNbr, inyNbr) != inProgressValue))
202 {
203 ++numValidNbrsOfNbr;
204 }
205 }
206 if(numValidNbrsOfNbr >= minNumValidNbrs){
207 cells.push_back(Cell(inx, iny, inProgressValue));
208 field.at(inx, iny) = inProgressValue;
209 }
210 }
211 }
212 }
213 }
214 }
215 }
216
217 return numInvalidCells == 0;
218}
219
220
221template <class T>
222void 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 const int xadd[numNbrs] = {0, -1, 1, 0};
232 const int yadd[numNbrs] = {-1, 0, 0, 1};
233
234// this field stores whether we already visited the given cell
235 Field<bool> visited(field.width(), field.height(), false);
236 vector<pair<int, int> > cells;
237
238 const int fwidth = (int)field.width();
239 const int fheight = (int)field.height();
240
241 for(int outerIy = 0; outerIy < fheight; ++outerIy){
242 for(int outerIx = 0; outerIx < fwidth; ++outerIx){
243 if(visited.at(outerIx, outerIy)
244 || (field.at(outerIx, outerIy) == noDataValue))
245 {
246 continue;
247 }
248
249 cells.clear();
250 cells.push_back(make_pair(outerIx, outerIy));
251 size_t curCell = 0;
252 while(curCell < cells.size()){
253 int ix = cells[curCell].first;
254 int iy = cells[curCell].second;
255
256 for(size_t inbr = 0; inbr < numNbrs; ++inbr){
257 int nx = ix + xadd[inbr];
258 int ny = iy + yadd[inbr];
259 if((nx >= 0 && nx < fwidth && ny >= 0 && ny < fheight)
260 && !visited.at(nx, ny))
261 {
262 visited.at(nx, ny) = true;
263 if(field.at(nx, ny) != noDataValue){
264 cells.push_back(make_pair(nx, ny));
265 }
266 }
267 }
268 ++curCell;
269 }
270
271 if(cells.size() < thresholdCellCount){
272 for(size_t i = 0; i < cells.size(); ++i){
273 int ix = cells[i].first;
274 int iy = cells[i].second;
275 field.at(ix, iy) = noDataValue;
276 }
277 }
278 }
279 }
280}
281
282}// end of namespace
283
284#endif //__H__UG_field_util_impl
Definition field.h:42
size_t width() const
Definition field.h:57
T & at(size_t x, size_t y)
Definition field_impl.hpp:113
size_t height() const
Definition field.h:58
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
double number
Definition types.h:124
Definition smart_pointer.h:814
the ug namespace
void BlurField(Field< T > &field, number alpha, size_t numIterations, const T &noDataValue)
Smoothens the field by adjusting the value of each pixel towards the average of its neighbours.
Definition field_util_impl.h:43
bool EliminateInvalidCells(Field< T > &field, const T &noDataValue)
eliminates invalid cells by repeatedly filling those cells with averages of neighboring cells
Definition field_util_impl.h:83
void InvalidateSmallLenses(Field< T > &field, size_t thresholdCellCount, const T &noDataValue)
invalidates cells that belong to a small lense
Definition field_util_impl.h:222
Definition field_util_impl.h:74
int x
Definition field_util_impl.h:76
int y
Definition field_util_impl.h:77
T value
Definition field_util_impl.h:78
Cell(int _x, int _y, T _val)
Definition field_util_impl.h:75