/*
 * Copyright (c) 2010 Dale Thomas http://www.inaneblabbering.com/dalethomas/
 *
 * Permission to use, copy, modify, distribute and sell this software
 * and its documentation for any purpose is hereby granted without fee,
 * provided that the above copyright notice appear in all copies.
 * Dale Thomas makes no representations about the suitability 
 * of this software for any purpose.  
 * It is provided "as is" without express or implied warranty.
 */

//#define SHOW_SAMPLING_GRID

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <vector>
using namespace std;

const int WIDTH = 80*8;
const int HEIGHT = 80*6;
const int WIDTH_PADDED = WIDTH+1;
const int HEIGHT_PADDED = HEIGHT+1;
const int MAX_DEPTH = 2;
const int NUM_SAMPLES = 200;
const int STEP = 8;
const float SCALE = 15;
const float SIGMA = 0.4;
const float INV_ASPECT = float(HEIGHT)/WIDTH;
const int NUM_FRAMES = 1000;


#define MYEPSILON 1e-5
#define RAND (float(rand())/RAND_MAX)
#define RRAND (RAND*2-1)

// Don't ask why it uses 3D vectors. IT JUST DOES, ALRIGHT!!!
struct Vec
{
	double x, y, z; // position, also color (r,g,b)
	Vec(double x_=0, double y_=0, double z_=0){ x=x_; y=y_; z=z_; }
	Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); }
	Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); }
	Vec operator*(double b) const { return Vec(x*b,y*b,z*b); }
	void operator*=(double b) { x*=b; y*=b; z*=b; }
	void operator+=(const Vec &b) { x+=b.x; y+=b.y; z+=b.z; }
	void operator*=(const Vec &b) { x*=b.x; y*=b.y; z*=b.z; }
	Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); }
	Vec& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); }
	double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross:
	Vec operator%(Vec&b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
};


vector<Vec> colours;
vector<unsigned int> flags;

struct Sphere
{
	Vec _pos, _col;
	float _rad, _rad2, _rfl;
	int _id;
	Sphere(Vec _P,float _R,Vec _C,float _Rfl,int _ID)
		: _pos(_P)
		, _rad(_R)
		, _rad2(_R*_R)
		, _col(_C)
		, _rfl(_Rfl)
		, _id(_ID)
	{ }
	bool intersect(Vec _Org,Vec _Dir,float& _Dist)
	{
		const Vec cam_to_sphere = _pos-_Org;
//		if(cam_to_sphere.dot(cam_to_sphere)<_rad2) { _Dist = 0.f; return true; }
		const float b = cam_to_sphere.dot(_Dir);
		float det = b*b-cam_to_sphere.dot(cam_to_sphere)+_rad2, t1, t2;
		if( det<MYEPSILON ) return false;
		det = sqrtf(det);
		if((t1=b-det)>_Dist) return false;
		if(t1>MYEPSILON) { _Dist = t1; return true; }
		if((t2=b+det)>_Dist) return false;
		if(t2>MYEPSILON) { _Dist = t2; return true; }
		return false;
	}
};

vector<Sphere> spheres;

inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }

void setup()
{
	colours.resize(WIDTH_PADDED*HEIGHT_PADDED);
	flags.resize(WIDTH_PADDED*HEIGHT_PADDED,0);

	spheres.push_back(Sphere(Vec( 0, 0, 0),1.5,Vec(0.2, 0.0, 0.0),0, 7));
	spheres.push_back(Sphere(Vec( 0, 0, 0),1.5,Vec(0.0, 0.2, 0.0),0, 8));
	spheres.push_back(Sphere(Vec( 0, 0, 0),1.5,Vec(0.0, 0.0, 0.2),0, 9));
	spheres.push_back(Sphere(Vec( 0, 0, 0),1.5,Vec(0.2, 0.2, 0.0),0,10));
	spheres.push_back(Sphere(Vec(-7, 0, 0),2.0,Vec(1.0, 0.0, 0.0),0, 1));
	spheres.push_back(Sphere(Vec( 0, 0, 0),2.0,Vec(1.0, 1.0, 1.0),0, 2));
	spheres.push_back(Sphere(Vec( 7, 0, 0),2.0,Vec(0.0, 1.0, 0.0),0, 3));
	spheres.push_back(Sphere(Vec(-7, 0, 0),7.0,Vec(0.0, 0.0, 0.0),0, 4));
	spheres.push_back(Sphere(Vec( 0, 0, 0),7.0,Vec(0.0, 0.0, 0.0),0, 5));
	spheres.push_back(Sphere(Vec( 7, 0, 0),7.0,Vec(0.0, 0.0, 0.0),0, 6));
}

Vec raytrace(int _Depth,Vec _Org,Vec _Dir,int _Index)
{
	Vec col(0,0,0);
	if(_Depth++>MAX_DEPTH) return col;
	Sphere* hit_obj = 0;
	float min_dist = 1e10;
	for(unsigned int i=0; i<spheres.size(); i++)
	{
		float temp_min_dist = 1e10;
		if( spheres[i].intersect(_Org,_Dir,temp_min_dist) )
		{
			if( min_dist>temp_min_dist )
			{
                min_dist = temp_min_dist;
				hit_obj = &(spheres[i]);
			}
		}
	}
	
	if( hit_obj )
	{
		if(_Depth==1) flags[_Index] += hit_obj->_id;
		col += hit_obj->_col*exp(-min_dist*SIGMA);
		if(hit_obj->_rfl>MYEPSILON)
		{

			Vec hit_pos = _Org+_Dir*min_dist;
			Vec hit_nrm = hit_pos-hit_obj->_pos; hit_nrm.norm();
			if( hit_nrm.dot(_Dir)>0.f ) hit_nrm *= -1.f;
			Vec rfl_dir = _Dir-hit_nrm*(_Dir.dot(hit_nrm)*2);
			col += raytrace(_Depth,hit_pos,rfl_dir,_Index)*hit_obj->_rfl;
		}
	}

	return col;
}

