// (C) 2002 Peter Horsfield; N-Dimensional Perlin Noise
// The perlin class defines it's methods here but could
// would normally be #included, or made into a non-template
// class. I have used this class to generate terrain and 
// ascii clouds.

// Template parameters:

// NOISE:	The source of the values to be interpolated between; 
//			I usually use the Mersenne Twister since it is fast
//			and very random. Must look like a PCTable with a seed
//			function.
		
// N:		The number of dimensions in the hypercube.

// M:		If you consider that the Perlin noise algorithm allows
//   an interpolation between two points, then the coordinates
//   of those points should be considered integral, and the 
//   distance between them and the point to be sampled, fractional.
//   Defining them in this way shows how we could store them
//   in a fixed-point arithmetic value: I<<M + F

//   The horsepower of which I am particularly proud is used lower 
//   down in the axiom::Perlin::operator[] function, and described
//   in more detail.

/*
	// Load precomputed vertex value into the stack
	stack.push(hypercube.pop());

	// This is key: determine which bits(dimensions) overflowed					
	looplogic = (dimIt ^ (dimIt+1)) & ~(dimIt+1);
	
	// Store next
	++dimIt;

	// Perform interpolations
	while(looplogic) stack.push(interpolate(stack.pop(),stack.pop()));
*/

namespace axiom {

template<class NOISE, int M, int N>
class Perlin {
private:
	// Possibly the noise sources may be precomputed, and hence may
	// be quite large. Avoid the possibility of them being on the
	// stack by dynamically constructing them.
	std::auto_ptr<NOISE> d[N];

	// Define limits using a handy template power class
	enum { WRAP = PWR2<M>::VALUE };
	enum { TRUNC = WRAP-1 };

	// Handy inner class used to compute the ease curve
	class EaseCompute { public: 
		inline float MINVAL() { return 0.0f + Epsilon; }
		inline float MAXVAL() { return 1.0f - Epsilon; }
		float operator() (unsigned int n) { 
			float x = (float)n / (float)(WRAP);
			return 3*x*x - 2*x*x*x; 
		}
	};

	// Precomputable table containing the Ease curve.
	PCTable<float, PWR2<16>::VALUE, -2.0f, EaseCompute> ease;

public:
	// The parameter indicates precomputation of the tables
	Perlin(bool pc = false) : ease(pc) { 
		for(int i = 0; i < N; ++i) d[i] = std::auto_ptr<NOISE>(new NOISE(pc));
	}

	// Allow for each dimension to be individually seeded and possibly recomputed
	void seed(unsigned int dm, unsigned int sd, bool rc = false) { 
		d[dm]->seed(sd);
		if(rc)d[dm]->precompute();
		else d[dm]->nulltable();
	}

	// Performs the interpolation in one dimension
	inline float interpolate(float flow, float fhi, unsigned int adj)
	{
		float offs = ease[adj];				
		return (1.0f - offs) * flow + offs * fhi;
	}

	// Return a dimension
	NOISE& dim(int n) { return *(d[n]); }

	// Compute a Perlin noise value at the given coordinates.
	float operator[](unsigned int coords[N]) {

		// This routine performs a stack based interpolation of the 
		// vertexes of the hypercube surrounding coords.
		// Each vertex of that hypercube corresponds to a random noise value
		// The perlin computation interpolates between the vertices to produce
		// a smoothed value at any point within that hypercube.
		// Our computation does this by using the results of two 1d interpolations
		// as input to a 2d interpolation, similarly two n-d interpolations as
		// input to an (n+1)-d interpolation.

		// Temporary results are stored on the following stack
		float		  stackbuf[N+1]; 
		float*		  stack = stackbuf - 1;
		// Precomputation buffers of the hypercube within which we'll interpolate
		float hyperCubeBuf[1<<N];
		float *hyperCube = hyperCubeBuf + (1<<N);

		// Variables used to control whether an interpolation must be performed
		unsigned int looplogic;
		// Our interpolation is based on the coordinate in the current dimension
		unsigned int *curdim;
		// And on the current dimension on which we are operating.
		// bit0 = dim-1, bit1 = dim-2, bit n-1 = dim n
		// When a bit carries, then an interpolation can be performed
		unsigned int  dimIt;  

		// Calculate the value of the vertex of the N dimensional cube
		// We implemented this by looking up in N 1d noise vectors and 
		// then averaging. This is the killer since we are doing mem 
		// lookups into 64kb arrays, so we're precomputing: (Note loop termination occurs when the unsigned variable underflows) 
		for((unsigned)dimIt = (1<<N) - 1; dimIt < (1<<N) ; --dimIt)
		{
			*--hyperCube = 0;
			curdim = coords;
			for(unsigned int i = 0; i < N; ++i)
			{
				*hyperCube += dim(i)[(*curdim>>M) + ((dimIt>>i)&1)];
				++curdim;
			}
			*hyperCube /= N;
		}
		
		// Examples where N=3
		// dimIt = 000b, endmask = 111b, hyperCube = hyperCubeBuf-1, stack = stackBuf-1
		unsigned int  endmask = (1<<N)-1;
		++dimIt;
		--hyperCube;

		// Begin processing
		do {
			// Load precomputed vertex value into the stack
			*++stack = *++hyperCube;					

			// This is key: determine which bits(dimensions) overflowed					
			looplogic = (dimIt ^ (dimIt+1)) & ~(dimIt+1);

			// Store new value
			++dimIt;

			// Reset our input coord coordstor
			curdim = coords;

			// Now loop through the rolled dimensions, and interpolate the points
			// i.e.  1-2			1-2				1--
			//		3...4->5   =>   3-4->5    =>	  --4
			while(looplogic)
			{
				// since we will combine the top two points on the cube into one
				--stack;
				// perform interpolation
				*stack = interpolate(*stack, *(stack+1),*curdim);
				// move to next test
				looplogic >>= 1;
				// increment which coordinate value to check
				++curdim;
			}
		}
		while(dimIt & endmask);
		// Loop exits when we overflow in the N'th dimension
		// At this point all interpolations are done and we return
		// the bottom of the stack (stack == stackbuf)
		return *stack;
	}

	// My architecture indicates these should exist for testing purposes
	inline float MINVAL() { return d[0]->MINVAL(); }
	inline float MAXVAL() { return d[0]->MAXVAL(); }
}; /*Perlin*/
}  /*axiom*/