Unity《游戏特效编程》作业1：纹理扭曲动画+UV坐标扰动
本文禁止转载
B站：Heskey0
《图形特效编程》第一次课程作业
本次作业包括两个部分:
- 课堂作业1：通过时间参数扰动纹理坐标，实现纹理的扭曲动画
- 课堂作业2：通过时间参数扰动纹理坐标，使用水面纹理实现水流动画
最终效果
作业一：（效果）
https://www.bilibili.com/video/BV12g411S7m3
作业一：（代码）
https://github.com/Heskey0/Unity-Shader-Graphics/tree/main/ProcedureAnim-WaterCaustics
作业二：（效果）
https://www.bilibili.com/video/BV12B4y1g7xV
作业二：（代码）
https://github.com/Heskey0/Unity-Shader-Graphics/tree/main/FluidSim-VotexConfinement
作业一：纹理扭曲动画
作业二：UV坐标扰动
1. 函数波动现象
首先来看一些函数图像：
\(\sin(x)\):
\(\cos(sin(x))\):
\(cos(sin(cos(x)))\):
很容易发现：sin与cos互相嵌套，函数值域会变为 \((0,1)\) 并且其导数的频率也会降低。
如果在sin与cos相互嵌套的过程中，加入偏移(例如 \(cos(sin(cos(x))+t)\))，则随着t的增加整个函数图像会在 \((0,1)\) 之间周期性波动。例如：
我们可以利用这个波动现象，来模拟水的波纹。
2. 思路分析：
我们先设计一个sin和cos不断嵌套的函数，例如：
\[\cos(\sin(x)+\cos(x))+\sin(\sin(x)+\cos(x))
\]
然后再嵌套过程中加入偏移 \(t\)，例如：
\[\cos(\sin(x+t)+\cos(x+t)+t)+\sin(\sin(x+t)+\cos(x+t)+t)
\]
这样一来，随着 \(t\) 的增长，函数图像就会在 \((0,1)\) 之间一直波动。但要模拟水的波纹，一个波是不够的，我们可以将多个波叠加起来。那么该如何设计多个波？
这里我使用的波，它们sin和cos嵌套的层数不同，但是嵌套的方式不变。考虑函数的自变量是二维向量（这里我使用的是uv坐标） \(p=(x,y)\)，然后时间记为 \(t\)，则函数组 \(i=(x,y)\) 可以写成：
\[i.x=p.x+cos(t-i.x)+sin(t+i.y)\\
i.y=p.y+sin(t-i.y)+cos(t+i.x)
\]
只需要不断更新 \(i\) 的值，就会产生sin和cos不断嵌套的波。
但这样设计函数，会导致函数的值域超出 \((-1,1)\)，所以我们引入一个新的向量：
\[v=\pmatrix{
\frac{p.x}{sin(i.x+t)}\\
\frac{p.y}{cos(i.y+t)}
}
\]
然后取其模长的倒数，累加称为最终的值：
\[c=\sum_{n=0}^{5}\frac{1}{inten*||v||}
\]
3. 代码实现：
Shader "Homework/01_1"
{
Properties
{
_WaveFrequency("Frequency", Range(0.1, 10)) = 3.5
_Speed("Speed", Range(0.1, 3)) = 0.5
_TAU("Scale", Range(1,10)) = 6.28318530718
_Inten("Brightness", Range(0.0005, 0.02)) = 0.005
_FOO("FOO", Range(1,500)) = 250.0
}
SubShader
{
Tags
{
"RenderPipeline" = "UniversalPipeline"
"LightMode"="UniversalForward"
}
Pass
{
HLSLPROGRAM
#pragma vertex vert
#pragma fragment frag
#include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/Core.hlsl"
#include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/Lighting.hlsl"
#include "Packages/com.unity.render-pipelines.core/ShaderLibrary/SpaceTransforms.hlsl"
#include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/Shadows.hlsl"
CBUFFER_START(UnityPerMaterial)
float _FOO;
float _WaveFrequency;
float _Speed;
float _TAU;
float _Inten;
CBUFFER_END
struct Attributes
{
float4 positionOS : SV_POSITION;
float3 normalOS : NORMAL;
float2 uv : TEXCOORD0;
};
struct Varyings
{
float3 positionWS : TEXCOORD0;
float4 positionCS : SV_POSITION;
float2 uv : TEXCOORD2;
float3 normalWS : TEXCOORD3;
float4 screenPos : TEXCOORD4;
};
Varyings vert(Attributes i)
{
Varyings o;
o.positionWS = TransformObjectToWorld(i.positionOS);
o.positionCS = TransformWorldToHClip(o.positionWS);
o.uv = i.uv;
o.normalWS = TransformObjectToWorldNormal(i.normalOS);
o.screenPos = ComputeScreenPos(o.positionCS);
return o;
}
float4 frag(Varyings v) : SV_Target
{
float3 o = float3(1,1,1);
float time = _Time.g * _Speed;
float2 uv = v.screenPos.xy / v.screenPos.w;
float2 p = uv*_TAU - _FOO;
float2 i = p;
float c = 1.0;
float inten = _Inten;
float MAX_ITER = 5;
float n = 0; //0~5
//0
float t = time * (1.0 - _WaveFrequency/float(n+1));
i = p + float2(cos(t-i.x)+sin(t+i.y), sin(t-i.y)+cos(t+i.x));
c += 1.0/(inten*length(float2(p.x / sin(i.x+t), p.y / cos(i.y+t))));
n++;
// 1
t = time * (1.0 - _WaveFrequency/float(n+1));
i = p + float2(cos(t-i.x)+sin(t+i.y), sin(t-i.y)+cos(t+i.x));
c += 1.0/length(float2(p.x / (sin(i.x+t)/inten), p.y / (cos(i.y+t)/inten)));
n++;
//2
t = time * (1.0 - _WaveFrequency/float(n+1));
i = p + float2(cos(t-i.x)+sin(t+i.y), sin(t-i.y)+cos(t+i.x));
c += 1.0/length(float2(p.x / (sin(i.x+t)/inten), p.y / (cos(i.y+t)/inten)));
n++;
//3
t = time * (1.0 - _WaveFrequency/float(n+1));
i = p + float2(cos(t-i.x)+sin(t+i.y), sin(t-i.y)+cos(t+i.x));
c += 1.0/length(float2(p.x / (sin(i.x+t)/inten), p.y / (cos(i.y+t)/inten)));
n++;
//4
t = time * (1.0 - _WaveFrequency/float(n+1));
i = p + float2(cos(t-i.x)+sin(t+i.y), sin(t-i.y)+cos(t+i.x));
c += 1.0/length(float2(p.x / (sin(i.x+t)/inten), p.y / (cos(i.y+t)/inten)));
c /= float(MAX_ITER);
c = 1.17-pow(c, 1.4);
o = float3(pow(abs(c), 8.0)*float3(1,1,1));
o = clamp(o + float3(0.0, 0.35, 0.5), 0.0, 1.0);
return float4(o,1);
}
ENDHLSL
}
}
}
作业2：UV坐标扰动
1. 思路分析
目标：使用无散流场的解算结果（速度场）来扭曲纹理
参考：
- 代码：胡大佬的主页 https://www.zhihu.com/people/hua-la-la-53-15
- 理论：Bridson和Muller的course：https://www.cs.ubc.ca/~rbridson/
步骤：
- 使用半欧拉法，求解对流速度（简单处理流固耦合）
- 施加外力（重力，粘性力等）
- 计算速度场的散度
- 修复涡量
- 求解压强的泊松方程
- 使用压强的梯度，将速度场投影为无散
- 使用速度场扭曲纹理
2. 代码实现
2.1 半欧拉法求解对流速度
\[q^{n+1}=advect(v^n,\Delta t,q^n)
\]
将RT的纹理查询设置为双线性插值，就不需要自己写插值算法：
VelocityRT.filterMode = FilterMode.Bilinear;
使用半欧拉法：（在backtrace的过程中，还可以使用Runge-Kutta进一步调高时间积分的精度）
// backtrace
float2 last_pos = i.uv - dt*tex2D(VelocityTex, i.uv);
// semi-Lagrangian: bi-linear interpolation
col.xy = tex2D(QuantityTex, last_pos).xy; // advection: q = q
简单处理下流固耦合，固体速度场直接设置为0：
if(tex2D(BlockTex, i.uv).x > 0.99f)col.xy = float2(0.0f, 0.0f);
2.2 施加外力
对于烟雾模拟，可以忽略重力。对于水的模拟，可以忽略粘性（在模拟过程中存在数值耗散）
2.3 速度场的散度
\[d=\nabla\cdot u
\]
float Top = tex2D(_VelocityTex, i.uv + float2(0.0f, _VelocityTex_TexelSize.y)).y;
float Bottom = tex2D(_VelocityTex, i.uv + float2(0.0f, -_VelocityTex_TexelSize.y)).y;
float Right = tex2D(_VelocityTex, i.uv + float2(_VelocityTex_TexelSize.x, 0.0f)).x;
float Left = tex2D(_VelocityTex, i.uv + float2(-_VelocityTex_TexelSize.x, 0.0f)).x;
float divergence = 0.5f * (Right - Left + Top - Bottom); // delta x = 1
2.4 修复涡量
注：这一部分是自己实现的
由于模拟过程存在数值耗散，所以额外施加一个生成涡量的力，推动涡以维持涡的存在。
\[f_{\textrm{conf}}=\epsilon\Delta x(\vec N\times\vec\omega)
\]
需要注意的是，\(\vec N\) 的计算过程中需要对涡量的梯度做归一化，归一化时要在分母上添加一个值防止除0：
\[\vec N_{i,j,k}=\frac{\nabla|\vec\omega|_{i,j,k}}{||\nabla|\vec\omega|_{i,j,k}||+10^{-20}}
\]
计算涡量的代码：
float4 frag (v2f i) : SV_Target
{
float4 col = float4(0,0,0,1);
float2 L = tex2D(VelocityTex, i.uv - float2(VelocityTex_TexelSize.x, 0.0)).xy;
float2 R = tex2D(VelocityTex, i.uv + float2(VelocityTex_TexelSize.x, 0.0)).xy;
float2 B = tex2D(VelocityTex, i.uv - float2(0.0, VelocityTex_TexelSize.y)).xy;
float2 T = tex2D(VelocityTex, i.uv + float2(0.0, VelocityTex_TexelSize.y)).xy;
col.xy = float2(0, 0);
col.z = (R.y-L.y - (T.x-B.x)) * 0.5f; // along the z axis
return col;
}
速度更新的代码：
float4 frag (v2f i) : SV_Target
{
float4 col = float4(0,0,0,1);
float L = tex2D(VorticityTex, i.uv - float2(VorticityTex_TexelSize.x, 0.0)).z;
float R = tex2D(VorticityTex, i.uv + float2(VorticityTex_TexelSize.x, 0.0)).z;
float B = tex2D(VorticityTex, i.uv - float2(0.0, VorticityTex_TexelSize.y)).z;
float T = tex2D(VorticityTex, i.uv + float2(0.0, VorticityTex_TexelSize.y)).z;
float C = tex2D(VorticityTex, i.uv).z;
float2 N = float2(T-B, L-R);
float2 force = curl_strength * C * N/(length(N)+0.01);
float2 velocity = tex2D(VelocityTex, i.uv).xy;
col.xy = velocity + force * dt;
col.z = 0.0;
return col;
}
2.5 解泊松方程
因为使用 Graphics.Blit() 在RT之间传递物理场，所以在Unity里面不方便做多层循环。尝试了MGPCG方法，写到一半发现CPU和GPU之间传递的变量太多，且不方便切换solver，于是沿用了胡大佬的Jacobi作为standalone Solver。
2.6 速度场投影
\[\vec u^{n+1}=\vec u-\Delta t\frac1\rho\nabla p
\]
代码如下
velocity.xy -= float2(R - L,T - B);
2.7 显示
为了简便，没有传递Dye场，而直接展示速度场
float4 frag (v2f i) : SV_Target
{
float4 col = tex2D(_MainTex, i.uv);
float x = (col.x + 1.0) / 2;//value的值在-1到1之间，x的值在0到1之间
if (x < 0.25) col = float4(0.0, 4.0 * x, 1.0, 1.0f);
else if (x < 0.5) col = float4(0.0, 1.0, 1.0 + 4.0 * (0.25 - x), 1.0f);
else if (x < 0.75) col = float4(4.0 * (x - 0.5), 1.0, 0.0, 1.0f);
else col = float4(1.0, 1.0 + 4.0 * (0.75 - x), 0.0,1.0f);
if (tex2D(BlockTex, i.uv).x > 0.9f)col = float4(0.0f, 0.0f, 0.0f, 1.0f);
return col;
}
