球谐光照总结与实现

实习需要学习了一下球谐光照,总结如下:

1. 首先是下面的rendering equation

第二部分是一个半球面上的积分。

2. 球谐变换

类比傅里叶变换[采用定义在圆上的三角函数],球谐变换采用定义在球面上的一组球谐基函数。[Spherical Harmonic Lighting: The Gritty Details]有这些函数的介绍。

现在只需要知道

  • 一个定义在球面上的函数可以近似分解成一组正交基函数的加权和:F(s) ≈ ∑ fi Yi(s).
  • 两个函数f(s),g(s)乘积的积分可以用两个基向量的点积近似表示:∫f(s)g(s)d≈∑ fi gi.
  • 坐标fi的求解,类别向量的点积,函数的点积为积分:fi=∫F(s)Yi(s)d

这样我们就把一个球面上的积分转换成了计算向量的点积,问题是如何求解坐标fi

3. 蒙特卡罗积分

解fi需要用到神奇的蒙特卡罗积分

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

下面可以写代码了,看一下最简单的“dot product lighting”,也就是漫反射模型,对应的渲染方程里面的brdf部分f(p,wo,wi)是一个常数,可以拿出积分号外,现在积分号里面剩下了∫Li(p,wi)*max(n.dot(wi),0)dwi。根据2.2,我们只需要分别投影分解Li(p,wi)和max(n.dot(wi),0)dwi两个函数,利用蒙特卡洛方法求积分的时候还可以加上可见性的判断。注意需要对模型的每一个顶点都需要计算一组基向量fi,如果加上自遮挡的判断会耗费大量时间,例如下面的dragon用了6000s来处理所有的模型顶点。

 

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

贴一下代码

1. SHLightingMesh类 负责模型顶点的处理并打包绑定到gpu.

