#include "geoworld.h"

#include <mars_calc.h>
#include <mars_util.h>
#include <mars_lookup.h>
#include <math.h>
#include <stdio.h>

#include "geomorph.h"
#include "lodtools.h"

namespace geoworld
{
	void LockedGWSection::blitFrom( const FileSource::MainDEM::Block & block, const unsigned short nTileX, const unsigned short nTileY )
	{ 
		block.copyTo(_hm, nTileX, nTileY); 
	}

	void LockedGWSection::blitTo( FileSource::MainDEM::Block & block, const unsigned short nTileX, const unsigned short nTileY ) const
		{ block.copyFrom (_hm, nTileX, nTileY); }

	void LockedGWSection::operator >> ( PNGDEM & pngdem ) const
		{ pngdem << _hm; }

	void LockedGWSection::operator>>( LockedGWSection::View && view ) const
	{
		if (view.getHeightMap() ->width == _hm.width && view.getHeightMap() ->height == _hm.height)
		{
			assert(view.getHeightMap() ->width == view.getStratum() ->width && view.getHeightMap() ->height == view.getStratum() ->height);
			_hm >> *view.getHeightMap();
			_stratum >> *view.getStratum();
		} else
		{
			assert(abs(bbox.left) < MAX_LOCK_DIM && abs(bbox.top) < MAX_LOCK_DIM);

			const signed int 
				x = static_cast< signed int > (bbox.left), 
				y = static_cast< signed int > (bbox.top);

			_hm >> GeoHeightMap::View::Offset(view.getHeightMap(), x, y);
			_stratum >> Stratum::View::Offset(view.getStratum(), x, y);
		}
	}

	void LockedGWSection::operator << ( const PNGDEM & pngdem )
	{
		assert(pngdem.getWidth() <= _hm.width && pngdem.getHeight() <= _hm.height);
		GeoHeightMap::Precision * pnData = new GeoHeightMap::Precision[pngdem.getWidth()];
		for (unsigned int i = 0; i < pngdem.getHeight(); ++i)
		{
			pngdem.retrieveLine(pnData, i);
			_hm.setn(0, i, pnData, pngdem.getWidth());
		}
	}

	void LockedGWSection::operator << ( const LockedGWSection::View & view )
	{
		if (view.getHeightMap()->width == _hm.width && view.getHeightMap()->height == _hm.height)
		{
			assert(view.getHeightMap() ->width == view.getStratum() ->width && view.getHeightMap() ->height == view.getStratum() ->height);
			_hm << *view.getHeightMap();
			_stratum << *view.getStratum();
		} else
		{
			assert(abs(bbox.left) < MAX_LOCK_DIM && abs(bbox.top) < MAX_LOCK_DIM);

			const signed int 
				x = static_cast< signed int > (bbox.left), 
				y = static_cast< signed int > (bbox.top);

			_hm << GeoHeightMap::View::Offset(view.getHeightMap(), x, y);
			_stratum << Stratum::View::Offset(view.getStratum(), x, y);
		}
	}

	mars::BBox< long > LockedGWSection::computeLockedSection( const mars::BBox< long > & bbox, const unsigned short nTileDim, const bool bPad /*= true*/ )
	{
		assert(((nTileDim - 1) & nTileDim) == 0); // Must be a power of 2

		using namespace mars;
		const unsigned short nPadHalfTile = (bPad ? nTileDim / 2 : 0);

		// DEPS: LGWSBBA
		return mars::BBox< long > (
			ALIGN(bbox.left - nPadHalfTile, static_cast< long > (nTileDim)),
			ALIGN(bbox.top - nPadHalfTile, static_cast< long > (nTileDim)),
			ALIGN(bbox.right + nPadHalfTile + nTileDim, static_cast< long > (nTileDim)) - 1,
			ALIGN(bbox.bottom + nPadHalfTile + nTileDim, static_cast< long > (nTileDim)) - 1
		);
	}

	bool LockedGWSection::isValid( const mars::BBox< long > & bbox, const unsigned short nTileDim, const bool bPad /*= true */ )
	{
		assert(((nTileDim - 1) & nTileDim) == 0); // Must be a power of 2
		if (bbox.empty())
			return false;

		const mars::BBox< long > bboxActual = computeLockedSection(bbox, nTileDim, bPad);
		return bboxActual.getWidth() <= static_cast< long > (MAX_LOCK_DIM) && bboxActual.getHeight() <= static_cast< long > (MAX_LOCK_DIM);
	}

	LockedGWSection::LockedGWSection( const mars::BBox< long > & bbox, const long nDepth, const unsigned short nDEMLOD, const size_t nLayers, const unsigned short nStratumLOD ) 
		: _hm(static_cast< unsigned short > (bbox.getWidth()), static_cast< unsigned short > (bbox.getHeight()), mars::Normal), 
		_stratum(
			nLayers, 
			static_cast< WorldUnit > (bbox.getWidth()) << nDEMLOD, 
			static_cast< WorldUnit > (bbox.getHeight()) << nDEMLOD, 
			static_cast< WorldUnit > (nDepth) << nDEMLOD,
			nStratumLOD, 
			mars::Normal
		), 
		bbox(bbox),
		_ssop(SSync_Perfect), _nDEMLOD(nDEMLOD)
	{
		assert(bbox.getWidth() <= static_cast< long > (MAX_LOCK_DIM) && bbox.getHeight() <= static_cast< long > (MAX_LOCK_DIM));
	}

