Ускоряем метод Виолы-Джонса (Viola-Jones)

В последнее время метод Виолы-Джонса, который долгое время был основным способом детектирования объектов на изображении, отступает под натиском более новых и совершенных алгоритмов. Тем не менее, актуальность этого метода еще сохраняется и в настоящем времени.

Да, каскадный классификатор основанный на признаках Хаара (метод Виолы-Джонса) уступает в скорости работы каскадному LBP классификатору. Он менее точен, чем детектор, основанный на HOG признаках, и тем более детектор, базирующийся на сверточных нейронных сетях. И все же у него есть определенная ниша, когда требуется точность выше, чем у LBP каскада, но скорость работы более точных детекторов недостаточна высока. Не менее важным фактором является то, что для каскадного Хаар классификатора существует большое количество уже обученных каскадов, в том числе в стандартной поставке библиотеки OpenCV. Поэтому скорость работы этого алгоритма весьма важна. Что и побудило автора в свое время занятся его оптимизацией.


Ну и какая статья об детектировании лиц, может обойтись без фотографии Лены?

Способы повышения производительности

Существует довольно много подробных описаний того, как работает метод Виолы-Джонса, в том числе и на Хабрахабре: источник 1, источник 2. Поэтому я не буду подробно останавливаться на его описании, кому будет интересно прочтут это в выше перечисленных источниках. Хотелось бы сказать пару слов о возможных путях, с помощью которых метод Виолы-Джонса может быть ускорен.

1. Метод Виолы-Джонса не является каким-то жестким фиксированным алгоритмом. Его каскады формируются в результатет обучения. Потому скорость работы может сильно варьироваться в зависимости от сложности объекта и размера объекта, который мы хотим обнаружить. Простейший способ ускорения его работы — это задать ограничение на минимальный объект или увеличить шаг пирамиды. Тем не менее эти способы имеют побочные эффекты в виде уменьшения точности работы или в увеличении числа ложных срабатываний. Переобучение каскадов может иногда дать очень значительный эффект в плане точности и скорости работы работы. Однако это достаточно затратный процесс, связанный со сбором большой обучающей выборки. Да и само обучение не сказать, чтобы очень быстрый процесс.
   
2. Достаточно очевидным способом ускорения, является обработка различных частей изображения в параллельных потоках. Данный метод уже реализован в OpenCV, потому здесь делать ничего не надо.
3. Использовать для параллельной обработки графический ускоритель. Здесь тоже уже есть реализация в OpenCV. Минусом можно считать требование к наличию графического ускорителя.
4. Использовать для ускорения работы алгоритма векторные инструкции процессора SIMD. Векторные инструкции есть практически во всех современных процессорах. К сожалению в разных процессорах они разные — потому приходится иметь разные реализации под разные процессоры. Кроме того, эффективное применение векторных инструкций требует значительной переделки алгоритма. Данному методу и посвящена настоящая статья.

Исходная (скалярная) версия алгоритма

Если взять исходники данного алгоритма из OpenCV, то видны многочисленные попытки его ускорить при помощи различных SIMD (SSE, AVX, NEON) инструкций. У меня нет точных численных оценок, на сколько данные попытки были успешными. Однако, судя по исходным кодам, ускорение вряд ли было значительным, так как или в оптимизациях вообще не используются векторные инструкции (только их скалярные аналоги), или загрузка данных в векторные регистры осуществляется поэлементно, что фатальным образом сказывается на общей производительности алгоритма.

И так, рассмотрим скалярную версию алгоритма каскадного Хаар-классфикатора (далее приводится его упрощенная версия) для заданного масштаба:

void HaarDetect(const HaarCascade & hid, const Image & mask, const Rect & rect, int step, Image & dst)
{
  for (ptrdiff_t row = rect.top; row < rect.bottom; row += step)
  {
    size_t p_offset = row * hid.sum.stride / sizeof(uint32_t);
    size_t pq_offset = row * hid.sqsum.stride / sizeof(uint32_t);
    for (ptrdiff_t col = rect.left; col < rect.right; col += step)
    {
      if (mask.At<uint32_t>(col, row) == 0)
        continue;
      float norm = Norm(hid, pq_offset + col);
      if (Detect(hid, p_offset + col, 0, norm) > 0)
        dst.At<uint32_t>(col, row) = 1;
    }
  }
}

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

