DirectX11进阶8_计算着色器(入门、流体模拟)

一、计算着色器：入门

GPU通常被设计为从一个位置或连续的位置读取并处理大量的内存数据(即流操作)，而CPU则被设计为专门处理随机内存的访问。

由于顶点数据和像素数据可以分开处理，GPU架构使得它能够高度并行，在处理图像上效率非常高。但是一些非图像应用程序也能够利用GPU强大的并行计算能力以获得效益。GPU用在非图像用途的应用程序可以称之为：通用GPU(GPGPU)编程。

GPU需要数据并行的算法才能从GPU的并行架构中获得优势，并不是所有的算法都适合用GPU来实现。对于大量的数据，我们需要保证它们都进行相似的操作以确保并行处理。比如顶点着色器都是对大量的顶点数据进行处理，而像素着色器也是对大量的像素片元进行处理。

对于GPGPU编程，用户通常需要从显存中获取运算结果，将其传回CPU。这需要从显存将结果复制到内存中，这样虽然速度会慢一些，但起码还是比直接在CPU运算会快很多。如果是用于图形编程的话倒是可以省掉数据传回CPU的时间，比如说我们要对渲染好的场景再通过计算着色器来进行一次模糊处理。

在Direct3D中，计算着色器也是一个可编程着色器，它并不属于渲染管线的一个直接过程。我们可以通过它对GPU资源进行读写操作，运行的结果通常会保存在Direct3D的资源中，我们可以将它作为结果显示到屏幕，可以给别的地方作为输入使用，甚至也可以将它保存到本地。

1.1 线程和线程组

在GPU编程中，我们编写的着色器程序会同时给大量的线程运行，可以将这些线程按网格来划分成线程组。一个线程组由一个多处理器来执行，如果你的GPU有16个多处理器，你会想要把问题分解成至少16个线程组以保证每个多处理器都工作。为了获取更好的性能，让每个多处理器来处理至少2个线程组是一个比较不错的选择，这样当一个线程组在等待别的资源时就可以先去考虑完成另一个线程组的工作。

每个线程组都会获得共享内存，这样每个线程都可以访问它。但是不同的线程组不能相互访问对方获得的共享内存。

线程同步操作可以在线程组中的线程之间进行，但处于不同线程组的两个线程无法被同步。事实上，我们没有办法控制不同线程组的处理顺序，毕竟线程组可以在不同的多处理器上执行。

一个线程组由N个线程组成。硬件实际上会将这些线程划分成一系列warps(一个warp包含32个线程)，并且一个warp由SIMD32中的多处理器进行处理(32个线程同时执行相同的指令)。在Direct3D中，你可以指定一个线程组不同维度下的大小使得它不是32的倍数，但是出于性能考虑，最好还是把线程组的维度大小设为warp的倍数。

将线程组的大小设为256看起来是个比较好的选择，它适用于大量的硬件情况。修改线程组的大小意味着你还需要修改需要调度的线程组数目。

注意：NVIDIA硬件中，每个warp包含32个线程。而ATI则是每个wavefront包含64个线程。warp或者wavefront的大小可能随后续硬件的升级有所修改。

1.1.1 调度线程组执行计算着色器程序

ID3D11DeviceContext::Dispatch方法–调度线程组执行计算着色器程序：

void ID3D11DeviceContext::Dispatch(
	UINT ThreadGroupCountX,		// [In]X维度下线程组数目
	UINT ThreadGroupCountY,		// [In]Y维度下线程组数目
	UINT ThreadGroupCountZ);	// [In]Z维度下线程组数目

可以看到上面列出了X, Y, Z三个维度，说明线程组本身是可以3维的。当前例子的一个线程组包含了8x8x1个线程，而线程组数目为3x2x1，即我们进行了这样的调用：

m_pd3dDeviceContext->Dispatch(3, 2, 1);

1.2 第一份计算着色器程序

现在我们加载两张图片，并将它混合并将结果输出到一张图片：

下面的这个着色器负责对两个纹理的像素颜色进行分量乘法运算。

Texture2D g_TexA : register(t0);
Texture2D g_TexB : register(t1);

RWTexture2D<float4> g_Output : register(u0);

