// vec3 constructNormal(float depth, vec2 texcoords, sampler2D depthtex) {
//     const vec2 offset1 = vec2(0.0,0.0001);
//     const vec2 offset2 = vec2(0.0001,0.0);
  
//     float depth1 = texture2D(depthtex, texcoords + offset1).r;
//     float depth2 = texture2D(depthtex, texcoords + offset2).r;
  
//     vec3 p1 = vec3(offset1, depth1 - depth);
//     vec3 p2 = vec3(offset2, depth2 - depth);
  
//     vec3 normal = cross(p1, p2);
//     normal.z = -normal.z;
  
//     normal = normalize(normal);
//     return vec3(normal.x, normal.z, normal.y);
// }

vec2 hammersley(uint i, float numSamples) {
    uint bits = i;
    bits = (bits << 16) | (bits >> 16);
    bits = ((bits & 0x55555555) << 1) | ((bits & 0xAAAAAAAA) >> 1);
    bits = ((bits & 0x33333333) << 2) | ((bits & 0xCCCCCCCC) >> 2);
    bits = ((bits & 0x0F0F0F0F) << 4) | ((bits & 0xF0F0F0F0) >> 4);
    bits = ((bits & 0x00FF00FF) << 8) | ((bits & 0xFF00FF00) >> 8);
    return vec2(i / numSamples, bits / exp2(32));
}


float map(float n, float start1, float stop1, float start2, float stop2) {
    return ((n - start1) / (stop1 - start1)) * (stop2 - start2) + start2;
}

vec3 blackbody(float t) {//return vec3(1.0, 0.949, 0.937);
    // http://en.wikipedia.org/wiki/Planckian_locus

    vec4 vx = vec4(-0.2661239e9,-0.2343580e6,0.8776956e3,0.179910);
    vec4 vy = vec4(-1.1063814,-1.34811020,2.18555832,-0.20219683);
    float it = 1./t;
    float it2= it*it;
    float x = dot(vx,vec4(it*it2,it2,it,1.));
    float x2 = x*x;
    float y = dot(vy,vec4(x*x2,x2,x,1.));
    float z = 1. - x - y;

    // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
    mat3 xyzToSrgb = mat3(
    3.2404542,-1.5371385,-0.4985314,
    -0.9692660, 1.8760108, 0.0415560,
    0.0556434,-0.2040259, 1.0572252
    );

    vec3 srgb = vec3(x/y,1.,z/y) * xyzToSrgb;
    return max(srgb,0.);
}

float signum(float x) {
	return uintBitsToFloat((floatBitsToUint(x) & 0x80000000u) | floatBitsToUint(1.0));
}

mat3 getRotMat(vec3 from, vec3 to) {
	float cosine = dot(from, to);

	float tmp = signum(cosine);
	    tmp = 1.0 / (tmp + cosine);

	vec3 axis = cross(to, from);
	vec3 tmpv = axis * tmp;

	return mat3(
		axis.x * tmpv.x + cosine, axis.x * tmpv.y - axis.z, axis.x * tmpv.z + axis.y,
		axis.y * tmpv.x + axis.z, axis.y * tmpv.y + cosine, axis.y * tmpv.z - axis.x,
		axis.z * tmpv.x - axis.y, axis.z * tmpv.y + axis.x, axis.z * tmpv.z + cosine
	);
}


// Returns the Spherical Harmonics base coefficients for the 0th and 1st order
// vec4 getSHBASE(vec3 normal) {
//     float Y00     = 0.282095;
//     float Y11     = 0.488603 * normal.x;
//     float Y10     = 0.488603 * normal.z;
//     float Y1_1    = 0.488603 * normal.y;
//     return vec4(Y00, Y11, Y10, Y1_1);
// }

// float computeProjectedSHFactor(vec4 origSHCoeff, vec3 normal) {
//     vec4 newBase = getSHBASE(normal);
//     vec3 A = vec3(3.1415926, 2.094395, 0.785398);
//     return A.x*newBase.x*origSHCoeff.x + A.y*newBase.w*origSHCoeff.w + A.y*newBase.z*origSHCoeff.z + A.y*newBase.y*origSHCoeff.y;
// }

// struct SH
// {
//     vec4 shY;
//     vec2 CoCg;
// };

// SH initSH() {
//     SH result;
//     result.shY = vec4(0);
//     result.CoCg = vec2(0);
//     return result;
// }

