球谐函数实验

实验参考这篇全面讲解的球谐光照的pdf
spherical-harmonic-lighting.pdf
以前几次看到球谐的地方就停下来了，因为最讨厌的就是这些文章直接告诉你有无数个正交的球谐基，但怎么来的不知道。也有一些文章说这些基是来自于解拉普拉斯方程，但为什么解一个偏微分方程就能得到正交的球谐基呢，也没有哪里给出明确的解释。或许这需要对这个方程解空间的结构进行学习才能了解。对于我来说超纲了，先接受这个事实，在接受这个事实的情况下，进行实验。
实验要达到的结果就是
环境辐照度图–>分解球谐系数–>通过系统重建辐照度图
说白了就是一个傅立叶变换和反变换的过程。

一. 通过ComputeShader实现SH系数计算

这里的SH系数计算我打算采用CS的方式来计算，将CubeMap传入CS，然后将预计算的采样点也传入CS。然后开8个线程，将计算结果写入一个Vector3的数组中。

public class Reconstruction : MonoBehaviour
{
    const int n_bands = 3;

    [StructLayout(LayoutKind.Explicit)]
    public unsafe struct SHSample
    {
        //public Vector3 sph;

        [FieldOffset(0)]
        public Vector3 vec;
        [FieldOffset(12)]
        public fixed float coeff[n_bands * n_bands];
    }    
    SHSample[] samples;
	....

这里先定义预计算采样点的数据结构SHSample， 因为要传给cs，所以结构体必须是指定内存布局的。SHSample的vec指明了采样的方向向量，coeff存放了在这个方向上，前9个球谐基函数的返回值。n_bands就是球谐基的层数，3层球谐基一共有9个。

接下来照搬前文提到的pdf文献中的所有函数。其中接口函数就是SH_setup_spherical_samples，它就是传入采样点个数和球谐层数，计算出来的结果放入数组，代码中的采样是cosinWeight的采样，具体可见[此处文献]

    unsafe void SH_setup_spherical_samples(SHSample[] samples, int sqrt_n_samples, int n_bands)
    {
        // fill an N*N*2 array with uniformly distributed
        // samples across the sphere using jittered stratification
        int i = 0; // array index
        float oneoverN = 1.0f / sqrt_n_samples;
        for (int a = 0; a < sqrt_n_samples; a++)
        {
            for (int b = 0; b < sqrt_n_samples; b++)
            {
                // generate unbiased distribution of spherical coords
                float x = (a + Random.Range(0f, 1.0f)) * oneoverN; // do not reuse results
                float y = (b + Random.Range(0f, 1.0f)) * oneoverN; // each sample must be random
                float theta = 2.0f * Mathf.Acos(Mathf.Sqrt(1.0f - x));
                float phi = 2.0f * Mathf.PI * y;
                //sph = new Vector3(theta, phi, 1.0f);
                // convert spherical coords to unit vector
                samples[i].vec = new Vector3(Mathf.Sin(theta) * Mathf.Cos(phi), Mathf.Sin(theta) * Mathf.Sin(phi), Mathf.Cos(theta));


                // precompute all SH coefficients for this sample
                for (int l = 0; l < n_bands; ++l)
                {
                    for (int m = -l; m <= l; ++m)
                    {
                        int index = l * (l + 1) + m;

                        fixed (float* pCoeef = samples[i].coeff)
                        {
                            pCoeef[index] = SH(l, m, theta, phi);
                        }
                    }
                }

                ++i;
            }
        }
    }


    float P(int l, int m, float x)
    {
        // evaluate an Associated Legendre Polynomial P(l,m,x) at x
        float pmm = 1.0f;
        if (m > 0)
        {
            float somx2 = Mathf.Sqrt((1.0f - x) * (1.0f + x));

            float fact = 1.0f;
            for (int i = 1; i <= m; i++)
            {
                pmm *= (-fact) * somx2;
                fact += 2.0f;
            }
        }
        if (l == m) return pmm;
        float pmmp1 = x * (2.0f * m + 1.0f) * pmm;
        if (l == m + 1) return pmmp1;
        float pll = 0.0f;
        for (int ll = m + 2; ll <= l; ++ll)
        {
            pll = ((2.0f * ll - 1.0f) * x * pmmp1 - (ll + m - 1.0f) * pmm) / (ll - m);
            pmm = pmmp1;
            pmmp1 = pll;
        }
        return pll;
    }