// 一个线程组中的线程数目。线程可以1维展开，也可以
// 2维或3维排布
[numthreads(16, 16, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
    g_Output[DTid.xy] = g_TexA[DTid.xy] * g_TexB[DTid.xy];
}

上面的代码有如下需要注意的：

• Texture2D仅能作为输入，但RWTexture2D类型支持读写，在本样例中主要是用于输出

• RWTexture2D使用时也需要指定寄存器，u说明使用的是无序访问视图寄存器

• [numthreads(X, Y, Z)]修饰符指定了一个线程组包含的线程数目，以及在3D网格中的布局

• 每个线程都会执行一遍该函数

• SV_DispatchThreadID是当前线程在3D网格中所处的位置，每个线程都有独立的SV_DispatchThreadID

• Texture2D除了使用Sample方法来获取像素外，还支持通过索引的方式来指定像素

如果使用1D纹理，线程修饰符通常为[numthreads(X, 1, 1)]或[numthreads(1, Y, 1)]

如果使用2D纹理，线程修饰符通常为[numthreads(X, Y, 1)]，即第三维度为1

2D纹理X和Y的值会影响你在调度线程组时填充的参数

注意：
 - 在cs_4_x下，一个线程组的最大线程数为768，且Z的最大值为1.
 - 在cs_5_0下，一个线程组的最大线程数为1024，且Z的最大值为64.

1.2.1 纹理输出与无序访问视图

留意上面着色器代码中的类型RWTexture2D，你可以对他进行像素写入，也可以从中读取像素。不过模板参数类型填写就比较讲究了。我们需要保证纹理的数据格式和RWTexture2D的模板参数类型一致，这里使用下表来描述比较常见的纹理数据类型和HLSL类型的对应关系：

此外，UAV不支持DXGI_FORMAT_B8G8R8A8_UNORM

其中unorm float表示的是一个32位无符号的，规格化的浮点数，可以表示范围0到1
而与之对应的snorm float表示的是32位有符号的，规格化的浮点数，可以表示范围-1到1

从上表可以得知DXGI_FORMAT枚举值的后缀要和HLSL的类型对应(浮点型对应浮点型，整型对应整型，规格化浮点型对应规格化浮点型)，否则可能会引发下面的错误(这里举DXGI_FORMAT为unorm，HLSL类型为float的例子)：

D3D11 ERROR: ID3D11DeviceContext::Dispatch: The resource return type for component 0 declared in the shader code (FLOAT) is not compatible with the resource type bound to Unordered Access View slot 0 of the Compute Shader unit (UNORM). This mismatch is invalid if the shader actually uses the view (e.g. it is not skipped due to shader code branching). [ EXECUTION ERROR #2097372: DEVICE_UNORDEREDACCESSVIEW_RETURN_TYPE_MISMATCH]

由于DXGI_FORMAT的部分格式比较紧凑，HLSL中能表示的最小类型通常又比较大。比如DXGI_FORMAT_R16G16B16A16_FLOAT和float4，个人猜测HLSL的类型为了能传递给DXGI_FORMAT，允许做丢失精度的同类型转换。

现在我们回到C++代码，现在需要创建一个2D纹理，然后在此基础上再创建无序访问视图作为着色器输出。

bool InitResource()
{

	HR(CreateDDSTextureFromFile(m_pd3dDevice.Get(), L"Texture\\flare.dds",
		nullptr, m_pTextureInputA.GetAddressOf()));
	HR(CreateDDSTextureFromFile(m_pd3dDevice.Get(), L"Texture\\flarealpha.dds",
		nullptr, m_pTextureInputB.GetAddressOf()));
	
	// 创建用于UAV的纹理，必须是非压缩格式
	D3D11_TEXTURE2D_DESC texDesc;
	texDesc.Width = 512;
	texDesc.Height = 512;
	texDesc.MipLevels = 1;
	texDesc.ArraySize = 1;
	texDesc.Format = DXGI_FORMAT_R32G32B32A32_FLOAT;
	texDesc.SampleDesc.Count = 1;
	texDesc.SampleDesc.Quality = 0;
	texDesc.Usage = D3D11_USAGE_DEFAULT;
	texDesc.BindFlags = D3D11_BIND_SHADER_RESOURCE |
		D3D11_BIND_UNORDERED_ACCESS;
	texDesc.CPUAccessFlags = 0;
	texDesc.MiscFlags = 0;

	HR(m_pd3dDevice->CreateTexture2D(&texDesc, nullptr, m_pTextureOutputA.GetAddressOf()));
	
	texDesc.Format = DXGI_FORMAT_R8G8B8A8_UNORM;
	HR(m_pd3dDevice->CreateTexture2D(&texDesc, nullptr, m_pTextureOutputB.GetAddressOf()));

	// 创建无序访问视图
	D3D11_UNORDERED_ACCESS_VIEW_DESC uavDesc;
	uavDesc.Format = DXGI_FORMAT_R32G32B32A32_FLOAT;
	uavDesc.ViewDimension = D3D11_UAV_DIMENSION_TEXTURE2D;
	uavDesc.Texture2D.MipSlice = 0;
	HR(m_pd3dDevice->CreateUnorderedAccessView(m_pTextureOutputA.Get(), &uavDesc,
		m_pTextureOutputA_UAV.GetAddressOf()));

	uavDesc.Format = DXGI_FORMAT_R8G8B8A8_UNORM;
	HR(m_pd3dDevice->CreateUnorderedAccessView(m_pTextureOutputB.Get(), &uavDesc,
		m_pTextureOutputB_UAV.GetAddressOf()));

	// 创建计算着色器
	ComPtr<ID3DBlob> blob;
	HR(CreateShaderFromFile(L"HLSL\\TextureMul_R32G32B32A32_CS.cso",
		L"HLSL\\TextureMul_R32G32B32A32_CS.hlsl", "CS", "cs_5_0", blob.GetAddressOf()));
	HR(m_pd3dDevice->CreateComputeShader(blob->GetBufferPointer(), blob->GetBufferSize(), nullptr, m_pTextureMul_R32G32B32A32_CS.GetAddressOf()));

	HR(CreateShaderFromFile(L"HLSL\\TextureMul_R8G8B8A8_CS.cso",
		L"HLSL\\TextureMul_R8G8B8A8_CS.hlsl", "CS", "cs_5_0", blob.GetAddressOf()));
	HR(m_pd3dDevice->CreateComputeShader(blob->GetBufferPointer(), blob->GetBufferSize(), nullptr, m_pTextureMul_R8G8B8A8_CS.GetAddressOf()));

	return true;
}

观察上面的代码，如果我们想要让纹理绑定到无序访问视图，就需要提供D3D11_BIND_UNORDERED_ACCESS绑定标签。

注意：如果你还为纹理创建了着色器资源视图，那么UAV和SRV不能同时绑定到渲染管线上。

1.2.2 计算着色阶段设置无序访问视图

ID3D11DeviceContext::CSSetUnorderedAccessViews–计算着色阶段设置无序访问视图

void ID3D11DeviceContext::CSSetUnorderedAccessViews(
	UINT                      StartSlot,						// [In]起始槽，值与寄存器对应
	UINT                      NumUAVs,							// [In]UAV数目
	ID3D11UnorderedAccessView * const *ppUnorderedAccessViews,	// [In]UAV数组
	const UINT                *pUAVInitialCounts				// [In]忽略
);

调度过程实现如下：

void Compute()
{
	assert(m_pd3dImmediateContext);

	m_pd3dImmediateContext->CSSetShaderResources(0, 1, m_pTextureInputA.GetAddressOf());
	m_pd3dImmediateContext->CSSetShaderResources(1, 1, m_pTextureInputB.GetAddressOf());

	// DXGI Format: DXGI_FORMAT_R32G32B32A32_FLOAT
	// Pixel Format: A32B32G32R32
	m_pd3dImmediateContext->CSSetShader(m_pTextureMul_R32G32B32A32_CS.Get(), nullptr, 0);
	m_pd3dImmediateContext->CSSetUnorderedAccessViews(0, 1, m_pTextureOutputA_UAV.GetAddressOf(), nullptr);
	m_pd3dImmediateContext->Dispatch(32, 32, 1);

	// DXGI Format: DXGI_FORMAT_R8G8B8A8_SNORM
	// Pixel Format: A8B8G8R8
	m_pd3dImmediateContext->CSSetShader(m_pTextureMul_R8G8B8A8_CS.Get(), nullptr, 0);
	m_pd3dImmediateContext->CSSetUnorderedAccessViews(0, 1, m_pTextureOutputB_UAV.GetAddressOf(), nullptr);
	m_pd3dImmediateContext->Dispatch(32, 32, 1);

	HR(SaveDDSTextureToFile(m_pd3dImmediateContext.Get(), m_pTextureOutputA.Get(), L"Texture\\flareoutputA.dds"));
	HR(SaveDDSTextureToFile(m_pd3dImmediateContext.Get(), m_pTextureOutputB.Get(), L"Texture\\flareoutputB.dds"));
	
	MessageBox(nullptr, L"请打开Texture文件夹观察输出文件flareoutputA.dds和flareoutputB.dds", L"运行结束", MB_OK);
}

由于我们的位图是512x512x1大小，一个线程组的线程布局为16x16x1，线程组的数目自然就是32x32x1了。如果调度的线程组宽度或高度不够，输出的位图也不完全。而如果提供了过宽或过高的线程组并不会影响运行结果，只是提供的线程组资源过多有些浪费而已。

最后通过ScreenGrab库将纹理保存到文件，就可以结束程序了。

1.3 纹理子资源的采样和索引

从上面的例子可以看到，我们能够使用2D索引来指定纹理的某一像素。如果2D索引越界访问，在计算着色器中是拥有良好定义的：读取越界资源将返回0，尝试写入越界资源的操作将不会执行。

但是上面这种采样只针对mip等级为0的纹理子资源，如果我们想指定其它mip等级的纹理子资源，可以使用mip.operator[][]方法：

R mips.Operator[][](
  in uint mipSlice,		// [In]mip切片值
  in uint2 pos			// [In]2D索引
);

返回值R视纹理数据类型而定。

用法如下：

g_Output.mip[g_MipSlice][DTid.xy] = 
	(unorm float4)(g_TexA.mip[gMipSlice][DTid.xy] * g_TexB.mip[g_MipSlice][DTid.xy]);

不过我们的演示程序用到的纹理Mip等级都为1，这里就不在代码端演示了。

纹理的Sample方法通常情况下你是不知道它具体选择的是哪些Mip等级的纹理子资源来进行采样，具体的行为交给采样器状态来决定。但是我们可以使用SampleLevel方法来指定要对纹理的哪个mip等级的子资源进行采样：

R SampleLevel(
	in SamplerState S,		// [In]采样器状态
	in float2 Location,		// [In]纹理坐标
	in float LOD			// [In]mip等级
);

当LOD为整数时，指定的就是具体某个mip等级的纹理，但如果LOD为浮点数，如3.3f，则会对mip等级为3和4的纹理子资源都进行一次采样，然后根据小数部分进行线性插值求得最终的插值颜色。

用法如下：

float4 texColor = g_Tex.SampleLevel(g_Sam, pIn.Tex, 0.0f);	// 使用第一个mip等级的纹理子资源

还有一种办法则是使用Load方法，它需要传递int3类型的参数，其中xy分量分别对应x和y方向基于0的索引（而不是0.0-1.0），z分量则为要使用的mipmap等级：

float4 texColor = g_Tex.Load(int3(DTid.xy, g_MipSlice))

二、计算着色器：共享内存与同步

2.1 DirectCompute 内存模型

DirectCompute提供了三种内存模型：基于寄存器的内存、设备内存和组内共享内存。不同的内存模型在内存大小、速度、访问方式等地方有所区别。

基于寄存器的内存：它的访问速度非常快，但是寄存器不仅有数目限制，寄存器指向的资源内部也有大小限制。如纹理寄存器（t#），常量缓冲区寄存器（b#），无序访问视图寄存器（u#），临时寄存器（r#或x#）等。而且我们在使用寄存器时也不是直接指定某一个寄存器来使用，而是通过着色器对象（例如tbuffer，它是在GPU内部的内存，因此访问速度特别快）对应到某一寄存器，然后使用该着色器对象来间接使用该寄存器的。而且这些寄存器是随着着色器编译出来后定死了的，因此寄存器的可用情况取决于当前使用的着色器代码。

下面的代码展示了如何声明一个基于寄存器的内存：

tbuffer tb : register(t0)
{
    float weight[256];        // 可以从CPU更新，只读
}

设备内存：通常指的是D3D设备创建出来的资源（如纹理、缓冲区），这些资源可以长期存在，只要引用计数不为0。你可以给这些资源创建很大的内存空间来使用，并且你还可以将它们作为着色器资源或者无序访问资源绑定到寄存器中供使用。当然，这种作为着色器资源的访问速度还是没有直接在寄存器上创建的内存对象来得快，因为它是存储在GPU外部的显存中。尽管这些内存可以通过非常高的内部带宽来访问，但是在请求值和返回值之间也有一个相对较高的延迟。尽管无序访问视图可以用于在设备内存中实现与基于寄存器的内存相同的操作，但当执行频繁的读写操作时，性能将会收到严重影响。此外，由于每个线程都可以通过无序访问视图读取或写入资源中的任何位置，这需要手动同步对资源的访问，也可以使用原子操作，又或者定义一个合理的访问方式避免出现多个线程访问到设备内存的同一个数据。

组内共享内存：前面两种内存模型是所有可编程着色阶段都可使用的，但是group shared memory只能在计算着色器使用。它的访问速度比设备内存资源快些，比寄存器慢，但是也有明显的内存限制——每个线程组最多只能分配32KB内存，供内部所有线程使用。组内共享的内存必须确定线程将如何与内存交互和使用内存，因此它还必须同步对该内存的访问。这将取决于正在实现的算法，但它通常涉及到前面描述的线程寻址。

这三种类型的内存提供了不同的访问速度和可用的大小，使得它们可以用于与其能力相匹配的不同情况，这也给计算着色器提供了更大的内存操作灵活性。下表则是对内存模型的总结：

2.2 线程标识符与线程寻址机制

对于线程组(大小(ThreadDimX, ThreadDimY, ThreadDimZ))中的每一个线程，它们都有一个唯一的线程ID值。我们可以使用系统值SV_GroupThreadID来取得，它的索引范围为(0, 0, 0)到(ThreadDimX - 1, ThreadDimY - 1, ThreadDimZ - 1)。

而对于整个线程组来说，由于线程组集合也是在3D空间中排布，它们也有一个唯一的线程组ID值。我们可以使用系统值SV_GroupID来取得，线程组的索引范围取决于调用ID3D11DeviceContext::Dispatch时提供的线程组(大小(GroupDimX, GroupDimY, GroupDimZ))，范围为(0, 0, 0)到(GroupDimX - 1, GroupDimY - 1, GroupDimZ - 1)。

紧接着就是系统值SV_GroupIndex，它是单个线程组内的线程三维索引的一维展开。若已知线程组的大小为(ThreadDimX, ThreadDimY, ThreadDimZ)，则可以确定SV_GroupIndex与SV_GroupThreadID满足下面关系：

SV_GroupIndex = SV_GroupThreadID.z * ThreadDimX * ThreadDimY + SV_GroupThreadID.y * ThreadDimX + SV_GroupThreadID.x;

最后就是系统值SV_DispatchThreadID，线程组中的每一个线程在ID3D11DeviceContext::Dispatch提供的线程组集合中都有其唯一的线程ID值。若已知线程组的大小为 (ThreadDimX, ThreadDimY, ThreadDimZ)，则可以确定SV_DispatchThreadID，SV_GroupThreadID和SV_GroupID满足以下关系：

SV_DispatchThreadID.xyz = SV_GroupID.xyz * float3(ThreadDimX, ThreadDimY, ThreadDimZ) + SV_GroupThreadID.xyz;

2.3 共享内存和线程同步

由于大量线程同时运行，并且线程能够通过组内共享内存或通过无序访问视图对应的资源进行交互，因此需要能够同步线程之间的内存访问。与传统的多线程编程一样，许多线程可用读取和写入相同的内存位置，存在写后读（Read After Write，简称RAW）导致内存损坏的危险。如何在不损失GPU并行性带来的性能的情况下还能够高效地同步这么多线程？幸运的是，有几种不同的机制可用用于同步线程组内的线程。

2.3.1 内存屏障（Memory Barriers）

这是一种最高级的同步技术。HLSL提供了许多内置函数，可用于同步线程组中所有线程的内存访问。需要注意的是，它只同步线程组中的线程，而不是整个调度。这些函数有两个不同的属性。第一个是调用函数时线程正在同步的内存类别（设备内存、组内共享内存，还是两者都有），第二个则指定给定线程组中的所有线程是否同步到其执行过程中的同一处。根据这两个属性，衍生出了下面这些不同版本的内置函数：

这些函数中都会阻止线程继续，直到满足该函数的特定条件位置。其中第一个函数GroupMemoryBarrior()阻塞线程的执行，直到线程组中的所有线程对组内共享内存的所有写入都完成。这用于确保当线程在组内共享内存中彼此共享数据时，所需的值在被其他线程读取之前有机会写入组内共享内存。这里有一个很重要的区别，即着色器核心执行一个写指令，而那个指令实际上是由GPU的内存系统执行的，并且写入内存中，然后在内存中它将再次对其他线程可用。从开始写入值到完成写入到目标位置有一个可变的时间量，这取决于硬件实现。通过执行阻塞操作，直到这些写操作被保证已经完成，开发人员可以确定不会有任何写后读错误引发的问题。

不过话说了那么多，总得实践一下。个人将双调排序项目中BitonicSort_CS.hlsl第15行的GroupMemoryBarrierWithGroupSync()修改为GroupMemoryBarrier()，执行后发现多次运行程序会出现一例排序结果不一致的情况。因此可以这样判断：GroupMemoryBarrier()仅在线程组内的所有线程组存在线程写入操作时阻塞，因此可能会出现阻塞结束时绝大多数线程完成了共享数据写入，但仍有少量线程甚至还没开始写入共享数据。因此实际上很少能够见到他出场的机会。

然后是GroupMemoryBarriorWithGroupSync()函数，相比上一个函数，他还阻止那些先到该函数的线程执行，直到所有的线程都到达该函数才能继续。很明显，在所有组内共享内存都加载之前，我们不希望任何线程前进，这使它成为完美的同步方法。

而第二对同步函数也执行类似的操作，只不过它们是在设备内存池上操作。这意味着在继续执行着色器程序前，可以同步通过无序访问视图写入资源的所有挂起内存的写入操作。这对于同步更大数目的内存更为有用，如果所需的共享存储器的大小太大不适合用组内共享内存，则可以将数据存在更大的设备内存的资源中。

第三对同步函数则是一起执行前面两种类型的同步，用于同时存在共享内存和设备内存的访问和同步上。

在一个线程组内，允许设置一片共享内存区域，使得当前线程组内的所有线程都可以访问当前的共享内存。一旦设置，那么每个线程都会各自配备一份共享内存。共享内存的访问速度非常快，就像寄存器访问CPU缓存那样。

共享内存的声明方式如下：

groupshared float4 g_Cache[256];

对于每个线程组来说，它所允许分配的总空间最大为32kb(即8192个标量，或2048个向量)。内部线程通常应该使用SV_ThreadGroupID来写入共享内存，这样以保证每个线程不会出现重复写入操作，而读取共享内存一般是线程安全的。

分配太多的共享内存会导致性能问题。假如一个多处理器支持32kb的共享内存，然后你的计算着色器需要20kb的共享内存，这意味着一个多处理器只适合处理一个线程组，因为剩余的共享内存不足以给新的线程组运行，这也会限制GPU的并行运算，当该线程组因为某些原因需要等待，会导致当前的多处理器处于闲置状态。因此保证一个多处理器至少能够处理两个或以上的线程组(比如每个线程组分配16kb以下的共享内存)，以尽可能减少该多处理器的闲置时间。

现在来考虑下面的代码：

Texture2D g_Input : register(t0);
RWTexture2D<float4> g_Output : register(u0);

groupshared float4 g_Cache[256];

[numthreads(256, 1, 1)]
void CS(uint3 GTid : SV_GroupThreadID,
	uint3 DTid : SV_DispatchThreadID)
{
	// 将纹理像素值缓存到共享内存
	g_Cache[GTid.x] = g_Input[DTid.xy];
	
	// 取出共享内存的值进行计算
	
	// 注意！！相邻的两个线程可能没有完成对纹理的采样
	// 以及存储到共享内存的操作
	float left = g_Cache[GTid.x - 1];
	float right = g_Cache[GTid.x + 1];
	
	// ...
}

因为多个线程同时运行，同一时间各个线程当前执行的指令有所偏差，有的线程可能已经完成了共享内存的赋值操作，有的线程可能还在进行纹理采样操作。如果当前线程正在读取相邻的共享内存片段，结果将是未定义的。为了解决这个问题，我们必须在读取共享内存之前让当前线程等待线程组内其它的所有线程完成写入操作。这里我们可以使用GroupMemoryBarrierWithGroupSync函数：

Texture2D g_Input : register(t0);
RWTexture2D<float4> g_Output : register(u0);

groupshared float4 g_Cache[256];

[numthreads(256, 1, 1)]
void CS(uint3 GTid : SV_GroupThreadID,
	uint3 DTid : SV_DispatchThreadID)
{
	// 将纹理像素值缓存到共享内存
	g_Cache[GTid.x] = g_Input[DTid.xy];
	
	// 等待所有线程完成写入
	GroupMemoryBarrierWithGroupSync();
	
	// 现在读取操作是线程安全的，可以开始进行计算
	float left = g_Cache[GTid.x - 1];
	float right = g_Cache[GTid.x + 1];
	
	// ...
}

三、水波实现

3.1 生成三角形栅格表示地形或水面

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

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

然后我们再运用函数y=f(x,z) 来为每个栅格点获取对应的y坐标。再利用*(x,f(x,z),z)* 的所有顶点构造出地形栅格

3.1.1 生成栅格顶点

由以上分析可知，我们首先需要完成的就是来构建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 = depth / 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平面，然后再得到对应的未经标准化的切线向量，如下图所示：

3.1.2 山峰函数

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

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

其中y对x和对z的偏导分别为：

山峰体面mesh如下：

因此创建上述地形的函数可以写成：

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) }; })	// 法向量函数