	LockedGWSection::LockedGWSection( const GeoHeightMap & hm, const unsigned short nDEMLOD, const Stratum & stratum ) 
		: _hm(hm), _stratum(stratum.layers, hm.width << stratum.lod, hm.height << stratum.lod, stratum.absdepth, stratum.lod, stratum.getBase().gridtype), bbox(hm.left, hm.top, hm.right, hm.bottom),
		_ssop(SSync_Perfect), _nDEMLOD(nDEMLOD)
	{
		_stratum << Stratum::View::Offset(&stratum, bbox.left, bbox.top);
	}

	LockedGWSection::View * LockedGWSection::createView( const long nLeft, const long nTop, const long nRight, const long nBottom )
	{
		assert (bbox.inside(nLeft, nTop) && bbox.inside(nRight, nBottom));

		const unsigned short
			nLeft2 = static_cast< unsigned short > (nLeft - bbox.left),
			nTop2 = static_cast< unsigned short > (nTop - bbox.top),
			nRight2 = static_cast< unsigned short > (nRight - bbox.left),
			nBottom2 = static_cast< unsigned short > (nBottom - bbox.top);

		View 
			* pView = new View(
				_hm.createView (nLeft2, nTop2, nRight2, nBottom2),
				_stratum.createView(nLeft2, nTop2, nRight2, nBottom2)
			);

		LOG(mars::Log::Debug) << bbox << ", " << pView->getHeightMap();
		return pView;
	}

	const LockedGWSection::View * LockedGWSection::createView( const long nLeft, const long nTop, const long nRight, const long nBottom ) const
	{
		assert (bbox.inside(nLeft, nTop) && bbox.inside(nRight, nBottom));

		const unsigned short
			nLeft2 = static_cast< unsigned short > (nLeft - bbox.left),
			nTop2 = static_cast< unsigned short > (nTop - bbox.top),
			nRight2 = static_cast< unsigned short > (nRight - bbox.left),
			nBottom2 = static_cast< unsigned short > (nBottom - bbox.top);

		const View
			* pView = new View(
				_hm.createView (nLeft2, nTop2, nRight2, nBottom2),
				_stratum.createView(nLeft2, nTop2, nRight2, nBottom2)
			);

		LOG(mars::Log::Debug) << bbox << ", " << pView->getHeightMap();
		return pView;
	}

	PagedGeoWorld::PagedGeoWorld(DepositSynthesizer * pDepSynth, FileSource * pFile, const unsigned short nPasses)
		: _pFile(pFile), _pAgentMaps(NULL), _passes(nPasses), _nNumProcs(DEFAULT_NUM_PROCS), _plsbboxLockedMutableRegions(NULL),
		_pDepSynth(pDepSynth)
	{
		const unsigned long
			nWorldWidth = pFile->getVirtualWidth(),
			nWorldHeight = pFile->getVirtualHeight();

		assert(_plsbboxLockedMutableRegions == NULL && _pAgentMaps == NULL);
		assert(_pDepSynth->getStratum().getBase().gridtype == mars::Wrap);

		_pFile = pFile;
		for (PassList::iterator i = _passes.begin(); i != _passes.end(); ++i)
		{
			delete *i;
			*i = new Pass(nWorldWidth, nWorldHeight);
		}
		_plsbboxLockedMutableRegions = new RegionContainer(nWorldWidth, nWorldHeight);
		_pAgentMaps = new AllAgentMaps(nWorldWidth, nWorldHeight);
	}

	PagedGeoWorld::~PagedGeoWorld()
	{
		GW_SCOPED_LOCK(mutex);
		{
			delete _plsbboxLockedMutableRegions;
			delete _pAgentMaps;

			for (PassList::iterator i = _passes.begin(); i != _passes.end(); ++i)
				delete *i;
		}
	}

	template< PagedGeoWorld::Mode MODE, typename TYPE_LockedGWSection >
	void PagedGeoWorld::performIOOperation( TYPE_LockedGWSection * pSection ) const
	{
		assert(_pFile != NULL);
		FileSource::MainDEM::Block block (_pFile->tilemap.getVirtualTileDim());

		const unsigned int 
			nTX0 = static_cast< unsigned int > (mars::WRAP(pSection->bbox.left, _pFile->tilemap.getVirtualWidth()) / _pFile->tilemap.getVirtualTileDim()),
			nTY0 = static_cast< unsigned int > (mars::WRAP(pSection->bbox.top, _pFile->tilemap.getVirtualHeight()) / _pFile->tilemap.getVirtualTileDim()),
			nTW = (pSection->getActualWidth() - 1) / _pFile->tilemap.getVirtualTileDim(),	// - 1 accounts for non-even redundant tile-blocking
			nTH = (pSection->getActualHeight() - 1) / _pFile->tilemap.getVirtualTileDim();

		if (nTX0 >= 0 && nTX0 + nTW <= _pFile->tilemap.getTilesWide() &&
			nTY0 >= 0 && nTY0 + nTH <= _pFile->tilemap.getTilesHigh())
		{
			seekTileOp<MODE>( nTX0, nTY0, std::ios::beg );

			unsigned short j = 0;

			while (j < nTH)
			{
				for (unsigned short i = 0; i < nTW; ++i)
					blockOp<MODE>(pSection, block, i, j);

				if (++j < nTH)
					seekTileOp<MODE>(-static_cast< signed int > (nTW), 1);
			}
		} else
		{
			unsigned int 
				nTXI, nTYI;
			bool bThumb;

			seekTileOp<MODE>(nTX0, nTY0, std::ios::beg);
			nTYI = 0;
			while (nTY0 + nTYI < _pFile->tilemap.getTilesHigh() && nTYI < nTH)
			{
				nTXI = 0;
				while (nTXI < nTW && nTX0 + nTXI < _pFile->tilemap.getTilesWide())
				{
					blockOp<MODE>(pSection, block, nTXI++, nTYI);
				}

				if (bThumb = (nTXI < nTW))
				{
					seekTileOp<MODE>(0, -1);
					do
					{
						blockOp<MODE>(pSection, block, nTXI++, nTYI);
					} while (nTXI < nTW);
				}

				++nTYI;
				if (nTY0 + nTYI < _pFile->tilemap.getTilesHigh() && nTYI < nTH)
					seekTileOp<MODE>(-static_cast< signed int > (nTW), 1 + (bThumb ? 1 : 0));
				else
					break;
			}
			seekTileOp<MODE>(nTX0, 0, std::ios::beg);
			while (nTYI < nTH)
			{
				nTXI = 0;
				while (nTXI < nTW && nTX0 + nTXI < _pFile->tilemap.getTilesWide())
				{
					blockOp<MODE>(pSection, block, nTXI++, nTYI);
				}
				if (bThumb = (nTXI < nTW))
				{
					seekTileOp<MODE>(0, -1);
					do
					{
						blockOp<MODE>(pSection, block, nTXI++, nTYI);
					} while (nTXI < nTW);
				}

				if (++nTYI < nTH)
					seekTileOp<MODE>(-static_cast< signed int > (nTW), 1 + (bThumb ? 1 : 0));
				else
					break;
			}
		}
	}