// vec3 projectSHIrradiance(SH sh, vec3 N) {
//     float d = dot(sh.shY.xyz, N);
//     float Y = 2.0 * (1.023326 * d + 0.886226 * sh.shY.w);
//     Y = max(Y, 0.0);

//     sh.CoCg *= Y * 0.282095 / (sh.shY.w + 1e-6);

//     float   T       = Y - sh.CoCg.y * 0.5;
//     float   G       = sh.CoCg.y + T;
//     float   B       = T - sh.CoCg.x * 0.5;
//     float   R       = B + sh.CoCg.x;

//     return max(vec3(R, G, B), vec3(0.0));
// }

// SH irradianceToSH(vec3 color, vec3 dir) {
//     SH result;

//     float   Co      = color.r - color.b;
//     float   t       = color.b + Co * 0.5;
//     float   Cg      = color.g - t;
//     float   Y       = max(t + Cg * 0.5, 0.0);

//     result.CoCg = vec2(Co, Cg);

//     float   L00     = 0.282095;
//     float   L1_1    = 0.488603 * dir.y;
//     float   L10     = 0.488603 * dir.z;
//     float   L11     = 0.488603 * dir.x;

//     result.shY = vec4 (L11, L1_1, L10, L00) * Y;

//     return result;
// }

// vec3 SHToIrradiance(SH sh) {
//     float   Y       = sh.shY.w / 0.282095;

//     float   T       = Y - sh.CoCg.y * 0.5;
//     float   G       = sh.CoCg.y + T;
//     float   B       = T - sh.CoCg.x * 0.5;
//     float   R       = B + sh.CoCg.x;

//     return max(vec3(R, G, B), vec3(0.0));
// }

// void accumulateSH(inout SH accum, SH new, float scale) {
//     accum.shY += new.shY * scale;
//     accum.CoCg += new.CoCg * scale;
// }

// SH mixSH(SH a, SH b, float alpha) {
//     SH result;
//     result.shY = mix(a.shY, b.shY, vec4(alpha));
//     result.CoCg = mix(a.CoCg, b.CoCg, vec2(alpha));
//     return result;
// }

struct SH
{
    vec3 L00;
    vec3 L11;
    vec3 L10;
    vec3 L1_1;
    vec3 L21;
    vec3 L2_1;
    vec3 L2_2;
    vec3 L20;
    vec3 L22;
};

SH initSH() {
    SH result;
    result.L00 = vec3(0.0);
    result.L11 = vec3(0.0);
    result.L10 = vec3(0.0);
    result.L1_1 = vec3(0.0);
    result.L21 = vec3(0.0);
    result.L2_1 = vec3(0.0);
    result.L2_2 = vec3(0.0);
    result.L20 = vec3(0.0);
    result.L22 = vec3(0.0);
    return result;
}

vec3 projectSHIrradiance(SH sh, vec3 normal) {
    //Multipliers for convolution with a cosine lobe
    float A0 = 3.141593;
    float A1 = 2.094395;
    float A2 = 0.785398;

    // Re-compute basis for projecting normal
    float Y00     = 0.282095;
    float Y11     = 0.488603 * normal.x;
    float Y10     = 0.488603 * normal.z;
    float Y1_1    = 0.488603 * normal.y;
    float Y21     = 1.092548 * normal.x*normal.z;
    float Y2_1    = 1.092548 * normal.y*normal.z;
    float Y2_2    = 1.092548 * normal.y*normal.x;
    float Y20     = 0.946176 * normal.z * normal.z - 0.315392;
    float Y22     = 0.546274 * (normal.x*normal.x - normal.y*normal.y);

    return     A0*Y00*sh.L00
             + A1*Y1_1*sh.L1_1 + A1*Y10*sh.L10 + A1*Y11*sh.L11 
             + A2*Y2_2*sh.L2_2 + A2*Y2_1*sh.L2_1 + A2*Y20*sh.L20 + A2*Y21*sh.L21 + A2*Y22*sh.L22;
}