3.2 流体模拟

在许多游戏中，你可能会看到有水面流动的场景，实际上他们不一定是真正的流体，而有可能只是一个正在运动的水体表面。出于运行效率的考虑，这些效果的实现或多或少涉及到数学公式或者物理公式。
本节我们只讨论龙书所实现的方法，即使用波动方程来表现局部位置激起的水波。它的公式推导比较复杂，我们此处仅给出推导结果并运用。

3.2.1 波动方程

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

一维的波动方程(c^2=T/ρ)：

同理，二维的波动方程可以通过添加一个y项得到：

常数c具有单位时间距离的尺度，因此可以表示速度。事实上我们也不会证明c实际上就是波沿着绳索或表面传递的速度。这是有意义的，因为波的速度随介质所受张力T的变大而增加，随介质密度μ的减小而减小。

满足一维波动方程的解有无穷多个，例如一种常见的波函数形式为z(x,t)=Asin(ω(t−xv))，函数随着时间的推移图像如下：

然而二维的波动方程仅包含了张力，没有考虑到其他阻力因素，这导致波的平均振幅并不会有任何损耗。我们可以给方程组添加一个与张力方向相反，且与点的运动速度有关的粘性阻尼力：

其中非负实数μ代表了液体的粘性，用来控制水面的波什么时候能够平静下来。μ越大，水波消亡的速度越快。对于水来说，通常它的μ值会比较小，使得水波能够存留比较长的时间；但对于油来说，它的μ值会比较大一些，因此水波消亡的速度会比较快。