	void PagedGeoWorld::unlock( LockedGWSection * pSection )
	{
		using namespace std;
		using namespace mars;

		GW_SCOPED_LOCK(mutex);
		{
			assert(_pFile != NULL);
			assert(_plsbboxLockedMutableRegions != NULL);

			if (!_plsbboxLockedMutableRegions->removeRegion(pSection->bbox))
				throw UnlockEx(pSection->bbox);

			pSection->processSyncOp();
			performIOOperation<Write>(pSection);

			delete pSection;
		}
	}

	void PagedGeoWorld::unlock( const LockedGWSection * pSection ) const
	{
		using namespace std;
		using namespace mars;

		GW_SCOPED_LOCK(mutex);
		{
			assert(_pFile != NULL);
			assert(_plsbboxLockedMutableRegions != NULL);

			// TODO: This doesn't quite work, will fail when a r/w section was requested while this region was locked and the regions overlap/intersect
			if (_plsbboxLockedMutableRegions->intersects(pSection->bbox))
				throw UnlockEx(pSection->bbox);

			delete pSection;
		}
	}

	geoworld::LockedGWSection * PagedGeoWorld::doLock( const mars::BBox< long > & bbox, const bool bPad /*= true*/, const bool bReadOnly /*= false*/ ) const
	{
		using namespace std;
		using namespace mars;

		assert(static_cast< unsigned long > (bbox.getWidth()) <= getWorldWidth() && static_cast< unsigned long > (bbox.getHeight()) < getWorldHeight());

		GW_SCOPED_LOCK(mutex);
		{
			assert(_pFile != NULL);
			assert(_plsbboxLockedMutableRegions != NULL);

			const unsigned short nTileDim = _pFile->tilemap.getVirtualTileDim() >> getLOD();
			const mars::BBox< long > 
				bboxActual = LockedGWSection::computeLockedSection(bbox, nTileDim, bPad);

			if (_plsbboxLockedMutableRegions->intersects(bboxActual))
				return NULL;

			if (!bReadOnly)
				_plsbboxLockedMutableRegions->addRegion(bboxActual);

			LockedGWSection * pSection = new LockedGWSection(bboxActual, _pFile->getWorldDepth(), getLOD(), _pFile->geohost.count, _pFile->geohost.getLOD());
			performIOOperation<Read>(pSection);
			*pSection << Stratum::Offset(&_pDepSynth->getStratum(), bboxActual.left, bboxActual.top);
			return pSection;
		}
	}

	void PagedGeoWorld::processWorkQueue( PagedGeoWorld::WorkQueue & q, PagedGeoWorld::IWorkQueueListener * pListener /*= NULL*/ )
	{
		bool bOuterRunning = !q.inbox.empty();
		bool bInnerRunning = false;

		if (!bOuterRunning)
			return;

		#pragma omp parallel default(none) firstprivate(bOuterRunning, bInnerRunning, pListener) shared(q)
		do
		{
			#pragma omp critical(work)
			{
				#pragma omp master
				{
					WorkQueue::ItemQueue inbox, ready, collisions;

					#pragma omp critical(inbox)
					{ inbox = q.inbox; }

					while (!inbox.empty() && ready.size() < _nNumProcs)
					{
						WorkItem item = inbox.front();
						inbox.pop_front();

						item.section = lock(item.picks->getBBox());

						if (item.section != NULL)
							ready.push_back(item);
						else
							collisions.push_back(item);
					}

					#pragma omp critical(ready)
					{ 
						q.ready.insert(q.ready.end(), ready.begin(), ready.end()); 
						if (pListener != NULL && !ready.empty())
							pListener->readyQueue(ready);
					}

					#pragma omp critical(inbox)
					{ 
						q.inbox.clear();
						q.inbox.insert(q.inbox.end(), inbox.begin(), inbox.end());
						q.inbox.insert(q.inbox.end(), collisions.begin(), collisions.end()); 

						if (pListener != NULL && !collisions.empty())
							pListener->collisions(collisions);
					}
				}
			}

			do 
			{
				WorkItem item;

				#pragma omp critical (ready)
				{
					bInnerRunning = !q.ready.empty();
					if (bInnerRunning)
					{
						item = q.ready.front();
						q.ready.pop_front();
					}
				}

				if (bInnerRunning)
				{
					item.picks->consume(_pDepSynth, item.section, item.minerals);

					#pragma omp critical (outbox)
					{ q.outbox.push_back(item); }
				}
			} while (bInnerRunning);

			#pragma omp master
			{
				WorkQueue::ItemQueue outbox;

				#pragma omp critical(outbox)
				{
					outbox = q.outbox;
					q.outbox.clear();
				}

				if (pListener != NULL && !outbox.empty())
					pListener->unlockQueue(outbox);

				flushOutbox(outbox);
			}

			#pragma omp critical (inbox)
			{
				bOuterRunning = !q.inbox.empty();
			}
		} while (bOuterRunning);
		
		// Outside of parallel region, this ensures everything is released
		if (pListener != NULL && !q.outbox.empty())
			pListener->unlockQueue(q.outbox);
		flushOutbox(q.outbox);

		assert(_plsbboxLockedMutableRegions->empty());
	}

