DirectX11 Windows Windows SDK--28 计算着色器：波浪(水波)

时间 2020-01-22

标签 directx11 directx windows sdk 计算着色波浪水波栏目 Windows 繁體版

原文原文链接

前言

有关计算着色器的基础其实并非不少。接下来继续讲解如何使用计算着色器实现水波效果，即龙书中所实现的水波。可是光看代码但是彻底看不出来是在作什么的。我的根据书中所给的参考书籍找到了对应的实现原理，可是里面涉及到比较多的物理公式。所以，要看懂这一章须要有高数功底（求导、偏导、微分方程），我会把推导过程给列出来。html

本章演示项目还用到了其余的一些效果，在学习本章以前建议先了解以下内容：git

章节内容
11 混合状态
17 利用几何着色器实现公告板效果（雾效部分）
26 计算着色器：入门
27 计算着色器：双调排序（非双调排序部分）

学习目标：github

熟悉如何利用三角形栅格来表示地形和水
熟悉如何用CPU或者计算着色器生成水波

DirectX11 With Windows SDK完整目录算法

Github项目源码数组

欢迎加入QQ群: 727623616 能够一块儿探讨DX11，以及有什么问题也能够在这里汇报。缓存

用三角形栅格表示地形

咱们可使用一个函数\(y=f(x, z)\)来表示一个曲面，在xz平面内构造一个栅格来近似地表示这个曲面。其中的每一个四边形都是由两个三角形所构成的，接下来再利用该函数计算出每一个栅格点处的高度便可。app

以下图所示，咱们先在xz平面内“铺设”一层栅格ide

而后咱们再运用函数\(y=f(x, z)\)来为每一个栅格点获取对应的y坐标。再利用\((x, f(x, z), z)\)的全部顶点构造出地形栅格函数

生成栅格顶点

由以上分析可知，咱们首先须要完成的就是来构建xz平面内的栅格。若规定一个m行n列的栅格，那么栅格包含m * n个四边形，即对应2 * m * n个三角形，顶点数为(m + 1) * (n + 1)。学习

在Geometry中建立地形的方法比较复杂，后三行的形参是根据(x, z)坐标来分别肯定出高度y、法向量n和颜色c的函数：

template<class VertexType = VertexPosNormalTex, class IndexType = DWORD>
MeshData<VertexType, IndexType> CreateTerrain(
    const DirectX::XMFLOAT2& terrainSize,           // 地形宽度与深度
    const DirectX::XMUINT2& slices,                 // 行栅格数与列栅格数
    const DirectX::XMFLOAT2 & maxTexCoord,          // 最大纹理坐标(texU, texV)
    const std::function<float(float, float)>& heightFunc,               // 高度函数y(x, z)
    const std::function<DirectX::XMFLOAT3(float, float)>& normalFunc,   // 法向量函数n(x, z)
    const std::function<DirectX::XMFLOAT4(float, float)>& colorFunc);   // 颜色函数c(x, z)

template<class VertexType = VertexPosNormalTex, class IndexType = DWORD>
    MeshData<VertexType, IndexType> CreateTerrain(
    float width, float depth,               // 地形宽度与深度
    UINT slicesX, UINT slicesZ,             // 行栅格数与列栅格数
    float texU, float texV,                 // 最大纹理坐标(texU, texV)
    const std::function<float(float, float)>& heightFunc,   // 高度函数y(x, z)
    const std::function<DirectX::XMFLOAT3(float, float)>& normalFunc,   // 法向量函数n(x, z)
    const std::function<DirectX::XMFLOAT4(float, float)>& colorFunc);   // 颜色函数c(x, z)

咱们的代码是从左下角开始，逐渐向右向上地计算出其他顶点坐标。须要注意的是，纹理坐标是以左上角为坐标原点，U轴朝右，V轴朝下。

下面的代码展现了如何生成栅格顶点：

template<class VertexType, class IndexType>
    MeshData<VertexType, IndexType> CreateTerrain(
    float width, float depth, 
    UINT slicesX, UINT slicesZ, 
    float texU, float texV, 
    const std::function<float(float, float)>& heightFunc, 
    const std::function<DirectX::XMFLOAT3(float, float)>& normalFunc, 
    const std::function<DirectX::XMFLOAT4(float, float)>& colorFunc)
{
    using namespace DirectX;

    MeshData<VertexType, IndexType> meshData;
    UINT vertexCount = (slicesX + 1) * (slicesZ + 1);
    UINT indexCount = 6 * slicesX * slicesZ;
    meshData.vertexVec.resize(vertexCount);
    meshData.indexVec.resize(indexCount);

    Internal::VertexData vertexData;
    UINT vIndex = 0;
    UINT iIndex = 0;

    float sliceWidth = width / slicesX;
    float sliceDepth = width / slicesZ;
    float leftBottomX = -width / 2;
    float leftBottomZ = -depth / 2;
    float posX, posZ;
    float sliceTexWidth = texU / slicesX;
    float sliceTexDepth = texV / slicesZ;

    XMFLOAT3 normal;
    XMFLOAT4 tangent;
    // 建立网格顶点
    //  __ __
    // | /| /|
    // |/_|/_|
    // | /| /| 
    // |/_|/_|
    for (UINT z = 0; z <= slicesZ; ++z)
    {
        posZ = leftBottomZ + z * sliceDepth;
        for (UINT x = 0; x <= slicesX; ++x)
        {
            posX = leftBottomX + x * sliceWidth;
            // 计算法向量并归一化
            normal = normalFunc(posX, posZ);
            XMStoreFloat3(&normal, XMVector3Normalize(XMLoadFloat3(&normal)));
            // 计算法平面与z=posZ平面构成的直线单位切向量，维持w份量为1.0f
            XMStoreFloat4(&tangent, XMVector3Normalize(XMVectorSet(normal.y, -normal.x, 0.0f, 0.0f)) + g_XMIdentityR3);

            vertexData = { XMFLOAT3(posX, heightFunc(posX, posZ), posZ),
                          normal, tangent, colorFunc(posX, posZ), XMFLOAT2(x * sliceTexWidth, texV - z * sliceTexDepth) };
            Internal::InsertVertexElement(meshData.vertexVec[vIndex++], vertexData);
        }
    }
    // 放入索引
    for (UINT i = 0; i < slicesZ; ++i)
    {
        for (UINT j = 0; j < slicesX; ++j)
        {
            meshData.indexVec[iIndex++] = i * (slicesX + 1) + j;
            meshData.indexVec[iIndex++] = (i + 1) * (slicesX + 1) + j;
            meshData.indexVec[iIndex++] = (i + 1) * (slicesX + 1) + j + 1;

            meshData.indexVec[iIndex++] = (i + 1) * (slicesX + 1) + j + 1;
            meshData.indexVec[iIndex++] = i * (slicesX + 1) + j + 1;
            meshData.indexVec[iIndex++] = i * (slicesX + 1) + j;
        }
    }

    return meshData;
}

其中须要额外了解的是切线向量的产生。因为要让切线向量可以与xOy平面平行，须要先让法向量投影到xOy平面，而后再获得对应的未经标准化的切线向量，以下图所示：

一种表示山峰的函数

正弦函数适合用于表示起伏不定的山坡，一种二维山川地形的函数为：

\[ y(x,z)=\frac{3}{10}(zsin(\frac{1}{10}x)+xcos(\frac{1}{10}z)) \]

其在(x, y, z)处未经标准化的法向量为：

\[ \mathbf{n}=(-\frac{\partial{y}}{\partial{x}}, 1, -\frac{\partial{y}}{\partial{z}}) \]

