const float LUT_SIZE  = 64.0;
const float LUT_SCALE = (LUT_SIZE - 1.0)/LUT_SIZE;
const float LUT_BIAS  = 0.5/LUT_SIZE;
const float PI = 3.14159265;

vec3 mul(mat3 m, vec3 v)
{
	return m * v;
}

mat3 mul(mat3 m1, mat3 m2)
{
	return m1 * m2;
}

vec3 rotation_y(vec3 v, float a)
{
	vec3 r;
	r.x =  v.x*cos(a) + v.z*sin(a);
	r.y =  v.y;
	r.z = -v.x*sin(a) + v.z*cos(a);
	return r;
}

vec3 rotation_z(vec3 v, float a)
{
	vec3 r;
	r.x =  v.x*cos(a) - v.y*sin(a);
	r.y =  v.x*sin(a) + v.y*cos(a);
	r.z =  v.z;
	return r;
}

vec3 rotation_yz(vec3 v, float ay, float az)
{
	return rotation_z(rotation_y(v, ay), az);
}

vec3 rotation_yz_inv(vec3 v, float ay, float az)
{
	return rotation_y(rotation_z(v, -az), -ay);
}

float Fpo(float d, float l)
{
	return l/(d*(d*d + l*l)) + atan(l/d)/(d*d);
}

float Fwt(float d, float l)
{
	return l*l/(d*(d*d + l*l));
}

float ltc_R = 0.15;
vec3 cylinderCenter(in LightProperties light)  {return light.position.xyz;}
vec3 cylinderTangent(in LightProperties light) {return light.up.xyz;}
vec3 cylinderP1(in LightProperties light)      {return cylinderCenter(light) - 0.5 * light.dimensions.x * cylinderTangent(light);}
vec3 cylinderP2(in LightProperties light)      {return cylinderCenter(light) + 0.5 * light.dimensions.x * cylinderTangent(light);}

mat3 ltc_Minv;
#if 0
float D(vec3 w)
{
	vec3 wo = ltc_Minv * w;
	float lo = length(wo);
	float res = 1.0/PI * max(0.0, wo.z/lo) * abs(determinant(ltc_Minv)) / (lo*lo*lo);
	return res;
}

float I_ltc_disks(vec3 p1, vec3 p2, float R)
{
	float A = PI * R * R;
	vec3 wt  = normalize(p2 - p1);
	vec3 wp1 = normalize(p1);
	vec3 wp2 = normalize(p2);
	float Idisks = A * (
	D(wp1) * max(0.0, dot(+wt, wp1)) / dot(p1, p1) +
	D(wp2) * max(0.0, dot(-wt, wp2)) / dot(p2, p2));
	return Idisks;
}
#endif


float I_diffuse_line(vec3 p1, vec3 p2)
{
	// tangent
	vec3 wt = normalize(p2 - p1);

	// clamping
	if (p1.z <= 0.0 && p2.z <= 0.0) return 0.0;
	if (p1.z < 0.0) p1 = (+p1*p2.z - p2*p1.z) / (+p2.z - p1.z);
	if (p2.z < 0.0) p2 = (-p1*p2.z + p2*p1.z) / (-p2.z + p1.z);

	// parameterization
	float l1 = dot(p1, wt);
	float l2 = dot(p2, wt);

	// shading point orthonormal projection on the line
	vec3 po = p1 - l1*wt;

	// distance to line
	float d = length(po);

	// integral
	float I = (Fpo(d, l2) - Fpo(d, l1)) * po.z +
			  (Fwt(d, l2) - Fwt(d, l1)) * wt.z;
	return I / PI;
}

float I_ltc_line(vec3 p1, vec3 p2)
{
	// transform to diffuse configuration
	vec3 p1o = ltc_Minv * p1;
	vec3 p2o = ltc_Minv * p2;
	float I_diffuse = I_diffuse_line(p1o, p2o);

	// width factor
	vec3 ortho = normalize(cross(p1, p2));
	float w =  1.0 / length(inverse(transpose(ltc_Minv)) * ortho);

	return w * I_diffuse;
}


// LTC Support (N - surface normal, V - ray direction from the camera, P - shaded point position)
vec3 LTC_Evaluate(in LightProperties light, vec3 N, vec3 V, vec3 P)
{
	bool endCaps = true;

	// construct orthonormal basis around N
	vec3 T1, T2;
	T1 = normalize(V - N*dot(V, N));
	T2 = cross(N, T1);

	mat3 B = transpose(mat3(T1, T2, N));

	vec3 p1 = mul(B, cylinderP1(light) - P);
	vec3 p2 = mul(B, cylinderP2(light) - P);

	//if (analytic) // analytic integration
	{
		float Iline = ltc_R * I_ltc_line(p1, p2);
		//float Idisks = endCaps ? I_ltc_disks(p1, p2, R) : 0.0;
		float Idisks = 0.0; // endCaps ? I_ltc_disks(p1, p2, R) : 0.0;
		return vec3(min(1.0, Iline + Idisks));
	}
	//else // numerical integration
	//{
	//	float Icylinder = I_cylinder_numerical(p1, p2, R);
	//	float Idisks = endCaps ? I_disks_numerical(p1, p2, R) : 0.0;
	//	return vec3(Icylinder + Idisks);
	//}
}
