65.9K
CodeProject 正在变化。 阅读更多。
Home

GPGPU Papyrus 演示

John Michael Hauck
starIconstarIconstarIconstarIconstarIcon

5.00/5 (12投票s)

2013年4月24日

CPOL

9分钟阅读

viewsIcon

25994

对于 C# 桌面开发人员来说，以前从未如此容易地编写利用现代图形卡惊人的计算性能的代码。在这篇文章中，我将分享一些解决一个简单（但仍然有趣）的图像分析问题的技术。源代码 https://www.assembla.com/co

对于 C# 桌面开发人员来说，以前从未如此容易地编写利用现代图形卡惊人的计算性能的代码。在这篇文章中，我将分享一些解决一个简单（但仍然有趣）的图像分析问题的技术。

ScreenShot

https://www.assembla.com/code/cudafypapyrus/subversion/nodes

构建应用程序

要重现以下工作，您需要一张支持 CUDA 或 OpenCL 的显卡。我在 GTX Titan (1000 美元) 和 GT 230 (100 美元) 上都可以很好地运行此代码，并且支持 CUDA 和 OpenCL 模式。

  1. 获取软件
    1. 对于 Nvidia GPU 的 CUDA，请从 Nvidia 获取 CUDA：https://developer.nvidia.com/cuda-downloads  
    2. 对于 GPU 的 OpenCL，您的视频驱动程序应该就足够了。
    3. 对于 Intel CPU 的 OpenCL，请获取此 SDK：http://software.intel.com/en-us/vcsource/tools/opencl-sdk
    4. 对于 CUDAfy（C# 开发），请获取此 DLL：http://cudafy.codeplex.com/。截至撰写本文时，如果您想在 CUDAfy 中获得 OpenCL 支持，您需要获取 1.23 版（alpha）。
    5. 为了开发，请获取 Visual Studio。如果您在进行 CUDA 工作，您需要在您的硬盘上安装 2010 C++ 编译器。您可以使用 VS2010 或 VS2012 来完成本文中的所有工作。
  2. 对于 CUDA 工作，请确保您的环境变量 PATH 包含指向 VS2010 C++ 编译器的链接。我的包含此字符串：C:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\bin
  3. 向应用程序添加对 Cudafy.net 的引用。

问题

恢复一张被分成 9 个图块的图像，这些图块已被随机重新排列。这是一段被打乱的卷轴（包含以赛亚书 6:10）的片段。

Scrambled Isaiah_6_10

缩放图像

第一步是将图像缩放到可管理的尺寸。在这种情况下，我决定将缩放后的图像尺寸设置为 192x192 像素。我还决定去除颜色分量，直接处理灰度图像。

Scaled

在 CPU 上运行的 C# 代码如下所示：

// load image from disk
var cpuImage = new MyImageReader(sourceFileName);
 
// allocate image memory on the GPU
var gpuImage = Gpu.Allocate<uint>(cpuImage.Pixels.Length);
 
// copy the image pixels to the GPU
Gpu.CopyToDevice(cpuImage.Pixels, gpuImage);
 
// allocate scaled image memory on the GPU
var gpuScaledImage = Gpu.Allocate<byte>(64 * 3, 64 * 3);
 
// rescale the image using the GPU
Gpu.Launch(new dim3(12, 12), new dim3(16, 16),
    ScaleImageKernel, gpuImage, cpuImage.Width, cpuImage.Height, gpuScaledImage);

最后一行代码在 GPU 上启动一个名为 ScaleImageKernel 的函数。Gpu.Launch 的前两个参数指定了块的数量（12x12）和每个块中的线程数（16x16）。这样，我们就得到了每个像素一个线程来处理缩放后的图像。

ScaleDiag

您可能会想，为什么我们不启动一个（无论是什么）块，并包含 192x192 个线程呢？答案很简单。CUDA 每个块最多只允许 1024 个线程。那么，您可能会想，为什么我们不启动 192x192 个块，每个块只有一个线程呢？答案也很简单。每个块都会带来开销，您希望将其最小化。您可能会想，对于这个问题，最佳的块数和每个块的线程数是多少。Nvidia 提供了一个名为 CUDA_Occupancy_Calculator.xls 的精美电子表格，可以提供帮助。另一个选择是尝试不同的尺寸，看看什么效果最好。