class SHLightingMesh : public glMeshWrapper {
            GLuint* vbos_4_shcoefs_gpu;
            float9* buff_4_shcoefs_cpu;
        public:
            SHLightingMesh(const redips::Triangles* model, const SamplesOnUnitSphere& samples, redips::ShaderSource shaderSource = redips::ShaderSource()) 
                : glMeshWrapper(model, shaderSource, WrapOption::_genMipmap_){
                bindVaoAttribData(0, 1, 2);
                if (shaderSource.sourceType == redips::ShaderSource::SourceType::_null_) {
                    char strbuf[512];
                    sprintf_s(strbuf, "%s/OpenglWrappers/DemoMeshes/SHLighting", _REDIPS_ROOT_PATH_);
                    useShader(strbuf);
                }

                //allocate memory
                vbos_4_shcoefs_gpu = new GLuint[mesh->groups.size()]; {
                    glGenBuffers(mesh->groups.size(), vbos_4_shcoefs_gpu);
                    int maxFaceCnt = 0; {
                        for (auto& g : mesh->groups) maxFaceCnt = std::max(maxFaceCnt, g.faceCnt);
                        _RUNTIME_ASSERT_(maxFaceCnt > 0, "assert [maxFaceCnt>0] failed");

                        buff_4_shcoefs_cpu = new float9[maxFaceCnt * 3];
                    }
                }
                
                //build-tree, accelerte intersection calculation
                const_cast<Triangles*>(model_ptr())->buildTree();

                //calculate sh coefficients
                const std::vector<float3>& normals = mesh->normals;
                const std::vector<float3>& vertices = mesh->vertices;
                for (int i = 0; i < mesh->groups.size(); ++i) {
                    auto& g = mesh->groups[i];
                    if (g.faceCnt <= 0) continue;
                    _RUNTIME_ASSERT_(g.faceType >= GROUP_FACE_TYPE::_withnormal_, "assert [g.faceType>= GROUP_FACE_TYPE::_withnormal_] failed");

                    //fill cpu array
                    float9* ptr_f9 = buff_4_shcoefs_cpu;
                    int nid = g.nsid, fid = g.fsid, faceCnt = g.faceCnt;
                    for (int neid = nid + faceCnt; nid < neid; ++nid,++fid) {
                        const int3& indices_n = mesh->faces_vn[nid];
                        const int3& indices_v = mesh->faces_v [fid];
                        CalculateSHCoefficients(vertices[indices_v.x], normals[indices_n.x].unit(), samples, *ptr_f9++);
                        CalculateSHCoefficients(vertices[indices_v.y], normals[indices_n.y].unit(), samples, *ptr_f9++);
                        CalculateSHCoefficients(vertices[indices_v.z], normals[indices_n.z].unit(), samples, *ptr_f9++);
                    }

                    //copy to gpu and bind to vertex
                    glBindVertexArray(vaos[i]);
                    glBindBuffer(GL_ARRAY_BUFFER, vbos_4_shcoefs_gpu[i]);
                    glBufferData(GL_ARRAY_BUFFER, sizeof(float) * 9 * faceCnt * 3, buff_4_shcoefs_cpu, GL_STATIC_DRAW);
                    
                    glVertexAttribPointer(3, 3, GL_FLOAT, GL_FALSE, sizeof(float) * 9, NULL);
                    glEnableVertexAttribArray(3);
                    glVertexAttribPointer(4, 3, GL_FLOAT, GL_FALSE, sizeof(float) * 9, (const void*)(sizeof(float)*3));
                    glEnableVertexAttribArray(4);
                    glVertexAttribPointer(5, 3, GL_FLOAT, GL_FALSE, sizeof(float) * 9, (const void*)(sizeof(float)*6));
                    glEnableVertexAttribArray(5);
                }

delete buff_4_shcoefs_cpu;
}; SHLightingMesh(
const glMeshWrapper& another, redips::ShaderSource shaderSource = redips::ShaderSource()) = delete; ~SHLightingMesh() { glDeleteBuffers(mesh->groups.size(), vbos_4_shcoefs_gpu); } void draw() { if (!m_shader) { std::cerr << "shader error" << std::endl; return; }; m_shader->Use(); for (int i = 0; i < meshCnt; i++) { unsigned int flags = 0; if (meshMtls[i]->texture_ka != NULL && meshFaceTypes[i] == redips::GROUP_FACE_TYPE::_withtex_) { glActiveTexture(GL_TEXTURE0); glBindTexture(GL_TEXTURE_2D, mtlTextureHandle[meshMtls[i]->texture_ka]); glUniform1i(glGetUniformLocation(m_shader->Program, shaderAmbientTextureUniformStr), shaderAmbientTextureLocation); flags |= 1u; } { redips::float3 color = meshMtls[i]->ambient; glUniform3f(glGetUniformLocation(m_shader->Program, shaderAmbientColorUniformStr), color.x, color.y, color.z); } if (meshMtls[i]->texture_kd != NULL && meshFaceTypes[i] == redips::GROUP_FACE_TYPE::_withtex_) { glActiveTexture(GL_TEXTURE1); glBindTexture(GL_TEXTURE_2D, mtlTextureHandle[meshMtls[i]->texture_kd]); glUniform1i(glGetUniformLocation(m_shader->Program, shaderDiffuseTextureUniformStr), shaderDiffuseTextureLocation); flags |= 2u; } { redips::float3 color = meshMtls[i]->diffuse; glUniform3f(glGetUniformLocation(m_shader->Program, shaderDiffuseColorUniformStr), color.x, color.y, color.z); } glUniform1ui(glGetUniformLocation(m_shader->Program, shaderSurfaceTypeUniformStr), flags); glBindVertexArray(vaos[i]); glDrawArrays(GL_TRIANGLES, 0, 3 * meshFaceCnt[i]); } } private: void CalculateSHCoefficients(const float3& vertice, const float3& normal, const SamplesOnUnitSphere& samples, float9& result) { result = 0; redips::HitRecord hitRecord; for (int i = 0; i < samples.size(); ++i) { auto fs_i = std::max(samples[i].dot(normal.unit())*-1, 0.0f); if (fs_i > 0.0) { Ray ray(vertice, normal); if (!const_cast<Triangles*>(model_ptr())->intersect(ray, hitRecord)) { result = result + samples.samples_vec9[i] * fs_i; } } } result = result * (4 * PI / samples.size()); } const GLchar* shaderSurfaceTypeUniformStr = "material.flags"; const GLchar* shaderAmbientColorUniformStr = "material.ambient"; const GLchar* shaderDiffuseColorUniformStr = "material.diffuse"; const GLchar* shaderAmbientTextureUniformStr = "material.ambientTexture"; const GLchar* shaderDiffuseTextureUniformStr = "material.diffuseTexture"; const GLuint shaderAmbientTextureLocation = 0; const GLuint shaderDiffuseTextureLocation = 1; };

