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_raster_layer_util
34 : #define __H__UG_raster_layer_util
35 :
36 : #include <utility>
37 : #include <string>
38 : #include <vector>
39 : #include "lib_grid/algorithms/heightfield_util.h"
40 : #include "common/boost_serialization.h"
41 :
42 : namespace ug{
43 :
44 0 : struct RasterLayerDesc{
45 0 : RasterLayerDesc(const std::string& filename, number minHeight) :
46 0 : m_filename(filename), m_minHeight(minHeight) {}
47 :
48 0 : const std::string& filename() const {return m_filename;}
49 0 : number min_height() const {return m_minHeight;}
50 :
51 : private:
52 : std::string m_filename;
53 : number m_minHeight;
54 : };
55 :
56 : typedef SmartPtr<RasterLayerDesc> SPRasterLayerDesc;
57 :
58 :
59 0 : class RasterLayers{
60 : public:
61 0 : struct layer_t{
62 0 : layer_t() : heightfield(), minHeight(SMALL) {}
63 :
64 : Heightfield heightfield;
65 : number minHeight;
66 :
67 : private:
68 : friend class boost::serialization::access;
69 : template <class Archive>
70 0 : void serialize( Archive& ar, const unsigned int version)
71 : {
72 0 : ar & minHeight;
73 0 : ar & heightfield;
74 0 : }
75 : };
76 :
77 :
78 : typedef RasterLayerDesc LayerDesc;
79 : typedef SPRasterLayerDesc SPLayerDesc;
80 :
81 :
82 : /// loads raster data from a list of .asc files.
83 : /** layerDescs.front() represents the bottom of the lowest layer.
84 : * layerDescs.top() represents the terrain surface. All data inbetween
85 : * is interpreted as sorted layer-bottoms from bottom to top.
86 : * \{ */
87 : void load_from_files(const std::vector<LayerDesc>& layerDescs);
88 :
89 : void load_from_files(const std::vector<SPLayerDesc>& layerDescs);
90 : /** \} */
91 :
92 : /// loads raster data from a list of .asc files.
93 : /** filenames.front() represents the bottom of the lowest layer.
94 : * filenames.top() represents the terrain surface. All data inbetween
95 : * is interpreted as sorted layer-bottoms from bottom to top.
96 : *
97 : * \param minLayerHeight If the height of the layer at a given point is
98 : * smaller than this value, then the layer is considered
99 : * to be non-existant at this point (i.e. has a hole)*/
100 : void load_from_files(const std::vector<std::string>& filenames,
101 : number minLayerHeight);
102 :
103 : // void save_to_files(const char* filenamePrefix);
104 : // void save_rel_to_glob_table_to_files(const char* filenamePrefix);
105 :
106 :
107 : void resize(size_t newSize);
108 :
109 : size_t size() const {return m_layers.size();}
110 : size_t num_layers() const {return m_layers.size();}
111 : bool empty() const {return m_layers.empty();}
112 :
113 : layer_t& operator[] (size_t i) {return *m_layers[i];}
114 : const layer_t& operator[] (size_t i) const {return *m_layers[i];}
115 :
116 : layer_t& layer (size_t i) {return *m_layers[i];}
117 : const layer_t& layer (size_t i) const {return *m_layers[i];}
118 :
119 : const layer_t& top () const {return *m_layers.back();}
120 :
121 0 : Heightfield& heightfield (size_t i) {return layer(i).heightfield;}
122 0 : const Heightfield& heightfield (size_t i) const {return layer(i).heightfield;}
123 :
124 0 : void set_min_height (size_t i, number h) {layer(i).minHeight = h;}
125 0 : number min_height (size_t i) const {return layer(i).minHeight;}
126 :
127 :
128 : /// invalidates cells in lower levels which are too close to valid cells in higher levels
129 : void invalidate_flat_cells();
130 :
131 : /// invalidates cells that belong to a small lense regarding its horizontal area
132 : void invalidate_small_lenses(number minArea);
133 :
134 : /// removes small holes by expanding the layer in those regions to the specified height
135 : void remove_small_holes(number maxArea);
136 :
137 : /// sets invalid or flat cells to the value of the corresponding cell in the level above
138 : /** This method is somehow antithetical to 'invalidate_flat_cells', since it reassigns
139 : * values to invalid cells which are shadowed by valid cells.*/
140 : void snap_cells_to_higher_layers();
141 :
142 : /// eliminates invalid cells by filling those cells with averages of neighboring valid cells
143 : void eliminate_invalid_cells();
144 :
145 : /// smoothens the values in each layer by averaging with neighboured values
146 : void blur_layers (number alpha, size_t numIterations);
147 :
148 : /// finds the first valid value at the given x-y-coordinate starting at the specified layer moving downwards.
149 : /** returns a pair containing the layer-index in which the valid value was found (first)
150 : * and the value at the given coordinate in that layer. Returns -1 if no such layer
151 : * was found.*/
152 : std::pair<int, number> trace_line_down (const vector2& c, size_t firstLayer) const;
153 :
154 : /// finds the first valid value at the given x-y-coordinate starting at the specified layer moving downwards.
155 : /** returns a pair containing the layer-index in which the valid value was found (first)
156 : * and the value at the given coordinate in that layer. Returns -1 if no such layer
157 : * was found.*/
158 : std::pair<int, number> trace_line_up (const vector2& c, size_t firstLayer) const;
159 :
160 : /// returns an index-pair of the layers above and below the specified point
161 : /** If there is no layer above or below, the associated component of the
162 : * returned is set to -1.*/
163 : std::pair<int, int> get_layer_indices (const vector3& c) const;
164 :
165 : /// transforms a relative height to an absolute height for a given x-y-coordinate.
166 : /** relative height is a value between 0 and #numLayers-1. if it is an integer
167 : * value the returned height will match the height of the associated layer.
168 : * If not or if the value would be invalid, it the non-integer fraction is used
169 : * to interpolate between the next higher and the next lower level.
170 : *
171 : * \note if 'construct_relative_to_global_height_table' was called before
172 : * this method, a more sophisticated height-value computation is
173 : * performed in inner invalid cells. This is especially useful if
174 : * a pure prism geometry was constructed from layers with holes.*/
175 : number relative_to_global_height (const vector2& c, number relHeight) const;
176 :
177 : /// transforms a relative height to an absolute height for a given x-y-coordinate.
178 : /** This method works similar to the original 'relative_to_global_height',
179 : * however, it always works on the orignal layer data and follows an equal
180 : * distances approach for no-data-cells. It thus ignores tables constructed
181 : * through 'construct_relative_to_global_height_table'. If the former method
182 : * hasn't been called or if 'invalidate_relative_to_global_height_table' has been
183 : * called, this method will do the same as 'relative_to_global_height'.*/
184 : number relative_to_global_height_simple (const vector2& c, number relHeight) const;
185 :
186 : /// Prepares a table for better 'relative_to_global_height' values in invalid inner regions.
187 : /** Constructs a table in which interior no-data-values are replaced by a
188 : * a relaxed value, computed through smoothing the relative distances to
189 : * the upper and lower layers of local neighbor cells.
190 : *
191 : * \note if the underlying layers have been changed or new ones have been added,
192 : * this method has to be called again to reflect those changes in
193 : * the constructed table.
194 : *
195 : * \param iterations the number of relaxation iterations (e.g. 1000)
196 : * \param alpha the relative amount of how much a value may change in
197 : * each iteration (between 0 and 1, e.g. 0.5)*/
198 : void construct_relative_to_global_height_table (size_t iterations, number alpha);
199 :
200 : /// invalidates the table construced by 'construct_relative_to_global_height_table'
201 : /** Use this method if you want to make sure that no special table is used
202 : * during 'relative_to_global_height'.*/
203 : void invalidate_relative_to_global_height_table ();
204 :
205 : private:
206 : /// returns dist(middle, upper, ix, iy) / dist(lower, upper, ix, iy)
207 : /** If dist(lower, upper, ix, iy) == 0, 0 is returned.*/
208 : number upper_lower_dist_relation ( Field<number>&lower,
209 : Field<number>& middle,
210 : Field<number>& upper,
211 : size_t ix,
212 : size_t iy);
213 :
214 : // BEGIN SERIALIZATION
215 : friend class boost::serialization::access;
216 :
217 : template <class Archive>
218 0 : void save( Archive& ar, const unsigned int version) const
219 : {
220 0 : size_t numLayers = m_layers.size();
221 : ar & numLayers;
222 0 : for(size_t i = 0; i < m_layers.size(); ++i){
223 : ar & *m_layers[i];
224 : }
225 :
226 0 : size_t numRelToGlob = m_relativeToGlobalHeights.size();
227 : ar & numRelToGlob;
228 0 : for(size_t i = 0; i < m_relativeToGlobalHeights.size(); ++i){
229 : ar & *m_relativeToGlobalHeights[i];
230 : }
231 0 : }
232 :
233 : template <class Archive>
234 0 : void load( Archive& ar, const unsigned int version)
235 : {
236 0 : size_t numLayers = 0;
237 : ar & numLayers;
238 0 : m_layers.resize(numLayers);
239 0 : for(size_t i = 0; i < numLayers; ++i){
240 0 : m_layers[i] = make_sp(new layer_t);
241 : ar & *m_layers[i];
242 : }
243 :
244 0 : size_t numRelToGlob = 0;
245 : ar & numRelToGlob;
246 0 : m_relativeToGlobalHeights.resize(numRelToGlob);
247 0 : for(size_t i = 0; i < numRelToGlob; ++i){
248 0 : m_relativeToGlobalHeights[i] = make_sp(new Heightfield);
249 : ar & *m_relativeToGlobalHeights[i];
250 : }
251 0 : }
252 :
253 : BOOST_SERIALIZATION_SPLIT_MEMBER()
254 : // END SERIALIZATION
255 :
256 : std::vector<SmartPtr<layer_t> > m_layers;
257 :
258 : /** In this array, interior pixels with no-data-value are replaced by a
259 : * relaxed value, computed through smoothing the relative distances to
260 : * the upper and lower layers of local neighbor cells.
261 : * This array has to be constructed explicitely through a call to
262 : * 'construct_relative_to_global_height_table'.*/
263 : std::vector<SmartPtr<Heightfield> > m_relativeToGlobalHeights;
264 : };
265 :
266 : typedef SmartPtr<RasterLayers> SPRasterLayers;
267 :
268 : }// end of namespace
269 :
270 : #endif //__H__UG_raster_layer_util
|