vec3 projectSHIrradianceBentCone(SH sh, vec3 normal, float ao) {
    float t = 1.0 - ao;
    float t2 = t*t;
    float t3 = t2*t;
    float t4 = t2*t2;

    //Multipliers for convolution with a cosine lobe
    float A0 = 3.141593 * (1 - t2);
    float A1 = 2.094395 * (1 - t3);
    float A2 = 0.785398 * (3*(1 - t4) - 2*(1 - t2));

    // Re-compute basis for projecting normal
    float Y00     = 0.282095;
    float Y11     = 0.488603 * normal.x;
    float Y10     = 0.488603 * normal.z;
    float Y1_1    = 0.488603 * normal.y;
    float Y21     = 1.092548 * normal.x*normal.z;
    float Y2_1    = 1.092548 * normal.y*normal.z;
    float Y2_2    = 1.092548 * normal.y*normal.x;
    float Y20     = 0.946176 * normal.z * normal.z - 0.315392;
    float Y22     = 0.546274 * (normal.x*normal.x - normal.y*normal.y);

    return     A0*Y00*sh.L00
             + A1*Y1_1*sh.L1_1 + A1*Y10*sh.L10 + A1*Y11*sh.L11 
             + A2*Y2_2*sh.L2_2 + A2*Y2_1*sh.L2_1 + A2*Y20*sh.L20 + A2*Y21*sh.L21 + A2*Y22*sh.L22;
}

SH irradianceToSH(vec3 color, vec3 dir) {
    SH result;

    float Y00     = 0.282095;
    float Y11     = 0.488603 * dir.x;
    float Y10     = 0.488603 * dir.z;
    float Y1_1    = 0.488603 * dir.y;
    float Y21     = 1.092548 * dir.x*dir.z;
    float Y2_1    = 1.092548 * dir.y*dir.z;
    float Y2_2    = 1.092548 * dir.y*dir.x;
    float Y20     = 0.946176 * dir.z * dir.z - 0.315392;
    float Y22     = 0.546274 * (dir.x*dir.x - dir.y*dir.y);

    result.L00   = color * Y00;
    result.L11   = color * Y11;
    result.L10   = color * Y10;
    result.L1_1  = color * Y1_1;
    result.L21   = color * Y21;
    result.L2_1  = color * Y2_1;
    result.L2_2  = color * Y2_2;
    result.L20   = color * Y20;
    result.L22   = color * Y22;

    return result;
}

// vec3 SHToIrradiance(SH sh) {
//     return vec3( sh.red.w, sh.green.w, sh.blue.w ) / 0.282095;
// }

void accumulateSH(inout SH accum, SH new, float scale) {
    accum.L00 += new.L00 * scale;
    accum.L11 += new.L11 * scale;
    accum.L10 += new.L10 * scale;
    accum.L1_1 += new.L1_1 * scale;
    accum.L21 += new.L21 * scale;
    accum.L2_1 += new.L2_1 * scale;
    accum.L2_2 += new.L2_2 * scale;
    accum.L20 += new.L20 * scale;
    accum.L22 += new.L22 * scale;
}

// SH mixSH(SH a, SH b, float alpha) {
//     SH result;
//     result.red = mix(a.red, b.red, vec4(alpha));
//     result.green = mix(a.green, b.green, vec4(alpha));
//     result.blue = mix(a.blue, b.blue, vec4(alpha));
//     return result;
// }

// void divideSH(inout SH a, float b) {
//     a.red /= b;
//     a.green /= b;
//     a.blue /= b;
// }

// SH minSH(SH a, SH b) {
//     SH result;
//     result.red = min(a.red, b.red);
//     result.green = min(a.green, b.green);
//     result.blue = min(a.blue, b.blue);
//     return result;
// }

// SH maxSH(SH a, SH b) {
//     SH result;
//     result.red = max(a.red, b.red);
//     result.green = max(a.green, b.green);
//     result.blue = max(a.blue, b.blue);
//     return result;
// }

SH scaleSH(SH a, float b) {
    SH result = a;
    result.L00 *= b;
    result.L11 *= b;
    result.L10 *= b;
    result.L1_1 *= b;
    result.L21 *= b;
    result.L2_1 *= b;
    result.L2_2 *= b;
    result.L20 *= b;
    result.L22 *= b;
    return result;
}


//// LOBE AWARE FILTERING ////

// SPHERICAL GAUSSIAN FUNCTION
float SG(vec3 dir1, vec3 dir2, float lambda) {
    return exp(lambda*( clamp(dot(dir1, dir2), 0.0, 1.0) - 1));
}

float lobe_sharpness(float roughness) {
    return 2.0 / (roughness * roughness);
}

float lobe_weight(vec3 centerNorm, vec3 norm, float centerRoughness, float roughness, float beta) {
    float lamba_i = lobe_sharpness(centerRoughness);
    float lamba_j = lobe_sharpness(roughness);
    return pow(2*sqrt(lamba_i*lamba_j)/(lamba_i + lamba_j), beta) * SG(centerNorm, norm, beta*lamba_i*lamba_j/(lamba_i + lamba_j));
}