其中y对x和对z的偏导分别为：
\[ \frac{\partial{y}}{\partial{x}}=\frac{3}{10}(\frac{1}{10}zcos(\frac{1}{10}x)+cos(\frac{1}{10}z))\\ \frac{\partial{y}}{\partial{z}}=\frac{3}{10}(sin(\frac{1}{10}x)-\frac{1}{10}xsin(\frac{1}{10}z)) \]

所以建立上述地形的函数能够写成：

MeshData<VertexPosNormalTex, DWORD> landMesh = Geometry::CreateTerrain(XMFLOAT2(160.0f, 160.0f),
    XMUINT2(50, 50), XMFLOAT2(10.0f, 10.0f),
    [](float x, float z) { return 0.3f*(z * sinf(0.1f * x) + x * cosf(0.1f * z)); },    // 高度函数
    [](float x, float z) { return XMFLOAT3{ -0.03f * z * cosf(0.1f * x) - 0.3f * cosf(0.1f * z), 1.0f, -0.3f * sinf(0.1f * x) + 0.03f * x * sinf(0.1f * z) }; })    // 法向量函数

固然，Lambda函数不理解，你也能够写成这样：

float GetHillsHeight(float x, float z)
{
    return 0.3f * (z * sinf(0.1f * x) + x * cosf(0.1f * z));
}

XMFLOAT3 GetHillsNormal(float x, float z)
{
    return XMFLOAT3(-0.03f * z * cosf(0.1f * x) - 0.3f * cosf(0.1f * z), 1.0f,
        -0.3f * sinf(0.1f * x) + 0.03f * x * sinf(0.1f * z));
}

MeshData<VertexPosNormalTex, DWORD> landMesh = Geometry::CreateTerrain(XMFLOAT2(160.0f, 160.0f), 
    XMUINT2(50, 50), XMFLOAT2(10.0f, 10.0f), GetHillsHeight, GetHillsNormal);

流体模拟

在许多游戏中，你可能会看到有水面流动的场景，实际上他们不必定是真正的流体，而有可能只是一个正在运动的水体表面。出于运行效率的考虑，这些效果的实现或多或少涉及到数学公式或者物理公式。

本节咱们只讨论龙书所实现的方法，即便用波动方程来表现局部位置激起的水波。它的公式推导比较复杂，可是用心看下去的会应该仍是能看得懂的。

话很少说，接下来咱们须要啃硬骨头了。

波动方程的推导

波动方程是一个描述在持续张力做用下的一维绳索或二维表面的某一点处的运动。在一维状况下，咱们能够考虑将一个富有弹性的绳索牢牢绑在两端（有轻微拉伸），让绳索落在x轴上来派生出一维的波动方程。咱们假定当前的绳索在每一点处的线密度(单位长度下的质量)都是恒等的ρ，而且沿着绳索在该点处受到沿着切线的方向持续的张力T。

令函数\(z(x,t)\)表明绳索在x点处，t时刻的垂直位移量。当绳索在z方向产生位移时，代表绳子正在被拉伸。由牛顿第二定律，咱们能够知道在t时刻内，位于\(x=s\)和\(x=s+\Delta x\)之间的绳索段在协力\(\mathbf{F}(x, t)\)的做用下，加速度为：

\[ \mathbf{a}(x,t)=\frac{\mathbf{F}(x,t)}{\rho\Delta x} \tag{28.1} \]
在下面的公式中，咱们能够将位于\(x=s\)和\(x=s+\Delta x\)之间的绳索段的两个端点所受的力分解成水平方向和竖直方向的力，获得\(H(x,t)\)和\(V(x,t)\)。让θ表示绳索在\(x=s\)处切向量与x轴方向的夹角。因为张力T沿着切线方向做用，水平份量\(H(s,t)\)和垂直份量\(V(s,t)\)能够表示为：
\[ H(s,t)=Tcos\theta \\ V(s,t)=Tsin\theta \tag{28.2} \]

让\(\theta+\Delta\theta\)表示另外一个端点\(x=s+\Delta x\)处切线向量与x轴的夹角。这样做用在该端点上的张力的水平份量\(H(s+\Delta x,t)\)和垂直份量\(V(s+\Delta x,t)\)则表示为：
\[ H(s+\Delta x,t)=Tcos(\theta + \Delta\theta) \\ V(s+\Delta x,t)=Tsin(\theta + \Delta\theta) \tag{28.3} \]
对于小幅度的运动，咱们假定绳索段的水平协力为0，这样绳索段的加速度仅仅包含垂直方向。所以，对于\(x=s\)与\(x=s+\Delta x\)之间的绳索段，咱们有下面的公式：
\[ H(s+\Delta x, t) - H(s, t) = 0 \tag{28.4} \]
这样H函数就不须要依赖x了，咱们能够用\(H(t)\)取代\(H(x,t)\)。

做用在\(x=s\)和\(x=s+\Delta x\)之间的绳索段的垂直协力会在z方向产生一个加速度。因为垂直加速度等于位置函数\(z(x,t)的二阶导\)，所以有：
\[ a_z(s,t)=\frac{{\partial}^2}{{\partial}t^2}z(s,t)=\frac{V(s+\Delta x,t)-V(s,t)}{\rho\Delta x} \tag{28.5} \]
等式两边乘上ρ，而后让\(\Delta x\)趋向于0，此时等式右边正好为偏导数的定义：
\[ \rho\frac{{\partial}^2}{{\partial}t^2}z(s,t)=\lim_{\Delta x \to 0}\frac{V(s+\Delta x,t)-V(s,t)}{\Delta x} \tag{28.6} \]
所以又有：
\[ \rho\frac{{\partial}^2}{{\partial}t^2}z(s,t)=\frac{\partial}{{\partial}x}V(s,t) \tag{28.7} \]

将方程组(28.2)联立，咱们能够写成：
\[ V(s, t)=H(t)tan\theta \tag{28.8} \]
由于θ是绳索在X轴方向与切线之间的夹角，\(\tan\theta\)正好也是函数z(x, t)在s处的斜率，所以：
\[ V(s, t)=H(t)\frac{\partial}{\partial x}z(s,t) \tag{28.9} \]
即：
\[ \rho\frac{{\partial}^2}{{\partial}t^2}z(s,t)=\frac{\partial}{\partial x}[H(t)\frac{\partial}{\partial x}z(s,t)] \tag{28.10} \]

因为H(t)并不依赖于x，咱们又能够写成：
\[ \rho\frac{{\partial}^2}{{\partial}t^2}z(s,t)=H(t)\frac{{\partial}^2}{\partial x^2}z(s,t) \tag{28.11} \]
对于小幅度的运动，\(cos\theta\)接近于1（此时水平份量的力为0），所以咱们用\(H(t)\)来近似表示张力T。让\(c^2=T/\rho\)，咱们如今获得了一维的波动方程：
\[ \frac{{\partial}^2 z}{{\partial}t^2}=c^2\frac{{\partial}^2 z}{\partial x^2} \tag{28.12} \]
同理，二维的波动方程能够经过添加一个y项获得：
\[ \frac{{\partial}^2 z}{{\partial}t^2}=c^2(\frac{{\partial}^2 z}{\partial x^2}+\frac{{\partial}^2 z}{\partial y^2}) \tag{28.13} \]
常数c具备单位时间距离的尺度，所以能够表示速度。事实上咱们也不会证实c实际上就是波沿着绳索或表面传递的速度。这是有意义的，由于波的速度随介质所受张力T的变大而增长，随介质密度μ的减少而减少。

