Line data Source code
1 : /*
2 : * Copyright (c) 2016: 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_impl
34 : #define __H__UG_raster_impl
35 :
36 : #include <limits>
37 : #include <cstring>
38 : #include <algorithm>
39 : #include <fstream>
40 : #include "common/error.h"
41 : #include "common/util/file_util.h"
42 : #include "common/util/string_util.h"
43 : #include "raster_kernels.h"
44 :
45 : namespace ug{
46 :
47 : ////////////////////////////////////////////////////////////////////////////////
48 : // Raster::MultiIndex
49 :
50 : template <class T, int TDIM>
51 0 : Raster<T, TDIM>::MultiIndex::
52 : MultiIndex()
53 : {}
54 :
55 : template <class T, int TDIM>
56 0 : Raster<T, TDIM>::MultiIndex::
57 : MultiIndex(size_t i)
58 : {
59 : set(i);
60 : }
61 :
62 : template <class T, int TDIM>
63 : int Raster<T, TDIM>::MultiIndex::
64 : dim () const
65 : {
66 : return TDIM;
67 : }
68 :
69 : template <class T, int TDIM>
70 : void Raster<T, TDIM>::MultiIndex::
71 : set (size_t i)
72 : {
73 0 : for(int d = 0; d < TDIM; ++d)
74 0 : m_ind[d] = i;
75 : }
76 :
77 : template <class T, int TDIM>
78 : size_t& Raster<T, TDIM>::MultiIndex::
79 : operator[] (int d)
80 : {
81 : return m_ind[d];
82 : }
83 :
84 : template <class T, int TDIM>
85 : size_t Raster<T, TDIM>::MultiIndex::
86 : operator[] (int d) const
87 : {
88 0 : return m_ind[d];
89 : }
90 :
91 :
92 : ////////////////////////////////////////////////////////////////////////////////
93 : // Raster::Coordinate
94 :
95 : template <class T, int TDIM>
96 0 : Raster<T, TDIM>::Coordinate::
97 : Coordinate()
98 : {}
99 :
100 : template <class T, int TDIM>
101 0 : Raster<T, TDIM>::Coordinate::
102 : Coordinate(number c)
103 : {
104 : set(c);
105 : }
106 :
107 : template <class T, int TDIM>
108 : Raster<T, TDIM>::Coordinate::
109 : Coordinate(const MathVector<TDIM, number>& v)
110 : {
111 : for(int d = 0; d < TDIM; ++d)
112 : m_coord[d] = v[d];
113 : }
114 :
115 : template <class T, int TDIM>
116 : int Raster<T, TDIM>::Coordinate::
117 : dim () const
118 : {
119 : return TDIM;
120 : }
121 :
122 : template <class T, int TDIM>
123 : void Raster<T, TDIM>::Coordinate::
124 : set (number c)
125 : {
126 0 : for(int d = 0; d < TDIM; ++d)
127 0 : m_coord[d] = c;
128 : }
129 :
130 :
131 : template <class T, int TDIM>
132 : number& Raster<T, TDIM>::Coordinate::
133 : operator[] (int d)
134 : {
135 : return m_coord[d];
136 : }
137 :
138 : template <class T, int TDIM>
139 : number Raster<T, TDIM>::Coordinate::
140 : operator[] (int d) const
141 : {
142 0 : return m_coord[d];
143 : }
144 :
145 : template <class T, int TDIM>
146 : typename Raster<T, TDIM>::Coordinate& Raster<T, TDIM>::Coordinate::
147 : operator+= (const Coordinate& c)
148 : {
149 : for(int d = 0; d < TDIM; ++d)
150 : m_coord[d] += c[d];
151 : }
152 :
153 : template <class T, int TDIM>
154 : typename Raster<T, TDIM>::Coordinate& Raster<T, TDIM>::Coordinate::
155 : operator-= (const Coordinate& c)
156 : {
157 : for(int d = 0; d < TDIM; ++d)
158 : m_coord[d] -= c[d];
159 : }
160 :
161 : template <class T, int TDIM>
162 : typename Raster<T, TDIM>::Coordinate& Raster<T, TDIM>::Coordinate::
163 : operator*= (number s)
164 : {
165 : for(int d = 0; d < TDIM; ++d)
166 : m_coord[d] *= s;
167 : }
168 :
169 :
170 :
171 : ////////////////////////////////////////////////////////////////////////////////
172 : // Raster - public
173 :
174 : template <class T, int TDIM>
175 0 : Raster<T, TDIM>::
176 : Raster () :
177 0 : m_data(NULL),
178 : m_numNodes(0),
179 : m_selNode(0),
180 : m_minCorner(0),
181 : m_extension(1),
182 : m_cellExtension(1),
183 : m_cursor(0),
184 : m_numNodesTotal(0),
185 0 : m_noDataValue(std::numeric_limits<T>::max())
186 : {
187 : update_num_nodes_total();
188 : update_cell_extension();
189 0 : }
190 :
191 :
192 : template <class T, int TDIM>
193 : Raster<T, TDIM>::
194 : Raster (const Raster<T, TDIM>& raster) :
195 : m_data(NULL),
196 : m_numNodes(raster.m_numNodes),
197 : m_selNode(raster.m_selNode),
198 : m_minCorner(raster.m_minCorner),
199 : m_extension(raster.m_extension),
200 : m_cellExtension(raster.m_cellExtension),
201 : m_cursor(raster.m_cursor),
202 : m_numNodesTotal(raster.m_numNodesTotal),
203 : m_noDataValue(raster.m_noDataValue)
204 : {
205 : update_num_nodes_total();
206 : update_cell_extension();
207 :
208 : create();
209 :
210 : memcpy(m_data, raster.m_data, num_nodes_total() * sizeof(T));
211 : }
212 :
213 : template <class T, int TDIM>
214 : Raster<T, TDIM>::
215 : Raster (const MultiIndex& numNodes) :
216 : m_data(NULL),
217 : m_numNodes(numNodes),
218 : m_selNode(0),
219 : m_minCorner(0),
220 : m_extension(1),
221 : m_cellExtension(1),
222 : m_cursor(0),
223 : m_numNodesTotal(0),
224 : m_noDataValue(std::numeric_limits<T>::max())
225 : {
226 : update_num_nodes_total();
227 :
228 : for(int d = 0; d < TDIM; ++d)
229 : m_extension[d] = numNodes[d] - 1;
230 :
231 : update_cell_extension();
232 :
233 : create();
234 : }
235 :
236 : template <class T, int TDIM>
237 : Raster<T, TDIM>::
238 : Raster (const MultiIndex& numNodes,
239 : const Coordinate& extension,
240 : const Coordinate& minCorner) :
241 : m_data(NULL),
242 : m_numNodes(numNodes),
243 : m_selNode(0),
244 : m_minCorner(minCorner),
245 : m_extension(extension),
246 : m_cellExtension(1),
247 : m_cursor(0),
248 : m_numNodesTotal(0),
249 : m_noDataValue(std::numeric_limits<T>::max())
250 : {
251 : update_num_nodes_total();
252 : update_cell_extension();
253 : create();
254 : }
255 :
256 :
257 : template <class T, int TDIM>
258 : Raster<T, TDIM>::
259 : ~Raster ()
260 : {
261 0 : if(m_data)
262 0 : delete[] m_data;
263 0 : }
264 :
265 : template <class T, int TDIM>
266 : Raster<T, TDIM>& Raster<T, TDIM>::
267 : operator= (const Raster& raster)
268 : {
269 : m_numNodes = raster.m_numNodes;
270 : m_selNode = raster.m_selNode;
271 : m_minCorner = raster.m_minCorner;
272 : m_extension = raster.m_extension;
273 : m_cursor = raster.m_cursor;
274 : m_noDataValue = raster.m_noDataValue;
275 :
276 : update_num_nodes_total();
277 : update_cell_extension();
278 :
279 : create();
280 :
281 : memcpy(m_data, raster.m_data, num_nodes_total() * sizeof(T));
282 :
283 : return *this;
284 : }
285 :
286 : template <class T, int TDIM>
287 0 : int Raster<T, TDIM>::
288 : dim () const
289 : {
290 0 : return TDIM;
291 : }
292 :
293 : template <class T, int TDIM>
294 0 : void Raster<T, TDIM>::
295 : set_num_nodes (int dim, size_t num)
296 : {
297 0 : m_numNodes[dim] = num;
298 : update_num_nodes_total();
299 : update_cell_extension(dim);
300 0 : }
301 :
302 : template <class T, int TDIM>
303 0 : void Raster<T, TDIM>::
304 : set_num_nodes (const typename Raster<T, TDIM>::MultiIndex& mi)
305 : {
306 0 : m_numNodes = mi;
307 : update_num_nodes_total();
308 : update_cell_extension();
309 0 : }
310 :
311 : template <class T, int TDIM>
312 : size_t Raster<T, TDIM>::
313 : num_nodes_total () const
314 : {
315 0 : return m_numNodesTotal;
316 : }
317 :
318 : template <class T, int TDIM>
319 0 : size_t Raster<T, TDIM>::
320 : num_nodes (int dim) const
321 : {
322 0 : return m_numNodes[dim];
323 : }
324 :
325 : template <class T, int TDIM>
326 : const typename Raster<T, TDIM>::MultiIndex& Raster<T, TDIM>::
327 : num_nodes () const
328 : {
329 : return m_numNodes;
330 : }
331 :
332 :
333 : template <class T, int TDIM>
334 0 : void Raster<T, TDIM>::
335 : create ()
336 : {
337 0 : if(m_data){
338 0 : delete[] m_data;
339 0 : m_data = NULL;
340 : }
341 :
342 : update_num_nodes_total(); // this isn't strictly necessary if everything works right.
343 :
344 : const size_t num = num_nodes_total();
345 0 : if(num){
346 0 : m_data = new T[num];
347 : }
348 0 : }
349 :
350 :
351 : template <class T, int TDIM>
352 : T& Raster<T, TDIM>::
353 : node_value (const MultiIndex& mi)
354 : {
355 0 : return m_data[data_index(mi)];
356 : }
357 :
358 : template <class T, int TDIM>
359 : T Raster<T, TDIM>::
360 : node_value (const MultiIndex& mi) const
361 : {
362 0 : return m_data[data_index(mi)];
363 : }
364 :
365 :
366 : template <class T, int TDIM>
367 0 : void Raster<T, TDIM>::
368 : set_min_corner (int dim, number coord)
369 : {
370 0 : m_minCorner[dim] = coord;
371 0 : }
372 :
373 : template <class T, int TDIM>
374 : void Raster<T, TDIM>::
375 : set_min_corner (const Coordinate& coord)
376 : {
377 0 : m_minCorner = coord;
378 : }
379 :
380 : template <class T, int TDIM>
381 : const typename Raster<T, TDIM>::Coordinate& Raster<T, TDIM>::
382 : min_corner () const
383 : {
384 : return m_minCorner;
385 : }
386 :
387 : template <class T, int TDIM>
388 0 : number Raster<T, TDIM>::
389 : min_corner (int dim) const
390 : {
391 0 : return m_minCorner[dim];
392 : }
393 :
394 : template <class T, int TDIM>
395 0 : void Raster<T, TDIM>::
396 : set_extension (int dim, number ext)
397 : {
398 0 : m_extension[dim] = ext;
399 : update_cell_extension(dim);
400 0 : }
401 :
402 : template <class T, int TDIM>
403 0 : void Raster<T, TDIM>::
404 : set_extension (const Coordinate& ext)
405 : {
406 0 : m_extension = ext;
407 : update_cell_extension();
408 0 : }
409 :
410 : template <class T, int TDIM>
411 : const typename Raster<T, TDIM>::Coordinate& Raster<T, TDIM>::
412 : extension () const
413 : {
414 : return m_extension;
415 : }
416 :
417 : template <class T, int TDIM>
418 0 : number Raster<T, TDIM>::
419 : extension (int dim) const
420 : {
421 0 : return m_extension[dim];
422 : }
423 :
424 :
425 : template <class T, int TDIM>
426 0 : T Raster<T, TDIM>::
427 : interpolate (const Coordinate& coord, int order) const
428 : {
429 0 : switch(order){
430 : case 0: {
431 : MultiIndex mi;
432 0 : for(size_t d = 0; d < TDIM; ++d)
433 : {
434 0 : mi[d] = static_cast<int>(0.5 + (coord[d] - m_minCorner[d]) / m_cellExtension[d]);
435 : if(mi[d] < 0) mi[d] = 0;
436 0 : else if(mi[d] >= num_nodes(d)) mi[d] = num_nodes(d) - 1;
437 : }
438 : return node_value(mi);
439 : } break;
440 :
441 : case 1: {
442 : MultiIndex mi;
443 : Coordinate lc;
444 0 : for(size_t d = 0; d < TDIM; ++d)
445 : {
446 0 : mi[d] = static_cast<int>((coord[d] - m_minCorner[d]) / m_cellExtension[d]);
447 : if(mi[d] < 0){
448 : mi[d] = 0;
449 : lc[d] = 0;
450 : }
451 0 : else if(mi[d] + 1 >= num_nodes(d)){
452 0 : mi[d] = num_nodes(d) - 2;
453 0 : lc[d] = 1;
454 : }
455 : else{
456 0 : lc[d] = ( coord[d]
457 0 : - ((number)mi[d] * m_cellExtension[d]
458 0 : + m_minCorner[d]))
459 0 : / m_cellExtension[d];
460 : }
461 : }
462 :
463 0 : return interpolate_linear(mi, lc);
464 : } break;
465 :
466 0 : default:
467 0 : UG_THROW("Raster::interpolate(): Unsupported interpolation order: " << order);
468 : }
469 : return m_noDataValue;
470 : }
471 :
472 :
473 : template <class T, int TDIM>
474 0 : void Raster<T, TDIM>::
475 : set_no_data_value(T val)
476 : {
477 0 : m_noDataValue = val;
478 0 : }
479 :
480 : template <class T, int TDIM>
481 0 : T Raster<T, TDIM>::
482 : no_data_value() const
483 : {
484 0 : return m_noDataValue;
485 : }
486 :
487 :
488 : template <class T, int TDIM>
489 0 : void Raster<T, TDIM>::
490 : blur(T alpha, size_t iterations)
491 : {
492 : raster_kernels::Blur<T, TDIM> blurKernel (alpha);
493 : const MultiIndex start(0);
494 :
495 0 : for(size_t iiter = 0; iiter < iterations; ++iiter)
496 : run_on_all (blurKernel);
497 0 : }
498 :
499 :
500 : template <class T, int TDIM>
501 : template <class TKernel>
502 : typename TKernel::result_t Raster<T, TDIM>::
503 : run_on_all()
504 : {
505 : TKernel kernel;
506 : run_on_all (MultiIndex(0), kernel, TDIM - 1);
507 : return kernel.result();
508 : }
509 :
510 : template <class T, int TDIM>
511 : template <class TKernel>
512 : void Raster<T, TDIM>::
513 : run_on_all(TKernel& kernel)
514 : {
515 0 : run_on_all (MultiIndex(0), kernel, TDIM - 1);
516 : }
517 :
518 :
519 : template <class T, int TDIM>
520 : template <class TKernel>
521 0 : void Raster<T, TDIM>::
522 : run_on_all(const MultiIndex& start, TKernel& kernel, int curDim)
523 : {
524 0 : if(curDim > 0) {
525 : const size_t numNodes = num_nodes(curDim);
526 0 : for(MultiIndex cur = start; cur[curDim] < numNodes; ++cur[curDim]){
527 0 : run_on_all(cur, kernel, curDim - 1);
528 : }
529 : }
530 : else {
531 : const size_t numNodes = num_nodes(0);
532 0 : for(MultiIndex cur = start; cur[0] < numNodes; ++cur[0]){
533 0 : kernel (*this, cur);
534 : }
535 : }
536 0 : }
537 :
538 :
539 : template <class T, int TDIM>
540 : template <class TKernel>
541 : typename TKernel::result_t Raster<T, TDIM>::
542 : run_on_nbrs(const MultiIndex& center)
543 : {
544 : TKernel kernel;
545 0 : run_on_nbrs(center, kernel, TDIM - 1);
546 : return kernel.result();
547 : }
548 :
549 :
550 : template <class T, int TDIM>
551 : template <class TKernel>
552 : void Raster<T, TDIM>::
553 : run_on_nbrs(const MultiIndex& center, TKernel& kernel)
554 : {
555 : run_on_nbrs(center, kernel, TDIM - 1);
556 : }
557 :
558 :
559 : template <class T, int TDIM>
560 : template <class TKernel>
561 0 : void Raster<T, TDIM>::
562 : run_on_nbrs(const MultiIndex& center, TKernel& kernel, int curDim)
563 : {
564 0 : if(curDim > 0)
565 0 : run_on_nbrs(center, kernel, curDim - 1);
566 :
567 0 : if(center[curDim] > 0){
568 0 : MultiIndex c = center;
569 0 : --c[curDim];
570 0 : kernel(*this, c);
571 : }
572 :
573 0 : if(center[curDim] + 1 < num_nodes(curDim)){
574 0 : MultiIndex c = center;
575 0 : ++c[curDim];
576 0 : kernel(*this, c);
577 : }
578 0 : }
579 :
580 :
581 : template <class T, int TDIM>
582 0 : void Raster<T, TDIM>::
583 : select_node (int dim, size_t index)
584 : {
585 0 : m_selNode[dim] = index;
586 0 : }
587 :
588 : template <class T, int TDIM>
589 : void Raster<T, TDIM>::
590 : select_node (const MultiIndex& mi)
591 : {
592 : m_selNode = mi;
593 : }
594 :
595 :
596 : template <class T, int TDIM>
597 0 : void Raster<T, TDIM>::
598 : set_selected_node_value (T val)
599 : {
600 0 : node_value(m_selNode) = val;
601 0 : }
602 :
603 : template <class T, int TDIM>
604 0 : T Raster<T, TDIM>::
605 : selected_node_value () const
606 : {
607 0 : return node_value(m_selNode);
608 : }
609 :
610 :
611 : template <class T, int TDIM>
612 0 : void Raster<T, TDIM>::
613 : set_cursor (int dim, number coord)
614 : {
615 0 : m_cursor[dim] = coord;
616 0 : }
617 :
618 : template <class T, int TDIM>
619 : void Raster<T, TDIM>::
620 : set_cursor (const Coordinate& coord)
621 : {
622 : m_cursor = coord;
623 : }
624 :
625 :
626 : template <class T, int TDIM>
627 0 : T Raster<T, TDIM>::
628 : interpolate_at_cursor (int order) const
629 : {
630 0 : return interpolate(m_cursor, order);
631 : }
632 :
633 :
634 : template <class T, int TDIM>
635 0 : void Raster<T, TDIM>::
636 : load_from_asc (const char* filename)
637 : {
638 : using namespace std;
639 :
640 : #define LFA_ERR_WHERE "Error in Raster::load_from_asc('" << filename << "'): "
641 : // this macro helps with error-checks for bad dimensions
642 : #define LFA_CHECK_DIM(d, line)\
643 : UG_COND_THROW(d >= TDIM, LFA_ERR_WHERE << "Bad dimension '" << d <<\
644 : "' in line " << line << " of file " << filename << "," <<\
645 : "while trying to read a " << TDIM << "d raster.");
646 :
647 0 : std::string fullFileName = FindFileInStandardPaths(filename);
648 0 : UG_COND_THROW(fullFileName.empty(),
649 : LFA_ERR_WHERE << "Couldn't find the specified file in any of the standard paths.");
650 :
651 0 : ifstream in(fullFileName.c_str());
652 0 : UG_COND_THROW(!in, LFA_ERR_WHERE << "Couldn't access file.");
653 :
654 : MultiIndex numNodes(0);
655 :
656 : // indicate whether minCoord was specified as cell-center
657 : MultiIndex minCoordIsCenter(0);
658 :
659 : Coordinate minCoord(0);
660 : Coordinate cellSize(0);
661 :
662 : T noDataValue = T();
663 :
664 : // parse header
665 : // the header lenght varies between different asc files. This depends on the
666 : // dimension and whether equlateral cells are specified or not, i.e., whether
667 : // 'cellsize' or whether 'xcellsize', 'ycellsize', ... was specified.
668 : // We're trying to guess the correct length here.
669 : int headerLen = 0;
670 : if(TDIM == 1) headerLen = 4;
671 : else if(TDIM == 2) headerLen = 6;
672 : else if(TDIM == 3) headerLen = 8;
673 : else{
674 : UG_THROW("Raster::load_from_asc only supports 1, 2, and 3 dimensions\n");
675 : }
676 :
677 0 : for(int i = 0; i < headerLen; ++i){
678 : string name;
679 : double value;
680 0 : in >> name >> value;
681 0 : UG_COND_THROW(!in, LFA_ERR_WHERE <<
682 : "Couldn't parse expected name-value pair in row " << i);
683 :
684 0 : name = ToLower(name);
685 :
686 0 : if(name.compare("ncols") == 0){
687 : LFA_CHECK_DIM(0, i);
688 0 : numNodes[0] = (int)value;
689 : }
690 0 : else if(name.compare("nrows") == 0){
691 0 : LFA_CHECK_DIM(1, i);
692 0 : numNodes[1] = (int)value;
693 : }
694 0 : else if(name.compare("nstacks") == 0){
695 0 : LFA_CHECK_DIM(2, i);
696 0 : numNodes[2] = (int)value;
697 : }
698 :
699 0 : else if(name.compare("xllcenter") == 0){
700 : LFA_CHECK_DIM(0, i);
701 0 : minCoord[0] = value;
702 0 : minCoordIsCenter[0] = 1;
703 : }
704 0 : else if(name.compare("yllcenter") == 0){
705 0 : LFA_CHECK_DIM(1, i);
706 0 : minCoord[1] = value;
707 0 : minCoordIsCenter[1] = 1;
708 : }
709 0 : else if(name.compare("zllcenter") == 0){
710 0 : LFA_CHECK_DIM(2, i);
711 0 : minCoord[2] = value;
712 0 : minCoordIsCenter[2] = 1;
713 : }
714 0 : else if(name.compare("xllcorner") == 0){
715 : LFA_CHECK_DIM(0, i);
716 0 : minCoord[0] = value;
717 : }
718 0 : else if(name.compare("yllcorner") == 0){
719 0 : LFA_CHECK_DIM(1, i);
720 0 : minCoord[1] = value;
721 : }
722 0 : else if(name.compare("zllcorner") == 0){
723 0 : LFA_CHECK_DIM(2, i);
724 0 : minCoord[2] = value;
725 : }
726 :
727 0 : else if(name.compare("cellsize") == 0){
728 0 : for(int d = 0; d < TDIM; ++d)
729 0 : cellSize[d] = value;
730 : }
731 :
732 0 : else if(name.compare("xcellsize") == 0){
733 : LFA_CHECK_DIM(0, i);
734 : // we have to read additional cell-sizes for the other dimensions
735 0 : headerLen += (TDIM - 1);
736 0 : cellSize[0] = value;
737 : }
738 0 : else if(name.compare("ycellsize") == 0){
739 0 : LFA_CHECK_DIM(1, i);
740 0 : cellSize[1] = value;
741 : }
742 0 : else if(name.compare("zcellsize") == 0){
743 0 : LFA_CHECK_DIM(2, i);
744 0 : cellSize[2] = value;
745 : }
746 :
747 0 : else if(name.compare("nodata_value") == 0){
748 0 : noDataValue = value;
749 : }
750 :
751 : else{
752 0 : UG_THROW(LFA_ERR_WHERE << "unknown identifier in header: " << name);
753 : }
754 : }
755 :
756 0 : for(int d = 0; d < TDIM; ++d){
757 0 : if(minCoordIsCenter[d])
758 0 : minCoord[d] -= 0.5 * cellSize[d];
759 : }
760 :
761 : // check validity
762 0 : for(int d = 0; d < TDIM; ++d){
763 0 : UG_COND_THROW(numNodes[d] == 0, LFA_ERR_WHERE << "Num nodes may not be 0 for dim " << d);
764 0 : UG_COND_THROW(cellSize[d] <= 0, LFA_ERR_WHERE << "cell-size must be bigger than 0 for dim " << d);
765 : }
766 :
767 0 : set_num_nodes(numNodes);
768 : set_min_corner(minCoord);
769 0 : Coordinate extension = cellSize;
770 0 : for(int d = 0; d < TDIM; ++d)
771 0 : extension[d] *= (number)(numNodes[d] - 1);
772 0 : set_extension(extension);
773 : set_no_data_value(noDataValue);
774 :
775 0 : create();
776 :
777 : // parse values
778 : // y and z are inverted
779 0 : size_t num[3] = {0, 1, 1};
780 0 : for(size_t i = 0; i < TDIM; ++i)
781 0 : num[i] = m_numNodes[i];
782 :
783 0 : for(size_t iz = 0; iz < num[2]; ++iz){
784 0 : for(size_t iy = 0; iy < num[1]; ++iy){
785 0 : for(size_t ix = 0; ix < num[0]; ++ix)
786 : {
787 0 : const size_t ty = num[1] - 1 - iy;
788 0 : const size_t tz = num[2] - 1 - iz;
789 0 : in >> m_data[ix + num[0] * (ty + num[1] * tz)];
790 0 : UG_COND_THROW(!in, LFA_ERR_WHERE << "Couldn't read value for at ("
791 : << ix << ", " << iy << ", " << iz << ")");
792 : }
793 : }
794 : }
795 0 : }
796 :
797 : template <class T, int TDIM>
798 0 : void Raster<T, TDIM>::
799 : save_to_asc (const char* filename) const
800 : {
801 : using namespace std;
802 : #define STA_ERR_WHERE "Error in Raster::save_to_asc('" << filename << "'): "
803 :
804 0 : UG_COND_THROW(!m_data, STA_ERR_WHERE << "Can't write an unitinialized raster."
805 : "Please call 'create' or 'load_from_asc' first.");
806 :
807 0 : ofstream out(filename);
808 0 : UG_COND_THROW(!out, STA_ERR_WHERE << "Couldn't open file for writing.");
809 :
810 : out << "ncols " << num_nodes(0) << endl;
811 : if(TDIM > 1)
812 : out << "nrows " << num_nodes(1) << endl;
813 : if(TDIM > 2)
814 : out << "nstacks " << num_nodes(2) << endl;
815 :
816 : out << "xllcorner " << setprecision(16) << min_corner(0) << endl;
817 : if(TDIM > 1)
818 : out << "yllcorner " << setprecision(16) << min_corner(1) << endl;
819 : if(TDIM > 2)
820 : out << "zllcorner " << setprecision(16) << min_corner(2) << endl;
821 :
822 : bool equlateralCells = true;
823 0 : for(int d = 1; d < TDIM; ++d){
824 0 : if(m_cellExtension[d] != m_cellExtension[0]){
825 : equlateralCells = false;
826 : break;
827 : }
828 : }
829 :
830 0 : if(equlateralCells)
831 : out << "cellsize " << m_cellExtension[0] << endl;
832 : else{
833 : out << "xcellsize " << m_cellExtension[0] << endl;
834 : if(TDIM > 1)
835 : out << "ycellsize " << m_cellExtension[1] << endl;
836 : if(TDIM > 2)
837 : out << "zcellsize " << m_cellExtension[2] << endl;
838 : }
839 :
840 : out << "NODATA_value " << no_data_value() << endl;
841 :
842 : // write values
843 : // y and z are inverted
844 0 : size_t num[3] = {0, 1, 1};
845 0 : for(size_t i = 0; i < TDIM; ++i)
846 0 : num[i] = num_nodes(i);
847 :
848 0 : for(size_t iz = 0; iz < num[2]; ++iz){
849 0 : for(size_t iy = 0; iy < num[1]; ++iy){
850 0 : for(size_t ix = 0; ix < num[0]; ++ix)
851 : {
852 0 : if(ix > 0)
853 0 : out << " ";
854 0 : const size_t ty = num[1] - 1 - iy;
855 0 : const size_t tz = num[2] - 1 - iz;
856 0 : out << m_data[ix + num[0] * (ty + num[1] * tz)];
857 : }
858 : out << endl;
859 : }
860 : out << endl;
861 : }
862 : out << endl;
863 0 : }
864 :
865 : ////////////////////////////////////////////////////////////////////////////////
866 : // Raster - private
867 : template <class T, int TDIM>
868 : size_t Raster<T, TDIM>::
869 : data_index (const MultiIndex& mi, int curDim, size_t curVal) const
870 : {
871 0 : if(curDim == 0)
872 0 : return curVal + mi[0];
873 : else{
874 0 : return data_index(mi, curDim - 1, m_numNodes[curDim - 1] * (mi[curDim] + curVal));
875 : }
876 : }
877 :
878 : template <class T, int TDIM>
879 : void Raster<T, TDIM>::
880 : update_num_nodes_total()
881 : {
882 0 : m_numNodesTotal = 1;
883 0 : for(int d = 0; d < TDIM; ++d)
884 0 : m_numNodesTotal *= num_nodes(d);
885 : }
886 :
887 : template <class T, int TDIM>
888 : void Raster<T, TDIM>::
889 : update_cell_extension()
890 : {
891 0 : for(int d = 0; d < TDIM; ++d)
892 : update_cell_extension(d);
893 : }
894 :
895 : template <class T, int TDIM>
896 : void Raster<T, TDIM>::
897 : update_cell_extension(int dim)
898 : {
899 0 : if(m_numNodes[dim] > 1 && m_extension[dim] > 0)
900 0 : m_cellExtension[dim] = m_extension[dim] / (m_numNodes[dim] - 1);
901 : else
902 0 : m_cellExtension[dim] = 1;
903 : }
904 :
905 :
906 : template <class T, int TDIM>
907 0 : T Raster<T, TDIM>::
908 : interpolate_linear (
909 : const MultiIndex& minNodeInd,
910 : Coordinate& localCoord,
911 : int curDim) const
912 : {
913 0 : if(curDim == 0)
914 0 : return node_value(minNodeInd);
915 :
916 0 : MultiIndex miMax = minNodeInd;
917 0 : miMax[curDim - 1] += 1;
918 :
919 0 : T val0 = interpolate_linear(minNodeInd, localCoord, curDim - 1);
920 0 : T val1 = interpolate_linear(miMax, localCoord, curDim - 1);
921 :
922 : // perform linear interpolation
923 0 : val0 *= (1. - localCoord[curDim - 1]);
924 0 : val1 *= localCoord[curDim - 1];
925 0 : val0 += val1;
926 0 : return val0;
927 : }
928 :
929 : }// end of namespace
930 :
931 : #endif //__H__UG_raster_impl
|