3.2.2 波动方程近似导数

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

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

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

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


同理，我们取顶点[i][j]的两个y方向上的相邻顶点[i][j-1]和[i][j+1]来近似计算z对于y的偏导数：

至于时间，我们可以通过计算顶点在当前时刻分别与上一个时刻和下一个时刻的平均位移差来定义z对时间t偏导Δz/Δt：

二阶偏导也可以使用和一阶偏导相同的方法来计算。假如我们已经计算出了顶点[i-1][j]和[i+1][j]的偏移z对x的一阶偏导数，那么我们就可以得到两者平均差值：

将方程带入上式，可以得到：

除以d使得我们可以得到Δ(∂z/∂x)/Δx，对应二阶偏导：

该公式要求我们使用x轴距离为2的顶点[i+2][j]和[i-2][j]来计算二阶偏导。不过相邻的两个顶点就没有用到了，我们可以基于顶点[i][j]将x轴缩小一半，使用距离最近的两个相邻点[i+1][j]和[i-1][j]来计算二阶偏导：

同理可得：

3.2.3 求出水面位移

联立z对t的一阶偏导以及二阶偏导，带粘性阻尼力的二维波动方程可以表示为：

我们想要能够在传递模拟间隔t来确定下一次模拟位移量z(i,j,k+1)，现在我们已经知道当前位移量z(i,j,k)和上一次模拟的位移量z(i,j,k−1)。因此z(i,j,k+1)的解为：