	void PagedGeoWorld::flushOutbox( WorkQueue::ItemQueue &outbox )
	{
		while(!outbox.empty())
		{
			unlock(outbox.front().section);
			flushMinerals(outbox.front().minerals);
			outbox.pop_front();
		}
	}

	void PagedGeoWorld::runMorphs( const unsigned short nPass, MorphList & rejects )
	{
		size_t c = 0;
		Pass & pass = *_passes[nPass];

		const long
			nWW = _pFile->tilemap.getVirtualWidth(),
			nWH = _pFile->tilemap.getVirtualHeight();
		const unsigned short
			nLW = static_cast< unsigned short > (std::min (static_cast< long > (LockedGWSection::MAX_LOCK_DIM), nWW / 4)),
			nLH = static_cast< unsigned short > (std::min (static_cast< long > (LockedGWSection::MAX_LOCK_DIM), nWH / 4));

		MasterSet master (&pass.quadtree, pass.list.begin(), pass.list.end());
		const size_t nTotal = master.size();

		class QueueListener : public IWorkQueueListener
		{
		private:
			PagedGeoWorld * _pParent;
			size_t _nTotal;

		public:
			size_t counter;

			QueueListener (PagedGeoWorld * pParent, const size_t nTotal)
				: _pParent(pParent), _nTotal(nTotal), counter(0) {}

			void readyQueue (const WorkQueue::ItemQueue & readyq)
			{
				MorphList mrphs;
				for (WorkQueue::ItemQueue::const_iterator i = readyq.begin(); i != readyq.end(); ++i)
					for (Picks::iterator j = i->picks->begin(); j != i->picks->end(); ++j)
						mrphs.push_back(*j);

				_pParent->fire(&IGeoWorldListener::runningMorphs, mrphs.begin(), mrphs.end(), counter, mrphs.size(), _nTotal);
				counter += mrphs.size();
			}
		} qListener(this, nTotal);

		long i, j;
		WorkQueue q;

		for (j = 0; j < nWW; j += nLH)
		{
			for (i = 0; i < nWH; i += nLW)
			{
				mars::ptr< Picks > pPicks;
				
				pPicks = 
					master.collect(
						QTMorphs::Region (
							QTMorphs::VECTOR(i, j), 
							QTMorphs::VECTOR(i + nLW - 1, j + nLH - 1)
						)
					);
				
				if (isPickValid(*pPicks))
					q.inbox.push_back(WorkItem(pPicks));
			}
		}

		processWorkQueue(q, &qListener);

		if (!master.empty())
		{
			for (MorphSet::const_iterator i = master.begin(); i != master.end(); )
			{
				mars::ptr< GeoMorph > pp = *i;
				BoundedGeoMorph * pBGM = mars::ptr_dynamic_cast< BoundedGeoMorph > (pp);

				if (pBGM != NULL)
				{
					mars::BBox< long > bbox = pBGM->getBBox();
					QTMorphs::VECTOR ptCenter = (bbox.getMaximum() - bbox.getMinimum()) / 2 + bbox.getMinimum();
					mars::ptr< Picks > pPicks = master.collect(
						QTMorphs::Region(
							QTMorphs::VECTOR(ptCenter.x - nLW / 2, ptCenter.y - nLH / 2), 
							QTMorphs::VECTOR(ptCenter.x + nLW / 2, ptCenter.y + nLH / 2)
						)
					);

					if (isPickValid(*pPicks))
						q.inbox.push_back(WorkItem(pPicks));

					if (!pPicks->empty())
						i = master.begin();	// Iterator is potentially invalid after calling "collect" and getting some picks, restart the list iteration
					else
						++i;
				} else
					++i;
			}

			processWorkQueue(q, &qListener);
		}

		if (!master.empty())
		{
			const unsigned short 
				nMax2LockWidth = 
					static_cast< unsigned short > (
						std::min(
							static_cast< unsigned long > (getWorldWidth() / 4), 
							static_cast< unsigned long > (1 << mars::LOG2(getMaxGeoMorphWidth()))
						)
					),
				nMax2LockHeight = 
					static_cast< unsigned short > (
						std::min(
							static_cast< unsigned long > (getWorldHeight() / 4), 
							static_cast< unsigned long > (1 << mars::LOG2(getMaxGeoMorphHeight()))
						)
					);

			const unsigned short
				n2iX = static_cast< unsigned short > (getWorldWidth() / nMax2LockWidth),
				n2iY = static_cast< unsigned short > (getWorldHeight() / nMax2LockHeight),
				n2Count = n2iX * n2iY;

			for (MorphSet::const_iterator i = master.begin(); i != master.end(); ++i)
			{
				mars::ptr< GeoMorph > pp = *i;
				BoundedGeoMorph * pBGM = mars::ptr_dynamic_cast< BoundedGeoMorph > (pp);

				if (pBGM != NULL)	// TODO: Support all GM types
				{
					const mars::BBox< long > bbox = pBGM->getBBox();

					if (isValidSectionBBox(bbox))
						q.inbox.push_back(WorkItem(pp));
					else
					{
						fire(&IGeoWorldListener::rejectedMorph, *i);
						rejects.push_back(*i);
					}
				} else
				{
					HeightMapGeoMorph * pGM = mars::ptr_dynamic_cast< HeightMapGeoMorph > (pp);

					if (pGM != NULL)
					{
						long x, y;
					
						for (unsigned short j = 0; j < n2iY; ++j)
						{
							y = static_cast< long > ((j * 2) % n2iY + (j * 2 / n2iY)) * nMax2LockHeight;
							for (unsigned short i = 0; i < n2iX; ++i)
							{
								x = static_cast< long > ((i * 2) % n2iX + (i * 2 / n2iX)) * nMax2LockWidth;
								q.inbox.push_back(
									WorkItem(
										pp, 
										mars::BBox< long > (
											x, 
											y, 
											x + nMax2LockWidth - 1, 
											y + nMax2LockHeight - 1
										)
									)
								);
							}
						}

						this->fire(&IGeoWorldListener::runningMorph, pp, qListener.counter++, nTotal);
						processWorkQueue(q);
					}
				}
			}
			processWorkQueue(q, &qListener);
		}
	}

