#include "pbrt.h"
#include "texture.h"
#include "paramset.h"

#include <iostream>
#include <cmath>
#include <list>
#include "cheap_rand.c"


static int Poisson_count[256]=
{4,3,1,1,1,2,4,2,2,2,5,1,0,2,1,2,2,0,4,3,2,1,2,1,3,2,2,4,2,2,5,1,2,3,2,2,2,2,2,3,
 2,4,2,5,3,2,2,2,5,3,3,5,2,1,3,3,4,4,2,3,0,4,2,2,2,1,3,2,2,2,3,3,3,1,2,0,2,1,1,2,
 2,2,2,5,3,2,3,2,3,2,2,1,0,2,1,1,2,1,2,2,1,3,4,2,2,2,5,4,2,4,2,2,5,4,3,2,2,5,4,3,
 3,3,5,2,2,2,2,2,3,1,1,4,2,1,3,3,4,3,2,4,3,3,3,4,5,1,4,2,4,3,1,2,3,5,3,2,1,3,1,3,
 3,3,2,3,1,5,5,4,2,2,4,1,3,4,1,5,3,3,5,3,4,3,2,2,1,1,1,1,1,2,4,5,4,5,4,2,1,5,1,1,
 2,3,3,3,2,5,2,3,3,2,0,2,1,1,4,2,1,3,2,1,2,2,3,2,5,5,3,4,5,5,2,4,4,5,3,2,2,2,1,4,
 2,3,3,4,2,5,4,2,4,2,2,2,4,5,3,2};


class DistMetric {
public:
  virtual float operator()(const Point &p1, const Point &p2) = 0;
  virtual ~DistMetric() {}
};

class EuclideanDist: public DistMetric {
public:
  float operator()(const Point &p1, const Point &p2) {
	return Distance(p1, p2);
  }

};

class ManhattanDist: public DistMetric {
public:
  float operator()(const Point &p1, const Point &p2) {
	 return fabs(p1.x - p2.x) + fabs(p1.y - p2.y) + fabs(p1.z - p2.z);
  }
};


class MinkovskyDist: public DistMetric {
public:
  MinkovskyDist(float _e): e(_e) {}; 
  float operator()(const Point &p1, const Point &p2) {
	return pow( pow(fabs(p1.x - p2.x), e) + 
				pow(fabs(p1.y - p2.y), e) +
				pow(fabs(p1.z - p2.z), e), 1.0/e);

  }
private:
  float e;
};

class ChebyshevDist: public DistMetric {
public:
  float operator()(const Point &p1, const Point &p2) {
  return max(fabs(p1.x - p2.x), max(fabs(p1.y - p2.y), fabs(p1.z - p2.z)));
  }

};



void insert_sorted_point(std::list<float> &l, const Point &p, const Point &center, DistMetric *metric) {

	float Pdist = (*metric)(p, center);
	std::list<float>::iterator pos;

	for(pos = l.begin(); pos != l.end(); pos++) {
		if(	*pos >= Pdist ) {
			break;
		}
	}

	l.insert(pos, Pdist);

	if (l.size() > 4 ) {
		l.pop_back();
	}
}


void add_cube(long cube_x, long cube_y, long cube_z, 
			  std::list<float> &l, const Point &at,
			  DistMetric *metric) {

	// get seed for cell
	unsigned long seed = (541*cube_x+79*cube_y +31*cube_z) % int(pow(2, 23));
	cheap_srand(seed);
	
	// find number of points in this cube
	int cell_count = Poisson_count[int(floor(cheap_rand() * 255))];
	if (cell_count == 0) {
		cell_count = 1;
	}
	
	for ( int i = 0; i < cell_count; i++) {
		Point p(cube_x + cheap_rand(), 
						cube_y + cheap_rand(), 
						cube_z + cheap_rand());

		insert_sorted_point(l, p, at, metric);
	}
}


float worley(Point at, const std::vector<float> &combination, DistMetric *metric) {
	std::list<float> point_list;

	for ( int off_x = -1; off_x < 2; off_x++ ) {
	  for ( int off_y = -1; off_y < 2; off_y++ ) {
		for ( int off_z = -1; off_z < 2; off_z++ ) {
		  long cube_x, cube_y, cube_z;
		  cube_x = long(floor(at.x) + off_x);
		  cube_y = long(floor(at.y) + off_y);
		  cube_z = long(floor(at.z) + off_z);
		  
		  add_cube(cube_x, cube_y, cube_z, point_list, at, metric);
		  
		}
	  }
	}

	float ret = 0;
	
	std::list<float>::iterator f = point_list.begin();
	for (unsigned int i = 1; i < combination.size(); i++) {
	  ret += combination[i-1] * sqrt(*f);
	  f++;
	}

	return ret;
}


// WorleyTexture Declarations
template <class T> class WorleyTexture : public Texture<T> {
public:
	// WorleyTexture Public Methods
	~WorleyTexture() {
		delete mapping;
	}
	WorleyTexture(vector<float> c, float t_scale, 
				  std::string metric_name, TextureMapping3D *map) {
		mapping = map;

		if( metric_name == "euclidean") {
		  metric = new EuclideanDist;
		} else if (metric_name == "manhattan") {
		  metric = new ManhattanDist;
		} else if (metric_name == "chebyshev") {
		  metric = new ChebyshevDist;
		} else if (std::string::npos != metric_name.find("minkovsky")) {
		  metric_name.erase(metric_name.find("minkovsky"), 9);
		  metric = new MinkovskyDist(atof(metric_name.c_str()));
		} else {
		  metric = new EuclideanDist;
		}

		scale = t_scale;
		combination = c;
	}
	T Evaluate(const DifferentialGeometry &dg) const {
		Vector dpdx, dpdy;
		Point P = mapping->Map(dg, &dpdx, &dpdy);

		return worley(P*scale, combination, metric);
	}
private:
  // WorleyTexture Private Data
  TextureMapping3D *mapping;
  vector<float> combination;
  float scale;
  DistMetric* metric;
};

// WorleyTexture Method Definitions
extern "C" DLLEXPORT Texture<float> * CreateFloatTexture(const Transform &tex2world,
		const TextureParams &tp) {
	// Initialize 3D texture mapping _map_ from _tp_
	TextureMapping3D *map = new IdentityMapping3D(tex2world);

	std::vector<float> c;

	c.push_back(tp.FindFloat("F1", 1));
	c.push_back(tp.FindFloat("F2", 0));
	c.push_back(tp.FindFloat("F3", 0));
	c.push_back(tp.FindFloat("F4", 0));
	
	return new WorleyTexture<float>(c, tp.FindFloat("texscale", 1), 
									tp.FindString("metric"), map);
}

extern "C" DLLEXPORT Texture<Spectrum> * CreateSpectrumTexture(const Transform &tex2world,
		const TextureParams &tp) {
	// Initialize 3D texture mapping _map_ from _tp_
	TextureMapping3D *map = new IdentityMapping3D(tex2world);

	std::vector<float> c;

	c.push_back(tp.FindFloat("F1", 1));
	c.push_back(tp.FindFloat("F2", 0));
	c.push_back(tp.FindFloat("F3", 0));
	c.push_back(tp.FindFloat("F4", 0));

	return new WorleyTexture<Spectrum>(c, tp.FindFloat("texscale", 1), 
									   tp.FindString("metric"), map);
}