【体绘制】raymarching算法

最近研究云的渲染接触到了体绘制，在shadertoy上看了一些例子，不过它们都是全场景体绘制的，场景中的所有物体，地形均由数学方法生成。然后用一个quad覆盖整个屏幕，对每个片段进行raymarching，甚至是raytracing。好像有一个64kb的比赛，就是用最小的程序大小，看谁能够渲染出最漂亮的图像，令人不得不感叹数学的神奇。

要想系统的了解体绘制，推荐有一本非常好的书，Real-Time Volume Graphics，里面不仅介绍了用于科学领域的体绘制，而且针对游戏开者，非常详细的讲述了如何在GPU上实现并优化，如何将体绘制与传统的基于几何面片的渲染方式集成起来。

我自己也粗陋的实现了一下raymarching算法，绘制了一个3D柏林噪声，可以下载下来试试，效果我还挺满意的：

3D纹理

体绘制用到的纹理是3D的，第一次接触3D纹理的我还有点小懵逼…其实也很好理解，就是给2D纹理再加上一个维度。存储方式通常就是n张w*h的图片啦，不过在OpenGL中载入的时候会稍麻烦一些：

GLuint load3DPerlinNoise()
{
    GLuint perlinNoise;
    //perlinNoise = -1;
    glGenTextures(1, &perlinNoise);
    glBindTexture(GL_TEXTURE_3D, perlinNoise);

    char fileName[100];
    unsigned char* noiseData = new unsigned char[256*256*256*3];

    for(int i = 0; i < 256; i++){
        sprintf(fileName, "Perlin_Noise_BMP/Noise_Perlin%04d.bmp",i);
        int width, height, channel;
        unsigned char *data = SOIL_load_image(fileName, &width, &height, &channel, SOIL_LOAD_RGB);
        memcpy(noiseData + i*65536*3, data, sizeof(unsigned char)*width*height*3);
        SOIL_free_image_data(data);
    }

    glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB, 256, 256, 256, 0, GL_RGB, GL_UNSIGNED_BYTE, (GLvoid*)noiseData);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_BORDER);
    delete noiseData;
    return perlinNoise;
}

一个是纹理数据存入buffer的时候，要加上一个维度的offset。再一个就是生成和绑定纹理的时候，都改成了GL_TEXTURE_3D，其他都是大同小异了。最后传入顶点纹理坐标时也变成了3维，从(0,0,0)到(1,1,1)。

raymarching 算法

raymarching 算法思想很直观：首先有一个3D的体纹理，然后从相机发射n条射线，射线有一个采样的步长。当射线处在体纹理中时，每个步长采一次样，获取纹理值（实际上表示该点的密度值），计算光照，然后和该条射线当前累积的颜色值进行混合。

为什么这样就可以渲染出正确的图案呢？因为光路是可逆的，从光源射出的光线经过散射，最终进入摄像机的效果等同于从摄像机发出的射线进行着色和采样，这个raytracing的道理是一样的。

这种算法很适合在GPU上实现，因为每条射线的计算都是独立并行的，GPU在大量并行计算上有先天的优势。为了在GPU上实现，我们需要解决的问题主要有2个：

1. 哪些片段需要raymarching。
2. raymarching的方向和终点在哪里。

下面来逐一进行解决。

确定raymarching的片段

体绘制首先需要一个载体(proxy geometry)，也就是为了确定屏幕上的哪些像素是属于某个体纹理的。这个问题很容易就让人联想到包围盒，问题也就引刃而解。

我们只需将体纹理的包围盒绘制出来，那么包围盒在屏幕上覆盖的片段自然就是需要进行raymarching的了。如下图所示：

随后只需要执行raymarching的片断着色器即可。

raymarching的方向和终止点

在使用包围盒作为体绘制的载体时，起/终点就是每根ray进出包围盒时的两个交点。关于如何得到这两个点的坐标，有一种2个pass的算法：

• 绘制包围盒的背面，即将OpenGL背面剔除设置为GL_FRONT，并将每个片段的世界坐标保存在纹理缓存中。
• 绘制正面，将每个片段的坐标和上一个pass中的每个片段的坐标相减，即可的到ray的方向和长度，然后进行raymarching算法，达到长度终止即可（采样时要转换为纹理坐标）。

