MISTを用いた実践的なプログラム例その２として、テンプレートマッチングによる物体の検出を行います。ここでは位置と大きさが不明な物体を入力画像中から精度良く検出するために、参照画像の位置と入力画像の大きさを変化させながら照合を行います。

output

検出結果イメージ画像 (実際の出力は、入力画像中の検出領域と参照画像が同じ大きさになるように、入力画像を拡縮しています)

以下に、SSDAとアクティブ探索というアルゴリズムを用いた高速検出手法を紹介します。

SSDA (Sequential Similarity Detection Algorithm)とは

画素間の相違度を用いてマッチングを行う検出手法です。これまでに求められた最小の相違度を暫定値として、ある位置のテンプレートの相違度を計算し、これが明らかに暫定値を超える場合は途中で計算を打ち切るという手法です。【参考文献：D. I. Barnea and H. F. Silverman: A class of algorithms for fast digital image registration, IEEE Trans. Computers，21，pp.179-186，1972】 #アクティブ探索とは

色ヒストグラム間類似度を用いて、テンプレートを大きくスキップさせながらマッチングを行う検出手法です。テンプレートをある位置に置いたときの類似度から、その近傍の類似度の上限値を計算します。上限値が探索値よりも小さければ、その領域での探索が省略できるため、照合回数を極端に低減できます。【参考文献：村瀬洋, V.V.Vinod: 局所色情報を用いた高速物体探索---アクティブ探索法---, 電子情報通信学会論文誌，Vol.J81-DII, No.9，pp.2035-2042，1998】 #手法の比較

各手法の比較実験として、テンプレートをラスタスキャンして検出を行う関数を用意します。用いる手法は以下の4つです。

SSDA
アクティブ探索
ピクセル輝度差を用いたラスタスキャン
ヒストグラムインターセクションを用いたラスタスキャン

SSDAとアクティブ探索は、MISTを用いて以下のように書けます。各関数には、入力画像input、参照画像referenceを渡します。 x、y、scaleには、検出結果（入力画像のスケールとそのx座標、y座標）が入ります。

各手法とも、まず入力画像を参照画像と同じくらいの大きさに縮小した後、定数SCALE_FACTORだけ拡大しながら探索を行います。

1.定義

はじめに，プログラム内で定義した定数を紹介します．

#define Q 16                　//量子化を16階調で行う
#define SCALE_FACTOR 1.25   //入力画像の拡大率を設定
#define MAX_SCALE 5.0       //最大の拡大率を設定

//RGB画素型の設定
typedef mist::rgb < unsigned char > pixel_type;
typedef mist::array2 < pixel_type > image_type;

2.SSDA

SSDAでは，画像間相違度の尺度として，画素間の差異を使用します．以下に，それらを計算する関数を用意します．

diff_SDA() : 画素間の差の絶対値（SDA:Sum of Absolute Difference）
diff_SSD() : 差の二乗（SSD:Sum of Squared Difference）

double diff_SAD( const pixel_type &p1,
                 const pixel_type &p2 )
{
    double ret = 0;
    ret += std::abs( static_cast< int >( p1.r ) - p2.r );
    ret += std::abs( static_cast< int >( p1.g ) - p2.g );
    ret += std::abs( static_cast< int >( p1.b ) - p2.b );
    return ret;
}

inline double _square( const int v )
{
    return v * v;
}

double diff_SSD( const pixel_type &p1,
                 const pixel_type &p2 )
{
    double ret = 0;
    ret += _square( static_cast< int >( p1.r ) - p2.r );
    ret += _square( static_cast< int >( p1.g ) - p2.g );
    ret += _square( static_cast< int >( p1.b ) - p2.b );
    return ret;
}

diff_SAD()またはdiff_SSD()を使ったSSDAは以下のようになります．

void ssda_search( const image_type  &input,
                  const image_type  &reference,
                  size_t            &ret_x,
                  size_t            &ret_y,
                  double            &ret_scale )
{   
    //入力画像を参照画像のサイズにするときの縦横を計算（縦横比は固定で、縦横の大きい方に合わせる）
    double scale;
    if( input.width( ) * reference.height( ) <  reference.width( ) * input.height( ) )
    {
        scale = static_cast< double >( reference.width( ) ) / input.width( );
    }
    else
    {
        scale = static_cast< double >( reference.height( ) ) / input.height( );
    }

    //暫定の最小相違度、位置
    double min_diff = 255 * 255 * 3;
    double diff;
    size_t width, height, size;
    image_type ref;

    //参照画像の大きさを変えて照合
    while( scale < MAX_SCALE ) //拡大率の上限を指定
    {
        std::cout << scale << " 倍)\n";
        std::cout << "暫定の最小相違度:" << min_diff << "\n";

        width = static_cast< size_t >( reference.width( ) / scale );
        height = static_cast< size_t >( reference.height( ) / scale );
        size = width * height;
        
        //参照画像のサイズを変更
        mist::linear::interpolate( reference, ref, width, height );
        
        min_diff *= size;
        //参照画像を1画素ずつずらして照合
        for( size_t j = 0 ; j < input.height( ) - height + 1 ; j ++ )
        {
            for( size_t i = 0 ; i < input.width( ) - width + 1 ; i ++ )
            {
                //重なった画像内で1画素ずつ差異を計算
                diff = 0;
                for( size_t l = 0 ; l < height ; l ++ )
                {
                    for( size_t k = 0 ; k < width ; k ++ )
                    {
                        diff += diff_SSD( input( i + k , j + l ) , ref( k , l ) );
                        if( diff > min_diff )
                        {
                            goto SSDA_loop_exit;
                        }
                    }
                }
                //差異が最小のとき、返す座標を更新
                min_diff = diff;
                ret_scale = scale;
                ret_x = static_cast< unsigned int >( i * ret_scale );
                ret_y = static_cast< unsigned int >( j * ret_scale );
SSDA_loop_exit:;
            }
        }
        min_diff /= size;

        //スケールを変更する
        scale *= SCALE_FACTOR;
    }
    std::cout << width << " " << height << " " << min_diff << std::endl;
}

