source: Revenant/marslib/include/mars_grids.h

port/mars-tycoon
Last change on this file was 80a6a52, checked in by Jonathan Neufeld <support@…>, 3 years ago

Get to a compile state for terrain procedural generation
  • Property mode set to 100755
File size: 54.2 KB
Line 
1#ifndef __MARSGRIDS_H__
2#define __MARSGRIDS_H__
3
4#include <exception>
5#include <iterator>
6#include <list>
7#include <limits>
8#include <istream>
9#include <ostream>
10
11#include <math.h>
12#include <memory.h>
13
14#ifdef _OPENMP
15#include<omp.h>
16#endif
17
18#include "coroutines.h"
19#include "mars_util.h"
20#include "mars_calc.h"
21#include "mars_log.h"
22
23namespace mars
24{
25 template <typename T>
26 class GridCoords
27 {
28 public:
29 typedef T Precision;
30
31 T x, y, z;
32
33 inline GridCoords ()
34 : x(0), y(0), z(0) {}
35 inline GridCoords (T x, T y, T z)
36 : x(x), y(y), z(z) {}
37 };
38
39 template <typename T>
40 inline bool operator < (const GridCoords <T> & left, const GridCoords <T> & right)
41 {
42 return left.z < right.z || (left.z == right.z && left.y < right.y || (left.z == right.z && left.y == right.y && left.x < right.x));
43 }
44
45 template <typename T>
46 class GridEx : public std::exception
47 {
48 public:
49 const GridCoords <T> coords;
50
51 inline GridEx (const GridCoords <T> & coords)
52 : coords(coords) {}
53 };
54
55 struct HexGrid_TAG {};
56 struct BasicGrid_TAG {};
57
58 template <typename T>
59 inline GridCoords <T> tx2basic (HexGrid_TAG, const GridCoords <T> & coordsHexa)
60 {
61 return GridCoords <T>
62 (
63 coordsHexa.x / 2,
64 coordsHexa.y / 2,
65 coordsHexa.z
66 );
67 }
68
69 template <typename T>
70 inline GridCoords <T> tx2basic (BasicGrid_TAG, const GridCoords <T> & coords)
71 {
72 return coords;
73 }
74
75 template <typename T>
76 inline GridCoords <T> txFbasic (HexGrid_TAG, const GridCoords <T> & coords)
77 {
78 return GridCoords <T>
79 (
80 coords.x * 2 + (coords.y + coords.z) % 2,
81 coords.y * 2 + (coords.z % 2),
82 coords.z
83 );
84 }
85
86 template <typename T>
87 inline GridCoords <T> txFbasic (BasicGrid_TAG, const GridCoords <T> & coords)
88 {
89 return coords;
90 }
91
92 template <typename T>
93 void populate_neighbors (HexGrid_TAG, std::list <GridCoords <T> > & list)
94 {
95 list.push_back(GridCoords<T> (-2, -1, 0));
96 list.push_back(GridCoords<T> (-2, +1, 0));
97 list.push_back(GridCoords<T> (+2, -1, 0));
98 list.push_back(GridCoords<T> (+2, +1, 0));
99 list.push_back(GridCoords<T> (0, -2, 0));
100 list.push_back(GridCoords<T> (0, +2, 0));
101 list.push_back(GridCoords<T> (0, -2, -1));
102 list.push_back(GridCoords<T> (0, +2, -1));
103 list.push_back(GridCoords<T> (-2, 0, -1));
104 list.push_back(GridCoords<T> (+2, 0, -1));
105 list.push_back(GridCoords<T> (0, -2, +1));
106 list.push_back(GridCoords<T> (0, +2, +1));
107 list.push_back(GridCoords<T> (-2, 0, +1));
108 list.push_back(GridCoords<T> (+2, 0, +1));
109 }
110
111 template <typename T>
112 void populate_neighbors (BasicGrid_TAG, std::list <GridCoords <T> > & list)
113 {
114 list.push_back(GridCoords<T> (-1, -1, -1));
115 list.push_back(GridCoords<T> (-1, 0, -1));
116 list.push_back(GridCoords<T> (-1, +1, -1));
117 list.push_back(GridCoords<T> (0, -1, -1));
118 list.push_back(GridCoords<T> (0, 0, -1));
119 list.push_back(GridCoords<T> (0, +1, -1));
120 list.push_back(GridCoords<T> (+1, -1, -1));
121 list.push_back(GridCoords<T> (+1, 0, -1));
122 list.push_back(GridCoords<T> (+1, +1, -1));
123
124
125 list.push_back(GridCoords<T> (-1, -1, 0));
126 list.push_back(GridCoords<T> (-1, 0, 0));
127 list.push_back(GridCoords<T> (-1, +1, 0));
128 list.push_back(GridCoords<T> (0, -1, 0));
129
130 list.push_back(GridCoords<T> (0, +1, 0));
131 list.push_back(GridCoords<T> (+1, -1, 0));
132 list.push_back(GridCoords<T> (+1, 0, 0));
133 list.push_back(GridCoords<T> (+1, +1, 0));
134
135
136 list.push_back(GridCoords<T> (-1, -1, +1));
137 list.push_back(GridCoords<T> (-1, 0, +1));
138 list.push_back(GridCoords<T> (-1, +1, +1));
139 list.push_back(GridCoords<T> (0, -1, +1));
140 list.push_back(GridCoords<T> (0, 0, +1));
141 list.push_back(GridCoords<T> (0, +1, +1));
142 list.push_back(GridCoords<T> (+1, -1, +1));
143 list.push_back(GridCoords<T> (+1, 0, +1));
144 list.push_back(GridCoords<T> (+1, +1, +1));
145 }
146
147 template <typename E, typename T = short, typename Ty = BasicGrid_TAG>
148 class BaseGrid;
149
150 template <typename E, typename T, typename Ty>
151 class BaseGrid
152 {
153 public:
154 class NeighborIterator : public std::iterator<std::input_iterator_tag, GridCoords <T> >
155 {
156 private:
157 const BaseGrid & _parent;
158 GridCoords <T> _start, _curr;
159 typename std::list <GridCoords <T> >::const_iterator _itneigh;
160
161 CR_CONTEXT;
162
163 private:
164 inline NeighborIterator (const BaseGrid & parent, const GridCoords <T> & start)
165 : _start (start), _parent(parent)
166 {
167 CR_INIT();
168 process();
169 }
170 void process ()
171 {
172 CR_START();
173
174 for (_itneigh = _parent._neighbors.begin(); _itneigh != _parent._neighbors.end(); _itneigh++)
175 {
176 _curr.x = _start.x + _itneigh->x;
177 _curr.y = _start.y + _itneigh->y;
178 _curr.z = _start.z + _itneigh->z;
179
180 if (_parent.checkBounds(tx2basic(_parent.TYPE, _curr)) && _parent.check(_curr))
181 {
182 CR_RETURN_VOID;
183 }
184 }
185
186 CR_END();
187 }
188
189 public:
190 /*
191 NeighborIterator (const NeighborIterator & copy)
192 : _parent (copy._parent), _start (copy._start), _curr (copy._curr), _itneigh (copy._itneigh), CR_COPY(copy) {}
193 */
194
195 inline GridCoords <T> & operator*() { return _curr; }
196 inline GridCoords <T> * operator->() { return &_curr; }
197 inline operator bool () { return _parent._neighbors.end() != _itneigh; }
198 inline NeighborIterator & operator++()
199 {
200 process();
201 return *this;
202 }
203 inline NeighborIterator operator++(int)
204 {
205 NeighborIterator old = *this;
206
207 process();
208 return old;
209 }
210
211 friend class BaseGrid;
212 };
213
214 private:
215 typedef std::list <GridCoords <T> > NeighborsList;
216 T _w, _h, _d;
217 E ** _grid;
218 Ty TYPE;
219 float _interval;
220 NeighborsList _neighbors;
221
222 protected:
223 template <typename TyA>
224 inline GridCoords<T> fixNC (const GridCoords <T> & coords) const
225 {
226 TyA TYPE;
227 const GridCoords <T> & t = tx2basic (TYPE, coords);
228
229 return t;
230 }
231 inline GridCoords<T> fixNC (const GridCoords <T> & coords) const
232 {
233 return fixNC <Ty> (coords);
234 }
235
236 template <typename TyA>
237 inline GridCoords<T> fix (const GridCoords <T> & coords) const
238 {
239 TyA TYPE;
240 const GridCoords <T> & t = tx2basic (TYPE, coords);
241
242 if (!checkBounds(t))
243 throw GridEx<T>(t);
244
245 return t;
246 }
247 template <typename TyA>
248 inline bool checkCoords (const GridCoords <T> & coords) const
249 {
250 TyA TYPE;
251 return checkBounds(tx2basic(TYPE, coords));
252 }
253
254 inline GridCoords<T> fix (const GridCoords <T> & coords) const
255 {
256 return fix <Ty> (coords);
257 }
258 inline bool checkBounds (const GridCoords <T> & coords) const
259 {
260 return (coords.x >= 0 && coords.y >= 0 && coords.z >= 0 && coords.x < _w && coords.y < _h && coords.z < _d);
261 }
262
263 template <typename TyA>
264 inline const GridCoords <T> dim () const
265 {
266 TyA TYPE;
267 return txFbasic(TYPE, GridCoords <T> (_w, _h, _d));
268 }
269
270 inline const long index (const GridCoords <T> & gc) const
271 {
272 return gc.z * _w * _h + gc.y * _w + gc.x;
273 }
274
275 public:
276 inline BaseGrid(T w, T h, T d)
277 : _w(w), _h(h), _d(d), _grid(new E*[_w *_h * _d])
278 {
279 memset(_grid, 0, _w * _h * _d * sizeof (E *));
280 populate_neighbors(TYPE, _neighbors);
281
282 float sum = 0.0f;
283
284 for (typename NeighborsList::iterator it = _neighbors.begin(); it != _neighbors.end(); it++)
285 {
286 GridCoords <T> & gc = (*it);
287
288 sum += sqrt(static_cast<float> (gc.x * gc.x + gc.y * gc.y + gc.z * gc.z));
289 }
290 _interval = sum / _neighbors.size();
291 }
292
293 inline const T width () const { return dim <Ty> ().x; }
294 inline const T height () const { return dim <Ty> ().y; }
295 inline const T depth () const { return dim <Ty> ().z; }
296
297 inline bool valid (const GridCoords <T> & gc) const { return checkCoords <Ty> (gc); }
298
299 template <typename TyA>
300 inline const T width () const { return dim <TyA> ().x; }
301
302 template <typename TyA>
303 inline const T height () const { return dim <TyA> ().y; }
304
305 template <typename TyA>
306 inline const T depth () const { return dim <TyA> ().z; }
307
308 inline NeighborIterator neighbors (const GridCoords <T> & coords)
309 {
310 return NeighborIterator(*this, coords);
311 }
312 inline bool move (const GridCoords <T> & fromA, const GridCoords <T> & toA)
313 {
314 const long & ifrom = index(fix(fromA));
315 const GridCoords <T> & to = fixNC(toA);
316
317 if (!checkBounds(to))
318 {
319 delete _grid[ifrom];
320 _grid[ifrom] = NULL;
321 return false;
322 } else
323 {
324 const long & ito = index(to);
325
326 if (ifrom == ito)
327 return true;
328
329 if (_grid[ito] != NULL)
330 delete _grid[ito];
331
332 _grid[ito] = _grid[ifrom];
333 _grid[ifrom] = NULL;
334
335 return true;
336 }
337 }
338 inline void erase (const GridCoords <T> & at)
339 {
340 const long & i = index(fix(at));
341
342 delete _grid[i];
343 _grid[i] = NULL;
344 }
345 inline float interval () const { return _interval; }
346
347 inline void set (const GridCoords <T> & coordsA, const E & elem)
348 {
349 const long & i = index(fix(coordsA));
350
351 if (_grid[i] == NULL)
352 _grid[i] = new E (elem);
353 else
354 {
355 delete _grid[i];
356 _grid[i] = new E (elem);
357 }
358 }
359
360 template <typename TyA>
361 inline void set (const GridCoords <T> & coordsA)
362 {
363 const long & i = index( fix <TyA> (coordsA) );
364
365 if (_grid[i] == NULL)
366 _grid[i] = new E ();
367 else
368 {
369 delete _grid[i];
370 _grid[i] = new E();
371 }
372 }
373 inline E * getp (const GridCoords <T> & coordsA)
374 {
375 return _grid[index( fix(coordsA) )];
376 }
377 inline E & get (const GridCoords <T> & coordsA)
378 {
379 const long & i = index( fix(coordsA) );
380
381 if (_grid[i] == NULL)
382 _grid[i] = new E ();
383
384 return *_grid[i];
385 }
386 inline E & operator [] (const GridCoords <T> & coordsA)
387 {
388 return get(coordsA);
389 }
390 inline bool check (const GridCoords <T> & coordsA) const
391 {
392 const GridCoords <T> gc = fix(coordsA);
393 int idx = index(gc);
394 return _grid[idx] != NULL;
395 }
396 inline const GridCoords <T> & oob (const GridCoords <T> & from)
397 {
398 const GridCoords<T> & coords = fix (from);
399 GridCoords <T> result;
400
401 if (coords.x < 0)
402 result.x = -1;
403 else if (coords.x >= _w)
404 result.x = 1;
405 else
406 result.x = 0;
407
408 if (coords.y < 0)
409 result.y = -1;
410 else if (coords.y >= _h)
411 result.y = 1;
412 else
413 result.x = 0;
414
415 if (coords.z < 0)
416 result.z = -1;
417 else if (coords.z >= _d)
418 result.z = 1;
419 else
420 result.z = 0;
421
422 return result;
423 }
424
425 inline virtual ~BaseGrid()
426 {
427 for (T zi = 0; zi < _d; zi++)
428 for (T yi = 0; yi < _h; yi++)
429 for (T xi = 0; xi < _w; xi++)
430 delete _grid[index(GridCoords<T> (xi, yi, zi))];
431
432 delete [] _grid;
433 }
434 };
435
436 enum GridType
437 {
438 Normal,
439 Wrap
440 };
441 enum StitchEdge
442 {
443 StitchLeft = 1 << 0,
444 StitchTop = 1 << 1,
445 StitchRight = 1 << 2,
446 StitchBottom = 1 << 3
447 };
448
449 template< typename T >
450 class MemoryBlock
451 {
452 public:
453 T * data;
454 const size_t scanline, numlines, length, bytesize, linebytesize;
455
456 explicit MemoryBlock (const MemoryBlock & copy)
457 : data(new T[copy.length]),
458 scanline(copy.scanline), numlines(copy.numlines), length(copy.length), bytesize(copy.bytesize), linebytesize(copy.linebytesize)
459 {
460 // TODO: memcpy might be faster in the case of MPDG if using 32-bit "int" as opposed to 16-bit "short"
461 memcpy(data, copy.data, bytesize);
462 }
463
464 MemoryBlock (const size_t nScanLine, const size_t nNumScanLines)
465 : data (new T[nScanLine * nNumScanLines]), scanline(nScanLine), numlines(nNumScanLines),
466 length(nScanLine * nNumScanLines), bytesize(nScanLine * nNumScanLines * sizeof(T)), linebytesize(nScanLine * sizeof(T))
467 {}
468
469 inline T * offset (const unsigned int x, const unsigned int y)
470 {
471 assert (x >= 0 && x < scanline && y >= 0 && y < numlines);
472 return &data[y * scanline + x];
473 }
474 inline void memcopy (T * dest, const MemoryBlock< T > * pBlockSrc, const T * src, const unsigned short count)
475 {
476 assert(dest >= data && dest - data + count <= static_cast< long > (length));
477 assert(src >= pBlockSrc->data && src - pBlockSrc->data + count <= static_cast< long > (pBlockSrc->length));
478 memcpy(dest, src, count * sizeof(T));
479 }
480 inline void memcopy (T * dest, const T * src, const unsigned short count)
481 {
482 assert(dest >= data && dest - data + count < static_cast< long > (length));
483 memcpy(dest, src, count * sizeof(T));
484 }
485
486 inline std::istream & operator << (std::istream & ins)
487 {
488 ins.read(reinterpret_cast< char * > (data), bytesize);
489 return ins;
490 }
491 inline std::ostream & operator >> (std::ostream & outs) const
492 {
493 outs.write(reinterpret_cast< char * > (data), bytesize);
494 return outs;
495 }
496
497 ~MemoryBlock()
498 {
499 LOGD << "Delete " << data;
500 delete [] data;
501 }
502 };
503
504 template <typename T, typename ST = Numerics< T > >
505 class ScalarFieldView
506 {
507 private:
508#define KERNEL_OPERATION(k,op) KERNEL_OPERATION_T(T, k, op)
509#define CONST_KERNEL_OPERATION(k,op) KERNEL_OPERATION_T(const T, k, op)
510#define KERNEL_OPERATION_T(TYPE,k,op) \
511 if (_bContiguous) \
512 { \
513 int i, len = static_cast< int > (length); \
514 PARALLEL_FOR(i) \
515 for (i = 0; i < len; ++i) \
516 { \
517 TYPE & k = operator[] (i); \
518 op; \
519 } \
520 } \
521 else \
522 { \
523 int i,j; \
524 PARALLEL_FOR2(i,j) \
525 for (j = 0; j < height; ++j) \
526 { \
527 TYPE * pLine = line(j); \
528 for (i = 0; i < width; ++i) \
529 { \
530 TYPE & k = pLine[i]; \
531 op; \
532 } \
533 } \
534 }
535
536#define CHUNKED_KERNEL_OPERATION(k,l,op) \
537 if (_bContiguous) \
538 { \
539 T * k = _map; \
540 const unsigned int l = length; \
541 op; \
542 } else \
543 { \
544 int j; \
545 \
546 PARALLEL_FOR(j) \
547 for (j = 0; j < height; ++j) \
548 { \
549 T * k = &_map[j * _pBlock->scanline]; \
550 const unsigned int l = width; \
551 op; \
552 } \
553 }
554
555#define DUAL_CHUNKED_KERNEL_OPERATION(hm,a,b,l,op) DUAL_CHUNKED_KERNEL_OPERATION_T(ScalarFieldView< T, ST >, hm, a, b, l, op)
556#define DUAL_CHUNKED_KERNEL_OPERATION_T(J,hm,a,b,l,op) \
557 if (_bContiguous && hm._bContiguous) \
558 { \
559 T * a = _map; \
560 const J::Precision * b = hm._map; \
561 const unsigned int l = length; \
562 op; \
563 } else \
564 { \
565 int j; \
566 \
567 PARALLEL_FOR(j) \
568 for (j = 0; j < height; ++j) \
569 { \
570 T * a = line(j); \
571 const J::Precision * b = hm.line(j); \
572 const unsigned int l = width; \
573 op; \
574 } \
575 }
576
577#define DUAL_KERNEL_OPERATION(hm,k,l,op) DUAL_KERNEL_OPERATION_T(ScalarFieldView< T, ST >, hm, k, l, op)
578#define DUAL_KERNEL_OPERATION_T(J,hm,k,l,op) \
579 if (_bContiguous) \
580 { \
581 int i, len = static_cast< int > (length); \
582 PARALLEL_FOR(i) \
583 for (i = 0; i < len; ++i) \
584 { \
585 T & k = operator[] (i); \
586 const typename J::Precision & l = hm[i]; \
587 op; \
588 } \
589 } \
590 else \
591 { \
592 int i,j; \
593 PARALLEL_FOR2(i,j) \
594 for (j = 0; j < height; ++j) \
595 { \
596 T * pDestLine = line(j); \
597 const typename J::Precision * pSrcLine = hm.line(j); \
598 for (i = 0; i < width; ++i) \
599 { \
600 T & k = pDestLine[i]; \
601 const typename J::Precision & l = pSrcLine[i]; \
602 op; \
603 } \
604 } \
605 }
606
607 mutable MemoryBlock< T > * _pBlock;
608 T * _map;
609
610 const bool _bContiguous;
611 mutable unsigned short _nViewRefCount;
612
613 inline void assertDimensions () const
614 {
615 assert(length < (1 << 31) && "width and height must be less than 32768");
616 assert(right >= left && bottom >= top);
617 assert(height == bottom - top + 1 && width == right - left + 1);
618 assert(&_pBlock->data[top * _pBlock->scanline + left] == _map);
619 assert(&_map[(height - 1) * _pBlock->scanline + width] == &_pBlock->data[bottom * _pBlock->scanline + right + 1]);
620 assert(top >= 0 && left >= 0);
621 assert(_pBlock->numlines >= static_cast< size_t > (height + top));
622 assert(_pBlock->scanline >= static_cast< size_t > (width + left));
623 }
624
625 protected:
626 ScalarFieldView (MemoryBlock< T > * pBlock, const unsigned short nLeft, const unsigned short nTop, const unsigned short nWidth, const unsigned short nHeight, const GridType engt = Normal)
627 : _pBlock(pBlock), _map(pBlock->offset(nLeft, nTop)),
628 _bContiguous(pBlock->scanline == nWidth), _nViewRefCount(0), width(nWidth), height(nHeight), left(nLeft), top(nTop), right(nLeft + nWidth - 1),
629 bottom(nTop + nHeight - 1),
630 length(pBlock->length), gridtype(engt)
631 { assertDimensions(); }
632 ScalarFieldView (MemoryBlock< T > * pBlock, const BBox< unsigned short > & bbox, const GridType engt = Normal)
633 : _pBlock(pBlock), _map(pBlock->offset(bbox.left, bbox.top)), left(bbox.left), top(bbox.top), width(bbox.getWidth()), length(pBlock->length), height(bbox.getHeight()), right(bbox.right), bottom(bbox.bottom),
634 _bContiguous(pBlock->scanline == bbox.getWidth()),
635 gridtype(engt), _nViewRefCount(0)
636 { assertDimensions(); }
637
638 // Reserved for use by the ScalarField and subclasses, used to duplicate the memory block
639 ScalarFieldView (const ScalarFieldView & copy, MemoryBlock< T > * const pBlock)
640 : _pBlock(pBlock), _map(pBlock->offset(copy.left, copy.top)), left(copy.left), top(copy.top), width(copy.width), height(copy.height), length(copy.length), right(copy.right), bottom(copy.bottom),
641 _bContiguous(copy._bContiguous),
642 gridtype(copy.gridtype), _nViewRefCount(0)
643 {
644 assert(pBlock != copy._pBlock); // Use the copy constructor instead
645 assert(pBlock->scanline == copy._pBlock->scanline && pBlock->numlines == copy._pBlock->numlines);
646 assertDimensions();
647 }
648
649 // Copy constructor
650 ScalarFieldView (const ScalarFieldView & copy)
651 : _pBlock(copy._pBlock), _map(copy._map), left(copy.left), top(copy.top), width(copy.width), height(copy.height), length(copy.length), right(copy.right), bottom(copy.bottom),
652 _bContiguous(copy._bContiguous),
653 gridtype(copy.gridtype), _nViewRefCount(0)
654 {}
655
656 template< typename VT, typename ED >
657 VT * createTypedView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom, ED & ed)
658 {
659 assert (nLeft >= 0 && nRight < width && nTop >= 0 && nBottom < height);
660 assert (nRight < width && nBottom < height);
661 assert (nRight >= nLeft && nBottom >= nTop);
662
663 ++_nViewRefCount;
664 return new VT (
665 _pBlock,
666 BBox< unsigned short > (
667 nLeft + left,
668 nTop + top,
669 nRight + left,
670 nBottom + top
671 ),
672 ed
673 );
674 }
675
676 template< typename VT, typename ED >
677 const VT * createTypedView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom, const ED & ed) const
678 {
679 assert (nLeft >= 0 && nRight < width && nTop >= 0 && nBottom < height);
680 assert (nRight < width && nBottom < height);
681 assert (nRight >= nLeft && nBottom >= nTop);
682
683 ++_nViewRefCount;
684 return new VT (
685 _pBlock,
686 BBox< unsigned short > (
687 nLeft + left,
688 nTop + top,
689 nRight + left,
690 nBottom + top
691 ),
692 ed
693 );
694 }
695
696 template< typename VT >
697 VT * createTypedView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom)
698 {
699 assert (nLeft >= 0 && nRight < width && nTop >= 0 && nBottom < height);
700 assert (nRight < width && nBottom < height);
701 assert (nRight >= nLeft && nBottom >= nTop);
702
703 ++_nViewRefCount;
704 return new VT (
705 _pBlock,
706 BBox< unsigned short > (
707 nLeft + left,
708 nTop + top,
709 nRight + left,
710 nBottom + top
711 )
712 );
713 }
714
715 template< typename VT >
716 const VT * createTypedView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom) const
717 {
718 assert (nLeft >= 0 && nRight < width && nTop >= 0 && nBottom < height);
719 assert (nRight < width && nBottom < height);
720 assert (nRight >= nLeft && nBottom >= nTop);
721
722 ++_nViewRefCount;
723 return new VT (
724 _pBlock,
725 BBox< unsigned short > (
726 nLeft + left,
727 nTop + top,
728 nRight + left,
729 nBottom + top
730 )
731 );
732 }
733
734 inline MemoryBlock< T > * block() { return _pBlock; }
735 inline const MemoryBlock< T > * block() const { return _pBlock; }
736
737 inline const T * data() const { return _map; }
738 inline T * data() { return _map; }
739 inline const T * line(const unsigned short j) const
740 {
741 assert(j >= 0 && j < height);
742 return &_map[j * _pBlock->scanline];
743 }
744 inline T * line(const unsigned short j)
745 {
746 assert(j >= 0 && j < height);
747 return &_map[j * _pBlock->scanline];
748 }
749
750 static inline void determineOrder (const ScalarFieldView & sfv1, const ScalarFieldView & sfv2, const ScalarFieldView * & psfvSmaller, const ScalarFieldView * & psfvBigger)
751 {
752 if (sfv1.width < sfv2.width || sfv1.height < sfv2.height)
753 {
754 psfvBigger = &sfv2;
755 psfvSmaller = &sfv1;
756 } else
757 {
758 psfvBigger = &sfv1;
759 psfvSmaller = &sfv2;
760 }
761 }
762
763 friend class ScalarFieldView< long double >;
764 friend class ScalarFieldView< double >;
765 friend class ScalarFieldView< float >;
766 friend class ScalarFieldView< signed long >;
767 friend class ScalarFieldView< unsigned long >;
768 friend class ScalarFieldView< signed int >;
769 friend class ScalarFieldView< unsigned int >;
770 friend class ScalarFieldView< signed short >;
771 friend class ScalarFieldView< unsigned short >;
772 friend class ScalarFieldView< signed char >;
773 friend class ScalarFieldView< unsigned char >;
774
775 public:
776 typedef T Precision;
777 typedef ST ScalarTraits;
778
779 class Offset
780 {
781 public:
782 ScalarFieldView * const view;
783 const ScalarFieldView * const_view;
784 const signed int x, y;
785
786 Offset (ScalarFieldView * const pView, const signed int x, const signed int y)
787 : view(pView), const_view(pView), x(x), y(y) {}
788 Offset (const ScalarFieldView * pView, const signed int x, const signed int y)
789 : view(NULL), const_view(pView), x(x), y(y) {}
790
791 void operator >> (ScalarFieldView & other) const
792 { const_view->blit(other, x, y); }
793 void operator << (const ScalarFieldView & other)
794 { other.blit(*view, x, y); }
795 };
796 typedef ScalarFieldView < T, ST > View;
797
798 const unsigned short width, height, left, top, right, bottom;
799 const unsigned int length;
800 const GridType gridtype;
801
802 ScalarFieldView (ScalarFieldView && move)
803 : _pBlock(move._pBlock), _map(move._map), left(move.left), top(move.top), width(move.width), height(move.height), length(move.length), right(move.right), bottom(move.bottom),
804 _bContiguous(move._bContiguous),
805 gridtype(move.gridtype), _nViewRefCount(move._nViewRefCount)
806 {}
807
808 virtual ~ScalarFieldView ()
809 { assert(_nViewRefCount == 0); }
810
811 inline mars::BBox< unsigned short > bbox () const
812 { return mars::BBox< unsigned short > (left, top, left + width - 1, top + height - 1); }
813
814 View * createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom)
815 { return createTypedView< View > (nLeft, nTop, nRight, nBottom); }
816
817 const View * createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom) const
818 { return createTypedView< View > (nLeft, nTop, nRight, nBottom); }
819
820 inline static View createHosted(const unsigned short nWidth, const unsigned short nHeight, const GridType engt)
821 { return View (new MemoryBlock< T >(nWidth, nHeight), 0, 0, nWidth, nHeight, engt); }
822
823 virtual void releaseView (const View * pView) const
824 {
825 delete pView;
826 --_nViewRefCount;
827 }
828
829 const MemoryBlock< T > & getMemoryBlock () const { return *_pBlock; }
830
831 virtual mars::Range< T, ST > range () const
832 {
833 if (length > 0)
834 {
835 T max = ST::min();
836 T min = ST::max();
837
838 CONST_KERNEL_OPERATION(k,{
839 if (k > max)
840 max = k;
841 if (k < min)
842 min = k;
843 });
844
845 return mars::Range< T, ST > (min, max);
846 } else
847 {
848 return mars::Range< T, ST > ();
849 }
850 }
851
852 virtual void clear ()
853 { CHUNKED_KERNEL_OPERATION(k,l, memset(k, 0, l * sizeof(T))); }
854
855 virtual void fill (const T init)
856 { KERNEL_OPERATION(k, k = init); }
857
858 template< typename J, typename JST >
859 void add (const ScalarFieldView< J, JST > & hm)
860 {
861 assert(width == hm.width && height == hm.height);
862 typedef ScalarFieldView< J, JST > OtherView;
863 DUAL_KERNEL_OPERATION_T(OtherView, hm, k, l, k += static_cast< T > (l));
864 }
865 template< typename J, typename JST >
866 void sub (const ScalarFieldView< J, JST > & hm)
867 {
868 assert(width == hm.width && height == hm.height);
869 typedef ScalarFieldView< J, JST > OtherView;
870 DUAL_KERNEL_OPERATION_T(OtherView, hm, k, l, k -= static_cast< T > (l));
871 }
872
873 void add (const T & n)
874 { KERNEL_OPERATION(val, val += n); }
875 void sub (const T & n)
876 { KERNEL_OPERATION(val, val -= n); }
877
878 void mult (const T nFactor)
879 { KERNEL_OPERATION(val, val *= nFactor); }
880 void div (const T nFactor)
881 { KERNEL_OPERATION(val, val /= nFactor); }
882
883 template< typename J, typename JST >
884 void unionn (const ScalarFieldView< J, JST > & hm)
885 {
886 assert(width == hm.width && height == hm.height);
887 typedef ScalarFieldView< J, JST > OtherView;
888 DUAL_KERNEL_OPERATION_T(OtherView, hm, k, l, {
889 if (k < static_cast< T > (l))
890 k = static_cast< T > (l);
891 });
892 }
893 template< typename J, typename JST >
894 void intersection (const ScalarFieldView< J, JST > & hm)
895 {
896 assert(width == hm.width && height == hm.height);
897 typedef ScalarFieldView< J, JST > OtherView;
898 DUAL_KERNEL_OPERATION_T(OtherView, hm, k, l, {
899 if (k > static_cast< T > (l))
900 k = static_cast< T > (l);
901 });
902 }
903 inline View & operator += (const View & hm)
904 {
905 add(hm);
906 return *this;
907 }
908 inline View & operator += (const T & n)
909 {
910 add(n);
911 return *this;
912 }
913 template< typename J, typename JST >
914 inline View & operator += (const ScalarFieldView< J, JST > & hm)
915 {
916 add< J, JST >(hm);
917 return *this;
918 }
919 inline View & operator -= (const View & hm)
920 {
921 sub(hm);
922 return *this;
923 }
924 inline View & operator -= (const T & n)
925 {
926 sub(n);
927 return *this;
928 }
929 template< typename J, typename JST >
930 inline View & operator -= (const ScalarFieldView< J, JST > & hm)
931 {
932 sub< J, JST >(hm);
933 return *this;
934 }
935 template< typename J, typename JST >
936 inline View & operator |= (const ScalarFieldView< J, JST > & hm)
937 {
938 unionn< J, JST > (hm);
939 return *this;
940 }
941 template< typename J, typename JST >
942 inline View & operator &= (const ScalarFieldView< J, JST > & hm)
943 {
944 intersection< J, JST > (hm);
945 return *this;
946 }
947 inline View & operator *= (const T nFactor)
948 {
949 mult(nFactor);
950 return *this;
951 }
952 inline View & operator /= (const T nFactor)
953 {
954 div(nFactor);
955 return *this;
956 }
957
958 inline T & operator [] (const unsigned int index)
959 {
960 assert(index < length);
961 return _map[index];
962 }
963 inline const T & operator [] (const unsigned int index) const
964 {
965 assert(index < length);
966 return _map[index];
967 }
968
969 inline T & at (const unsigned short x, const unsigned short y)
970 {
971 assert(x >= 0 && x < width && y >= 0 && y < height);
972 return _map[y * _pBlock->scanline + x];
973 }
974 inline const T & at (const unsigned short x, const unsigned short y) const
975 {
976 assert(x >= 0 && x < width && y >= 0 && y < height);
977 return _map[y * _pBlock->scanline + x];
978 }
979
980 inline T & atw (const signed int x, const signed int y)
981 { return at(WRAP(x, width), WRAP(y, height)); }
982 inline const T & atw (const signed int x, const signed int y) const
983 { return at(WRAP(x, width), WRAP(y, height)); }
984
985 inline T & setn (const unsigned short x, const unsigned short y, const T & val)
986 { return at(x, y) = val; }
987 inline T & setw (const signed int x, const signed int y, const T & val)
988 { return atw(x, y) = val; }
989
990 inline T & set (const signed int x, const signed int y, const T & val)
991 {
992 switch (gridtype)
993 {
994 case Wrap: return setw(x, y, val);
995 case Normal:
996 default:
997 return setn(x, y, val);
998 }
999 }
1000
1001 inline void setw (const signed int x, const signed int y, const T * pnData, const unsigned short nCount)
1002 {
1003 const unsigned short
1004 x0 = WRAP(x, width),
1005 y0 = WRAP(y, height),
1006 xN = x0 + nCount;
1007
1008 if (xN > width)
1009 {
1010 setn(0, y0, pnData, xN - width);
1011 setn(x0, y0, pnData, (width - 1) - x0);
1012 }
1013 else
1014 setn(x0, y0, pnData, nCount);
1015 }
1016 inline void setn (const unsigned short x, const unsigned short y, const T * pnData, const unsigned short nCount)
1017 {
1018 assert(y >= 0 && y < height && x >= 0 && x <= width);
1019 assert(x + nCount <= width);
1020 _pBlock->memcopy(&line(y)[x], pnData, nCount);
1021 }
1022
1023 inline void set (const signed int x, const signed int y, const T * pnData, const unsigned short nCount)
1024 {
1025 switch (gridtype)
1026 {
1027 case Wrap: return setw(x, y, pnData, nCount);
1028 case Normal:
1029 default:
1030 return setn(x, y, pnData, nCount);
1031 }
1032 }
1033
1034 inline const T & getn (const unsigned short x, const unsigned short y) const
1035 { return at (x, y); }
1036 inline const T & getw (const signed int x, const signed int y) const
1037 { return atw (x, y); }
1038
1039 inline T & getn (const unsigned short x, const unsigned short y)
1040 { return at (x, y); }
1041 inline T & getw (const signed int x, const signed int y)
1042 { return atw (x, y); }
1043
1044 inline T & get(const signed int x, const signed int y)
1045 {
1046 switch (gridtype)
1047 {
1048 case Wrap: return getw(x, y);
1049 case Normal:
1050 default:
1051 return getn(x, y);
1052 }
1053 }
1054 inline const T & get(const signed int x, const signed int y) const
1055 {
1056 switch (gridtype)
1057 {
1058 case Wrap: return getw(x, y);
1059 case Normal:
1060 default:
1061 return getn(x, y);
1062 }
1063 }
1064
1065 inline const T & operator () (const signed int x, const signed int y) const
1066 { return get(x, y); }
1067
1068 inline T & operator () (const signed int x, const signed int y)
1069 { return get(x, y); }
1070
1071 inline void stitch(const int nStitch)
1072 { stitch(nStitch, this); }
1073 inline void stitch(const int nStitch, const View * pView)
1074 { stitchAll(nStitch, pView, pView, pView, pView); }
1075 void stitchAll(const int nStitch)
1076 { stitchAll(nStitch, this, this, this, this); }
1077 void stitchAll(const int nStitch, const View * pLeft, const View * pTop, const View * pRight, const View * pBottom)
1078 {
1079 assert(!(nStitch & StitchTop) || pTop != NULL);
1080 assert(!(nStitch & StitchBottom) || pBottom != NULL);
1081 assert(!(nStitch & StitchLeft) || pLeft != NULL);
1082 assert(!(nStitch & StitchRight) || pRight != NULL);
1083
1084 if (nStitch & StitchTop)
1085 {
1086 assert(pTop->width == width && pTop->height == height);
1087 _pBlock->memcopy(line(0), pTop->_pBlock, pTop->line(bottom), width);
1088 }
1089 if (nStitch & StitchBottom)
1090 {
1091 assert(pBottom->width == width && pBottom->height == height);
1092 _pBlock->memcopy(line(bottom), pBottom->_pBlock, pBottom->line(0), width);
1093 }
1094
1095 if (nStitch & StitchLeft)
1096 {
1097 assert(pLeft->width == width && pLeft->height == height);
1098 int j;
1099
1100#ifdef _OPENMP
1101#pragma omp parallel for private(j)
1102#endif
1103 for (j = 0; j < height; ++j)
1104 set(0, j, pLeft->get(right, j));
1105 }
1106 if (nStitch & StitchRight)
1107 {
1108 assert(pRight->width == width && pRight->height == height);
1109 int j;
1110
1111#ifdef _OPENMP
1112#pragma omp parallel for private(j)
1113#endif
1114 for (j = 0; j < height; ++j)
1115 set(right, j, pRight->get(0, j));
1116 }
1117 }
1118
1119 void copyEdges(const int nEdges, const View * pOther)
1120 {
1121 assert(width == pOther->width && height == pOther->height);
1122
1123 if (nEdges & StitchTop)
1124 _pBlock->memcopy(line(0), pOther->_pBlock, pOther->line(0), width);
1125 if (nEdges & StitchBottom)
1126 _pBlock->memcopy(line(bottom), pOther->_pBlock, pOther->line(bottom), width);
1127
1128 if (nEdges & StitchLeft)
1129 {
1130 int j;
1131
1132#ifdef _OPENMP
1133#pragma omp parallel for private(j)
1134#endif
1135 for (j = 0; j < height; ++j)
1136 set(0, j, pOther->get(0, j));
1137 }
1138 if (nEdges & StitchRight)
1139 {
1140 int j;
1141
1142#ifdef _OPENMP
1143#pragma omp parallel for private(j)
1144#endif
1145 for (j = 0; j < height; ++j)
1146 set(right, j, pOther->get(right, j));
1147 }
1148 }
1149
1150 inline bool isInside (const signed int x, const signed int y) const
1151 {
1152 return x >= left && y >= top && x <= right && y <= bottom;
1153 }
1154
1155 // Samples elevation in a rotated Von Neumann neighborhood
1156 inline T meanVN (const signed int x, const signed int y) const
1157 {
1158 typedef typename ST::UP R;
1159
1160 return static_cast <T> (
1161 (
1162 static_cast <R> (get(x - 1, y - 1)) +
1163 static_cast <R> (get(x + 1, y - 1)) +
1164 static_cast <R> (get(x - 1, y + 1)) +
1165 static_cast <R> (get(x + 1, y + 1))
1166 )
1167 / 4);
1168 }
1169
1170 inline T mean5 (const signed int x, const signed int y) const
1171 {
1172 typedef typename ST::UP R;
1173 return static_cast <T> (
1174 (
1175 static_cast <R> (get(x - 1, y - 1)) +
1176 static_cast <R> (get(x, y - 1)) +
1177 static_cast <R> (get(x + 1, y - 1)) +
1178 static_cast <R> (get(x - 1, y)) +
1179 static_cast <R> (get(x, y)) +
1180 static_cast <R> (get(x + 1, y)) +
1181 static_cast <R> (get(x - 1, y + 1)) +
1182 static_cast <R> (get(x, y + 1)) +
1183 static_cast <R> (get(x + 1, y + 1)) +
1184
1185 static_cast <R> (get(x - 2, y - 2)) +
1186 static_cast <R> (get(x - 2, y - 1)) +
1187 static_cast <R> (get(x - 2, y)) +
1188 static_cast <R> (get(x - 2, y + 1)) +
1189 static_cast <R> (get(x - 2, y + 2)) +
1190
1191 static_cast <R> (get(x + 2, y - 2)) +
1192 static_cast <R> (get(x + 2, y - 1)) +
1193 static_cast <R> (get(x + 2, y)) +
1194 static_cast <R> (get(x + 2, y + 1)) +
1195 static_cast <R> (get(x + 2, y + 2)) +
1196
1197 static_cast <R> (get(x - 1, y - 2)) +
1198 static_cast <R> (get(x, y - 2)) +
1199 static_cast <R> (get(x + 1, y - 2)) +
1200
1201 static_cast <R> (get(x - 1, y + 2)) +
1202 static_cast <R> (get(x, y + 2)) +
1203 static_cast <R> (get(x + 1, y + 2))
1204 ) / 25
1205 );
1206 }
1207
1208 inline T mean (const signed int x, const signed int y) const
1209 {
1210 typedef typename ST::UP R;
1211
1212 R a = 0;
1213
1214 a += get(x-1, y-1); a += get(x+0, y-1); a += get(x+1, y-1);
1215 a += get(x-1, y-0); a += get(x+0, y+0); a += get(x+1, y-0);
1216 a += get(x-1, y+1); a += get(x+0, y+1); a += get(x+1, y+1);
1217
1218 return static_cast <T> (a / 9);
1219 }
1220
1221 static void blitw (
1222 const ScalarFieldView & sfvSrc, ScalarFieldView & sfvDest,
1223 const signed int x = 0, const signed int y = 0
1224 )
1225 {
1226 const ScalarFieldView
1227 * psfvBig,
1228 * psfvSmall;
1229
1230 determineOrder(sfvSrc, sfvDest, psfvSmall, psfvBig);
1231
1232 assert(psfvBig->width >= psfvSmall->width && psfvBig->height >= psfvSmall->height);
1233 assert(
1234 sfvSrc.left >= 0 && sfvDest.left >= 0 && sfvSrc.top >= 0 && sfvDest.top >= 0 &&
1235 sfvSrc.right >= 0 && sfvDest.right >= 0 && sfvSrc.bottom >= 0 && sfvDest.bottom >= 0
1236 );
1237
1238 const unsigned short
1239 w = psfvSmall->width,
1240 h = psfvSmall->height;
1241
1242 const unsigned short
1243 x0 = WRAP(x, psfvBig->width),
1244 y0 = WRAP(y, psfvBig->height),
1245
1246 xN = x0 + w - 1,
1247 yN = y0 + h - 1;
1248
1249 if (xN >= psfvBig->width || yN >= psfvBig->height)
1250 {
1251 const unsigned short
1252 xH = std::min(xN, static_cast< unsigned short > (psfvBig->width - 1)),
1253 yH = std::min(yN, static_cast< unsigned short > (psfvBig->height - 1)),
1254
1255 w0 = xH - x0 + 1,
1256 h0 = yH - y0 + 1,
1257
1258 wN = psfvSmall->width - w0,
1259 hN = psfvSmall->height - h0;
1260
1261 assert(psfvBig->gridtype == mars::Wrap);
1262 if (psfvBig == &sfvDest)
1263 {
1264 if (w0 > 0 && h0 > 0) blitn (sfvSrc, sfvDest, w0, h0, 0, 0, x0, y0);
1265 if (wN > 0 && h0 > 0) blitn (sfvSrc, sfvDest, wN, h0, w0, 0, 0, y0);
1266 if (w0 > 0 && hN > 0) blitn (sfvSrc, sfvDest, w0, hN, 0, h0, x0, 0);
1267 if (wN > 0 && hN > 0) blitn (sfvSrc, sfvDest, wN, hN, w0, h0, 0, 0);
1268 } else
1269 {
1270 if (w0 > 0 && h0 > 0) blitn (sfvSrc, sfvDest, w0, h0, x0, y0, 0, 0);
1271 if (wN > 0 && h0 > 0) blitn (sfvSrc, sfvDest, wN, h0, 0, y0, w0, 0);
1272 if (w0 > 0 && hN > 0) blitn (sfvSrc, sfvDest, w0, hN, x0, 0, 0, h0);
1273 if (wN > 0 && hN > 0) blitn (sfvSrc, sfvDest, wN, hN, 0, 0, w0, h0);
1274 }
1275 } else
1276 {
1277 if (psfvBig == &sfvDest)
1278 blitn(sfvSrc, sfvDest, psfvSmall->width, psfvSmall->height, 0, 0, x0, y0);
1279 else
1280 blitn(sfvSrc, sfvDest, psfvSmall->width, psfvSmall->height, x0, y0, 0, 0);
1281 }
1282 }
1283
1284 static void blitn (
1285 const ScalarFieldView & sfvSrc, ScalarFieldView & sfvDest,
1286 const unsigned short w, const unsigned short h,
1287 const unsigned short xs = 0, const unsigned short ys = 0,
1288 const unsigned short xd = 0, const unsigned short yd = 0
1289 )
1290 {
1291 assert(w > 0 && h > 0);
1292 assert(w <= sfvSrc.width && h <= sfvSrc.height);
1293 assert(w <= sfvDest.width && h <= sfvDest.height);
1294
1295 assert(xs >= 0 && ys >= 0);
1296 assert(xd >= 0 && yd >= 0);
1297 assert(
1298 (xs + w <= sfvSrc.width && ys + h <= sfvSrc.height) &&
1299 (xd + w <= sfvDest.width && yd + h <= sfvDest.height)
1300 );
1301
1302 T * pnDestOffset = sfvDest._pBlock->offset(xd + sfvDest.left, yd + sfvDest.top);
1303 const T * pnSrcOffset = sfvSrc._pBlock->offset(xs + sfvSrc.left, ys + sfvSrc.top);
1304 const size_t nCount = w;
1305
1306 int j;
1307#ifdef _OPENMP
1308#pragma omp parallel for private(j)
1309#endif
1310 for (j = 0; j < h; ++j)
1311 sfvDest.block() ->memcopy(
1312 &pnDestOffset[j * sfvDest._pBlock->scanline],
1313 sfvSrc._pBlock,
1314 &pnSrcOffset[j * sfvSrc._pBlock->scanline],
1315 nCount
1316 );
1317 }
1318
1319 inline void blit (ScalarFieldView & sfvDest, const signed int x = 0, const signed int y = 0) const
1320 {
1321 if (gridtype == Wrap || sfvDest.gridtype == Wrap) // TODO: Confirm bounds too
1322 blitw(*this, sfvDest, x, y);
1323 else
1324 {
1325 const ScalarFieldView
1326 * psfvBig,
1327 * psfvSmall;
1328
1329 determineOrder(*this, sfvDest, psfvSmall, psfvBig);
1330 if (psfvBig == &sfvDest)
1331 blitn(*this, sfvDest, psfvSmall->width, psfvSmall->height, 0, 0, x, y);
1332 else
1333 blitn(*this, sfvDest, psfvSmall->width, psfvSmall->height, x, y, 0, 0);
1334 }
1335 }
1336
1337 void operator >> (std::vector< std::vector< T > > & stdfield) const
1338 {
1339 stdfield.resize(height);
1340
1341 size_t j;
1342 PARALLEL_FOR(j);
1343 for (j = 0; j < height; ++j)
1344 {
1345 stdfield[j].resize(width);
1346 const T * pLine = line(j);
1347 for (size_t i = 0; i < width; ++i)
1348 stdfield[j][i] = pLine[i];
1349 }
1350 }
1351
1352 inline void operator >> (ScalarFieldView & dest) const
1353 {
1354 assert(width == dest.width && height == dest.height);
1355 blit(dest);
1356 }
1357 inline void operator << (const ScalarFieldView & other)
1358 {
1359 assert(width == other.width && height == other.height);
1360 other.blit(*this);
1361 }
1362 inline void operator >> (Offset & dest) const
1363 {
1364 blit(*dest.view, dest.x, dest.y);
1365 }
1366 inline void operator << (const Offset & other)
1367 {
1368 other.const_view->blit(*this, other.x, other.y);
1369 }
1370
1371 template< typename J, typename JST >
1372 void operator << (const ScalarFieldView< J, JST > & other)
1373 {
1374 assert(width == other.width && height == other.height);
1375 typedef ScalarFieldView< J, JST > OtherView;
1376 DUAL_KERNEL_OPERATION_T(OtherView, other,nDest,nSrc,{
1377 nDest = static_cast< T > (nSrc);
1378 });
1379 }
1380
1381 template< typename J, typename JST >
1382 inline void operator >> (ScalarFieldView< J, JST > & dest) const
1383 { dest << *this; }
1384 };
1385
1386 template< typename T, typename ST = Numerics< T > >
1387 class ScalarField : public ScalarFieldView< T, ST >
1388 {
1389 private:
1390 typedef ScalarFieldView< T, ST > Supa;
1391
1392 public:
1393 ScalarField (const unsigned short width, const unsigned short height, const ScalarField & copy, const unsigned short x, const unsigned short y, const GridType engt = Normal)
1394 : Supa(Supa::createHosted(width, height, engt))
1395 { blitFrom (copy, x, y); }
1396
1397 ScalarField (const unsigned short width, const unsigned short height, const GridType engt = Normal)
1398 : Supa(Supa::createHosted(width, height, engt))
1399 { memset(this->data(), 0, this->block() -> bytesize); }
1400 ScalarField (const Supa & copy)
1401 : Supa(copy, new MemoryBlock< T >(copy.getMemoryBlock()))
1402 {}
1403 ScalarField (const ScalarField & copy)
1404 : Supa(copy, new MemoryBlock< T >(copy.getMemoryBlock()))
1405 {}
1406
1407 inline std::istream & operator << (std::istream & ins)
1408 { return *this->block() << ins; }
1409 inline std::ostream & operator >> (std::ostream & outs) const
1410 { return *this->block() >> outs; }
1411 inline void operator >> (Supa & dest) const
1412 { Supa::operator >> (dest); }
1413 inline void operator << (const Supa & other)
1414 { Supa::operator << (other); }
1415 inline void operator >> (typename Supa::Offset dest) const
1416 { Supa::operator >> (dest); }
1417 inline void operator << (const typename Supa::Offset other)
1418 { Supa::operator << (other); }
1419 template< typename J, typename JST >
1420 inline void operator << (const ScalarFieldView< J, JST > & other)
1421 { Supa::operator << (other); }
1422 template< typename J, typename JST >
1423 inline void operator >> (ScalarFieldView< J, JST > & dest) const
1424 { Supa::operator >> (dest); }
1425
1426 virtual ~ScalarField() { delete this->block(); }
1427 };
1428
1429 template <typename T, typename ST = Numerics< T > >
1430 class MidpointDisplacementGrid : public ScalarFieldView< T, ST >
1431 {
1432 private:
1433 typedef MidpointDisplacementGrid< T, ST > Self;
1434 typedef ScalarFieldView< T, ST > Supa;
1435
1436 const enum StorageMode
1437 {
1438 Concrete,
1439 Virtual
1440 } _ensmStorageMode;
1441
1442 unsigned short _c;
1443 unsigned int _stride;
1444
1445 static inline unsigned short createMinDim (const ScalarFieldView< T, ST > & sfv)
1446 { return std::min(sfv.width, sfv.height); }
1447 static inline unsigned int createStride (const ScalarFieldView< T, ST > & sfv)
1448 {
1449 const unsigned short nMinDim = createMinDim(sfv);
1450 return nMinDim - nMinDim % 2;
1451 }
1452 static inline unsigned short createMinDim (const unsigned short nWidth, const unsigned short nHeight)
1453 { return std::min(nWidth, nHeight); }
1454 static inline unsigned int createStride (const unsigned short nWidth, const unsigned short nHeight)
1455 {
1456 const unsigned short nMinDim = createMinDim(nWidth, nHeight);
1457 return nMinDim - nMinDim % 2;
1458 }
1459
1460 inline void assertions () const
1461 {
1462 assert(this->width == dimensionFor(this->width, Supa::gridtype) && this->height == dimensionFor(this->height, Supa::gridtype));
1463 assert(_stride >= 2);
1464 }
1465
1466 protected:
1467 MidpointDisplacementGrid (const short nIterCount, const unsigned short nWidth, const unsigned short nHeight, const GridType engt = mars::Normal)
1468 : Supa (Supa::createHosted(nWidth, nHeight, engt)), _c (0),
1469 _stride(createStride(nWidth, nHeight)),
1470 itercount(nIterCount), lastX(nWidth - 1), lastY(nHeight - 1), mindim(createMinDim(nWidth, nHeight)), _ensmStorageMode(Virtual)
1471 { assertions(); }
1472 MidpointDisplacementGrid (const short nIterCount, ScalarFieldView< T, ST > & sfv)
1473 : Supa (sfv), _c (0), itercount(nIterCount), lastX(sfv.width - 1), lastY(sfv.height - 1), _stride(createStride(sfv)), mindim(createMinDim(sfv)), _ensmStorageMode(Concrete)
1474 { assertions(); }
1475
1476 inline static short calcIterCount (const unsigned short nMinDim)
1477 {
1478 return static_cast <unsigned short> (
1479 ceil(
1480 log(static_cast <float> (nMinDim)) / log(2.0f)
1481 )
1482 );
1483 }
1484 inline static unsigned short calcDimLast (const unsigned short nIterCount)
1485 {
1486 return 1 << nIterCount;
1487 }
1488
1489 public:
1490 template <typename S, unsigned dBp = 1>
1491 class AbstractFunctor
1492 {
1493 private:
1494 template <unsigned N>
1495 inline const unsigned int pow(const unsigned short c, quantity<N>) const
1496 { return c * pow(c, quantity<N-1>()); }
1497
1498 inline const unsigned int pow (const unsigned short c, quantity<0>) const
1499 { return 1; }
1500
1501 protected:
1502 const unsigned short _nFlux;
1503 const float _fCoarseness;
1504
1505 float _nRandAmount;
1506
1507 template <typename REAL>
1508 inline REAL applyFlux (const REAL & rMean) const
1509 {
1510 return rMean + static_cast <REAL> (mars::RANDf(_nRandAmount)) - static_cast <REAL> (_nRandAmount) / static_cast <REAL> (2);
1511 }
1512 inline float computeRandomAmount (const unsigned short c) const
1513 {
1514 return static_cast <float> (_nFlux) / static_cast< float > (pow<dBp>(c + 1, quantity<dBp>())) * _fCoarseness;
1515 }
1516
1517 public:
1518 typedef S StateType;
1519
1520 inline AbstractFunctor (const unsigned short nFlux, const float fCoarseness)
1521 : _nFlux (nFlux), _fCoarseness (fCoarseness) {}
1522
1523 inline const float & getRandomAmount () const { return _nRandAmount; }
1524
1525 inline void iteration (const unsigned short nStride, const unsigned short c)
1526 {
1527 //_nRandAmount = _nDim / (static_cast <int> (static_cast <float> (c) * (1.0f / _fCoarseness)) + 1);
1528 //_nRandAmount = static_cast <float> (_nStride) * _fCoarseness;
1529 // XTODO: Struggling with finding an appropriate noise generating algorithm here.
1530 // The problem is it's either too noisy at the finest level resulting in what look like innumerable spires everywhere
1531 // or the crater doesn't look noisy enough and looks like a perfectly uniform ring
1532 _nRandAmount = computeRandomAmount(c);
1533 }
1534
1535 inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c, const T & d) const;
1536 inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c) const;
1537 };
1538
1539 class PinkFn : public AbstractFunctor<int>
1540 {
1541 public:
1542 inline PinkFn (const unsigned short nFlux, const float fCoarseness)
1543 : AbstractFunctor<int>(nFlux, fCoarseness) {}
1544
1545 inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c, const T & d) const
1546 {
1547 return static_cast <T> (AbstractFunctor<int>::applyFlux(static_cast <float> (a + b + c + d) / 4.0f));
1548 }
1549 inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c) const
1550 {
1551 return static_cast <T> (AbstractFunctor<int>::applyFlux(static_cast <float> (a + b + c) / 3.0f));
1552 }
1553 };
1554
1555 class BrownianFn : public AbstractFunctor<int, 2>
1556 {
1557 public:
1558 inline BrownianFn (const unsigned short nFlux, const float fCoarseness)
1559 : AbstractFunctor<int, 2>(nFlux, fCoarseness) {}
1560
1561 inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c, const T & d) const
1562 {
1563 return static_cast <T> (AbstractFunctor<int, 2>::applyFlux(static_cast <float> (a + b + c + d) / 4.0f));
1564 }
1565 inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c) const
1566 {
1567 return static_cast <T> (AbstractFunctor<int, 2>::applyFlux(static_cast <float> (a + b + c) / 3.0f));
1568 }
1569 };
1570
1571 class EdgeSeedIterator : public std::iterator <std::input_iterator_tag, T>
1572 {
1573 private:
1574 MidpointDisplacementGrid< T > & _grid;
1575 unsigned short _j, _i;
1576 unsigned short _halfstride;
1577 bool _done;
1578
1579 CR_CONTEXT;
1580
1581 void process ()
1582 {
1583 // NOTE: Redundant code in MidpointDisplacementGrid::iterate(...);
1584 CR_START();
1585
1586 _halfstride = _grid.stride() / 2;
1587 assert(_halfstride > 0);
1588
1589 for (_j = 0; _j <= _grid.lastY; _j += _grid.stride())
1590 for (_i = 0; _i <= _grid.lastX; _i += _grid.stride())
1591 {
1592 CR_RETURN_VOID;
1593 }
1594
1595 for (_j = 0; _j <= _grid.lastY; _j += _halfstride)
1596 for (_i = (_halfstride + _j) % _grid.stride(); _i <= _grid.lastX; _i += _grid.stride())
1597 {
1598 CR_RETURN_VOID;
1599 }
1600
1601 _done = true;
1602 CR_END();
1603 }
1604
1605 public:
1606 inline EdgeSeedIterator (MidpointDisplacementGrid & grid)
1607 : _grid(grid), _done(false)
1608 {
1609 CR_INIT();
1610 process();
1611 }
1612
1613 inline EdgeSeedIterator & operator++ ()
1614 {
1615 process();
1616 return *this;
1617 }
1618 inline EdgeSeedIterator operator++ (int)
1619 {
1620 EdgeSeedIterator old = *this;
1621
1622 process();
1623 return old;
1624 }
1625
1626 inline T & assign (const T & a) { return _grid(_i, _j) = a; }
1627 inline T & ref () { return _grid(_i, _j); }
1628 inline const T & ref () const { return _grid(_i, _j); }
1629
1630 inline T & operator = (const T & a) { return assign(a); }
1631
1632 inline const unsigned short getX() const { return _i; }
1633 inline const unsigned short getY() const { return _j; }
1634
1635 inline operator bool () const { return !_done; }
1636 };
1637
1638 public:
1639 const unsigned short itercount;
1640 const unsigned short mindim, lastX, lastY;
1641
1642 virtual ~MidpointDisplacementGrid()
1643 {
1644 switch (_ensmStorageMode)
1645 {
1646 case Concrete: delete this->block(); break;
1647 }
1648 }
1649
1650 const unsigned short index () const { return _c; }
1651 const unsigned short stride () const { return _stride; }
1652
1653 inline static Self * createInstance (const unsigned short nMinWidth, const unsigned short nMinHeight, const GridType engt)
1654 {
1655 assert(nMinHeight > 0 && nMinWidth > 0);
1656 return new Self (
1657 calcIterCount(std::min(nMinHeight, nMinWidth)),
1658 dimensionFor(nMinWidth, engt),
1659 dimensionFor(nMinHeight, engt),
1660 engt
1661 );
1662 }
1663 inline static Self * createInstance (Supa & sfv)
1664 {
1665 return new Self (
1666 calcIterCount(std::min(sfv.width, sfv.height) - (sfv.gridtype == Normal ? 1 : 0)),
1667 sfv
1668 );
1669 }
1670
1671 inline static unsigned short dimensionFor (const unsigned short n, const GridType engt)
1672 {
1673 switch (engt)
1674 {
1675 case Wrap:
1676 return calcDimLast (calcIterCount (n));
1677 case Normal:
1678 default:
1679 return calcDimLast (calcIterCount (n - n % 2)) + 1;
1680 }
1681 }
1682
1683 inline void run (const float fCoarseNess, const int nStitch = 0)
1684 { run (BrownianFn (mindim, fCoarseNess), nStitch); }
1685
1686 template <typename F>
1687 void run (F functor, const int nStitch = 0)
1688 {
1689 iterations (itercount, functor, nStitch);
1690 }
1691 template <typename F>
1692 void iterations (const unsigned short nCount, F & functor, const int nStitch = 0)
1693 {
1694 unsigned short nLastExcl = _c + nCount;
1695
1696 while (_c <= nLastExcl && _stride >= 2)
1697 {
1698 iterate (functor, nStitch);
1699 }
1700 }
1701
1702 inline EdgeSeedIterator edgeSeedIterator () { return EdgeSeedIterator (*this); }
1703
1704 void seedFrom (const ScalarFieldView< T, ST > & other, const unsigned short xOffs, const unsigned short yOffs, const unsigned short nLOD)
1705 {
1706 int i, j, x, y;
1707
1708 _c = itercount - nLOD;
1709 _stride = (1 << nLOD);
1710
1711 assert(this->height / _stride + yOffs <= other.height);
1712 assert(this->width / _stride + xOffs <= other.width);
1713
1714#ifdef _OPENMP
1715#pragma omp parallel for private(i, j, x, y)
1716#endif
1717 for (j = 0; j < this->height; j += _stride)
1718 {
1719 y = yOffs + (j >> nLOD);
1720 for (i = 0, x = xOffs; i < this->width; i += _stride, ++x)
1721 Self::operator () (i, j) = other(x, y) << nLOD;
1722 }
1723 }
1724
1725 inline void reset ()
1726 { step(-_c); }
1727
1728 inline void step (const signed int nSteps = 1)
1729 {
1730 _c += nSteps;
1731 if (nSteps < 0)
1732 _stride *= (1 << -nSteps);
1733 else
1734 _stride /= (1 << nSteps);
1735 }
1736
1737 template < typename Functor > void iterate (Functor & func, const int nStitch = 0)
1738 {
1739 const int
1740 nHalfStride = _stride / 2;
1741
1742 int i, j;
1743 typename Functor::StateType state;
1744
1745 assert(nHalfStride > 0);
1746 func.iteration (_stride, _c);
1747
1748 switch (this->gridtype)
1749 {
1750 case Wrap:
1751 // Square phase
1752 // NOTE: Redundant code in EdgeSeedIterator
1753#ifdef _OPENMP
1754#pragma omp parallel for private(i, j, state)
1755#endif
1756 for (j = 0; j <= lastY; j += _stride)
1757 for (i = 0; i <= lastX; i += _stride)
1758 Self::operator()(i + nHalfStride, j + nHalfStride) =
1759 func(state,
1760 i + nHalfStride, j + nHalfStride,
1761 Self::operator() (i, j),
1762 Self::operator() ((i + _stride) % this->width, j),
1763 Self::operator() (i, (j + _stride) % this->height),
1764 Self::operator() ((i + _stride) % this->width, (j + _stride) % this->height)
1765 );
1766
1767 // Diamond phase
1768 // NOTE: Redundant code in EdgeSeedIterator and other iterator
1769#ifdef _OPENMP
1770#pragma omp parallel for private(i, j, state)
1771#endif
1772 for (j = 0; j <= lastY; j += nHalfStride)
1773 for (i = (nHalfStride + j) % _stride; i <= lastX; i += _stride)
1774 Self::operator() (i, j) =
1775 func(state,
1776 i, j,
1777 Self::operator()(i, (j - nHalfStride + this->height) % this->height),
1778 Self::operator()(i, (j + nHalfStride) % this->height),
1779 Self::operator()((i - nHalfStride + this->width) % this->width, j),
1780 Self::operator()((i + nHalfStride) % this->width, j)
1781 );
1782 break;
1783 case Normal:
1784 default:
1785 {
1786 const int
1787 nCtxLastX = lastX - nHalfStride,
1788 nCtxLastY = lastY - nHalfStride;
1789
1790 // Square-phase
1791 // NOTE: Redundant code in EdgeSeedIterator
1792#ifdef _OPENMP
1793#pragma omp parallel for private(i, j, state)
1794#endif
1795 for (j = 0; j < nCtxLastY; j += _stride)
1796 for (i = 0; i < nCtxLastX; i += _stride)
1797 Self::operator()(i + nHalfStride, j + nHalfStride) =
1798 func(state,
1799 i + nHalfStride, j + nHalfStride,
1800 Self::operator() (i, j),
1801 Self::operator() (i + _stride, j),
1802 Self::operator() (i, j + _stride),
1803 Self::operator() (i + _stride, j + _stride)
1804 );
1805
1806#ifdef _OPENMP
1807#pragma omp parallel sections private(i, j, state)
1808#endif
1809 {
1810 // Next 4 blocks explicitly perform the square-step (second-step) of the diamond-
1811 // square algorithm along all four edges of the grid so that we don't have to do
1812 // any bounds checking for the rest of the iteration
1813
1814 // Top edge
1815#ifdef _OPENMP
1816#pragma omp section
1817#endif
1818 if ((nStitch ^ StitchTop) & StitchTop)
1819 {
1820 const int jj = 0;
1821 for (i = nHalfStride; i < lastX; i += _stride)
1822 Self::operator() (i, jj) =
1823 func(state,
1824 i, jj,
1825 Self::operator()(i - nHalfStride, jj),
1826 Self::operator()(i + nHalfStride, jj),
1827 Self::operator()(i, jj + nHalfStride)
1828 );
1829 }
1830
1831 // Bottom edge
1832#ifdef _OPENMP
1833#pragma omp section
1834#endif
1835 if ((nStitch ^ StitchBottom) & StitchBottom)
1836 {
1837 const int j = lastY;
1838 for (i = nHalfStride; i < lastX; i += _stride)
1839 Self::operator() (i, j) =
1840 func(state,
1841 i, j,
1842 Self::operator()(i - nHalfStride, j),
1843 Self::operator()(i + nHalfStride, j),
1844 Self::operator()(i, j - nHalfStride)
1845 );
1846 }
1847
1848 // Left edge
1849#ifdef _OPENMP
1850#pragma omp section
1851#endif
1852 if ((nStitch ^ StitchLeft) & StitchLeft)
1853 {
1854 const int i = 0;
1855 for (j = nHalfStride; j < lastY; j += _stride)
1856 Self::operator() (i, j) =
1857 func(state,
1858 i, j,
1859 Self::operator()(i, j - nHalfStride),
1860 Self::operator()(i, j + nHalfStride),
1861 Self::operator()(i + nHalfStride, j)
1862 );
1863 }
1864
1865 // Right edge
1866#ifdef _OPENMP
1867#pragma omp section
1868#endif
1869 if ((nStitch ^ StitchRight) & StitchRight)
1870 {
1871 const int i = lastX;
1872 for (j = nHalfStride; j < lastY; j += _stride)
1873 Self::operator() (i, j) =
1874 func(state,
1875 i, j,
1876 Self::operator()(i, j - nHalfStride),
1877 Self::operator()(i, j + nHalfStride),
1878 Self::operator()(i - nHalfStride, j)
1879 );
1880 }
1881 }
1882
1883 // Perform the second phase (square phase) of the diamond-square algorithm on the rest of
1884 // the grid excluding its edges.
1885 // NOTE: Redundant code in EdgeSeedIterator
1886#ifdef _OPENMP
1887#pragma omp parallel for private(i, j, state)
1888#endif
1889 for (j = nHalfStride; j < lastY; j += nHalfStride)
1890 for (i = nHalfStride + j % _stride; i < lastX; i += _stride)
1891 Self::operator() (i, j) =
1892 func(state,
1893 i, j,
1894 Self::operator()(i, j - nHalfStride),
1895 Self::operator()(i, j + nHalfStride),
1896 Self::operator()(i - nHalfStride, j),
1897 Self::operator()(i + nHalfStride, j)
1898 );
1899 }
1900 break;
1901 }
1902
1903 ++_c;
1904 _stride /= 2;
1905 }
1906
1907 inline std::istream & operator << (std::istream & ins)
1908 { return *this->block() << ins; }
1909 inline std::ostream & operator >> (std::ostream & outs) const
1910 { return *this->block() >> outs; }
1911
1912 };
1913}
1914
1915#endif
Note: See TracBrowser for help on using the repository browser.