GerstnerWave
我们海洋模拟在那个GerstnerWave和FFT是两种常用的方法,这次我参考《GEM GEMS1》实现了一套GerstnerWave,大概简单总结下
顶点位置计算:
(1) x,y也就是坐标点在XZ平面的位置分量
(2) Q为波浪的陡峭参数
(3) D控制波浪的移动方向,为float2
(4) A为波浪的振幅
(5) t 为游戏开始到目前运行的时间
(6) 假设波长为λ, 则 频率 w = 2* π / λ,
(7) 假设我们的波长每秒移动的速度为 S,则相位常量 δ = w * S
总体上计算GerstnerWave的顶点位置的方式是,初始化一个大网格,对网格的每个顶点初始化XZ屏幕的值,代入上面的公式计算。不过从上面公式看
struct WaveParam
{
public:
float fAmplitude;
float fSpeed;
float fWaveLength;
WaveParam::WaveParam(float amplitude, float speed, float waveLength):
fAmplitude(amplitude), fSpeed(speed), fWaveLength(waveLength)
{
}
};
struct GerstnerWaveParam : public DirectionSineWaveParam
{
public:
float fSteepness;
GerstnerWaveParam(XMFLOAT2 direction, float amplitude, float speed,
float waveLength, float steepness) :
DirectionSineWaveParam(direction, amplitude, speed, waveLength),
fSteepness(steepness)
{
}
};
//计算位置
void GerstnerWave::CalculateVertexPos(float time)
{
UINT index = 0;
//从左到右,从下到上
for (int posZ = -mWaveHeight / 2; posZ <= mWaveHeight / 2; ++posZ)
{
for (int posX = -mWaveWidth / 2; posX <= mWaveWidth / 2; ++posX)
{
float fPositionX = (float)posX * mWaveGridSize;
float fPositionZ = (float)posZ * mWaveGridSize;
float fGerstnerX = fPositionX + GetPosXOffset(fPositionX, fPositionZ, time);
float fGerstnerZ = fPositionZ + GetPosZOffset(fPositionX, fPositionZ, time);
//pos
mWaveVertexData[index].pos.x = fGerstnerX;
mWaveVertexData[index].pos.z = fGerstnerZ;
mWaveVertexData[index].pos.y = (float)GetWaveVertexHeight(fPositionX, fPositionZ, time);
++index;
}
}
}
float GerstnerWave::GetWaveVertexHeight(int x, int z, float time)
{
float fHeight = 0.0f;
for (UINT iParamIndex = 0; iParamIndex < m_vecGerstnerWaveParam.size(); ++iParamIndex)
{
DirectionSineWaveParam& dirWaveParam = m_vecGerstnerWaveParam[iParamIndex];
float fDdotXZ = dirWaveParam.fDrection.x * x + dirWaveParam.fDrection.y * z;
//频率
float w = 2.0f * XM_PI / dirWaveParam.fWaveLength;
//相位差常量
float phaseConstant = dirWaveParam.fSpeed * w;
//posZ
fHeight += (float)dirWaveParam.fAmplitude * sin(fDdotXZ * w + time * phaseConstant);
}
return fHeight;
}
float GerstnerWave::GetPosXOffset(float x, float z, float time)
{
float fPosXOffset = 0.0f;
for (UINT iParamIndex = 0; iParamIndex < m_vecGerstnerWaveParam.size(); ++iParamIndex)
{
GerstnerWaveParam& gerstnerWaveParam = m_vecGerstnerWaveParam[iParamIndex];
float fDdotXZ = gerstnerWaveParam.fDrection.x * x + gerstnerWaveParam.fDrection.y * z;
//频率
float w = 2.0f * XM_PI / gerstnerWaveParam.fWaveLength;
//相位差常量
float phaseConstant = gerstnerWaveParam.fSpeed * w;
fPosXOffset += (float)gerstnerWaveParam.fAmplitude * gerstnerWaveParam.fSteepness *
gerstnerWaveParam.fDrection.x * cos(fDdotXZ * w + time * phaseConstant);
}
return fPosXOffset;
}
float GerstnerWave::GetPosZOffset(float x, float z, float time)
{
float fPosZOffset = 0.0f;
for (UINT iParamIndex = 0; iParamIndex < m_vecGerstnerWaveParam.size(); ++iParamIndex)
{
GerstnerWaveParam& gerstnerWaveParam = m_vecGerstnerWaveParam[iParamIndex];
float fDdotXZ = gerstnerWaveParam.fDrection.x * x + gerstnerWaveParam.fDrection.y * z;
//频率
float w = 2.0f * XM_PI / gerstnerWaveParam.fWaveLength;
//相位差常量
float phaseConstant = gerstnerWaveParam.fSpeed * w;
fPosZOffset += (float)gerstnerWaveParam.fAmplitude * gerstnerWaveParam.fSteepness *
gerstnerWaveParam.fDrection.y * cos(fDdotXZ * w + time * phaseConstant);
}
return fPosZOffset;
}
注意一下这里采用的是多组GerstnerWave参数叠加,因为仅仅1组效果不是很好,我参考了下网上的资料以及自己的实验结果,大概三组以及以上效果才能不错。
3组参数:
vector<GerstnerWaveParam> vecGerstnerWaveParam;
vecGerstnerWaveParam.push_back(GerstnerWaveParam(XMFLOAT2(0.7f, 0.7), 2.0f, 4.0f, 20.0f,0.3f));
vecGerstnerWaveParam.push_back(GerstnerWaveParam(XMFLOAT2(0.7f, -0.7), 1.5f, 5.0f, 30.0f, 0.4f));
vecGerstnerWaveParam.push_back(GerstnerWaveParam(XMFLOAT2(-0.7f, -0.3), 1.0f, 6.0f, 40.0f, 0.5f));
网格数是80 * 80, 每个格子大小为 1.0
大致得到的网格运算效果如下:
帧数只有17,有点 惨
顶点法线计算:
每个顶点法线的计算是计算临接的四个三角形的法线的规格向量的平均值。
void Wave::CalculateVertexNormal()
{
for (int row = 0; row < mWaveHeight + 1; ++row)
{
for (int column = 0; column < mWaveWidth + 1; ++column)
{
int index = row * (mWaveWidth + 1) + column;
XMFLOAT3 vertex1 = mWaveVertexData[index].pos;
//左下角
if ((0 == row && 0 == column))
{
XMFLOAT3 vertex2 = mWaveVertexData[(row + 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 vertex3 = mWaveVertexData[row * (mWaveWidth + 1) + column + 1].pos;
mWaveVertexData[index].normal = CrossNormal(vertex1, vertex2, vertex3);
}
//右下角
else if (0 == row && column == mWaveWidth)
{
XMFLOAT3 vertex2 = mWaveVertexData[row * (mWaveWidth + 1) + column - 1].pos;
XMFLOAT3 vertex3 = mWaveVertexData[(row + 1) * (mWaveWidth + 1) + column].pos;
mWaveVertexData[index].normal = CrossNormal(vertex1, vertex2, vertex3);
}
//左上角
else if (mWaveHeight == row && 0 == column)
{
XMFLOAT3 vertex2 = mWaveVertexData[row * (mWaveWidth + 1) + column + 1].pos;
XMFLOAT3 vertex3 = mWaveVertexData[(row - 1) * (mWaveWidth + 1) + column].pos;
mWaveVertexData[index].normal = CrossNormal(vertex1, vertex2, vertex3);
}
//右上角
else if (mWaveHeight == row && column == mWaveWidth)
{
XMFLOAT3 vertex2 = mWaveVertexData[(row - 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 vertex3 = mWaveVertexData[row * (mWaveWidth + 1) + column - 1].pos;
mWaveVertexData[index].normal = CrossNormal(vertex1, vertex2, vertex3);
}
//左边
else if (0 == column)
{
XMFLOAT3 vertex2 = mWaveVertexData[(row + 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 vertex3 = mWaveVertexData[row * (mWaveWidth + 1) + column + 1].pos;
XMFLOAT3 vertex4 = mWaveVertexData[(row - 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 normal1 = CrossNormal(vertex1, vertex2, vertex3);
XMFLOAT3 normal2 = CrossNormal(vertex1, vertex3, vertex4);
mWaveVertexData[index].normal =
XMFLOAT3((normal1.x + normal2.x) / 2.0f, (normal1.y + normal2.y) / 2.0f, (normal1.z + normal2.z) / 2.0f);
}
//右边
else if (column == mWaveWidth)
{
XMFLOAT3 vertex2 = mWaveVertexData[(row - 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 vertex3 = mWaveVertexData[row * (mWaveWidth + 1) + column - 1].pos;
XMFLOAT3 vertex4 = mWaveVertexData[(row + 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 normal1 = CrossNormal(vertex1, vertex2, vertex3);
XMFLOAT3 normal2 = CrossNormal(vertex1, vertex3, vertex4);
mWaveVertexData[index].normal =
XMFLOAT3((normal1.x + normal2.x) / 2.0f, (normal1.y + normal2.y) / 2.0f, (normal1.z + normal2.z) / 2.0f);
}
//下边
else if (0 == row)
{
XMFLOAT3 vertex2 = mWaveVertexData[row * (mWaveWidth + 1) + column - 1].pos;
XMFLOAT3 vertex3 = mWaveVertexData[(row + 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 vertex4 = mWaveVertexData[row * (mWaveWidth + 1) + column + 1].pos;
XMFLOAT3 normal1 = CrossNormal(vertex1, vertex2, vertex3);
XMFLOAT3 normal2 = CrossNormal(vertex1, vertex3, vertex4);
mWaveVertexData[index].normal =
XMFLOAT3((normal1.x + normal2.x) / 2.0f, (normal1.y + normal2.y) / 2.0f, (normal1.z + normal2.z) / 2.0f);
}
//上边
else if (mWaveHeight == row)
{
XMFLOAT3 vertex2 = mWaveVertexData[row * (mWaveWidth + 1) + column + 1].pos;
XMFLOAT3 vertex3 = mWaveVertexData[(row - 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 vertex4 = mWaveVertexData[row * (mWaveWidth + 1) + column - 1].pos;
XMFLOAT3 normal1 = CrossNormal(vertex1, vertex2, vertex3);
XMFLOAT3 normal2 = CrossNormal(vertex1, vertex3, vertex4);
mWaveVertexData[index].normal =
XMFLOAT3((normal1.x + normal2.x) / 2.0f,
(normal1.y + normal2.y) / 2.0f, (normal1.z + normal2.z) / 2.0f);
}
//里面的
else
{
XMFLOAT3 vertex2 = mWaveVertexData[(row + 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 vertex3 = mWaveVertexData[row * (mWaveWidth + 1) + column + 1].pos;
XMFLOAT3 vertex4 = mWaveVertexData[(row - 1) * (mWaveWidth + 1) + column].pos;
XMFLOAT3 vertex5 = mWaveVertexData[row * (mWaveWidth + 1) + column - 1].pos;
XMFLOAT3 normal1 = CrossNormal(vertex1, vertex2, vertex3);
XMFLOAT3 normal2 = CrossNormal(vertex1, vertex3, vertex4);
XMFLOAT3 normal3 = CrossNormal(vertex1, vertex4, vertex5);
XMFLOAT3 normal4 = CrossNormal(vertex1, vertex5, vertex2);
float normalX = (normal1.x + normal2.x + normal3.x + normal4.x) / 4.0f;
float normalY = (normal1.y + normal2.y + normal3.y + normal4.y) / 4.0f;
float normalZ = (normal1.z + normal2.z + normal3.z + normal4.z) / 4.0f;
mWaveVertexData[index].normal = XMFLOAT3(normalX, normalY, normalZ);
}
NormalizeFloat3(mWaveVertexData[index].normal);
}
}
}
因为有些边缘或者角落的点并没有四个临接面,所以要做出一定的条件判断。
顶点切线的计算:
方法和之前的教程 Directx11教程十二之NormalMap(法线贴图)是一样的,利用位置边向量和纹理UV边向量
//计算切线
void Wave::CalculateVertexTangent()
{
for (int row = 0; row < mWaveHeight + 1; ++row)
{
for (int column = 0; column < mWaveWidth + 1; ++column)
{
int index = row * (mWaveWidth + 1) + column;
VertexPCNTT vertex1;
VertexPCNTT vertex2;
VertexPCNTT vertex3;
//除开最上面行和最右行
if (row >= 0 && row <= mWaveHeight - 1 && column >= 0 && column <= mWaveWidth - 1)
{
vertex1 = mWaveVertexData[index];
vertex2 = mWaveVertexData[(row + 1) * (mWaveWidth + 1) + column];
vertex3 = mWaveVertexData[row * (mWaveWidth + 1) + column + 1];
}
//除开右上角的最右行
else if (row != mWaveHeight &&column == mWaveWidth)
{
vertex1 = mWaveVertexData[index];
vertex2 = mWaveVertexData[(row + 1) * (mWaveWidth + 1) + column];
vertex3 = mWaveVertexData[row * (mWaveWidth + 1) + column - 1];
}
//最上面一行
else
{
mWaveVertexData[index].tangent = mWaveVertexData[index - mWaveWidth - 1].tangent;
continue;
}
float Edge1[3], Edge2[3];
float TexEdge1[2], TexEdge2[2];
//计算面的两个向量
//边向量1
Edge1[0] = vertex2.pos.x - vertex1.pos.x; //E0X
Edge1[1] = vertex2.pos.y - vertex1.pos.y; //E0Y
Edge1[2] = vertex2.pos.z - vertex1.pos.z; //E0Z
//边向量2
Edge2[0] = vertex3.pos.x - vertex1.pos.x; //E1X
Edge2[1] = vertex3.pos.y - vertex1.pos.y; //E1Y
Edge2[2] = vertex3.pos.z - vertex1.pos.z; //E1Z
//纹理边向量1
TexEdge1[0] = vertex2.uv.x - vertex1.uv.x; //U0
TexEdge1[1] = vertex2.uv.y - vertex1.uv.y; //V0
//纹理边向量2
TexEdge2[0] = vertex3.uv.x - vertex1.uv.x; //U1
TexEdge2[1] = vertex3.uv.y - vertex1.uv.y; //V1
//求出TB在模型空间坐标的方程系数
float den = 1.0f / (TexEdge1[0] * TexEdge2[1] - TexEdge1[1] * TexEdge2[0]);
//求出Tangent
XMFLOAT3 tangent;
tangent.x = den*(TexEdge2[1] * Edge1[0] - TexEdge1[1] * Edge2[0]);
tangent.y = den*(TexEdge2[1] * Edge1[1] - TexEdge1[1] * Edge2[1]);
tangent.z = den*(TexEdge2[1] * Edge1[2] - TexEdge1[1] * Edge2[2]);
NormalizeFloat3(tangent);
//求出垂直于Nomral的Tangent
XMFLOAT3& normal = mWaveVertexData[index].normal;
float dot = Dot(normal, tangent);
tangent.x =-(tangent.x - dot * normal.x);
tangent.y = -(tangent.y - dot * normal.y);
tangent.z = -(tangent.z - dot * normal.z);
NormalizeFloat3(tangent);
mWaveVertexData[index].tangent = tangent;
}
}
}
顶点UV
这个倒没什么好说的
void Wave::CalculateVertexUV()
{
for (int row = 0; row < mWaveHeight + 1; ++row)
{
for (int column = 0; column < mWaveWidth + 1; ++column)
{
int index = row * (mWaveWidth + 1) + column;
float u = ((float)(column) / (float)mWaveHeight);
float v = ((float)(row) / (float)mWaveWidth);
mWaveVertexData[index].uv.x = u;
mWaveVertexData[index].uv.y = v;
}
}
}
Shader代码:
Texture2D DiffuseTex:register(t0);
Texture2D NormalTex:register(t1);
SamplerState SampleWrapLinear:register(s0);
SamplerState SampleClampPoint:register(s1);
cbuffer CBMatrix:register(b0)
{
matrix World;
matrix View;
matrix Proj;
matrix WorldInvTranspose;
float3 cameraPos;
float pad1;
float4 dirLightColor;
float3 dirLightDir;
float pad2;
float3 ambientLight;
float pad3;
};
cbuffer CBEveryFrame:register(b1)
{
float4 surfaceColor;
}
struct VertexIn
{
float3 Pos:POSITION;
float3 Color:COLOR;
float3 Normal:NORMAL;
float3 Tangent:TANGENT;
float2 Tex:TEXCOORD;
};
struct VertexOut
{
float4 Pos:SV_POSITION;
float3 worldPos:POSITION0;
float3 WorldNormal:NORMAL;
float3 WorldTangent:TANGENT;
float2 Tex:TEXCOORD0;
};
VertexOut VS(VertexIn ina)
{
VertexOut outa;
outa.Pos = mul(float4(ina.Pos, 1.0f), World);
outa.worldPos = outa.Pos.xyz;
outa.Pos = mul(outa.Pos, View);
outa.Pos = mul(outa.Pos, Proj);
outa.WorldNormal = mul(ina.Normal, (float3x3)WorldInvTranspose); //此事世界逆转置矩阵的第四行本来就没啥用
outa.WorldNormal = normalize(outa.WorldNormal);
outa.WorldTangent = mul(ina.Tangent, (float3x3)World);
outa.WorldTangent = normalize(outa.WorldTangent);
outa.Tex = ina.Tex;
return outa;
}
float4 PS(VertexOut outa) : SV_Target
{
float4 color = float4(0.0, 0.0, 0.0, 1.0);
float3 normal = NormalTex.Sample(SampleWrapLinear, outa.Tex).xyz;
normal = normal * 2.0 - 1.0;
float3 N = normalize(outa.WorldNormal);
float3 T = normalize(outa.WorldTangent);
float3 B = normalize(cross(N, T));
float3 worldNormal;
worldNormal.x = normal.x * T.x + normal.y * B.x + normal.z * N.x;
worldNormal.y = normal.x * T.y + normal.y * B.y + normal.z * N.y;
worldNormal.z = normal.x * T.z + normal.y * B.z + normal.z * N.z;
worldNormal = normalize(worldNormal);
//float3 worldNormal = normalize(outa.WorldNormal);
//计算diffuseLight
float3 invLightDir = -normalize(dirLightDir);
float diffuseFactor = saturate(dot(worldNormal, invLightDir));
float3 diffuse = DiffuseTex.Sample(SampleWrapLinear, outa.Tex).xyz;
diffuse = pow(diffuse, float4(2.0f, 2.0f, 2.0f, 0.0f));
//SpecularLight
float3 viewDir = normalize(cameraPos - outa.worldPos);
float3 halfDir = normalize(invLightDir + viewDir);
float spec = pow(saturate(dot(halfDir, worldNormal)), 32);
float3 specular = float3(spec, spec, spec);
color = float4((ambientLight + dirLightColor.xyz * diffuseFactor) * diffuse + specular * 0.0f, 1.0);
float correctGamma = 1.0 / 2.0;
color = pow(color, float4(correctGamma, correctGamma, correctGamma, 0.0));
return color;
}
DiffuseMap:
NormalMap:
80 X 80的大小只有16帧 只能说这样的大量并行数据下 CPU计算性能太差了,所以还是得指望ComputeShader。
下篇博客:Directx11进阶之基于GerstnerWave模拟二之基于ComputeShader的实现
源码:
https://github.com/2047241149/SDEngine 的Wave.h, Wave.Cpp, GerstnerWave.h, GerstnerWave.cpp
参考资料:
【1】《GPU GEMS 1》的 Chapter 1. Effective Water Simulation from Physical Models