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 : #include <vector>
34 : #include "raster_layer_util.h"
35 : #include "field_util.h"
36 : #include "lib_grid/file_io/file_io_asc.h"
37 :
38 : using namespace std;
39 :
40 : namespace ug{
41 :
42 : // initialize all values with simple-values and record a list of all invalid ones
43 : struct CellIdx {
44 : CellIdx() {}
45 0 : CellIdx(size_t x, size_t y) : ix(x), iy(y) {}
46 : size_t ix;
47 : size_t iy;
48 : };
49 :
50 0 : void RasterLayers::
51 : load_from_files(const std::vector<LayerDesc>& layerDescs)
52 : {
53 0 : resize(layerDescs.size());
54 0 : for(size_t i = 0; i < layerDescs.size(); ++i){
55 0 : LoadHeightfieldFromASC(heightfield(i), layerDescs[i].filename().c_str());
56 : set_min_height(i, layerDescs[i].min_height());
57 : }
58 0 : }
59 :
60 0 : void RasterLayers::
61 : load_from_files(const std::vector<SPLayerDesc>& layerDescs)
62 : {
63 : vector<LayerDesc> descs;
64 0 : descs.reserve(layerDescs.size());
65 0 : for_each_in_vec(const SPLayerDesc& d, layerDescs){
66 0 : descs.push_back(*d);
67 : }end_for;
68 0 : load_from_files(descs);
69 0 : }
70 :
71 0 : void RasterLayers::
72 : load_from_files(const std::vector<std::string>& filenames, number minHeight)
73 : {
74 : vector<LayerDesc> descs;
75 0 : descs.reserve(filenames.size());
76 0 : for_each_in_vec(const std::string& fname, filenames){
77 0 : descs.push_back(LayerDesc(fname, minHeight));
78 : }end_for;
79 0 : load_from_files(descs);
80 0 : }
81 :
82 0 : void RasterLayers::
83 : resize(size_t newSize)
84 : {
85 : size_t oldSize = size();
86 0 : m_layers.resize(newSize);
87 0 : for(size_t i = oldSize; i < newSize; ++i){
88 0 : m_layers[i] = make_sp(new layer_t);
89 : }
90 0 : }
91 :
92 0 : void RasterLayers::
93 : invalidate_flat_cells()
94 : {
95 : // since the top-layer is considered to be the terrain surface, we'll start
96 : // one layer below the top.
97 0 : if(size() <= 1)
98 : return;
99 0 : for(int lvl = (int)size() - 2; lvl >= 0; --lvl){
100 0 : Heightfield& curHF = heightfield(lvl);
101 : Field<number>& innerField = curHF.field();
102 : number minHeight = min_height(lvl);
103 :
104 : // iterate over the cells of each heightfield and invalidate values
105 0 : for(int iy = 0; iy < (int)innerField.height(); ++iy){
106 0 : for(int ix = 0; ix < (int)innerField.width(); ++ix){
107 0 : number val = innerField.at(ix, iy);
108 0 : if(val != curHF.no_data_value()){
109 : bool gotDataVal = false;
110 0 : vector2 c = curHF.index_to_coordinate(ix, iy);
111 : // compare the value against values from higher layers
112 0 : for(int ulvl = lvl + 1; ulvl < (int)size(); ++ulvl){
113 0 : Heightfield& uHF = heightfield(ulvl);
114 : number uval = uHF.interpolate(c);
115 0 : if(uval != uHF.no_data_value()){
116 : gotDataVal = true;
117 0 : if(val + minHeight > uval){
118 : // the height in the current cell is too small
119 0 : innerField.at(ix, iy) = curHF.no_data_value();
120 : }
121 : break;
122 : }
123 : }
124 :
125 : if(!gotDataVal){
126 : // The cell lies outside of the domain which is specified
127 : // by the upmost field. The cell thus has to be invalidated
128 0 : innerField.at(ix, iy) = curHF.no_data_value();
129 : }
130 : }
131 : }
132 : }
133 : }
134 : }
135 :
136 :
137 0 : void RasterLayers::
138 : invalidate_small_lenses(number minArea)
139 : {
140 0 : for(size_t lvl = 0; lvl < size(); ++lvl){
141 : Heightfield& hf = heightfield(lvl);
142 0 : number cellSize = hf.cell_size().x() * hf.cell_size().y();
143 0 : if(cellSize > 0){
144 0 : InvalidateSmallLenses(hf.field(), minArea / cellSize,
145 0 : hf. no_data_value());
146 : }
147 : }
148 0 : }
149 :
150 0 : void RasterLayers::
151 : remove_small_holes(number maxArea)
152 : {
153 : using namespace std;
154 :
155 : const size_t numNbrs = 4;
156 0 : const int xadd[numNbrs] = {0, -1, 1, 0};
157 0 : const int yadd[numNbrs] = {-1, 0, 0, 1};
158 :
159 : Field<bool> visited;
160 : vector<pair<int, int> > cells;
161 :
162 0 : for(int lvl = (int)size() - 2; lvl >= 0; --lvl){
163 0 : Heightfield& hf = heightfield(lvl);
164 : Field<number>& field = hf.field();
165 : number noDataValue = hf.no_data_value();
166 : number minHeight = min_height(lvl);
167 :
168 0 : number cellSize = hf.cell_size().x() * hf.cell_size().y();
169 0 : if(cellSize <= 0)
170 0 : continue;
171 :
172 0 : size_t thresholdCellCount(maxArea / cellSize);
173 :
174 :
175 0 : visited.resize_no_copy(field.width(), field.height());
176 : visited.fill_all(false);
177 :
178 0 : const int fwidth = (int)field.width();
179 0 : const int fheight = (int)field.height();
180 :
181 0 : for(int outerIy = 0; outerIy < fheight; ++outerIy){
182 0 : for(int outerIx = 0; outerIx < fwidth; ++outerIx){
183 0 : if(visited.at(outerIx, outerIy)
184 0 : || (field.at(outerIx, outerIy) != noDataValue))
185 : {
186 0 : continue;
187 : }
188 :
189 : cells.clear();
190 0 : cells.push_back(make_pair(outerIx, outerIy));
191 : size_t curCell = 0;
192 0 : while(curCell < cells.size()){
193 0 : int ix = cells[curCell].first;
194 0 : int iy = cells[curCell].second;
195 :
196 0 : for(size_t inbr = 0; inbr < numNbrs; ++inbr){
197 0 : int nx = ix + xadd[inbr];
198 0 : int ny = iy + yadd[inbr];
199 0 : if((nx >= 0 && nx < fwidth && ny >= 0 && ny < fheight)
200 0 : && !visited.at(nx, ny))
201 : {
202 0 : visited.at(nx, ny) = true;
203 0 : if(field.at(nx, ny) == noDataValue){
204 0 : cells.push_back(make_pair(nx, ny));
205 : }
206 : }
207 : }
208 0 : ++curCell;
209 : }
210 0 : if(cells.size() < thresholdCellCount){
211 0 : for(size_t i = 0; i < cells.size(); ++i){
212 0 : int ix = cells[i].first;
213 0 : int iy = cells[i].second;
214 0 : vector2 c = hf.index_to_coordinate(ix, iy);
215 0 : pair<int, number> result = trace_line_up(c, (size_t)lvl);
216 0 : if(result.first > lvl){
217 : // we have to adjust not only the value in the current layer
218 : // but also possibly values in lower levels, to avoid
219 : // creating new thin layers
220 0 : number curVal = result.second - minHeight;
221 0 : field.at(ix, iy) = curVal;
222 0 : for(int lowerLvl = lvl - 1; lowerLvl >= 0 ; --lowerLvl){
223 0 : Heightfield& lhf = heightfield(lowerLvl);
224 : Field<number>& lfield = lhf.field();
225 0 : pair<int, int> index = lhf.coordinate_to_index(c.x(), c.y());
226 0 : int lx = index.first;
227 0 : int ly = index.second;
228 0 : if(lx >= 0 && lx < (int)lfield.width()
229 0 : && ly >= 0 && ly < (int)lfield.height())
230 : {
231 0 : number lval = lfield.at(lx, ly);
232 0 : if(lval != lhf.no_data_value()){
233 0 : if(curVal - lval < minHeight){
234 0 : curVal -= minHeight;
235 0 : lfield.at(lx, ly) = curVal;
236 : }
237 : else
238 : break;
239 : }
240 : }
241 : }
242 : }
243 : }
244 : }
245 : }
246 : }
247 : }
248 0 : }
249 :
250 :
251 0 : void RasterLayers::
252 : snap_cells_to_higher_layers()
253 : {
254 0 : if(size() <= 1)
255 : return;
256 :
257 0 : for(int lvl = (int)size() - 2; lvl >= 0; --lvl){
258 0 : Heightfield& curHF = heightfield(lvl);
259 : Field<number>& curField = curHF.field();
260 0 : Heightfield& upperHF = heightfield(lvl + 1);
261 : number minHeight = min_height(lvl);
262 :
263 0 : for(size_t iy = 0; iy < curField.height(); ++iy){
264 0 : for(size_t ix = 0; ix < curField.width(); ++ix){
265 0 : number curVal = curField.at(ix, iy);
266 0 : vector2 c = curHF.index_to_coordinate(ix, iy);
267 : number upperVal = upperHF.interpolate(c);
268 0 : if(upperVal == upperHF.no_data_value())
269 0 : continue;
270 :
271 0 : if((curVal == curHF.no_data_value()) || (curVal > upperVal)
272 0 : || (upperVal - curVal < minHeight))
273 : {
274 0 : curField.at(ix, iy) = upperVal;
275 : }
276 : }
277 : }
278 : }
279 : }
280 :
281 :
282 0 : void RasterLayers::
283 : eliminate_invalid_cells()
284 : {
285 0 : for(size_t lvl = 0; lvl < size(); ++lvl){
286 : Heightfield& hf = heightfield(lvl);
287 0 : EliminateInvalidCells(hf.field(), hf.no_data_value());
288 : }
289 0 : }
290 :
291 :
292 0 : void RasterLayers::
293 : blur_layers(number alpha, size_t numIterations)
294 : {
295 0 : for(size_t i = 0; i < size(); ++i){
296 0 : BlurField(heightfield(i).field(), alpha, numIterations, heightfield(i).no_data_value());
297 : }
298 0 : }
299 :
300 0 : std::pair<int, number> RasterLayers::
301 : trace_line_down(const vector2& c, size_t firstLayer) const
302 : {
303 0 : for(int i = (int)firstLayer; i >= 0; --i){
304 0 : number val = heightfield(i).interpolate(c);
305 0 : if(val != heightfield(i).no_data_value()){
306 : return make_pair(i, val);
307 : }
308 : }
309 :
310 0 : return make_pair<int, number>(-1, 0);
311 : }
312 :
313 0 : std::pair<int, number> RasterLayers::
314 : trace_line_up(const vector2& c, size_t firstLayer) const
315 : {
316 0 : for(size_t i = firstLayer; i < size(); ++i){
317 : number val = heightfield(i).interpolate(c);
318 0 : if(val != heightfield(i).no_data_value()){
319 0 : return make_pair((int)i, val);
320 : }
321 : }
322 :
323 0 : return make_pair<int, number>(-1, 0);
324 : }
325 :
326 0 : std::pair<int, int> RasterLayers::
327 : get_layer_indices(const vector3& c) const
328 : {
329 0 : vector2 xy(c.x(), c.y());
330 :
331 : int upperLayer = -1;
332 : int lowerLayer = -1;
333 : while(1){
334 0 : std::pair<int, number> cut = trace_line_down(xy, upperLayer + 1);
335 0 : if(cut.first == -1)
336 : break;
337 0 : else if(cut.second > c.z()){
338 : lowerLayer = cut.first;
339 : break;
340 : }
341 : else
342 : upperLayer = cut.first;
343 0 : }
344 :
345 0 : return make_pair(upperLayer, lowerLayer);
346 : }
347 :
348 0 : number RasterLayers::
349 : relative_to_global_height(const vector2& c, number relHeight) const
350 : {
351 : static const int order = 1;
352 :
353 0 : if(m_relativeToGlobalHeights.empty())
354 0 : return relative_to_global_height_simple(c, relHeight);
355 : else{
356 0 : const int relHeightLow = (int)(relHeight + SMALL);
357 0 : const int topLvl = (int)m_relativeToGlobalHeights.size() - 1;
358 0 : if((relHeightLow >= 0)
359 0 : && (relHeightLow + 1 <= topLvl))
360 : {
361 0 : number rel = relHeight - (number)relHeightLow;
362 0 : return (1. - rel) * m_relativeToGlobalHeights[relHeightLow]
363 : ->interpolate(c, order)
364 0 : + rel * m_relativeToGlobalHeights[relHeightLow + 1]
365 0 : ->interpolate(c, order);
366 : }
367 0 : else if(relHeightLow >= topLvl)
368 0 : return m_relativeToGlobalHeights[topLvl]->interpolate(c, order);
369 : else
370 0 : return m_relativeToGlobalHeights[0]->interpolate(c, order);
371 : }
372 : }
373 :
374 0 : number RasterLayers::
375 : relative_to_global_height_simple(const vector2& c, number relHeight) const
376 : {
377 0 : const int relHeightLower = max((int)(relHeight + SMALL), 0);
378 0 : const int relHeightUpper = min(relHeightLower + 1, (int)num_layers() - 1);
379 0 : const std::pair<int, number> lower = trace_line_down(c, relHeightLower);
380 0 : const std::pair<int, number> upper = trace_line_up(c, relHeightUpper);
381 :
382 0 : UG_COND_THROW(lower.first < 0, "Invalid lower layer for coordinate " << c
383 : << " and relative height " << relHeight);
384 0 : UG_COND_THROW(upper.first < 0, "Invalid upper layer for coordinate " << c
385 : << " and relative height " << relHeight);
386 :
387 0 : if(upper.second < lower.second)
388 : return upper.second;
389 :
390 0 : number layerDiff = max(1, upper.first - lower.first);
391 0 : number rel = (clip<number>(relHeight - (number)relHeightLower, 0, 1)
392 0 : + relHeightLower - lower.first)
393 0 : / layerDiff;
394 0 : return (1. - rel) * lower.second + rel * upper.second;
395 : }
396 :
397 0 : void RasterLayers::
398 : construct_relative_to_global_height_table(size_t iterations, number alpha)
399 : {
400 : m_relativeToGlobalHeights.clear();
401 0 : for(size_t i = 0; i < m_layers.size(); ++i){
402 0 : m_relativeToGlobalHeights.push_back(
403 0 : make_sp(new Heightfield(m_layers[i]->heightfield)));
404 : }
405 :
406 :
407 : // make sure that the bottom-layer has no holes (compared to the top layer).
408 : // we'll set it to the same height as its local upper layer.
409 : {
410 : Heightfield& baseHF = *m_relativeToGlobalHeights[0];
411 : Field<number>& baseField = baseHF.field();
412 :
413 0 : for(size_t iy = 0; iy < baseField.height(); ++iy){
414 0 : for(size_t ix = 0; ix < baseField.width(); ++ix){
415 0 : if(baseField.at(ix, iy) == baseHF.no_data_value()){
416 0 : vector2 c = baseHF.index_to_coordinate(ix, iy);
417 0 : std::pair<int, number> upper = trace_line_up(c, 0);
418 0 : if(upper.first != -1)
419 0 : baseField.at(ix, iy) = upper.second;
420 : }
421 : }
422 : }
423 : }
424 :
425 :
426 : std::vector<std::vector<CellIdx> > allCells;
427 0 : allCells.resize(m_relativeToGlobalHeights.size());
428 :
429 0 : for(int lvl = 1; lvl + 1 < (int)m_relativeToGlobalHeights.size(); ++lvl){
430 0 : Heightfield& curHF = *m_relativeToGlobalHeights[lvl];
431 : Field<number>& curField = curHF.field();
432 :
433 : std::vector<CellIdx>& cells = allCells[lvl];
434 :
435 0 : for(size_t iy = 0; iy < curField.height(); ++iy){
436 0 : for(size_t ix = 0; ix < curField.width(); ++ix){
437 0 : if(curField.at(ix, iy) == curHF.no_data_value()){
438 0 : vector2 c = curHF.index_to_coordinate(ix, iy);
439 0 : if(trace_line_up(c, lvl).first != -1){
440 0 : curField.at(ix, iy) = relative_to_global_height_simple(c, lvl);
441 0 : cells.push_back(CellIdx(ix, iy));
442 : }
443 : }
444 : }
445 : }
446 : }
447 :
448 : // now iterate over all recorded cells and perform relaxation
449 0 : for(size_t iteration = 0; iteration < iterations; ++iteration){
450 0 : for(int lvl = 1; lvl + 1 < (int)m_relativeToGlobalHeights.size(); ++lvl){
451 0 : Field<number>& cur = m_relativeToGlobalHeights[lvl]->field();
452 0 : Field<number>& lower = m_relativeToGlobalHeights[lvl-1]->field();
453 0 : Field<number>& upper = m_relativeToGlobalHeights[lvl+1]->field();
454 :
455 : std::vector<CellIdx>& cells = allCells[lvl];
456 0 : for(size_t icell = 0; icell < cells.size(); ++icell){
457 0 : const CellIdx ci = cells[icell];
458 : // average dist-relations of neighbors
459 : number nbrVal = 0;
460 : number numNbrs = 0;
461 :
462 0 : if(ci.ix > 0){
463 0 : nbrVal += upper_lower_dist_relation(
464 : lower, cur,
465 0 : upper, ci.ix - 1, ci.iy);
466 : ++numNbrs;
467 : }
468 0 : if(ci.ix + 1 < cur.width()){
469 :
470 0 : nbrVal += upper_lower_dist_relation(
471 : lower, cur,
472 0 : upper, ci.ix + 1, ci.iy);
473 0 : ++numNbrs;
474 : }
475 0 : if(ci.iy > 0){
476 0 : nbrVal += upper_lower_dist_relation(
477 : lower, cur,
478 : upper, ci.ix, ci.iy - 1);
479 0 : ++numNbrs;
480 : }
481 0 : if(ci.iy + 1 < cur.height()){
482 0 : nbrVal += upper_lower_dist_relation(
483 : lower, cur,
484 : upper, ci.ix, ci.iy + 1);
485 0 : ++numNbrs;
486 : }
487 :
488 0 : if(numNbrs > 0){
489 0 : number locVal = upper_lower_dist_relation(
490 : lower, cur,
491 : upper, ci.ix, ci.iy);
492 :
493 0 : locVal = (1. - alpha) * locVal + alpha * nbrVal / numNbrs;
494 :
495 : // reconstruct height value from dist-relation
496 0 : const number upperVal = upper.at(ci.ix, ci.iy);
497 0 : const number lowerVal = lower.at(ci.ix, ci.iy);
498 :
499 0 : number dirTotal = upperVal - lowerVal;
500 0 : cur.at(ci.ix, ci.iy) = upperVal - locVal * dirTotal;
501 :
502 : // enforce minWidth constraints
503 0 : const number upperDist = upperVal - cur.at(ci.ix, ci.iy);
504 0 : const number lowerDist = cur.at(ci.ix, ci.iy) - lowerVal;
505 0 : if(upperDist < layer(lvl).minHeight)
506 0 : if(lowerDist < layer(lvl-1).minHeight)
507 : cur.at(ci.ix, ci.iy)
508 0 : = 0.5 * (upperVal - upperDist + lowerVal + lowerDist);
509 : else
510 0 : cur.at(ci.ix, ci.iy) = upperVal - upperDist;
511 0 : else if(lowerDist < layer(lvl-1).minHeight)
512 0 : cur.at(ci.ix, ci.iy) = lowerVal + lowerDist;
513 : }
514 : }
515 : }
516 : }
517 0 : }
518 :
519 :
520 0 : void RasterLayers::
521 : invalidate_relative_to_global_height_table ()
522 : {
523 : m_relativeToGlobalHeights.clear();
524 0 : }
525 :
526 0 : number RasterLayers::
527 : upper_lower_dist_relation ( Field<number>&lower,
528 : Field<number>& middle,
529 : Field<number>& upper,
530 : size_t ix,
531 : size_t iy)
532 : {
533 0 : number totalDist = fabs(upper.at(ix, iy) - lower.at(ix, iy));
534 0 : if(totalDist > 0)
535 0 : return (upper.at(ix, iy) - middle.at(ix, iy)) / totalDist;
536 : return 0;
537 : }
538 :
539 : // void RasterLayers::
540 : // save_to_files(const char* filenamePrefix_cstr)
541 : // {
542 : // std::string filenamePrefix = filenamePrefix_cstr;
543 : // for(size_t i = 0; i < m_layers.size(); ++i){
544 :
545 : // }
546 : // }
547 :
548 : // void RasterLayers::
549 : // save_rel_to_glob_table_to_files(const char* filenamePrefix_cstr)
550 : // {
551 :
552 : // }
553 :
554 :
555 :
556 : }// end of namespace
|
