//
// $Id: main.cxx,v 1.12 2002/07/27 06:21:53 keijiro Exp $
//
#include <iostream>
#include <vector>
#include <boost/array.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_01.hpp>
#include <tvmet/Vector.h>
#include <IL/il.h>

// pi
const double PI = 3.14159265358979323846;

// real number class
typedef double Real;

// 3-vector class
typedef tvmet::Vector<Real, 3> Vector;

// vector normalization
Vector normalize(const Vector &v)
{   
	return Vector(v / sqrt(tvmet::dot(v, v)));
}

//
// uniform random number generator (mersenne twister method)
//
class Random
{
	static boost::mt19937 rand_gen;
	static boost::uniform_01<boost::mt19937> uni;

	const Real min, max;

public:

	Random(Real min, Real max) : min(min), max(max) {}

	Random(Real max = 1.0) : min(0.0), max(max) {}

	Real operator()()
	{
		return (max - min) * uni() + min;
	}
};

boost::mt19937 Random::rand_gen(std::time(NULL));
boost::uniform_01<boost::mt19937> Random::uni(rand_gen);

//
// scattered vector class
//
class ScatterVector : public Vector
{
public:

	ScatterVector(const Vector &normal)
	{
		Vector tan1, tan2;

		if (fabs(normal[1]) < 0.9)
		{
			// tan1 = normalize(normal x (0,1,0))
			tan1 = normalize(Vector(-normal[2], 0.0, normal[0]));
			tan2 = normalize(tvmet::cross(tan1, normal));
		}
		else
		{
			// tan1 = normalize(normal x (1,0,0))
			tan1 = normalize(Vector(0.0, normal[2], -normal[1]));
			tan2 = normalize(tvmet::cross(tan1, normal));
		}

		Real vert = sqrt(Random()());
		Real hori = sqrt(1.0 - vert);

		Real phi = Random(-PI, PI)();
		Real x = hori * sin(phi);
		Real z = hori * cos(phi);

		(*this)[0] = normal[0] * vert + tan1[0] * x + tan2[0] * z;
		(*this)[1] = normal[1] * vert + tan1[1] * x + tan2[1] * z;
		(*this)[2] = normal[2] * vert + tan1[2] * x + tan2[2] * z;
	}
};

//
// funtion object for intersection of ray and object
//
// TODO: consider using static polymorphism
//
class RayIntersection
{
protected:

	Vector inter, normal;

public:

	RayIntersection() {}
	virtual ~RayIntersection() {}

	// get coordinate of the previous intersection
	const Vector& prevInterCoord() const
	{
		return inter;
	}

	// get normal vector of the previous intersection
	const Vector& prevInterNormal() const
	{
		return normal;
	}

	//
	// intersection test
	//
	// Returns a negative value if there is no intersection,
	// or parameter (t) at the intersection.
	//
	virtual Real test(const Vector &org, const Vector &vec) = NULL;
};

//
// derivation: ray-sphere intersection
//
class SphereRayIntersection : public RayIntersection
{
	const Vector center;
	const Real radius;

public:

	SphereRayIntersection(const Vector &center, Real radius) :
		center(center), radius(radius) {}

	Real test(const Vector &org, const Vector &vec)
	{
		Vector diff(org - center);

		Real b = tvmet::dot(vec, diff);
		Real b_sq = b * b;

		Real c = tvmet::dot(diff, diff) - radius * radius;

		if (b_sq < c)
		{
			// no intersection
			return -1.0;
		}

		Real t = -b - sqrt(b_sq - c);

		if (t < 0.0)
		{
			// no intersection (org is inside the sphere)
			return -1.0;
		}

		// get the intersection
		inter = org + vec * t;
		normal = normalize(Vector(inter - center));

		return t;
	}
};

//
// derivation: ray and half infinite wall intersection
//
template <int axis, int sign>
class HalfWallIntersection : public RayIntersection
{
	const Real pos, height;

public:

	HalfWallIntersection(Real pos, Real height) : pos(pos), height(height)
	{
		normal = 0;
		normal[axis] = sign;
	}

	Real test(const Vector &org, const Vector &vec)
	{
		Real t = pos - org[axis];

		if (t * sign > 0 || vec[axis] * sign > 0.0)
		{
			// no intersection
			return -1.0;
		}

		t /= vec[axis];
		inter = org + vec * t;

		if (inter[1] > height)
		{
			// no intersection
			return -1.0;
		}

		return t;
	}
};

//
// Scene class
//
class Scene
{
	// object class
	struct Object
	{
		RayIntersection *inter;
		const Vector albedo;

		Object(RayIntersection *inter, const Vector &albedo) :
			inter(inter), albedo(albedo) {}

		~Object()
		{
			delete inter;
		}
	};

	typedef std::vector<Object*> ObjectList;

	ObjectList objs;
	Object *inter;

public:

