FC2ブログ

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

SSE によるテンプレートマッチングの高速化

テンプレートマッチングは,ある画像中から特定のパターンに類似した領域を探すポピュラーな方法です.原理がシンプルで使いやすい手法ですが,画像全体にわたって反復的に計算を行うため,計算量は膨大になります.計算速度を向上するには,計算を並列化することが有効であり,SSE (Streaming SIMD Extensions) などの利用が考えられます.SSE は,CPU の 128 bit レジスタを用いて 32 bit 単精度浮動小数点演算を 4 並列で行うことのできる拡張命令です.浮動小数点演算を大量に行う処理で高速化が期待できます.本記事では,SSE によるテンプレートマッチングの高速化を行います.

Keywords: テンプレートマッチング, SSE


テンプレートとの相似度を評価する指標には,SSD (Sum of Squared Difference) や SAD (Sum of Absolute Difference) ,正規化相互相関 (NCC, Normalized Cross-Correlation) といったものがあり,検出したい特徴や取得される画像の条件に応じて使い分けます.本記事では,ゼロ平均正規化相互相関によるテンプレートマッチングを実装します.

ゼロ平均正規化相互相関 (Zero-mean Normalized Cross-Correlation) は以下の式で表されます.

simd-template_match-eq1.png (1)

計算効率を向上するため,本記事では式(2)を用いて相関係数を計算します.式(2)はこちらの記事を参考にしました.

simd-template_match-eq2.png (2)


通常の浮動小数点演算と,SSE をそれぞれ利用して,式(2)によるテンプレートマッチングを実装しました.ソースコードはそれぞれ template_match.c,template_match_sse.c となります.
#ifndef TEMPLATE_MATCH_H__
#define TEMPLATE_MATCH_H__

typedef struct _BitmapData
{
int dim0;
int dim1;

int stride0;
int stride1;

float *data;

} BitmapData;

int create_bitmap_data(int dim0, int dim1, BitmapData **dst);

void dispose_bitmap_data(BitmapData *obj);

//! Calculate similarity using template matching technique.
/*!
@param[in] src Source image of template matching. Image size must
be larger than size of template.
@param[in] tmp Template image illustrating characteristics of your
target.
@param[out] dst Destination buffer for storing evaluated similarity.

@return true/false indicates whether process has been succeeded or
not.
*/
int template_match(BitmapData *src, BitmapData *tmp, BitmapData **dst);

#endif /* TEMPLATE_MATCH_H__ */
#include "template_match.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <math.h>

#define MYAPP_TRUE 1
#define MYAPP_FALSE 0

#ifdef _MSC_VER
typedef __int32 int32_t;
typedef unsigned __int32 uint32_t;
#else
#include <stdint.h>
#endif

#ifdef _MSC_VER
#define MYAPP_INLINE __inline
#else
#define MYAPP_INLINE inline
#endif

#ifndef INFINITY
#define INFINITY myapp_inf()

static MYAPP_INLINE float myapp_inf(void)
{
// Bit field
// 0 | 111 1111 1 | 000 0000 -- 0000 0000
// sign | exp | fraction
// 31 | 30 23 | 22 0
static uint32_t value = 0x7F800000;

return *((float*)&value);
}

#endif /* INFINITY */

typedef struct _NCCRegion
{
/// @{
int dim0;
int dim1;
/// @}

int scanline;

float ptp;

float scaled_mean;

float scaled_length;

float *data;

} NCCRegion;

static int
create_ncc_region(int dim0, int dim1, NCCRegion **dst)
{
NCCRegion *p = NULL;

if(!(p = malloc(sizeof(NCCRegion))) ||
!(p->data = malloc(sizeof(float) * dim0 * dim1)))
{
free(p);
return MYAPP_FALSE;
}

memset(p->data, 0, sizeof(float) * dim0 * dim1);

p->dim0 = dim0;
p->dim1 = dim1;
p->scanline = dim1;

p->scaled_mean = p->scaled_length = 0;

*dst = p;

return MYAPP_TRUE;
}

static void dispose_ncc_region(NCCRegion *p)
{
if(p != NULL) {
free(p->data);
free(p);
}
}