剩余的参数将传递给 GPU 上的 ScaleImageKernel 函数。需要明确的是，ScaleImageKernel 函数会针对 192x192 个线程中的每一个被调用一次。经过 CUDAfy 转换后在 GPU 上运行的 C# 代码如下所示：

[Cudafy]
public static void ScaleImageKernel(GThread gThread, uint[] sourceImage, int sourceWidth, int sourceHeight, byte[,] scaledImage)
{
    var scaledX = gThread.blockIdx.x * gThread.blockDim.x + gThread.threadIdx.x;
    var scaledY = gThread.blockIdx.y * gThread.blockDim.y + gThread.threadIdx.y;
 
    var sourceX = scaledX * sourceWidth / scaledImage.GetLength(0);
    var sourceY = scaledY * sourceHeight / scaledImage.GetLength(1);
 
    var sourcePixel = sourceImage[sourceX + sourceY * sourceWidth];
    var blue = sourcePixel & 0xFF;
    var green = (sourcePixel >> 8) & 0xFF;
    var red = (sourcePixel >> 16) & 0xFF;
 
    scaledImage[scaledX, scaledY] = (byte)(red * 0.3f + green * 0.6f + blue * 0.1f);
}

对于每一次调用，都可以通过当前的块和线程索引直接获得缩放后的目标像素。这些值通过一个特殊的 GThread 类型参数传递给函数。您应该能够轻松理解代码，包括转换为灰度。

您可能注意到，这个缩减内核并没有对各个源像素进行平均，而是简单地复制了其中一个源像素而忽略了其他像素。好吧，最好从简单开始，但现在我们可以介绍一个名为 EnhancedScaleImageKernel 的增强版内核。

[Cudafy]
public static void EnhancedScaleImageKernel(GThread gThread, uint[] sourceImage, int sourceWidth, int sourceHeight, byte[,] scaledImage)
{
    var scaledX = gThread.blockIdx.x * gThread.blockDim.x + gThread.threadIdx.x;
    var scaledY = gThread.blockIdx.y * gThread.blockDim.y + gThread.threadIdx.y;
 
    EnhancedScaleImagePixel(sourceImage, sourceWidth, sourceHeight, scaledImage, scaledX, scaledY);
}
 
[Cudafy]
private static float EnhancedScaleImagePixel(uint[] sourceImage, int sourceWidth, int sourceHeight, byte[,] scaledImage, int scaledX, int scaledY)
{
    var startX = scaledX * sourceWidth / scaledImage.GetLength(0);
    var startY = scaledY * sourceHeight / scaledImage.GetLength(1);
 
    var endX = (scaledX + 1) * sourceWidth / scaledImage.GetLength(0);
    var endY = (scaledY + 1) * sourceHeight / scaledImage.GetLength(1);
 
    var sum = 0f;
    var count = 0;
 
    for (var sourceX = startX; sourceX < endX; sourceX++)
    {
        for (var sourceY = startY; sourceY < endY; sourceY++)
        {
            var sourcePixel = sourceImage[sourceX + sourceY * sourceWidth];
            var blue = sourcePixel & 0xFF;
            var green = (sourcePixel >> 8) & 0xFF;
            var red = (sourcePixel >> 16) & 0xFF;
            sum += red * 0.3f + green * 0.6f + blue * 0.1f;
            count++;
        }
    }
 
    scaledImage[scaledX, scaledY] = (byte)(sum / count);
 
    return 0;
}
出于稍后将解释的原因，我决定将内核的实现放在一个名为 EnhancedScaleImagePixel 的单独函数中。无论如何，这就完成了……一个通过平均源像素并将颜色分量去除来将图像缩放到 192x192 的内核。

与 CPU 对比

