十六，镜面IBL--预滤波环境贴图

又到了开心的公式时刻了。
先看看渲染方程

现在关注第二部分，镜面反射。

其中
这里很棘手，与输入wi和输出w0都有关系，所以，再近似

其中第一部分，就是预滤波环境贴图，形式上与前面的辐照度图很相似，那么能不能用同样的方法呢？
先看看镜面反射和漫反射的图

可以看到，镜面反射是绕着出射向量的一个范围（成为波瓣），而漫反射是绕着法线方向均匀分布的。
再想想积分辐照度图时，是以法线向量为中心，进行积分的。
那很自然的想到，积分镜面反射的预滤波环境贴图可以以出射向量为中心，在波瓣范围内积分。

然而， 波瓣有大有小，是因为粗糙度不同，

所以，不能只积分一次，而是多次，按照不同粗糙度积分后写到mipmap，或者单独的纹理中。这里为了方便，分别写到不同的纹理中。

那么该如何积分呢？辐照度图是在经度0到360，纬度0到90内均匀积分。

而镜面反射中，给定入射方向，波瓣指向方向就是微平面半向量的反射方向。所以，只在波瓣内积分就可以了，即重要性采样。
这时就可以使用蒙特卡洛积分，即在大数定律基础上，采取N样本即可。N越大越准。pdf为概率密度函数。

比如

采样样本越多，越靠近中间范围。因为中间范围概率大。

以上为均匀采样，

如果采样样本有偏，则会更快收敛。比如通过低差异序列获取样本。
float RadicalInverse_VdC(uint bits)
{
bits = (bits << 16u) | (bits >> 16u);
bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
return float(bits) * 2.3283064365386963e-10; // / 0x100000000
}

vec2 Hammersley(uint i, uint N)
{
return vec2(float(i)/float(N), RadicalInverse_VdC(i));
}
或者无位运算的
"float VanDerCorpus(uint n, uint base) "
"{ "
" float invBase = 1.0 / float(base); "
" float denom = 1.0; "
" float result = 0.0; "
" for (uint i = 0u; i < 32u; ++i) "
" { "
" if (n > 0u) "
" { "
" denom = mod(float(n), 2.0); "
" result += denom * invBase; "
" invBase = invBase / 2.0; "
" n = uint(float(n) / 2.0); "
" } "
" } "
"return result; "
"} "
" "
"vec2 HammersleyNoBitOps(uint i, uint N) "
"{ "
" return vec2(float(i) / float(N), VanDerCorpus(i, 2u)); "
"} "

然后根据法线方向，粗糙度和低差异序列生成采样向量，该向量大体围绕着预估的波瓣方向。
“vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness)”
“{”
“float a = roughness * roughness;”
“float phi = 2.0 * PI * Xi.x;”
“float cosTheta = sqrt((1.0 - Xi.y)/(1.0+(a*a-1.0) * Xi.y));”
“float sinTheta = sqrt(1.0 - cosTheta * cosTheta);”
“vec3 H;”
“H.x = cos(phi) * sinTheta;”
“H.y = sin(phi) * sinTheta;”
“H.z = cosTheta;”
“vec3 up = abs(N.z) < 0.999 ? vec3(0.0,0.0,1.0) : vec3(1.0,0.0,0.0);”
“vec3 tangent = normalize(cross(up,N));”
“vec3 bitangent = cross(N,tangent);”
“vec3 sampleVec = tangent * H.x + bitangent * H.y + N * H.z;”
“return normalize(sampleVec);”
“}”
运行结果如下

代码如下
#include <osg/TextureCubeMap>
#include <osg/TexGen>
#include <osg/TexEnvCombine>
#include <osgUtil/ReflectionMapGenerator>
#include <osgDB/ReadFile>
#include <osgViewer/Viewer>
#include <osg/NodeVisitor>
#include <osg/ShapeDrawable>

static const char * vertexShader =
{
//“#version 120 core\n”
“in vec3 aPos;\n”
“varying vec3 localPos;\n”
“void main(void)\n”
“{\n”
“localPos = aPos;\n”
" gl_Position = ftransform();\n"
//“gl_Position = view * view * vec4(aPos,1.0);”
“}\n”
};