static int
get_ncc_region(BitmapData *src, int r0, int r1, NCCRegion *dst)
{
int i, j;
float pmin, pmax, *p, *q, *origin;

if(r0 < 0 || r0 + dst->dim0 > src->dim0 ||
r1 < 0 || r1 + dst->dim1 > src->dim1)
{
// out of bounds error
return MYAPP_FALSE;
}

origin = src->data + r0 * src->stride0 +
r1 * src->stride1;

pmin = INFINITY;
pmax = -INFINITY;

for(i=0; i<dst->dim0; i++) {
p = origin + i * src->stride0;
q = dst->data + i * dst->scanline;

for(j=0; j<dst->dim1; j++) {
if(*p < pmin) pmin = *p;
if(*p > pmax) pmax = *p;

*q = *p;

p += src->stride1;
q += 1;
}
}

dst->ptp = pmax - pmin;

return MYAPP_TRUE;
}

static void build_template(NCCRegion *win)
{
int i, j;
float sum, sumsq, *p;

sum = sumsq = 0;

for(i=0; i<win->dim0; i++) {
p = win->data + i * win->scanline;

for(j=0; j<win->dim1; j++) {
sum += (*p);
sumsq += (*p) * (*p);

p += 1;
}
}

win->scaled_mean = sum;
win->scaled_length = sumsq * win->dim0 * win->dim1 - sum * sum;
}

static float cross_correlation(NCCRegion *src, NCCRegion *tmp)
{
int i, j;
float a0, a1, sum, sumsq, dot, *p, *q;

if(src->ptp == 0) {
// not enough characteristics
return 0;
}

sum = sumsq = dot = 0;

for(i=0; i<src->dim0; i++) {
p = src->data + i * src->scanline;
q = tmp->data + i * tmp->scanline;

for(j=0; j<src->dim1; j++) {
sum += (*p);
sumsq += (*p) * (*p);

dot += (*p) * (*q);

p += 1;
q += 1;
}
}

src->scaled_mean = sum;
src->scaled_length = sumsq * src->dim0 * src->dim1 - sum * sum;

a0 = dot * src->dim0 * src->dim1 - src->scaled_mean * tmp->scaled_mean;
a1 = sqrtf(src->scaled_length * tmp->scaled_length);

return a0 / a1;
}

int template_match(BitmapData *src, BitmapData *tmp, BitmapData **dst)
{
int i, j, imax, jmax, ret;
BitmapData *out;
NCCRegion *ncc_tmp, *ncc_src;

ret = MYAPP_FALSE;

imax = src->dim0 - tmp->dim0 + 1;
jmax = src->dim1 - tmp->dim1 + 1;

if(imax <= 0 || jmax <= 0) {
// Template is larger than source image, exceptional parameter
return MYAPP_FALSE;
}

ncc_src = ncc_tmp = NULL; out = NULL;

if(!create_ncc_region(tmp->dim0, tmp->dim1, &ncc_src) ||
!create_ncc_region(tmp->dim0, tmp->dim1, &ncc_tmp) ||
!create_bitmap_data(imax, jmax, &out))
{
// memory exception
goto func_error;
}

if(!get_ncc_region(tmp, 0, 0, ncc_tmp)) {
// improbable case, fatal error
goto func_error;
}

if(ncc_tmp->ptp > 0) {
// Template has enough characteristics

build_template(ncc_tmp);

for(i=0; i<out->dim0; i++) {
for(j=0; j<out->dim1; j++) {
if(!get_ncc_region(src, i, j, ncc_src)) {
// improbable case, return immediately
goto func_error;
}

out->data[i * out->stride0 +
j * out->stride1]
= cross_correlation(ncc_src, ncc_tmp);
}
}
}

*dst = out; ret = MYAPP_TRUE;

goto func_final;

func_error:
dispose_bitmap_data(out);

func_final:
dispose_ncc_region(ncc_src);
dispose_ncc_region(ncc_tmp);
return MYAPP_TRUE;
}

int create_bitmap_data(int dim0, int dim1, BitmapData **dst)
{
BitmapData *obj = NULL;

if(dim0 < 0 || dim1 < 0) { // parameter check
return MYAPP_FALSE;
}

if(!(obj = malloc(sizeof(BitmapData))) ||
!(obj->data = malloc(sizeof(float) * dim0 * dim1)))
{
free(obj);
return MYAPP_FALSE;
}

memset(obj->data, 0, sizeof(float) * dim0 * dim1);

obj->dim0 = dim0;
obj->dim1 = dim1;

obj->stride0 = dim1;
obj->stride1 = 1;

*dst = obj;

return MYAPP_TRUE;
}

void dispose_bitmap_data(BitmapData *obj)
{
if(obj != NULL) {
free(obj->data);
free(obj);
}
}
#include "template_match.h"

#include <xmmintrin.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <math.h>

#define MYAPP_TRUE 1
#define MYAPP_FALSE 0

#ifdef _MSC_VER
typedef __int32 int32_t;
typedef unsigned __int32 uint32_t;
#else
#include <stdint.h>
#endif

