#pragma once
#include "includes.h"
#include "bvh.h"
#include "file_helpers.h"
#include "points.h"
#include "config.h"
#include "texture.h"

#ifdef SPEED
//#include "poly_math.h"
inline double fast_pow(double a, double b) {
	return pow((float)a,(float)b);
}
#endif

double mod(double a, double b) {
	double result = fmod(a,b);
	if (result<0.0) result += b;
	return result;
}

struct Tile {
	int x,y,w,h;
};

struct ray_info {
	double distance;
	bool empty;
	bool valid;
};

void hsv_to_rgb(double hsv[3], double rgb[3]) {
	double chroma = hsv[2] * hsv[1];
	double hprime = hsv[0] / 60.0;
	double x = chroma  *  ( 1.0 - fabs(fmod(hprime,2.0)-1.0) );
	if      ( hprime>=0.0 && hprime<1.0 ) { rgb[0] += chroma; rgb[1] +=      x; rgb[2] +=      0; }
	else if ( hprime>=1.0 && hprime<2.0 ) { rgb[0] +=      x; rgb[1] += chroma; rgb[2] +=      0; }
	else if ( hprime>=2.0 && hprime<3.0 ) { rgb[0] +=      0; rgb[1] += chroma; rgb[2] +=      x; }
	else if ( hprime>=3.0 && hprime<4.0 ) { rgb[0] +=      0; rgb[1] +=      x; rgb[2] += chroma; }
	else if ( hprime>=4.0 && hprime<5.0 ) { rgb[0] +=      x; rgb[1] +=      0; rgb[2] += chroma; }
	else if ( hprime>=5.0 && hprime<6.0 ) { rgb[0] += chroma; rgb[1] +=      0; rgb[2] +=      x; }
	else {
		throw glLibError("Something went wrong in converting from HSV!");
	}
}
void add_color_heatmap(double color[3], double intensity) {
	intensity = glLibMathMisc::clamp(intensity,0.0,1.0);
	double hsv[3] = {(1.0-intensity)*240.0,1.0,1.0};
	hsv_to_rgb(hsv,color);
}
void map_vec_to_rgb(double vec[3], double color[3]) {
	#ifdef OUTPUT_ABS
	color[0] = fabs(vec[0]);
	color[1] = fabs(vec[1]);
	color[2] = fabs(vec[2]);
	#else
	color[0] = (vec[0]+1.0)*0.5;
	color[1] = (vec[1]+1.0)*0.5;
	color[2] = (vec[2]+1.0)*0.5;
	#endif
}

double fresnel_rs(double n1, double n2, double theta_i, double theta_t) {
	double a = n1*cos(theta_i), b = n2*cos(theta_t);
	double term = (a-b) / (a+b);
	return term * term;
}
double fresnel_rp(double n1, double n2, double theta_i, double theta_t) {
	double a = n1*cos(theta_t), b = n2*cos(theta_i);
	double term = (a-b) / (a+b);
	return term * term;
}

template <typename type_float> void get_disk(type_float point[2], type_float radius=1) {
	//http://mathworld.wolfram.com/DiskPointPicking.html
	type_float rand_angle = glLibMathRandom::get_uniform(0.0,2.0*glLibMathMisc::pi);
	double rand_radius_sqrt = sqrt(glLibMathRandom::get_uniform(0.0,radius*radius));
	point[0] = rand_radius_sqrt * cos(rand_angle);
	point[1] = rand_radius_sqrt * sin(rand_angle);
}
/*template <typename type_float> void get_hemisphere() {
	glLibMathRandom::get_sphere(type_float point[3]); //TODO;
}*/
template <typename type_float> void get_hemisphere_cosine(type_float point[3], type_float radius=1) {
	//Global Illumination Compendium, pg. 19
	type_float rand_angle = glLibMathRandom::get_uniform(0.0,2.0*glLibMathMisc::pi);
	double rand_radius = glLibMathRandom::get_uniform(0.0,radius*radius);
	double rand_radius_sqrt = sqrt(rand_radius);
	point[0] = rand_radius_sqrt * cos(rand_angle);
	point[1] = sqrt(1.0-rand_radius); //i.e., sqrt(1.0-rand_radius_sqrt*rand_radius_sqrt)
	point[2] = rand_radius_sqrt * sin(rand_angle);
}
template <typename type_float, typename type_int> void get_hemisphere_cosine_n(type_float point[3], type_int n, type_float radius=1) {
	//Global Illumination Compendium, pg. 19
	double rand_angle = glLibMathRandom::get_uniform(0.0,2.0*glLibMathMisc::pi);
	double rand_radius = glLibMathRandom::get_uniform(0.0,1.0);

	#ifdef SPEED
		type_float radical = sqrt(   1   -   fast_pow( rand_radius, 2.0/(n+1.0) )   );
	#else
		type_float radical = sqrt(   1   -   pow( rand_radius, 2.0/(n+1.0) )   );
	#endif

	point[0] = radius * cos(rand_angle) * radical;

	#ifdef SPEED
		point[1] = radius * fast_pow(rand_radius,1.0/(n+1.0));
	#else
		point[1] = radius * pow(rand_radius,1.0/(n+1.0));
	#endif

	point[2] = radius * sin(rand_angle) * radical;
}