	Scene() : inter(NULL)
	{
		// left sphere
		objs.push_back(new Object(
							new SphereRayIntersection(Vector(-0.5, 0.5, 0.0), 0.5),
							Vector(0.9, 0.9, 0.9)));

		// right sphere
		objs.push_back(new Object(
							new SphereRayIntersection(Vector(0.5, 0.5, 0.0), 0.5),
							Vector(0.9, 0.9, 0.9)));

		// top sphere
		objs.push_back(new Object(
							new SphereRayIntersection(Vector(0.0, 1.37, 0.0), 0.5),
							Vector(0.9, 0.9, 0.9)));

		// ground
		objs.push_back(new Object(
							new HalfWallIntersection<1, 1>(0.0, 1.0),
							Vector(0.75, 0.75, 0.75)));

		// left wall
		objs.push_back(new Object(
							new HalfWallIntersection<0, 1>(-1.2, 2.0),
							Vector(0.5, 0.5, 1.0)));

		// right wall
		objs.push_back(new Object(
							new HalfWallIntersection<0, -1>(1.2, 2.0),
							Vector(0.5, 1.0, 0.5)));

		// behind wall
		objs.push_back(new Object(
							new HalfWallIntersection<2, -1>(1.2, 2.0),
							Vector(0.75, 0.75, 0.75)));
	}

	~Scene()
	{
		for (ObjectList::iterator i = objs.begin(); i != objs.end(); i++)
		{
			delete *i;
		}
	}

	// returns the coordinate of the previous intersection
	const Vector& coord() const
	{
		return inter->inter->prevInterCoord();
	}

	// returns the normal vector of the previous intersection
	const Vector& normal() const
	{
		return inter->inter->prevInterNormal();
	}

	// return the albedo of the previous intersection
	const Vector& albedo() const
	{
		return inter->albedo;
	}

	// shoot a ray into the scene
	bool shoot(const Vector &org, const Vector &vec)
	{
		const float INF = 1e+10;
		Real dist = INF;

		for (ObjectList::iterator i = objs.begin(); i != objs.end(); i++)
		{
			Real t = (*i)->inter->test(org, vec);

			if (t > 0.0 && t < dist)
			{
				inter = *i;
				dist = t;
			}
		}

		return dist < INF;
	}
};

//
// renderer class
//
class Renderer
{
	Scene scene;
	Vector eye;

public:

	Renderer() : eye(0.0, 1.0, -4.0) {}

	Vector render(const Vector &ray)
	{
		Vector org(eye);
		Vector vec(ray);
		Vector rad(1.0, 1.0, 1.0);

		while (true)
		{
			// shoot a ray into the scene
			if (!scene.shoot(org, vec))
			{
				return rad;
			}

			// random shoot
			org = scene.coord();
			vec = ScatterVector(scene.normal());

			if (tvmet::sum(rad) > 1.5)
			{
				// hight importance: normal reflection
				rad *= scene.albedo();
			}
			else
			{
				// low importance: russian roulette
				if (Random(3.0)() > tvmet::sum(scene.albedo()))
				{
					return Vector(0.0);
				}
				else
				{
					rad = rad * scene.albedo() * 3.0 / tvmet::sum(scene.albedo());
				}
			}
		}
	}
};

//
// screen image class
//
template <int width, int height>
class Image
{
	const char *name;
	unsigned char raw[height * width * 3];

public:

	typedef const char* Arg;

	Image(Arg name) : name(name) {}

	// set a pixel
	void set(int x, int y, const Vector &c)
	{
		unsigned char *ptr = raw + 3 * (y * width + x);
		ptr[0] = (unsigned char)(c[0] * 255.0);
		ptr[1] = (unsigned char)(c[1] * 255.0);
		ptr[2] = (unsigned char)(c[2] * 255.0);
	}

	// save the image into a file
	void save()
	{
		ILuint id;

		ilGenImages(1, &id);
		ilBindImage(id);

		ilTexImage(width, height, 1, 3, IL_RGB, IL_UNSIGNED_BYTE, raw);
		ilEnable(IL_FILE_OVERWRITE);
		ilSaveImage(const_cast<char*>(name));

		ilDeleteImages(1, &id);
	}
};

//
// screen scanner class
//
template <typename Renderer, int width, int height, int itr>
class Scanner
{
	typedef Image<width, height> Image;

public:

	Real fov, dist;

	Scanner(Real fov, Real dist) : fov(fov), dist(dist) {}

	void scan(typename Image::Arg img_arg)
	{
		Renderer rend;
		Image img(img_arg);

		// screen width / height
		Real scr_w = tan(fov) * dist;
		Real scr_h = scr_w * height / width;

		// dot width / height
		const Real dot_w = scr_w / (width - 1.0);
		const Real dot_h = scr_h / (height - 1.0);

		// jitter generator
		Random dx(-0.5 * dot_w, 0.5 * dot_w);
		Random dy(-0.5 * dot_h, 0.5 * dot_h);

		for (int row = 0; row < height; row++)
		{
			Real y = scr_h * (row / (height - 1.0) - 0.5);

			for (int col = 0; col < width; col++)
			{
				Real x = scr_w * (col / (width - 1.0) - 0.5);
				Vector color(0.0);

				for (int i = 0; i < itr; i++)
				{
					color += rend.render(normalize(Vector(x + dx(), y + dy(), dist)));
				}
				color /= itr;

				img.set(col, row, color);
			}

			std::cout << "*" << std::flush;
		}

		img.save();

		std::cout << " done." << std::endl;
	}
};

//
// main procedure
//
int main()
{
	ilInit();

	Scanner<Renderer, 160, 120, 200>(PI / 4, 1.0).scan("out.png");

	return 0;
}