uint32_t Sum(uint32_t * const ptr[4], size_t offset)
{
  return ptr[0][offset] - ptr[1][offset] - ptr[2][offset] + ptr[3][offset];
}

float Norm(const HaarCascade & hid, size_t offset)
{
  float sum = float(Sum(hid.p, offset));
  float sqsum = float(Sum(hid.pq, offset));
  float q = sqsum*hid.windowArea - sum *sum;
  return q < 0.0f ? 1.0f : sqrtf(q);
}

А также непосредственно сама классификация:

float WeightedSum(const WeightedRect & rect, size_t offset)
{
  int sum = rect.p0[offset] - rect.p1[offset] - rect.p2[offset] + rect.p3[offset];
  return rect.weight*sum;
}

int Detect(const HaarCascade & hid, size_t offset, int startStage, float norm)
{
  const HaarCascade::Stage * stages = hid.stages.data();
  const HaarCascade::Node * node = hid.nodes.data() + stages[startStage].first;
  const float * leaves = hid.leaves.data() + stages[startStage].first * 2;
  for (int i = startStage, n = (int)hid.stages.size(); i < n; ++i)
  {
    const HaarCascade::Stage & stage = stages[i];
    const HaarCascade::Node * end = node + stage.ntrees;
    float stageSum = 0.0;
    for (; node < end; ++node, leaves += 2)
    {
      const HaarCascade::Feature & feature = hid.features[node->featureIdx];
      float sum = WeightedSum(feature.rect[0], offset) + 
        WeightedSum(feature.rect[1], offset);
      if (feature.rect[2].p0)
        sum += WeightedSum(feature.rect[2], offset);
      stageSum += leaves[sum >= node->threshold*norm];
    }
    if (stageSum < stage.threshold)
      return -i;
  }
  return 1;
}

Здесь признаки Хаара расчитываются при помощи предварительно рассчитанных интегральных изображений.

Способы векторизации алгоритма

Анализ приведенного алгоритма показывает, что основные вычислительные ресурсы тратятся на вычисление взвешенных интегральных сумм (функция WeightedSum). Поэтому их и нужно прежде всего векторизовать. Способов это сделать всего 2:

• Пытаемся параллельно вычислять сразу несколько Хаар-признаков в одном SIMD векторе.
• Пытаемся параллельно обсчитывать сразу несколько точек изображения в одном SIMD векторе.

Векторизация по признакам

На первый взгляд наиболее оптимально производить векторизацию по вычислению Хаар-признаков (к стати, так и сделано в OpenCV). Однако есть одно но: данные, над которыми нужно проводить вычисления, случайным образом разбросаны в памяти. Потому их нужно сначала собирать оттуда при помощи скалярных операций чтения, а затем результаты вычисления раскидывать при помощи скалярных операций записи. Можно конечно воспользоваться gather-scatter операциями (они есть в AVX2 и AVX-512), но скорость их работы сопоставима, а иногда и медленнее, чем если просто использовать скалярный код.

Векторизация по точкам

Векторизация по точкам изображения не так очевидна (фактически мы проводим векторизацию практически по самому внешнему циклу). Однако, данные в интегральном изображении для соседних точек лежат рядом, что позволяет их эффективно загружать при помощи векторных операций. Кроме того метод легко масштабируется для разной длины SIMD вектора. На этом плюсы заканчиваются, и дальше список проблем:

1. Эффективность алгоритма Виолы-Джонса основана на том факте, что на каждой стадии отбрасываются заведомо неправильные точки. Доля отбрасываемых точек зависит от параметров обучения, и по по умолчанию равна 50%. Если мы обрабатываем SIMD-вектор с несколькими точками, то мы будем вынуждены продолжать обсчет стадий до тех пор, пока не отбросятся все его элементы. Это снижает эффективность метода при расчете последующих стадий. К счастью, стадии, на которых отбрасываются соседние точки, достаточно сильно коррелируют. Кроме того, если в векторе осталась только одна точка для проверки, то мы можем перейти к скалярному варианту алгоритма, который в данном случае будет эффективнее.
2. Оригинальный алгоритм из OpenCV в целях оптимизации для небольшого масштаба проверяет точки с шагом 2. Это ускоряет скалярный код, но снижает эффективность векторных инструкций, так как нам уже на первой стадии приходится половину вычислений выполнять впустую. К счастью у данной проблемы есть довольно элегантное решение. Но мы его подробнее распишем чуть позже.

