SPH（光滑粒子流体动力学）流体模拟实现七：屏幕空间流体渲染（SSF）

之前都是用的Marching Cube重建流体表面的方法。近来为了做对比实验，尝试了屏幕空间流体渲染的方法。发现屏幕空间的方法不仅在效率上有极大的提升，效果也还是非常可观的，先上两张图：

本文方法来自NYIDIA的文章“Screen Space Fluid Rendering for Games”。本文的流体粒子运动模拟方法使用的是前文的PBF方法，我们已经根据PBF计算完成粒子运动。利用粒子完成屏幕空间渲染的主题思路如下：

利用点精灵的方式所有粒子。
根据绘制的点精灵记录其深度，其深度计算将点精灵视为球体，保存至深度缓存。
根据绘制的点精灵记录其厚度，其厚度计算同样将点精灵视为球体，保存至厚度缓存。厚度不能和深度使用MRT同时输出，因为厚度需要累加，因此我们要关闭深度测试，利用alpha加法混合，累加厚度。
利用双边高斯滤波对深度缓存进行滤波，用普通高斯滤波对厚度缓存滤波。并将其用MRT输出至两张帧缓存中。
利用过滤完的深度重建法线，并使用一张帧缓存保存法线。
使用法线以及厚度缓存，进行水的表面渲染。

其流程图如下：

该图为论文中的原图，我个人认为Thickness Image到Surface Shader这部之间应该还需要一步高斯滤波，并且可以放在Depth Smoothing步骤中同时进行。

接下来我们看具体实现步骤。

绘制点精灵

我们在渲染粒子时，使用几何着色器，在眼空间中将粒子渲染成正方形方片。我们之后所有的计算都在眼空间中进行。然后在片段着色器中，裁剪掉超出圆半径的片段，示例如下：

我们的顶点着色器和几何着色器的代码如下：

//顶点着色器----------------------------------------------------
#version 410 core
layout (location = 0) in vec3 position;uniform mat4 view;
uniform mat4 projection;
uniform mat4 model;
uniform float pointRadius;out float PointRadiusView;
out mat4 Projection;void main()
{vec4 viewPos=view*model*vec4(position,1.0f);PointRadiusView=length(view*model*vec4(pointRadius,0,0,0));Projection=projection;gl_Position=viewPos;
}//几何着色器----------------------------------------------------
#version 410 corelayout (points) in;
layout (triangle_strip,max_vertices = 6) out;in float PointRadiusView[];
in mat4 Projection[];
in mat4 View[];out float disCP;
out vec3 viewCenterPos;
out vec3 fragPos;
out mat4 Projectionf;
void main()
{Projectionf=Projection[0];disCP=PointRadiusView[0];viewCenterPos=gl_in[0].gl_Position.xyz;fragPos=vec3((gl_in[0].gl_Position+vec4(-PointRadiusView[0],-PointRadiusView[0],0,0)));gl_Position = Projection[0]*vec4(fragPos,1.0);EmitVertex();fragPos=vec3((gl_in[0].gl_Position+vec4(PointRadiusView[0],-PointRadiusView[0],0,0)));gl_Position = Projection[0]*vec4(fragPos,1.0);EmitVertex();fragPos=vec3((gl_in[0].gl_Position+vec4(-PointRadiusView[0],PointRadiusView[0],0,0)));gl_Position = Projection[0]*vec4(fragPos,1.0);EmitVertex();fragPos=vec3((gl_in[0].gl_Position+vec4(PointRadiusView[0],PointRadiusView[0],0,0)));gl_Position = Projection[0]*vec4(fragPos,1.0);EmitVertex();EndPrimitive();
}

生成深度缓存和厚度缓存

我们之前说过，由于厚度缓存需要关闭深度缓存，并采用加法混合的方式计算，因此厚度和深度缓存需要分布计算输出。在输出深度缓存时，我们需要注意圆片的不同位置深度是不同的，中间厚两边薄，因此中心的深度值较小。我们生成深度缓存片段着色器的代码如下：

