#ifndef FRACTALMOUNTAINS_HXX
#define FRACTALMOUNTAINS_HXX

#include "Object.hxx"
#include "Vec3f.hxx"
#include "Noise.hxx"

/* This class creates mountains on a given plane using fractal geometry.
** For details of algorithm see "Texturing & Modeling: A Procedural Approach" page 500.
*/
class FractalMountains 
{

private:
  
    Vec3f** points;   // stores computed points
    int resX;
    int resY;

    // using wrap-around boundary conditions 
    Vec3f GetPoint(int x, int y)
    {
        if(x == -1)
	    x = resX;
	else if(x == (resX + 1))
	   x = 0;

	if(y == -1)
	   y = resY;
	else if(y == (resY + 1))
	    y = 0;

	return points[x][y];
    }

public:

    FractalMountains(Object* obj, Vec3f p1, Vec3f p2, Vec3f p3, int resX, int resY, float H, float lacunarity, int octaves, float height)
      : resX(resX),resY(resY)
    {
        // initialize points matrix
        points = new Vec3f* [resX+1];
	for(int i = 0; i <= resX; ++i)
	{
	    points[i] = new Vec3f[resY+1];
	}
	
	// compute the two direction vectors, spanning the mountains plane
	// p1 is origin, p2 lies on x-axis, p3 lies on y-axis
	Vec3f dirX = p2 - p1;
	Vec3f dirY = p3 - p1;

	// compute normal of plan (growth direction of mountains)
	Vec3f normal = Cross(dirX,dirY);
	Normalize(normal);
      
	// compute step sizes
	float stepX = Length(dirX) / static_cast<float>(resX);
	float stepY = Length(dirY) / static_cast<float>(resY);
	
	Normalize(dirX);
	Normalize(dirY);
	
	// precompute weights array
	float exponent_array[octaves+1];
	float frequency = 1.0;
	for(int i = 0; i <= octaves; ++i)
	{
	    exponent_array[i] = powf(frequency, -H);
	    frequency *= lacunarity;
	}	
	
	// compute value for each point according algorithm in book
	for(int x = 0; x <= resX; ++x)
	{  
	    for (int y = 0; y <= resY; ++y)
	    {
		points[x][y] = p1 + dirX * stepX * static_cast<float>(x) + dirY * stepY * static_cast<float>(y);
		Vec3f point = Vec3f(x,y,0) / resX * 5.0f;
	      
		float value = PerlinNoise3D_new(point.x(), point.y(), point.z());
		
		point *= lacunarity;
	      
		for(int i = 1; i < octaves; ++i)
		{
		    value += PerlinNoise3D_new(point.x(), point.y(), point.z()) * exponent_array[i] * value;
		    point *= lacunarity; 
		} 

		points[x][y] += normal * height * value; 
	    }
	}
	
	// compute a normal for each point
	Vec3f** normals = new Vec3f* [resX+1];
	for(int i = 0; i <= resX; ++i)
	    normals[i] = new Vec3f[resY+1];
	
	for(int x = 0; x <= resX; ++x )
	{
	    for(int y = 0; y <= resY; ++y)
	    {
	        Vec3f n1 = Cross(GetPoint(x,y+1) - GetPoint(x,y), GetPoint(x-1,y+1) - GetPoint(x,y));
		Normalize(n1);
		
		Vec3f n2 = Cross(GetPoint(x+1,y) - GetPoint(x,y), GetPoint(x,y+1) - GetPoint(x,y));
		Normalize(n2);
		
		Vec3f n3 = Cross(GetPoint(x+1,y-1) - GetPoint(x,y), GetPoint(x+1,y) - GetPoint(x,y));
		Normalize(n3);
		
		Vec3f n4 = Cross(GetPoint(x,y-1) - GetPoint(x,y), GetPoint(x+1,y-1) - GetPoint(x,y));
		Normalize(n4);
		
		Vec3f n5 = Cross(GetPoint(x-1,y) - GetPoint(x,y), GetPoint(x,y-1) - GetPoint(x,y));
		Normalize(n5);
		
		Vec3f n6 = Cross(GetPoint(x-1,y+1) - GetPoint(x,y), GetPoint(x-1,y) - GetPoint(x,y));
		Normalize(n6);
		
		normals[x][y] = n1 + n2 + n3 + n4 + n5 + n6;
		Normalize(normals[x][y]);
	    }
	}

	// add triangles to object
	Object::SmoothTriangleFactory factory;
	
	for(int x = 0; x < resX; ++x )
	{
	    for(int y = 0; y < resY; ++y)
	    {
		obj->Add(factory.Create(points[x][y], points[x+1][y], points[x][y+1], normals[x][y], normals[x+1][y], normals[x][y+1], 0, 0, 0));
		obj->Add(factory.Create(points[x+1][y], points[x+1][y+1], points[x][y+1], normals[x+1][y], normals[x+1][y+1], normals[x][y+1], 0, 0, 0));
	    }
	}
	
	// delete points matrix
	for(int i = 0; i <= resX; ++i)
	    delete[] points[i];
      
	delete[] points;

	// delete normals matrix
	for(int i = 0; i <= resX; ++i)
	    delete[] normals[i];

	delete[] normals;
    };
};

#endif