知足一维波动方程的解有无穷多个，例如一种常见的波函数形式为\(z(x,t)=Asin(\omega(t-\frac{x}{v}))\)，函数随着时间的推移图像以下：

然而方程(28.13)仅包含了张力，没有考虑到其余阻力因素，这致使波的平均振幅并不会有任何损耗。咱们能够给方程组添加一个与张力方向相反，且与点的运动速度有关的粘性阻尼力：
\[ \frac{{\partial}^2 z}{{\partial}t^2}=c^2(\frac{{\partial}^2 z}{\partial x^2}+\frac{{\partial}^2 z}{\partial y^2})-\mu\frac{\partial z}{\partial t} \tag{28.14} \]
其中非负实数μ表明了液体的粘性，用来控制水面的波何时可以平静下来。μ越大，水波消亡的速度越快。对于水来讲，一般它的μ值会比较小，使得水波可以存留比较长的时间；但对于油来讲，它的μ值会比较大一些，所以水波消亡的速度会比较快。

近似导数

带粘性阻尼力的二维波动方程(28.14)能够经过可分离变量的形式解出来。然而它的解十分复杂，须要大规模的实时模拟演算。取而代之的是，咱们将使用一种数值上的技术来模拟波在流体表面的传播。

假定咱们的流体表面能够表示成一个n×m的矩形栅格，以下图所示。

其中d为两个邻近顶点在x方向的距离及y方向的距离（规定相等），t为时间间隔。咱们用\(z(i, j, k)\)来表示顶点的位移量，其中i和j分别要知足\(0\leq i<n\)及\(0\leq j<m\)，表明世界位置(id, jd)的顶点。k则用来表示时间。所以，\(z(i, j, k)\)等价于顶点(id, jd)在t时刻的z方向位移。

此外，咱们须要施加边界条件，让处在边缘的顶点位移量固定为0.内部顶点的偏移则可使用(28.14)方程，采用近似导数的方法计算出来。以下图所示，咱们能够经过在x方向上计算顶点[i][j]分别和它的相邻的两个顶点[i-1][j]和[i+1][j]的微分的平均值\(\frac{\Delta z}{\Delta x}\)来逼近具备坐标[i][j]的顶点上与x轴对齐的切线。这样就有\(\Delta x = d\)，咱们将偏导数\(\frac{\partial z}{\partial x}\)定义为：

\[ \begin{equation} \begin{split} \frac{\partial}{\partial x}z(i,j,k) &= \frac{\frac{z(i,j,k)-z(i-1,j,k)}{d}+\frac{z(i+1,j,k)-z(i,j,k)}{d}}{2} \\ &=\frac{z(i+1,j,k)-z(i-1,j,k)}{2d} \\ \end{split} \end{equation} \tag{28.15} \]

同理，咱们取顶点[i][j]的两个y方向上的相邻顶点[i][j-1]和[i][j+1]来近似计算z对于y的偏导数：
\[ \frac{\partial}{\partial y}z(i,j,k) =\frac{z(i,j+1,k)-z(i,j-1,k)}{2d} \tag{28.16} \]
至于时间，咱们能够经过计算顶点在当前时刻分别与上一个时刻和下一个时刻的平均位移差来定义z对时间t偏导\(\frac{\Delta z}{\Delta t}\)：
\[ \frac{\partial}{\partial t}z(i,j,k) =\frac{z(i,j,k+1)-z(i,j,k-1)}{2t} \tag{28.17} \]
二阶偏导也可使用和一阶偏导相同的方法来计算。假如咱们已经计算出了顶点[i-1][j]和[i+1][j]的偏移z对x的一阶偏导数，那么咱们就能够获得二者平均差值：
\[ \begin{equation} \begin{split} \Delta[\frac{\partial}{\partial x}z(i,j,k)] &= \frac{\frac{\partial}{\partial x}z(i+1,j,k)-\frac{\partial}{\partial x}z(i,j,k)}{2} + \frac{\frac{\partial}{\partial x}z(i,j,k)-\frac{\partial}{\partial x}z(i-1,j,k)}{2} \\ &=\frac{\frac{\partial}{\partial x}z(i+1,j,k)-\frac{\partial}{\partial x}z(i-1,j,k)}{2} \\ \end{split} \end{equation} \tag{28.18} \]
将方程(28.15)带入上式，能够获得：
\[ \begin{equation} \begin{split} \Delta[\frac{\partial}{\partial x}z(i,j,k)] &= \frac{\frac{z(i+2,j,k)-z(i,j,k)}{2d} - \frac{z(i,j,k)-z(i-2,j,k)}{2d}}{2} \\ &=\frac{z(i+2,j,k)-2z(i,j,k)+z(i-2,j,k)}{4d} \\ \end{split} \end{equation} \tag{28.19} \]
除以d使得咱们能够获得\(\Delta(\frac{\partial z}{\partial x})/\Delta x\)，对应二阶偏导：
\[ \frac{{\partial}^2}{{\partial}x^2}z(i,j,k)=\frac{z(i+2,j,k)-2z(i,j,k)+z(i-2,j,k)}{4d^2} \tag{28.20} \]

该公式要求咱们使用x轴距离为2的顶点[i+2][j]和[i-2][j]来计算二阶偏导。不过相邻的两个顶点就没有用到了，咱们能够基于顶点[i][j]将x轴缩小一半，使用距离最近的两个相邻点[i+1][j]和[i-1][j]来计算二阶偏导：
\[ \frac{{\partial}^2}{{\partial}x^2}z(i,j,k)=\frac{z(i+1,j,k)-2z(i,j,k)+z(i-1,j,k)}{d^2} \tag{28.21} \]
同理可得：
\[ \frac{{\partial}^2}{{\partial}y^2}z(i,j,k)=\frac{z(i,j+1,k)-2z(i,j,k)+z(i,j-1,k)}{d^2} \tag{28.22} \]

\[ \frac{{\partial}^2}{{\partial}t^2}z(i,j,k)=\frac{z(i,j,k+1)-2z(i,j,k)+z(i,j,k-1)}{t^2} \tag{28.23} \]

求出水面位移

联立z对t的一阶偏导（公式28.17）以及二阶偏导（公式28.2一、28.2二、28.23），带粘性阻尼力的二维波动方程能够表示为：
\[ \frac{z(i,j,k+1)-2z(i,j,k)+z(i,j,k-1)}{t^2}=c^{2}\frac{z(i+1,j,k)-2z(i,j,k)+z(i-1,j,k)}{d^2}\\ +c^{2}\frac{z(i,j+1,k)-2z(i,j,k)+z(i,j-1,k)}{d^2}-\mu\frac{z(i,j,k+1)-z(i,j,k-1)}{2t} \tag{28.24} \]

咱们想要可以在传递模拟间隔t来肯定下一次模拟位移量\(z(i,j,k+1)\)，如今咱们已经知道当前位移量\(z(i,j,k)\)和上一次模拟的位移量\(z(i,j,k-1)\)。所以\(z(i,j,k+1)\)的解为：
\[ z(i,j,k+1)=\frac{4-8c^2t^2/d^2}{\mu t+2}z(i,j,k)+\frac{\mu t-2}{\mu t+2}z(i,j,k-1)\\ +\frac{2c^2t^2/d^2}{\mu t+2}[z(i+1,j,k)+z(i-1,j,k)+z(i,j+1,k)+z(i,j-1,k)] \tag{28.25} \]
这条公式正是咱们想要的。其中常量部分能够预先计算出来，只剩下3个带t的因式和4个加法须要给网格中每一个顶点进行运算。