3.アクティブ探索

アクティブ探索は，上記の通り，ヒストグラムインターセクションを類似度基準とした探索手法です．アクティブ探索本体の説明の前に，それに用いる関数を説明します．下の4つです．

(1). 階調数を落とす (2). ヒストグラムを作成する (3). ヒストグラムインターセクションを計算する (4). 周囲の上限値を計算し，類似度の計算不要な画素をパスする

(1)階調数を落とす．

img配列を，tone階調へ落とした結果の配列をoutputに格納します．但しimg配列に格納した画像は 256階調であることを仮定しています．

void tone_down( const image_type &in, 
        image_type &out, 
        const unsigned char base )
{
    out.resize( in.width( ), in.height() );

    for( size_t i = 0 ; i < out.size( ) ; i ++ )
    {
        out[ i ].r  = in[ i ].r / base;
        out[ i ].g  = in[ i ].g / base;
        out[ i ].b  = in[ i ].b / base;
    }
}

(2)ヒストグラムを作成する．

img配列中の，横方向の位置pi，縦方向の位置pjを窓の左上のピクセルとして，サイズw_width × w_heightの窓を対象領域としたヒストグラムを histogramへ格納します．但し，histogramは，この関数へ渡す前にリサイズを行う必要があります．

void mk_histogram( const image_type &img, 
           const size_t pi, const size_t pj, 
           const size_t w_width, const size_t w_height, 
           mist::array3< double > &histogram )
{
    histogram.fill( 0.0 );

    for( size_t j = pj ; j < pj + w_height ; j ++ )
    {
        for(size_t i = pi ; i < pi + w_width ; i ++ )
        {
          histogram( static_cast< unsigned int >( img( i, j ).r ), static_cast< unsigned int >( img( i, j ).g ), static_cast< unsigned int >( img( i, j ).b) ) ++;
        }
    }
    const size_t w_size = w_width * w_height;
    for( size_t i = 0 ; i < histogram.size( ) ; i ++ )
    {
        histogram[ i ] /= w_size;
    }
}

(3)ヒストグラムインターセクションを計算する．

2つのヒストグラム，histogram_aとhistogram_bのヒストグラムインターセクションを計算します．

double histogram_intersection( const mist::array3<double> &histogram_a, 
                   const mist::array3<double> &histogram_b)
{
    double intersection = 0;

    for( size_t i = 0 ; i < histogram_a.size( ) ; i ++ )
    {
        intersection += ( histogram_a[ i ] < histogram_b[ i ] ) ? histogram_a[ i ] : histogram_b[ i ];
    }
    return intersection;
}

(4). 周囲の上限値を計算し，類似度の計算不要な画素をパスする

今回は，アクティブ探索時に上限値を超えない，ヒストグラムインターセクションの計算が不要なピクセルの位置情報を，bool型の2次元配列culling_tableを用いて管理します．この配列は入力画像と同じサイズの配列で，culling_tableがfalseの位置のピクセルでは，ヒストグラムインターセクションの計算が不要であることを示しています．

inline double _minimum( const double v0, const double v1 )
{
    return ( v0 < v1 ) ? v0 : v1;
}

inline double _upper_bound( const size_t a_and_b, const size_t window_size, const double &similarity )
{
    return ( _minimum( similarity * window_size, a_and_b ) + window_size - a_and_b ) / window_size;
}

void _culling_out( const size_t pi, const size_t pj, const size_t w_width, const size_t w_height, 
           const double &sim, const double &max_sim, 
           mist::array2< bool > &culling_table )
{
    const size_t w_size = w_width * w_height;
    for( size_t j = 0; j < w_height ; j ++ )
    {
        if( _upper_bound( w_width * ( w_height - j ), w_size, sim ) >= max_sim )
        {
            break;
        }
        for( size_t i = 0 ; i < w_width ; i ++ )
        {
            //上限値が目標類似値より小さい場合，ヒストグラムインターセクションを計算する必要が無いことをマークする
            if( _upper_bound( ( w_width - i ) * ( w_height - j ), w_size, sim ) < max_sim )
            {
                culling_table( pi + i, pj + j ) = true;
                if( i < pi )
                {
                    culling_table( pi - i, pj + j ) = true;
                }
            }
            else
            {
                break;
            }
        }
    }
}