现在是“稍后”了。我想比较 GPU 和 CPU 的速度。由于其编写方式，我现在可以从 CPU 的并行循环中调用它来测量差异。
Parallel.For(0, 192 * 192, p => {
    var x = p % 192;
    var y = p / 192;
    EnhancedScaleImagePixel(sourceImage, sourceWidth, sourceHeight, scaledImage, x, y);});
令人惊讶的是，对于我的配置，GPU 完成这项任务的速度比 CPU 快 100 多倍……我的配置是 Dell Precision T3600，16GB RAM，Intel Xeon E5-2665 0 @ 2.40GHz，Nvidia GTX Titan。

提取边缘

出于我的明智考虑，我决定解决这个问题的最佳方法是比较 9 个图块中每个图块的边缘像素，以确定图块的最佳排列。这意味着我可以省去大部分缩放图像的工作。这有点令人遗憾，因为它很有趣。
Edges
我们有 9 个图块，每个图块有 4 个边缘，每个边缘有 64 个像素。在 CPU 上运行的 C# 代码如下所示：
// allocate edges memory on the GPU
var gpuEdges = Gpu.Allocate<byte>(9, 4, 64);
 
// extract edge information using the GPU
Gpu.Launch(new dim3(9, 4), 64, ExtractEdgeKernel, gpuScaledImage, gpuEdges);
正如您所见，我们在 GPU 上创建了 9x4x64 个线程来执行复制这些像素的任务。您可能会错过的是，上一步缩放的图像仍然在 GPU 内存中。这意味着我们不需要传输缩放的图像。EdgeDiag 在 GPU 上运行的 C# 代码如下所示：
[Cudafy]
public static void ExtractEdgeKernel(GThread gThread, byte[,] image, byte[, ,] edges)
{
    var tileIndex = gThread.blockIdx.x;
    var tileX = tileIndex % 3;
    var tileY = tileIndex / 3;
    var edgeIndex = gThread.blockIdx.y;
    var pixelIndex = gThread.threadIdx.x;
 
    var sourceX = tileX * 64;
    var sourceY = tileY * 64;
 
    switch (edgeIndex)
    {
        case 0: sourceY += pixelIndex; break; // left
        case 1: sourceX += pixelIndex; break; // top
        case 2: sourceY += pixelIndex; sourceX += 63; break; // right
        case 3: sourceX += pixelIndex; sourceY += 63; break; // bottom
    }
 
    edges[tileIndex, edgeIndex, pixelIndex] = image[sourceX, sourceY];

 
}
这个内核中最难的部分是弄清楚源坐标映射到哪个边缘坐标。

常量内存

内核完成后，我们就知道会有许多后续的边缘数据操作：事实证明，在 CUDA 中有许多类型的内存。到目前为止，我们一直在使用所谓的“设备内存”。这就像 CPU RAM。它速度很快，但远不及 GPU 上其他类型的内存快。Edges 2
还有另一种类型的内存称为“常量内存”。GPU 可以比设备内存更快地读取这种内存，因为它被缓存了。GPU 无法写入它。此外，常量内存供应稀少——只有 64KB。这是一个使用常量内存的绝佳时机。我们将把上面计算出的边缘数据移到常量内存中，以使后续的内核调用速度更快。要在 CUDAfy 中声明常量内存，您可以输入类似以下内容：
[Cudafy]
public static byte[, ,] Edges = new byte[9, 4, 64];
这里奇怪的是，这实际上在 CPU 可访问的 RAM 和 GPU 设备内存中分配了内存。到目前为止，我们会在内存变量前面加上“gpu”或“cpu”以保持清晰。如果您从 CPU 引用 Edges，它会引用 CPU RAM；如果您从 GPU 引用 Edges，它会引用 GPU 上的常量内存。这应该足够清晰了。
Edges 5
现在，将数据从设备内存移动到常量内存需要两个步骤。首先，我们需要将数据从 GPU 的设备内存复制到 CPU RAM，然后从 CPU RAM 复制回 GPU 上的常量内存。开始吧：
// copy edges to GPU constant memory
Gpu.CopyFromDevice(gpuEdges, Edges);
Gpu.CopyToConstantMemory(Edges, Edges);