若是波速c过快，或者时间段t太长，上述式子的偏移量有可能趋向于无穷。为了保持结果是有穷的，咱们须要给上式肯定额外的条件，而且要保证在咱们抬高一个顶点并释放后可以确保水波逐渐远离顶点位置。

假定咱们拥有一个n×m顶点数组(其任意\(z(i,j,0)=0\)和\(z(i,j,1)=0\))，如今让某处顶点被抬高使得\(z(i_0, j_0, 0)=h\)和\(z(i_0, j_0, 1)=h\)，h是一个非0位移值。若该处顶点被释放了2t时间，此时式子(28.25)中第三个加法部分的值为0，故有：
\[ \begin{equation} \begin{split} z(i,j,2)&=\frac{4-8c^2t^2/d^2}{\mu t+2}z(i,j,1)+\frac{\mu t-2}{\mu t+2}z(i,j,0)\\ &=\frac{2-8c^2t^2/d^2+\mu t}{\mu t+2}h \end{split} \end{equation} \tag{28.26} \]
为了让顶点向周围平坦的水面移动，它在2t时刻的位移必需要比在t时刻的更小一些。所以就有：
\[ \left| z(i_0, j_0, 2)\right| < \left| z(i_0, j_0, 1) \right| = \left| h \right| \tag{28.27} \]
代入方程(28.26)，有：
\[ \left| \frac{2-8c^2t^2/d^2 + \mu t}{\mu t + 2} \right|\cdot\left| h\right| < \left| h \right| \tag{28.28} \]
所以，
\[ -1 < \frac{2-8c^2t^2/d^2 + \mu t}{\mu t + 2} < 1 \tag{28.29} \]
把速度c解出来，咱们能够获得：
\[ 0<c<\frac{d}{2t}\sqrt{\mu t + 2} \tag{28.30} \]
这告诉咱们对于公式(28.25)，给定两个顶点的距离d以及时间间隔t，波速c必须小于上式的最大值

又或者，给定距离d和波速c，咱们可以计算出最大的时间间隔t。对不等式(28.29)乘上\((-\mu t + 2)\)来简化获得：
\[ 0<\frac{4c^2}{d^2}t^2<\mu t+2 \tag{28.31} \]
因为t>0，中间部分恒大于0，去掉左半部分解一元二次不等式，并舍去t<=0的部分，得：
\[ 0<t<\frac{\mu+\sqrt{\mu ^2+32c^2/d^2}}{8c^2/d^2} \tag{28.32} \]
使用c区间和t区间外的值会致使位移z的结果呈指数级爆炸。

代码实现

如今咱们须要准备两个二维顶点位置数组，其中一个数组表示的是当前模拟全部顶点的位置，而另外一个数组则表示的是上一次模拟全部顶点的位置。当咱们计算新的位移时，将下一次模拟的结果直接写在存放上一次模拟顶点的数组便可（能够观察公式28.25，咱们须要的是当前顶点上一次模拟、当前模拟的数据，以及当前模拟相邻四个顶点的数据，所以其余顶点在计算位移量的时候不会有影响）。

为了执行光照计算，咱们还须要为每一个顶点获取正确的法向量以及可能正确的切线向量。对于顶点坐标(i, j)，未标准化的，与x轴对齐的切向量T和与y轴对齐的切线向量B以下：
\[ \mathbf{T}=(1,0,\frac{\partial}{\partial x}z(i,j,k))\\ \mathbf{B}=(0,1,\frac{\partial}{\partial y}z(i,j,k)) \tag{28.33} \]
用公式28.15和28.16代入28.33，可得：
\[ \mathbf{T}=(1,0,\frac{z(i+1,j,k)-z(i-1,j,k)}{2d})\\ \mathbf{B}=(0,1,\frac{z(i,j+1,k)-z(i,j-1,k)}{2d}) \tag{28.34} \]
通过叉乘后能够获得未经标准化的法向量：
\[ \begin{equation} \begin{split} \mathbf{N} &= \mathbf{T}\times\mathbf{B} \\ &=\begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ 1 & 0 & \frac{z(i+1,j,k)-z(i-1,j,k)}{2d} \\ 0 & 1 & \frac{z(i,j+1,k)-z(i,j-1,k)}{2d} \\ \end{vmatrix} \\ &= (-\frac{z(i+1,j,k)-z(i-1,j,k)}{2d},-\frac{z(i,j+1,k)-z(i,j-1,k)}{2d}, 1) \end{split} \end{equation} \tag{28.35} \]
对上述向量乘上2d的倍数并不会改变它的方向，但能够消除除法：
\[ \mathbf{T}=(2d,0,z(i+1,j,k)-z(i-1,j,k))\\ \mathbf{B}=(0,2d,z(i,j+1,k)-z(i,j-1,k))\\ \mathbf{N}=(z(i-1,j,k)-z(i+1,j,k),z(i,j-1,k)-z(i,j+1,k), 2d) \tag{28.36} \]
注意这里T和B并无相互正交。

要意识到两次模拟期间的时间间隔必须是恒定的，而不是依据帧间隔。不一样游戏运行帧数不一样，所以在帧速率较高的状况下，必须采用必定的机制保证在通过足够多的时间后才更新位置。

CPU计算法

利用上述方式实现的波浪有两种代码实现形式：

使用动态顶点缓冲区，模拟过程在CPU完成，而后结果写入顶点缓冲区
经过计算着色器，将结果写入位移贴图（即纹理数据存放的是位移值，而不是颜色），而后在渲染的时候再利用它

其中前者的效率通常不如后者（在Debug模式下差异过于明显，而在Release模式下好不少），但实现起来比较简单。

这里两种方法都已经实现了，按顺序讲解。

基类WavesRender

基类WavesRender定义及Init过程以下：

class WavesRender
{
public:
    template<class T>
    using ComPtr = Microsoft::WRL::ComPtr<T>;

    void SetMaterial(const Material& material);
    void XM_CALLCONV SetWorldMatrix(DirectX::FXMMATRIX world);

    UINT RowCount() const;
    UINT ColumnCount() const;

protected:
    // 不容许直接构造WavesRender，请从CpuWavesRender或GpuWavesRender构造
    WavesRender() = default;
    ~WavesRender() = default;
    // 不容许拷贝，容许移动
    WavesRender(const WavesRender&) = delete;
    WavesRender& operator=(const WavesRender&) = delete;
    WavesRender(WavesRender&&) = default;
    WavesRender& operator=(WavesRender&&) = default;

    void Init(
        UINT rows,          // 顶点行数
        UINT cols,          // 顶点列数
        float texU,         // 纹理坐标U方向最大值
        float texV,         // 纹理坐标V方向最大值
        float timeStep,     // 时间步长
        float spatialStep,  // 空间步长
        float waveSpeed,    // 波速
        float damping,      // 粘性阻尼力
        float flowSpeedX,   // 水流X方向速度
        float flowSpeedY);  // 水流Y方向速度

protected:
    UINT m_NumRows;                     // 顶点行数
    UINT m_NumCols;                     // 顶点列数

    UINT m_VertexCount;                 // 顶点数目
    UINT m_IndexCount;                  // 索引数目

    DirectX::XMFLOAT4X4 m_WorldMatrix;  // 世界矩阵
    DirectX::XMFLOAT2 m_TexOffset;      // 纹理坐标偏移
    float m_TexU;                       // 纹理坐标U方向最大值
    float m_TexV;                       // 纹理坐标V方向最大值
    Material m_Material;                // 水面材质

    float m_FlowSpeedX;                 // 水流X方向速度
    float m_FlowSpeedY;                 // 水流Y方向速度
    float m_TimeStep;                   // 时间步长
    float m_SpatialStep;                // 空间步长
    float m_AccumulateTime;             // 累积时间