    int factorial(int n)
    {
        int result = 1;
        for (int i = 2; i <= n; ++i)
        {
            result *= i;
        }
        return result;

    }

    float K(int l, int m)
    {
        // renormalisation constant for SH function
        float temp = ((2.0f * l + 1.0f) * factorial(l - m)) / (4.0f * Mathf.PI * factorial(l + m));
        return Mathf.Sqrt(temp);
    }
    float SH(int l, int m, float theta, float phi)
    {
        // return a point sample of a Spherical Harmonic basis function
        // l is the band, range [0..N]
        // m in the range [-l..l]
        // theta in the range [0..Pi]
        // phi in the range [0..2*Pi]
        float sqrt2 = Mathf.Sqrt(2.0f);
        if (m == 0) return K(l, 0) * P(l, m, Mathf.Cos(theta));
        else if (m > 0) return sqrt2 * K(l, m) * Mathf.Cos(m * phi) * P(l, m, Mathf.Cos(theta));
        else return sqrt2 * K(l, -m) * Mathf.Sin(-m * phi) * P(l, -m, Mathf.Cos(theta));
    }

接下来写Start函数。在启动的时候，先预计算出32*32=1024个球面采样点，每个点的9个球谐基系数。然后它传入compute shader进行计算，computer shader开一个线程组，每个组含9个线程，每个线程计算一个投影系数，最终将反馈的结果存入coeff数组中，代码如下：

    public ComputeShader myshader;
    public Cubemap cube;
    Vector4[] coeff = new Vector4[n_bands * n_bands];
    unsafe void RunShader()
    {
        int k = myshader.FindKernel("CSMain");
        myshader.SetTexture(k, "_CubeTex", cube);

        ComputeBuffer sampleDataBuffer = new ComputeBuffer(samples.Length, Marshal.SizeOf(typeof(SHSample)));

        sampleDataBuffer.SetData(samples);
        myshader.SetBuffer(k, "sampleData", sampleDataBuffer);
        myshader.SetInt("n_samples", samples.Length);

        ComputeBuffer resultBuffer = new ComputeBuffer(n_bands * n_bands, sizeof(Vector4));
        myshader.SetBuffer(k, "Result", resultBuffer);

        myshader.Dispatch(k, 1, 1, 1);

        resultBuffer.GetData(coeff);

        for (int i = 0; i < coeff.Length; ++i)
        {
            Debug.Log(string.Format("output{0}: {1}, {2}, {3}", i, coeff[i].x, coeff[i].y, coeff[i].z));
        }
    }

    void Start()
    {
        int sqrt_n_samples = 32;
        samples = new SHSample[sqrt_n_samples * sqrt_n_samples];
        SH_setup_spherical_samples(samples, sqrt_n_samples, n_bands);

        RunShader();
    }

接下来是Compute Shader的写法， 其中开的线程是一维的，一共有9个。然后对CubeMap的采样必须用_CubeTex.SampleLevel， 用texCUBE或Sample都是编译不过的。最后采用求积分的方式将结果存在Result数组中

#pragma kernel CSMain


struct SampleStruct
{
    float3 vec;
    float coeff[9];
};

RWStructuredBuffer<SampleStruct> sampleData;
RWStructuredBuffer<float4> Result;

TextureCube<float4> _CubeTex; 
SamplerState my_Linear_clamp_sampler;
int n_samples;



[numthreads(9, 1, 1)]
void CSMain(uint3 id : SV_DispatchThreadID)
{ 
    const double weight = 4.0f * 3.14159265f;

    float4 product = 0;
    for (int i = 0; i<n_samples; ++i)
    {
        float4 color = _CubeTex.SampleLevel(my_Linear_clamp_sampler, sampleData[i].vec, 0);
        product += (color * sampleData[i].coeff[id.x]);
    }

    // divide the result by weight and number of samples
    float factor =  weight / n_samples;
    Result[id.x] = product *factor;
}

然后运行场景，在控制台得到如下日志

日志输出了，九组投影系数的RGB分量，那么这个结果对不对，需要验证一下，我们就利用这组值，来还原整个cubemap

二. 通过SH系数重建CubeMap

还原的和上一篇博客搭的框架基本一致
场景辐照度练习
主要方式就是在场景中放一个单位正方体，然后用摄像机拍6个表面，然后将系数传入Shader，在Shader中进行重建

