#ifndef SCHLICKPHOTONSHADER_HXX
 #define SCHLICKPHOTONSHADER_HXX
 
 #include "Shader.hxx"
 #include "Material.hxx"
 #include "PhotonKDTree.hxx"
 
 #include "defines.h"
  
 // This shader is equivalent to the photon shader, but uses photon maps for indirect illumination
 class SchlickPhotonShader : public Shader
 {
 	Randomizer randomizer;
 public:
 	Vec3f        ambient,             // ambient radiance 
 	             light_ray_org_std;   // template for light_ray.org, to calculate it only once
 	float        p12,                 // square of the isotropy factor for primary surface
 	             p22;                 // square of the isotropy factor for primary surface
 	float        mirror,              // Material based mirroring    percentage
 	             glossy_mirror,       // Material based glossy reflection percentage
 	             glossy_transparency, // Material based glossy refraction percentage
 	             transparency;        // Material based transparency percentage
 	PhotonKDTree *photon_tree;
 
 	SchlickPhotonShader(Scene *scene, Material *material, PhotonKDTree *photon_tree)
 		: Shader(scene, material), photon_tree(photon_tree)
 	{
 		// Ambient illumination part
 		ambient = material->ambient * Product(material->color, La);
 
 		// Initialization of reflection and refraction
 		mirror              = (material->mirror              == UNDEF) ? 0.0f : material->mirror;
 		glossy_mirror       = (material->glossy_mirror       == UNDEF) ? 0.0f : material->glossy_mirror;
 		glossy_transparency = (material->glossy_transparency == UNDEF) ? 0.0f : material->glossy_transparency;
 		transparency        = (material->transparency        == UNDEF) ? 0.0f : material->transparency;
 
 		if (mirror + glossy_transparency + glossy_mirror + transparency > 1.0)
 			cerr << "WARNING: Material '" << material << "' has wrong reflection/refraction percentages (> 1.0)!" << endl;
 
 		// Precomputation of p^2
 		p12 = material->p1 * material->p1;
 		p22 = material->p2 * material->p2;
 	};
 
 	virtual ~SchlickPhotonShader()
 	{};
 
 #include "SchlickTerm.include"
 
 	virtual Vec3f FancyColor(const Ray &origray, const Vec3f &normal, const bool &into)
 	{
 		Ray   ray;
 		float m = mirror,
 		      t = transparency,
 		      gm = glossy_mirror,
 		      gt = glossy_transparency;
 
 		// Project the ray direction onto the normal vector...
 		Vec3f projected   = Dot(normal, -origray.dir) * normal,
 		      std_ray_org = origray.org + origray.t * origray.dir; 
 
 
 		// Start with the reflection calculation (which is needed anyway for t/m)
 		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 		Vec3f reflective_color        = Vec3f(0),
 		      refractive_color        = Vec3f(0),
 		      glossy_reflection_color = Vec3f(0),
 		      glossy_refraction_color = Vec3f(0);
 
 		if (t || m)
 		{
 			ray.org = std_ray_org;
 			ray.t   = Infinity;
 			ray.hit = NULL;
 
 			// ...and mirror away
 			ray.dir = origray.dir + 2 * projected;
 			if (Nonzero(ray.dir))
 			{
 				Normalize(ray.dir);
 
 				ray.refrac_index = origray.refrac_index;
 				reflective_color = scene->RayTrace(ray);
 			}
 			else
 				reflective_color = Vec3f(0);
 		}
 
 		// Now perform the transparency part
 		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 		if (t)
 		{
 			ray.org = std_ray_org;
 			ray.t   = Infinity;
 			ray.hit = NULL;
 
 			// Cosine of the incident angle
 			float cosa       = Dot(normal, -origray.dir),
 			      new_refrac_index;
 
 			if (into)
 				new_refrac_index = origray.refrac_index + (material->refrac_index - 1.0f);
 			else
 				new_refrac_index = origray.refrac_index - (material->refrac_index - 1.0f);
 
 			float n        = origray.refrac_index / new_refrac_index,
 
 			/*
 			 * According to Snell's law we know that
 			 *       sina/sinb = n2/n1 = 1/n,
 			 * hence
 			 *       (sina/sinb)^2 = (1/n)^2
 			 * or
 			 *       sin^2 b = sin^2 a * n^2.
 			 * As
 			 *       sin^2 a + cos^2 a = 1,
 			 * we know that
 			 *       sin^2 b = n^2 * (1 - cos^2 a)
 			 * and thus
 			 *       cos^2 b = 1 - n^2 * (1 - cos^2 a).
 			 */
 
 			      cosb2 = 1.0f - n * n * (1.0f - cosa*cosa);
 
 			if (cosb2 < 0.0f)
 				// Total reflection
 				m += t;
 			else
 			{
 				// Here again, we use the Schlick approximation
 				// (Hey, we are in the Schlick shader ;-) !)
 
 				float cosb    = sqrtf(cosb2),
 				      c       = (n - 1) / (n + 1),
 				      schlick;
 
 				c *= c;
 
 				// We use the angle in the less dense medium for Schlick's term
 				if (into)
 					schlick = c + (1 - c) * powf((1 - cosa), 5);
 				else
 					schlick = c + (1 - c) * powf((1 - cosb), 5);
 
 				schlick *= t;
 
 				m += schlick;
 				t -= schlick;
 
 				ray.dir = (n * cosa - cosb) * normal + n * origray.dir;
 
 				// Apply refraction index now as we are sure to have changed the medium
 				ray.refrac_index = new_refrac_index;
 			}
 
 			refractive_color = scene->RayTrace(ray);
 		}
 
 		if (gm || gt)
 		{
 			// Glossy reflection... think what? Yes... Same as in the photon distributor!
 			// Have fun...
 
 			// Get a new random direction
 			// ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 			// Randomize the parametrized direction
 			float incident_part = randomizer.drand() * 2 - 1,
 			      normal_part   = randomizer.drand();
 
 			// Span the incident plane
 			Vec3f incident_direction,
 			      new_direction;
 			if (projected == normal)
 				new_direction = normal;
 			else
 			{
 				incident_direction = origray.dir + projected;
 				Normalize(incident_direction);
 
 				// Deparametrize this direction
 				new_direction = incident_part * incident_direction + normal_part * normal;
 				if (Nonzero(new_direction))
 					Normalize(new_direction);
 			}
 
 			// Calculate the inverse BRDF and act accordingly
 			// We do this in ways of the Schlick model (Compare Schlick Shader for further details)
 			// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 			float specular;
 
 			// Get the halfway vector in between incoming and outgoing ray
 			Vec3f H = new_direction - origray.dir;
 			if (Nonzero(H))
 			{
 				Normalize(H);
 
 				float t, u, v, x, w;
 
 				// Calculation of various scaling factors
 				t = fabs(Dot(H,             normal));
 				u = fabs(Dot(-origray.dir,  H));
 				v = fabs(Dot(-origray.dir,  normal));
 				x = fabs(Dot(new_direction, normal));
 				w = fabs(Dot(H - normal,    normal));
 
 				// Precomputation of w^2
 				w = w * w;
 
 				specular = SchlickTerm(t, u, v, x, w);
 			}
 			else
 				specular = 0.0f;
 
 			if (specular)
 			{
 				if (gt)
 				{
 					ray.org = std_ray_org;
 					ray.t   = Infinity;
 					ray.hit = NULL;
 					ray.dir = new_direction - 2 * Dot(new_direction, normal);
 
 					// Cosine of the incident angle
 					float cosa  = Dot(normal, ray.dir),
 					      new_refrac_index;
 
 					if (into)
 						new_refrac_index = origray.refrac_index + (material->refrac_index - 1.0f);
 					else
 						new_refrac_index = origray.refrac_index - (material->refrac_index - 1.0f);
 
 					if (new_refrac_index)
 					{
 						float n     = origray.refrac_index / new_refrac_index,
 						      cosb2 = 1.0f - n * n * (1.0f - cosa*cosa);
 
 						if (cosb2 < 0.0f)
 						{
 							// Total reflection
 							gm += gt;
 							gt  = 0;
 						}
 						else
 						{
 							// Here again, we use the Schlick approximation
 							// (Hey, we are in the Schlick shader ;-) !)
 
 							float cosb    = sqrtf(cosb2),
 							      c       = (n - 1) / (n + 1),
 							      schlick;
 
 							c *= c;
 
 							// We use the angle in the less dense medium for Schlick's term
 							if (into)
 								schlick = c + (1 - c) * powf((1 - cosa), 5);
 							else
 								schlick = c + (1 - c) * powf((1 - cosb), 5);
 
 							schlick *= gt;
 
 							gm += schlick;
 							gt -= schlick;
 
 							ray.dir = (n * cosa - cosb) * normal + n * ray.dir;
 
 							// Apply refraction index now as we are sure to have changed the medium
 							ray.refrac_index = new_refrac_index;
 
 							glossy_refraction_color += (15.0f * specular) * scene->RayTrace(ray);
 
 							// Apply optional color filter
 							if (material->glossy_transparency_filter != UNDEF)
 							{
 								glossy_refraction_color = (material->glossy_transparency_filter
 											* (glossy_refraction_color.x + glossy_refraction_color.y + glossy_refraction_color.z) / 3.0f)
 											* material->color
 											+ (1.0f - material->glossy_transparency_filter) * glossy_refraction_color;
 							}
 						}
 					}
 				}
 
 				if (gm)
 				{
 					ray.org = std_ray_org;
 					ray.t   = Infinity;
 					ray.hit = NULL;
 					ray.dir = new_direction;
 
 					// What is this 15.0f about? Well... we look out in all directions where there is only little
 					// irradiance according to the BRDF. This value compensates for that.
 					glossy_reflection_color += (15.0f * specular) * scene->RayTrace(ray);
 
 					// Apply optional color filter
 					if (material->glossy_mirror_filter != UNDEF)
 					{
 						glossy_reflection_color = (material->glossy_mirror_filter
 									* (glossy_reflection_color.x + glossy_reflection_color.y + glossy_reflection_color.z) / 3.0f)
 									* material->color
 									+ (1.0f - material->glossy_mirror_filter) * glossy_reflection_color;
 					}
 				}
 			}
 			else
 			{
 				glossy_reflection_color = Vec3f(0);
 				glossy_refraction_color = Vec3f(0);
 			}
 		}
 
 		// Apply optional color filter
 		if (material->refrac_filter != UNDEF)
 		{
 			refractive_color = (material->refrac_filter * (refractive_color.x + refractive_color.y + refractive_color.z) / 3.0f)
 					                                    * material->color
 			                 + (1.0f - material->refrac_filter) * refractive_color;
 		}
 
 		// Apply optional color filter
 		if (material->mirror_filter != UNDEF)
 		{
 			reflective_color = (material->mirror_filter * (reflective_color.x + reflective_color.y + reflective_color.z) / 3.0f)
 					                                    * material->color
 			                 + (1.0f - material->mirror_filter) * reflective_color;
 		}
 
 		return t * refractive_color + m * reflective_color + gm * glossy_reflection_color + gt * glossy_refraction_color;
 	}
 
 
 	virtual Vec3f Shade(Ray &ray)
 	{
 		Vec3f Lr(0,0,0),
 		      Ll(0.0f),
 		      H,
 		      N = ray.hit->GetNormal(ray),
 		      pixelColor;
 		Ray   light_ray;
 		float Lr_fac,            // Factor of Lr, for time improvement
 		      t,                 // +
 		      u,                 //  |
 		      v,                 //   > cosines of various angles between normal vectors
 		      x,                 //  |
 		      w;                 // +
 		bool  into;
 
 		into = true;
 		if (Dot(N,ray.dir) > 0)
 		{
 			N = -N;
 			into = false;
 		}
 
 		light_ray_org_std = ray.org + ray.t * ray.dir;
 
 		for (vector<Procedural*>::iterator it = material->procedural.begin(); it != material->procedural.end(); ++it)
 			if ((*it)->ProvidesBumpMap())
 				(*it)->BumpMap(N, light_ray_org_std);
 
 		#define NEIGHBORS 1000
 		#define MAXDIST2  1.0f
 		int         count = 0;
 
 		Normalize(ray.dir);
 
 		// Integrate over a neighborhood of NEIGHBORS elements within MAXDIST2
 		PhotonHeap heap(light_ray_org_std, NEIGHBORS, MAXDIST2);
 		photon_tree->GetNextPhotons(heap);
 
 		if (material->c1 != UNDEF)
 		{
 			for (int i = 0; i < NEIGHBORS; ++i)
 				if (heap[i])
 				{
 					// Check whether pure photon mapping or simple subsurface scattering should be applied
 					Vec3f photon_dir = heap[i]->point - light_ray_org_std,
 					      photon_len = photon_dir;
 
 					if (Nonzero(photon_dir))
 						Normalize(photon_dir);
 
 					if (Dot(N, photon_dir) < -0.15)
 					{
 						const float maxdist = sqrtf(MAXDIST2);
 
 						// We are in a frustrum of 2*80 degrees wrt the surface normal, this is
 						// most probably a problem for subsurface scattering!
 						// Subsurface scattering? Err... well... sort of. This is actually a quite risky approximation
 						// and does only hold for rather thick objects. But hey, we don't want any points for it, okay?
 						// It is simply for the sake of cool images!
 						float weight = 0;
 
 						if (Nonzero(photon_len))
 							weight = Length(photon_len)/maxdist;
 
 						Lr += (1.0f - weight) * material->subsurface_scat * Product(heap[i]->energy, material->color); 
 					}
 					else
 					{
 						// These photons were gathered around the hitpoint and will be responsible for
 						// "classical" indirect illumination
 						light_ray.dir = -heap[i]->dir;
 						Ll  = heap[i]->energy;
 
 						count++;
 
 						if (Nonzero(light_ray.dir))
 						{
 							Normalize(light_ray.dir);
 
 // This was a test we did: Lambertian illumination. We got better results when we weighted the
 // incoming illumination directions and energies, why we discarded this one
 #if 0
 							float weight = Dot(light_ray.dir, ray.dir);
 							if (weight > 0.0f)
 								Ll *= weight;
 							else
 								Ll = Vec3f(0.0f);
 #endif
 
 							Lr_fac = 0;
 
 							// Get the halfway vector in between incoming and outgoing ray
 							H = light_ray.dir - ray.dir;
 							Normalize(H);
 
 							// Calculation of various scaling factors
 							t = fabs(Dot(H,             N));
 							u = fabs(Dot(light_ray.dir, H));
 							v = fabs(Dot(light_ray.dir, N));
 							x = fabs(Dot(-ray.dir,      N));
 							w = fabs(Dot(H - N,         N));
 
 							// Precomputation of w^2
 							w = w * w;
 
 							Lr_fac = SchlickTerm(t, u, v, x, w);
 
 							pixelColor = Product(material->color, Ll);
 
 							// This is the one "bitmap" texture
 							if (material->texture)
 							{
 								float u = 0,
 								      v = 0;
 
 								ray.hit->GetUV(ray, u, v);
 								v = 1.0f - v;
 
 								pixelColor = Product(pixelColor, material->texture->GetTexel(u,v));
 							}
 
 							// Check for (multiple) procedural textures to be applied...
 							for (vector<Procedural*>::iterator it = material->procedural.begin(); it != material->procedural.end(); ++it)
 								if ((*it)->ProvidesTexture())
 									pixelColor = Product(pixelColor, (*it)->GetColor(light_ray_org_std));
 
 							Lr += Lr_fac * pixelColor;
 						}
 					}
 				}
 
 			float scale = heap.maxdist2 * M_PI;
 
 			if (count && scale)
 				Lr *= 1.0f / (scale * static_cast<float>(count));
 		}
 
 		if (mirror || transparency || glossy_mirror || glossy_transparency)
 		{
 			Lr *= (1.0f - mirror - transparency - glossy_mirror - glossy_transparency);
 			Lr += FancyColor(ray, N, into);
 		}
 
 		Lr = Min(Vec3f(1.0f), Max(Vec3f(0.0), Lr));
 
 		Lr += ambient;
 		return Lr;
 	};
 };
 
 #endif