    //
    // 咱们能够预先计算出来的常量
    //

    float m_K1;
    float m_K2;
    float m_K3;

    
};

void WavesRender::Init(UINT rows, UINT cols, float texU, float texV,
    float timeStep, float spatialStep, float waveSpeed, float damping, float flowSpeedX, float flowSpeedY)
{
    XMStoreFloat4x4(&m_WorldMatrix, XMMatrixIdentity());

    m_NumRows = rows;
    m_NumCols = cols;

    m_TexU = texU;
    m_TexV = texV;
    m_TexOffset = XMFLOAT2();

    m_VertexCount = rows * cols;
    m_IndexCount = 6 * (rows - 1) * (cols - 1);

    m_TimeStep = timeStep;
    m_SpatialStep = spatialStep;
    m_FlowSpeedX = flowSpeedX;
    m_FlowSpeedY = flowSpeedY;
    m_AccumulateTime = 0.0f;

    float d = damping * timeStep + 2.0f;
    float e = (waveSpeed * waveSpeed) * (timeStep * timeStep) / (spatialStep * spatialStep);
    m_K1 = (damping * timeStep - 2.0f) / d;
    m_K2 = (4.0f - 8.0f * e) / d;
    m_K3 = (2.0f * e) / d;
}

Init方法将公式(28.25)的三个重要参数先初始化，以供后续计算使用。

CPUWavesRender类

CPUWavesRender的定义以下：

class CpuWavesRender : public WavesRender
{
public:
    CpuWavesRender() = default;
    ~CpuWavesRender() = default;
    // 不容许拷贝，容许移动
    CpuWavesRender(const CpuWavesRender&) = delete;
    CpuWavesRender& operator=(const CpuWavesRender&) = delete;
    CpuWavesRender(CpuWavesRender&&) = default;
    CpuWavesRender& operator=(CpuWavesRender&&) = default;

    HRESULT InitResource(ID3D11Device* device,
        const std::wstring& texFileName,    // 纹理文件名
        UINT rows,          // 顶点行数
        UINT cols,          // 顶点列数
        float texU,         // 纹理坐标U方向最大值
        float texV,         // 纹理坐标V方向最大值
        float timeStep,     // 时间步长
        float spatialStep,  // 空间步长
        float waveSpeed,    // 波速
        float damping,      // 粘性阻尼力
        float flowSpeedX,   // 水流X方向速度
        float flowSpeedY);  // 水流Y方向速度

    void Update(float dt);
    
    // 在顶点[i][j]处激起高度为magnitude的波浪
    // 仅容许在1 < i < rows和1 < j < cols的范围内激起
    void Disturb(UINT i, UINT j, float magnitude);
    // 绘制水面
    void Draw(ID3D11DeviceContext* deviceContext, BasicEffect& effect);

    void SetDebugObjectName(const std::string& name);
private:
    

    std::vector<VertexPosNormalTex> m_Vertices;         // 保存当前模拟结果的顶点二维数组的一维展开
    std::vector<VertexPos> m_PrevSolution;              // 保存上一次模拟结果的顶点位置二维数组的一维展开

    ComPtr<ID3D11Buffer> m_pVertexBuffer;               // 当前模拟的顶点缓冲区
    ComPtr<ID3D11Buffer> m_pIndexBuffer;                // 当前模拟的索引缓冲区

    ComPtr<ID3D11ShaderResourceView> m_pTextureDiffuse; // 水面纹理
    bool m_isUpdated;                                   // 当前是否有顶点数据更新
};

其中顶点的位置咱们只保留了两个副本，即当前模拟和上一次模拟的。而因为在计算出下一次模拟的结果后，咱们就能够抛弃掉上一次模拟的结果。所以咱们能够直接把结果写在存放上一次模拟的位置，而后再进行交换便可。此时本来是当前模拟的数据则变成了上一次模拟的数据，而下一次模拟的结果则变成了当前模拟的数据。顶点的法向量只须要在完成了下一次模拟后再更新，所以也不须要多余的副本了。

在使用CpuWavesRender以前当然是要调用InitResource先进行初始化的，但如今咱们跳过这部分代码，直接看和算法相关的几个方法。

CpuWavesRender::Disturb方法--激起波浪

因为咱们施加了边界0的条件，所以不能对边界区域激起波浪。在修改高度时，咱们还对目标点的相邻四个顶点也修改了高度使得一开始的波浪不会看起来太突兀：

void CpuWavesRender::Disturb(UINT i, UINT j, float magnitude)
{
    // 不要对边界处激起波浪
    assert(i > 1 && i < m_NumRows - 2);
    assert(j > 1 && j < m_NumCols - 2);

    float halfMag = 0.5f * magnitude;

    // 对顶点[i][j]及其相邻顶点修改高度值
    size_t curr = i * (size_t)m_NumCols + j;
    m_Vertices[curr].pos.y += magnitude;
    m_Vertices[curr - 1].pos.y += halfMag;
    m_Vertices[curr + 1].pos.y += halfMag;
    m_Vertices[curr - m_NumCols].pos.y += halfMag;
    m_Vertices[curr + m_NumCols].pos.y += halfMag;

    m_isUpdated = true;
}

CpuWavesRender::Update方法--更新波浪

以前提到，两次模拟期间的时间间隔必须是恒定的，而不是依据帧间隔。所以在设置好初始的时间步长后，每当经历了大于时间步长的累积时间就能够进行更新了。一样在更新过程当中咱们要始终限制边界值为0。虽然公式复杂，但好在实现过程并不复杂。详细见代码：

void CpuWavesRender::Update(float dt)
{
    m_AccumulateTime += dt;
    m_TexOffset.x += m_FlowSpeedX * dt;
    m_TexOffset.y += m_FlowSpeedY * dt;

    // 仅仅在累积时间大于时间步长时才更新
    if (m_AccumulateTime > m_TimeStep)
    {
        m_isUpdated = true;
        // 仅仅对内部顶点进行更新
        for (size_t i = 1; i < m_NumRows - 1; ++i)
        {
            for (size_t j = 1; j < m_NumCols - 1; ++j)
            {
                // 在此次更新以后，咱们将丢弃掉上一次模拟的数据。
                // 所以咱们将运算的结果保存到Prev[i][j]的位置上。
                // 注意咱们可以使用这种原址更新是由于Prev[i][j]
                // 的数据仅在当前计算Next[i][j]的时候才用到
                m_PrevSolution[i * m_NumCols + j].pos.y =
                    m_K1 * m_PrevSolution[i * m_NumCols + j].pos.y +
                    m_K2 * m_Vertices[i * m_NumCols + j].pos.y +
                    m_K3 * (m_Vertices[(i + 1) * m_NumCols + j].pos.y +
                        m_Vertices[(i - 1) * m_NumCols + j].pos.y +
                        m_Vertices[i * m_NumCols + j + 1].pos.y +
                        m_Vertices[i * m_NumCols + j - 1].pos.y);
            }
        }

        // 因为把下一次模拟的结果写到了上一次模拟的缓冲区内，
        // 咱们须要将下一次模拟的结果与当前模拟的结果交换
        for (size_t i = 1; i < m_NumRows - 1; ++i)
        {
            for (size_t j = 1; j < m_NumCols - 1; ++j)
            {
                std::swap(m_PrevSolution[i * m_NumCols + j].pos, m_Vertices[i * m_NumCols + j].pos);
            }
        }

        m_AccumulateTime = 0.0f;    // 重置时间

        // 使用有限差分法计算法向量
        for (size_t i = 1; i < m_NumRows - 1; ++i)
        {
            for (size_t j = 1; j < m_NumCols - 1; ++j)
            {
                float left = m_Vertices[i * m_NumCols + j - 1].pos.y;
                float right = m_Vertices[i * m_NumCols + j + 1].pos.y;
                float top = m_Vertices[(i - 1) * m_NumCols + j].pos.y;
                float bottom = m_Vertices[(i + 1) * m_NumCols + j].pos.y;
                m_Vertices[i * m_NumCols + j].normal = XMFLOAT3(-right + left, 2.0f * m_SpatialStep, bottom - top);
                XMVECTOR nVec = XMVector3Normalize(XMLoadFloat3(&m_Vertices[i * m_NumCols + j].normal));
                XMStoreFloat3(&m_Vertices[i * m_NumCols + j].normal, nVec);
            }
        }   
    }
}