	bool PagedGeoWorld::isValidSectionBBox (const mars::BBox< long > & bbox) const
	{
		if (bbox.empty())
			return false;

		if (bbox.getWidth() > getMaxGeoMorphWidth() || bbox.getHeight() > getMaxGeoMorphHeight())
			return false;

		return LockedGWSection::isValid(bbox, _pFile->getVirtualTileDim());
	}
	bool PagedGeoWorld::isPickValid (const Picks & p) const
	{
		return isValidSectionBBox(p.getBBox()) && !p.empty();
	}

	unsigned short PagedGeoWorld::getMaxGeoMorphWidth() const
	{ 
		return static_cast< unsigned short > (
			std::min(
				static_cast< unsigned long > (
					LockedGWSection::getMaxViewportDim(_pFile->tilemap.getVirtualTileDim())
				), 
				static_cast< unsigned long > (
					LockedGWSection::getActualViewportDim(
						static_cast< unsigned short > (
							std::min(
								_pFile->tilemap.getVirtualWidth(), 
								static_cast< unsigned long > (std::numeric_limits< unsigned short >::max())
							)
						), 
						_pFile->tilemap.getVirtualTileDim()
					)
				)
			) >> getLOD()
		); 
	}

	unsigned short PagedGeoWorld::getMaxGeoMorphHeight() const
	{ 
		return static_cast< unsigned short > (
			std::min(
				static_cast< unsigned long > (
					LockedGWSection::getMaxViewportDim(_pFile->tilemap.getVirtualTileDim())
				), 
				static_cast< unsigned long > (
					LockedGWSection::getActualViewportDim(
						static_cast< unsigned short > (
							std::min(
								_pFile->tilemap.getVirtualHeight(), 
								static_cast< unsigned long > (std::numeric_limits< unsigned short >::max())
							)
						), 
						_pFile->tilemap.getVirtualTileDim()
					)
				)
			) >> getLOD()
		); 
	}

	void PagedGeoWorld::addMorphs( MorphList & morphs, MorphList & rejects, const unsigned short nPass )
	{
		Pass & pass = *_passes[nPass];

		for (MorphList::const_iterator i = morphs.begin(); i != morphs.end(); ++i)
		{
			if (!checkGeoMorph(*i))
			{
				rejects.push_back(*i);
				i = morphs.erase(i);

				if (i == morphs.end())
					break;
			}
			else
				pass.quadtree.addMorph(*i);
		}
		pass.list.insert(pass.list.end(), morphs.begin(), morphs.end()); 
	}

	bool PagedGeoWorld::addMorph( mars::ptr <GeoMorph> & gm, const unsigned short nPass )
	{
		if (checkGeoMorph(gm))
		{
			Pass & pass = *_passes[nPass];
			pass.list.push_back(gm); 
			pass.quadtree.addMorph(gm);
			return true;
		}

		return false;
	}

	void PagedGeoWorld::readBlock( LockedGWSection * pSection, FileSource::MainDEM::Block & block, const unsigned short i, const unsigned short j ) const
	{
		block << _pFile->tilemap;
		pSection->blitFrom(block, i, j);
	}

	void PagedGeoWorld::writeBlock( const LockedGWSection * pSection, FileSource::MainDEM::Block & block, const unsigned short i, const unsigned short j ) const
	{
		pSection->blitTo(block, i, j);
		block >> _pFile->tilemap;
	}

	void PagedGeoWorld::flushMinerals( MineralsContainer & ctrMins )
	{
		GW_SCOPED_LOCK(mutex);
		{
			assert(_pFile != NULL);
			ctrMins.flushTo(_pFile->minerals);
		}
	}