#version 410 coreout float FVDepth;in float disCP;
in vec3 viewCenterPos;
in vec3 fragPos;
in mat4 Projectionf;uniform float NEAR;
uniform float FAR;void main(){float discp=distance(viewCenterPos,fragPos);vec4 fColor;if(discp>disCP){discard;}float height=sqrt(disCP*disCP-discp*discp);//深度float depthView=(fragPos.z+height);vec4 clip_space_pos=Projectionf*vec4(vec3(depthView),1.0);gl_FragDepth=(clip_space_pos.z/clip_space_pos.w)*0.5+0.5;FVDepth=-depthView;}

我们裁剪掉超出半径的部分，并利用圆公式，计算深度。我们这里记录眼空间下的深度，因为眼空间z值为负，我们为其取正，方便计算和可视化。我们这里用FVDepth记录线性深度，但仍给gl_FragDepth赋值透视投影变化后的深度，因为这样才能保证深度测试的进行。

厚度缓存需要关闭深度测试，并开启混合，我们的设置代码如下：

glDisable(GL_DEPTH_TEST);
glEnable(GL_BLEND);
glBlendFunc(GL_ONE, GL_ONE);//加法混合
render(shader,wiThick);//将厚度绘制在厚度缓存中
glDisable(GL_BLEND);
glEnable(GL_DEPTH_TEST);

片段着色器代码为：

#version 410 coreout float FragColor;in float disCP;
in vec3 viewCenterPos;
in vec3 fragPos;uniform float NEAR;
uniform float FAR;void main(){float discp=distance(viewCenterPos,fragPos);vec4 fColor;if(discp>disCP){discard;}float height=sqrt(disCP*disCP-discp*discp);FragColor=2.0*height;
}

我们分别可视化深度缓存和厚度缓存，效果如下：

我们实际的深度和厚度是线性的，且大于1（屏幕输出的颜色会被钳制为0～1），我们根据摄像机远平面，对深度和厚度进行缩放后才能可视化。我们这里仅仅为了观察才对其缩放，实际帧缓存中是未缩放的值，帧缓存中保存的是和HDR一致的浮点型帧缓存，用于保存大于1的值。

我们可以发现，其具有强烈的颗粒感，因此我们之后需要对其做平滑处理。我们观察一下未平滑处理生成的最后效果：

滤波处理

使用MRT（多渲染目标）方法可以在一个着色器内平滑并输出深度和厚度缓存。对厚度缓存，我们直接采用高斯模糊的方法。但是深度缓存并不能直接使用高斯模糊，因为我们需要用深度信息重建法线，如果直接使用高斯模糊过滤深度信息会产生如下效果：

浮于上方的粒子会被平滑至同一表面内，因此我们需要值域（深度）信息。我们这里就使用双边过滤（Bilateral Filter）。双边滤波思维导图如下：

双边滤波公式如下：

其中 $W_{p}$ 是权重和， $G_{\sigma }_{s}$ 为空域高斯权值， $G_{\sigma }_{r}$ 为值域高斯权值。其不仅在空域对图像进行滤波，还考虑了值域的影响，我们的深度值正是如此，深度差距较大的值，经过双边滤波，依旧能保存其原本的深度信息。我们的滤波代码如下：

#version 410 corelayout (location = 0) out float Fdepth;
layout (location = 1) out float Fthickness;in vec2 aTexCoords;uniform sampler2D depthMap;
uniform sampler2D thickMap;
uniform vec2 blurDir;
uniform float filterRadius;
uniform float spatialScale;
uniform float rangeScale;const float NEAR=0.1f;
const float FAR=200.0f;void main(){vec2 depTexSize=textureSize(depthMap,0);vec2 thickTexSize=textureSize(thickMap,0);float depth=texture(depthMap,aTexCoords).r;float thickness=texture(thickMap,aTexCoords).r;float sumDep=0.0f,sumThick=0.0f;float wsumDep=0.0f,wsumThick=0.0f;for (float x=-filterRadius; x<=filterRadius; x+=1.0) {float sampleDepth=texture(depthMap,aTexCoords+x*blurDir).r;float sampleThick=texture(thickMap,aTexCoords+x*blurDir).r;//空域float r=x*spatialScale;float w=exp(-r*r);//值域float r2Dep=(sampleDepth-depth)*rangeScale;float gDep=exp(-r2Dep*r2Dep);//深度采用双边滤波，厚度只使用高斯滤波sumDep+=sampleDepth*w*gDep;sumThick+=sampleThick*w;wsumDep+=w*gDep;wsumThick+=w;}sumDep/=wsumDep;sumThick/=wsumThick;Fdepth=sumDep;Fthickness=(sumThick);
}

我们滤波后的效果如下：

为了提升速度，我们这里还采用了1D的多次滤波方法，代码如下：

bool firstBlur=1;
for (int i=0; i<blurNumber;i++) {//绘制模糊glBindFramebuffer(GL_FRAMEBUFFER,blurFrame[0]);glViewport(0, 0, downSample*winWidth, downSample*winHeight);glClearColor(0.0f, 0.0f, 0.0f, 1.0f);glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);wiBlur.use();if(firstBlur){wiBlur.setInt("depthMap", 1);wiBlur.setInt("thickMap", 2);firstBlur=0;}else{wiBlur.setInt("depthMap", 6);wiBlur.setInt("thickMap", 7);}wiBlur.setVec2("blurDir", glm::vec2(1.0/(downSample*winWidth),0.0f));wiBlur.setMat4("invP", glm::inverse(projection));renderWinPlane(wiBlur);wiBlur.close();glBindFramebuffer(GL_FRAMEBUFFER,blurFrame[1]);glViewport(0, 0, downSample*winWidth, downSample*winHeight);glClearColor(0.0f, 0.0f, 0.0f, 1.0f);glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);wiBlur.use();wiBlur.setInt("depthMap", 4);wiBlur.setInt("thickMap", 5);wiBlur.setVec2("blurDir", glm::vec2(0.0f,1.0/(downSample*winHeight)));renderWinPlane(wiBlur);wiBlur.close();
}

重构法线

我们已经得到双边滤波后的深度缓存了，我们接下来重构法线，我们这里的方法也很简单，直接用深度还原片眼空间的位置，然后利用位置求解法线，代码如下：

#version 410 coreout vec4 FragColor;in vec2 aTexCoords;uniform sampler2D depthMap;
uniform float NEAR;
uniform float FAR;
uniform float aspect;//近平面高:宽
uniform float nearHeight;//近平面高const float epo=1e-7;//裁剪空间转换为眼空间
vec3 uvToEye(vec2 texCoord,float depth){
//    vec4 clipSpacePos;
//    clipSpacePos.xy = texCoord.xy * 2.0f - 1.0f;
//    clipSpacePos.z = depth*2.0-1.0;
//    clipSpacePos.w = 1;
//    clipSpacePos = invMVP * clipSpacePos;
//    return clipSpacePos.xyz / clipSpacePos.w;vec2 deltaUV=(2.0*texCoord-vec2(1.0))*vec2(aspect,1.0);//计算近平面的平移向量vec2 deltaView=0.5*nearHeight*deltaUV*depth/NEAR;return vec3(vec2(deltaView),-depth);
}vec3 getEyePos(sampler2D depthTex,vec2 texCoord){float depth=(texture(depthTex,texCoord).r);return uvToEye(texCoord,depth);
}void main(){vec2 texSize=vec2(1.0/textureSize(depthMap,0).s,1.0/textureSize(depthMap,0).t);//深度float depthV=(texture(depthMap,aTexCoords).r);if (depthV>=FAR-1.0) {discard;}//获得当前点的view空间位置vec3 posEye=uvToEye(aTexCoords,depthV);//计算微分vec3 ddu=getEyePos(depthMap,aTexCoords+vec2(texSize.x,0))-posEye;vec3 ddub=posEye-getEyePos(depthMap,aTexCoords-vec2(texSize.x,0));if (abs(ddu.z)>abs(ddub.z)) {ddu=ddub;}vec3 ddv=getEyePos(depthMap,aTexCoords+vec2(0,texSize.y))-posEye;vec3 ddvb=posEye-getEyePos(depthMap,aTexCoords-vec2(0,texSize.y));if (abs(ddv.z)>abs(ddvb.z)) {ddv=ddvb;}//计算法线vec3 N=cross(ddu,ddv);N=normalize(N);//N=vec3(inverse(transpose(inverse(mat3(invMVP))))*N);FragColor=vec4(N,1.0);
}

这里由于我保存的深度是线性眼空间的深度，非透视投影计算后的深度，因此我没用投影矩阵的逆矩阵求解片段眼空间的位置，而是使用三角形的相似去求解眼空间的位置：

根据上图，我们只需要近平面的高度和宽高比就能计算片段在眼空间的位置了。近平面的高度为：

float nearHeight=2*NEAR_PLANE*glm::tan(glm::radians(m_camera.Zoom)/2.0);

之后我们对u和v方向求微分，两个微分向量的叉乘就是我们的法线了，值得注意的是我们使用两个相反的方向进行微分计算，选择绝对值分量更小的。这是因为如果使用单一的微分，会导致法线在物体边缘突变，而产生错误，示意图如下：

我们看到未经处理的左图，粒子和粒子之间有明显的黑色边缘，而处理后的右图则明显改善。

我们可视化法线，效果如下：

我们发现双边滤波对表面有着极大的平滑提升。

着色计算

所有我们需要的信息我们都获取了，我们接下来便可以进行着色计算了。我们这里只考虑反射部分和折射部分。折射部分又考虑厚度对透射率的影响。我们使用近似菲涅尔方程，计算反射分量：

$F_{Schlick}=F_{0}+(1-F_{0})(1-h\cdot v)^5$

反射部分使用反射向量索引天空盒的方式。

折射部分使用Beer定律计算透射率：

$T_{r}=e^{-kd}$

其中d为厚度，k为透射率控制常量。透射部分，我们根据折射采样2D背景图（该图被我们事先渲染至帧缓存内），折射采样公式为：

texture(bg2D,texcoords-N.xy*thickness*k)

其中N为法线，thickness为厚度，k为折射标量（越大折射越明显）。未能透射的部分，我们假定返回默认水体颜色，水体的着色代码如下：

#version 410 coreout vec4 FragColor;in vec2 aTexCoords;uniform float NEAR;
uniform float FAR;
uniform float aspect;//近平面高:宽
uniform float nearHeight;//近平面高
uniform float shininess;//高光参数
uniform sampler2D depthMap;
uniform sampler2D thicknessMap;
uniform sampler2D normalMap;
uniform samplerCube cubeMap;
uniform sampler2D cubeMap2D;
uniform vec3 lightPos;//指向光源（眼空间）
uniform mat4 view;const vec3 diffuseColor=vec3(0.5,0.5,0.5);
const vec3 specularColor=vec3(0.5f);
const vec3 waterRefrectColor=vec3(0.05,0.5,0.8);
const float epo=1e-2;
const vec3 waterF0=vec3(0.15f);
const float refrectScale=2.0;
const float waterK=0.05;//裁剪空间转换为眼空间
vec3 uvToEye(vec2 texCoord,float depth){vec2 deltaUV=(2.0*texCoord-vec2(1.0))*vec2(aspect,1.0);//计算近平面的平移向量vec2 deltaView=0.5*nearHeight*deltaUV*depth/NEAR;return vec3(vec2(deltaView),-depth);
}vec3 culFresnel(vec3 f0,float cosTheta){return f0+(1.0-f0)*pow(1.0-cosTheta,5.0);
}//眼空间光照计算
void main(){vec3 lightDir=normalize((mat3(view)*lightPos));float depth=texture(depthMap,aTexCoords).r;if (depth>=FAR-epo) {discard;}vec2 texSize=vec2(1.0/textureSize(normalMap,0).s,1.0/textureSize(normalMap,0).t);vec3 normal=texture(normalMap,aTexCoords).xyz;float thickness=texture(thicknessMap,aTexCoords).r;//phong漫反射（钳制后）vec3 ambient=diffuseColor*0.1;vec3 diffuse=ambient+diffuseColor*(max(dot(normal,lightDir),0.0));vec3 eyePos=uvToEye(aTexCoords,depth);vec3 veiwDir=normalize(-eyePos);//blinn-phong镜面反射vec3 halfDir=normalize(veiwDir+lightDir);vec3 specular=specularColor*pow(max(dot(normal, halfDir), 0.0), shininess);//fresnel反射分量vec3 Rfresnel=culFresnel(waterF0,max(dot(normal,veiwDir),0.0));vec3 RreflectDir=normalize(mat3(inverse(view))*reflect(-veiwDir,normal));//世界空间vec3 cubeReflect=texture(cubeMap,RreflectDir).rgb;//折射分量（t=p-b(N*P)）vec3 Refrect=texture(cubeMap2D,aTexCoords-normal.xy*thickness/FAR*refrectScale).xyz;//透射率float transparency=exp(-thickness*waterK);Refrect=mix(waterRefrectColor,Refrect,transparency);FragColor=vec4(mix(Refrect,cubeReflect,Rfresnel),1.0);}

完成渲染效果如下：

整体效果还是比较好的，不过双边滤波还有有一些比较明显的缺点。比如单体粒子的边缘会过于扁平，导致我们看到上面的视频中飞溅的水滴边界明显。对比我们Marching Cube的重建表面效果如下：

SSF方法：

各向异性后的Marching Cube方法：

我们可以看到屏幕空间已经能带来极大的平滑效果，然而其不足就是单一粒子的表现不足。我接下来可能会去尝试更好的滤波方式，使用自适应卷积核半径，以及边界处理方式来得到更好的屏幕空间平滑效果