## 主要内容一览

1. Gerstner(Trochoidal) Wave
2. Phillips Spectrum
3. 基于GPU的快速傅立叶变换
4. 海面浪端白沫(Whitecap)检测

## 目录

1. 引入 – Sinusoid正弦波.
2. 进化 – Gerstner Wave与其实现.
3. 介绍 – 海洋学统计模型.
4. 引入 – 完全基于CPU运算的海洋学统计模型实现.
5. 进化 – GPU上的FFT, 随机数与Implementation Notes.
6. 搬家 – 基于GPU的海洋学统计模型实现.
7. 后记
8. Renferences

## 进化 – Gerstner Wave与其实现

Gerstner Wave的解如下:

```inline void GerstnerLevelOne(out half3 offsets, out half3 normal, half3 vertex, half3 sVertex,
half amplitude, half frequency, half steepness, half speed,
half4 directionAB, half4 directionCD)
{
half3 offs = 0;

for(int i = 0;i < 5;i++)
{
offs.x += steepness * amplitude * steeps[i] * amps[i] * dir[i].x * cos(frequency * fs[i] * dot(sVertex.xz, dir[i]) + speeds[i] * frequency * fs[i] * _Time.y);
offs.z += steepness * amplitude * steeps[i] * amps[i] * dir[i].y * cos(frequency * fs[i] * dot(sVertex.xz, dir[i]) + speeds[i] * frequency * fs[i] * _Time.y);
offs.y += amplitude * amps[i] * sin(frequency * fs[i] * dot(sVertex.xz, dir[i].xy) + speeds[i] * frequency * fs[i] * _Time.y);
}

offsets = offs;

normal = half3(0, 1, 0);
for(int i = 0;i < 5;i++)
{
normal.x -= dir[i].x * frequency * fs[i] * amplitude * amps[i] * cos(dot(offs, frequency * fs[i] * dir[i]) + speeds[i] * frequency * fs[i] * _Time.y);
normal.z -= dir[i].y * frequency * fs[i] * amplitude * amps[i] * cos(dot(offs, frequency * fs[i] * dir[i]) + speeds[i] * frequency * fs[i] * _Time.y);
normal.y -= steepness * steeps[i] * frequency * fs[i] * amps[i] * amplitude * sin(dot(offs, frequency * fs[i] * dir[i]) + speeds[i] * frequency * fs[i] * _Time.y);
}
}```

Gerstner Wave在游戏中比较常用. 网络上有一个UE4的实现教程 [3], 写得非常好, 大家可以参考一下.

Gerstner Wave还有一个变种 — Skewed Trochoidal Wave, 可以用来模拟涌向岸边的潮水. 其原理可以在这里了解(似乎需要梯子).

## 引入 – 完全基于CPU运算的海洋学统计模型实现

```#region MonoBehaviours

private void Update()
{
timer += Time.deltaTime / tDivision;
EvaluateWaves(timer);
}

private void Awake()
{
filter = GetComponent();
mesh = new Mesh();
filter.mesh = mesh;
SetParams();
GenerateMesh();
}

#endregion```