这条公式正是我们想要的。其中常量部分可以预先计算出来，只剩下3个带t的因式和4个加法需要给网格中每个顶点进行运算。

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

假定我们拥有一个n×m顶点数组(其任意z(i,j,0)=0和z(i,j,1)=0)，现在让某处顶点被抬高使得z(i0,j0,0)=h和z(i0,j0,1)=h，h是一个非0位移值。若该处顶点被释放了2t时间，此时式子(28.25)中第三个加法部分的值为0，故有：

为了让顶点向周围平坦的水面移动，它在2t时刻的位移必须要比在t时刻的更小一些。因此就有：

代入方程，有：

因此，

把速度c解出来，我们可以得到：

这告诉我们对于公式，给定两个顶点的距离d以及时间间隔t，波速c必须小于上式的最大值

又或者，给定距离d和波速c，我们能够计算出最大的时间间隔t。对不等式乘上(−μt+2)来简化得到：

由于t>0，中间部分恒大于0，去掉左半部分解一元二次不等式，并舍去t<=0的部分，得：

tip： 使用c区间和t区间外的值会导致位移z的结果呈指数级爆炸。

3.2.4 求出顶点法向量和切线向量

现在我们需要准备两个二维顶点位置数组，其中一个数组表示的是当前模拟所有顶点的位置，而另一个数组则表示的是上一次模拟所有顶点的位置。当我们计算新的位移时，将下一次模拟的结果直接写在存放上一次模拟顶点的数组即可（可以观察位移公式，我们需要的是当前顶点上一次模拟、当前模拟的数据，以及当前模拟相邻四个顶点的数据，因此其他顶点在计算位移量的时候不会有影响）。

