Canny算法中的非极大值抑制：用SSE优化

Question

我有“荣幸”改善以下其他人的代码的运行时间。 （这是来自精明算法的非最大抑制”）。我的第一个想法是使用 SSE 内在代码，我在这方面很新，所以我的问题是。

有没有机会这样做？ 如果是这样，有人可以给我一些提示吗？

void vNonMaximumSupression(
          float* fpDst, 
          float const*const fpMagnitude, 
          unsigned char  const*const ucpGradient,                                                                           ///< [in] 0 -> 0°, 1 -> 45°, 2 -> 90°, 3 -> 135°
int iXCount, 
int iXOffset, 
int iYCount, 
int ignoreX, 
int ignoreY)
{
    memset(fpDst, 0, sizeof(fpDst[0]) * iXCount * iXOffset);

    for (int y = ignoreY; y < iYCount - ignoreY; ++y)
    {
        for (int x = ignoreX; x < iXCount - ignoreX; ++x)
        {
            int idx = iXOffset * y + x;
            unsigned char dir = ucpGradient[idx];
            float fMag = fpMagnitude[idx];

            if (dir == 0 && fpMagnitude[idx - 1]           < fMag && fMag > fpMagnitude[idx + 1] ||
                dir == 1 && fpMagnitude[idx - iXCount + 1] < fMag && fMag > fpMagnitude[idx + iXCount - 1] ||
                dir == 2 && fpMagnitude[idx - iXCount]     < fMag && fMag > fpMagnitude[idx + iXCount] ||
                dir == 3 && fpMagnitude[idx - iXCount - 1] < fMag && fMag > fpMagnitude[idx + iXCount + 1]
                )
                    fpDst[idx] = fMag;
            else
                fpDst[idx] = 0;
        }
    }
}

Answer 1

讨论

正如@harold 所指出的，这里矢量化的主要问题是该算法对每个像素使用不同的偏移量（由方向矩阵指定）。我可以想到几种可能的矢量化方式：

1. @nikie：一次评估所有分支，即将每个像素与其所有邻居进行比较。这些比较的结果会根据方向值进行混合。
2. @PeterCordes：将大量像素加载到 SSE 寄存器中，然后使用 _mm_shuffle_epi8 仅选择给定方向上的邻居。然后执行两个向量化比较。
3. (me)：使用标量指令沿方向加载正确的两个相邻像素。然后将四个像素的这些值组合到 SSE 寄存器中。最后，在 SSE 中进行两次比较。

第二种方法很难有效实施，因为对于一组 4 个像素，有 18 个相邻像素可供选择。我认为这需要太多的洗牌。

第一种方法看起来不错，但它每像素执行的操作要多四倍。我想向量指令的加速会被太多的计算所淹没。

我建议使用第三种方法。您可以在下面看到有关提高性能的提示。

混合方法

首先，我们希望尽可能快地制作标量代码。您提供的代码包含太多分支。其中大多数是不可预测的，例如按方向切换。

为了删除分支，我建议创建一个数组delta = {1, stride - 1, stride, stride + 1}，它给出了方向的索引偏移量。通过使用此数组，您可以找到要与之比较的相邻像素的索引（没有分支）。然后你做两个比较。最后，你可以写一个像res = (isMax ? curr : 0);这样的三元运算符，希望编译器可以为它生成无分支代码。

不幸的是，编译器（至少 MSVC2013）不够聪明，无法避免 isMax 分支。这就是为什么我们可以从使用标量 SSE 内在函数重写内部循环中受益。查找the guide 以供参考。您主要需要以 _ss 结尾的内部函数，因为代码是完全标量的。

最后，我们可以对除了加载相邻像素之外的所有内容进行矢量化。为了加载相邻像素，我们可以使用带有标量参数的_mm_setr_ps 内在函数，要求编译器为我们生成一些好的代码 =)

__m128 forw = _mm_setr_ps(src[idx+0 + offset0], src[idx+1 + offset1], src[idx+2 + offset2], src[idx+3 + offset3]);
__m128 back = _mm_setr_ps(src[idx+0 - offset0], src[idx+1 - offset1], src[idx+2 - offset2], src[idx+3 - offset3]);

我刚刚自己实现了它。在 Ivy Bridge 3.4Ghz 上进行单线程测试。使用 1024 x 1024 分辨率的随机图像作为源。结果（以毫秒为单位）为：

original: 13.078     //your code
branchless: 8.556    //'branchless' code
scalarsse: 2.151     //after rewriting to sse intrinsics
hybrid: 1.159        //partially vectorized code

他们确认了每一步的性能改进。最终的代码需要超过一毫秒来处理一个百万像素的图像。总加速约为11.3倍。确实，您可以在 GPU 上获得更好的性能 =)

我希望提供的信息足以让您重现这些步骤。如果您正在寻找可怕的剧透，请查看 here 我对所有这些阶段的实现。