	bool PagedGeoWorld::checkGeoMorph( const mars::ptr<GeoMorph> &gm )
	{
		const BoundedGeoMorph * pBGM = dynamic_cast< const BoundedGeoMorph * > (gm.pointer());

		if (pBGM != NULL)
		{
			if (!isValidSectionBBox(pBGM->getBBox()))
				return false;
		}

		return true;
	}

	void PagedGeoWorld::Picks::consume( const IMineralQuery * pQMin, LockedGWSection * pSection, MineralsContainer & ctrMins )
	{
		assert(!_bbox.empty());
		// TODO: Got an access violation in this loop, invalid iterator? list pulled-out from beneath by something else?
		for (iterator i = begin(); i != end(); ++i)
		{
			mars::ptr< GeoMorph > pp = *i;
			HeightMapGeoMorph * pSyGM = mars::ptr_dynamic_cast< HeightMapGeoMorph > (pp);

			if (pSyGM != NULL)
				pSyGM->doMorph(*pSection);

			MineralSpawnGeoMorph * pMSGM = mars::ptr_dynamic_cast< MineralSpawnGeoMorph > (pp);

			if (pMSGM != NULL)
				pMSGM->doSpawnMinerals(pQMin, ctrMins);
		}
	}

	void PagedGeoWorld::Picks::computeBBox()
	{
		_bbox = mars::BBox< long >();

		for (MorphSet::const_iterator k = begin(); k != end(); ++k)
		{
			const BoundedGeoMorph * pBGM = mars::ptr_dynamic_cast< BoundedGeoMorph > (*k);

			if (pBGM != NULL)
			{
				if (_bbox.empty())
					_bbox = pBGM->getBBox();
				else
					_bbox.add(pBGM->getBBox());
			}
		}
	}

	mars::ptr< PagedGeoWorld::Picks > PagedGeoWorld::MasterSet::collect( const QTMorphs::Region & rgBBox )
	{
		GW_SCOPED_LOCK(mutex);
		{
			mars::ptr< Picks > picks = new Picks();
			QTMorphs::MorphIterator mit = _pQTMorphs->morphs(rgBBox);

			while (mit)
			{
				mars::ptr< GeoMorph > pGM = static_cast< mars::ptr< GeoMorph > > (mit);

				if (find(pGM) != end() && (rgBBox ^ mit.region()))
				{
					picks->insert(pGM);
					erase(pGM);
				}
				++mit;
			}

			picks->computeBBox();
			return picks;
		}
	}

	bool PagedGeoWorld::MasterSet::empty()
	{
		GW_SCOPED_LOCK(mutex);
		{
			return MorphSet::empty();
		}
	}

	PagedGeoWorld::MasterSet::MasterSet( QTMorphs * pQTMorphs, MorphList::const_iterator iMorphs0, MorphList::const_iterator iMorphsN ) 
		: MorphSet(iMorphs0, iMorphsN), _pQTMorphs(pQTMorphs) {}

	void PagedGeoWorld::MasterSet::push( const mars::ptr< Picks > & pPicks )
	{
		GW_SCOPED_LOCK(mutex);
		{
			insert(pPicks->begin(), pPicks->end());
		}
	}

	GeoWorld::GeoWorld( DepositSynthesizer * pDepSynth, const unsigned short nLODLowest, const unsigned short nLODMid, FileSource * pFile )
		: _nLODLowest(nLODLowest), _nLODMid(nLODMid), _phmLowest(NULL), _nTileDim(pFile->getVirtualTileDim() >> (nLODLowest - nLODMid)), 
		_nLowWidth(pFile->getWorldWidth()), _nLowHeight(pFile->getWorldHeight()), _pDepSynth(pDepSynth), _pFile(pFile)
	{
		using namespace mars;

		assert(nLODLowest >= MIN_LOWEST_LOD);

		_phmLowest = new GeoHeightMap (
			static_cast< unsigned short > (pFile->getWorldWidth() >> (nLODLowest - nLODMid)),
			static_cast< unsigned short > (pFile->getWorldHeight() >> (nLODLowest - nLODMid)),
			mars::Wrap
		);
		_phmLowest->clear();
		assert(_pDepSynth->getStratum().getBase().gridtype == mars::Wrap);

		_plsbboxLockedMutableRegions = new RegionContainer(_phmLowest->width, _phmLowest->height);
		_pAllAgentMaps = new AllAgentMaps(_phmLowest->width, _phmLowest->height);
		assert((_nLowWidth % _nTileDim) == 0 && (_nLowHeight % _nTileDim) == 0);

		assert(
			_nLowWidth <= FileSource::MAX_AXIS_TILES * FileSource::MAX_TILE_DIM &&
			_nLowHeight <= FileSource::MAX_AXIS_TILES * FileSource::MAX_TILE_DIM
		);
	}

	GeoWorld::~GeoWorld()
	{
		delete _phmLowest;
		delete _plsbboxLockedMutableRegions;
		delete _pAllAgentMaps;
	}

	void GeoWorld::unlock( LockedGWSection * pSection )
	{
		using namespace std;
		using namespace mars;

		GW_SCOPED_LOCK(mutex);
		{
			if (!_plsbboxLockedMutableRegions->removeRegion(pSection->bbox))
				throw UnlockEx(pSection->bbox);

			pSection->processSyncOp();
			*pSection >> LockedGWSection::View(_phmLowest, &_pDepSynth->getStratum());
			delete pSection;
		}
	}