```private void SetParams()
{
vertices = new Vector3[resolution * resolution];
indices = new int[(resolution - 1) * (resolution - 1) * 6];
normals = new Vector3[resolution * resolution];
vertConj = new Vector2[resolution * resolution];
verttilde = new Vector2[resolution * resolution];
vertMeow = new Vector3[resolution * resolution];//Meow ~
uvs = new Vector2[resolution * resolution];
}

private void GenerateMesh()
{
int indiceCount = 0;
int halfResolution = resolution / 2;
for (int i = 0; i < resolution; i++)
{
float horizontalPosition = (i - halfResolution) * unitWidth;
for (int j = 0; j < resolution; j++)
{
int currentIdx = i * (resolution) + j;
float verticalPosition = (j - halfResolution) * unitWidth;
vertices[currentIdx] = new Vector3(horizontalPosition + (resolution % 2 == 0? unitWidth / 2f : 0f), 0f, verticalPosition + (resolution % 2 == 0? unitWidth / 2f : 0f));
normals[currentIdx] = new Vector3(0f, 1f, 0f);
verttilde[currentIdx] = htilde0(i, j);
Vector2 temp = htilde0(resolution - i, resolution - j);
vertConj[currentIdx] = new Vector2(temp.x, -temp.y);
uvs[currentIdx] = new Vector2(i * 1.0f / (resolution - 1), j * 1.0f / (resolution - 1));
if (j == resolution - 1)
continue;
if (i != resolution - 1)
{
indices[indiceCount++] = currentIdx;
indices[indiceCount++] = currentIdx + 1;
indices[indiceCount++] = currentIdx + resolution;
}
if (i != 0)
{
indices[indiceCount++] = currentIdx;
indices[indiceCount++] = currentIdx - resolution + 1;
indices[indiceCount++] = currentIdx + 1;
}
}
}
mesh.vertices = vertices;
mesh.SetIndices(indices, MeshTopology.Triangles, 0);
mesh.normals = normals;
mesh.uv = uvs;
filter.mesh = mesh;
}```

```private Vector3 Displacement(Vector2 x, float t, out Vector3 nor)
{
Vector2 h = new Vector2(0f, 0f);
Vector2 d = new Vector2(0f, 0f);
Vector3 n = Vector3.zero;
Vector2 c, htilde_c, k;
float kx, kz, k_length, kDotX;
for (int i = 0; i < resolution; i++)
{
kx = 2 * PI * (i - resolution / 2.0f) / length;
for (int j = 0; j < resolution; j++)
{
kz = 2 * PI * (j - resolution / 2.0f) / length;
k = new Vector2(kx, kz);
k_length = k.magnitude;
kDotX = Vector2.Dot(k, x);
c = new Vector2(Mathf.Cos(kDotX), Mathf.Sin(kDotX));
Vector2 temp = htilde(t, i, j);
htilde_c = new Vector2(temp.x * c.x - temp.y * c.y, temp.x * c.y + temp.y * c.x);
h += htilde_c;
n += new Vector3(-kx * htilde_c.y, 0f, -kz * htilde_c.y);
if (k_length < EPSILON)
continue;
d += new Vector2(kx / k_length * htilde_c.y, -kz / k_length * htilde_c.y);
}
}
nor = Vector3.Normalize(Vector3.up - n);
return new Vector3(d.x, h.x, d.y);
}```

```private void EvaluateWaves(float t)
{
hds = new Vector2[resolution * resolution];

for (int i = 0; i < resolution; i++)
{
for (int j = 0; j < resolution; j++)
{
int index = i * resolution + j;
vertMeow[index] = vertices[index];
}
}
for (int i = 0; i < resolution; i++)
{
for (int j = 0; j < resolution; j++)
{
int index = i * resolution + j;
Vector3 nor = new Vector3(0f, 0f, 0f);
Vector3 hd = Displacement(new Vector2(vertMeow[index].x, vertMeow[index].z), t, out nor);
vertMeow[index].y = hd.y;
vertMeow[index].z = vertices[index].z - hd.z * choppiness;
vertMeow[index].x = vertices[index].x - hd.x * choppiness;
normals[index] = nor;
hds[index] = new Vector2(hd.x, hd.z);
}
}

Color[] colors = new Color[resolution * resolution];

for (int i = 0; i < resolution; i++)//写得并不正确,
{
for (int j = 0; j < resolution; j++)
{
int index = i * resolution + j;
Vector2 dDdx = Vector2.zero;
Vector2 dDdy = Vector2.zero;
if (i != resolution - 1)
{
dDdx = 0.5f * (hds[index] - hds[index + resolution]);
}
if (j != resolution - 1)
{
dDdy = 0.5f * (hds[index] - hds[index + 1]);
}
float jacobian = (1 + dDdx.x) * (1 + dDdy.y) - dDdx.y * dDdy.x;
Vector2 noise = new Vector2(Mathf.Abs(normals[index].x), Mathf.Abs(normals[index].z)) * 0.3f;
float turb = Mathf.Max(1f - jacobian + noise.magnitude, 0f);
float xx = 1f + 3f * Mathf.SmoothStep(1.2f, 1.8f, turb);
xx = Mathf.Min(turb, 1.0f);
xx = Mathf.SmoothStep(0f, 1f, turb);
colors[index] = new Color(xx, xx, xx, xx);
}
}
mesh.vertices = vertMeow;
mesh.normals = normals;
mesh.colors = colors;
}```

