#include "winderode.h"

#include <omp.h>

#include "geoworld.h"

namespace geoworld
{
	void WindErosionGM::init(SynthesisSession & manager, const IGeoWorldAccessor & accessor, const unsigned short nMaxWidth, const unsigned short nMaxHeight)
	{
		_psfDensity = accessor.queryField(FT_Density);
	}

	void WindErosionGM::doMorph( LockedGWSection & section ) const
	{
		using namespace mars;
		using namespace geoworld;

		LockedGWSection::View
			* pView = section.createView(section.bbox.left, section.bbox.top, section.bbox.right, section.bbox.bottom);

		ScalarField< RationalType > 
			hmf(pView->width, pView->height, Wrap),
			hmdf(pView->width, pView->height, Wrap);

		const RangeX< GeoHeightMap::Precision > rnghm = pView->getHeightMap() ->range();
		GaussianFn< RationalType > 
			gauss(1.0, 1.0, rnghm.delta * _frDecayLevel, 0);
		signed int i, j, k, g;
		VectorTag< RationalType >::V2 vgr;

		ptr< NoiseGen::NoiseMap > pnm = _ngen.compute(pView->width, pView->height);
		const NoiseGen::NoiseMap & nm = *pnm;
		const RangeX< NoiseGen::NoiseMap::Precision > rngnm = nm.range();
		const signed int len = static_cast< signed int > (pView->getHeightMap() ->length);	// It's ok, HeightMap.length is never more than the positive maximum of a signed-integer

		// Copy integral-precision field to rational-precision field
		hmf << *pView->getHeightMap();

		for (unsigned int c = 0; c < _iterations; ++c)
		{
			for (j = 0; j < hmf.height; ++j)
				for (i = 0; i < hmf.width; ++i)
				{
					const RationalType mh = hmf.mean5(i, j);
					// This ensures that pits and gullies are ignored
					if (mh - hmf(i, j) < 0)
					{
						const RationalType	// Get the noise map ratio at i,j as a ratio to the noise map range
							nr = static_cast< RationalType > ((nm(i,j) - rngnm.minimum) / rngnm.delta);

						// Obtain gradient into vgr, sampling distance of 5
						grad5(hmf, i, j, vgr); 

						const RationalType
							// Sample density half-way down at current spot
							fp = static_cast< RationalType > (1) / 
								static_cast< RationalType > (
									(*_psfDensity) (i, j, static_cast< GeoHeightMap::Precision > (hmf(i,j) * 0.5f))
								),	
							pgr = MAG(vgr),

							// Amount of material to remove
							rdh = pgr * fp * nr;

						// If the displacement is zero or the volume of material to move is zero then skip
						if (rdh > 0.0)
						{
							// Remove a smooth pile of material from the source site
							for (k = i - 2; k <= i + 2; ++k)
								for (g = j - 2; g <= j + 2; ++g)
									hmdf(k, g) -= gaussblur(k - i, g - j) * rdh;

							const RationalType
								rddh = rdh * gauss.f(pgr);	// TODO: Export gauss fn to pre-init tables state

							// Scale-up the displacement by a random amount, makes a very pleasing result
							vgr *= RANDf(mh) + static_cast< RationalType > (1);

							const signed int 
								x0 = i - RNDi(vgr.x), 
								y0 = j - RNDi(vgr.y);

							// Deposits "alluvium" into a smooth pile
							for (k = x0 - 2; k <= x0 + 2; ++k)
								for (g = y0 - 2; g <= y0 + 2; ++g)
									hmdf(k, g) += gaussblur(k - x0, g - y0) * rddh;
						}
					}
				}

			hmf += hmdf;
			hmdf.clear();
		}

		if (_bSmoothPost)
		{
			#pragma omp parallel for private(i, j)
			for (j = 0; j < hmf.height; ++j)
				for (i = 0; i < hmf.width; ++i)
					(*pView->getHeightMap())(i, j) = static_cast< GeoHeightMap::Precision > (
						hmf(i,j) + (
							static_cast< RationalType > (hmf.meanVN(i, j) - hmf(i, j))
						) 
						* (nm(i, j) - rngnm.minimum) / rngnm.delta
					);
		} else
			*pView->getHeightMap() << hmf;

		section.releaseView(pView);
	}

	WindErosionGM::~WindErosionGM()
	{
		// HACK: Use releaseField
		delete _psfDensity;
	}

	WindErosionGM::WindErosionGM( const unsigned int nIterations /*= 8*/, const float frDecayLevel /*= 0.1f*/, const bool bSmoothPost /*= true */ ) : _iterations(nIterations), _frDecayLevel(frDecayLevel), _ngen(0.3f, 0.5f), _bSmoothPost(bSmoothPost), _pAccessor(NULL), _psfDensity(NULL)
	{
		using namespace mars;
		GaussianFn< RationalType > gauss(
			static_cast< RationalType > (1), 
			static_cast< RationalType > (1) / sqrt (
				static_cast< RationalType > (2) * static_cast< RationalType > (PI)
			), 
			1, 
			0
		);	// TODO: Export to pre-init tables API

		for (signed int x = -2; x <= +2; ++x)
			for (signed int y = -2; y <= +2; ++y)
				gaussblur(x, y) = gauss.f(static_cast< RationalType > (x), static_cast< RationalType > (y));
	}

	GeoMorph * WindErosionGTFactory::createRandomInstance() const
	{ return new WindErosionGM(_mmnIterations.next(), _frDecayLevel, _bSmoothPost); }

	void WindErosionGTFactory::configure( IConfigGTFactory * pFactoryConfig, const IConfigGTFactory::Settings & settings )
	{
		const geoworld::IConfigSection & section = *pFactoryConfig->getSection();

		*section["iterations"] >> _mmnIterations;
		*section["decay"] >> _frDecayLevel;
		*section["smooth"] >> _bSmoothPost;
	}
}