计算拟合度

这里的目标是预先计算每个图块边缘之间的拟合度。这意味着我们想知道图块 1 的左边缘与图块 2 的右边缘的匹配程度等等。FitsDiag
为了保存这些信息，我们创建了更多常量内存，如下所示：
[Cudafy]
public static float[,] LeftRightFit = new float[9, 9];
然后，我们使用以下 CPU 代码计算所有左-右拟合的可能性，并将结果存储到常量内存中：
// allocate fit memory on the GPU
var gpuFit = Gpu.Allocate<float>(9, 9);
 
// compare edge fitting using the GPU
Gpu.Launch(new dim3(9, 9), 64, ComputeFitsKernel, 2, 0, gpuFit);
 
// copy edges to GPU constant memory
Gpu.CopyFromDevice(gpuFit, LeftRightFit);
Gpu.CopyToConstantMemory(LeftRightFit, LeftRightFit);
这里没有新东西，只是需要意识到我们为每个边缘对启动了一个块。是的，我们最终会将一个图块的左边缘与自己的右边缘进行比较。这类优化可以通过启动 8x8 网格和一些棘手的内核数学来实现。另一种选择是启动 9x9，但然后让块中的线程立即返回一个糟糕的匹配（如果两个图块相同）。无论哪种方式都很好。本文中未显示任何优化。此时，我想让您期待内核中的一个新问题，将在下一节中展示。ComputeFitsKernel 以每块 64 个线程启动，但我们只想从每个块返回一个数字。在之前的内核中，每个线程都设置了一个特定的输出值。现在，我们有 64 个线程需要协同工作才能得到一个数字。不出所料，块有一些巧妙的功能使这项工作成为可能。执行实际拟合计算的 GPU 内核如下所示：
[Cudafy]
public static void ComputeFitsKernel(GThread gThread, int edgeIndexA, int edgeIndexB, float[,] fit)
{
    var sum = gThread.AllocateShared<float>("sum", 64);
 
    var tileIndexA = gThread.blockIdx.x;
    var tileIndexB = gThread.blockIdx.y;
    var pixelIndex = gThread.threadIdx.x;
 
    var diff = Edges[tileIndexA, edgeIndexA, pixelIndex] - Edges[tileIndexB, edgeIndexB, pixelIndex];
    sum[pixelIndex] = diff * diff;
 
    gThread.SyncThreads();
 
    for (var i = 64 / 2; i > 0; i /= 2)
    {
        if (pixelIndex < i)
        {
            sum[pixelIndex] += sum[pixelIndex + i];
        }
 
        gThread.SyncThreads();
    }
 
    if (pixelIndex == 0)
    {
        fit[tileIndexA, tileIndexB] = sum[0];
    }
}
第一行应该让您困惑。为什么我们要让 64 个线程中的每个线程分配 64 个浮点数的“共享”内存（那是什么）？我相信函数名称不正确。函数应该命名为“TheFirstThreadPerBlockThatGetHereShouldAllocateThisMemoryAndAllSusequentThreadsInTheBlockShouldObtainAReferenceToThisAllocatedMemory”。也许 ObtainShared 也可以。共享内存（每块 64 个浮点数）速度非常非常快。它的速度与 CPU 寄存器相当，比设备内存快约 100-200 倍。每个线程计算一个图块的一个边缘的一个像素与另一个图块的另一个边缘的另一个像素之间的差值。这个差值（平方）被放入共享内存的相应位置。到目前为止，一切都很好。SyncThreads 函数使块中的所有 64 个线程在此处互相等待。这意味着以下 for 循环在所有线程将数据放入共享数组之前不会开始。这是一个告诉您，尽管别人可能说了什么，但在不同块之间同步线程从来没有一个好方法。缩减循环非常简单。阅读代码并理解其优雅之处。还要理解，如果有人决定将 SyncThreads 放在“if”语句中，可能会发生的死锁。那永远是个坏主意，并且可能导致您的设备驱动程序锁定。以上所有内容都为顶部-底部边缘重复。比较排列 下一步是分析所有（9!）排列，并选择拟合度最好的（最小的度量）。
PermDiag
在这种情况下，我们意识到存在大量的排列。我们将（随意）选择每个块包含 256 个线程。然后，我们将创建足够的块（1418 个），以便每个排列都有自己的线程。然后，我们将编写内核，以便每个块返回 256 个线程评估的排列中的最佳排列。然后，我们将这 1418 个最佳候选者从 GPU 传输到 CPU。然后，CPU 将找到这 1418 个评估中的每一个的最佳排列。这是 CPU 代码：
// evaluate all permutations
const int threads = 256;
const int blocks = Permutations / threads + 1;
var cpuEvaluations = new Evaluation[blocks];
var gpuEvaluations = Gpu.Allocate(cpuEvaluations);
Gpu.Launch(blocks, threads, ExplorePermutationsKernel, gpuEvaluations);
Gpu.CopyFromDevice(gpuEvaluations, cpuEvaluations);
 
