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

最新推荐文章于 2024-06-24 20:31:48 发布

置顶 0小龙虾0

最新推荐文章于 2024-06-24 20:31:48 发布

阅读量4k

点赞数 9

分类专栏：流体模拟(SPH)

本文链接：https://blog.csdn.net/qq_39300235/article/details/110183675

版权

流体模拟(SPH) 专栏收录该内容

13 篇文章 194 订阅

订阅专栏

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 core

layout (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 core

out 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 core

out 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 core

layout (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 core

out 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 core

out 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);
    
}

完成渲染效果如下：