2个pass过程相对繁琐，我使用了1个pass：

• 获得片段的世界坐标，然后减去视点位置得到ray的方向。然后每次步进时都判断当前的纹理坐标是否超出了包围盒的边界，一旦超出，就停止算法。

确定ray的起始和终止点，是一个对资源消耗很大的过程，这方面也有许多内容值得研究，作为一个探索原理的程序，就以体现算法思想为主了。

实现

顶点着色器如下：

#version 450 core

layout(location = 0) in vec3 vertices;
layout(location = 1) in vec3 UVW;

uniform mat4 M;
uniform mat4 V;
uniform mat4 P;
uniform vec3 viewPos;

out vec3 textCoord;
out vec3 fragWorldPos;

void main()
{
    gl_Position = P*V*M*vec4(vertices, 1.0);
    // 把片段的世界坐标以及纹理坐标传递给下一个阶段
    fragWorldPos = vec3(M*vec4(vertices, 1.0));
    textCoord = UVW;
}

然后是片断着色器：

#version 450 core

uniform sampler3D noiseSampler;
out vec4 outColor;

in vec3 textCoord;
in vec3 fragWorldPos;
uniform vec3 viewPos;

uniform mat4 M;

struct Ray{
    vec3 o;
    vec3 d;
}eyeRay;

// 这个着色函数copy的shadertoy上iq大神的，关于体绘制的着色还没有开始深入的学习，日后再说哈哈
vec4 integrate( in vec4 sum, in float dif, in float den, in vec3 bgcol, in float t )
{
    // lighting
    vec3 lin = vec3(0.65,0.7,0.75)*1.4 + vec3(1.0, 0.6, 0.3)*dif;        
    vec4 col = vec4( mix( vec3(1.0,0.95,0.8), vec3(0.25,0.3,0.35), den ), den );
    col.xyz *= lin;
    col.xyz = mix( col.xyz, bgcol, 1.0-exp(-0.003*t*t) );
    // front to back blending    
    col.a *= 0.2;
    col.rgb *= col.a;
    return sum + col*(1.0-sum.a);
}

#define CHECK_IN_BOX(p) \
if(p.x < 0.0 || p.x > 1.0\
 ||p.y < 0.0 || p.y > 1.0\
 ||p.z < 0.0 || p.z > 1.0)\
 break;

void main()
{
    vec3 bgColor = vec3(0.8,0.0,0.4);
    vec3 lightColor = vec3(0.6, 0.8, 0.7);
    vec3 lightDir = normalize(vec3(5,5,5));

    // 射线的起点和方向
    eyeRay.o = viewPos;
    // 由于采样时使用的纹理坐标，
    // 而这个3D纹理坐标系和包围盒的本地坐标系是平行的（原点可能不一致），
    // 因此将ray的方向向量从世界坐标转换回本地坐标
    eyeRay.d = inverse(mat3(M)) * normalize(fragWorldPos - viewPos);

    float stepSize = 0.005;
    vec4 result = vec4(0.0);
    // 采样其实是从这个片段的纹理坐标开始的，然后沿着ray的方向步进
    vec3 p = textCoord;
    float steps = 0;
    for(int i = 0; i < 1000; i++, steps++){
        CHECK_IN_BOX(p);
        float dens = texture(noiseSampler, p).r;
        if(dens >= 0.01){
            float dif = clamp((dens - texture(noiseSampler, p + lightDir *0.5).r), 0.0, 1.0);
            result = integrate(result, dif, dens, bgColor, stepSize);
        }
        p += stepSize*eyeRay.d;
    }

    outColor = result;
    //outColor = vec4(fragWorldPos, 1.0);
}

shader代码和算法有非常直观的对应，也比较好理解。

接下来就是学习关于体绘制着色的算法了。相比于传统的几何面光照计算，体绘制复杂了许多，光线在粒子中的散射行为比较复杂，也有许多公式用来近似这个过程。看到书中的大把公式，我这数学白痴顿时有点腿软的感觉了…Anyway，体绘制是一个非常有意思的领域，能够绘制出很多美妙的效果。