SIMD версия алгоритма

В данной статье я приведу упрощенную версию алгоритма для SSE4.1, версии для AVX2, AVX-512 и NEON используют тот же самый подход.

Для начала рассмотрим случай, когда мы сканируем все точки изображения без пропусков:

void DetectionHaarDetect32fp(const HidHaarCascade & hid, const Image & mask, const Rect & rect, Image & dst)
{
  size_t rightAligned = rect.left + Simd::AlignLo(rect.Width(), 4);
  for (ptrdiff_t row = rect.top; row < rect.bottom; row += 1)
  {
    size_t p_offset = row * hid.sum.stride / sizeof(uint32_t);
    size_t pq_offset = row * hid.sqsum.stride / sizeof(uint32_t);
    size_t col = rect.left;
    for (; col < rightAligned; col += 4)
    {
      __m128i result = _mm_loadu_si128((__m128i*)(mask.Row<uint32_t>(row) + col));
      if (_mm_testz_si128(result, _mm_set1_epi32(1)))
        continue;
      __m128 norm = Norm32fp(hid, pq_offset + col);
      Detect32f(hid, p_offset + col, norm, result);
      _mm_storeu_si128((__m128i*)(dst.Row<uint32_t>(row) + col), result);
    }
    for (; col < rect.right; col += 1)
    {
      if (mask.At<uint32_t>(col, row) == 0)
        continue;
      float norm = Norm32f(hid, pq_offset + col);
      if (Detect(hid, p_offset + col, 0, norm) > 0)
        dst.At<uint32_t>(col, row) = 1;
    }
  }
}

Здесь в отличие от скалярной версии происходит сканирование изображения не по одиночным точкам, а блоками по 4 точки (лишь краевые точки обрабатываются скалярным образом). Нормирующий коэффициент для блока из четырех точек рассчитывается следующим образом:

inline __m128 ValidSqrt(__m128 value)
{
  return _mm_blendv_ps(_mm_set1_ps(1.0f), _mm_sqrt_ps(value), _mm_cmpgt_ps(value, _mm_set1_ps(0.0f)));
}

inline __m128i Sum32ip(uint32_t * const ptr[4], size_t offset)
{
  __m128i s0 = _mm_loadu_si128((__m128i*)(ptr[0] + offset));
  __m128i s1 = _mm_loadu_si128((__m128i*)(ptr[1] + offset));
  __m128i s2 = _mm_loadu_si128((__m128i*)(ptr[2] + offset));
  __m128i s3 = _mm_loadu_si128((__m128i*)(ptr[3] + offset));
  return _mm_sub_epi32(_mm_sub_epi32(s0, s1), _mm_sub_epi32(s2, s3));
}

inline __m128 Norm32fp(const HidHaarCascade & hid, size_t offset)
{
  __m128 area = _mm_set1_ps(hid.windowArea);
  __m128 sum = _mm_cvtepi32_ps(Sum32ip(hid.p, offset));
  __m128 sqsum = _mm_cvtepi32_ps(Sum32ip(hid.pq, offset));
  return ValidSqrt(_mm_sub_ps(_mm_mul_ps(sqsum, area), _mm_mul_ps(sum, sum)));
}

Сама классификация для блока из четырех точек:

inline int ResultCount(__m128i result)
{
  uint32_t SIMD_ALIGNED(16) buffer[4];
  _mm_store_si128((__m128i*)buffer, _mm_sad_epu8(result, _mm_setzero_si128()));
  return buffer[0] + buffer[2];
}