class Renderer;
Renderer* instance__renderer;
inline int wrapper__render_tile(void*);
class Renderer {
	friend int wrapper__render_tile(void* arg);
	public:
		bool finished;
		bool stop;
	private:
		SDL_Surface* _image;
		int _img_w,_img_h;
		int _num_tiles_x,_num_tiles_y;
		double img_aspect;
		SDL_PixelFormat* _pixel_format;
		int _numof_samples_pixel, _numof_samples_aperture;

		int num_tiles_rendered; //Rendered.  Not loaded from cache
		#ifdef SAVE_INTERMEDIARY
		char loading_or_rendering;
		#endif

		int _thread_ids[NUM_RENDER_THREADS];

		//Must have a separate matrix for each thread, and must have a separate reflection rotation matrix for each recursion level
		glLibMathMatrix<double,char>* _matrices_rotation[NUM_RENDER_THREADS*MAX_DEPTH];
		glLibMathMatrix<double,char>* _matrices_result[NUM_RENDER_THREADS];
		glLibMathMatrix<double,char>* _matrices_point[NUM_RENDER_THREADS];

		Scene* _scene;
		double _background_color[3];

		map<string,Light*>* _lights;

		deque<Tile>* schedule;

		SDL_sem* schedule_lock;

		list<SDL_Thread*>* render_threads;

		double vec_eye_x[3],vec_eye_y[3]; //x, y
		double vec_eye_v[3],vec_eye_z[3],vec_img_center[3],vec_camera_position[3]; //view, scaled view, image center

		glLibTimer* _timer_total;
		glLibTimer* _timer_render;

		WeightedBVH* _bvh;

		PoissonPoint** _poisson_samples_pixel;
		PoissonPoint** _poisson_samples_aperture;

	private:
		void init_image_attributes(void) {
			_img_w = _scene->image->width;
			_img_h = _scene->image->height;
			img_aspect = (double)_img_w / (double)_img_h;

			if (TILE_SIZE*(_img_w/TILE_SIZE)==_img_w) { _num_tiles_x = _img_w/TILE_SIZE; }
			else _num_tiles_x = _img_w/TILE_SIZE+1;
			if (TILE_SIZE*(_img_h/TILE_SIZE)==_img_h) { _num_tiles_y = _img_h/TILE_SIZE; }
			else _num_tiles_y = _img_h/TILE_SIZE+1;

			_pixel_format = _image->format;
			_numof_samples_pixel = _scene->image->samples;
			_numof_samples_aperture = _scene->camera->aperture_samples;
		}
		void init_scene_attributes(void) {
			for (int i=0;i<3;++i) _background_color[i] = _scene->background_color[i];

			_lights = _scene->lights;
		}
		void init_vectors(void) {
			glLibMathVec::copy3(_scene->camera->position,vec_camera_position);

			glLibMathVec::subtract3(_scene->camera->center,vec_camera_position,vec_eye_v);
			glLibMathVec::normalize3(vec_eye_v);

			glLibMathVec::cross(vec_eye_v,_scene->camera->up,vec_eye_x);
			glLibMathVec::normalize3(vec_eye_x);

			glLibMathVec::cross(vec_eye_x,vec_eye_v,vec_eye_y); //vec_eye_y is normalized

			glLibMathVec::scale3(_scene->camera->image_plane,vec_eye_v,vec_eye_z);
			glLibMathVec::subtract3(vec_camera_position,vec_eye_z,vec_img_center);
		}
		void init_matrices(void) {
			for (unsigned long i=0;i<NUM_RENDER_THREADS*MAX_DEPTH;++i) {
				_matrices_rotation[i] = new glLibMathMatrix<double,char>(3,3,GLLIB_MATRIX_EMPTY);
			}
			for (unsigned long i=0;i<NUM_RENDER_THREADS;++i) {
				_matrices_result[i] = new glLibMathMatrix<double,char>(3,1,GLLIB_MATRIX_EMPTY);
				_matrices_point[i] = new glLibMathMatrix<double,char>(3,1,GLLIB_MATRIX_EMPTY);
			}
		}
		void init_threads(void) {
			schedule_lock = SDL_CreateSemaphore(1);

			schedule = new deque<Tile>();
			for (int x=0;x<_num_tiles_x;++x) {
				for (int y=0;y<_num_tiles_y;++y) {
					Tile new_tile;
					new_tile.x = x*TILE_SIZE;
					new_tile.y = y*TILE_SIZE;
					new_tile.w = TILE_SIZE;
					new_tile.h = TILE_SIZE;
					if (new_tile.x+new_tile.w>_img_w) new_tile.w = _img_w - new_tile.x;
					if (new_tile.y+new_tile.h>_img_h) new_tile.h = _img_h - new_tile.y;
					schedule->push_back(new_tile);
				}
			}
			//random_shuffle(schedule->begin(),schedule->end());

			render_threads = new list<SDL_Thread*>();
		}
		void init_random_samples(void) {
			int i;
			list<PoissonPoint*>* s;

			_poisson_samples_pixel = new PoissonPoint*[_numof_samples_pixel];
			i = 0;
			s = generate_poisson_square(_numof_samples_pixel);
			for (list<PoissonPoint*>::const_iterator iter=s->begin();iter!=s->end();++iter) {
				_poisson_samples_pixel[i++] = *iter;
			}
			delete s;

			_poisson_samples_aperture = new PoissonPoint*[_numof_samples_aperture*NUM_RANDOM_DISTR];
			i = 0;
			for (int j=0;j<NUM_RANDOM_DISTR;++j) {
				s = generate_poisson_disk(_numof_samples_aperture);
				for (list<PoissonPoint*>::const_iterator iter=s->begin();iter!=s->end();++iter) {
					PoissonPoint* p = *iter;
					p->x = _scene->camera->aperture*(p->x-0.5);
					p->y = _scene->camera->aperture*(p->y-0.5);
					_poisson_samples_aperture[i++] = p;
				}
				delete s;
			}
		}
		