以上の定義・関数を利用して，アクティブ探索を実装します．

void active_search( const image_type &input, 
            const image_type &reference, 
            unsigned int &ret_x, unsigned int &ret_y, 
            double &ret_scale )
{

    //階調数を落とす
    image_type degraded_input, degraded_reference;
    tone_down( input, degraded_input, 256 / Q );
    tone_down( reference, degraded_reference, 256 / Q );

    //データを初期化する
    mist::array3< double > input_histogram( Q, Q, Q );
    mist::array3< double > reference_histogram( Q, Q, Q ); 

    //参照画像のヒストグラムを作成する
    mk_histogram( degraded_reference, 0, 0, reference.width( ), reference.height( ), reference_histogram );

    //アスペクト比を変えずに，画像サイズを変更
    double scale;
    if( input.width( ) * reference.height( ) <  reference.width( ) * input.height( ) )
    {
        scale = static_cast< double >( reference.width( ) ) / input.width( );
    }
    else
    {
        scale = static_cast< double >( reference.height( ) ) / input.height( );
    }
    
    size_t window_width = static_cast< size_t >( reference.width( ) / scale );
    size_t window_height = static_cast< size_t >( reference.height( ) / scale );
    size_t window_size = window_width * window_height;

    mist::array2< bool > culling_table( input.width( ), input.height( ) );

    double max_similarity = 0; 
    double similarity;

    //設定した倍率になるまで探索を行う
    while( scale < MAX_SCALE )
    {
        culling_table.fill( false );

        //image_type image( window_width, window_height );

        for( size_t j = 0 ; j < input.height( ) - window_height + 1 ; j ++ )
        {
            for( size_t i = 0 ; i < input.width( ) - window_width + 1 ; i ++ )
            {

                //未処理のピクセルである
                if( !culling_table( i, j ) )
                {
                    //ヒストグラムインターセクションを計算する
                    mk_histogram( degraded_input, i, j, window_width, window_height, input_histogram );
                    similarity = histogram_intersection( input_histogram, reference_histogram );


                    //最大類似度を更新する
                    if( similarity > max_similarity )
                    {
                        max_similarity = similarity;
                        ret_scale = scale;
                        ret_x = static_cast< unsigned int >( i * ret_scale );
                        ret_y = static_cast< unsigned int >( j * ret_scale );
                    }
                    
                    //周囲のピクセルを探索する
                    _culling_out( i, j, window_width, window_height, similarity, max_similarity, culling_table );
                }
            }
        }
        
        //スケールを変更する
        scale *= SCALE_FACTOR;
        window_width = static_cast< size_t >( reference.width( ) / scale );
        window_height = static_cast< size_t >( reference.height( ) / scale );
        window_size = window_width * window_height;
    }
}

4.探索スピードの比較

紹介した4つの探索方法の比較をします．プログラム本体はこちら探索方法の切り替えは61行目からの4つの関数を切り替えます．画像は，入力画像としてこれ，参照画像としてはこれを，それぞれ使用しました．使用した計算機は下の2台です．

	コンパイラ	CPU	メモリ
マシンA	Visual C++	Pentium IV 2.8GHz	1GB
マシンB	GCC3.4.4	Athlon Dual 2600+	2GB

それぞれのマシンでの計算時間．マシンAについて

SSDA	30.29 (sec)
アクティブ探索	7.98 (sec)
ピクセル輝度差	34.9 (sec)
ヒストグラムインターセクション	92.30 (sec)

マシンBについて

SSDA	34.47 (sec)
アクティブ探索	04.54 (sec)
ピクセル輝度差	44.43 (sec)
ヒストグラムインターセクション	146.26 (sec)

この結果からも，高速化が有効に働いていることが確認できます．どちらのマシンにおいても，アクティブ探索が最も高速に探索を終了しています． SSDAと比較しても，アクティブ探索が高速な理由は，類似度計算に，サイズに対して不変量である正規化ヒストグラムを用いているため，線形補間での画像の変形が不要であることが考えられます．

物体検出の高速化 - mist-team/mist GitHub Wiki

SSDA (Sequential Similarity Detection Algorithm)とは

1.定義

2.SSDA

3.アクティブ探索

(1)階調数を落とす．

(2)ヒストグラムを作成する．

(3)ヒストグラムインターセクションを計算する．

(4). 周囲の上限値を計算し，類似度の計算不要な画素をパスする

4.探索スピードの比較

⚠️ GitHub.com Fallback ⚠️

物体検出の高速化 - mist-team/mist GitHub Wiki

SSDA (Sequential Similarity Detection Algorithm)とは

1.定義

2.SSDA

3.アクティブ探索

(1)階調数を落とす．

(2)ヒストグラムを作成する．

(3)ヒストグラムインターセクションを計算する．

(4). 周囲の上限値を計算し，類似度の計算不要な画素をパスする

4.探索スピードの比較

⚠️ **GitHub.com Fallback** ⚠️

⚠️ GitHub.com Fallback ⚠️