#include <iostream>
#include <string>
#include <stdlib.h>
#include <time.h>


using namespace std;

#include "pbrt_wrapper.cpp"

#ifndef FLT_EPSILON
#define FLT_EPSILON 1.192092896e-07F
#endif

#ifndef FLT_MAX
#define FLT_MAX         3.402823466e+38F        /* max value */
#endif

float RandomFloat()
{
	return (float)rand()/(float)RAND_MAX;
}

// L_r(1): specular reflection on the first hit:
//		   when cos(wo)>0, i.e., entering, it is the contaminant-air interface
//		   when cos(wo)<0, i.e., outgoing, it is the contaminant-glass interface
//
class DustyGlassSpecularReflectionDust
{
public:
	// DustyGlassSpecularReflectionDust Public Methods
	DustyGlassSpecularReflectionDust(float r_glass, float r_dust,
		float ior_air, float ior_dust, float ior_glass):
		R_glass(r_glass), R_dust(r_dust),  
		fresnel_air_dust(ior_air, ior_dust), 
		fresnel_dust_glass(ior_dust, ior_glass){}

		float Sample_f(const Vector &wo, Vector *wi,
			float u1, float u2, float *pdf) const;

private:
	// DustyGlassSpecularReflectionDust Private Data
	float R_glass, R_dust;
	FresnelDielectric fresnel_air_dust;
	FresnelDielectric fresnel_dust_glass;
};

// L_r(2): attenuated specular reflection (on the second hit)
//		when cos(wo)>0, light path is air-->contaminant-->glass-->contaminant-->air
//		when cos(wo)<0, light path is glass-->contaminant-->air-->contaminant-->glass
//
//		so it involves two refractions, one reflection, and two attenuations.
//
class DustyGlassSpecularReflection 
{
public:
	// DustyGlassSpecularReflection Public Methods
	DustyGlassSpecularReflection(float r_glass, 
		float t_glass, float r_dust,float t_dust, 
		float ior_air, float ior_dust, float ior_glass,
		float tau_dust):
		R_glass(r_glass), T_glass(t_glass), R_dust(r_dust), T_dust(t_dust), 
		fresnel_air_dust(ior_air, ior_dust),
		fresnel_dust_glass(ior_dust, ior_glass),
		Tau(tau_dust){}
		float Sample_f(const Vector &wo, Vector *wi,
			float u1, float u2, float *pdf) const;
private:
	// DustyGlassSpecularReflection Private Data
	float R_glass, T_glass, R_dust, T_dust;
	FresnelDielectric fresnel_air_dust;
	FresnelDielectric fresnel_dust_glass;
	float Tau;
};

// L_t(1): attenuated transmission 
class DustyGlassSpecularTransmission 
{
public:
	// DustyGlassSpecularTransmission Public Methods
	DustyGlassSpecularTransmission(float t_glass, float t_dust,
		float ior_air, float ior_dust, float ior_glass,
		float tau_dust)
		: T_glass(t_glass), T_dust(t_dust), 
		fresnel_air_dust(ior_air, ior_dust), 
		fresnel_dust_glass(ior_dust, ior_glass),
		Tau(tau_dust){}
		float Sample_f(const Vector &wo, Vector *wi, float u1, float u2, float *pdf) const;

private:
	// DustyGlassSpecularTransmission Private Data
	float T_glass, T_dust;
	FresnelDielectric fresnel_air_dust;
	FresnelDielectric fresnel_dust_glass;
	float Tau;
};

// L_r(3)+L_r(4): scattering components in reflection
class DustyGlassDiffuseReflection
{
public:
	// DustyGlassDiffuseReflection Public Methods
	DustyGlassDiffuseReflection(float r_glass, float t_glass, 
		float r_dust, float t_dust, 
		float ior_air, float ior_dust, float ior_glass, 
		float rho_dust, float tau_dust, float g_dust, int b_hk93)
		: R_glass(r_glass), T_glass(t_glass), R_dust(t_dust), T_dust(t_dust), 
		fresnel_air_dust(ior_air, ior_dust),
		fresnel_dust_glass(ior_dust, ior_glass), 
		rho(rho_dust), tau(tau_dust), g(g_dust), bHK93(b_hk93){
		}
		float f(const Vector &, const Vector &) const;

private:
	// DustyGlassDiffuseReflection Private Data
	float R_glass, T_glass, R_dust, T_dust;
	FresnelDielectric fresnel_air_dust;
	FresnelDielectric fresnel_dust_glass;
	float rho, tau, g;
	int bHK93;
};