		void deinit_matrices(void) {
			for (int i=0;i<NUM_RENDER_THREADS*MAX_DEPTH;++i) {
				delete _matrices_rotation[i];
			}
			for (int i=0;i<NUM_RENDER_THREADS;++i) {
				delete _matrices_result[i];
				delete _matrices_point[i];
			}
		}
		void deinit_threads(void) {
			stop = true;
			for (list<SDL_Thread*>::iterator iter=render_threads->begin();iter!=render_threads->end();++iter) {
				int status;
				SDL_WaitThread(*iter,&status);
			}
			delete render_threads;

			delete schedule;

			SDL_DestroySemaphore(schedule_lock);
		}
		void deinit_random_samples(void) {
			for (int i=0;i<_numof_samples_pixel;++i) {
				delete _poisson_samples_pixel[i];
			}
			delete [] _poisson_samples_pixel;

			for (int i=0;i<_numof_samples_aperture*NUM_RANDOM_DISTR;++i) {
				delete _poisson_samples_aperture[i];
			}
			delete [] _poisson_samples_aperture;
		}

	public:
		Renderer(Scene* scene, SDL_Surface* image) {
			_scene = scene;
			_bvh = NULL;
			_image = image;

			init_image_attributes();
			init_scene_attributes();
			init_vectors();
			init_matrices();
			init_threads();
			init_random_samples();

			finished = false;
			stop = false;
			num_tiles_rendered = 0;
			#ifdef SAVE_INTERMEDIARY
			loading_or_rendering = 0;
			#endif

			_timer_total = new glLibTimer();
			_timer_render = new glLibTimer();

			instance__renderer = this;
		}
		~Renderer() {
			delete _timer_total;
			delete _timer_render;

			deinit_threads();

			if (_bvh!=NULL) delete _bvh;

			deinit_random_samples();

			deinit_matrices();
		}

		void construct_scene(void) {
			#ifdef USING_BVH
			_bvh = new WeightedBVH(_scene->objects,16);
			printf("  Bounding volume hierarchy has:\n");
			printf("    Number of objects:  %u\n",_bvh->get_num_objects());
			printf("      Triangles:        %u\n",_bvh->get_num_objects_triangles());
			printf("      Spheres:          %u\n",_bvh->get_num_objects_spheres());
			printf("      Planes:           %u\n",_bvh->get_num_objects_planes());
			printf("    Leaf depth:         max %u, average %2.3f, possible %u\n",_bvh->get_max_depth(),_bvh->get_avg_depth(),_bvh->get_max_possible_depth());
			printf("    Leaf contents:      max %u, average %2.3f\n",_bvh->get_max_in_leaf(),_bvh->get_avg_in_leaf());
			#endif
		}
		
		void begin(void) {
			_timer_total->mark();
			for (int i=0;i<NUM_RENDER_THREADS;++i) {
				_thread_ids[i] = i;
				render_threads->push_front( SDL_CreateThread(wrapper__render_tile,&(_thread_ids[i])) );
			}
		}

