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

Question

在Hough Line Transform中,对于每个边缘像素,我们在Hough参数空间中找到对应的Rho和Theta.Rho和Theta的累加器应该是全局的.如果我们想要并行化算法,那么拆分累加器空间的最佳方法是什么？

Answer 1

并行化算法的最佳方法可能取决于几个方面.一个重要的方面是您要定位的硬件.正如您使用"openmp"标记了您的问题,我认为,在您的情况下,目标是SMP系统.

为了回答你的问题,让我们首先看一下Hough变换的典型,直接的实现(我将使用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空间的角度分量的宽度,高度和目标分辨率,该函数hough返回Hough空间的累加器数组(组织的"角度") -wise")并在输出参数中存储其距离维度的上限rlimit.也就是说,返回的累加器数组中的元素数由下式给出res * (*rlimit).

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

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

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

1. 数据分解基于输入数据的结构,和
2. 基于输出数据结构的数据分解.

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

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

我们首先看输出数据分解,因为对于Hough变换,它导致最简单的并行算法.

输出数据分解

将输出数据分解成可以相对独立地生成的单元,实现了使最内层循环的迭代并行执行.

这样做,必须考虑到循环并行化的任何所谓的循环携带依赖性.在这种情况下,这是直截了当的,因为没有这样的循环携带依赖:循环的所有迭代都需要对共享数组的读写访问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四核英特尔至强E5520上进行的.所有加速都是五次运行的平均值.顺序实现的平均运行时间为2.66秒.)

虚假分享

注意,并行实现易受累加器阵列的错误共享的影响.对于基于输出数据的分解的实现,通过转置累加器阵列(即,通过将其"按距离"组织),可以在很大程度上避免这种错误共享.在我的实验中,这样做并测量了影响,并没有导致任何可观察到的进一步加速.

结论

回到你的问题,"什么是分割累加器空间的最佳方法？",答案似乎是最好不要分割累加器空间,而是分割输入空间.

如果由于某种原因,你已经决定分割累加器空间,你可以考虑改变算法的结构,以便最外面的循环迭代霍夫空间和内部循环,无论输入图片的最小尺寸是什么.这样,您仍然可以派生一个并行实现,该实现仅产生一次线程开销,并且没有数据同步开销.但是,在该方案中,三角函数不再是有条件的,因此,总的来说,每次循环迭代都必须比上面的方案做更多的工作.