## 进化 – GPU上的FFT, 随机数与Implementation Notes

Stockham FFT的伪代码如下:

```Shader "Hidden/Mistral Water/Helper/Vertex/Stockham"
{
Properties
{
_Input ("Input Sampler", 2D) = "black" {}
_TransformSize ("Transform Size", Float) = 256
_SubTransformSize ("Log Size", Float) = 8
}

{
Pass
{
Cull Off
ZWrite Off
ZTest Off

CGPROGRAM

#include "FFTCommon.cginc"

#pragma fragment frag
#pragma multi_compile _HORIZONTAL _VERTICAL

uniform sampler2D _Input;
uniform float _TransformSize;
uniform float _SubTransformSize;

float4 frag(FFTVertexOutput i) : SV_TARGET
{
float index;

#ifdef _HORIZONTAL
index = i.texcoord.x * _TransformSize - 0.5;
#else
index = i.texcoord.y * _TransformSize - 0.5;
#endif

float evenIndex = floor(index / _SubTransformSize) * (_SubTransformSize * 0.5) + fmod(index, _SubTransformSize * 0.5) + 0.5;

#ifdef _HORIZONTAL
float4 even = tex2D(_Input, float2((evenIndex), i.texcoord.y * _TransformSize) / _TransformSize);
float4 odd  = tex2D(_Input, float2((evenIndex + _TransformSize * 0.5), i.texcoord.y * _TransformSize) / _TransformSize);
#else
float4 even = tex2D(_Input, float2(i.texcoord.x * _TransformSize, (evenIndex)) / _TransformSize);
float4 odd  = tex2D(_Input, float2(i.texcoord.x * _TransformSize, (evenIndex + _TransformSize * 0.5)) / _TransformSize);
#endif

float twiddleV = GetTwiddle(index / _SubTransformSize);
float2 twiddle = float2(cos(twiddleV), sin(twiddleV));
float2 outputA = even.xy + MultComplex(twiddle, odd.xy);
float2 outputB = even.zw + MultComplex(twiddle, odd.zw);

return float4(outputA, outputB);
}

ENDCG
}
}
}

// FFTCommon.cginc
//inline float GetTwiddle(float ratio)
//{
//	return -2 * PI * ratio;
//}```

```inline float UVRandom(float2 uv, float salt, float random)
{
uv += float2(salt, random);
return frac(sin(dot(uv, float2(12.9898, 78.233))) * 43758.5453);
}```

```inline float2 hTilde0(float2 uv, float r1, float r2, float phi)
{
float2 r;
float rand1 = UVRandom(uv, 10.612, r1);
float rand2 = UVRandom(uv, 11.899, r2);
rand1 = clamp(rand1, EPSILON, 1);
rand2 = clamp(rand2, EPSILON, 1);
float x = sqrt(-2 * log(rand1));
float y = 2 * PI * rand2;
r.x = x * cos(y);
r.y = x * sin(y);
return r * sqrt(phi / 2);
}```

## 搬家 – 基于GPU的海洋学统计模型实现