2. SamplesOnUnitSphere 类 负责在单位球面上随机生成均匀的采样点

class SamplesOnUnitSphere {
        public:
            std::vector<float9> samples_vec9;
            std::vector<redips::float3> samples_vec3;
            int size() const{ return samples_vec3.size(); }
            SamplesOnUnitSphere(int sampleCnt) {
                std::default_random_engine generator;
                std::uniform_real_distribution<float> normal_floats(0, 1);
                samples_vec3.resize(sampleCnt);
                samples_vec9.resize(sampleCnt);
                for (int i = 0; i < sampleCnt; ++i) {
                    auto x = normal_floats(generator);
                    auto y = normal_floats(generator);
                    auto theta = acos(sqrt(x)) * 2 - PI * 0.5;
                    auto phi = (2 * y - 1)*PI;
                    samples_vec3[i] = redips::float3(cos(theta)*sin(phi), sin(theta), cos(theta)*cos(phi));
                    samples_vec9[i] = float9(samples_vec3[i]);
                }
            }

            redips::float3& operator[](size_t idx) {
                _RUNTIME_ASSERT_(idx < samples_vec3.size(), "assert [idx<samples_vec3.size()] failed");
                return samples_vec3[idx];
            }
            const redips::float3& operator[](size_t idx) const {
                _RUNTIME_ASSERT_(idx < samples_vec3.size(), "assert [idx<samples_vec3.size()] failed");
                return samples_vec3[idx];
            }
        };

3. Cubemap类 里面的computeSH函数负责将来自六个面的环境光编码成球谐基向量

        class Cubemap {
            redips::FImage* faces[6];
            //a cube with side-length 2 centered at (0,0,0)
            redips::Triangles cube = redips::float3(2, 2, 2);
        public:
            Cubemap(const char* picnames[6]) {
                for (int i = 0; i < 6; ++i) faces[i] = new redips::FImage(picnames[i]);
            }
            ~Cubemap() { for (int i = 0; i < 6; ++i) delete faces[i]; }
            redips::float3 texture(const redips::float3& direction) {
                using namespace redips;

                //find the intersection
                redips::Ray ray({ 0,0,0 }, direction);
                redips::HitRecord recorder;
                cube.intersect(ray, recorder);
                auto hitp = ray.ori + ray.dir * recorder.distance;

                //fabs(one-dim-of-$hitp$) must equal 1, because $hitp$ is on one face of $cube$
                int FaceId = -1;
                float2 texcoord;
                const float epsi = 1e-5;
                if (fabs(hitp.x - 1) < epsi) {
                    FaceId = 0;
                    texcoord = float2(hitp.z, hitp.y)*0.5 + float2{0.5f, 0.5f};
                }
                else if (fabs(hitp.x + 1) < epsi) {
                    FaceId = 1;
                    texcoord = float2(hitp.z*-1, hitp.y)*0.5 + float2{ 0.5f, 0.5f };
                }
                else if (fabs(hitp.y - 1) < epsi) {
                    FaceId = 2;
                    texcoord = float2(hitp.x, hitp.z)*0.5 + float2{ 0.5f, 0.5f };
                }
                else if (fabs(hitp.y + 1) < epsi) {
                    FaceId = 3;
                    texcoord = float2(hitp.x, hitp.z*-1)*0.5 + float2{ 0.5f, 0.5f };
                }
                else if (fabs(hitp.z - 1) < epsi) {
                    FaceId = 4;
                    texcoord = float2(hitp.x*-1, hitp.y)*0.5 + float2{ 0.5f, 0.5f };
                }
                else if (fabs(hitp.z + 1) < epsi) {
                    FaceId = 5;
                    texcoord = float2(hitp.x, hitp.y)*0.5 + float2{ 0.5f, 0.5f };
                }

                if (FaceId < 0) {
                    std::cerr << "[Cubemap] : Ops, something wrong in Cubemap.texture" << std::endl;
                    return redips::float3();
                }
                return faces[FaceId]->tex2d(texcoord.x, texcoord.y);
            }

            void computeSH(SamplesOnUnitSphere& samples, float9& channel_r, float9& channel_g, float9& channel_b) {
                channel_r = 0, channel_g = 0, channel_b = 0;
                for (int i = 0; i < samples.size(); ++i) {
                    auto Es_i = texture(samples[i]*-1);
                    auto& Ys = samples.samples_vec9[i];
                    channel_r = channel_r + Ys*Es_i.x;
                    channel_g = channel_g + Ys*Es_i.y;
                    channel_b = channel_b + Ys*Es_i.z;
                }
                auto scale = (4 * PI / samples.size()); {
                    channel_r = channel_r * scale;
                    channel_g = channel_g * scale;
                    channel_b = channel_b * scale;
                }
            }
        };