// get the best permutation
var bestEvaluation = cpuEvaluations[0];
foreach (var evaluation in cpuEvaluations)
{
    if (evaluation.Metric < bestEvaluation.Metric)
    {
        bestEvaluation = evaluation;
    }
}
这个最后的内核在技术上没有什么新的，但它无疑是本项目中最复杂的内核。
[Cudafy]
public static void ExplorePermutationsKernel(GThread gThread, Evaluation[] evaluations)
{
    var blockEvaluations = gThread.AllocateShared<Evaluation>("be", 256);
    var v = gThread.AllocateShared<byte>("v", 256, 9);
    var t = gThread.threadIdx.x;
 
    var permutation = gThread.blockIdx.x * gThread.blockDim.x + gThread.threadIdx.x;
 
    // 0 1 2
    // 3 4 5
    // 6 7 8
 
    TileOrderFromPermutation(Permutations, permutation, 9, v, t);
 
    var metric = 0f;
 
    metric += LeftRightFit[v[t, 0], v[t, 1]] + LeftRightFit[v[t, 1], v[t, 2]];
    metric += LeftRightFit[v[t, 3], v[t, 4]] + LeftRightFit[v[t, 4], v[t, 5]];
    metric += LeftRightFit[v[t, 6], v[t, 7]] + LeftRightFit[v[t, 7], v[t, 8]];
 
    metric += TopBottomFit[v[t, 0], v[t, 3]] + TopBottomFit[v[t, 3], v[t, 6]];
    metric += TopBottomFit[v[t, 1], v[t, 4]] + TopBottomFit[v[t, 4], v[t, 7]];
    metric += TopBottomFit[v[t, 2], v[t, 5]] + TopBottomFit[v[t, 5], v[t, 8]];
 
    blockEvaluations[t].Permutation = permutation;
    blockEvaluations[t].Metric = metric;
 
    gThread.SyncThreads();
 
    for (var i = 256 / 2; i > 0; i /= 2)
    {
        if (t < i)
        {
            if (blockEvaluations[t].Metric > blockEvaluations[t + i].Metric)
            {
                blockEvaluations[t] = blockEvaluations[t + i];
            }
        }
 
        gThread.SyncThreads();
    }
 
    if (gThread.threadIdx.x == 0)
    {
        evaluations[gThread.blockIdx.x] = blockEvaluations[0];
    }
}
好的，所以 TileOrderFromPermutation 函数需要稍作解释。这是一个经典的例子，说明为 GPU 编写代码与为 CPU 编写代码有何不同。在 CPU 上，我们会寻找一种可以顺序迭代排列的算法。在 GPU 上，我们更喜欢一种能够根据排列索引为我们提供排列顺序的算法。尽管这种方法成本更高，但由于 GPU 提供的巨大并行性，它在 GPU 上的速度远远快于 CPU。
© . All rights reserved.