【体绘制】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个：

哪些片段需要raymarching。

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，体绘制是一个非常有意思的领域，能够绘制出很多美妙的效果。