        camera.transform.position = Vector3.zero;
        camera.transform.up = Vector3.up;
        for (int i = 0; i < 6; ++i)
        {
            rt[i] = new RenderTexture(512, 512, 0, RenderTextureFormat.ARGBFloat);
            camera.targetTexture = rt[i];
            camera.transform.forward = forwards[i];
            Shader.SetGlobalVectorArray("_coeff", coeff); //系数传入shader
            camera.RenderWithShader(SHConstructionShader, "");
            camera.targetTexture = null;
        }

而Shader的处理也很简单，在VS中计算好世界坐标，作为每个像素的采样输入，然后将这个坐标代入九个球谐基函数（参考前面给出的PDF），最后和系数相乘。最后因为要保存到文件，所以需要gamma校正

Shader "Unlit/SHConstruction"
{
    Properties
    {
    }
    
    SubShader
    {
        Tags{ "RenderType" = "Opaque" }
        LOD 100

        Pass
        {
            Cull Off
            CGPROGRAM
    #pragma vertex vert
    #pragma fragment frag

    #include "UnityCG.cginc"

            struct appdata
            {
                float4 vertex : POSITION;
                float2 uv : TEXCOORD0;
            };

            struct v2f
            {
                float2 uv : TEXCOORD0;
                float3 worldPos :TEXCOORD1;
                float4 vertex : SV_POSITION;
            };


            v2f vert(appdata v)
            {
                v2f o;
                o.vertex = UnityObjectToClipPos(v.vertex);
                o.worldPos = mul(unity_ObjectToWorld, v.vertex);

                return o;
            }


            float4 _coeff[9];
            fixed4 frag(v2f i) : SV_Target
            {

                float3 N = normalize(i.worldPos);

                float kPI = 3.1415926f;
                float kSqrtPI = sqrt(kPI);
                float kSqrt3 = sqrt(3.0f);
                float kSqrt15 = sqrt(15.0f);

                float x = N.x;
                float y = N.y;
                float z = N.z;

                float outsh[9];
                outsh[0] = 1.0f / (2.0f * kSqrtPI);
                outsh[1] = -(kSqrt3 * y) / (2.0f * kSqrtPI);
                outsh[2] = (kSqrt3 * z) / (2.0f * kSqrtPI);
                outsh[3] = -(kSqrt3 * x) / (2.0f * kSqrtPI);
                outsh[4] = (kSqrt15 * x * y) / (2.0f * kSqrtPI);
                outsh[5] = -(kSqrt15 * y * z) / (2.0f * kSqrtPI);
                outsh[6] = (sqrt(5.0f) * (3.0f * z * z - 1.0f)) / (4.0f * kSqrtPI);
                outsh[7] = -(kSqrt15 * x * z) / (2.0f * kSqrtPI);
                outsh[8] = (kSqrt15 * (x * x - y * y)) / (4.0f * kSqrtPI);

                float4 FragColor = float4(0, 0, 0, 1);
                for (int i = 0; i <9 ; ++i)
                {
                    FragColor += (outsh[i] * _coeff[i]);
                }

                FragColor = float4(LinearToGammaSpace(FragColor.rgb), 1);

                return FragColor;
            }
            ENDCG
        }
    }
}

下面是重建结果和原图的对照，看起来像那么回事