inline __m128 WeightedSum32f(const WeightedRect & rect, size_t offset)
{
  __m128i s0 = _mm_loadu_si128((__m128i*)(rect.p0 + offset));
  __m128i s1 = _mm_loadu_si128((__m128i*)(rect.p1 + offset));
  __m128i s2 = _mm_loadu_si128((__m128i*)(rect.p2 + offset));
  __m128i s3 = _mm_loadu_si128((__m128i*)(rect.p3 + offset));
  __m128i sum = _mm_sub_epi32(_mm_sub_epi32(s0, s1), _mm_sub_epi32(s2, s3));
  return _mm_mul_ps(_mm_cvtepi32_ps(sum), _mm_set1_ps(rect.weight));
}

inline void StageSum32f(const float * leaves, float threshold, const __m128 & sum, const __m128 & norm, __m128 & stageSum)
{
  __m128 mask = _mm_cmplt_ps(sum, _mm_mul_ps(_mm_set1_ps(threshold), norm));
  stageSum = _mm_add_ps(stageSum, _mm_blendv_ps(_mm_set1_ps(leaves[1]), _mm_set1_ps(leaves[0]), mask));
}

void Detect32f(const HidHaarCascade & hid, size_t offset, const __m128 & norm, __m128i & result)
{
  typedef HidHaarCascade Hid;
  const float * leaves = hid.leaves.data();
  const Hid::Node * node = hid.nodes.data();
  const Hid::Stage * stages = hid.stages.data();
  for (int i = 0, n = (int)hid.stages.size(); i < n; ++i)
  {
    const Hid::Stage & stage = stages[i];
    if (stage.canSkip)
      continue;
    const Hid::Node * end = node + stage.ntrees;
    __m128 stageSum = _mm_setzero_ps();
    if (stage.hasThree)
    {
      for (; node < end; ++node, leaves += 2)
      {
        const Hid::Feature & feature = hid.features[node->featureIdx];
        __m128 sum = _mm_add_ps(WeightedSum32f(feature.rect[0], offset), WeightedSum32f(feature.rect[1], offset));
        if (feature.rect[2].p0)
          sum = _mm_add_ps(sum, WeightedSum32f(feature.rect[2], offset));
        StageSum32f(leaves, node->threshold, sum, norm, stageSum);
      }
    }
    else
    {
      for (; node < end; ++node, leaves += 2)
      {
        const Hid::Feature & feature = hid.features[node->featureIdx];
        __m128 sum = _mm_add_ps(WeightedSum32f(feature.rect[0], offset), WeightedSum32f(feature.rect[1], offset));
        StageSum32f(leaves, node->threshold, sum, norm, stageSum);
      }
    }
    result = _mm_andnot_si128(_mm_castps_si128(_mm_cmpgt_ps(_mm_set1_ps(stage.threshold), stageSum)), result);
    int resultCount = ResultCount(result);
    if (resultCount == 0) // Если в векторе не осталось ни одного кандидата, прекращаем поиск:
    {
      return;
    }
    else if (resultCount == 1) // Если в векторе осталась только 1 точка, переходим к скалярной версии алгоритма:
    {
      uint32_t SIMD_ALIGNED(16) _result[4];
      float SIMD_ALIGNED(16) _norm[4];
      _mm_store_si128((__m128i*)_result, result);
      _mm_store_ps(_norm, norm);
      for (int j = 0; j < 4; ++j)
      {
        if (_result[j])
        {
          _result[j] = Detect32f(hid, offset + j, i + 1, _norm[j]) > 0 ? 1 : 0;
          break;
        }
      }
      result = _mm_load_si128((__m128i*)_result);
      return;
    }
  }
}

Как видно, алгоритм практически такой же, как и для скалярной версии (конечно же с поправкой на использование векторов SSE вместо обычных вещественных чисел).

Теперь рассмотрим ситуацию, когда нам надо это же алгоритм выполнить с шагом 2. Самый простой вариант — это грузить из памяти два вектора по 4 числа, а потом из них формировать один (оставляя только четные числа). Однако он показывает совершенно неудовлетвлорительную производительность. К счастью можно сделать небольшой трюк с переупорядочиванием исходных данных — в исходных интегральных изображениях мы четные точки перемещаем в начало строчки, а нечетные в конец. Тогда загрузка векторов с четными/нечетными элементами превращается в тривиальную операцию. Кроме того мы можем использовать для сканирования один и тот же код.

Данные оптимизации легко аналогично реализовать и для другого размера вектора (8 для AVX2 и 16 для AVX-512 ) и для других платформ (например, NEON для платформы ARM).