void paintPixel(int _Index)
{
	Vec col(0,0,0);
	for(int sample=0; sample<NUM_SAMPLES; sample++)
	{
		const float fx = (RAND+(_Index%WIDTH_PADDED))/WIDTH;
		const float fy = (RAND+(_Index/WIDTH_PADDED))/HEIGHT;
		const float f = M_PI*2*float(sample)/NUM_SAMPLES;
		Vec dir(sinf(f),cosf(f),0);
		col += raytrace(0,Vec((fx*2.0-1.0),(fy*2.0-1.0)*INV_ASPECT,0)*SCALE,dir,_Index);
	}
	colours[_Index] = col*(1.f/NUM_SAMPLES);	
}

void sampleGrid(int _Index,int _Size)
{
	int i[4];
	i[0] = _Index;
	i[1] = i[0]+_Size;
	i[2] = i[1]+_Size*WIDTH_PADDED;
	i[3] = i[2]-_Size;

	unsigned int flag = flags[i[0]];
	if( flag==flags[i[1]] && flag==flags[i[2]] && flag==flags[i[3]] )
	{ // simply interpolate colour for this patch
		for(int y=0; y<_Size; ++y)
		{
			float fy = float(y)/(_Size-1), ffy = 1.0-fy;
			for(int x=0; x<_Size; ++x)
			{
				int index = _Index+x+y*WIDTH_PADDED;
				if( !flags[index] )
				{
					flags[index] = flag;
#ifdef SHOW_SAMPLING_GRID
					colours[index] = Vec(0,1,0);
#else
					float fx = float(x)/(_Size-1), ffx = 1.0-fx;

					colours[index] =
						(colours[i[0]]*ffx + colours[i[1]]*fx)*ffy +
						(colours[i[3]]*ffx + colours[i[2]]*fx)*fy;
#endif
				}
			}
		}
	}
	else
	{
		_Size >>= 1;
		int ii[5];

		ii[0] = _Index+_Size;
		ii[2] = ii[0]+_Size*WIDTH_PADDED;
		ii[1] = ii[2]-_Size;
		ii[3] = ii[2]+_Size;
		ii[4] = ii[2]+_Size*WIDTH_PADDED;

		if( !flags[ii[0]] ) paintPixel(ii[0]);
		if( !flags[ii[1]] ) paintPixel(ii[1]);
		if( !flags[ii[2]] ) paintPixel(ii[2]);
		if( !flags[ii[3]] ) paintPixel(ii[3]);
		if( !flags[ii[4]] ) paintPixel(ii[4]);

		if( _Size>1 )
		{ // recursively subdivide
			sampleGrid( i[0],_Size);
			sampleGrid(ii[0],_Size);
			sampleGrid(ii[1],_Size);
			sampleGrid(ii[2],_Size);
		}
	}
}

int main()
{
	setup();


	Vec freq[7],phase[7],pos[0];
	float dist[7];
    for(int s=0; s<7; s++)
    {
		freq[s] = Vec(RRAND,RRAND,0)*0.1f;
		phase[s] = Vec(RRAND,RRAND,0)*M_PI;
		pos[s] = spheres[s]._pos;

		dist[s] = (s<4) ? 3 : 1;
    }

	for(int frame=0; frame<NUM_FRAMES; frame++)
    {
  		flags.assign(WIDTH_PADDED*HEIGHT_PADDED,0);
  		
        // *** animate some objects ***
		{
			const float f = float(frame+1000);
			for(int s=0; s<7; s++)
			{
				spheres[s]._pos = pos[s]+Vec(sinf(f*freq[s].x+phase[s].x),cosf(f*freq[s].y+phase[s].y),0)*dist[s];
			}
		}
            
        
        // *** paint a coarse grid ***
    	for( int y=0; y<HEIGHT_PADDED; y+=STEP )
    	{
    		for( int x=0; x<WIDTH_PADDED; x+=STEP )
    		{
    			paintPixel(x+y*WIDTH_PADDED);
    
    		}
    	}
    
        // *** fill in gaps ***
    	for( int y=0; y<HEIGHT; y+=STEP )
    	{
    		for( int x=0; x<WIDTH; x+=STEP )
    		{
    			sampleGrid(x+y*WIDTH_PADDED,STEP);
    		}
    	}

        {    
            char buff[256];
            sprintf(buff,"image%.5d.ppm",frame);
        	FILE *f = fopen(buff, "w"); // Write image to PPM file.
        	fprintf(f, "P3\n%d %d\n%d\n", WIDTH, HEIGHT, 255);
        	for( int y=0; y<HEIGHT; y++ )
        	{
        		int i = y*WIDTH_PADDED;
        		for( int x=0; x<WIDTH; x++,i++ )
        		{
        			fprintf(f,"%d %d %d ", toInt(colours[i].x), toInt(colours[i].y), toInt(colours[i].z));
        		}
        	}
         }
		cout << "finished frame " << frame+1 << "/" << NUM_FRAMES << endl;
     }
	getchar();
	return 0;
}