实时体积云渲染(地平线)：三.云体渲染

实时体积云渲染(地平线)：三.云体渲染

 

体渲染

最常见的体可视化方法就是Ray-marching法。该方法如下图：

在像云这样的部分透明的介质的情况下，采样点将沿着视线进一步增加到场景中，并增加固定的步长，并且在每一步中都会累积介质的密度。 在每个云密度不为零的样本处，都会向太阳发出类似的射线行进，以累积样本点与太阳之间的密度。 当前渲染像素的最终alpha值仅取决于沿视线的密度，而颜色(RGB)还取决于太阳光的颜色以及太阳与采样点之间的密度。

在本文的方法中，我们将天空视为巨大的圆顶，如下图：

我们的Ray-marching将在圆顶的低层到高层进行积分计算。云层积分步长为：

向太阳积分步长的计算公式为：

其中n和Sn分别为云层积分步数和向太阳积分步数。(ch_stop-ch_start)为云层进入点到离开点到距离。圆顶的创建可以直接通过创建限定的球体，代码如下：

void renderSphere()
{
    if (sphereVAO == 0)
    {
        glGenVertexArrays(1, &sphereVAO);
        
        unsigned int vbo, ebo;
        glGenBuffers(1, &vbo);
        glGenBuffers(1, &ebo);
        
        std::vector<glm::vec3> positions;
        std::vector<glm::vec3> uv;
        std::vector<glm::vec3> normals;
        std::vector<unsigned int> indices;
        
        const unsigned int X_SEGMENTS = 64;
        const unsigned int Y_SEGMENTS = 64;
        const float PI = 3.14159265359;
        
        float percentSph=acos(planetSize/(planetSize+startHeight))/(PI);
        float oriUV=-1.0*std::sin(percentSph*PI);
        float desUV=std::sin(percentSph*PI);
        float denoUV=desUV-oriUV;

        for (unsigned int y = 0; y <= Y_SEGMENTS; ++y)
        {
            for (unsigned int x = 0; x <= X_SEGMENTS; ++x)
            {
                float xSegment = (float)x / (float)X_SEGMENTS;
                float ySegment = (float)y / (float)Y_SEGMENTS;
                float pSph=ySegment * PI * percentSph;
                float xPos = std::cos(xSegment * 2.0f * PI) * std::sin(pSph);
                float yPos = std::cos(pSph);
                float zPos = std::sin(xSegment * 2.0f * PI) * std::sin(pSph);
                
                positions.push_back(glm::vec3(xPos, yPos, zPos));
                uv.push_back(glm::vec3((xPos-oriUV)/denoUV,0.0,(zPos-oriUV)/denoUV));
                normals.push_back(glm::vec3(xPos, yPos, zPos));
            }
        }
        
        for (int y = 0; y < Y_SEGMENTS; ++y)
        {
            for (int x = 0; x < X_SEGMENTS; ++x)
            {
                indices.push_back(y       * (X_SEGMENTS + 1) + x);
                indices.push_back((y + 1) * (X_SEGMENTS + 1) + x);
                indices.push_back(y * (X_SEGMENTS + 1) + x + 1);
                indices.push_back((y + 1) * (X_SEGMENTS + 1) + x);
                indices.push_back((y + 1) * (X_SEGMENTS + 1) + x + 1);
                indices.push_back(y * (X_SEGMENTS + 1) + x + 1);
            }
            
        }
        indexCount = indices.size();
        
        std::vector<float> data;
        for (int i = 0; i < positions.size(); ++i)
        {
            data.push_back(positions[i].x);
            data.push_back(positions[i].y);
            data.push_back(positions[i].z);
            if (uv.size() > 0)
            {
                data.push_back(uv[i].x);
                data.push_back(uv[i].y);
                data.push_back(uv[i].z);
            }
            if (normals.size() > 0)
            {
                data.push_back(normals[i].x);
                data.push_back(normals[i].y);
                data.push_back(normals[i].z);
            }
        }
        glBindVertexArray(sphereVAO);
        glBindBuffer(GL_ARRAY_BUFFER, vbo);
        glBufferData(GL_ARRAY_BUFFER, data.size() * sizeof(float), &data[0], GL_STATIC_DRAW);
        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ebo);
        glBufferData(GL_ELEMENT_ARRAY_BUFFER, indices.size() * sizeof(unsigned int), &indices[0], GL_STATIC_DRAW);
        float stride = (3 + 3 + 3) * sizeof(float);
        glEnableVertexAttribArray(0);
        glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, stride, (void*)0);
        glEnableVertexAttribArray(1);
        glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, stride, (void*)(3 * sizeof(float)));
        glEnableVertexAttribArray(2);
        glVertexAttribPointer(2, 3, GL_FLOAT, GL_FALSE, stride, (void*)(5 * sizeof(float)));
    }
    
    glBindVertexArray(sphereVAO);
    glDrawElements(GL_TRIANGLES, indexCount, GL_UNSIGNED_INT, 0);
}