```/// Called in Start
private void RenderInitial()
{
Graphics.Blit(null, initialTexture, initialMat);
spectrumMat.SetTexture("_Initial", initialTexture);
heightMat.SetTexture("_Initial", initialTexture);
}

/// Called in Update
private void GenerateTexture()
{
float deltaTime = Time.deltaTime;
//使用两张RT确保参数修改前后连续
currentPhase = !currentPhase;
RenderTexture rt = currentPhase ? pingPhaseTexture : pongPhaseTexture;
dispersionMat.SetTexture("_Phase", currentPhase? pongPhaseTexture : pingPhaseTexture);
dispersionMat.SetFloat("_DeltaTime", deltaTime * mult);
Graphics.Blit(null, rt, dispersionMat);

spectrumMat.SetTexture("_Phase", currentPhase? pingPhaseTexture : pongPhaseTexture);
Graphics.Blit(null, spectrumTexture, spectrumMat);

fftMat.EnableKeyword("_HORIZONTAL");
fftMat.DisableKeyword("_VERTICAL");
int iterations = Mathf.CeilToInt((float)Math.Log(resolution * 8, 2)) * 2;
for (int i = 0; i < iterations; i++)
{
fftMat.SetFloat("_SubTransformSize", Mathf.Pow(2, (i % (iterations / 2)) + 1));
if (i == 0)
{
fftMat.SetTexture("_Input", spectrumTexture);
Graphics.Blit(null, pingTransformTexture, fftMat);
}
else if (i == iterations - 1)
{
fftMat.SetTexture("_Input", (iterations % 2 == 0) ? pingTransformTexture : pongTransformTexture);
Graphics.Blit(null, displacementTexture, fftMat);
}
else if (i % 2 == 1)
{
fftMat.SetTexture("_Input", pingTransformTexture);
Graphics.Blit(null, pongTransformTexture, fftMat);
}
else
{
fftMat.SetTexture("_Input", pongTransformTexture);
Graphics.Blit(null, pingTransformTexture, fftMat);
}
if (i == iterations / 2 - 1)
{
fftMat.DisableKeyword("_HORIZONTAL");
fftMat.EnableKeyword("_VERTICAL");
}
}

heightMat.SetTexture("_Phase", currentPhase? pingPhaseTexture : pongPhaseTexture);
Graphics.Blit(null, spectrumTexture, heightMat);
fftMat.EnableKeyword("_HORIZONTAL");
fftMat.DisableKeyword("_VERTICAL");
for (int i = 0; i < iterations; i++)
{
fftMat.SetFloat("_SubTransformSize", Mathf.Pow(2, (i % (iterations / 2)) + 1));
if (i == 0)
{
fftMat.SetTexture("_Input", spectrumTexture);
Graphics.Blit(null, pingTransformTexture, fftMat);
}
else if (i == iterations - 1)
{
fftMat.SetTexture("_Input", (iterations % 2 == 0) ? pingTransformTexture : pongTransformTexture);
Graphics.Blit(null, heightTexture, fftMat);
}
else if (i % 2 == 1)
{
fftMat.SetTexture("_Input", pingTransformTexture);
Graphics.Blit(null, pongTransformTexture, fftMat);
}
else
{
fftMat.SetTexture("_Input", pongTransformTexture);
Graphics.Blit(null, pingTransformTexture, fftMat);
}
if (i == iterations / 2 - 1)
{
fftMat.DisableKeyword("_HORIZONTAL");
fftMat.EnableKeyword("_VERTICAL");
}
}

Graphics.Blit(null, normalTexture, normalMat);
Graphics.Blit(null, whiteTexture, whiteMat);
}```

