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 <cmath>
34 : #include "heightfield_util.h"
35 : #include "field_util.h"
36 : #include "common/util/file_util.h"
37 : #include "lib_grid/file_io/file_io_asc.h"
38 : using namespace std;
39 :
40 : namespace ug{
41 :
42 : ////////////////////////////////////////////////////////////////////////////////
43 0 : Heightfield::Heightfield() :
44 : m_cellSize(1, 1),
45 : m_offset(0, 0),
46 0 : m_noDataValue(1e30)
47 : {
48 0 : }
49 :
50 0 : number Heightfield::
51 : interpolate(number x, number y, int interpOrder) const
52 : {
53 0 : switch(interpOrder){
54 0 : case 0: {
55 0 : pair<int, int> ind = coordinate_to_index(x, y);
56 0 : if( ind.first >= 0 && ind.first < (int)m_field.width() &&
57 0 : ind.second >= 0 && ind.second < (int)m_field.height())
58 : {
59 0 : return m_field.at(ind.first, ind.second);
60 : }
61 : }break;
62 :
63 0 : case 1:{
64 0 : int cx = int((x - m_offset.x()) / m_cellSize.x());
65 0 : int cy = int((y - m_offset.y()) / m_cellSize.y());
66 :
67 0 : const number vx = (x - ((number)cx * m_cellSize.x() + m_offset.x()))
68 0 : / m_cellSize.x();
69 0 : const number vy = (y - ((number)cy * m_cellSize.y() + m_offset.y()))
70 0 : / m_cellSize.y();
71 0 : const number ux = 1-vx;
72 0 : const number uy = 1-vy;
73 :
74 : if( (cx >= 0)
75 0 : && (cy >= 0)
76 0 : && (cx + 1 < (int)m_field.width())
77 0 : && (cy + 1 < (int)m_field.height()) )
78 : {
79 :
80 0 : return ux*uy*m_field.at(cx, cy) + vx*uy*m_field.at(cx+1,cy)
81 0 : + ux*vy*m_field.at(cx, cy+1) + vx*vy*m_field.at(cx+1, cy+1);
82 : }
83 : else{
84 0 : const int x0 = max(cx, 0);
85 0 : const int x1 = min(cx + 1, (int)m_field.width() - 1);
86 0 : const int y0 = max(cy, 0);
87 0 : const int y1 = min(cy + 1, (int)m_field.height() - 1);
88 :
89 0 : if(x0 >= (int)m_field.width() || x1 < 0 || y0 >= (int)m_field.height() || y1 < 0)
90 0 : return m_noDataValue;
91 :
92 0 : return ux*uy*m_field.at(x0, y0) + vx*uy*m_field.at(x1,y0)
93 0 : + ux*vy*m_field.at(x0, y1) + vx*vy*m_field.at(x1, y1);
94 :
95 : }
96 : }break;
97 :
98 0 : default:
99 0 : UG_THROW("Currently only interpolation orders 0 and 1 supported in Heightfield::interpolate")
100 : }
101 :
102 0 : return m_noDataValue;
103 : }
104 :
105 0 : std::pair<int, int> Heightfield::
106 : coordinate_to_index(number x, number y) const
107 : {
108 : // roundOffset 0: value is constant across each cell
109 : // roundOffset 0.5: value is constant around each node
110 : const number roundOffset = 0.5;
111 :
112 : pair<int, int> c;
113 0 : if(m_cellSize.x() != 0)
114 0 : c.first = (int)(roundOffset + (x - m_offset.x()) / m_cellSize.x());
115 : else
116 : c.first = 0;
117 :
118 0 : if(m_cellSize.y() != 0)
119 0 : c.second = (int)(roundOffset + (y - m_offset.y()) / m_cellSize.y());
120 : else
121 : c.second = 0;
122 0 : return c;
123 : }
124 :
125 0 : vector2 Heightfield::index_to_coordinate(int ix, int iy) const
126 : {
127 0 : return vector2( m_offset.x() + (number)ix * m_cellSize.x(),
128 0 : m_offset.y() + (number)iy * m_cellSize.y());
129 : }
130 :
131 0 : vector2 Heightfield::extent() const
132 : {
133 0 : return vector2(m_cellSize.x() * (number)m_field.width(),
134 0 : m_cellSize.y() * (number)m_field.height());
135 : }
136 :
137 0 : void Heightfield::blur(number alpha, size_t numIterations)
138 : {
139 0 : BlurField(field(), alpha, numIterations, no_data_value());
140 0 : }
141 :
142 :
143 0 : bool Heightfield::eliminate_invalid_cells()
144 : {
145 0 : return EliminateInvalidCells(field(), no_data_value());
146 : }
147 :
148 :
149 :
150 :
151 : ////////////////////////////////////////////////////////////////////////////////
152 0 : void CreateGridFromField(Grid& grid,
153 : const Heightfield& hfield,
154 : Grid::VertexAttachmentAccessor<APosition> aaPos)
155 : {
156 0 : CreateGridFromField(grid, hfield.field(), hfield.cell_size(), hfield.offset(),
157 : hfield.no_data_value(), aaPos);
158 0 : }
159 :
160 : ////////////////////////////////////////////////////////////////////////////////
161 0 : void CreateGridFromFieldBoundary(Grid& grid,
162 : const Heightfield& hfield,
163 : Grid::VertexAttachmentAccessor<APosition> aaPos)
164 : {
165 0 : CreateGridFromFieldBoundary(grid, hfield.field(), hfield.cell_size(), hfield.offset(),
166 : hfield.no_data_value(), aaPos);
167 0 : }
168 :
169 : ////////////////////////////////////////////////////////////////////////////////
170 0 : void LoadHeightfieldFromASC(Heightfield& hfield, const char* filename)
171 : {
172 0 : std::string name = FindFileInStandardPaths(filename);
173 0 : UG_COND_THROW(name.empty(), "Couldn't locate file " << filename);
174 :
175 0 : FileReaderASC reader;
176 0 : reader.set_field(&hfield.field());
177 0 : reader.load_file(name.c_str());
178 : hfield.set_cell_size(vector2(reader.cell_size(), reader.cell_size()));
179 : hfield.set_offset(vector2(reader.lower_left_corner_x(),
180 : reader.lower_left_corner_y()));
181 : hfield.set_no_data_value(reader.no_data_value());
182 0 : }
183 :
184 : }// end of namespace
|
