Блог

Я хотел бы выполнить вычисление суммы абсолютных разностей с переменным размером блока с помощью двумерного массива из 16-битных целых чисел в программе на C как можно эффективнее. Меня интересует код сопоставления блоков в реальном времени. Мне было интересно, есть ли какие-либо программные библиотеки, доступные для этого? Код выполняется в Windows XP, и я застрял с использованием Visual Studio 2010 для выполнения компиляции. Процессор представляет собой 2-ядерный AMD Athlon 64 x2 4850e.

Под вычислением суммы абсолютных разностей с переменным размером блока (SAD) я подразумеваю следующее.

У меня есть один двумерный массив меньшего размера, который я буду называть template_grid , и один двумерный массив большего размера, который я буду называть image . Я хочу найти область на изображении, которая минимизирует сумму абсолютных разностей между пикселями в шаблоне и пикселями в области на изображении.

Самый простой способ вычислить SAD в C if был бы следующим:

 for(int shiftY = 0; shiftY < rangeY; shiftY  ) {
    for(int shiftX = 0; shiftX < rangeX; shiftX  ) {
        for(int x = 0; x < lenTemplateX; x  ) {
            for(int y = 0; y < lenTemplateY; y  ) {
                SAD[shiftY][shiftX]=abs(template_grid[x][y] - image[y   shiftY][x   shiftX]);
            }
        }
    }
}
  

Вычисление SAD для определенных размеров массива было оптимизировано в библиотеке Intel performance primitives library. Однако массивы, с которыми я работаю, не соответствуют размерам в этих библиотеках.

Есть два диапазона поиска, с которыми я работаю,

большой диапазон: rangeY = 45, rangeX = 10

небольшой диапазон: rangeY = 4, rangeX = 2

Существует только один размер шаблона, и он равен: lenTemplateY = 61, lenTemplateX = 7

Небольшая оптимизация:

 for(int shiftY = 0; shiftY < rangeY; shiftY  ) {
  for(int shiftX = 0; shiftX < rangeX; shiftX  ) {
    // if you can assume SAD is already filled with 0-es, 
    // you don't need the next line
    SAD[shiftX][shiftY]=0;
    for(int tx = 0, imx=shiftX; x < lenTemplateX; tx  ,imx  ) {
      for(int ty = 0, imy=shiftY; y < lenTemplateY; ty  ,imy  ) {
        // two increments of imx/imy may be cheaper than 
        // two addition with offsets
        SAD[shiftY][shiftX] =abs(template_grid[tx][ty] - image[imx][imy]);
      }
    }
  }
}
  

Развертывание цикла с использованием шаблонов C

Может быть, это сумасшедшая идея для вашей конфигурации (меня беспокоит компилятор C ), но она может сработать. Я не даю никаких гарантий, но попробуйте.

Идея может сработать, потому что ваши template_grid размеры и диапазоны постоянны — таким образом, известны во время компиляции.
Кроме того, чтобы это сработало, ваши image и template_grid должны быть организованы с одинаковым расположением (сначала столбец или сначала строка) — способ, которым ваш «пример кода» изображен в вопросе, смешивает SAD x/y с template_grid y/x .
В последующих разделах я предполагаю организацию «сначала по столбцу», так что SAD[ix] это ix ^{обозначает} -й SAD** столбец вашей в вашей матрице. Код выполняется точно так же для «row first», за исключением того, что имена переменных не будут соответствовать значению ваших массивов значений.

Итак, давайте начнем:

 template <
  typename sad_type, typename val_type,
  size_t template_len
> struct sad1D_simple {
  void operator()(
    const val_type* img, const val_type* templ,
    sad_typeamp; result
  ) {
    // template specialization recursion, with one less element to add
    sad1D_simple<sad_type, val_type, template_len-1> one_shorter;
    // call it incrementing the img and template offsets
    one_shorter(img 1, templ 1, result);
    // the add the contribution of the first diff we skipped over above
    result =abs(*(img template_len-1)-*(templ template_len-1));
  }
};

// at len of 0, the result is zero. We need it to stop the
template <
  typename sad_type, typename val_type
>
struct sad1D_simple<sad_type, val_type, 0> {
  void operator()(
    const val_type* img, const val_type* templ,
    sad_typeamp; result
  ) {
    result=0;
  }
};
  

Почему структура функтора — struct with operator? C не допускает частичной специализации шаблонов функций.
Что sad1D_simple делает: запускает for цикл, который вычисляет SAD два массива во входных данных без какого-либо смещения, основываясь на том факте, что длина вашего template_grid массива является константой, известной во время компиляции. Это в том же духе, что и «вычисление факториала времени компиляции с использованием шаблонов C «