static const char psShader =
{
“varying vec3 localPos;\n”
“uniform samplerCube environmentMap;”
“uniform float roughness;”
“const float PI = 3.1415926;”
“float VanDerCorpus(uint n, uint base) "
“{ "
" float invBase = 1.0 / float(base); "
" float denom = 1.0; "
" float result = 0.0; "
" for (uint i = 0u; i < 32u; ++i) "
" { "
" if (n > 0u) "
" { "
" denom = mod(float(n), 2.0); "
" result += denom * invBase; "
" invBase = invBase / 2.0; "
" n = uint(float(n) / 2.0); "
" } "
" } "
“return result; "
“} "
" "
“vec2 HammersleyNoBitOps(uint i, uint N) "
“{ "
" return vec2(float(i) / float(N), VanDerCorpus(i, 2u)); "
“} "
//“float RadicalInverse_Vdc(uint bits)\n”
//”{”
//“bits = (bits << 16u) | (bits >> 16u);”
//“bits = ((bits & 0x55555555u) << 1u ) | (bits & 0xAAAAAAAAu) >> 1u);”
//“bits = ((bits & 0x33333333u) << 2u ) | (bits & 0xCCCCCCCCu) >> 2u);”
//“bits = ((bits & 0x0F0F0F0Fu) << 4u ) | (bits & 0xF0F0F0F0u) >> 4u);”
//“bits = ((bits & 0x00FF00FFu) << 8u ) | (bits & 0xFF00FF00u) >> 8u);”
//“return float(bits) * 2.3283064365386963e-10;”
//”}”
//“vec2 Hammersley(uint i, uint N)”
//”{”
//“return vec2(float(i) / float(N), RadicalInverse_Vdc(i));”
//”}"
“vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness)”
“{”
“float a = roughness * roughness;”
“float phi = 2.0 * PI * Xi.x;”
"float cosTheta = sqrt((1.0 - Xi.y)/(1.0+(aa-1.0) * Xi.y));"
“float sinTheta = sqrt(1.0 - cosTheta * cosTheta);”
“vec3 H;”
“H.x = cos(phi) * sinTheta;”
“H.y = sin(phi) * sinTheta;”
“H.z = cosTheta;”
“vec3 up = abs(N.z) < 0.999 ? vec3(0.0,0.0,1.0) : vec3(1.0,0.0,0.0);”
“vec3 tangent = normalize(cross(up,N));”
“vec3 bitangent = cross(N,tangent);”
“vec3 sampleVec = tangent * H.x + bitangent * H.y + N * H.z;”
“return normalize(sampleVec);”
“}”
"void main() "
"{ "
" vec3 N = normalize(localPos); "
" vec3 R = N; "
" vec3 V = R; "
" "
" const uint SAMPLE_COUNT = 1024u; "
" float totalWeight = 0.0; "
" vec3 prefilteredColor = vec3(0.0); "
" for (uint i = 0u; i < SAMPLE_COUNT; ++i) "
" { "
" vec2 Xi = HammersleyNoBitOps(i, SAMPLE_COUNT); "
" vec3 H = ImportanceSampleGGX(Xi, N, roughness); "
" vec3 L = normalize(2.0 * dot(V, H) * H - V); "
" "
" float NdotL = max(dot(N, L), 0.0); "
" if (NdotL > 0.0) "
" { "
" prefilteredColor += texture(environmentMap, L).rgb * NdotL; "
" totalWeight += NdotL; "
" } "
" } "
" prefilteredColor = prefilteredColor / totalWeight; "
" "
" gl_FragColor = vec4(prefilteredColor, 1.0); "
"} "
};
class MyNodeVisitor : public osg::NodeVisitor
{
public:
MyNodeVisitor() : osg::NodeVisitor(osg::NodeVisitor::TRAVERSE_ALL_CHILDREN)
{

}
void apply(osg::Geode& geode)
{int count = geode.getNumDrawables();for (int i = 0; i < count; i++){osg::ref_ptr<osg::Geometry> geometry = geode.getDrawable(i)->asGeometry();if (!geometry.valid()){continue;}osg::Array* vertexArray = geometry->getVertexArray();geometry->setVertexAttribArray(1, vertexArray);}traverse(geode);
}

};

