DirectX11海洋模拟实践

框架基于https://www.cnblogs.com/X-Jun/p/9028764.html

尽力用最精简的语言描述FFT海洋

总览

预计算阶段

计算静态海洋频谱值 h0

频谱计算阶段

在此阶段进行更新，引入时间变量

计算高度偏移频谱、水平偏移频谱、法线频谱

逆傅里叶变换阶段

进行二维傅里叶变换

得到偏移数据以及法线数据

最终处理阶段

对数据进行整合

理论

参考https://zhuanlan.zhihu.com/p/64414956

海面的IDFT模型

h(x⃗,t)=∑k⃗h~(k⃗,t)eik⃗⋅x⃗h(\vec{x}, t)=\sum_{\vec{k}} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}} h(x,t)=k∑h~(k,t)eik⋅x

此方程为二维逆傅里叶变换

其中k空间的两个分量定义为
kx=2πnL,n∈{−N2,−N2+1,…,N2−1}kz=2πmL,m∈{−N2,−N2+1,…,N2−1}\begin{array}{l} k_{x}=\frac{2 \pi n}{L}, n \in\left\{-\frac{N}{2},-\frac{N}{2}+1, \ldots, \frac{N}{2}-1\right\} \\ k_{z}=\frac{2 \pi m}{L}, m \in\left\{-\frac{N}{2},-\frac{N}{2}+1, \ldots, \frac{N}{2}-1\right\} \end{array} kx=L2πn,n∈{−2N,−2N+1,…,2N−1}kz=L2πm,m∈{−2N,−2N+1,…,2N−1}

海洋频谱

菲利普频谱（Phillips spectrum）
h~(k⃗,t)=h~0(k⃗)eiω(k)t+h~0∗(−k⃗)e−iω(k)t\tilde{h}(\vec{k}, t)=\tilde{h}_{0}(\vec{k}) e^{i \omega(k) t}+\tilde{h}_{0}^{*}(-\vec{k}) e^{-i \omega(k) t} h~(k,t)=h~0(k)eiω(k)t+h~0∗(−k)e−iω(k)t
其中
h~0∗=conj(h~0)\tilde{h}_{0}^{*} = conj(\tilde{h}_{0}) h~0∗=conj(h~0)

即其共轭复数
k=∣k⃗∣ω(k)=gkk = | \vec{k}| \\\omega(k)=\sqrt{g k} k=∣k∣ω(k)=gk
更准确来说是（h是海洋深度）
ω(k)=gktanh⁡(kh)\omega(k)=\sqrt{g k \tanh (k h)} ω(k)=gktanh(kh)

g是重力常数
h~0(k⃗)=12(ξr+iξi)Ph(k⃗)\tilde{h}_{0}(\vec{k})=\frac{1}{\sqrt{2}}\left(\xi_{r}+i \xi_{i}\right) \sqrt{P_{h}(\vec{k})} h~0(k)=21(ξr+iξi)Ph(k)
其中 ξr 和 ξi 是相互独立的随机数，均服从均值为0，标准差为1的正态分布。
Ph(k⃗)=Ae−1/(kL)2k4∣k⃗⋅w⃗∣2P_{h}(\vec{k})=A \frac{e^{-1 /(k L)^{2}}}{k^{4}}|\vec{k} \cdot \vec{w}|^{2} Ph(k)=Ak4e−1/(kL)2∣k⋅w∣2

风向：∣w⃗∣2=1风向：|\vec{w}|^{2} = 1 风向：∣w∣2=1

L=V2gL=\frac{V^{2}}{g} L=gV2

V为风速

具体含义如下

法线

因为求导的线性性质，h的梯度值仍然是IDFT形式，可以用IDFT计算
∇h(x⃗,t)=∑k⃗ik⃗h~(k⃗,t)eik⃗⋅x⃗\nabla h(\vec{x}, t)=\sum_{\vec{k}} i \vec{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}} ∇h(x,t)=k∑ikh~(k,t)eik⋅x

N⃗=normalize ((0,1,0)−(∇hx(x⃗,t),0,∇hz(x⃗,t)))=normalize (−∇hx(x⃗,t),1,−∇hz(x⃗,t))\begin{array}{l} \vec{N}=\text { normalize }\left((0,1,0)-\left(\nabla h_{x}(\vec{x}, t), 0, \nabla h_{z}(\vec{x}, t)\right)\right) \\ =\text { normalize }\left(-\nabla h_{x}(\vec{x}, t), 1,-\nabla h_{z}(\vec{x}, t)\right) \end{array} N= normalize ((0,1,0)−(∇hx(x,t),0,∇hz(x,t)))= normalize (−∇hx(x,t),1,−∇hz(x,t))

Gerstner wave

https://www.bilibili.com/video/BV1E64y1D78T?spm_id_from=333.880.my_history.page.click&vd_source=0cfebf737cf740b1e5610209f09e99a8Gerstner wave进行了水平偏移，压缩海浪形成尖角