	void GeoWorld::unlock( const LockedGWSection * pSection ) const
	{
		GW_SCOPED_LOCK(mutex);
		{
			if (_plsbboxLockedMutableRegions->intersects(pSection->bbox))
				throw UnlockEx(pSection->bbox);

			delete pSection;
		}
	}

	geoworld::LockedGWSection * GeoWorld::lock( const mars::BBox< long > & bbox, const bool bReadOnly ) const
	{
		using namespace std;
		using namespace mars;

		assert(static_cast< unsigned short > (bbox.getWidth()) <= _phmLowest->width && static_cast< unsigned short > (bbox.getHeight()) <= _phmLowest->height);

		GW_SCOPED_LOCK(mutex);
		{
			const mars::BBox< long > 
				bboxActual = LockedGWSection::computeLockedSection(bbox, _nTileDim);

			if (_plsbboxLockedMutableRegions->intersects(bboxActual))
				return NULL;

			if (!bReadOnly)
				_plsbboxLockedMutableRegions->addRegion(bboxActual);

			LockedGWSection * pSection = new LockedGWSection(bboxActual, _pFile->getWorldDepth(), getLOD(), _pFile->geohost.count, _pFile->geohost.getLOD());

			*pSection << LockedGWSection::View(_phmLowest, &_pDepSynth->getStratum());
			return pSection;
		}
	}

	LockedGWSection * GeoWorld::lock(const bool bReadOnly) const
	{
		using namespace mars;

		GW_SCOPED_LOCK(mutex);
		{
			const BBox< long >
				bboxAll = BBox< long >(
					0, 
					0, 
					(getWidth() >> getLOD()) - 1, 
					(getHeight() >> getLOD()) - 1
				);

			if (_plsbboxLockedMutableRegions->intersects(bboxAll))
				return NULL;

			if (!bReadOnly)
				_plsbboxLockedMutableRegions->addRegion(bboxAll);

			LockedGWSection * pSection = new LockedGWSection(bboxAll, _pFile->getWorldDepth(), getLOD(), _pFile->geohost.count, _pFile->geohost.getLOD());
			*pSection << LockedGWSection::View(_phmLowest, &_pDepSynth->getStratum());

			return pSection;
		}
	}

	void GeoWorld::runMorphs( const MineralsDB & mindb, MorphList & rejects )
	{
		size_t c = 0;

		for (MorphList::iterator i = _morphs.begin(); i != _morphs.end(); ++i, ++c)
		{
			BoundedGeoMorph * pBGM = dynamic_cast< BoundedGeoMorph * > (i->pointer());	// TODO: Lazily going through bounded geomorphs

			if (pBGM != NULL)
			{
				mars::BBox< long > bbox = pBGM->getBBox();
				const mars::BBox< long > 
					bboxActual = LockedGWSection::computeLockedSection(bbox, _nTileDim);

				if (bboxActual.getWidth() <= _phmLowest->width && 
					bboxActual.getHeight() <= _phmLowest->height && 
					LockedGWSection::isValid(bbox, _nTileDim))
				{
					fire(&IGeoWorldListener::runningMorph, *i, c, _morphs.size());

					LockedGWSection * pLockedSection = lock(
						bbox.left, 
						bbox.top, 
						static_cast< unsigned short > (bbox.getWidth()), 
						static_cast< unsigned short > (bbox.getHeight())
					);

					runMorph(mindb, *i, pLockedSection);

					unlock(pLockedSection);
				} else
				{
					fire(&IGeoWorldListener::rejectedMorph, *i);
					rejects.push_back(*i);
				}
			} else
			{
				fire(&IGeoWorldListener::runningMorph, *i, c, _morphs.size());

				LockedGWSection * pLockedWorld = lock();
				runMorph(mindb, *i, pLockedWorld);
				unlock(pLockedWorld);
			}
		}

		LOGN << "Macro world DEM range: " << _phmLowest->range();
	}

	bool GeoWorld::addMorph( mars::ptr <GeoMorph> & gm )
	{
		if (checkGeoMorph(gm))
		{
			_morphs.push_back(gm);
			return true;
		}

		return false;
	}

	void GeoWorld::addMorphs( MorphList & morphs, MorphList & rejects )
	{
		for (MorphList::iterator i = morphs.begin(); i != morphs.end(); ++i)
		{
			if (!checkGeoMorph(*i))
			{
				rejects.push_back(*i);
				i = morphs.erase(i);

				if (i == morphs.end())
					break;
			}
		}		
		
		_morphs.insert(_morphs.end(), morphs.begin(), morphs.end());
	}

	bool GeoWorld::checkGeoMorph( const mars::ptr<GeoMorph> &gm )
	{
		const BoundedGeoMorph * pBGM = dynamic_cast< const BoundedGeoMorph * > (gm.pointer());

		if (pBGM != NULL)
		{
			const mars::BBox< long > bbox = pBGM->getBBox();

			if (bbox.getWidth() >= getMaxGeoMorphWidth() || bbox.getHeight() >= getMaxGeoMorphHeight())
				return false;
		}

		return true;
	}

	unsigned short GeoWorld::getMaxGeoMorphWidth() const
	{
		return static_cast< unsigned short > ((_nLowWidth >> getDeltaLOD()) - _nTileDim * 2 - 1);	// - 1 accounts for non-event redundant tile-blocking
	}

	unsigned short GeoWorld::getMaxGeoMorphHeight() const
	{
		return static_cast< unsigned short > ((_nLowHeight >> getDeltaLOD()) - _nTileDim * 2 - 1);  // - 1 accounts for non-event redundant tile-blocking
	}

