// 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
{
    vec4 red;
    vec4 green;
    vec4 blue;
};

SH initSH() {
    SH result;
    result.red = vec4(0);
    result.green = vec4(0);
    result.blue = vec4(0);
    return result;
}

vec3 projectSHIrradiance(SH sh, vec3 N) {
    //Multipliers for convolution with a cosine lobe
    vec2 A        = vec2(3.141593, 2.094395);
    //  vec2 A        = vec2(3.141593, 2.094395);
    // Re-compute basis for projecting normal
    float Y00     = 0.282095;
    float Y1_1    = 0.488603 * N.y;
    float Y10     = 0.488603 * N.z;
    float Y11     = 0.488603 * N.x;

    vec3 L00      = vec3(sh.red.w, sh.green.w, sh.blue.w);
    vec3 L1_1     = vec3(sh.red.y, sh.green.y, sh.blue.y);
    vec3 L10      = vec3(sh.red.z, sh.green.z, sh.blue.z);
    vec3 L11      = vec3(sh.red.x, sh.green.x, sh.blue.x);

    return A[0]*Y00*L00 + A[1]*Y1_1*L1_1 + A[1]*Y10*L10 + A[1]*Y11*L11;
}

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

    float   Y00     = 0.282095;
    float   Y1_1    = 0.488603 * dir.y;
    float   Y10     = 0.488603 * dir.z;
    float   Y11     = 0.488603 * dir.x;

    result.red = vec4 (Y11, Y1_1, Y10, Y00) * color.r;
    result.green = vec4 (Y11, Y1_1, Y10, Y00) * color.g;
    result.blue = vec4 (Y11, Y1_1, Y10,Y00) * color.b;

    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.red += new.red * scale;
    accum.green += new.green * scale;
    accum.blue += new.blue * 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.red *= b;
    result.green *= b;
    result.blue *= 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));
}