CpuWavesRender::Draw方法--绘制波浪

这里的绘制跟以前用的BasicEffect是能够直接适配的，动态顶点缓冲区的更新只须要在数据发生变化时再进行，以减小CPU向GPU的数据传输次数：

void CpuWavesRender::Draw(ID3D11DeviceContext* deviceContext, BasicEffect& effect)
{
    // 更新动态缓冲区的数据
    if (m_isUpdated)
    {
        m_isUpdated = false;
        D3D11_MAPPED_SUBRESOURCE mappedData;
        deviceContext->Map(m_pVertexBuffer.Get(), 0, D3D11_MAP_WRITE_DISCARD, 0, &mappedData);
        memcpy_s(mappedData.pData, m_VertexCount * sizeof(VertexPosNormalTex),
            m_Vertices.data(), m_VertexCount * sizeof(VertexPosNormalTex));
        deviceContext->Unmap(m_pVertexBuffer.Get(), 0);
    }

    UINT strides[1] = { sizeof(VertexPosNormalTex) };
    UINT offsets[1] = { 0 };
    deviceContext->IASetVertexBuffers(0, 1, m_pVertexBuffer.GetAddressOf(), strides, offsets);
    deviceContext->IASetIndexBuffer(m_pIndexBuffer.Get(), DXGI_FORMAT_R32_UINT, 0);

    effect.SetMaterial(m_Material);
    effect.SetTextureDiffuse(m_pTextureDiffuse.Get());
    effect.SetWorldMatrix(XMLoadFloat4x4(&m_WorldMatrix));
    effect.SetTexTransformMatrix(XMMatrixScaling(m_TexU, m_TexV, 1.0f) * XMMatrixTranslationFromVector(XMLoadFloat2(&m_TexOffset)));
    effect.Apply(deviceContext);
    deviceContext->DrawIndexed(m_IndexCount, 0, 0);
    
}

GPU计算法

相比CPU计算法，GPU计算法的实现则更为复杂了。由于它不只须要用到计算着色器，还须要跑到绘制基本物体的特效那边作一些修改才能用。

HLSL代码

首先计算着色器部分完成的是激起波浪和更新波浪的部分，分为两个函数：

// Waves.hlsli
// 用于更新模拟
cbuffer cbUpdateSettings : register(b0)
{
    float g_WaveConstant0;
    float g_WaveConstant1;
    float g_WaveConstant2;
    float g_DisturbMagnitude;
    
    int2 g_DisturbIndex;
    float2 g_Pad;
}

RWTexture2D<float> g_PrevSolInput : register(u0);
RWTexture2D<float> g_CurrSolInput : register(u1);
RWTexture2D<float> g_Output : register(u2);

// WavesDisturb_CS.hlsl
#include "Waves.hlsli"

