并行化霍夫变换算法的最佳方法是什么？

Question

在霍夫线变换中，对于每个边缘像素，我们在霍夫参数空间中找到对应的 Rho 和 Theta。 Rho 和 Theta 的累加器应该是全局的。如果我们想并行化算法，那么分割累加器空间的最佳方法是什么？

Answer 1

并行化算法的最佳方法可能取决于几个方面。一个重要的方面是您所针对的硬件。由于您已用“openmp”标记您的问题，我假设在您的情况下，目标是SMP 系统。

为了回答您的问题，让我们先看看Hough transform 的一个典型、直接的实现（我将使用 C，但以下内容也适用于 C++ 和 Fortran）：

size_t *hough(bool *pixels, size_t w, size_t h, size_t res, size_t *rlimit)
{
  *rlimit = (size_t)(sqrt(w * w + h * h));

  double step = M_PI_2 / res;
  size_t *accum = calloc(res * *rlimit, sizeof(size_t));
  size_t x, y, t;

  for (x = 0; x < w; ++x)
    for (y = 0; y < h; ++y)
      if (pixels[y * w + x])
        for (t = 0; t < res; ++t)
        {
          double theta = t * step;
          size_t r = x * cos(theta) + y * sin(theta);
          ++accum[r * res + t];
        }

  return accum;
}

给定一个黑白像素数组（按行存储）、宽度、高度和霍夫空间角度分量的目标分辨率，函数hough 返回一个累加器数组，用于霍夫空间（按角度组织）并将其距离维度的上限存储在输出参数rlimit 中。也就是说，返回的累加器数组中的元素个数由res * (*rlimit) 给出。

函数体以三个嵌套循环为中心：两个最外层循环遍历输入像素，而有条件执行的最内层循环遍历霍夫空间的角度维度。

为了使算法并行化，我们必须以某种方式将其分解为可以并发执行的片段。通常，这种分解是由计算的结构或被操作的数据的结构引起的。

由于除了迭代之外，函数执行的唯一计算上有趣的任务是最内层循环主体中的三角函数，因此没有明显的基于计算结构的分解机会。因此，让我们专注于基于数据结构的分解，并让我们区分

1. 基于输入数据结构的数据分解，以及
2. 基于输出数据结构的数据分解。

在我们的例子中，输入数据的结构由作为参数传递给函数 hough 的像素数组给出，并由函数体中最外层的两个循环进行迭代。

输出数据的结构由返回的累加器数组的结构给出，并由函数体中最内层的循环迭代。

我们首先将输出数据分解视为，对于霍夫变换，它导致最简单的并行算法。

输出数据分解

将输出数据分解为可以相对独立生成的单元，具体化为让最内层循环的迭代并行执行。

这样做，必须考虑任何所谓的loop-carried dependencies 以使循环并行化。在这种情况下，这很简单，因为没有这种循环携带的依赖关系：循环的所有迭代都需要对共享数组accum 进行读写访问，但每次迭代都在数组的它自己的“私有”段上运行（即，那些具有索引i 和i % res == t 的元素）。

使用 OpenMP，这为我们提供了以下简单的并行实现：

size_t *hough(bool *pixels, size_t w, size_t h, size_t res, size_t *rlimit)
{
  *rlimit = (size_t)(sqrt(w * w + h * h));

  double step = M_PI_2 / res;
  size_t *accum = calloc(res * *rlimit, sizeof(size_t));
  size_t x, y, t;

  for (x = 0; x < w; ++x)
    for (y = 0; y < h; ++y)
      if (pixels[y * w + x])
        #pragma omp parallel for
        for (t = 0; t < res; ++t)
        {
          double theta = t * step;
          size_t r = x * cos(theta) + y * sin(theta);
          ++accum[r * res + t];
        }

  return accum;
}

输入数据分解

可以通过并行化最外层循环来获得遵循输入数据结构的数据分解。

但是，该循环确实具有循环携带的流依赖性，因为每次循环迭代都可能需要对共享累加器数组的每个单元格进行读写访问。因此，为了获得正确的并行实现，我们必须同步这些累加器访问。在这种情况下，这可以通过原子地更新累加器来轻松完成。

循环还带有两个所谓的反依赖。这些是由内部循环的归纳变量y 和t 诱导的，并且通过将它们设为并行外部循环的私有变量来轻松处理。

我们最终得到的并行实现如下所示：

size_t *hough(bool *pixels, size_t w, size_t h, size_t res, size_t *rlimit)
{
  *rlimit = (size_t)(sqrt(w * w + h * h));

  double step = M_PI_2 / res;
  size_t *accum = calloc(res * *rlimit, sizeof(size_t));
  size_t x, y, t;

  #pragma omp parallel for private(y, t)
  for (x = 0; x < w; ++x)
    for (y = 0; y < h; ++y)
      if (pixels[y * w + x])
        for (t = 0; t < res; ++t)
        {
          double theta = t * step;
          size_t r = x * cos(theta) + y * sin(theta);
          #pragma omp atomic
          ++accum[r * res + t];
        }

  return accum;
}

评估

评估这两种数据分解策略，我们观察到：

• 对于这两种策略，我们最终实现了并行化，其中算法的计算量大的部分（三角函数）很好地分布在线程上。
• 分解输出数据为我们提供了函数hough 中最内层循环的并行化。由于这个循环没有任何循环携带的依赖关系，我们不会产生任何数据同步开销。但是，由于对每个设置的输入像素都执行了最内层循环，因此由于重复形成一组线程等，我们确实会产生相当多的开销。
• 分解输入数据可实现最外层循环的并行化。这个循环只执行一次，因此线程开销很小。但是，另一方面，我们确实会产生一些数据同步开销来处理循环承载的流依赖关系。

通常可以假设 OpenMP 中的原子操作非常高效，而众所周知，线程开销相当大。因此，人们期望，对于霍夫变换，输入数据分解提供了一种更有效的并行算法。一个简单的实验证实了这一点。在这个实验中，我将这两个并行实现应用于随机生成的 1024x768 黑白图片，目标分辨率为 90（即每弧度 1 个累加器），并将结果与​​顺序实现进行比较。下表显示了不同线程数的两种并行实现所获得的相对加速：

# threads | OUTPUT DECOMPOSITION | INPUT DECOMPOSITION
----------+----------------------+--------------------
     2    |          1.2         |         1.9         
     4    |          1.4         |         3.7
     8    |          1.5         |         6.8

（该实验是在没有负载的双 2.2 GHz 四核 Intel Xeon E5520 上进行的。所有加速都是五次运行的平均值。顺序实现的平均运行时间为 2.66 秒。）

虚假分享

请注意，并行实现容易受到累加器数组的false sharing 的影响。对于基于输出数据分解的实现，这种错误共享可以通过转置累加器数组（即，通过“按距离”组织它）在很大程度上避免。这样做并衡量影响，在我的实验中，并没有导致任何可观察到的进一步加速。

结论

回到你的问题，“什么是拆分累加器空间的最佳方式？”，答案似乎是最好不要拆分累加器空间，而是拆分输入空间。

如果出于某种原因，您决定拆分累加器空间，您可以考虑更改算法的结构，以便最外层循环遍历霍夫空间，内层循环遍历输入图片的尺寸。这样，您仍然可以派生仅产生一次线程开销并且没有数据同步开销的并行实现。但是，在该方案中，三角函数不再是有条件的，因此总的来说，每次循环迭代都必须比上述方案做更多的工作。