为了执行光照计算，我们还需要为每个顶点获取正确的法向量以及可能正确的切线向量。对于顶点坐标(i, j)，未标准化的，与x轴对齐的切向量T和与y轴对齐的切线向量B如下：

用偏导公式代入，可得：

经过叉乘后可以得到未经标准化的法向量：

对上述向量乘上2d的倍数并不会改变它的方向，但可以消除除法：

注意这里T和B并没有相互正交。

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

3.3 代码实现

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

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

其中前者的效率一般不如后者，但实现起来比较简单。

3.3.1 CPU计算实现

3.3.1.1 基类WavesRender

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

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

	void SetMaterial(const Material& material);

	Transform& GetTransform();
	const Transform& GetTransform() const;

	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 = 0;					// 顶点行数
	UINT m_NumCols = 0;					// 顶点列数

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

	Transform m_Transform = {};			// 水面变换
	DirectX::XMFLOAT2 m_TexOffset = {};	// 纹理坐标偏移
	float m_TexU = 0.0f;				// 纹理坐标U方向最大值
	float m_TexV = 0.0f;				// 纹理坐标V方向最大值
	Material m_Material = {};			// 水面材质

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

	//
	// 我们可以预先计算出来的常量
	//

	float m_K1 = 0.0f;
	float m_K2 = 0.0f;
	float m_K3 = 0.0f;
};