#ifdef _MSC_VER
#define MYAPP_INLINE __inline
#else
#define MYAPP_INLINE inline
#endif

#ifndef INFINITY
#define INFINITY myapp_inf()

static MYAPP_INLINE float myapp_inf(void)
{
// Bit field
// 0 | 111 1111 1 | 000 0000 -- 0000 0000
// sign | exp | fraction
// 31 | 30 23 | 22 0
static uint32_t value = 0x7F800000;

return *((float*)&value);
}

#endif /* INFINITY */

#ifdef _MSC_VER
#define MEM_ALIGN(x) __declspec(align(x))
#else
#define MEM_ALIGN(x) __attribute__((aligned(x)))
#endif

static MYAPP_INLINE void* malloc_aligned(size_t size, size_t alignment)
{
#ifdef _MSC_VER
return _aligned_malloc(size, alignment);
#else
void *p;
return posix_memalign(&p, alignment, size) == 0 ? p : NULL;
#endif
}

static MYAPP_INLINE void free_aligned(void *p)
{
#ifdef _MSC_VER
_aligned_free(p);
#else
free(p);
#endif
}

typedef struct _NCCRegion
{
/// @{
int dim0;
int dim1;
/// @}

int scanline;

float ptp;

float scaled_mean;

float scaled_length;

float *data;

} NCCRegion;

static int
create_ncc_region(int dim0, int dim1, NCCRegion **dst)
{
int scanline;

NCCRegion *p = NULL;

if(dim1 % 4 == 0) {
scanline = dim1;
} else {
scanline = dim1 + 4 - dim1 % 4;
}

if(!(p = malloc(sizeof(NCCRegion))) ||
!(p->data = malloc_aligned(sizeof(float) * dim0 * scanline, 16)))
{
free(p);
return MYAPP_FALSE;
}

memset(p->data, 0, sizeof(float) * dim0 * scanline);

p->dim0 = dim0;
p->dim1 = dim1;
p->scanline = scanline;

p->scaled_mean = p->scaled_length = 0;

*dst = p;

return MYAPP_TRUE;
}

static void dispose_ncc_region(NCCRegion *p)
{
if(p != NULL) {
free_aligned(p->data);
free(p);
}
}

static int
get_ncc_region(BitmapData *src, int r0, int r1, NCCRegion *dst)
{
int i, j;
float pmin, pmax, *p, *q, *origin;

__m128 pmin_, pmax_, pix_;

MEM_ALIGN(16) float r[4], s[4];

if(r0 < 0 || r0 + dst->dim0 > src->dim0 ||
r1 < 0 || r1 + dst->dim1 > src->dim1)
{
// out of bounds error
return MYAPP_FALSE;
}

origin = src->data + r0 * src->stride0 +
r1 * src->stride1;

for(i=0; i<dst->dim0; i++) {
p = origin + i * src->stride0;

q = dst->data + i * dst->scanline;

for(j=0; j<dst->dim1; j++) {
*q = *p;

p += src->stride1;
q += 1;
}
}

pmin = INFINITY; pmin_ = _mm_set1_ps( INFINITY);
pmax = -INFINITY; pmax_ = _mm_set1_ps(-INFINITY);

for(i=0; i<dst->dim0; i++) {
p = dst->data + i * dst->scanline;

for(j=0; j < dst->dim1/4; j++) {
pix_ = _mm_load_ps(p);

pmin_ = _mm_min_ps(pmin_, pix_);
pmax_ = _mm_max_ps(pmax_, pix_);

p += 4;
}

for(j=0; j < dst->dim1%4; j++) {
if(*p < pmin) pmin = *p;
if(*p > pmax) pmax = *p;

p += 1;
}
}

_mm_store_ps(r, pmin_);
_mm_store_ps(s, pmax_);

for(i=0; i<4; i++) {
if(r[i] < pmin) pmin = r[i];
if(s[i] > pmax) pmax = s[i];
}

dst->ptp = pmax - pmin;

return MYAPP_TRUE;
}

static void build_template(NCCRegion *win)
{
int i, j;
float sum, sumsq, *p;

__m128 sum_, sumsq_, pix_;

MEM_ALIGN(16) float r[4], s[4];

sum_ = sumsq_ = _mm_setzero_ps();

for(i=0; i<win->dim0; i++) {
p = win->data + i * win->scanline;

for(j=0; j < win->scanline/4; j++) {
pix_ = _mm_load_ps(p);

sum_ = _mm_add_ps(sum_, pix_);
sumsq_ = _mm_add_ps(sumsq_, _mm_mul_ps(pix_, pix_));

p += 4;
}
}

_mm_store_ps(r, sum_);
_mm_store_ps(s, sumsq_);

sum = r[0] + r[1] + r[2] + r[3];
sumsq = s[0] + s[1] + s[2] + s[3];

win->scaled_mean = sum;
win->scaled_length = sumsq * win->dim0 * win->dim1 - sum * sum;
}

