/********************************************************
    © 2020 Continuum Graphics LLC. All Rights Reserved
 ********************************************************/

#if !defined _SPHERICALHARMONICS_
    #define _SPHERICALHARMONICS_

    // Generates a uniform distribution of directions over a unit sphere. 
    // Adapted from http://www.pbr-book.org/3ed-2018/Monte_Carlo_Integration/2D_Sampling_with_Multidimensional_Transformations.html#fragment-SamplingFunctionDefinitions-6
    // azimuthX and zenithY are both in [0, 1]. You can use random value, stratified, etc.
    // Top and bottom sphere pole (+-zenith) are along the Y axis.
    vec3 shGetUniformSphereSample(float azimuthX, float zenithY) {
        float phi = 2.0f * PI * azimuthX;
        float z = 1.0f - 2.0f * zenithY;
        float r = sqrt(max(0.0f, 1.0f - z * z));
        return vec3(r * cos(phi), z, r * sin(phi));
    }

    vec4 shZero() {
        return vec4(0.0);
    }

    // Evaluates spherical harmonics basis for a direction dir. (from [2] Appendix A2)
    // (evaluating the associated Legendre polynomials using the polynomial forms)
    vec4 shEvaluate(vec3 dir) {
        vec4 result;
        result.x = 0.28209479177387814347403972578039f;			// L=0 , M= 0
        result.y =-0.48860251190291992158638462283836f * dir.y;	// L=1 , M=-1
        result.z = 0.48860251190291992158638462283836f * dir.z;	// L=1 , M= 0
        result.w =-0.48860251190291992158638462283836f * dir.x;	// L=1 , M= 1
        return result;
    }

    // Recovers the value of a SH function in the direction dir.
    float shUnproject(vec4 functionSh, vec3 dir) {
        vec4 sh = shEvaluate(dir);
        return dot(functionSh, sh);
    }
    vec3 shUnproject(vec4 functionShX, vec4 functionShY, vec4 functionShZ, vec3 dir) {
        vec4 sh = shEvaluate(dir);
        return vec3(dot(functionShX, sh), dot(functionShY, sh), dot(functionShZ, sh));
    }

    // Projects a cosine lobe function, with peak value in direction dir, into SH. (from [4])
    // The integral over the unit sphere of the SH representation is PI.
    vec4 shEvaluateCosineLobe(vec3 dir) {
        vec4 result;
        result.x = 0.8862269254527580137f;			// L=0 , M= 0
        result.y =-1.0233267079464884885f * dir.y;	// L=1 , M=-1
        result.z = 1.0233267079464884885f * dir.z;	// L=1 , M= 0
        result.w =-1.0233267079464884885f * dir.x;	// L=1 , M= 1
        return result;
    }

    // Projects a Henyey-Greenstein phase function, with peak value in direction dir, into SH. (from [11])
    // The integral over the unit sphere of the SH representation is 1.
    vec4 evaluatePhaseHG(vec3 dir, float g) {
        vec4 result;
        const float factor = 0.48860251190291992158638462283836 * g;
        result.x = 0.28209479177387814347403972578039;	// L=0 , M= 0
        result.y =-factor * dir.y;						// L=1 , M=-1
        result.z = factor * dir.z;						// L=1 , M= 0
        result.w =-factor * dir.x;						// L=1 , M= 1
        return result;
    }

    float shUnprojectCosineLobe(vec4 functionSh, vec3 dir) {
        vec4 sh = shEvaluateCosineLobe(dir);
        return dot(functionSh, sh);
    }
    vec3 shUnprojectCosineLobe(vec4 functionShX, vec4 functionShY, vec4 functionShZ, vec3 dir) {
        vec4 sh = shEvaluateCosineLobe(dir);
        return vec3(dot(functionShX, sh), dot(functionShY, sh), dot(functionShZ, sh));
    }

    // Adds two SH functions together.
    vec4 shAdd(vec4 shL, vec4 shR) {
        return shL + shR;
    }

    // Scales a SH function uniformly by v.
    vec4 shScale(vec4 sh, float v) {
        return sh * v;
    }

    // Operates a rotation of a SH function.
    vec4 shRotate(vec4 sh, mat3 rotation) {
        // TODO verify and optimize
        vec4 result;
        result.x = sh.x;
        vec3 tmp = vec3(sh.w, sh.y, sh.z);		// undo direction component shuffle to match source/function space
        result.yzw = (tmp * rotation).yzx;		// apply rotation and re-shuffle
        return result;
    }

    // Integrates the product of two SH functions over the unit sphere.
    float shFuncProductIntegral(vec4 shL, vec4 shR) {
        return dot(shL, shR);
    }

    // Computes the SH coefficients of a SH function representing the result of the multiplication of two SH functions. (from [4])
    // If sources have N bands, this product will result in 2N*1 bands as signal multiplication can add frequencies (think about two lobes intersecting).
    // To avoid that, the result can be truncated to N bands. It will just have a lower frequency, i.e. less details. (from [2], SH Products p.7)
    vec4 shProduct(vec4 shL, vec4 shR) {
        const float factor = 1.0f / (2.0f * sqrt(PI));
        return factor * vec4(
            shL.x*shR.w + shL.w+shR.x,
            shL.y*shR.w + shL.w*shR.y,
            shL.z*shR.w + shL.w*shR.z,
            dot(shL,shR)
        );
    }

    // Convolves a SH function using a Hanning filtering. This helps reducing ringing and negative values. (from [2], Windowing p.16)
    // A lower value of w will reduce ringing (like the frequency of a filter)
    vec4 shHanningConvolution(vec4 sh, float w) {
        vec4 result = sh;
        float invW = 1.0 / w;
        float factorBand1 =(1.0 + cos( PI * invW )) / 2.0f;
        result.y *= factorBand1;
        result.z *= factorBand1;
        result.w *= factorBand1;
        return result;
    }

    // Convolves a SH function using a cosine lob. This is tipically used to transform radiance to irradiance. (from [3], eq.7 & eq.8)
    vec4 shDiffuseConvolution(vec4 sh) {
        vec4 result = sh;
        // L0
        result.x   *= PI;
        // L1
        result.yzw *= 2.0943951023931954923f;
        return result;
    }

    /*
    vec3 FromSH(vec4 cR, vec4 cG, vec4 cB, vec3 lightDir) {
        const float sqrt1OverPI = sqrt(rPI);
        const float sqrt3OverPI = sqrt(3.0 * rPI);
        const vec2 halfnhalf = vec2(0.5, -0.5);
        const vec2 sqrtOverPI = vec2(sqrt1OverPI, sqrt3OverPI);
        const vec4 foo = halfnhalf.xyxy * sqrtOverPI.xyyy;

        vec4 sh = foo * vec4(1.0, lightDir.yzx);

        return vec3(
            dot(sh,cR),
            dot(sh,cG),
            dot(sh,cB)
        );
    }

    vec4 ToSH(float value, vec3 dir) {
        const float transferl1 = 0.3849 * PI;
        const float sqrt1OverPI = sqrt(rPI);
        const float sqrt3OverPI = sqrt(rPI * 3.0);

        const vec2 halfnhalf = vec2(0.5, -0.5);
        const vec2 transfer = vec2(PI * sqrt1OverPI, transferl1 * sqrt3OverPI);

        const vec4 foo = halfnhalf.xyxy * transfer.xyyy;

        return foo * vec4(1.0, dir.yzx) * value;
    }
    */
#endif