Как это помогает?
Пример использования в приведенном ниже коде:

 typedef ulong SAD_t;
typedef int16_t pixel_val_t;

const size_t lenTemplateX = 7; // number of cols in the template_grid
const size_t lenTemplateY = 61;
const size_t rangeX=10, rangeY=45;

pixel_val_t **image, **template_grid;
SAD_t** SAD;
// assume those are initialized somehow


for(size_t tgrid_col=0; tgrid_col<lenTemplateX; tgrid_col  ) {
  pixel_val_t* template_col=template_grid[tgrid_col];
  // the X axis - horizontal - is the column axis, right?
  for(size_t shiftX=0; shiftX < rangeX; shiftX  ) {
    pixel_val_t* img_col=image[shiftX];
    for(size_t shiftY = 0; shiftY < rangeY; shiftY  ) {
      // the Y axis - vertical - is the "offset in a column"=row, isn't it?
      pixel_val_t* img_col_offsetted=img_col shiftY;

      // this functor is made by recursive specialization
      // there's no cycle inside it, it was unrolled into
      // lenTemplateY individual subtractions, abs-es and additions 
      sad1D_simple<SAD_t, pixel_val_t, lenTemplateY> calc;
      calc(img_col_offsetted, template_col, SAD[shiftX][shiftY]);
    }
  }
}
  

Мммм … можем ли мы сделать лучше? Нет, это не будет разворачивание по оси X, мы по-прежнему хотим оставаться в одномерной области, но… что ж, может быть, если мы создадим ranged sad1D и развернем еще один цикл по той же оси?
Это сработает, еслиf rangeX также является постоянным.

 template <
  typename sad_type, typename val_type,
  size_t range, size_t template_len
> struct sad1D_ranged {
  void operator()(
    const val_type* img, const val_type* templ,
    // result is assumed to have at least `range` slots
    sad_type* result
  ) {
    // we'll compute here the first slot of the result
    sad1D_simple<sad_type, val_type, template_len>
      calculator_for_first_sad;
    calculator_for_first_sad(img, templ, *(result));

    // now, ask for a recursive specialization for 
    // the next (range-1) sad-s
    sad1D_ranged<sad_type, val_type, range-1, template_len>
       one_less_in_range;
    // when calling, pass the shifted img and result
    one_less_in_range(img 1, templ, result 1);
  }
};

// for a range of 0, there's nothing to do, but we need it
// to stop the template specialization recursion
template <
  typename sad_type, typename val_type,
  size_t template_len
> struct sad1D_ranged<sad_type, val_type, 0, template_len> {
  void operator()(
    const val_type* img, const val_type* templ,
    // result is assumed to have at least `range` slots
    sad_type* result
  ) {
  }
};
  

И вот как вы это используете:

 for(size_t tgrid_col=0; tgrid_col<lenTemplateX; tgrid_col  ) {
  pixel_val_t* template_col=template_grid[tgrid_col];
  for(size_t shiftX=0; shiftX < rangeX; shiftX  ) {
    pixel_val_t* img_col=image[shiftX];
    SAD_t* sad_col=SAD[shiftX];

    sad1D_ranged<SAD_t, pixel_val_t, rangeY, lenTemplateY> calc;
    calc(img_col, template_col, sad_col);
  }
}
  

ДА… но вопрос в том, улучшит ли это производительность?
Черт возьми, если я знаю. Для небольшого количества циклов в цикле и для сильной локализации данных (значения близки друг к другу, так что они находятся в кэшах процессора), развертывание цикла должно улучшить производительность. При большем количестве циклов вы можете негативно повлиять на прогнозирование ветвления процессора и прочую чушь, которая, как я знаю, может повлиять на производительность, но я не знаю, как.

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

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

вы могли бы попробовать сопоставление шаблона OpenCV с параметром square difference, см. Руководство здесь. OpenCV оптимизирован с помощью OpenCL, но я не знаю для этой конкретной функции. Я думаю, вам стоит попробовать.

Я не уверен, насколько вы ограничены в использовании SAD, или если вы вообще заинтересованы в поиске области на изображении, которая наилучшим образом соответствует шаблону. В последнем случае вы можете использовать свертку вместо SAD. Это может быть решено в области Фурье за O (N log N), включая преобразование Фурье (FFT).

Короче говоря, вы можете использовать БПФ (например, с помощью http://www.fftw.org /) для преобразования шаблона и изображения в частотную область, затем их умножения и обратного преобразования во временную область.

Конечно, все это не имеет значения, если вы привязаны к использованию SAD.