static float cross_correlation(NCCRegion *src, NCCRegion *tmp)
{
int i, j;
float a0, a1, sum, sumsq, dot, *p, *q;

__m128 sum_, sumsq_, dot_, p_, q_;

MEM_ALIGN(16) float r[4], s[4];

if(src->ptp == 0) {
// not enough characteristics
return 0;
}

sum_ = sumsq_ = dot_ = _mm_setzero_ps();

for(i=0; i<src->dim0; i++) {
p = src->data + i * src->scanline;
q = tmp->data + i * tmp->scanline;

for(j=0; j < src->scanline/4; j++) {
p_ = _mm_load_ps(p);

sum_ = _mm_add_ps(sum_, p_);
sumsq_ = _mm_add_ps(sumsq_, _mm_mul_ps(p_, p_));

q_ = _mm_load_ps(q);

dot_ = _mm_add_ps(dot_, _mm_mul_ps(p_, q_));

p += 4;
q += 4;
}
}

_mm_store_ps(r, sum_);
_mm_store_ps(s, sumsq_);

sum = r[0] + r[1] + r[2] + r[3];
sumsq = s[0] + s[1] + s[2] + s[3];

src->scaled_mean = sum;
src->scaled_length = sumsq * src->dim0 * src->dim1 - sum * sum;

_mm_store_ps(r, dot_);

dot = r[0] + r[1] + r[2] + r[3];

a0 = dot * src->dim0 * src->dim1 - src->scaled_mean * tmp->scaled_mean;
a1 = sqrtf(src->scaled_length * tmp->scaled_length);

return a0 / a1;
}

int template_match(BitmapData *src, BitmapData *tmp, BitmapData **dst)
{
int i, j, imax, jmax, ret;
BitmapData *out;
NCCRegion *ncc_tmp, *ncc_src;

ret = MYAPP_FALSE;

imax = src->dim0 - tmp->dim0 + 1;
jmax = src->dim1 - tmp->dim1 + 1;

if(imax <= 0 || jmax <= 0) {
// Template is larger than source image, exceptional parameter
return MYAPP_FALSE;
}

ncc_src = ncc_tmp = NULL; out = NULL;

if(!create_ncc_region(tmp->dim0, tmp->dim1, &ncc_src) ||
!create_ncc_region(tmp->dim0, tmp->dim1, &ncc_tmp) ||
!create_bitmap_data(imax, jmax, &out))
{
// memory exception
goto func_error;
}

if(!get_ncc_region(tmp, 0, 0, ncc_tmp)) {
// improbable case, fatal error
goto func_error;
}

if(ncc_tmp->ptp > 0) {
// Template has enough characteristics

build_template(ncc_tmp);

for(i=0; i<out->dim0; i++) {
for(j=0; j<out->dim1; j++) {
if(!get_ncc_region(src, i, j, ncc_src)) {
// improbable case, return immediately
goto func_error;
}

out->data[i * out->stride0 +
j * out->stride1]
= cross_correlation(ncc_src, ncc_tmp);
}
}
}

*dst = out; ret = MYAPP_TRUE;

goto func_final;

func_error:
dispose_bitmap_data(out);

func_final:
dispose_ncc_region(ncc_src);
dispose_ncc_region(ncc_tmp);
return ret;
}

int create_bitmap_data(int dim0, int dim1, BitmapData **dst)
{
BitmapData *obj = NULL;

if(dim0 < 0 || dim1 < 0) { // parameter check
return MYAPP_FALSE;
}

if(!(obj = malloc(sizeof(BitmapData))) ||
!(obj->data = malloc(sizeof(float) * dim0 * dim1)))
{
free(obj);
return MYAPP_FALSE;
}

memset(obj->data, 0, sizeof(float) * dim0 * dim1);

obj->dim0 = dim0;
obj->dim1 = dim1;

obj->stride0 = dim1;
obj->stride1 = 1;

*dst = obj;

return MYAPP_TRUE;
}