P=(x⃗+∑QiAiDi⃗⋅x⃗cos⁡Φi∑Aisin⁡Φi)\mathbf{P}=\left(\begin{array}{l} \vec{x}+\sum\ Q_{i} A_{i} \vec{\mathbf{D}_{i}} \cdot \vec{x} \cos \Phi_{i} \\ \sum A_{i} \sin \Phi_{i} \end{array}\right)\\ P=(x+∑ QiAiDi⋅xcosΦi∑AisinΦi)
Φi=wiDi⃗⋅x⃗+φit\Phi_{i} = w_{i} \vec{\mathbf{D}_{i}} \cdot\vec{x}+\varphi_{i} t Φi=wiDi⋅x+φit

在IDFT中的形式为

D⃗(x⃗,t)=∑k⃗−ik⃗kh~(k⃗,t)eik⃗⋅x⃗x⃗′=x⃗+λD⃗(x⃗,t)\begin{array}{l} \vec{D}(\vec{x}, t)=\sum_{\vec{k}}-i \frac{\vec{k}}{k} \tilde{h}(\vec{k}, t) e^{i \vec{k} \cdot \vec{x}} \\ \vec{x}^{\prime}=\vec{x}+\lambda \vec{D}(\vec{x}, t) \end{array} D(x,t)=∑k−ikkh~(k,t)eik⋅xx′=x+λD(x,t)

快速傅里叶变化 FFT

参见https://www.bilibili.com/video/BV1Y7411W73U?spm_id_from=333.337.search-card.all.click

这方面的知识已经有很多讲解，这里不做详解

对图像进行横纵FFT/IFFT即可得到二维FFT
IFFT2=IFFT（IFFT（Matrixrol）col）IFFT2 = IFFT（IFFT（Matrix_{rol}）_{col}） IFFT2=IFFT（IFFT（Matrixrol）col）

HLSL实现

注：GPU上将CPU中一起执行的剪刀操作拆分到两个线程中并且使用ping-pong buffer加速计算

因为HLSL不支持float2类型的随机访问贴图，所以将两个频谱打包进行傅里叶变换（也加速了傅里叶变换的过程）

//#define FFT_SIZE_256#if defined(FFT_SIZE_512)
#define SIZE 512
#define LOG_SIZE 9
#elif defined(FFT_SIZE_256)
#define SIZE 256
#define LOG_SIZE 8
#elif defined(FFT_SIZE_128)
#define SIZE 128
#define LOG_SIZE 7
#else
#define SIZE 64
#define LOG_SIZE 6
#endifstatic uint Size = SIZE;#ifdef FFT_ARRAY_TARGET
RWTexture2DArray<float4> Target;
#else
RWTexture2D<float4> g_Target;
#endifcbuffer Params
{uint TargetsCount;uint Direction;//booluint Inverse; //booluint Scale; //booluint Permute; //bool
};//两组数据进行打包输入
groupshared float4 buffer[2][SIZE];float2 ComplexMult(float2 a, float2 b)
{return float2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}void ButterflyValues(uint step, uint index, out uint2 indices, out float2 twiddle)
{const float twoPi = 6.28318530718;uint b = Size >> (step + 1);uint w = b * (index / b);uint i = (w + index) % Size;sincos(-twoPi / Size * w, twiddle.y, twiddle.x);if (Inverse)twiddle.y = -twiddle.y;indices = uint2(i, i + b);
}float4 DoFft(uint threadIndex, float4 input)
{buffer[0][threadIndex] = input;GroupMemoryBarrierWithGroupSync();bool flag = false;[unroll(LOG_SIZE)]for (uint step = 0; step < LOG_SIZE; step++){uint2 inputsIndices;float2 twiddle;ButterflyValues(step, threadIndex, inputsIndices, twiddle);float4 v = buffer[flag][inputsIndices.y];//一次执行两个计算//ak = a0 + w*a1  buffer[!flag][threadIndex] = buffer[flag][inputsIndices.x]+ float4(ComplexMult(twiddle, v.xy), ComplexMult(twiddle, v.zw));flag = !flag;GroupMemoryBarrierWithGroupSync();}return buffer[flag][threadIndex];
}[numthreads(SIZE, 1, 1)]
void CS(uint3 id : SV_DispatchThreadID)
{uint threadIndex = id.x;uint2 targetIndex;if (Direction)targetIndex = id.yx;elsetargetIndex = id.xy;#ifdef FFT_ARRAY_TARGETfor (uint k = 0; k < TargetsCount; k++){Target[uint3(targetIndex, k)] = DoFft(threadIndex, Target[uint3(targetIndex, k)]);}
#elseg_Target[targetIndex] = DoFft(threadIndex, g_Target[targetIndex]);
#endif
}

项目GitHub
https://github.com/StellarWarp/FFT-Ocean