这样我们Ray-marching的起点就是片源着色器中每个片源的世界坐标了。

 

云光照模型

在射线行进过程中，对于每个非零密度的样本，将使用样本处的密度以及样本与太阳之间的密度来计算光。为了近似照明，主要组成部分是比尔定律的简化版本。 比尔定律根据介质的性能来计算光的衰减。 在这种情况下，比尔定律简化为函数：

其中ds∈[0，∞)是朝向太阳的累积密度，而b∈[0，∞)确定光的吸收量。 由于ds与实际密度有些差别，因此b是一个经验性术语，用于平衡解决方案。

为了产生云体边缘高亮效果，比尔定律与Henyey-Greenstein的相位函数进行了补充，函数如下：

该函数可用于计算介质的进/出散射。 g∈[-1，1]的范围从反向散射到各向同性散射再到内部散射。 术语θ∈[-1，1]是太阳的光方向与相机发出的光线方向之间的夹角。 此函数最重要的是用于向内散射，以在面对太阳的云具有高亮效果。

除了Henyey-Greenstein的相位函数外，为了能够增加围绕太阳的集中强度，可以使用第二个术语来实现更大的日落效果：

其中csi是要在太阳周围增加的额外强度，而cse是决定额外强度应在太阳周围的集中程度的指数。 我们最后的内外散射函数如下：

其中ins∈[0,1]是内散射量，outs∈[0,1]是外散射量 。 通过从HG(θ，ins)或ISextra(θ)中选择最大强度，可以在太阳周围添加更多强度。 此外，将外散射项HG(θ，-outs)与项ivo∈[0，1]相加，以便能够线性插值，从而可以更好地控制内散射与外散射。

上述函数代码如下：

float HG(float cos_angle,float g){
    float g2=g*g;
    return (1-g2)/(pow((1+g2-2*g*cos_angle),1.5))/(PI*4);
}
float ISextra(float cos_angle){
    const float cs_i=10;//高光内散射
    const float cs_e=20;
    return cs_i*pow(SAT(cos_angle),cs_e);
}

float IOS(float cos_angle){
    const float in_s=0.2;//内散射
    const float out_s=0.1;//外散射
    const float ivo=0.5;//反散射
    return Li(max(HG(cos_angle,in_s),ISextra(cos_angle)),HG(cos_angle,-1.0*out_s),ivo);
}

 

其他非物理照明更改

然而，朝着太阳迈出非常大的步伐会造成一些问题，这些问题似乎是无法解决的，这些问题基于物理学实现。 如下图所示：

上图左：由于第一个太阳步长位置位于云层之外的可能性很高，因此几乎没有轮廓，并且看起来像一个大白点。 这在视觉上是令人不愉快的，因为云还倾向于产生更暗的边缘，如右图。 这可以通过添加完全非物理的环境吸收函数来解决：

其中d是在当前采样点采样的云密度，而osa∈[0，1]是向外散射的量。 该函数使用高度百分比向顶部添加更多的环境散射。(ps：用了这个函数，我自己跑的效果还是左图- -，只有很少云朵出现右图效果，而且效果很差)。该函数代码为：

float OutScatterAmbient(float Loddensity,float percent_height){
    const float os_a=0.9;
    //深度方向吸收
    float depth=SAT(os_a*pow(Loddensity,R(percent_height,0.3,0.9,0.5,1.0)));
//    float depth= pow(Loddensity, R(percent_height, 0.3, 0.85, 0.5, 1.3));
    //垂直方向吸收
    float verticle=pow(SAT(R(percent_height,0,0.3,0.8,1.0)),0.8);
//    float verticle= pow(R(percent_height, 0.07, 0.14, 0.1, 1.0), 0.8);
    return 1.0-SAT(depth*verticle);
}

使用比尔定律在光线行进较长步长可能引起的另一个问题是明暗区域之间的对比度高，如下图(左)：

通过使用函数钳制最大衰减可以部分解决此问题：

其中ds∈[0，∞]是光线向太阳行进采样的密度，而ac∈[0，1]是用于钳制衰减的值。 在上图(右)中显示了将衰减钳位在b = 12和ac = 0.2的情况。

但是，钳位最大衰减会导致另一个问题。 会使云层底面的灰色部分看起来没有特征。 为了解决这个问题，可以使用云密度更改衰减，如下图(右)所示：

 具体来说，这是通过如下函数完成的：

其中amin∈[0，1]是云密度d对结果的影响量。衰减项A函数代码如下：