	void GeoWorld::write(const float fCoarseness)
	{ 
		_pFile->applyDEM(_phmLowest, fCoarseness); 
		*_pDepSynth >> _pFile->geohost;
		*_pDepSynth >> _pFile->minerals;
	}

	void GeoWorld::runMorph( const MineralsDB & mindb, mars::ptr< GeoMorph > pGM, LockedGWSection * pLockedSection )
	{
		HeightMapGeoMorph * pCGM = mars::ptr_dynamic_cast< HeightMapGeoMorph > (pGM);

		if (pCGM != NULL)
		{
			pCGM->doMorph(*pLockedSection);
			LOGN << "Macro world DEM range after morph " << pGM->getTypeName() << ": " << pLockedSection->getElevationRange();
		}

		MineralSpawnGeoMorph * pMSGM = mars::ptr_dynamic_cast< MineralSpawnGeoMorph > (pGM);

		if (pMSGM != NULL)
		{
			pMSGM->doSpawnMinerals(_pDepSynth, _ctrMins);
		}
	}

	bool QTMorphs::addMorph( mars::ptr <GeoMorph> gm )
	{
		BoundedGeoMorph * pBGM = dynamic_cast< BoundedGeoMorph * > (gm.pointer());

		if (pBGM != NULL)
		{
			const mars::BBox< long > bbox = pBGM->getBBox();

			assert(!bbox.empty());
			add(gm, Region(bbox.getMinimum(), bbox.getMaximum()));

			return true;
		}

		return false;
	}

	bool filesource::InfoHeader::validate() const
	{
		return strncmp(IDENT, "PGWM", sizeof(IDENT)) == 0 && 
			endian == 1 &&
			version == 1 && 
			tdim > 0 && 
			tx > 0 && ty > 0 && depth > 0;
	}

	filesource::InfoHeader::InfoHeader() 
		: endian(1), version(1), tdim(0), tx(0), ty(0), depth(0)
	{
		memcpy(IDENT, "PGWM", sizeof(IDENT));
	}

	void MineralsContainer::flushTo( FileSource::MineralMap & minmap )
	{
		for (SphPackList::const_iterator i = _sphpacks.begin(); i != _sphpacks.end(); ++i)
			for (GeoSpherePack::const_iterator j = i->second->iterate(); j; ++j)
				minmap.add(i->first, *j);

		_sphpacks.clear();
	}

	bool RegionContainer::removeRegion( const mars::BBox< long > & bbox )
	{
		bool bFound = false;
		std::list< mars::BBox< long > > ::iterator i = begin();

		while (i != end())
			if (bbox == *i)
			{
				i = erase(i);
				bFound = true;
			} else
				++i;

			return bFound;
	}

	bool RegionContainer::intersects( const mars::BBox< long > & bbox ) const
	{
		using namespace mars;

		for (list< mars::BBox< long > > ::const_iterator i = begin(); i != end(); ++i)
			if (bbox.intersects(*i))
				return true;

		// Now perform more rigorous testing, checks bounding boxes that wrap around the world edges
		const BBox< long > 
			bbox_2r = bbox + vector2D< long > (static_cast< long > (_nWWidth), 0),
			bbox_2b = bbox + vector2D< long > (0, static_cast< long > (_nWHeight)),
			bbox_2br = bbox + vector2D< long > (static_cast< long > (_nWWidth), static_cast< long > (_nWHeight));
		
		for (list< mars::BBox< long > > ::const_iterator i = begin(); i != end(); ++i)
			if (bbox_2r.intersects(*i) || bbox_2b.intersects(*i) || bbox_2br.intersects(*i))
				return true;

		return false;
	}

	float AgentMap::sample( const GWCoords & coords ) const
	{
		using namespace mars;

		double d = 1;

		for (QuadTreeSpherePacks::Explicit< QuadTreeSpherePacks::Tag3D >::const_EntityPointIterator
			i = _qtField.contents(PointRegion< QuadTreeSpherePacks::Tag3D > (coords)); 
			i; ++i)
		{
			const double
				fDistSQ = static_cast< double > (MAGSQ(i.region().p3() - coords)),
				rSQ = static_cast < double > (mars::SQ(i.region().r));

			// TODO: Check
			d *= fDistSQ / rSQ;
		}

		return static_cast< float > (1.0 - d + base);
	}

	mars::vector3Df AgentMap::gradient( const GWCoords & coords ) const
	{
		using namespace mars;

		vector3Df v;
		GWCoords temp;

		for (QuadTreeSpherePacks::Explicit< QuadTreeSpherePacks::Tag3D >::const_EntityPointIterator 
			i = _qtField.contents(PointRegion< QuadTreeSpherePacks::Tag3D > (coords)); 
			i; ++i)
		{
			// TODO: Check
			return v += CAST< float > (
				i.region().p3() - coords
			) / static_cast< float > (i.region().r);
		}

		return v;
	}


	AgentMap & AllAgentMaps::operator[]( const FieldType enftType )
	{
		switch (enftType)
		{
		case FT_Density:
			return density;
		}

		throw PGWIllegal ("Invalid field type requested for agent map");
	}

	const AgentMap & AllAgentMaps::operator[]( const FieldType enftType ) const
	{
		switch (enftType)
		{
		case FT_Density:
			return density;
		}

		throw PGWIllegal ("Invalid field type requested for agent map");
	}

	AllAgentMaps::AllAgentMaps( const unsigned long nWidth, const unsigned long nHeight ) 
		: density(FT_Density, nWidth, nHeight, 1.0f)
	{}

}