		inline SDL_Surface* get_image(void) {
			return _image;
		}
		
		inline void print_time(void) {
			print_time(_timer_total->get(),2);
		}
		void print_time(double seconds, int precision, const char* label_d="d", const char* label_h="h", const char* label_m="m", const char* label_s="s") {
			int minutes = 0, hours = 0, days = 0;
			if (seconds/60.0>=1.0) {
				minutes = floor(seconds/60.0);
				seconds = fmod(seconds,60.0);
			}
			if (minutes/60>=1) {
				hours = minutes / 60;
				minutes = minutes % 60;
			}
			if (hours/60>=1) {
				days = hours / 24;
				hours = hours % 24;
			}
			if (days>1) printf("%d%s",days,label_d);
			if (hours>1) printf("%d%s",hours,label_h);
			if (minutes>1) printf("%d%s",minutes,label_m);
			printf(("%."+to_string(precision)+"f"+label_s).c_str(),seconds);
		}

	public:
		inline unsigned int trace_get_intersect(double origin[3], double ray[3], double* distance, Object** object_hit) {
			#ifdef USING_BVH
			return _bvh->intersect(origin,ray, distance,object_hit);
			#else
			*object_hit = NULL;
			*distance = MAX_DOUBLE;
			for (list<Object*>::const_iterator iter=_scene->objects->begin(); iter!=_scene->objects->end(); ++iter) {
				double intersect = intersect_object(iter,ray,origin);
				if (intersect<=0.0) continue;
				if (intersect<*distance) {
					*distance = intersect;
					*object_hit = *iter;
				}
			}
			return _scene->objects->size();
			#endif
		}
		inline double trace_shadow_ray(double origin[3], double ray[3]) {
			Object* object_hit; double closest;
			trace_get_intersect(origin,ray,&closest,&object_hit);
			return closest;
		}
		struct ray_info trace_ray(double origin[3], double ray[3], double color[3], int depth, deque<Material*>* material_stack, int thread_number) {
			//Determine intersect
			Object* object_hit;
			double closest;
			//_bvh->print();
			unsigned int num_tests = trace_get_intersect(origin,ray,&closest,&object_hit);
			#ifdef OUTPUT_BVH_HEATMAP
			add_color_heatmap(color,num_tests/15.0);
			#endif
			//Shade
			if (object_hit!=NULL) {
				struct ray_info result_ray_information = { closest, false, true };

				#ifdef OUTPUT_BVH_HEATMAP
				return result_ray_information;
				#endif
				
				#ifdef OUTPUT_T
				double part = (closest-OUTPUT_T_NEAR)/(OUTPUT_T_FAR-OUTPUT_T_NEAR);
				color[0] = part*OUTPUT_T_R;
				color[1] = part*OUTPUT_T_G;
				color[2] = part*OUTPUT_T_B;
				return result_ray_information;
				#endif

				Material* material_hit = object_hit->material;

				#ifdef ALLOW_AMBIENT
				//Ambient
				glLibMathVec::add3(material_hit->ambient,color,color);
				#endif

				double position[3];
				glLibMathVec::scale3(closest,ray,position);
				glLibMathVec::add3(position,origin,position);

				//r is aliased as the reflected ray or the refracted ray in the relevant blocks
				double l[3],n[3],r[3],v[3];

				glLibMathVec::negate3(ray,v);

				//Get normal
				object_hit->set_normal(position,n);

				#ifdef OUTPUT_NORMAL
				map_vec_to_rgb(n,color);
				return result_ray_information;
				#endif

				#if defined ALLOW_DIFFUSE || defined ALLOW_SPECULAR// || defined ALLOW_SHADOWS
				for (map<string,Light*>::const_iterator iter=_lights->begin();iter!=_lights->end();++iter) {
					Light* light = iter->second;
					glLibMathVec::subtract3(light->primitive->position,position,l);
					glLibMathVec::normalize3(l);

					double n_dot_l = glLibMathVec::dot3(n,l);
					if (n_dot_l<0.0) {
						//if (double_sided) {
						//	glLibMathVec::negate3(n);
						//	n_dot_l = -n_dot_l;
						//} else {
						continue;
						//}
					}

					double light_dependent_color[3], light_dependent_scalar = 1.0;

					#ifdef ALLOW_DIFFUSE
					//Diffuse
					double diffuse_color[3]; glLibMathVec::copy3(material_hit->diffuse,diffuse_color);
					SDL_Surface* texture = material_hit->texture;
					switch (object_hit->type) {
						case SCENE_OBJECT_TYPE_PLANE: {
							double sc = 4.0;
							if (((Plane*)(object_hit))->normal[1]==1) {
								if ((mod(sc*position[0],2.0)>1.0)^(mod(sc*position[2],2.0)>1.0)) {
									n_dot_l *= 0.5;//*max(mod(sc*position[0],1.0),mod(sc*position[2],1.0));
								}
							}
							else {
								if ((mod(sc*position[0],2.0)>1.0)^(mod(sc*position[1],2.0)>1.0)) {
									n_dot_l *= 0.5;//9*max(mod(sc*position[0],1.0),mod(sc*position[1],1.0));
								}
							}
							break;
						}
						case SCENE_OBJECT_TYPE_SPHERE: {
							if (texture!=NULL) {
								double uv[2] = {0.5*texture->w,0.5*texture->h};

								double angle_x = atan2(n[2],n[0]);
								double angle_y = atan2(n[1],sqrt(n[0]*n[0]+n[2]*n[2]));
								uv[0] = texture->w * (angle_x+glLibMathMisc::pi)/(2.0*glLibMathMisc::pi);
								uv[1] = texture->h * (angle_y + glLibMathMisc::pi*0.5)/glLibMathMisc::pi;
								double texture_color[3];
								sample_default(uv[0],uv[1],texture,texture_color);
								glLibMathVec::multiply3(texture_color,diffuse_color,diffuse_color);
							}
							break;
						}
						case SCENE_OBJECT_TYPE_TRIANGLE: {
							if (texture!=NULL) {
								double uv[2] = {0.5*texture->w,0.5*texture->h};

								Triangle* t = (Triangle*)(object_hit);
								t->set_uv(position,uv);
								double texture_color[3];
								sample_default(uv[0]*texture->w,uv[1]*texture->h,texture,texture_color);
								glLibMathVec::multiply3(texture_color,diffuse_color,diffuse_color);
							}
							break;
						}
						default: throw glLibError("Unrecognized type!");
					}

					#ifdef OUTPUT_UV_DIFFUSE
					color[0] = uv[0]/texture->w;
					color[1] = uv[1]/texture->h;
					color[2] = 0.0;
					return result_ray_information;
					#endif

					double scaled_diffuse[3]; glLibMathVec::scale3(n_dot_l,diffuse_color,scaled_diffuse);
					glLibMathVec::copy3(scaled_diffuse,light_dependent_color);
					#endif

					#ifdef ALLOW_SPECULAR
					//Specular
					//	Not using glLibMathVec::reflect, because can take advantage of the fact that dot(n,l) has already been computed.
					double temp[3];
					glLibMathVec::scale3(2.0*n_dot_l,n,temp);
					glLibMathVec::subtract3(temp,l,r);

					double spec = glLibMathVec::dot3(r,v);
					spec = (spec<0.0) ? 0.0 : spec; //when raising small negative numbers to large powers, sometimes it becomes positive.
					/*#ifdef SPEED
					spec = fast_pow(spec,material_hit->shininess);
					#else
					spec = pow(spec,material_hit->shininess);
					#endif*/
					//21.85007974481659% faster
					int pow_index = (int)(spec*(NUM_POW_TABLE-1));
					spec = (material_hit->pow_table[pow_index]+material_hit->pow_table[pow_index+1])/2.0;

					double scaled_specular[3]; glLibMathVec::scale3(spec,material_hit->specular,scaled_specular);
					glLibMathVec::add3(scaled_specular,light_dependent_color,light_dependent_color);
					#endif

					#ifdef ALLOW_SHADOWS
					//Shadow
					map<Light*,int>::const_iterator location = material_hit->to_shadow->find(light);
					if (location!=material_hit->to_shadow->end()) {
						light_dependent_scalar = 0.0;
						for (int i=0;i<location->second;++i) {
							double light_point[3];
							light->get_random_point(light_point);

							double shadow_ray_direction[3];
							glLibMathVec::subtract3(light_point,position,shadow_ray_direction);

							double distance_to_light_point = glLibMathVec::length3(shadow_ray_direction);
							glLibMathVec::scale3(1.0/distance_to_light_point,shadow_ray_direction);

							double offset_ray[3], offset_position[3];
							glLibMathVec::scale3(EPSILON,shadow_ray_direction,offset_ray);
							glLibMathVec::add3(position,offset_ray,offset_position);
							double intersect_distance = trace_shadow_ray(offset_position,shadow_ray_direction);
							light_dependent_scalar += ((intersect_distance==0.0)||(intersect_distance>distance_to_light_point)) ? 1.0 : 0.0;
						}
						light_dependent_scalar /= (double)(location->second);
					}
					#endif

					//Add (possibly scaled) diffuse and specular
					glLibMathVec::scale3(light_dependent_scalar,light_dependent_color);
					glLibMathVec::multiply3(light_dependent_color,light->color,light_dependent_color);
					glLibMathVec::add3(light_dependent_color,color,color);
				}
				#endif

				#ifdef ALLOW_REFLECTIONS
				//Reflection
				if (material_hit->reflect&&depth<MAX_DEPTH) {
					double reflection[3] = {0.0,0.0,0.0};

					glLibMathVec::reflect3(v,n,r);

					#ifdef OUTPUT_REFLECTION
					map_vec_to_rgb(r,color);
					return result_ray_information;
					#endif

					double offset_ray[3], offset_position[3];
					glLibMathVec::scale3(EPSILON,r,offset_ray);
					glLibMathVec::add3(position,offset_ray,offset_position);

					if (material_hit->samples_reflection>1) {
						//TODO: define glLibMathVec::y_R3, x_R3, etc.
						double world_y[3] = {0.0,1.0,0.0};
						double tangent[3]; glLibMathVec::cross(world_y,r,tangent); glLibMathVec::normalize3(tangent);
						double distribution_angular_offset = acos(glLibMathVec::dot3(world_y,r));
						
						glLibMathMatrix<double,char>* rot_matrix = _matrices_rotation[thread_number*MAX_DEPTH+depth];
						glLibMathMatrix<double,char>::set_rotation_matrix_rad(rot_matrix,tangent,distribution_angular_offset);
						for (int i=0;i<material_hit->samples_reflection;++i) {
							double random_cos_vec[3], perturbed_reflection[3];

							get_hemisphere_cosine_n(random_cos_vec,material_hit->shininess,1.0);
							for (char j=0;j<3;++j) _matrices_point[thread_number]->values[j][0] = random_cos_vec[j];

							glLibMathMatrix<double,char>::multiply(
								rot_matrix,
								_matrices_point[thread_number],
								_matrices_result[thread_number]
							);
							for (char j=0;j<3;++j) perturbed_reflection[j] = _matrices_result[thread_number]->get_at(0,j);
							
							if (glLibMathVec::dot3(perturbed_reflection,n)<0.0) {
								glLibMathVec::negate3(perturbed_reflection);
								glLibMathVec::reflect3(perturbed_reflection,n,perturbed_reflection);
							}

							deque<Material*> material_stack_temp = *material_stack;
							trace_ray(offset_position,perturbed_reflection,reflection,depth+1,&material_stack_temp, thread_number);
						}

						glLibMathVec::scale3(1.0/(double)material_hit->samples_reflection,reflection);
					}
					else {
						deque<Material*> material_stack_temp = *material_stack;
						trace_ray(offset_position,r,reflection,depth+1,&material_stack_temp, thread_number);
					}

					glLibMathVec::scale3(material_hit->coefficient_reflection,reflection);
					glLibMathVec::add3(reflection,color,color);
				}
				#endif

				#ifdef ALLOW_REFRACTIONS
				//Refraction
				if (material_hit->refract&&depth<MAX_DEPTH) {
					double refraction[3] = {0.0,0.0,0.0};

					char entering_leaving = (glLibMathVec::dot3(v,n)>0.0) ? 1 : 2;

					if (entering_leaving==1) { //entering something
						if (material_stack->empty()) { //vacuum to object
							glLibMathVec::refract3(ray,n,r,1.0/material_hit->refractive_index);
						}
						else { //superobject to subobject
							Material* material_current = material_stack->front();
							glLibMathVec::refract3(ray,n,r,material_current->refractive_index/material_hit->refractive_index);
						}
						#ifdef OUTPUT_REFRACTION
						map_vec_to_rgb(r,color);
						return result_ray_information;
						#endif
						material_stack->push_front(material_hit);
					}
					else { //leaving something
						if (material_stack->empty()) {
							//Numerical error caused an entry/exit to be less than an epsilon?
							//printf("Uh oh!\n");
							//refraction[0] = 10000.0;
							//goto SET_REFRACTION;
							result_ray_information.valid = false;
							return result_ray_information;
						}
						Material* material_current = material_stack->front();
						material_stack->pop_front();

						double temp[3]; glLibMathVec::copy3(n,temp);
						glLibMathVec::negate3(n,temp);

						if (material_stack->empty()) { //object to vacuum
							glLibMathVec::refract3(ray,temp,r,material_current->refractive_index/1.0);
						}
						else { //subobject to superobject
							glLibMathVec::refract3(ray,temp,r,
								material_current->refractive_index/material_stack->front()->refractive_index);
						}
					}

					double offset_ray[3], offset_position[3];
					glLibMathVec::scale3(EPSILON,r,offset_ray);
					glLibMathVec::add3(position,offset_ray,offset_position);

					if (material_hit->samples_refraction>1) {
						double world_y[3] = {0.0,1.0,0.0};
						double tangent[3]; glLibMathVec::cross(world_y,r,tangent); glLibMathVec::normalize3(tangent);
						double distribution_angular_offset = acos(glLibMathVec::dot3(world_y,r));
						
						glLibMathMatrix<double,char>* rot_matrix = _matrices_rotation[thread_number*MAX_DEPTH+depth];
						glLibMathMatrix<double,char>::set_rotation_matrix_rad(rot_matrix,tangent,distribution_angular_offset);
						int num_valid_samples = 0;
						double avg_dist = 0.0;
						for (int i=0;i<material_hit->samples_refraction;++i) {
							double random_cos_vec[3], perturbed_refraction[3];

							get_hemisphere_cosine_n(random_cos_vec,material_hit->exponent_refractive,1.0);
							for (char j=0;j<3;++j) _matrices_point[thread_number]->values[j][0] = random_cos_vec[j];

							glLibMathMatrix<double,char>::multiply(
								rot_matrix,
								_matrices_point[thread_number],
								_matrices_result[thread_number]
							);
							for (char j=0;j<3;++j) perturbed_refraction[j] = _matrices_result[thread_number]->get_at(0,j);
							
							/*if (glLibMathVec::dot(perturbed_refraction,n,3)>0.0) {
								glLibMathVec::negate(perturbed_refraction,3);
								glLibMathVec::reflect(perturbed_refraction,n,perturbed_refraction,3);
							}*/

							deque<Material*> material_stack_temp = *material_stack;
							struct ray_info refract_result = trace_ray(
								offset_position,
								perturbed_refraction,refraction,depth+1,&material_stack_temp, thread_number);
							if (!refract_result.empty) {
								++num_valid_samples;
								avg_dist += refract_result.distance;
							}
						}
						avg_dist /= num_valid_samples;
						glLibMathVec::scale3(1.0/(5.0*avg_dist),refraction);
						glLibMathVec::scale3(1.0/(double)num_valid_samples,refraction);
					}
					else {
						deque<Material*> material_stack_temp = *material_stack;
						struct ray_info refract_result = trace_ray(
							offset_position,r,refraction,depth+1,&material_stack_temp, thread_number);
						refraction[0] *= exp(-0.30*refract_result.distance);
						refraction[1] *= exp(-0.15*refract_result.distance);
						refraction[2] *= exp(-0.20*refract_result.distance);
					}

					glLibMathVec::scale3(material_hit->coefficient_refraction,refraction);
					//SET_REFRACTION:
					//double n1 = material_hit->refractive_index;
					//double n2 = 1.0;
					//fresnel_rs(n1,n2,);
					glLibMathVec::add3(refraction,color,color);
				}
				#endif

				return result_ray_information;
			}
			struct ray_info result_ray_information = {0.0,true,true};
			return result_ray_information;
		}