int main()
{

osg::ref_ptr<osgViewer::Viewer> viewer = new osgViewer::Viewer;osg::ref_ptr<osg::TextureCubeMap> tcm = new osg::TextureCubeMap;
tcm->setTextureSize(128, 128);
tcm->setFilter(osg::Texture::MIN_FILTER, osg::Texture::LINEAR_MIPMAP_LINEAR);
tcm->setFilter(osg::Texture::MAG_FILTER, osg::Texture::LINEAR);
tcm->setWrap(osg::Texture::WRAP_S, osg::Texture::CLAMP_TO_EDGE);
tcm->setWrap(osg::Texture::WRAP_T, osg::Texture::CLAMP_TO_EDGE);
tcm->setWrap(osg::Texture::WRAP_R, osg::Texture::CLAMP_TO_EDGE);std::string strImagePosX = "D:/hdr/Right face camera.bmp";
osg::ref_ptr<osg::Image> imagePosX = osgDB::readImageFile(strImagePosX);
tcm->setImage(osg::TextureCubeMap::POSITIVE_X, imagePosX);
std::string strImageNegX = "D:/hdr/Left face camera.bmp";
osg::ref_ptr<osg::Image> imageNegX = osgDB::readImageFile(strImageNegX);
tcm->setImage(osg::TextureCubeMap::NEGATIVE_X, imageNegX);std::string strImagePosY = "D:/hdr/Front face camera.bmp";;
osg::ref_ptr<osg::Image> imagePosY = osgDB::readImageFile(strImagePosY);
tcm->setImage(osg::TextureCubeMap::POSITIVE_Y, imagePosY);
std::string strImageNegY = "D:/hdr/Back face camera.bmp";;
osg::ref_ptr<osg::Image> imageNegY = osgDB::readImageFile(strImageNegY);
tcm->setImage(osg::TextureCubeMap::NEGATIVE_Y, imageNegY);std::string strImagePosZ = "D:/hdr/Top face camera.bmp";
osg::ref_ptr<osg::Image> imagePosZ = osgDB::readImageFile(strImagePosZ);
tcm->setImage(osg::TextureCubeMap::POSITIVE_Z, imagePosZ);std::string strImageNegZ = "D:/hdr/Bottom face camera.bmp";
osg::ref_ptr<osg::Image> imageNegZ = osgDB::readImageFile(strImageNegZ);
tcm->setImage(osg::TextureCubeMap::NEGATIVE_Z, imageNegZ);
tcm->setUseHardwareMipMapGeneration(true);
float minMipMapLevel = 0.0;
float maxMipMapLevel = 4.0;
tcm->setMinLOD(minMipMapLevel);
tcm->setMaxLOD(maxMipMapLevel);osg::ref_ptr<osg::Box> box = new osg::Box(osg::Vec3(0, 0, 0), 10);
osg::ref_ptr<osg::ShapeDrawable> drawable = new osg::ShapeDrawable(box);
osg::ref_ptr<osg::Geode> geode = new osg::Geode;
geode->addDrawable(drawable);
MyNodeVisitor nv;
geode->accept(nv);
osg::ref_ptr<osg::StateSet> stateset = geode->getOrCreateStateSet();
stateset->setTextureAttributeAndModes(0, tcm, osg::StateAttribute::OVERRIDE | osg::StateAttribute::ON);//shaderosg::ref_ptr<osg::Shader> vs1 = new osg::Shader(osg::Shader::VERTEX, vertexShader);
osg::ref_ptr<osg::Shader> ps1 = new osg::Shader(osg::Shader::FRAGMENT, psShader);
osg::ref_ptr<osg::Program> program1 = new osg::Program;
program1->addShader(vs1);
program1->addShader(ps1);
program1->addBindAttribLocation("aPos", 1);osg::ref_ptr<osg::Uniform> environmentMapUniform = new osg::Uniform("environmentMap", 0);
stateset->addUniform(environmentMapUniform);
float theMip = 3.0;
float roughness = theMip / maxMipMapLevel;
osg::ref_ptr<osg::Uniform> roughnessUniform = new osg::Uniform("roughness", roughness);
stateset->addUniform(roughnessUniform);stateset->setAttribute(program1, osg::StateAttribute::ON);viewer->setSceneData(geode);
viewer->realize();
return viewer->run();

}