```Shader "Hidden/Mistral Water/Helper/Vertex/Initial Spectrum"
{
Properties
{
_RandomSeed1 ("Random Seed 1", Float) = 1.5122
_RandomSeed2 ("Random Seed 2", Float) = 6.1152
_Amplitude ("Phillips Ampitude", Float) = 1
_Length ("Wave Length", Float) = 256
_Resolution ("Ocean Resolution", int) = 256
_Wind ("Wind Direction (XY)", Vector) = (1, 1, 0, 0)
}

{
Pass
{
Cull Off
ZWrite Off
ZTest Off

CGPROGRAM

#include "FFTCommon.cginc"

#pragma fragment frag

uniform float _RandomSeed1;
uniform float _RandomSeed2;
uniform float _Amplitude;
uniform float _Length;
uniform int _Resolution;
uniform float2 _Wind;

float4 frag(FFTVertexOutput i) : SV_TARGET
{
float n = (i.texcoord.x * _Resolution);
float m = (i.texcoord.y * _Resolution);

float phi1 = Phillips(n, m, _Amplitude, _Wind, _Resolution, _Length);
float phi2 = Phillips(_Resolution - n, _Resolution - m, _Amplitude, _Wind, _Resolution, _Length);

float2 h0 = hTilde0(i.texcoord.xy, _RandomSeed1 / 2, _RandomSeed2 * 2, phi1);
float2 h0conj = Conj(hTilde0(i.texcoord.xy, _RandomSeed1, _RandomSeed2, phi2));

return float4(h0, h0conj);
}

ENDCG
}
}
}

/// FFTCommon.cginc
//inline float Phillips(float n, float m, float amp, float2 wind, float res, float len)
//{
//	float2 k = GetWave(n, m, len, res);
//	float klen = length(k);
//	float klen2 = klen * klen;
//	float klen4 = klen2 * klen2;
//	if(klen < EPSILON)
//		return 0;
//	float kDotW = dot(normalize(k), normalize(wind));
//	float kDotW2 = kDotW * kDotW;
//	float wlen = length(wind);
//	float l = wlen * wlen / G;
//	float l2 = l * l;
//	float damping = 0.01;
//	float L2 = l2 * damping * damping;
//	return amp * exp(-1 / (klen2 * l2)) / klen4 * kDotW2 * exp(-klen2 * L2);
//}

//inline float2 hTilde0(float2 uv, float r1, float r2, float phi)
//{
//	float2 r;
//	float rand1 = UVRandom(uv, 10.612, r1);
//	float rand2 = UVRandom(uv, 11.899, r2);
//	rand1 = clamp(rand1, 0.01, 1);
//	rand2 = clamp(rand2, 0.01, 1);
//	float x = sqrt(-2 * log(rand1));
//	float y = 2 * PI * rand2;
//	r.x = x * cos(y);
//	r.y = x * sin(y);
//	return r * sqrt(phi / 2);
//}```

```Shader "Hidden/Mistral Water/Helper/Vertex/Dispersion"
{
Properties
{
_Length ("Wave Length", Float) = 256
_Resolution ("Ocean Resolution", int) = 256
_DeltaTime ("Delta Time", Float) = 0.016
_Phase ("Last Phase", 2D) = "black" {}
}

{
Pass
{
Cull Off
ZWrite Off
ZTest Off

CGPROGRAM

#include "FFTCommon.cginc"

#pragma fragment frag

uniform float _DeltaTime;
uniform float _Length;
uniform int _Resolution;
uniform sampler2D _Phase;

float4 frag(FFTVertexOutput i) : SV_TARGET
{
float n = (i.texcoord.x * _Resolution);
float m = (i.texcoord.y * _Resolution);

float deltaPhase = CalcDispersion(n, m, _Length, _DeltaTime, _Resolution);
float phase = tex2D(_Phase, i.texcoord).r;

return float4(GetDispersion(phase, deltaPhase), 0, 0, 0);
}

ENDCG
}
}
}
/// FFTCommon.cginc
//inline float GetDispersion(float oldPhase, float newPhase)
//{
//	return fmod(oldPhase + newPhase, 2 * PI);
//}

//inline float CalcDispersion(float n, float m, float len, float dt, float res)
//{
//	float2 wave = GetWave(n, m, len, res);
//	float w = 2 * PI / len;
//	float wlen = length(wave);
//	return sqrt(G * length(wave) * (1 + wlen * wlen / 370 / 370)) * dt;
//	//return floor(sqrt(G * length(wave) / w)) * w * dt;
//}```

