#ifndef __MARSGRIDS_H__
#define __MARSGRIDS_H__

#include <exception>
#include <iterator>
#include <list>
#include <limits>
#include <istream>
#include <ostream>

#include <math.h>
#include <memory.h>

#ifdef _OPENMP
#include<omp.h>
#endif

#include "coroutines.h"
#include "mars_util.h"
#include "mars_calc.h"
#include "mars_log.h"

namespace mars
{
	template <typename T>
	class GridCoords
	{
	public:
		typedef T Precision;

		T x, y, z;

		inline GridCoords ()
			: x(0), y(0), z(0) {}
		inline GridCoords (T x, T y, T z)
			: x(x), y(y), z(z) {}
	};

	template <typename T>
	inline bool operator < (const GridCoords <T> & left, const GridCoords <T> & right)
	{
		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));
	}

	template <typename T>
	class GridEx : public std::exception
	{
	public:
		const GridCoords <T> coords;

		inline GridEx (const GridCoords <T> & coords)
			: coords(coords) {}
	};

	struct HexGrid_TAG {};
	struct BasicGrid_TAG {};

	template <typename T>
	inline GridCoords <T> tx2basic (HexGrid_TAG, const GridCoords <T> & coordsHexa)
	{
		return GridCoords <T>
			(
			coordsHexa.x / 2,
			coordsHexa.y / 2,
			coordsHexa.z
			);
	}

	template <typename T>
	inline GridCoords <T> tx2basic (BasicGrid_TAG, const GridCoords <T> & coords)
	{
		return coords;
	}

	template <typename T>
	inline GridCoords <T> txFbasic (HexGrid_TAG, const GridCoords <T> & coords)
	{
		return GridCoords <T>
			(
			coords.x * 2 + (coords.y + coords.z) % 2,
			coords.y * 2 + (coords.z % 2),
			coords.z
			);
	}

	template <typename T>
	inline GridCoords <T> txFbasic (BasicGrid_TAG, const GridCoords <T> & coords)
	{
		return coords;
	}

	template <typename T>
	void populate_neighbors (HexGrid_TAG, std::list <GridCoords <T> > & list)
	{
		list.push_back(GridCoords<T> (-2, -1, 0));
		list.push_back(GridCoords<T> (-2, +1, 0));
		list.push_back(GridCoords<T> (+2, -1, 0));
		list.push_back(GridCoords<T> (+2, +1, 0));
		list.push_back(GridCoords<T> (0, -2, 0));
		list.push_back(GridCoords<T> (0, +2, 0));
		list.push_back(GridCoords<T> (0, -2, -1));
		list.push_back(GridCoords<T> (0, +2, -1));
		list.push_back(GridCoords<T> (-2, 0, -1));
		list.push_back(GridCoords<T> (+2, 0, -1));
		list.push_back(GridCoords<T> (0, -2, +1));
		list.push_back(GridCoords<T> (0, +2, +1));
		list.push_back(GridCoords<T> (-2, 0, +1));
		list.push_back(GridCoords<T> (+2, 0, +1));
	}

	template <typename T>
	void populate_neighbors (BasicGrid_TAG, std::list <GridCoords <T> > & list)
	{
		list.push_back(GridCoords<T> (-1, -1, -1));
		list.push_back(GridCoords<T> (-1, 0, -1));
		list.push_back(GridCoords<T> (-1, +1, -1));
		list.push_back(GridCoords<T> (0, -1, -1));
		list.push_back(GridCoords<T> (0, 0, -1));
		list.push_back(GridCoords<T> (0, +1, -1));
		list.push_back(GridCoords<T> (+1, -1, -1));
		list.push_back(GridCoords<T> (+1, 0, -1));
		list.push_back(GridCoords<T> (+1, +1, -1));


		list.push_back(GridCoords<T> (-1, -1, 0));
		list.push_back(GridCoords<T> (-1, 0, 0));
		list.push_back(GridCoords<T> (-1, +1, 0));
		list.push_back(GridCoords<T> (0, -1, 0));

		list.push_back(GridCoords<T> (0, +1, 0));
		list.push_back(GridCoords<T> (+1, -1, 0));
		list.push_back(GridCoords<T> (+1, 0, 0));
		list.push_back(GridCoords<T> (+1, +1, 0));


		list.push_back(GridCoords<T> (-1, -1, +1));
		list.push_back(GridCoords<T> (-1, 0, +1));
		list.push_back(GridCoords<T> (-1, +1, +1));
		list.push_back(GridCoords<T> (0, -1, +1));
		list.push_back(GridCoords<T> (0, 0, +1));
		list.push_back(GridCoords<T> (0, +1, +1));
		list.push_back(GridCoords<T> (+1, -1, +1));
		list.push_back(GridCoords<T> (+1, 0, +1));
		list.push_back(GridCoords<T> (+1, +1, +1));
	}

	template <typename E, typename T = short, typename Ty = BasicGrid_TAG>
	class BaseGrid;

	template <typename E, typename T, typename Ty>
	class BaseGrid
	{
	public:
		class NeighborIterator : public std::iterator<std::input_iterator_tag, GridCoords <T> >
		{
		private:
			const BaseGrid & _parent;
			GridCoords <T> _start, _curr;
			typename std::list <GridCoords <T> >::const_iterator _itneigh;

			CR_CONTEXT;

		private:
			inline NeighborIterator (const BaseGrid & parent, const GridCoords <T> & start)
				: _start (start), _parent(parent)
			{
				CR_INIT();
				process();
			}
			void process ()
			{
				CR_START();

				for (_itneigh = _parent._neighbors.begin(); _itneigh != _parent._neighbors.end(); _itneigh++)
				{
					_curr.x = _start.x + _itneigh->x;
					_curr.y = _start.y + _itneigh->y;
					_curr.z = _start.z + _itneigh->z;

					if (_parent.checkBounds(tx2basic(_parent.TYPE, _curr)) && _parent.check(_curr))
					{
						CR_RETURN_VOID;
					}
				}

				CR_END();
			}

		public:
			/*
			NeighborIterator (const NeighborIterator & copy)
			: _parent (copy._parent), _start (copy._start), _curr (copy._curr), _itneigh (copy._itneigh), CR_COPY(copy) {}
			*/

			inline GridCoords <T> & operator*() { return _curr; }
			inline GridCoords <T> * operator->() { return &_curr; }
			inline operator bool () { return _parent._neighbors.end() != _itneigh; }
			inline NeighborIterator & operator++()
			{
				process();
				return *this;
			}
			inline NeighborIterator operator++(int)
			{
				NeighborIterator old = *this;

				process();
				return old;
			}

			friend class BaseGrid;
		};

	private:
		typedef std::list <GridCoords <T> > NeighborsList;
		T _w, _h, _d;
		E ** _grid;
		Ty TYPE;
		float _interval;
		NeighborsList _neighbors;

	protected:
		template <typename TyA>
		inline GridCoords<T> fixNC (const GridCoords <T> & coords) const
		{
			TyA TYPE;
			const GridCoords <T> & t = tx2basic (TYPE, coords);

			return t;
		}
		inline GridCoords<T> fixNC (const GridCoords <T> & coords) const
		{
			return fixNC <Ty> (coords);
		}

		template <typename TyA>
		inline GridCoords<T> fix (const GridCoords <T> & coords) const
		{
			TyA TYPE;
			const GridCoords <T> & t = tx2basic (TYPE, coords);

			if (!checkBounds(t))
				throw GridEx<T>(t);

			return t;
		}
		template <typename TyA>
		inline bool checkCoords (const GridCoords <T> & coords) const
		{
			TyA TYPE;
			return checkBounds(tx2basic(TYPE, coords));
		}

		inline GridCoords<T> fix (const GridCoords <T> & coords) const
		{
			return fix <Ty> (coords);
		}
		inline bool checkBounds (const GridCoords <T> & coords) const
		{
			return (coords.x >= 0 && coords.y >= 0 && coords.z >= 0 && coords.x < _w && coords.y < _h && coords.z < _d);
		}

		template <typename TyA>
		inline const GridCoords <T> dim () const
		{
			TyA TYPE;
			return txFbasic(TYPE, GridCoords <T> (_w, _h, _d));
		}

		inline const long index (const GridCoords <T> & gc) const
		{
			return gc.z * _w * _h + gc.y * _w + gc.x;
		}

	public:
		inline BaseGrid(T w, T h, T d)
			: _w(w), _h(h), _d(d), _grid(new E*[_w *_h * _d])
		{
			memset(_grid, 0, _w * _h * _d * sizeof (E *));
			populate_neighbors(TYPE, _neighbors);

			float sum = 0.0f;

			for (typename NeighborsList::iterator it = _neighbors.begin(); it != _neighbors.end(); it++)
			{
				GridCoords <T> & gc = (*it);

				sum += sqrt(static_cast<float> (gc.x * gc.x + gc.y * gc.y + gc.z * gc.z));
			}
			_interval = sum / _neighbors.size();
		}

		inline const T width () const { return dim <Ty> ().x; }
		inline const T height () const { return dim <Ty> ().y; }
		inline const T depth () const { return dim <Ty> ().z; }

		inline bool valid (const GridCoords <T> & gc) const { return checkCoords <Ty> (gc); }

		template <typename TyA>
		inline const T width () const { return dim <TyA> ().x; }

		template <typename TyA>
		inline const T height () const { return dim <TyA> ().y; }

		template <typename TyA>
		inline const T depth () const { return dim <TyA> ().z; }

		inline NeighborIterator neighbors (const GridCoords <T> & coords)
		{
			return NeighborIterator(*this, coords);
		}
		inline bool move (const GridCoords <T> & fromA, const GridCoords <T> & toA)
		{
			const long & ifrom = index(fix(fromA));
			const GridCoords <T> & to = fixNC(toA);

			if (!checkBounds(to))
			{
				delete _grid[ifrom];
				_grid[ifrom] = NULL;
				return false;
			} else
			{
				const long & ito = index(to);

				if (ifrom == ito)
					return true;

				if (_grid[ito] != NULL)
					delete _grid[ito];

				_grid[ito] = _grid[ifrom];
				_grid[ifrom] = NULL;

				return true;
			}
		}
		inline void erase (const GridCoords <T> & at)
		{
			const long & i = index(fix(at));

			delete _grid[i];
			_grid[i] = NULL;
		}
		inline float interval () const { return _interval; }

		inline void set (const GridCoords <T> & coordsA, const E & elem)
		{
			const long & i = index(fix(coordsA));

			if (_grid[i] == NULL)
				_grid[i] = new E (elem);
			else
			{
				delete _grid[i];
				_grid[i] = new E (elem);
			}
		}

		template <typename TyA>
		inline void set (const GridCoords <T> & coordsA)
		{
			const long & i = index( fix <TyA> (coordsA) );

			if (_grid[i] == NULL)
				_grid[i] = new E ();
			else
			{
				delete _grid[i];
				_grid[i] = new E();
			}
		}
		inline E * getp (const GridCoords <T> & coordsA)
		{
			return _grid[index( fix(coordsA) )];
		}
		inline E & get (const GridCoords <T> & coordsA)
		{
			const long & i = index( fix(coordsA) );

			if (_grid[i] == NULL)
				_grid[i] = new E ();

			return *_grid[i];
		}
		inline E & operator [] (const GridCoords <T> & coordsA)
		{
			return get(coordsA);
		}
		inline bool check (const GridCoords <T> & coordsA) const
		{
			const GridCoords <T> gc = fix(coordsA);
			int idx = index(gc);
			return _grid[idx] != NULL;
		}
		inline const GridCoords <T> & oob (const GridCoords <T> & from)
		{
			const GridCoords<T> & coords = fix (from);
			GridCoords <T> result;

			if (coords.x < 0)
				result.x  = -1;
			else if (coords.x >= _w)
				result.x = 1;
			else
				result.x = 0;

			if (coords.y < 0)
				result.y  = -1;
			else if (coords.y >= _h)
				result.y = 1;
			else
				result.x = 0;

			if (coords.z < 0)
				result.z  = -1;
			else if (coords.z >= _d)
				result.z = 1;
			else
				result.z = 0;

			return result;
		}

		inline virtual ~BaseGrid()
		{
			for (T zi = 0; zi < _d; zi++)
				for (T yi = 0; yi < _h; yi++)
					for (T xi = 0; xi < _w; xi++)
						delete _grid[index(GridCoords<T> (xi, yi, zi))];

			delete [] _grid;
		}
	};

	enum GridType
	{
		Normal,
		Wrap
	};
	enum StitchEdge
	{
		StitchLeft = 1 << 0,
		StitchTop = 1 << 1,
		StitchRight = 1 << 2,
		StitchBottom = 1 << 3
	};

	template< typename T >
	class MemoryBlock
	{
	public:
		T * data;
		const size_t scanline, numlines, length, bytesize, linebytesize;

		explicit MemoryBlock (const MemoryBlock & copy)
			: data(new T[copy.length]),
			scanline(copy.scanline), numlines(copy.numlines), length(copy.length), bytesize(copy.bytesize), linebytesize(copy.linebytesize)
		{
			// TODO: memcpy might be faster in the case of MPDG if using 32-bit "int" as opposed to 16-bit "short"
			memcpy(data, copy.data, bytesize);
		}

		MemoryBlock (const size_t nScanLine, const size_t nNumScanLines)
			: data (new T[nScanLine * nNumScanLines]), scanline(nScanLine), numlines(nNumScanLines),
			length(nScanLine * nNumScanLines), bytesize(nScanLine * nNumScanLines * sizeof(T)), linebytesize(nScanLine * sizeof(T))
		{}

		inline T * offset (const unsigned int x, const unsigned int y)
		{
			assert (x >= 0 && x < scanline && y >= 0 && y < numlines);
			return &data[y * scanline + x];
		}
		inline void memcopy (T * dest, const MemoryBlock< T > * pBlockSrc, const T * src, const unsigned short count)
		{
			assert(dest >= data && dest - data + count <= static_cast< long > (length));
			assert(src >= pBlockSrc->data && src - pBlockSrc->data + count <= static_cast< long > (pBlockSrc->length));
			memcpy(dest, src, count * sizeof(T));
		}
		inline void memcopy (T * dest, const T * src, const unsigned short count)
		{
			assert(dest >= data && dest - data + count < static_cast< long > (length));
			memcpy(dest, src, count * sizeof(T));
		}

		inline std::istream & operator << (std::istream & ins)
		{
			ins.read(reinterpret_cast< char * > (data), bytesize);
			return ins;
		}
		inline std::ostream & operator >> (std::ostream & outs) const
		{
			outs.write(reinterpret_cast< char * > (data), bytesize);
			return outs;
		}

		~MemoryBlock()
		{
			LOGD << "Delete " << data;
			delete [] data;
		}
	};

	template <typename T, typename ST = Numerics< T > >
	class ScalarFieldView
	{
	private:
#define KERNEL_OPERATION(k,op) KERNEL_OPERATION_T(T, k, op)
#define CONST_KERNEL_OPERATION(k,op) KERNEL_OPERATION_T(const T, k, op)
#define KERNEL_OPERATION_T(TYPE,k,op)								\
	if (_bContiguous)												\
	{																\
		int i, len = static_cast< int > (length);					\
		PARALLEL_FOR(i)												\
		for (i = 0; i < len; ++i)									\
		{															\
			TYPE & k = operator[] (i);								\
			op;														\
		}															\
	}																\
	else															\
	{																\
		int i,j;													\
		PARALLEL_FOR2(i,j)											\
		for (j = 0; j < height; ++j)								\
		{															\
			TYPE * pLine = line(j);									\
			for (i = 0; i < width; ++i)								\
			{														\
				TYPE & k = pLine[i];								\
				op;													\
			}														\
		}															\
	}

#define CHUNKED_KERNEL_OPERATION(k,l,op)							\
	if (_bContiguous)												\
	{																\
		T * k = _map;												\
		const unsigned int l = length;								\
		op;															\
	} else															\
	{																\
		int j;														\
																	\
		PARALLEL_FOR(j)												\
		for (j = 0; j < height; ++j)								\
		{															\
			T * k = &_map[j * _pBlock->scanline];					\
			const unsigned int l = width;							\
			op;														\
		}															\
	}

#define DUAL_CHUNKED_KERNEL_OPERATION(hm,a,b,l,op) DUAL_CHUNKED_KERNEL_OPERATION_T(ScalarFieldView< T, ST >, hm, a, b, l, op)
#define DUAL_CHUNKED_KERNEL_OPERATION_T(J,hm,a,b,l,op)				\
	if (_bContiguous && hm._bContiguous)							\
	{																\
		T * a = _map;												\
		const J::Precision * b = hm._map;							\
		const unsigned int l = length;								\
		op;															\
	} else															\
	{																\
		int j;														\
																	\
		PARALLEL_FOR(j)												\
		for (j = 0; j < height; ++j)								\
		{															\
			T * a = line(j);										\
			const J::Precision * b = hm.line(j);					\
			const unsigned int l = width;							\
			op;														\
		}															\
	}

#define DUAL_KERNEL_OPERATION(hm,k,l,op) DUAL_KERNEL_OPERATION_T(ScalarFieldView< T, ST >, hm, k, l, op)
#define DUAL_KERNEL_OPERATION_T(J,hm,k,l,op)						\
	if (_bContiguous)												\
	{																\
		int i, len = static_cast< int > (length);					\
		PARALLEL_FOR(i)												\
		for (i = 0; i < len; ++i)									\
		{															\
			T & k = operator[] (i);									\
			const typename J::Precision & l = hm[i];							\
			op;														\
		}															\
	}																\
	else															\
	{																\
		int i,j;													\
		PARALLEL_FOR2(i,j)											\
		for (j = 0; j < height; ++j)								\
		{															\
			T * pDestLine = line(j);								\
			const typename J::Precision * pSrcLine = hm.line(j);				\
			for (i = 0; i < width; ++i)								\
			{														\
				T & k = pDestLine[i];								\
				const typename J::Precision & l = pSrcLine[i];				\
				op;													\
			}														\
		}															\
	}

		mutable MemoryBlock< T > * _pBlock;
		T * _map;

		const bool _bContiguous;
		mutable unsigned short _nViewRefCount;

		inline void assertDimensions () const
		{
			assert(length < (1 << 31) && "width and height must be less than 32768");
			assert(right >= left && bottom >= top);
			assert(height == bottom - top + 1 && width == right - left + 1);
			assert(&_pBlock->data[top * _pBlock->scanline + left] == _map);
			assert(&_map[(height - 1) * _pBlock->scanline + width] == &_pBlock->data[bottom * _pBlock->scanline + right + 1]);
			assert(top >= 0 && left >= 0);
			assert(_pBlock->numlines >= static_cast< size_t > (height + top));
			assert(_pBlock->scanline >= static_cast< size_t > (width + left));
		}

	protected:
		ScalarFieldView (MemoryBlock< T > * pBlock, const unsigned short nLeft, const unsigned short nTop, const unsigned short nWidth, const unsigned short nHeight, const GridType engt = Normal)
			: _pBlock(pBlock), _map(pBlock->offset(nLeft, nTop)), 
			_bContiguous(pBlock->scanline == nWidth), _nViewRefCount(0), width(nWidth), height(nHeight), left(nLeft), top(nTop), right(nLeft + nWidth - 1),
			bottom(nTop + nHeight - 1), 
			length(pBlock->length), gridtype(engt)
			{ assertDimensions(); }
		ScalarFieldView (MemoryBlock< T > * pBlock, const BBox< unsigned short > & bbox, const GridType engt = Normal)
			: _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),
			_bContiguous(pBlock->scanline == bbox.getWidth()),
			gridtype(engt), _nViewRefCount(0)
			{ assertDimensions(); }

		// Reserved for use by the ScalarField and subclasses, used to duplicate the memory block
		ScalarFieldView (const ScalarFieldView & copy, MemoryBlock< T > * const pBlock)
			: _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),
			_bContiguous(copy._bContiguous),
			gridtype(copy.gridtype), _nViewRefCount(0)
		{
			assert(pBlock != copy._pBlock);	// Use the copy constructor instead
			assert(pBlock->scanline == copy._pBlock->scanline && pBlock->numlines == copy._pBlock->numlines);
			assertDimensions();
		}

		// Copy constructor
		ScalarFieldView (const ScalarFieldView & copy)
			: _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),
			_bContiguous(copy._bContiguous),
			gridtype(copy.gridtype), _nViewRefCount(0)
		{}

		template< typename VT, typename ED >
		VT * createTypedView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom, ED & ed)
		{
			assert (nLeft >= 0 && nRight < width && nTop >= 0 && nBottom < height);
			assert (nRight < width && nBottom < height);
			assert (nRight >= nLeft && nBottom >= nTop);

			++_nViewRefCount;
			return new VT (
				_pBlock,
				BBox< unsigned short > (
					nLeft + left,
					nTop + top,
					nRight + left,
					nBottom + top
				),
				ed
			);
		}

		template< typename VT, typename ED >
		const VT * createTypedView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom, const ED & ed) const
		{
			assert (nLeft >= 0 && nRight < width && nTop >= 0 && nBottom < height);
			assert (nRight < width && nBottom < height);
			assert (nRight >= nLeft && nBottom >= nTop);

			++_nViewRefCount;
			return new VT (
				_pBlock,
				BBox< unsigned short > (
					nLeft + left,
					nTop + top,
					nRight + left,
					nBottom + top
				),
				ed
			);
		}

		template< typename VT >
		VT * createTypedView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom)
		{
			assert (nLeft >= 0 && nRight < width && nTop >= 0 && nBottom < height);
			assert (nRight < width && nBottom < height);
			assert (nRight >= nLeft && nBottom >= nTop);

			++_nViewRefCount;
			return new VT (
				_pBlock,
				BBox< unsigned short > (
					nLeft + left,
					nTop + top,
					nRight + left,
					nBottom + top
				)
			);
		}

		template< typename VT >
		const VT * createTypedView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom) const
		{
			assert (nLeft >= 0 && nRight < width && nTop >= 0 && nBottom < height);
			assert (nRight < width && nBottom < height);
			assert (nRight >= nLeft && nBottom >= nTop);

			++_nViewRefCount;
			return new VT (
				_pBlock,
				BBox< unsigned short > (
					nLeft + left,
					nTop + top,
					nRight + left,
					nBottom + top
				)
			);
		}

		inline MemoryBlock< T > * block() { return _pBlock; }
		inline const MemoryBlock< T > * block() const { return _pBlock; }

		inline const T * data() const { return _map; }
		inline T * data() { return _map; }
		inline const T * line(const unsigned short j) const
		{
			assert(j >= 0 && j < height);
			return &_map[j * _pBlock->scanline];
		}
		inline T * line(const unsigned short j)
		{
			assert(j >= 0 && j < height);
			return &_map[j * _pBlock->scanline];
		}

		static inline void determineOrder (const ScalarFieldView & sfv1, const ScalarFieldView & sfv2, const ScalarFieldView * & psfvSmaller, const ScalarFieldView * & psfvBigger)
		{
			if (sfv1.width < sfv2.width || sfv1.height < sfv2.height)
			{
				psfvBigger = &sfv2;
				psfvSmaller = &sfv1;
			} else
			{
				psfvBigger = &sfv1;
				psfvSmaller = &sfv2;
			}
		}

		friend class ScalarFieldView< long double >;
		friend class ScalarFieldView< double >;
		friend class ScalarFieldView< float >;
		friend class ScalarFieldView< signed long >;
		friend class ScalarFieldView< unsigned long >;
		friend class ScalarFieldView< signed int >;
		friend class ScalarFieldView< unsigned int >;
		friend class ScalarFieldView< signed short >;
		friend class ScalarFieldView< unsigned short >;
		friend class ScalarFieldView< signed char >;
		friend class ScalarFieldView< unsigned char >;

	public:
		typedef T Precision;
		typedef ST ScalarTraits;

		class Offset
		{
		public:
			ScalarFieldView * const view;
			const ScalarFieldView * const_view;
			const signed int x, y;

			Offset (ScalarFieldView * const pView, const signed int x, const signed int y)
				: view(pView), const_view(pView), x(x), y(y) {}
			Offset (const ScalarFieldView * pView, const signed int x, const signed int y)
				: view(NULL), const_view(pView), x(x), y(y) {}

			void operator >> (ScalarFieldView & other) const
				{ const_view->blit(other, x, y); }
			void operator << (const ScalarFieldView & other)
				{ other.blit(*view, x, y); }
		};
		typedef ScalarFieldView < T, ST > View;

		const unsigned short width, height, left, top, right, bottom;
		const unsigned int length;
		const GridType gridtype;
        
		ScalarFieldView (ScalarFieldView && move)
			: _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),
			_bContiguous(move._bContiguous),
			gridtype(move.gridtype), _nViewRefCount(move._nViewRefCount)
		{}

		virtual ~ScalarFieldView ()
			{ assert(_nViewRefCount == 0); }

		inline mars::BBox< unsigned short > bbox () const
			{ return mars::BBox< unsigned short > (left, top, left + width - 1, top + height - 1); }

		View * createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom)
			{ return createTypedView< View > (nLeft, nTop, nRight, nBottom); }

		const View * createView (const unsigned short nLeft, const unsigned short nTop, const unsigned short nRight, const unsigned short nBottom) const
			{ return createTypedView< View > (nLeft, nTop, nRight, nBottom); }

        inline static View createHosted(const unsigned short nWidth, const unsigned short nHeight, const GridType engt)
            { return View (new MemoryBlock< T >(nWidth, nHeight), 0, 0, nWidth, nHeight, engt); }

		virtual void releaseView (const View * pView) const
		{
			delete pView;
			--_nViewRefCount;
		}

		const MemoryBlock< T > & getMemoryBlock () const { return *_pBlock; }

		virtual mars::Range< T, ST > range () const
		{
			if (length > 0)
			{
				T max = ST::min();
				T min = ST::max();

				CONST_KERNEL_OPERATION(k,{
					if (k > max)
						max = k;
					if (k < min)
						min = k;
				});

				return mars::Range< T, ST > (min, max);
			} else
			{
				return mars::Range< T, ST > ();
			}
		}

		virtual void clear ()
			{ CHUNKED_KERNEL_OPERATION(k,l, memset(k, 0, l * sizeof(T))); }

		virtual void fill (const T init)
			{ KERNEL_OPERATION(k, k = init); }

		template< typename J, typename JST >
		void add (const ScalarFieldView< J, JST > & hm)
		{
			assert(width == hm.width && height == hm.height);
			typedef ScalarFieldView< J, JST > OtherView;
			DUAL_KERNEL_OPERATION_T(OtherView, hm, k, l, k += static_cast< T > (l));
		}
		template< typename J, typename JST >
		void sub (const ScalarFieldView< J, JST > & hm)
		{
			assert(width == hm.width && height == hm.height);
			typedef ScalarFieldView< J, JST > OtherView;
			DUAL_KERNEL_OPERATION_T(OtherView, hm, k, l, k -= static_cast< T > (l));
		}

		void add (const T & n)
			{ KERNEL_OPERATION(val, val += n); }
		void sub (const T & n)
			{ KERNEL_OPERATION(val, val -= n); }

		void mult (const T nFactor)
			{ KERNEL_OPERATION(val, val *= nFactor); }
		void div (const T nFactor)
			{ KERNEL_OPERATION(val, val /= nFactor); }

		template< typename J, typename JST >
		void unionn (const ScalarFieldView< J, JST > & hm)
		{
			assert(width == hm.width && height == hm.height);
			typedef ScalarFieldView< J, JST > OtherView;
			DUAL_KERNEL_OPERATION_T(OtherView, hm, k, l, {
				if (k < static_cast< T > (l))
					k = static_cast< T > (l);
			});
		}
		template< typename J, typename JST >
		void intersection (const ScalarFieldView< J, JST > & hm)
		{
			assert(width == hm.width && height == hm.height);
			typedef ScalarFieldView< J, JST > OtherView;
			DUAL_KERNEL_OPERATION_T(OtherView, hm, k, l, {
				if (k > static_cast< T > (l))
					k = static_cast< T > (l);
			});
		}
		inline View & operator += (const View & hm)
		{
			add(hm);
			return *this;
		}
		inline View & operator += (const T & n)
		{
			add(n);
			return *this;
		}
		template< typename J, typename JST >
		inline View & operator += (const ScalarFieldView< J, JST > & hm)
		{
			add< J, JST >(hm);
			return *this;
		}
		inline View & operator -= (const View & hm)
		{
			sub(hm);
			return *this;
		}
		inline View & operator -= (const T & n)
		{
			sub(n);
			return *this;
		}
		template< typename J, typename JST >
		inline View & operator -= (const ScalarFieldView< J, JST > & hm)
		{
			sub< J, JST >(hm);
			return *this;
		}
		template< typename J, typename JST >
		inline View & operator |= (const ScalarFieldView< J, JST > & hm)
		{
			unionn< J, JST > (hm);
			return *this;
		}
		template< typename J, typename JST >
		inline View & operator &= (const ScalarFieldView< J, JST > & hm)
		{
			intersection< J, JST > (hm);
			return *this;
		}
		inline View & operator *= (const T nFactor)
		{
			mult(nFactor);
			return *this;
		}
		inline View & operator /= (const T nFactor)
		{
			div(nFactor);
			return *this;
		}

		inline T & operator [] (const unsigned int index)
		{
			assert(index < length);
			return _map[index];
		}
		inline const T & operator [] (const unsigned int index) const
		{
			assert(index < length);
			return _map[index];
		}

		inline T & at (const unsigned short x, const unsigned short y)
		{
			assert(x >= 0 && x < width && y >= 0  && y < height);
			return _map[y * _pBlock->scanline + x];
		}
		inline const T & at (const unsigned short x, const unsigned short y) const
		{
			assert(x >= 0 && x < width && y >= 0  && y < height);
			return _map[y * _pBlock->scanline + x];
		}

		inline T & atw (const signed int x, const signed int y)
			{ return at(WRAP(x, width), WRAP(y, height)); }
		inline const T & atw (const signed int x, const signed int y) const
			{ return at(WRAP(x, width), WRAP(y, height)); }

		inline T & setn (const unsigned short x, const unsigned short y, const T & val)
			{ return at(x, y) = val; }
		inline T & setw (const signed int x, const signed int y, const T & val)
			{ return atw(x, y) = val; }

		inline T & set (const signed int x, const signed int y, const T & val)
		{
			switch (gridtype)
			{
			case Wrap: return setw(x, y, val);
			case Normal:
			default:
				return setn(x, y, val);
			}
		}

		inline void setw (const signed int x, const signed int y, const T * pnData, const unsigned short nCount)
		{
			const unsigned short
				x0 = WRAP(x, width),
				y0 = WRAP(y, height),
				xN = x0 + nCount;

			if (xN > width)
			{
				setn(0, y0, pnData, xN - width);
				setn(x0, y0, pnData, (width - 1) - x0);
			}
			else
				setn(x0, y0, pnData, nCount);
		}
		inline void setn (const unsigned short x, const unsigned short y, const T * pnData, const unsigned short nCount)
		{
			assert(y >= 0 && y < height && x >= 0 && x <= width);
			assert(x + nCount <= width);
			_pBlock->memcopy(&line(y)[x], pnData, nCount);
		}

		inline void set (const signed int x, const signed int y, const T * pnData, const unsigned short nCount)
		{
			switch (gridtype)
			{
			case Wrap: return setw(x, y, pnData, nCount);
			case Normal:
			default:
				return setn(x, y, pnData, nCount);
			}
		}

		inline const T & getn (const unsigned short x, const unsigned short y) const
			{ return at (x, y); }
		inline const T & getw (const signed int x, const signed int y) const
			{ return atw (x, y); }

		inline T & getn (const unsigned short x, const unsigned short y)
			{ return at (x, y); }
		inline T & getw (const signed int x, const signed int y)
			{ return atw (x, y); }

		inline T & get(const signed int x, const signed int y)
		{
			switch (gridtype)
			{
			case Wrap: return getw(x, y);
			case Normal:
			default:
				return getn(x, y);
			}
		}
		inline const T & get(const signed int x, const signed int y) const
		{
			switch (gridtype)
			{
			case Wrap: return getw(x, y);
			case Normal:
			default:
				return getn(x, y);
			}
		}

		inline const T & operator () (const signed int x, const signed int y) const
			{ return get(x, y); }

		inline T & operator () (const signed int x, const signed int y)
			{ return get(x, y); }

		inline void stitch(const int nStitch)
			{ stitch(nStitch, this); }
		inline void stitch(const int nStitch, const View * pView)
			{ stitchAll(nStitch, pView, pView, pView, pView); }
		void stitchAll(const int nStitch)
			{ stitchAll(nStitch, this, this, this, this); }
		void stitchAll(const int nStitch, const View * pLeft, const View * pTop, const View * pRight, const View * pBottom)
		{
			assert(!(nStitch & StitchTop) || pTop != NULL);
			assert(!(nStitch & StitchBottom) || pBottom != NULL);
			assert(!(nStitch & StitchLeft) || pLeft != NULL);
			assert(!(nStitch & StitchRight) || pRight != NULL);

			if (nStitch & StitchTop)
			{
				assert(pTop->width == width && pTop->height == height);
				_pBlock->memcopy(line(0), pTop->_pBlock, pTop->line(bottom), width);
			}
			if (nStitch & StitchBottom)
			{
				assert(pBottom->width == width && pBottom->height == height);
				_pBlock->memcopy(line(bottom), pBottom->_pBlock, pBottom->line(0), width);
			}

			if (nStitch & StitchLeft)
			{
				assert(pLeft->width == width && pLeft->height == height);
				int j;

#ifdef _OPENMP
#pragma omp parallel for private(j)
#endif
				for (j = 0; j < height; ++j)
					set(0, j, pLeft->get(right, j));
			}
			if (nStitch & StitchRight)
			{
				assert(pRight->width == width && pRight->height == height);
				int j;

#ifdef _OPENMP
#pragma omp parallel for private(j)
#endif
				for (j = 0; j < height; ++j)
					set(right, j, pRight->get(0, j));
			}
		}

		void copyEdges(const int nEdges, const View * pOther)
		{
			assert(width == pOther->width && height == pOther->height);

			if (nEdges & StitchTop)
				_pBlock->memcopy(line(0), pOther->_pBlock, pOther->line(0), width);
			if (nEdges & StitchBottom)
				_pBlock->memcopy(line(bottom), pOther->_pBlock, pOther->line(bottom), width);

			if (nEdges & StitchLeft)
			{
				int j;

#ifdef _OPENMP
#pragma omp parallel for private(j)
#endif
				for (j = 0; j < height; ++j)
					set(0, j, pOther->get(0, j));
			}
			if (nEdges & StitchRight)
			{
				int j;

#ifdef _OPENMP
#pragma omp parallel for private(j)
#endif
				for (j = 0; j < height; ++j)
					set(right, j, pOther->get(right, j));
			}
		}

		inline bool isInside (const signed int x, const signed int y) const
		{
			return x >= left && y >= top && x <= right && y <= bottom;
		}

		// Samples elevation in a rotated Von Neumann neighborhood
		inline T meanVN (const signed int x, const signed int y) const
		{
			typedef typename ST::UP R;

			return static_cast <T> (
				(
					static_cast <R> (get(x - 1, y - 1)) +
					static_cast <R> (get(x + 1, y - 1)) +
					static_cast <R> (get(x - 1, y + 1)) +
					static_cast <R> (get(x + 1, y + 1))
				)
				/ 4);
		}

		inline T mean5 (const signed int x, const signed int y) const
		{
			typedef typename ST::UP R;
			return static_cast <T> (
				(
					static_cast <R> (get(x - 1, y - 1)) +
					static_cast <R> (get(x, y - 1)) +
					static_cast <R> (get(x + 1, y - 1)) +
					static_cast <R> (get(x - 1, y)) +
					static_cast <R> (get(x, y)) +
					static_cast <R> (get(x + 1, y)) +
					static_cast <R> (get(x - 1, y + 1)) +
					static_cast <R> (get(x, y + 1)) +
					static_cast <R> (get(x + 1, y + 1)) +

					static_cast <R> (get(x - 2, y - 2)) +
					static_cast <R> (get(x - 2, y - 1)) +
					static_cast <R> (get(x - 2, y)) +
					static_cast <R> (get(x - 2, y + 1)) +
					static_cast <R> (get(x - 2, y + 2)) +

					static_cast <R> (get(x + 2, y - 2)) +
					static_cast <R> (get(x + 2, y - 1)) +
					static_cast <R> (get(x + 2, y)) +
					static_cast <R> (get(x + 2, y + 1)) +
					static_cast <R> (get(x + 2, y + 2)) +

					static_cast <R> (get(x - 1, y - 2)) +
					static_cast <R> (get(x, y - 2)) +
					static_cast <R> (get(x + 1, y - 2)) +

					static_cast <R> (get(x - 1, y + 2)) +
					static_cast <R> (get(x, y + 2)) +
					static_cast <R> (get(x + 1, y + 2))
				) / 25
			);
		}

		inline T mean (const signed int x, const signed int y) const
		{
			typedef typename ST::UP R;

			R a = 0;

			a += get(x-1, y-1);	a += get(x+0, y-1);	a += get(x+1, y-1);
			a += get(x-1, y-0);	a += get(x+0, y+0);	a += get(x+1, y-0);
			a += get(x-1, y+1);	a += get(x+0, y+1);	a += get(x+1, y+1);

			return static_cast <T> (a / 9);
		}

		static void blitw (
			const ScalarFieldView & sfvSrc,		ScalarFieldView & sfvDest,
			const signed int x = 0,				const signed int y = 0
		)
		{
			const ScalarFieldView
				* psfvBig,
				* psfvSmall;

			determineOrder(sfvSrc, sfvDest, psfvSmall, psfvBig);

			assert(psfvBig->width >= psfvSmall->width && psfvBig->height >= psfvSmall->height);
			assert(
				sfvSrc.left >= 0 && sfvDest.left >= 0 && sfvSrc.top >= 0 && sfvDest.top >= 0 &&
				sfvSrc.right >= 0 && sfvDest.right >= 0 && sfvSrc.bottom >= 0 && sfvDest.bottom >= 0
			);

			const unsigned short
				w = psfvSmall->width,
				h = psfvSmall->height;

			const unsigned short
				x0 = WRAP(x, psfvBig->width),
				y0 = WRAP(y, psfvBig->height),

				xN = x0 + w - 1,
				yN = y0 + h - 1;

			if (xN >= psfvBig->width || yN >= psfvBig->height)
			{
				const unsigned short
					xH = std::min(xN, static_cast< unsigned short > (psfvBig->width - 1)),
					yH = std::min(yN, static_cast< unsigned short > (psfvBig->height - 1)),

					w0 = xH - x0 + 1,
					h0 = yH - y0 + 1,

					wN = psfvSmall->width - w0,
					hN = psfvSmall->height - h0;

				assert(psfvBig->gridtype == mars::Wrap);
				if (psfvBig == &sfvDest)
				{
					if (w0 > 0 && h0 > 0)	blitn (sfvSrc, sfvDest, w0, h0, 0,	0,	x0,	y0);
					if (wN > 0 && h0 > 0)	blitn (sfvSrc, sfvDest, wN, h0, w0, 0,	0,	y0);
					if (w0 > 0 && hN > 0)	blitn (sfvSrc, sfvDest,	w0, hN, 0,	h0, x0,	0);
					if (wN > 0 && hN > 0)	blitn (sfvSrc, sfvDest,	wN, hN, w0, h0, 0,	0);
				} else
				{
					if (w0 > 0 && h0 > 0)	blitn (sfvSrc, sfvDest,	w0, h0, x0,	y0, 0,	0);
					if (wN > 0 && h0 > 0)	blitn (sfvSrc, sfvDest,	wN, h0, 0,	y0, w0, 0);
					if (w0 > 0 && hN > 0)	blitn (sfvSrc, sfvDest, w0, hN, x0,	0,	0,	h0);
					if (wN > 0 && hN > 0)	blitn (sfvSrc, sfvDest, wN, hN, 0,	0,	w0, h0);
				}
			} else
			{
				if (psfvBig == &sfvDest)
					blitn(sfvSrc, sfvDest, psfvSmall->width, psfvSmall->height, 0, 0, x0, y0);
				else
					blitn(sfvSrc, sfvDest, psfvSmall->width, psfvSmall->height, x0, y0, 0, 0);
			}
		}

		static void blitn (
			const ScalarFieldView & sfvSrc,		ScalarFieldView & sfvDest,
			const unsigned short w,				const unsigned short h,
			const unsigned short xs = 0,		const unsigned short ys = 0,
			const unsigned short xd = 0,		const unsigned short yd = 0
		)
		{
			assert(w > 0 && h > 0);
			assert(w <= sfvSrc.width && h <= sfvSrc.height);
			assert(w <= sfvDest.width && h <= sfvDest.height);

			assert(xs >= 0 && ys >= 0);
			assert(xd >= 0 && yd >= 0);
			assert(
				(xs + w <= sfvSrc.width && ys + h <= sfvSrc.height) &&
				(xd + w <= sfvDest.width && yd + h <= sfvDest.height)
			);

			T *				pnDestOffset = sfvDest._pBlock->offset(xd + sfvDest.left, yd + sfvDest.top);
			const T *		pnSrcOffset = sfvSrc._pBlock->offset(xs + sfvSrc.left, ys + sfvSrc.top);
			const size_t	nCount = w;

			int j;
#ifdef _OPENMP
#pragma omp parallel for private(j)
#endif
			for (j = 0; j < h; ++j)
				sfvDest.block() ->memcopy(
					&pnDestOffset[j * sfvDest._pBlock->scanline],
					sfvSrc._pBlock,
					&pnSrcOffset[j * sfvSrc._pBlock->scanline],
					nCount
				);
		}

		inline void blit (ScalarFieldView & sfvDest, const signed int x = 0, const signed int y = 0) const
		{
			if (gridtype == Wrap || sfvDest.gridtype == Wrap)	// TODO: Confirm bounds too
				blitw(*this, sfvDest, x, y);
			else
			{
				const ScalarFieldView
					* psfvBig,
					* psfvSmall;

				determineOrder(*this, sfvDest, psfvSmall, psfvBig);
				if (psfvBig == &sfvDest)
					blitn(*this, sfvDest, psfvSmall->width, psfvSmall->height, 0, 0, x, y);
				else
					blitn(*this, sfvDest, psfvSmall->width, psfvSmall->height, x, y, 0, 0);
			}
		}

		void operator >> (std::vector< std::vector< T > > & stdfield) const
		{
			stdfield.resize(height);

			size_t j;
			PARALLEL_FOR(j);
			for (j = 0; j < height; ++j)
			{
				stdfield[j].resize(width);
				const T * pLine = line(j);
				for (size_t i = 0; i < width; ++i)
					stdfield[j][i] = pLine[i];
			}
		}

		inline void operator >> (ScalarFieldView & dest) const
		{
			assert(width == dest.width && height == dest.height);
			blit(dest);
		}
		inline void operator << (const ScalarFieldView & other)
		{
			assert(width == other.width && height == other.height);
			other.blit(*this);
		}
		inline void operator >> (Offset & dest) const
		{
			blit(*dest.view, dest.x, dest.y);
		}
		inline void operator << (const Offset & other)
		{
			other.const_view->blit(*this, other.x, other.y);
		}

		template< typename J, typename JST >
		void operator << (const ScalarFieldView< J, JST > & other)
		{
			assert(width == other.width && height == other.height);
			typedef ScalarFieldView< J, JST > OtherView;
			DUAL_KERNEL_OPERATION_T(OtherView, other,nDest,nSrc,{
				nDest = static_cast< T > (nSrc);
			});
		}

		template< typename J, typename JST >
		inline void operator >> (ScalarFieldView< J, JST > & dest) const
			{ dest << *this; }
	};

	template< typename T, typename ST = Numerics< T > >
	class ScalarField : public ScalarFieldView< T, ST >
	{
	private:
		typedef ScalarFieldView< T, ST > Supa;

	public:
		ScalarField (const unsigned short width, const unsigned short height, const ScalarField & copy, const unsigned short x, const unsigned short y, const GridType engt = Normal)
			: Supa(Supa::createHosted(width, height, engt))
			{ blitFrom (copy, x, y); }

		ScalarField (const unsigned short width, const unsigned short height, const GridType engt = Normal)
			: Supa(Supa::createHosted(width, height, engt))
			{ memset(this->data(), 0, this->block() -> bytesize); }
		ScalarField (const Supa & copy)
			: Supa(copy, new MemoryBlock< T >(copy.getMemoryBlock()))
			{}
		ScalarField (const ScalarField & copy)
			: Supa(copy, new MemoryBlock< T >(copy.getMemoryBlock()))
			{}

		inline std::istream & operator << (std::istream & ins)
			{ return *this->block() << ins; }
		inline std::ostream & operator >> (std::ostream & outs) const
			{ return *this->block() >> outs; }
		inline void operator >> (Supa & dest) const
			{ Supa::operator >> (dest); }
		inline void operator << (const Supa & other)
			{ Supa::operator << (other); }
		inline void operator >> (typename Supa::Offset dest) const
			{ Supa::operator >> (dest); }
		inline void operator << (const typename Supa::Offset other)
			{ Supa::operator << (other); }
		template< typename J, typename JST >
		inline void operator << (const ScalarFieldView< J, JST > & other)
			{ Supa::operator << (other); }
		template< typename J, typename JST >
		inline void operator >> (ScalarFieldView< J, JST > & dest) const
			{ Supa::operator >> (dest); }

		virtual ~ScalarField() { delete this->block(); }
	};

	template <typename T, typename ST = Numerics< T > >
	class MidpointDisplacementGrid : public ScalarFieldView< T, ST >
	{
	private:
	    typedef MidpointDisplacementGrid< T, ST > Self;
	    typedef ScalarFieldView< T, ST > Supa;

		const enum StorageMode
		{
			Concrete,
			Virtual
		} _ensmStorageMode;

		unsigned short _c;
		unsigned int _stride;

		static inline unsigned short createMinDim (const ScalarFieldView< T, ST > & sfv)
			{ return std::min(sfv.width, sfv.height); }
		static inline unsigned int createStride (const ScalarFieldView< T, ST > & sfv)
		{
			const unsigned short nMinDim = createMinDim(sfv);
			return nMinDim - nMinDim % 2;
		}
		static inline unsigned short createMinDim (const unsigned short nWidth, const unsigned short nHeight)
			{ return std::min(nWidth, nHeight); }
		static inline unsigned int createStride (const unsigned short nWidth, const unsigned short nHeight)
		{
			const unsigned short nMinDim = createMinDim(nWidth, nHeight);
			return nMinDim - nMinDim % 2;
		}

		inline void assertions () const
		{
			assert(this->width == dimensionFor(this->width, Supa::gridtype) && this->height == dimensionFor(this->height, Supa::gridtype));
			assert(_stride >= 2);
		}

	protected:
		MidpointDisplacementGrid (const short nIterCount, const unsigned short nWidth, const unsigned short nHeight, const GridType engt = mars::Normal)
			: Supa (Supa::createHosted(nWidth, nHeight, engt)), _c (0),
			_stride(createStride(nWidth, nHeight)),
			itercount(nIterCount), lastX(nWidth - 1), lastY(nHeight - 1),  mindim(createMinDim(nWidth, nHeight)), _ensmStorageMode(Virtual)
		{ assertions(); }
		MidpointDisplacementGrid (const short nIterCount, ScalarFieldView< T, ST > & sfv)
			: Supa (sfv), _c (0), itercount(nIterCount), lastX(sfv.width - 1), lastY(sfv.height - 1), _stride(createStride(sfv)), mindim(createMinDim(sfv)), _ensmStorageMode(Concrete)
		{ assertions(); }

		inline static short calcIterCount (const unsigned short nMinDim)
		{
			return static_cast <unsigned short> (
				ceil(
					log(static_cast <float> (nMinDim)) / log(2.0f)
				)
			);
		}
		inline static unsigned short calcDimLast (const unsigned short nIterCount)
		{
			return 1 << nIterCount;
		}

	public:
		template <typename S, unsigned dBp = 1>
		class AbstractFunctor
		{
		private:
			template <unsigned N>
			inline const unsigned int pow(const unsigned short c, quantity<N>) const
			{ return c * pow(c, quantity<N-1>()); }

			inline const unsigned int pow (const unsigned short c, quantity<0>) const
			{ return 1; }

		protected:
			const unsigned short _nFlux;
			const float _fCoarseness;

			float _nRandAmount;

			template <typename REAL>
			inline REAL applyFlux (const REAL & rMean) const
			{
				return rMean + static_cast <REAL> (mars::RANDf(_nRandAmount)) - static_cast <REAL> (_nRandAmount) / static_cast <REAL> (2);
			}
			inline float computeRandomAmount (const unsigned short c) const
			{
				return static_cast <float> (_nFlux) / static_cast< float > (pow<dBp>(c + 1, quantity<dBp>())) * _fCoarseness;
			}

		public:
			typedef S StateType;

			inline AbstractFunctor (const unsigned short nFlux, const float fCoarseness)
				: _nFlux (nFlux), _fCoarseness (fCoarseness) {}

			inline const float & getRandomAmount () const { return _nRandAmount; }

			inline void iteration (const unsigned short nStride, const unsigned short c)
			{
				//_nRandAmount = _nDim / (static_cast <int> (static_cast <float> (c) * (1.0f / _fCoarseness)) + 1);
				//_nRandAmount = static_cast <float> (_nStride) * _fCoarseness;
				// XTODO: Struggling with finding an appropriate noise generating algorithm here.
				// The problem is it's either too noisy at the finest level resulting in what look like innumerable spires everywhere
				// or the crater doesn't look noisy enough and looks like a perfectly uniform ring
				_nRandAmount = computeRandomAmount(c);
			}

			inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c, const T & d) const;
			inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c) const;
		};

		class PinkFn : public AbstractFunctor<int>
		{
		public:
			inline PinkFn (const unsigned short nFlux, const float fCoarseness)
				: AbstractFunctor<int>(nFlux, fCoarseness) {}

			inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c, const T & d) const
			{
				return static_cast <T> (AbstractFunctor<int>::applyFlux(static_cast <float> (a + b + c + d) / 4.0f));
			}
			inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c) const
			{
				return static_cast <T> (AbstractFunctor<int>::applyFlux(static_cast <float> (a + b + c) / 3.0f));
			}
		};

		class BrownianFn : public AbstractFunctor<int, 2>
		{
		public:
			inline BrownianFn (const unsigned short nFlux, const float fCoarseness)
				: AbstractFunctor<int, 2>(nFlux, fCoarseness) {}

			inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c, const T & d) const
			{
				return static_cast <T> (AbstractFunctor<int, 2>::applyFlux(static_cast <float> (a + b + c + d) / 4.0f));
			}
			inline T operator () (int & state, const int x, const int y, const T & a, const T & b, const T & c) const
			{
				return static_cast <T> (AbstractFunctor<int, 2>::applyFlux(static_cast <float> (a + b + c) / 3.0f));
			}
		};

		class EdgeSeedIterator : public std::iterator <std::input_iterator_tag, T>
		{
		private:
			MidpointDisplacementGrid< T > & _grid;
			unsigned short _j, _i;
			unsigned short _halfstride;
			bool _done;

			CR_CONTEXT;

			void process ()
			{
				// NOTE: Redundant code in MidpointDisplacementGrid::iterate(...);
				CR_START();

				_halfstride = _grid.stride() / 2;
				assert(_halfstride > 0);

				for (_j = 0; _j <= _grid.lastY; _j += _grid.stride())
					for (_i = 0; _i <= _grid.lastX; _i += _grid.stride())
					{
						CR_RETURN_VOID;
					}

				for (_j = 0; _j <= _grid.lastY; _j += _halfstride)
					for (_i = (_halfstride + _j) % _grid.stride(); _i <= _grid.lastX; _i += _grid.stride())
					{
						CR_RETURN_VOID;
					}

				_done = true;
				CR_END();
			}

		public:
			inline EdgeSeedIterator (MidpointDisplacementGrid & grid)
				: _grid(grid), _done(false)
			{
				CR_INIT();
				process();
			}

			inline EdgeSeedIterator & operator++ ()
			{
				process();
				return *this;
			}
			inline EdgeSeedIterator operator++ (int)
			{
				EdgeSeedIterator old = *this;

				process();
				return old;
			}

			inline T & assign (const T & a) { return _grid(_i, _j) = a; }
			inline T & ref () { return _grid(_i, _j); }
			inline const T & ref () const { return _grid(_i, _j); }

			inline T & operator = (const T & a) { return assign(a); }

			inline const unsigned short getX() const { return _i; }
			inline const unsigned short getY() const { return _j; }

			inline operator bool () const { return !_done; }
		};

	public:
		const unsigned short itercount;
		const unsigned short mindim, lastX, lastY;

		virtual ~MidpointDisplacementGrid()
		{
			switch (_ensmStorageMode)
			{
			case Concrete: delete this->block(); break;
			}
		}

		const unsigned short index () const { return _c; }
		const unsigned short stride () const { return _stride; }

		inline static Self * createInstance (const unsigned short nMinWidth, const unsigned short nMinHeight, const GridType engt)
		{
			assert(nMinHeight > 0 && nMinWidth > 0);
			return new Self (
				calcIterCount(std::min(nMinHeight, nMinWidth)),
				dimensionFor(nMinWidth, engt),
				dimensionFor(nMinHeight, engt),
				engt
			);
		}
		inline static Self * createInstance (Supa & sfv)
		{
			return new Self (
				calcIterCount(std::min(sfv.width, sfv.height) - (sfv.gridtype == Normal ? 1 : 0)),
				sfv
			);
		}

		inline static unsigned short dimensionFor (const unsigned short n, const GridType engt)
		{
			switch (engt)
			{
			case Wrap:
				return calcDimLast (calcIterCount (n));
			case Normal:
			default:
				return calcDimLast (calcIterCount (n - n % 2)) + 1;
			}
		}

		inline void run (const float fCoarseNess, const int nStitch = 0)
		{ run (BrownianFn (mindim, fCoarseNess), nStitch); }

		template <typename F>
		void run (F functor, const int nStitch = 0)
		{
			iterations (itercount, functor, nStitch);
		}
		template <typename F>
		void iterations (const unsigned short nCount, F & functor, const int nStitch = 0)
		{
			unsigned short nLastExcl = _c + nCount;

			while (_c <= nLastExcl && _stride >= 2)
			{
				iterate (functor, nStitch);
			}
		}

		inline EdgeSeedIterator edgeSeedIterator () { return EdgeSeedIterator (*this); }

		void seedFrom (const ScalarFieldView< T, ST > & other, const unsigned short xOffs, const unsigned short yOffs, const unsigned short nLOD)
		{
			int i, j, x, y;

			_c = itercount - nLOD;
			_stride = (1 << nLOD);

			assert(this->height / _stride + yOffs <= other.height);
			assert(this->width / _stride + xOffs <= other.width);

#ifdef _OPENMP
#pragma omp parallel for private(i, j, x, y)
#endif
			for (j = 0; j < this->height; j += _stride)
			{
				y = yOffs + (j >> nLOD);
				for (i = 0, x = xOffs; i < this->width; i += _stride, ++x)
					Self::operator () (i, j) = other(x, y) << nLOD;
			}
		}

		inline void reset ()
			{ step(-_c); }

		inline void step (const signed int nSteps = 1)
		{
			_c += nSteps;
			if (nSteps < 0)
				_stride *= (1 << -nSteps);
			else
				_stride /= (1 << nSteps);
		}

		template < typename Functor > void iterate (Functor & func, const int nStitch = 0)
		{
			const int
				nHalfStride = _stride / 2;

			int i, j;
			typename Functor::StateType state;

			assert(nHalfStride > 0);
			func.iteration (_stride, _c);

			switch (this->gridtype)
			{
			case Wrap:
				// Square phase
				// NOTE: Redundant code in EdgeSeedIterator
#ifdef _OPENMP
#pragma omp parallel for private(i, j, state)
#endif
				for (j = 0; j <= lastY; j += _stride)
					for (i = 0; i <= lastX; i += _stride)
						Self::operator()(i + nHalfStride, j + nHalfStride) =
							func(state,
								i + nHalfStride, j + nHalfStride,
								Self::operator() (i, j),
								Self::operator() ((i + _stride) % this->width, j),
								Self::operator() (i, (j + _stride) % this->height),
								Self::operator() ((i + _stride) % this->width, (j + _stride) % this->height)
							);

				// Diamond phase
				// NOTE: Redundant code in EdgeSeedIterator and other iterator
#ifdef _OPENMP
#pragma omp parallel for private(i, j, state)
#endif
				for (j = 0; j <= lastY; j += nHalfStride)
					for (i = (nHalfStride + j) % _stride; i <= lastX; i += _stride)
						Self::operator() (i, j) =
							func(state,
								i, j,
								Self::operator()(i, (j - nHalfStride + this->height) % this->height),
								Self::operator()(i, (j + nHalfStride) % this->height),
								Self::operator()((i - nHalfStride + this->width) % this->width, j),
								Self::operator()((i + nHalfStride) % this->width, j)
							);
				break;
			case Normal:
			default:
				{
					const int
						nCtxLastX = lastX - nHalfStride,
						nCtxLastY = lastY - nHalfStride;

					// Square-phase
					// NOTE: Redundant code in EdgeSeedIterator
#ifdef _OPENMP
#pragma omp parallel for private(i, j, state)
#endif
					for (j = 0; j < nCtxLastY; j += _stride)
						for (i = 0; i < nCtxLastX; i += _stride)
							Self::operator()(i + nHalfStride, j + nHalfStride) =
								func(state,
									i + nHalfStride, j + nHalfStride,
									Self::operator() (i, j),
									Self::operator() (i + _stride, j),
									Self::operator() (i, j + _stride),
									Self::operator() (i + _stride, j + _stride)
								);

#ifdef _OPENMP
#pragma omp parallel sections private(i, j, state)
#endif
					{
						// Next 4 blocks explicitly perform the square-step (second-step) of the diamond-
						// square algorithm	along all four edges of the grid so that we don't have to do
						// any bounds checking for the rest of the iteration

						// Top edge
#ifdef _OPENMP
#pragma omp section
#endif
						if ((nStitch ^ StitchTop) & StitchTop)
						{
							const int jj = 0;
							for (i = nHalfStride; i < lastX; i += _stride)
								Self::operator() (i, jj) =
									func(state,
										i, jj,
										Self::operator()(i - nHalfStride, jj),
										Self::operator()(i + nHalfStride, jj),
										Self::operator()(i, jj + nHalfStride)
									);
						}

						// Bottom edge
#ifdef _OPENMP
#pragma omp section
#endif
						if ((nStitch ^ StitchBottom) & StitchBottom)
						{
							const int j = lastY;
							for (i = nHalfStride; i < lastX; i += _stride)
								Self::operator() (i, j) =
									func(state,
										i, j,
										Self::operator()(i - nHalfStride, j),
										Self::operator()(i + nHalfStride, j),
										Self::operator()(i, j - nHalfStride)
									);
						}

						// Left edge
#ifdef _OPENMP
#pragma omp section
#endif
						if ((nStitch ^ StitchLeft) & StitchLeft)
						{
							const int i = 0;
							for (j = nHalfStride; j < lastY; j += _stride)
								Self::operator() (i, j) =
									func(state,
										i, j,
										Self::operator()(i, j - nHalfStride),
										Self::operator()(i, j + nHalfStride),
										Self::operator()(i + nHalfStride, j)
									);
						}

						// Right edge
#ifdef _OPENMP
#pragma omp section
#endif
						if ((nStitch ^ StitchRight) & StitchRight)
						{
							const int i = lastX;
							for (j = nHalfStride; j < lastY; j += _stride)
								Self::operator() (i, j) =
									func(state,
										i, j,
										Self::operator()(i, j - nHalfStride),
										Self::operator()(i, j + nHalfStride),
										Self::operator()(i - nHalfStride, j)
									);
						}
					}

					// Perform the second phase (square phase) of the diamond-square algorithm on the rest of
					// the grid excluding its edges.
					// NOTE: Redundant code in EdgeSeedIterator
#ifdef _OPENMP
#pragma omp parallel for private(i, j, state)
#endif
					for (j = nHalfStride; j < lastY; j += nHalfStride)
						for (i = nHalfStride + j % _stride; i < lastX; i += _stride)
							Self::operator() (i, j) =
								func(state,
									i, j,
									Self::operator()(i, j - nHalfStride),
									Self::operator()(i, j + nHalfStride),
									Self::operator()(i - nHalfStride, j),
									Self::operator()(i + nHalfStride, j)
								);
				}
				break;
			}

			++_c;
			_stride /= 2;
		}

		inline std::istream & operator << (std::istream & ins)
			{ return *this->block() << ins; }
		inline std::ostream & operator >> (std::ostream & outs) const
			{ return *this->block() >> outs; }

	};
}

#endif