//beer衰减项(利用clamp钳制)
float Attenuation(float Density2Sun,float density,float cos_angle){
    const float beer=12;
    float CouldAttuentionClampVal=0.2;
    float a_min=0.2;
    float prim=exp(-beer*Density2Sun);
    float scnd=exp(-beer*CouldAttuentionClampVal);
    float Aclamp=max(prim,scnd);
//    float checkval=R(cos_angle,0.0,1.0,scnd,scnd*0.5);

//    return max(checkval,prim);
    return max(density*a_min,Aclamp);
}

最终照明函数是所有三个照明函数(衰减，入/外散射和外散射环境)的组合：

代码如下：

//最后光量
vec3 Lfinal(float density,float Lodden,float density2Sun,float cos_angle,float percent_height){
    float attenuation_prob=Attenuation(density2Sun,density,cos_angle);
    float ambient_out_scatter=OutScatterAmbient(Lodden,percent_height);
    //并不需要每次计算因为matching过程中cos_angle恒定
    float sun_highlight=IOS(cos_angle);
    float attenuation = attenuation_prob*sun_highlight*ambient_out_scatter;
    attenuation=max(density * 0.4, attenuation);
    return vec3(attenuation*sun_Color);
}

ray-matching代码如下：

vec4 raymatching(){
    float fp=0.0;
    vec3 rayDir=normalize(L_view);
    vec3 samplePos=WorldPos;
    //透射率，云颜色
    float transmittance=1.0;
    vec3 cloudColor=vec3(0.0);
    float cloudDepth=0.0;
    //远云层交点
    vec3 ch_stop=cul_sphR(viewPos,rayDir,planetCenter,planetSize+startHeight+thickness);
    //步数
    int mStep=numStep;
    float ch_th=distance(WorldPos,ch_stop);//射线进出云层总长度
    float stepLen=ch_th/mStep;
    float sunStepLen=0.5*ch_th/sunStep;
    vec2 curDen=vec2(0.0);
    for (int i=0; i<mStep; i++) {
        //当前采样云高在云中高度的百分比
        float p_h=SAT((distance(samplePos,planetCenter)-planetSize-startHeight)/(thickness));
        //获取密度
        curDen=getDensity(samplePos,p_h);
        //透射率
        cloudDepth+=curDen.r;
        transmittance = SAT(exp(-1.0*cloudDepth));

        //向太阳积分
        float sunDen=0.0;
        vec3 sunSamplePos=samplePos;
        for(uint is=0;is<sunStep;is++){
            sunSamplePos+=L_sun*sunStepLen;
            float sunp_h=(distance(sunSamplePos,planetCenter)-planetSize-startHeight)/(thickness);
            sunDen+=getDensity(sunSamplePos,sunp_h).g;
        }
        cloudColor+=Lfinal(curDen.r,curDen.g,sunDen,cos_angle,p_h)*transmittance;
        //fp+=0.0000002*distance(wordSamPos,planetCenter);

        samplePos+=stepLen*rayDir*Jitter(samplePos);
    }

    return vec4(vec3(cloudColor),1.0-(transmittance));
}

因为beer定律决定了透射率，因此我们的alpha值为1-transmittance。这里我们在光线积分中加入了抖动函数Jitter，我们可以看一下有无抖动的差别(右图为均匀采样，左图为加过抖动函数)：

抖动函数代码为：

vec3 frac(vec3 p){
    return vec3(p.x-floor(p.x),p.y-floor(p.y),p.z-floor(p.z));
}

float frac(float p){
    return p-floor(p);
}

float hash(vec3 p){
    p = frac(p * 0.3183099 + 0.1);
    p *= 17.0;
    return frac(p.x * p.y * p.z * (p.x + p.y + p.z));
}

float Noise(vec3 x){
    vec3 p = vec3(floor(x.x),floor(x.y),floor(x.z));
    vec3 f = frac(x);
    f = f * f * (3.0 - 2.0 * f);

    return Li(Li(Li(hash(p + vec3(0, 0, 0)),
                        hash(p + vec3(1, 0, 0)), f.x),
                   Li(hash(p + vec3(0, 1, 0)),
                        hash(p + vec3(1, 1, 0)), f.x), f.y),
               Li(Li(hash(p + vec3(0, 0, 1)),
                        hash(p + vec3(1, 0, 1)), f.x),
                   Li(hash(p + vec3(0, 1, 1)),
                        hash(p + vec3(1, 1, 1)), f.x), f.y), f.z);
}

float map2(vec3 p){
    vec3 q = p;
    float f = 0.0f;
    f = 0.65000 * Noise(q);
    q = q * 2.02;
    f += 0.35000 * Noise(q);
    return clamp(f, 0.0, 1.0);
}

float Jitter(vec3 p){
    return R(map2(p * 300), 0, 1, 0.5, 1);
}

最后实现效果：

云体这一块只能算开个坑，效果很奇怪，没内味，最近有别的事要干，先放一放了，以后项目需要在深究吧(- -可恶)。