[numthreads(1, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
    // 咱们不须要进行边界检验，由于:
    // --读取超出边界的区域结果为0，和咱们对边界处理的需求一致
    // --对超出边界的区域写入并不会执行
    uint x = g_DisturbIndex.x;
    uint y = g_DisturbIndex.y;
    
    float halfMag = 0.5f * g_DisturbMagnitude;
    
    // RW型资源容许读写，因此+=是容许的
    g_Output[uint2(x, y)] += g_DisturbMagnitude;
    g_Output[uint2(x + 1, y)] += halfMag;
    g_Output[uint2(x - 1, y)] += halfMag;
    g_Output[uint2(x, y + 1)] += halfMag;
    g_Output[uint2(x, y - 1)] += halfMag;
}

// WavesUpdate.hlsl
#include "Waves.hlsli"

[numthreads(16, 16, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
    // 咱们不须要进行边界检验，由于:
    // --读取超出边界的区域结果为0，和咱们对边界处理的需求一致
    // --对超出边界的区域写入并不会执行
    uint x = DTid.x;
    uint y = DTid.y;
    
    g_Output[uint2(x, y)] =
        g_WaveConstant0 * g_PrevSolInput[uint2(x, y)].x +
        g_WaveConstant1 * g_CurrSolInput[uint2(x, y)].x +
        g_WaveConstant2 * (
            g_CurrSolInput[uint2(x, y + 1)].x +
            g_CurrSolInput[uint2(x, y - 1)].x +
            g_CurrSolInput[uint2(x + 1, y)].x +
            g_CurrSolInput[uint2(x - 1, y)].x);
    
}

因为所有过程交给了GPU完成，如今咱们须要有三个UAV，两个用于输入，一个用于输出。而且因为咱们指定了线程组内部包含16x16个线程，在C++初始化GpuWavesRender时，咱们也应该指定行顶点数和列顶点数都为16的倍数。

此外，由于GPU对边界外的良好定义，这使得咱们不须要约束调用Disturb的索引条件。

紧接着就是要修改BasicEffect里面用到的顶点着色器以支持计算着色器的位移贴图：

#include "LightHelper.hlsli"

Texture2D g_DiffuseMap : register(t0);          // 物体纹理
Texture2D g_DisplacementMap : register(t1);     // 位移贴图
SamplerState g_SamLinearWrap : register(s0);    // 线性过滤+Wrap采样器
SamplerState g_SamPointClamp : register(s1);    // 点过滤+Clamp采样器

cbuffer CBChangesEveryInstanceDrawing : register(b0)
{
    matrix g_World;
    matrix g_WorldInvTranspose;
    matrix g_TexTransform;
}

cbuffer CBChangesEveryObjectDrawing : register(b1)
{
    Material g_Material;
}

cbuffer CBChangesEveryFrame : register(b2)
{
    matrix g_View;
    float3 g_EyePosW;
    float g_Pad;
}

cbuffer CBDrawingStates : register(b3)
{
    float4 g_FogColor;
    int g_FogEnabled;
    float g_FogStart;
    float g_FogRange;
    int g_TextureUsed;
    
    int g_WavesEnabled;                     // 开启波浪绘制
    float2 g_DisplacementMapTexelSize;      // 位移贴图两个相邻像素对应顶点之间的x,y方向间距
    float g_GridSpatialStep;                // 栅格空间步长
}

cbuffer CBChangesOnResize : register(b4)
{
    matrix g_Proj;
}

cbuffer CBChangesRarely : register(b5)
{
    DirectionalLight g_DirLight[5];
    PointLight g_PointLight[5];
    SpotLight g_SpotLight[5];
}

struct VertexPosNormalTex
{
    float3 PosL : POSITION;
    float3 NormalL : NORMAL;
    float2 Tex : TEXCOORD;
};

struct InstancePosNormalTex
{
    float3 PosL : POSITION;
    float3 NormalL : NORMAL;
    float2 Tex : TEXCOORD;
    matrix World : World;
    matrix WorldInvTranspose : WorldInvTranspose;
};

struct VertexPosHWNormalTex
{
    float4 PosH : SV_POSITION;
    float3 PosW : POSITION; // 在世界中的位置
    float3 NormalW : NORMAL; // 法向量在世界中的方向
    float2 Tex : TEXCOORD;
};

而后是修改BasicObject_VS.hlsl。由于水面是单个物体，不须要改到BasicInstance_VS.hlsl里面：

// BasicObject_VS.hlsl
#include "Basic.hlsli"

// 顶点着色器
VertexPosHWNormalTex VS(VertexPosNormalTex vIn)
{
    VertexPosHWNormalTex vOut;
    
    // 绘制水波时用到
    if (g_WavesEnabled)
    {
        // 使用映射到[0,1]x[0,1]区间的纹理坐标进行采样
        vIn.PosL.y += g_DisplacementMap.SampleLevel(g_SamLinearWrap, vIn.Tex, 0.0f).r;
        // 使用有限差分法估算法向量
        float du = g_DisplacementMapTexelSize.x;
        float dv = g_DisplacementMapTexelSize.y;
        float left = g_DisplacementMap.SampleLevel(g_SamPointClamp, vIn.Tex - float2(du, 0.0f), 0.0f).r;
        float right = g_DisplacementMap.SampleLevel(g_SamPointClamp, vIn.Tex + float2(du, 0.0f), 0.0f).r;
        float top = g_DisplacementMap.SampleLevel(g_SamPointClamp, vIn.Tex - float2(0.0f, dv), 0.0f).r;
        float bottom = g_DisplacementMap.SampleLevel(g_SamPointClamp, vIn.Tex + float2(0.0f, dv), 0.0f).r;
        vIn.NormalL = normalize(float3(-right + left, 2.0f * g_GridSpatialStep, bottom - top));
    }
    
    matrix viewProj = mul(g_View, g_Proj);
    vector posW = mul(float4(vIn.PosL, 1.0f), g_World);

    vOut.PosW = posW.xyz;
    vOut.PosH = mul(posW, viewProj);
    vOut.NormalW = mul(vIn.NormalL, (float3x3) g_WorldInvTranspose);
    vOut.Tex = mul(float4(vIn.Tex, 0.0f, 1.0f), g_TexTransform).xy;
    return vOut;
}

能够看到，顶点y坐标的值和法向量的计算都移步到了顶点着色器上。

由于咱们对位移贴图的采样是要取出与当前顶点相邻的4个顶点对应的4个像素，故不能使用含有线性插值法的采样器来采样。所以咱们还须要在RenderStates.h中添加点过滤+Clamp采样：

ComPtr<ID3D11SamplerState> RenderStates::SSPointClamp = nullptr;    // 采样器状态：点过滤与Clamp模式

// ...
D3D11_SAMPLER_DESC sampDesc;
ZeroMemory(&sampDesc, sizeof(sampDesc));

// 点过滤与Clamp模式
sampDesc.Filter = D3D11_FILTER_MIN_MAG_MIP_POINT;
sampDesc.AddressU = D3D11_TEXTURE_ADDRESS_CLAMP;
sampDesc.AddressV = D3D11_TEXTURE_ADDRESS_CLAMP;
sampDesc.AddressW = D3D11_TEXTURE_ADDRESS_CLAMP;
sampDesc.ComparisonFunc = D3D11_COMPARISON_NEVER;
sampDesc.MinLOD = 0;
sampDesc.MaxLOD = D3D11_FLOAT32_MAX;
HR(device->CreateSamplerState(&sampDesc, SSPointClamp.GetAddressOf()));

GpuWavesRender类

GpuWavesRender类定义以下：

class GpuWavesRender : public WavesRender
{
public:
    GpuWavesRender() = default;
    ~GpuWavesRender() = default;
    // 不容许拷贝，容许移动
    GpuWavesRender(const GpuWavesRender&) = delete;
    GpuWavesRender& operator=(const GpuWavesRender&) = delete;
    GpuWavesRender(GpuWavesRender&&) = default;
    GpuWavesRender& operator=(GpuWavesRender&&) = default;

    // 要求顶点行数和列数都能被16整除，以保证不会有多余
    // 的顶点被划入到新的线程组当中
    HRESULT InitResource(ID3D11Device* device,
        const std::wstring& texFileName,    // 纹理文件名
        UINT rows,          // 顶点行数
        UINT cols,          // 顶点列数
        float texU,         // 纹理坐标U方向最大值
        float texV,         // 纹理坐标V方向最大值
        float timeStep,     // 时间步长
        float spatialStep,  // 空间步长
        float waveSpeed,    // 波速
        float damping,      // 粘性阻尼力
        float flowSpeedX,   // 水流X方向速度
        float flowSpeedY);  // 水流Y方向速度

    void Update(ID3D11DeviceContext* deviceContext, float dt);

    // 在顶点[i][j]处激起高度为magnitude的波浪
    // 仅容许在1 < i < rows和1 < j < cols的范围内激起
    void Disturb(ID3D11DeviceContext* deviceContext, UINT i, UINT j, float magnitude);
    // 绘制水面
    void Draw(ID3D11DeviceContext* deviceContext, BasicEffect& effect);

    void SetDebugObjectName(const std::string& name);

private:
    struct {
        DirectX::XMFLOAT4 waveInfo;
        DirectX::XMINT4 index;
    } m_CBUpdateSettings;                                   // 对应Waves.hlsli的常量缓冲区

private:
    ComPtr<ID3D11Texture2D> m_pNextSolution;                // 缓存下一次模拟结果的y值二维数组
    ComPtr<ID3D11Texture2D> m_pCurrSolution;                // 保存当前模拟结果的y值二维数组
    ComPtr<ID3D11Texture2D> m_pPrevSolution;                // 保存上一次模拟结果的y值二维数组

    ComPtr<ID3D11ShaderResourceView> m_pNextSolutionSRV;    // 缓存下一次模拟结果的y值着色器资源视图
    ComPtr<ID3D11ShaderResourceView> m_pCurrSolutionSRV;    // 缓存当前模拟结果的y值着色器资源视图
    ComPtr<ID3D11ShaderResourceView> m_pPrevSolutionSRV;    // 缓存上一次模拟结果的y值着色器资源视图

    ComPtr<ID3D11UnorderedAccessView> m_pNextSolutionUAV;   // 缓存下一次模拟结果的y值无序访问视图
    ComPtr<ID3D11UnorderedAccessView> m_pCurrSolutionUAV;   // 缓存当前模拟结果的y值无序访问视图
    ComPtr<ID3D11UnorderedAccessView> m_pPrevSolutionUAV;   // 缓存上一次模拟结果的y值无序访问视图

    ComPtr<ID3D11Buffer> m_pVertexBuffer;                   // 当前模拟的顶点缓冲区
    ComPtr<ID3D11Buffer> m_pIndexBuffer;                    // 当前模拟的索引缓冲区
    ComPtr<ID3D11Buffer> m_pConstantBuffer;                 // 当前模拟的常量缓冲区

    ComPtr<ID3D11ComputeShader> m_pWavesUpdateCS;           // 用于计算模拟结果的着色器
    ComPtr<ID3D11ComputeShader> m_pWavesDisturbCS;          // 用于激起水波的着色器

    ComPtr<ID3D11ShaderResourceView> m_pTextureDiffuse;     // 水面纹理
};

其中m_pNextSolution、m_pCurrSolution、m_pPrevSolution都为2D位移贴图，它们不只可能会做为计算着色器的输入、输出（UAV），还可能会做为提供给顶点着色器的位移y输入用于计算（SRV）。

GpuWavesRender::Disturb方法--激起波浪

计算工做都交给GPU了，这里CPU也就负责提供所需的内容，而后再调度计算着色器便可。最后必定要把绑定到CS的UAV撤下来，避免资源同时做为一个地方的输入和另外一个地方的输出。

void GpuWavesRender::Disturb(ID3D11DeviceContext* deviceContext, UINT i, UINT j, float magnitude)
{
    // 更新常量缓冲区
    D3D11_MAPPED_SUBRESOURCE mappedData;
    m_CBUpdateSettings.waveInfo = XMFLOAT4(0.0f, 0.0f, 0.0f, magnitude);
    m_CBUpdateSettings.index = XMINT4(j, i, 0, 0);
    deviceContext->Map(m_pConstantBuffer.Get(), 0, D3D11_MAP_WRITE_DISCARD, 0, &mappedData);
    memcpy_s(mappedData.pData, sizeof m_CBUpdateSettings, &m_CBUpdateSettings, sizeof m_CBUpdateSettings);
    deviceContext->Unmap(m_pConstantBuffer.Get(), 0);
    // 设置计算所需
    deviceContext->CSSetShader(m_pWavesDisturbCS.Get(), nullptr, 0);
    ID3D11UnorderedAccessView* m_UAVs[1] = { m_pCurrSolutionUAV.Get() };
    deviceContext->CSSetConstantBuffers(0, 1, m_pConstantBuffer.GetAddressOf());
    deviceContext->CSSetUnorderedAccessViews(2, 1, m_UAVs, nullptr);

    deviceContext->Dispatch(1, 1, 1);

    // 清除绑定
    m_UAVs[0] = nullptr;
    deviceContext->CSSetUnorderedAccessViews(2, 1, m_UAVs, nullptr);
}

GpuWavesRender::Update方法--更新波浪

须要注意的是，这三个位移贴图是循环使用的。调度完成以后，本来是上一次模拟的纹理将用于等待下一次模拟的输出，而当前模拟的纹理则变成上一次模拟的纹理，下一次模拟的纹理则变成了当前模拟的纹理。这种循环交换方式称之为Ping-Pong交换：

void GpuWavesRender::Update(ID3D11DeviceContext* deviceContext, float dt)
{
    m_AccumulateTime += dt;
    m_TexOffset.x += m_FlowSpeedX * dt;
    m_TexOffset.y += m_FlowSpeedY * dt;

    // 仅仅在累积时间大于时间步长时才更新
    if (m_AccumulateTime > m_TimeStep)
    {
        // 更新常量缓冲区
        D3D11_MAPPED_SUBRESOURCE mappedData;
        m_CBUpdateSettings.waveInfo = XMFLOAT4(m_K1, m_K2, m_K3, 0.0f);
        deviceContext->Map(m_pConstantBuffer.Get(), 0, D3D11_MAP_WRITE_DISCARD, 0, &mappedData);
        memcpy_s(mappedData.pData, sizeof m_CBUpdateSettings, &m_CBUpdateSettings, sizeof m_CBUpdateSettings);
        deviceContext->Unmap(m_pConstantBuffer.Get(), 0);
        // 设置计算所需
        deviceContext->CSSetShader(m_pWavesUpdateCS.Get(), nullptr, 0);
        deviceContext->CSSetConstantBuffers(0, 1, m_pConstantBuffer.GetAddressOf());
        ID3D11UnorderedAccessView* pUAVs[3] = { m_pPrevSolutionUAV.Get(), m_pCurrSolutionUAV.Get(), m_pNextSolutionUAV.Get() };
        deviceContext->CSSetUnorderedAccessViews(0, 3, pUAVs, nullptr);
        // 开始调度
        deviceContext->Dispatch(m_NumCols / 16, m_NumRows / 16, 1);

        // 清除绑定
        pUAVs[0] = pUAVs[1] = pUAVs[2] = nullptr;
        deviceContext->CSSetUnorderedAccessViews(0, 3, pUAVs, nullptr);

        //
        // 对缓冲区进行Ping-pong交换以准备下一次更新
        // 上一次模拟的缓冲区再也不须要，用做下一次模拟的输出缓冲
        // 当前模拟的缓冲区变成上一次模拟的缓冲区
        // 下一次模拟的缓冲区变换当前模拟的缓冲区
        //
        auto resTemp = m_pPrevSolution;
        m_pPrevSolution = m_pCurrSolution;
        m_pCurrSolution = m_pNextSolution;
        m_pNextSolution = resTemp;

        auto srvTemp = m_pPrevSolutionSRV;
        m_pPrevSolutionSRV = m_pCurrSolutionSRV;
        m_pCurrSolutionSRV = m_pNextSolutionSRV;
        m_pNextSolutionSRV = srvTemp;
        
        auto uavTemp = m_pPrevSolutionUAV;
        m_pPrevSolutionUAV = m_pCurrSolutionUAV;
        m_pCurrSolutionUAV = m_pNextSolutionUAV;
        m_pNextSolutionUAV = uavTemp;

        m_AccumulateTime = 0.0f;        // 重置时间
    }
}

GpuWavesRender::Draw方法--绘制波浪

跟CPU的绘制区别基本上就在注释部分了：

void GpuWavesRender::Draw(ID3D11DeviceContext* deviceContext, BasicEffect& effect)
{
    UINT strides[1] = { sizeof(VertexPosNormalTex) };
    UINT offsets[1] = { 0 };
    deviceContext->IASetVertexBuffers(0, 1, m_pVertexBuffer.GetAddressOf(), strides, offsets);
    deviceContext->IASetIndexBuffer(m_pIndexBuffer.Get(), DXGI_FORMAT_R32_UINT, 0);

    effect.SetWavesStates(true, 1.0f / m_NumCols, 1.0f / m_NumCols, m_SpatialStep); // 开启波浪绘制
    effect.SetMaterial(m_Material);
    effect.SetTextureDiffuse(m_pTextureDiffuse.Get());
    effect.SetTextureDisplacement(m_pCurrSolutionSRV.Get());    // 须要额外设置位移贴图
    effect.SetWorldMatrix(XMLoadFloat4x4(&m_WorldMatrix));
    effect.SetTexTransformMatrix(XMMatrixScaling(m_TexU, m_TexV, 1.0f) * XMMatrixTranslationFromVector(XMLoadFloat2(&m_TexOffset)));
    effect.Apply(deviceContext);
    deviceContext->DrawIndexed(m_IndexCount, 0, 0);

    effect.SetTextureDisplacement(nullptr); // 解除占用
    effect.SetWavesStates(false);   // 关闭波浪绘制
    effect.Apply(deviceContext);
}

演示

从这一章开始带着个人光追显卡来跑渲染效果了（相比之前的集成显卡），帧数会有明显的提高

下图演示了GPU绘制的水波效果。因为限制了10M上传，只能勉强看到这一小段了。

测试帧数以下：

	GPU通用计算模式	CPU动态更新模式
Debug x64模式	3500	27
Release x64模式	4200	3900

Debug的话又是vector那边惹出的问题，忽略不计。

只不过当前的样例跟龙书的样例都没有处理好视角带来的过分混合问题，我的暂时也没有解决的办法：

练习题

不要在山体高于水面的顶点激起波浪。
修改波浪的初始参数，观察效果上有什么不一样。

参考资料

本文的算法实现参考的是 Mathematics for 3D Game Programming and Computer Graphics, Third Edition，书内页码443开始。群内能够下载。

DirectX11 With Windows SDK完整目录

Github项目源码

欢迎加入QQ群: 727623616 能够一块儿探讨DX11，以及有什么问题也能够在这里汇报。