```Shader "Hidden/Mistral Water/Helper/Vertex/Spectrum"
{
Properties
{
_Length ("Wave Length", Float) = 256
_Resolution ("Ocean Resolution", int) = 256
_Phase ("Last Phase", 2D) = "black" {}
_Initial ("Intial Spectrum", 2D) = "black" {}
_Choppiness ("Choppiness", Float) = 1
}

{
Pass
{
Cull Off
ZWrite Off
ZTest Off

CGPROGRAM

#include "FFTCommon.cginc"

#pragma fragment frag

uniform sampler2D _Phase;
uniform sampler2D _Initial;
uniform float _Length;
uniform int _Resolution;
uniform float _Choppiness;

float4 frag(FFTVertexOutput i) : SV_TARGET
{
float n = (i.texcoord.x * _Resolution);
float m = (i.texcoord.y * _Resolution);
float2 wave = GetWave(n, m, _Length, _Resolution);
float w = length(wave);
float phase = tex2D(_Phase, i.texcoord).r;
float2 pv = float2(cos(phase), sin(phase));
float2 h0 = tex2D(_Initial, i.texcoord).rg;
float2 h0conj = tex2D(_Initial, i.texcoord).ba;
//h0conj = MultByI(h0conj);
float2 h = MultComplex(h0, pv) + MultComplex(h0conj, Conj(pv));
//return float4(h, h);
w = max(0.0001, w);
float2 hx = -MultByI(h * wave.x / w) * _Choppiness;
float2 hz = -MultByI(h * wave.y / w) * _Choppiness;
return float4(hx, hz);
}

ENDCG
}
}
}```

```Shader "Hidden/Mistral Water/Helper/Map/Ocean Normal"
{
Properties
{
_Resolution ("Resolution", Float) = 512
_Length ("Sea Width", Float) = 512
_Displacement ("Displacement", 2D) = "black" {}
_Bump ("_Bump", 2D) = "bump" {}
}

{
Pass
{
Cull Off
ZWrite Off
ZTest Off

CGPROGRAM

#include "FFTCommon.cginc"

#pragma fragment frag

uniform float _Resolution;
uniform float _Length;
uniform sampler2D _Displacement;
uniform float4 _Displacement_TexelSize;
uniform sampler2D _Bump;

float4 frag(FFTVertexOutput i) : SV_TARGET
{
float texelSize = 1 / _Length;
float2 dDdy = -0.5 * (tex2D(_Displacement, i.texcoord + float4(0, -texelSize, 0, 0)).rb - tex2D(_Displacement, i.texcoord + float4(0, texelSize, 0, 0)).rb) / 8;
float2 dDdx = -0.5 * (tex2D(_Displacement, i.texcoord + float4(-texelSize, 0, 0, 0)).rb - tex2D(_Displacement, i.texcoord + float4(texelSize, 0, 0, 0)).rb) / 8;
float2 noise = 0.3 * tex2D(_Bump, i.texcoord).xz;
float jacobian = (1 + dDdx.x) * (1 + dDdy.y) - dDdx.y * dDdy.x;
float turb = max(0, 1 - jacobian + length(noise));
float xx = 1 + 3 * smoothstep(1.2, 1.8, turb);
xx = min(turb, 1);
xx = smoothstep(0, 1, turb);
return float4(xx, xx, xx, 1);
}

ENDCG
}
}
}```

## References

[1]GPU Gems

[4]Tessendorf, Jerry. “Simulating ocean water.” Simulating nature: realistic and interactive techniques. SIGGRAPH 1.2 (2001): 5.

[5]Dupuy, Jonathan, and Eric Bruneton. “Real-time animation and rendering of ocean whitecaps.” SIGGRAPH Asia 2012 Technical Briefs . ACM, 2012.

[6] Cooley, James W., and John W. Tukey. “An algorithm for the machine calculation of complex Fourier series.” Mathematics of computation 19.90 (1965): 297-301.

[7]Lloyd, D. Brandon, Chas Boyd, and Naga Govindaraju. “Fast computation of general Fourier transforms on GPUs.” Multimedia and Expo, 2008 IEEE International Conference on . IEEE, 2008.

[9] Pseudorandom number generator. (2017, December 18). In Wikipedia, The Free Encyclopedia . Retrieved 13:38, December 29, 2017, from https://en.wikipedia.org/w/index.php?title=Pseudorandom_number_generator&amp;amp;amp;amp;amp;oldid=816057825

[10]Box-Muller transform

|