Результаты оптимизации

Ниже в таблице приведены времена сканирования (в миллисекундах) изображения для разных реализиций алгоритма. Для тестирования использвалось изображение размером 1920x1080. Сканирование производилось без масштабирования исходного изображения:

Function	Scalar	SSE4.1	AVX2	AVX-512
Common	227.499	111.509	74.952	54.579
DetectionHaarDetect32fi(0)	84.792	45.771	32.716	25.961
DetectionHaarDetect32fi(1)	162.779	79.007	50.996	36.841
DetectionHaarDetect32fp(0)	329.904	159.355	109.725	80.862
DetectionHaarDetect32fp(1)	588.270	268.298	172.396	114.735

Здесь цифрами в круглых скобочках (0) и (1) — помечены результаты для разных каскадов, использованных для тестирования. DetectionHaarDetect32fp — функция, сканирующая изображение с шагом 1, DetectionHaarDetect32fi — с шагом 2. Common — среднее геометрическое.

В следующей таблице приведены ускорение, достигнутое относительно базовой (скаляной) версии:

Function	SSE4.1/Scalar	AVX2/Scalar	AVX-512/Scalar
Common	2.04	3.04	4.17
DetectionHaarDetect32fi(0)	1.85	2.59	3.27
DetectionHaarDetect32fi(1)	2.06	3.19	4.42
DetectionHaarDetect32fp(0)	2.07	3.01	4.08
DetectionHaarDetect32fp(1)	2.19	3.41	5.13

Ну и на закуску относительное ускорение, при увеличении размера вектора в два раза:

Function	SSE4.1/Scalar	AVX2/SSE4.1	AVX-512/AVX2
Common	2.04	1.49	1.37
DetectionHaarDetect32fi(0)	1.85	1.40	1.26
DetectionHaarDetect32fi(1)	2.06	1.55	1.38
DetectionHaarDetect32fp(0)	2.07	1.45	1.36
DetectionHaarDetect32fp(1)	2.19	1.56	1.50

Из таблицы наглядно видно убывающую полезность от каждого последующего удвоения вектора.

Программная реализация

В заключение, хотелось бы пару слов сказать о программной реализации вышеуказанных алгоритмов. В рамка проекта Simd был реализован С++ класс Simd::Detection. В нем под капотом спрятана низкоуровневая работа с сишным API билиотеки Simd. Одним из немаловажных достоинств является совместимость используемых каскадов (HAAR и LBP) с каскадами из OpenCV, а также простота использования. Ниже приведен пример детектирования лиц с использованием этого фреймворка:

#include <iostream>
#include <string>

#include "opencv2/opencv.hpp"
#define SIMD_OPENCV_ENABLE
#include "Simd/SimdDetection.hpp"
#include "Simd/SimdDrawing.hpp"

int main(int argc, char * argv[])
{
  if (argc < 2)
  {
    std::cout << "You have to set video source! It can be 0 for camera or video file name." << std::endl;
    return 1;
  }
  std::string source = argv[1];

  cv::VideoCapture capture;
  if (source == "0")
    capture.open(0);
  else
    capture.open(source);
  if (!capture.isOpened())
  {
    std::cout << "Can't capture '" << source << "' !" << std::endl;
    return 1;
  }

  typedef Simd::Detection<Simd::Allocator> Detection;
  Detection detection;
  detection.Load("../../data/cascade/haar_face_0.xml");
  bool inited = false;

  const char * WINDOW_NAME = "FaceDetection";
  cv::namedWindow(WINDOW_NAME, 1);
  for (;;)
  {
    cv::Mat frame;
    capture >> frame;

    Detection::View image = frame;

    if (!inited)
    {
      detection.Init(image.Size(), 1.2, image.Size() / 20);
      inited = true;
    }

    Detection::Objects objects;
    detection.Detect(image, objects);

    for (size_t i = 0; i < objects.size(); ++i)
      Simd::DrawRectangle(image, objects[i].rect, Simd::Pixel::Bgr24(0, 255, 255));

    cv::imshow(WINDOW_NAME, frame);
    if (cvWaitKey(1) == 27)// "press 'Esc' to break video";
      break;
  }
  return 0;
}