// L_t(2): scattering component in transmission
class DustyGlassDiffuseTransmission
{
public:
	// DustyGlassDiffuseTransmission Public Methods
	DustyGlassDiffuseTransmission(float t_glass, float t_dust,
		float ior_air, float ior_dust, float ior_glass,
		float rho_dust, float tau_dust, float g_dust)
		: T_glass(t_glass), T_dust(t_dust), 
		fresnel_air_dust(ior_air, ior_dust),
		fresnel_dust_glass(ior_dust, ior_glass),
		rho(rho_dust), tau(tau_dust), g(g_dust){
		}
		float f(const Vector &, const Vector &) const;

private:
	// DustyGlassDiffuseTransmission Private Data
	float T_glass, T_dust;
	FresnelDielectric fresnel_air_dust;
	FresnelDielectric fresnel_dust_glass;
	float rho, tau, g;
};

// L_r(1)
float DustyGlassSpecularReflectionDust::Sample_f(const Vector &wo,
					Vector *wi, float u1, float u2, float *pdf) const 
{
	bool entering = CosTheta(wo) > 0.f;
	if(entering) // mirror reflection on the air-contaminant interface
	{
		// if no interface between air and contaminant, this term is zero
		if(fresnel_air_dust.eta_i == fresnel_air_dust.eta_t)
		{
			return (0.f);
		}

		*wi = Vector(-wo.x, -wo.y, wo.z);
		*pdf = 1.f;
		float cosThetai = fabsf(CosTheta(*wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;
		return R_dust * fresnel_air_dust.Evaluate(CosTheta(wo)) / cosThetai;
	}
	else	// mirror reflection on the contaminant-glass interface
	{
		// if no interface between glass and contaminant, this term is zero
		if(fresnel_dust_glass.eta_i == fresnel_dust_glass.eta_t)
		{
			return (0.f);
		}

		*wi = Vector(-wo.x, -wo.y, wo.z);		
		*pdf = 1.f;
		float cosThetai = fabsf(CosTheta(*wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;
		return R_glass * fresnel_dust_glass.Evaluate(CosTheta(wo)) / cosThetai;
	}
}

// L_r(2)
float DustyGlassSpecularReflection::Sample_f(const Vector &wo,
						Vector *wi, float u1, float u2, float *pdf) const 
{	
	bool entering = CosTheta(wo) > 0.f;
	if(entering)
	{
		// if no interface between glass and contaminant, this term is zero
		if(fresnel_dust_glass.eta_i == fresnel_dust_glass.eta_t)
		{
			return (0.f);
		}

		if(T_dust==0 || R_glass==0)
		{
			return (0.f);
		}

		float weight(1.f), F;

		float ei = fresnel_air_dust.eta_i, et = fresnel_air_dust.eta_t;
		float sini2 = SinTheta2(wo);
		float eta = ei / et;
		float sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;		
		float cost = sqrtf(std::max(0.f, 1.f - sint2));
		F = fresnel_air_dust.Evaluate(CosTheta(wo));
		float TT = T_dust * ((1.f)-F);
		weight *= TT*TT; // Fresnel (the two Fresnel transmissions should be the same)
		float sintOverSini = eta;
		Vector wo1 = Vector(sintOverSini * -wo.x,
			sintOverSini * -wo.y,
			-cost);

		float cosThetaO1 = fabsf(CosTheta(wo1));
		if(cosThetaO1<FLT_EPSILON) cosThetaO1=FLT_EPSILON;
		weight *= exp(-2*Tau/cosThetaO1); // dust attenuation
		F = fresnel_dust_glass.Evaluate(CosTheta(wo1)); // Fresnel 
		weight *= R_glass * F; 

		*wi = Vector(-wo.x, -wo.y, wo.z);
		*pdf = 1.f;		
		float cosThetai = fabsf(CosTheta(*wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;
		return weight / cosThetai;
	}
	else
	{	
		// if no interface between air and contaminant, this term is zero
		if(fresnel_air_dust.eta_i == fresnel_air_dust.eta_t)
		{
			return (0.f);
		}

		if(T_glass==0 || R_dust==0)
		{
			return (0.f);
		}

		float weight(1.f), F;

		float ei = fresnel_dust_glass.eta_i, et = fresnel_dust_glass.eta_t;
		float sini2 = SinTheta2(wo);
		float eta = et / ei;
		float sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;
		float cost = sqrtf(std::max(0.f, 1.f - sint2));
		float sintOverSini = eta;
		Vector wo1(sintOverSini * -wo.x,
			sintOverSini * -wo.y,
			cost);

		// Fresnel
		F = fresnel_dust_glass.Evaluate(CosTheta(wo));
		float TT = T_glass * ((1.f)-F); 
		weight *= TT * TT; // The two Fresnel transmissions are the same.
		float cosThetaO1 = fabsf(CosTheta(wo1));
		if(cosThetaO1<FLT_EPSILON) cosThetaO1=FLT_EPSILON;
		weight *= exp(-2*Tau/cosThetaO1); // dust attenuation
		F = fresnel_air_dust.Evaluate(CosTheta(wo1)); // Fresnel 
		weight *= R_dust * F;

		*wi = Vector(-wo.x, -wo.y, wo.z);
		*pdf = 1.f;
		float cosThetai = fabsf(CosTheta(*wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;
		return weight/cosThetai; 
	}
}

// L_t(1)
float DustyGlassSpecularTransmission::Sample_f(const Vector &wo,
					Vector *wi, float u1, float u2, float *pdf) const 
{
	bool entering = CosTheta(wo) > 0.f;
	if(entering)
	{
		float weight(1.f), F;

		float ei = fresnel_air_dust.eta_i, et = fresnel_air_dust.eta_t;
		float sini2 = SinTheta2(wo);
		float eta = ei / et;
		float sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;
		F = fresnel_air_dust.Evaluate(CosTheta(wo));
		weight *= T_dust * ((1.f)-F); // Fresnel		
		float cost = sqrtf(std::max(0.f, 1.f - sint2));
		float sintOverSini = eta;
		Vector wo1 = Vector(sintOverSini * -wo.x,
			sintOverSini * -wo.y,
			-cost);

		float cosThetao1 = fabsf(CosTheta(wo1));
		if(cosThetao1<FLT_EPSILON) cosThetao1=FLT_EPSILON;
		weight *= exp(-Tau/cosThetao1); // dust attenuation

		// change for contaminant-glass interface 
		wo1 = -wo1;		
		F = fresnel_dust_glass.Evaluate(CosTheta(wo1));
		weight *= T_glass * ((1.f) - F); // Fresnel		

		// get the wi
		ei = fresnel_dust_glass.eta_i; et = fresnel_dust_glass.eta_t;
		sini2 = SinTheta2(wo1);
		eta   = ei / et;
		sint2 = eta * eta * sini2;
		if(sint2 > 1.) return 0.f;
		cost = sqrtf(std::max(0.f, 1.f-sint2));
		sintOverSini = eta;
		*wi = Vector(sintOverSini * -wo1.x,
			sintOverSini * -wo1.y,
			-cost);
		*pdf = 1.f;
		float cosThetai = fabsf(CosTheta(*wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;

		// additional term in Fresnel
		ei = fresnel_air_dust.eta_i;
		et = fresnel_dust_glass.eta_t;
		weight *= (ei*ei)/(et*et);		

		return weight / cosThetai;
	}
	else
	{
		float weight(1.f), F;

		// get the wo inside the dust
		float ei = fresnel_dust_glass.eta_t, et = fresnel_dust_glass.eta_i;
		float sini2 = SinTheta2(wo);
		float eta = ei / et;
		float sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;
		F = fresnel_dust_glass.Evaluate(CosTheta(wo));
		weight *= T_glass * ((1.f)-F); // Fresnel

		float cost = sqrtf(std::max(0.f, 1.f - sint2));
		float sintOverSini = eta;
		Vector wo1 = Vector(sintOverSini * -wo.x,
			sintOverSini * -wo.y,
			cost);

		float cosThetao1 = fabsf(CosTheta(wo1));
		if(cosThetao1<FLT_EPSILON) cosThetao1=FLT_EPSILON;
		weight *= exp(-Tau/cosThetao1); // dust attenuation

		// change for air-contaminant interface 
		wo1 = -wo1;
		F = fresnel_air_dust.Evaluate(CosTheta(wo1));
		weight *= T_dust * ((1.f) - F); // Fresnel

		// get the wi
		ei = fresnel_air_dust.eta_t; et = fresnel_air_dust.eta_i;
		sini2 = SinTheta2(wo1);
		eta   = ei / et;
		sint2 = eta * eta * sini2;
		if(sint2 > 1.) return 0.f;
		cost = sqrtf(std::max(0.f, 1.f-sint2));
		sintOverSini = eta;
		*wi = Vector(sintOverSini * -wo1.x,
			sintOverSini * -wo1.y,
			cost);
		*pdf = 1.f;
		float cosThetai = fabsf(CosTheta(*wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;

		// additional term in Fresnel
		ei = fresnel_air_dust.eta_i;
		et = fresnel_dust_glass.eta_t;
		weight *= (et*et)/(ei*ei);		

		return weight / cosThetai;
	}
}

// L_r(3)+L_r(4)
float DustyGlassDiffuseReflection::f(const Vector &wo, const Vector &wi) const 
{
	bool entering = CosTheta(wo) > 0.f;
	if (entering)
	{
		if(T_dust==0)
			return (0.f);

		float weight(1.f), F;

		float ei = fresnel_air_dust.eta_i, et = fresnel_air_dust.eta_t;
		float sini2 = SinTheta2(wo);
		float eta = ei / et;
		float sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;
		float cost = sqrtf(std::max(0.f, 1.f - sint2));
		F = fresnel_air_dust.Evaluate(CosTheta(wo));
		weight *= T_dust * ((1.f)-F); // Fresnel
		float sintOverSini = eta;
		Vector wo1 = Vector(sintOverSini * -wo.x,
			sintOverSini * -wo.y,
			-cost);

		sini2 = SinTheta2(wi);
		sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;
		cost = sqrtf(max(0.f, 1.f-sint2));
		F = fresnel_air_dust.Evaluate(CosTheta(wi));
		weight *= T_dust * (float(1.f)-F); // Fresnel 1
		sintOverSini = eta;
		Vector wi1 = Vector(sintOverSini * -wi.x,
			sintOverSini * -wi.y,
			-cost);

		// now, actual dust attenuation and scattering computation
		// first change to the contaminant-glass interface
		wo1 = -wo1;
		wi1 = -wi1; 
		float cosThetaO = fabsf(CosTheta(wo1));
		if(cosThetaO<FLT_EPSILON) cosThetaO=FLT_EPSILON;
		float cosThetaI = fabsf(CosTheta(wi1));
		if(cosThetaI<FLT_EPSILON) cosThetaI=FLT_EPSILON;
		float A  = exp(-tau*(1.0f/cosThetaI+1.0f/cosThetaO));		
		float p  = PhaseHG(-wi1, wo1, g);
		float fd = rho*p/(cosThetaI+cosThetaO);
		float Lr3 = (1-A)*fd*cosThetaI*weight;

		// now compute L_r(4) (and also L_r(5), similar as L_r(4). see dustyglass.cpp)
		float Lr4(0.f);
		
		// Lr4 exists only if there is an interface 
		// between the contaminant and the glass
		if((fresnel_dust_glass.eta_i!=fresnel_dust_glass.eta_t) && !R_glass==0)
		{						
			float Br4, Br5, B1, B2, B;
			if(cosThetaI==cosThetaO)
			{
				Br4  = exp(-2*tau/cosThetaI)*tau/cosThetaI;
				Br5  = Br4;
			}
			else
			{
				B1  = exp(-tau/cosThetaI);
				B2  = exp(-tau/cosThetaO);
				B   = (B1-B2)/(cosThetaI-cosThetaO)*cosThetaI;
				Br4   = B*exp(-tau/cosThetaI);
				Br5   = B*exp(-tau/cosThetaO);
			}
			Vector wo1_reflected=Vector(-wo1.x, -wo1.y, wo1.z);
			p = PhaseHG(wi1, wo1_reflected, g); // the same for L_r(4) and L_r(5)
			float Fr4 = fresnel_dust_glass.Evaluate(CosTheta(wi1));
			float Fr5 = fresnel_dust_glass.Evaluate(CosTheta(wo1));
		//	Lr4 = (Br4*Fr4+Br5*Fr5)*weight*rho*p*R_glass;
			Lr4 = Br4*Fr4*weight*rho*p*R_glass;
		}
		
		if(bHK93==1) Lr4 = (0.f);

		// finally done
		float cosThetai = fabsf(CosTheta(wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;
		return (Lr3+Lr4)/cosThetai;
	}
	else
	{
		if(T_glass==0)
			return (0.f);

		float weight(1.f), F;

		float ei = fresnel_dust_glass.eta_i, et = fresnel_dust_glass.eta_t;

		// get the wo and wi inside the dust
		float sint2 = SinTheta2(wi);
		float eta = et / ei;
		float sini2 = eta * eta * sint2;
		if (sini2 > 1.) return 0.f;

		float sint2o = SinTheta2(wo);
		float sini2o = eta * eta * sint2o;
		if (sini2o > 1.) return 0.f;

		float cosi = sqrtf(std::max(0.f, 1.f - sini2));
		float siniOverSint = eta;
		Vector wio(siniOverSint * -wi.x,
			siniOverSint * -wi.y,
			cosi);
		F  = fresnel_dust_glass.Evaluate(CosTheta(wi));
		weight *= T_glass*((1.)-F);

		float cosio = sqrtf(std::max(0.f, 1.f - sini2o));
		float sinioOverSint = eta;
		Vector woi(sinioOverSint * -wo.x,
			sinioOverSint * -wo.y,
			cosio);
		F  = fresnel_dust_glass.Evaluate(CosTheta(wo));
		weight *= T_glass*((1.)-F);


		// Lr3: dust attenuation and scattering computation
		woi = -woi;
		wio = -wio; 
		float cosThetaO = fabsf(CosTheta(woi));
		if(cosThetaO<FLT_EPSILON) cosThetaO=FLT_EPSILON;
		float cosThetaI = fabsf(CosTheta(wio));
		if(cosThetaI<FLT_EPSILON) cosThetaI=FLT_EPSILON;
		float A  = exp(-tau*(1.0f/cosThetaI+1.0f/cosThetaO));
		float p  = PhaseHG(wio, -woi, g);
		float fd = rho*p/(cosThetaI+cosThetaO);
		float Lr3 = (1-A) * fd * cosThetaI * weight;

		// now compute Lr4 (and Lr5)
		float Lr4(0.f);

		// Lr4 exists only if there is an interface 
		// between the contaminant and the air
		if((fresnel_air_dust.eta_i!=fresnel_air_dust.eta_t) && !R_dust==0)
		{						
			float Br4, Br5, B1, B2, B;
			if(cosThetaI==cosThetaO)
			{
				Br4  = exp(-2*tau/cosThetaI)*tau/cosThetaI;
				Br5  = Br4;
			}
			else
			{
				B1  = exp(-tau/cosThetaI);
				B2  = exp(-tau/cosThetaO);
				B   = (B1-B2)/(cosThetaI-cosThetaO)*cosThetaI;
				Br4   = B*exp(-tau/cosThetaI);
				Br5   = B*exp(-tau/cosThetaO);
			}
			Vector woi_reflected=Vector(-woi.x, -woi.y, woi.z);
			p = PhaseHG(wio, woi_reflected, g); // the same for L_r(4) and L_r(5)
			float Fr4 = fresnel_air_dust.Evaluate(CosTheta(wio));
			float Fr5 = fresnel_air_dust.Evaluate(CosTheta(woi));
		//	Lr4 = (Br4*Fr4+Br5*Fr5)*weight*rho*p*R_dust;
			Lr4 = Br4*Fr4*weight*rho*p*R_dust;
		}

		if(bHK93==1) Lr4 = (0.f);

		// finally done
		float cosThetai = fabsf(CosTheta(wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;

		return (Lr3+Lr4)/cosThetai;
	}
}

// L_t(2)
float DustyGlassDiffuseTransmission::f(const Vector &wo, const Vector &wi) const 
{
	bool entering = CosTheta(wo) > 0.f;
	if (entering)
	{
		float weight(1.f), F;

		// get the wi inside the dust
		float ei = fresnel_dust_glass.eta_i, et = fresnel_dust_glass.eta_t;
		float sint2 = SinTheta2(wi);
		float eta = et / ei;
		float sini2 = eta * eta * sint2;
		if (sini2 > 1.) return 0.f;
		float cosi = sqrtf(std::max(0.f, 1.f - sini2));
		float siniOverSint = eta;
		Vector wio(siniOverSint * -wi.x,
			siniOverSint * -wi.y,
			cosi);
		F  = fresnel_dust_glass.Evaluate(CosTheta(wi));
		weight *= T_glass * ((1.f)-F); // Fresnel


		// get the wo inside the dust
		ei = fresnel_air_dust.eta_i, et = fresnel_air_dust.eta_t;
		sini2 = SinTheta2(wo);
		eta = ei / et;
		sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;
		float cost = sqrtf(std::max(0.f, 1.f - sint2));		
		float sintOverSini = eta;
		Vector wo1 = Vector(sintOverSini * -wo.x,
			sintOverSini * -wo.y,
			-cost);
		F = fresnel_air_dust.Evaluate(CosTheta(wo));
		weight *= T_dust * ((1.f)-F); // Fresnel

		// now, actual dust attenuation and scattering computation
		wo1 = -wo1;	
		float    p  = PhaseHG(wio,wo1,g);
		float	 fd = rho*p;
		float    A1, A2, A;
		float cosThetaO = fabsf(CosTheta(wo1));
		if(cosThetaO<FLT_EPSILON) cosThetaO=FLT_EPSILON;
		float cosThetaI = fabsf(CosTheta(wio));
		if(cosThetaI<FLT_EPSILON) cosThetaI=FLT_EPSILON;
		float cosThetai = fabsf(CosTheta(wi));
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;

		if(cosThetaI==cosThetaO)
		{
			A  = exp(-tau/cosThetaI)*tau/cosThetaI;
		}
		else
		{
			A1  = exp(-tau/cosThetaI);
			A2  = exp(-tau/cosThetaO);
			A   = (A1-A2)/(cosThetaI-cosThetaO)*cosThetaI;
		}

		ei = fresnel_air_dust.eta_i; 
		et = fresnel_dust_glass.eta_t;
		weight *= (ei*ei)/(et*et);

		return  weight * fd * A / cosThetai;
	}	
	else
	{
		float weight(1.f), F;

		// get the wo inside the dust
		float ei = fresnel_dust_glass.eta_i, et = fresnel_dust_glass.eta_t;
		float sini2 = SinTheta2(wo);
		float eta = et / ei;
		float sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;
		float cost = sqrtf(std::max(0.f, 1.f - sint2));
		float sintOverSini = eta;
		Vector woi(sintOverSini * -wo.x,
			sintOverSini * -wo.y,
			cost);
		F = fresnel_dust_glass.Evaluate(CosTheta(wo));
		weight *= T_glass * ((1.f)-F); // Fresnel

		// get the wi inside the dust
		ei = fresnel_air_dust.eta_i; et = fresnel_air_dust.eta_t;
		sini2 = SinTheta2(wi);
		eta = ei / et;
		sint2 = eta * eta * sini2;
		if (sint2 > 1.) return 0.f;
		cost = sqrtf(std::max(0.f, 1.f-sint2));
		sintOverSini = eta;
		Vector wi1 = Vector(sintOverSini * -wi.x,
			sintOverSini * -wi.y,
			-cost);
		F = fresnel_air_dust.Evaluate(CosTheta(wi));
		weight *= T_dust * ((1.f)-F); // Fresnel 1


		// now actual dust attenuation and scattering computation
		woi = -woi;
		float    p  = PhaseHG(woi,wi1,g);
		float	 fd = rho*p;
		float    A1, A2, A;
		float cosThetaI = fabsf(CosTheta(wi1));
		if(cosThetaI<FLT_EPSILON) cosThetaI=FLT_EPSILON;
		float cosThetaO = fabsf(CosTheta(woi));
		if(cosThetaO<FLT_EPSILON) cosThetaO=FLT_EPSILON;
		if(cosThetaI==cosThetaO)
		{
			A  = exp(-tau/cosThetaI)*tau/cosThetaI;
		}
		else
		{
			A1  = exp(-tau/cosThetaI);
			A2  = exp(-tau/cosThetaO);
			A   = (A1-A2)/(cosThetaI-cosThetaO)*cosThetaI;
		}

		float cosThetai = fabsf(CosTheta(wi)); 
		if(cosThetai<FLT_EPSILON) cosThetai=FLT_EPSILON;

		ei = fresnel_air_dust.eta_i; 
		et = fresnel_dust_glass.eta_t;
		weight *= (et*et)/(ei*ei);
		return  weight * fd * A / cosThetai;
	}
}

// DustyGlassOil Class Declarations
//
// add a surface (interface) for the dust layer
// one extra parameter compared with DustyGlass:
//	1. the index of the dust layer
//
class DustyGlassOil 
{
public:
	// DustyGlassOil Public Methods
	DustyGlassOil(float r_glass, float t_glass,
		float i_glass,  float i_air, float i_dust, 
		float r_dust, float t_dust, float rho_dust, 
		float tau_dust, float g_dust, int bHK93)
	{
		Kr_glass = r_glass;
		Kt_glass = t_glass;
		index_glass = i_glass;
		index_air = i_air;
		index_dust = i_dust;
		Kr_dust = r_dust;
		Kt_dust = t_dust;
		Rho_dust = rho_dust;
		Tau_dust = tau_dust;
		G_dust   = g_dust;	
		B_HK93	 = bHK93;

		m_Lr1=NULL;
		m_Lr2=NULL;
		m_Lr345=NULL;
		m_Lt1=NULL;
		m_Lt2=NULL;
	}

	void GetBSDF();
	void ReleaseBSDF();
	float my_f(const Vector &wo, const Vector &wi);

	float SpecularReflection(const Vector &wo, Vector *wi) ;
	float SpecularTransmission(const Vector &wo, Vector *wi) ;

	bool HasSpecularComponents()
	{
		return true;
	}
private:
	// Glass Private Data
	float Kr_glass, Kt_glass, Kr_dust, Kt_dust, Rho_dust;
	float index_glass, index_air, index_dust, Tau_dust, G_dust;
	int B_HK93;

	DustyGlassSpecularReflectionDust*	m_Lr1;
	DustyGlassSpecularReflection*		m_Lr2;
	DustyGlassSpecularTransmission*		m_Lt1;
	DustyGlassDiffuseReflection*		m_Lr345;
	DustyGlassDiffuseTransmission*		m_Lt2;
};

// Glass Method Definitions
void DustyGlassOil::GetBSDF() 
{
	// Allocate _BSDF_, possibly doing bump-mapping with _bumpMap_
	float R_glass = Kr_glass;
	float T_glass = Kt_glass;
	float ior_glass = index_glass;
	float ior_air = index_air;

	float ior_dust = index_dust;
	float R_dust = Kr_dust;
	float T_dust = Kt_dust;
	float rho_dust = Rho_dust;
	float tau_dust    = Tau_dust;
	float g_dust      = G_dust;
	
	// Lr(1) specular reflection on contaminant or glass (1 specular reflection
	m_Lr1 = new DustyGlassSpecularReflectionDust(R_glass, R_dust, 
								ior_air, ior_dust, ior_glass);

	// Lr(2) attenuation parts (2 specular transmission, 1 specular reflection)
	m_Lr2 = new DustyGlassSpecularReflection(R_glass, 
								T_glass, R_dust, T_dust,
								ior_air, ior_dust, ior_glass, tau_dust);

	// Lt(1) 2 specular transmissions
	m_Lt1 = new DustyGlassSpecularTransmission(T_glass, T_dust, 
								ior_air, ior_dust, ior_glass, tau_dust);

	// single scattering parts 
	// diffuse reflection (Lr3, Lr4, Lr5)
	m_Lr345 = new DustyGlassDiffuseReflection(R_glass, 
								T_glass, R_dust, T_dust, 
								ior_air, ior_dust, ior_glass, rho_dust, tau_dust, g_dust, B_HK93);

	// diffuse transmission (Lt2)
	m_Lt2 = new DustyGlassDiffuseTransmission(T_glass, T_dust,
								ior_air, ior_dust, ior_glass, rho_dust, tau_dust, g_dust);
	
}

void DustyGlassOil::ReleaseBSDF()
{
	delete m_Lr1;
	delete m_Lr2;
	delete m_Lr345;
	delete m_Lt1;
	delete m_Lt2;
	m_Lr1=NULL;
	m_Lr2=NULL;
	m_Lr345=NULL;
	m_Lt1=NULL;
	m_Lt2=NULL;
}

float DustyGlassOil::my_f(const Vector &wo, const Vector &wi)
{
	Vector nn = Vector(0,0,1);
	if (Dot(wi, nn) * Dot(wo, nn) > 0) 
		return m_Lr345->f(wo,wi);
	else
		return m_Lt2->f(wo,wi);
}

float DustyGlassOil::SpecularReflection(const Vector &wo, Vector *wi) 
{
	Vector wii;
	float pdf;
	// This function will be called as (wi, *wo). So we first use this to get
	// the outgoing direction (wo). But then when computing the BRDF, we should
	// switch back because our model is not symmetric. 
	float f = m_Lr1->Sample_f(wo, wi, RandomFloat(), RandomFloat(), &pdf);
	if(f==0)
	{
		f = m_Lr2->Sample_f(wo, wi, RandomFloat(), RandomFloat(), &pdf);
	}

	f  = m_Lr1->Sample_f(*wi, &wii, RandomFloat(), RandomFloat(), &pdf);
	f += m_Lr2->Sample_f(*wi, &wii, RandomFloat(), RandomFloat(), &pdf);
	return f;
}

float DustyGlassOil::SpecularTransmission(const Vector &wo, Vector *wi) 
{
	Vector wii;
	float pdf;
	float f = m_Lt1->Sample_f(wo, wi, RandomFloat(), RandomFloat(), &pdf);
	if(f!=0)
		f = m_Lt1->Sample_f(*wi, &wii, RandomFloat(), RandomFloat(), &pdf);
	return f;
}

void printHelp()
{
	cout << "Tool for evaluation of the BRDF/BTDF model " << endl;
	cout << "for dirty glass (glass covered with a layer of scattering media) "<< endl;
	cout << "The illumination is collimated light (from top or from bottom)" << endl;
	cout << "" << endl;
	cout << "Usage: evaluator.exe dust_tau dust_g dust_albedo dust_ior theta_i ngrids outputfile [bHK93=0]" << endl;
	cout << "" << endl;
	cout << "dust_tau\t optical thickness of the contaminant. "<< endl;
	cout << "        \t typical in (0,0.4) for good approximation" << endl;
	cout << "dust_g\t\t scattering parameter of the contaminant. in (-1,1)" << endl;
	cout << "dust_albedo\t albedo of single scattering. in (0,1). " << endl;
	cout << "           \t typical less than 0.5 for good approximation" << endl;
	cout << "dust_ior\t refraction index of the contaminant." << endl;
	cout << "theta_i(wi)\t direction of the illumination. in (0,pi). " << endl;
	cout << "           \t if less than pi/2, air2glass; otherwise glass2air " << endl;
	cout << "ngrids\t\t sampling points. wo will be span in (0,pi/2) and (pi/2, pi) " << endl;
	cout << "      \t\t into ngrids" << endl;
	cout << "" << endl;
	cout << "Example: " << endl;
	cout << "\tevaluator.exe 0.4 0.7 0.5 1.0 0 200 result.txt" << endl;
	cout << "" << endl;
	cout << "Format of the output: " << endl;
	cout << "" << endl;
	cout << "There are 4 parts in the output file." << endl; 
	cout << "The values are reflectance or transmittance" << endl;
	cout << "(i.e. physically corresponds to radiance (W/m2/sr)"  << endl;
	cout << "So they are BRDF*cos(theta_i) or BTDF*cos(theta_i) " << endl;
	cout << "" << endl;
	cout << "  Specular Reflectance (angle, value)" << endl;
	cout << "  Diffuse Reflectance  (angles, values)" << endl;
	cout << "  Specular Transmission (angle, value)" << endl;
	cout << "  Diffuse Transmission (angles values)" << endl;
	cout << "" << endl;
	cout << "Jinwei GU (jwgu@cs.columbia.edu). 2007/6/5" << endl;
}

int main(int argc, char* argv[])
{
	float Kr_glass = float(1.f);
	float Kt_glass = float(1.f);
	float Kr_dust  = float(1.f);
	float Kt_dust  = float(1.f);

	float index_glass = 1.5f;
	float index_air   = 1.f;

	if(argc<8)
	{
		printHelp();
		return 0;
	}
	
	srand(time(0));

	float tau_dust    = Clamp(atof(argv[1]), 0.f, FLT_MAX);
	float g_dust      = Clamp(atof(argv[2]), -1.f, 1.f);
	float albedo_dust = Clamp(atof(argv[3]), 0.f, 1.f);
	float rho_dust = float(albedo_dust);
	float index_dust  = Clamp(atof(argv[4]), 0.f, FLT_MAX);
	float theta_i     = Clamp(atof(argv[5]), 0.f, M_PI);
	int   ngrids	  = (int)Clamp(atoi(argv[6]), 1, INT_MAX);
	string outputfile = string(argv[7]);
	FILE* fp = fopen(outputfile.c_str(), "wt");
	if(fp==NULL)
	{
		cout << "***Error: can not open file " << outputfile << endl;
		return 0;
	}
	
	int bHK93 = 0;
	if(argc>8) bHK93 = atoi(argv[8]);

	DustyGlassOil* dustyglass_oil = new DustyGlassOil(Kr_glass, Kt_glass, 
									index_glass, index_air, index_dust, 
									Kr_dust, Kt_dust, rho_dust, tau_dust, g_dust, bHK93);
	dustyglass_oil->GetBSDF();
	
	float sini = sin(theta_i);
	float cosi = cos(theta_i);
	bool entering = cosi > 0.f;
	Vector wi = Vector(sini,0,cosi);
	Vector wo;
	float f;

	// Part 1: compute the specular reflectance
	f = dustyglass_oil->SpecularReflection(wi, &wo);
	f = f*fabs(cosi);
	fprintf(fp, "SpecularReflectance\n");
	fprintf(fp, "wo %f %f %f\n", wo[0], wo[1], wo[2]);
	fprintf(fp, "R  %f\n\n", f);
	
	// Part 2: compute the specular transmittance
	f = dustyglass_oil->SpecularTransmission(wi, &wo);
	f = f*fabs(cosi);
	fprintf(fp, "SpecularTransmittance\n");
	fprintf(fp, "wo %f %f %f\n", wo[0], wo[1], wo[2]);
	fprintf(fp, "T  %f\n\n", f);

	// Part 3: compute the diffuse reflecance
	fprintf(fp, "DiffuseReflectance\n");
	for(int i=0; i<ngrids; i++)
	{
		float thetao = M_PI/2.0f-(0.5f+i)*M_PI/ngrids;
		if(entering)
			wo = Vector(sin(thetao), 0, cos(thetao));
		else
			wo = Vector(sin(thetao), 0, -cos(thetao));

		float f = dustyglass_oil->my_f(wo, wi);
		f = f*fabs(cosi);
		fprintf(fp, "%f\t%f\n", thetao, f);
	}
	fprintf(fp, "\n");


	// Part 4: compute the diffuse transmittance
	fprintf(fp, "DiffuseTransmittance\n");
	for(int i=0; i<ngrids; i++)
	{
		float thetao = M_PI/2.0f-(0.5f+i)*M_PI/ngrids;
		if(entering)
			wo = Vector(sin(thetao), 0, -cos(thetao));
		else
			wo = Vector(sin(thetao), 0, cos(thetao));

		float f = dustyglass_oil->my_f(wo, wi);
		f = f*fabs(cosi);
		fprintf(fp, "%f\t%f\n", thetao, f);
	}
	fprintf(fp, "\n");
	
	dustyglass_oil->ReleaseBSDF();
	delete dustyglass_oil;
	dustyglass_oil = NULL;

	fclose(fp); fp = NULL;
}