Каков наилучший способ распараллеливания алгоритма преобразования Хафа?

Лучший способ распараллеливания алгоритма может зависеть от нескольких аспектов. Важным таким аспектом является оборудование, на которое вы нацелены. Поскольку вы отметили свой вопрос "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 возвращает массив-накопитель для пространства Хоу (организованного "под углом") и сохраняет верхнюю границу для его измерения расстояния в выходном аргументе rlimit. То есть число элементов в возвращаемом массиве накопителей задается через res * (*rlimit).

Тело функции сосредоточено на трех вложенных циклы: два крайних цикла повторяются по входным пикселям, в то время как условно выполняемый внутренний цикл повторяется по угловому измерению пространства Хафа.

Чтобы распараллелить алгоритм, мы должны каким-то образом разложить его на части, которые могут выполняться одновременно. Обычно такая декомпозиция индуцируется либо структурой вычислений , либо структурой данных , которые оперируются.

Поскольку, помимо итерации, единственной вычислительно интересной задачей, которую выполняет функция, является тригонометрия в теле самой внутренней петли, при этом нет очевидных возможностей для декомпозиции на основе структуры вычислений. Поэтому давайте сосредоточимся на разложениях, основанных на структуре данных, и будем различать

1. декомпозиции данных, основанные на структуре входных данных, и
2. декомпозиции данных, основанные на структуре выходные данные.

Структура входных данных в нашем случае задается массивом пикселей, который передается в качестве аргумента функции 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 обычно можно считать довольно эффективными, в то время как накладные расходы на потоковую обработку, как известно, довольно велики. Следовательно, можно ожидать, что для преобразования Хоу декомпозиция входных данных дает более эффективный параллельный алгоритм. Это подтверждается простым экспериментом. Для этого экспериментируя, я применил две параллельные реализации к случайно сгенерированному черно-белому изображению 1024х768 с целевым разрешением 90 (то есть 1 аккумулятор на градус дуги) и сравнил результаты с последовательной реализацией. В этой таблице показаны относительные ускорения, полученные двумя параллельными реализациями для различных чисел потоков:

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

(эксперимент проводился на двухъядерном процессоре Intel Xeon E5520 с частотой 2,2 ГГц без нагрузки. Все результаты являются в среднем за пять пробегов. Среднее время выполнения последовательной реализации составило 2,66 с.)

Оно

Обратите внимание, что параллельные реализации подвержены ложному совместному использованию массива накопителей. Для реализации, основанной на декомпозиции выходных данных, этого ложного совместного использования можно в значительной степени избежать путем транспонирования массива накопителей (т. е. путем организации его "по расстоянию"). Делая это и измеряя воздействие, сделал, по-моему эксперименты, не приводящие к каким-либо наблюдаемым дальнейшим ускорениям.

Заключение

Возвращаясь к вашему вопросу: "как лучше всего разделить пространство аккумулятора?", ответ, по-видимому, заключается в том, что лучше вообще не разделять пространство аккумулятора, а вместо этого разделить пространство ввода.

Если по какой-то причине вы решили разбить пространство накопителя, вы можете рассмотреть возможность изменения структуры алгоритма таким образом, чтобы самые внешние циклы повторялись над Пространство Хафа и внутренняя петля над тем, что является наименьшим из размеров входного изображения. Таким образом, вы все еще можете получить параллельную реализацию, которая выполняет потоковую нагрузку только один раз и которая освобождается от накладных расходов на синхронизацию данных. Однако в этой схеме тригонометрия больше не может быть условной, и поэтому в общей сложности каждая итерация цикла должна будет выполнять больше работы, чем в приведенной выше схеме.