4. main.cpp负责组装
#include <common/glfwApp.h>
#include <OpenglWrappers/glHelper.h>
#include <OpenglWrappers/DemoMeshes/Skybox.h>
#include <OpenglWrappers/DemoMeshes/SHLightingMesh.h>

//#define MODEL_PATH "E:/Documents/models/spheres/sphere0.obj"
#define MODEL_PATH "E:/Documents/models/dragon.obj"

//window size
const redips::int2 Winsize{ 960,960 };

//glfw setup
auto application = redips::glfw::getInstance(Winsize.x, Winsize.y);

//screen capture
auto screenCapture = redips::glScreenCapture::getInstance(Winsize);

//skybox
const char* skybox_faces[] = {
    "D:/Documents/resources/skyboxes/pure/black.jpg",
    "D:/Documents/resources/skyboxes/pure/blue.jpg",
    "D:/Documents/resources/skyboxes/pure/black.jpg",
    "D:/Documents/resources/skyboxes/pure/red.jpg",
    "D:/Documents/resources/skyboxes/pure/black.jpg",
    "D:/Documents/resources/skyboxes/pure/green.jpg"
};
auto skybox = redips::SkyBox(skybox_faces);

//standard pinhole camera
redips::PhC phc(60, 1.0f, 0.1f, 10000);

//load a obj and then wrap into a glMesh。
redips::SphericalHarmonic::SHLightingMesh * shmesh = nullptr;

//mesh center
redips::float3 heart = redips::float3{ 0,0,0 };

//keyboard event
void movement() {
    if (application->keydown(GLFW_KEY_M)) {
        screenCapture->capture("capture.bmp", 0);
    }
}

//need register a display-callback function to tell glfwApp what to render
void display() {
    using namespace redips;
    
    shmesh->uniformMat44f("model", shmesh->model_ptr()->Transform().transpose().ptr());
    shmesh->uniformMat44f("projection_view", phc.glProjectionView().ptr());
    shmesh->draw();

    skybox.render(&phc);

    movement();
}

void initialize() {
    using namespace redips;

    //setup camera
    phc.lookAt(heart + float3(0, 0, 20), heart, float3(0, 1, 0));
    phc.setResolution(Winsize.x, Winsize.y);

    //generate samples
    SphericalHarmonic::SamplesOnUnitSphere samples(500);

    //calculate light
    SphericalHarmonic::Cubemap cubemap(skybox_faces);
    SphericalHarmonic::float9 light_r, light_g, light_b;
    cubemap.computeSH(samples, light_r, light_g, light_b);

    //warp mesh
    shmesh = new redips::SphericalHarmonic::SHLightingMesh(new Triangles(MODEL_PATH), samples);
    shmesh ->littedBy(light_r, light_g, light_b);

//transform mesh const_cast<Triangles*>(shmesh->model_ptr())->setTransform( Mat44f::scale(float3(60,60,60))* //Mat44f::rotatey(RAD(90))* Mat44f::translation(shmesh->model_ptr()->aabb_R().heart()*-1) ); } int main() { initialize(); application->registerDisplayCallback(display); application->bindCamera(&phc); application->loop(); return 0; }



redips
是之前封装的基础操作的类库
原文地址:https://www.cnblogs.com/redips-l/p/9221982.html