Init方法将水波位移公式的三个重要参数先初始化，以供后续计算使用。

void WavesRender::Init(UINT rows, UINT cols, float texU, float texV,
	float timeStep, float spatialStep, float waveSpeed, float damping, float flowSpeedX, float flowSpeedY)
{
	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;
}

3.3.1.2 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 = false;							// 当前是否有顶点数据更新
};

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

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

3.3.1.3 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;
}

3.3.1.4 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);
			}
		}	
	}
}

3.3.1.5 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(m_Transform.GetLocalToWorldMatrixXM());
	effect.SetTexTransformMatrix(XMMatrixScaling(m_TexU, m_TexV, 1.0f) * XMMatrixTranslationFromVector(XMLoadFloat2(&m_TexOffset)));
	effect.Apply(deviceContext);
	deviceContext->DrawIndexed(m_IndexCount, 0, 0);
	
}

3.3.2 GPU计算实现

相比CPU计算法，GPU计算法的实现则更为复杂了。因为它不仅需要用到计算着色器，还需要利用位移贴图处理水面效果。

3.3.2.1 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()));

3.3.2.2 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）。

3.3.2.3 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);
}

3.3.2.4 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;		// 重置时间
	}
}

3.3.2.5 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(m_Transform.GetLocalToWorldMatrixXM());
	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);
}