void dispose_bitmap_data(BitmapData *obj)
{
if(obj != NULL) {
free(obj->data);
free(obj);
}
}
これらを用いて SSE の使用による計算精度への影響と,高速化の効果を検討します.精度の評価では,既知の座標にテンプレート画像と等しい部分領域を有する画像を入力として,テンプレートマッチングを行います.相関値が 1 となれば正解です.速度の評価は,ランダム画像を入力としてテンプレートマッチングを行い,消費時間を計測します.
計測条件は,入力画像サイズを 1000x1000,テンプレートサイズを 32x32 とします.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <time.h> // For performance test

#include "template_match.h"

#define SRC_DIM0 1000
#define SRC_DIM1 1000

#define TMP_DIM0 32
#define TMP_DIM1 32

#define SIGMA 2.0

int main(int argc, char **argv)
{
int i, j, k;
int pos[][2] = {
{273, 293}, {934, 862}, {780, 366}, {516, 948}, {234, 523},
{854, 932}, {747, 757}, { 50, 337}, {480, 527}, { 92, 608}
};
BitmapData *tmp, *src, *dst0, *dst1;

clock_t t0, t1;

tmp = src = dst0 = dst1 = NULL;

if(!create_bitmap_data(TMP_DIM0, TMP_DIM1, &tmp)) {
// memory error
exit(1);
}

for(i=0; i<TMP_DIM0; i++) {
for(j=0; j<TMP_DIM1; j++) {
double x, y;

x = (double)j - (double)(TMP_DIM1 + 1) / 2;
y = (double)i - (double)(TMP_DIM0 + 1) / 2;

tmp->data[i * tmp->stride0 + j]
= (float)exp(-0.5 * (x*x + y*y) / (SIGMA * SIGMA));
}
}

if(!create_bitmap_data(SRC_DIM0, SRC_DIM1, &src)) {
// memory error
exit(1);
}

for(k=0; k<10; k++) {
int r0, r1;

r0 = pos[k][0];
r1 = pos[k][1];

for(i=0; i<TMP_DIM0; i++) {
for(j=0; j<TMP_DIM1; j++) {
float *p, *q;

p = src->data + (r0 + i) * src->stride0 +
(r1 + j) * src->stride1;
q = tmp->data + i * tmp->stride0 +
j * tmp->stride1;

*p = *q;
}
}
}

template_match(src, tmp, &dst0);

for(i=0; i<SRC_DIM0; i++) {
for(j=0; j<SRC_DIM1; j++) {
src->data[i * src->stride0 +
j * src->stride1]
= (float)rand() / 32768.f;
}
}

t0 = clock(); template_match(src, tmp, &dst1); t1 = clock();

for(k=0; k<10; k++) {
int r0, r1;

r0 = pos[k][0];
r1 = pos[k][1];

printf("%f\n", dst0->data[r0 * dst0->stride0 + r1 * dst0->stride1]);
}

printf("Elapsed time: %.2f\n", (double)(t1 - t0) / CLOCKS_PER_SEC);

dispose_bitmap_data(src);
dispose_bitmap_data(tmp);
dispose_bitmap_data(dst0);
dispose_bitmap_data(dst1);

return 0;
}
私の実行環境 (2.5GHz Mac mini) での出力は,それぞれ以下の通りです.
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
Elapsed time: 2.33
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
Elapsed time: 1.33
以上の検討から,SSE を使用することにより約二倍の高速化が見込まれ,本記事での画像およびテンプレートの条件では計算精度も十分得られることが分かりました.

最後に,ビルドに使用した CMakeLists.txt を掲載します.
cmake_minimum_required (VERSION 2.6)
project (TEMPLATE_MATCHING_SAMPLE)

if (${CMAKE_SIZEOF_VOID_P} EQUAL 4)
if (MSVC)
set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /arch:SSE2")
else ()
set (CMAKE_C_FLAGS ${CMAKE_C_FLAGS} "-msse2 -mfpmath=sse")
endif ()
endif ()

set (CMAKE_BUILD_TYPE release)

add_executable (tm main.c template_match.c)
add_executable (tm_sse main.c template_match_sse.c)

if (UNIX)
target_link_libraries(tm m)
target_link_libraries(tm_sse m)
endif ()
スポンサーサイト

テーマ : プログラミング
ジャンル : コンピュータ

コメントの投稿

非公開コメント

プロフィール

Ishida Akihiko

Author:Ishida Akihiko
FC2ブログへようこそ!

免責事項
当サイトに掲載する記事内容は,必ずしも正確性,信頼性,妥当性,有用性,完成度などを保証しません.記事の利用はすべて自己責任でお願いします.当サイトに掲載された内容によって発生したいかなる損害に対しても,管理人は一切の責任を負いかねます.
最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
アクセスカウンター
検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

QRコード
QR
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。