		int render_tile(void* data) {
			static int thread_number = *((int*)(data));

			Tile to_render;
			deque<Material*> material_stack;

			#ifdef SAVE_INTERMEDIARY
			string tile_name;
			SDL_Surface* current_tile = SDL_CreateRGBSurface(SDL_SWSURFACE,TILE_SIZE,TILE_SIZE,24,0,0,0,0);
			#endif
			while (true) {
				SDL_SemWait(schedule_lock);
				#ifdef SAVE_INTERMEDIARY
				do {
					if (schedule->size()==0||stop) {
						finished = true;
						break;
					}
					to_render = schedule->front();
					schedule->pop_front();
					tile_name = RENDER_TEMPDIR_NAME + "/" + to_string(to_render.x) + "_" + to_string(to_render.y) + ".bmp";
					//printf("Checking %d,%d = %s\n",to_render.x,to_render.y,tile_name.c_str());
					if (exists_file(tile_name)) {
						printf("\r  Loading tile (%d,%d) from cache.                ",to_render.x,to_render.y);
						SDL_Surface* temp = SDL_LoadBMP(tile_name.c_str());
						SDL_SemPost(schedule_lock);
						for (int y=0;y<to_render.h;++y) {
							for (int x=0;x<to_render.w;++x) {
								set_pixel_rgb24(
									to_render.x+x,
									_img_h-(to_render.y+y)-1,
									get_pixel_rgb24(x,TILE_SIZE-y-1,temp)
								);
							}
						}
						SDL_SemWait(schedule_lock);
						SDL_FreeSurface(temp);
					}
					else {
						if (loading_or_rendering==0) {
							loading_or_rendering = 1;
							_timer_render->mark();
						}
						break;
					}
					//printf("  Exists!  %d %d\n",is_directory(tile_name.c_str()),is_file(tile_name.c_str()));
				} while (true);
				#else
				if (schedule->size()==0) finished = true;
				else {
					to_render = schedule->front();
					schedule->pop_front();
				}
				#endif
				if (!(finished||stop)) {
					printf("\r  Rendering tile (%d,%d); expected: ",to_render.x,to_render.y);
					print_time(  _timer_total->get() / (double)(num_tiles_rendered) * schedule->size(),  0  );
					printf("                ");
				}
				SDL_SemPost(schedule_lock);
				if (finished||stop) {
					printf("\n");
					break;
				}
				
				#ifdef ALLOW_RENDER
				for (int pix_x=to_render.x;pix_x<to_render.x+to_render.w;pix_x+=PIXEL_STEP) {
					for (int pix_y=to_render.y;pix_y<to_render.y+to_render.h;pix_y+=PIXEL_STEP) {

						double color[3] = {0.0,0.0,0.0};
						for (int i=0;i<_numof_samples_pixel;++i) {
							material_stack.clear();

							double x_part = ((double)(pix_x)+_poisson_samples_pixel[i]->x) / (double)_img_w * -2.0 + 1.0;
							double y_part = ((double)(pix_y)+_poisson_samples_pixel[i]->y) / (double)_img_h * -2.0 + 1.0;
							y_part /= img_aspect;

							double ray_position[3] = {
								x_part*vec_eye_x[0] + y_part*vec_eye_y[0] + vec_img_center[0],
								x_part*vec_eye_x[1] + y_part*vec_eye_y[1] + vec_img_center[1],
								x_part*vec_eye_x[2] + y_part*vec_eye_y[2] + vec_img_center[2]
							};

							//Calculate focal point of this sample in the scene
							// 1/s1 + 1/s2 = 1/f
							//s2 is _scene->camera->image_plane
							//s1 is the camera's z distance from the lens to the focal point
							double central_ray[3], central_ray_prime[3]; //Rays through the lens' center
							glLibMathVec::subtract3(vec_camera_position,ray_position,central_ray);
							static double s2 = _scene->camera->image_plane;
							static double s1 = _scene->camera->recip_focal_length - 1.0/s2;
							//central_ray_prime is parallel to central_ray, but goes exactly to the focal point
							// ||central_ray||/s1 = ||central_ray_prime||/(s1+s2)
							glLibMathVec::scale3((s1+s2)*glLibMathVec::length3(central_ray)/s1,central_ray,central_ray_prime);
							double focal_point[3]; glLibMathVec::add3(ray_position,central_ray_prime,focal_point);

							int distr_index = _numof_samples_aperture*(rand()%NUM_RANDOM_DISTR);
							for (int j=0;j<_numof_samples_aperture;++j) {
								double ray_direction[3]; //ray_direction is also used as a temp
								glLibMathVec::scale3(_poisson_samples_aperture[distr_index+j]->x,vec_eye_x,ray_direction); //ray_direction is used as temp
								glLibMathVec::add3(vec_camera_position,ray_direction,ray_position);
								glLibMathVec::scale3(_poisson_samples_aperture[distr_index+j]->y,vec_eye_y,ray_direction); //ray_direction is used as temp
								glLibMathVec::add3(ray_direction,ray_position,ray_position);

								//Ray starts on the lens, heads toward focal point
								glLibMathVec::subtract3(focal_point,ray_position,ray_direction); //ray_direction gets set for real
								glLibMathVec::normalize3(ray_direction);

								#ifndef ALLOW_BACKGROUND
								trace_ray(ray_position,ray_direction,color,0,&material_stack, thread_number);
								#else
								if (trace_ray(ray_position,ray_direction,color,0,&material_stack, thread_number).empty) {
									color[0] += _background_color[0];
									color[1] += _background_color[1];
									color[2] += _background_color[2];
								}
								#endif
							}
						}
						glLibMathVec::scale3(255.0/(_numof_samples_pixel*_numof_samples_aperture),color);
						glLibMathVec::clamp3(color,0.0,255.0);

						Uint32 final_color = SDL_MapRGB(_pixel_format,(int)(color[0]),(int)(color[1]),(int)(color[2]));
						set_pixel_rgb24(pix_x,_img_h-pix_y-1,final_color,_image);
						#ifdef SAVE_INTERMEDIARY
						set_pixel_rgb24(pix_x-to_render.x,TILE_SIZE-(pix_y-to_render.y)-1,final_color,current_tile);
						#endif
					}
				}
				#else
				SDL_Delay(RENDER_TILE_DELAY);
				#endif

				#ifdef SAVE_INTERMEDIARY
				save(tile_name.c_str(),current_tile);
				#endif
				++num_tiles_rendered;
			}

			#ifdef SAVE_INTERMEDIARY
			SDL_FreeSurface(current_tile);
			#endif

			finished = true;

			return 1;
		}
};
inline int wrapper__render_tile(void* arg) {
	return instance__renderer->render_tile(arg);
}