画像処理

(6) Scale-Invariant Feature Transform ( SIFT )

画像認識においては、入力した画像がそのまま利用されることは少なく、その特徴を抽出したデータ ( 特徴量 ) への変換が行われます。特徴量は、例えば "りんご" を表すときに、それがどんな画像の中の "りんご" であっても不変であることが求められます。しかし、"りんご" はあらゆる方向を向いているし、大きさもバラバラでしょう。さらに明るさやノイズなどの因子などによっても特徴量が影響を受けることが考えられます。
特徴量が大きさや向き、明るさやノイズなどの因子に対して不変であることは、画像認識において非常に重要です。そのような特徴量の代表として、今回は「Scale-Invariant Feature Transform ( SIFT )」を紹介します。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) スケール・スペース ( Scale Space )

以前、「エッジ抽出」を行う手法の一つである「LoG フィルタ ( Laplacian of Gaussian Filter )」を紹介しました。このフィルタは、二階微分 ( ラプラス演算子 ) とガウシアン・フィルタを組み合わせたもので、ガウシアン・フィルタによって誤判定の原因となるノイズが除去されるという効果があります。
LoG フィルタには別の効果もあります。矩形波に対して、様々な標準偏差のガウシアン・フィルタと LoG フィルタを掛けた時の波形の変化を以下に示します。

図 1-1. 矩形波に対するガウシアン・フィルタと LoG フィルタ
ガウシアン・フィルタなし ( 二階微分のみ )
矩形波LoG フィルタ
矩形波矩形波のLoGフィルタ
σ = 1
ガウシアン・フィルタLoG フィルタ
σ=1のガウシアン・フィルタσ=1のLoGフィルタ
σ = 3
ガウシアン・フィルタLoG フィルタ
σ=3のガウシアン・フィルタσ=3のLoGフィルタ
σ = 5
ガウシアン・フィルタLoG フィルタ
σ=5のガウシアン・フィルタσ=5のLoGフィルタ

左側のグラフはガウシアン・フィルタのみを行った結果、右側のグラフはさらに二階微分を行った結果です。また、最も上にあるグラフはガウシアン・フィルタを掛けず、二階微分だけを行ったもの ( つまり左側は処理前の矩形波 ) で、下になるほどガウシアン・フィルタの標準偏差 ( σ ) が大きくなっていきます。離散データを処理した結果なので、滑らかな曲線にはなっていませんが、標準偏差が 3 以上になるとかなり滑らかに見えるようになります。
ここで注目したいのが σ = 3 や σ = 5 での LoG フィルタの処理結果で、どちらも矩形波の中央付近で極値となっています。エッジ部分の極値が、ガウシアン・フィルタによってぼやけた形になり、両端のエッジが重ね合わされることで極値が生じるわけです。ある小領域が、その周囲の領域とは状況が異なっている場合、その領域を「ブロブ ( Blob )」と呼びます。LoG フィルタは、ブロブを検出するのにも利用できることをこの結果は示しています。

しかし、ブロブの形状や大きさによって、極値が得られる σ の値は変わります。そこで、複数の σ で LoG フィルタを掛けて、その中から極値を検出します。このとき、二次元の画像 ( x, y ) に σ を加えた ( x, y, σ ) の三次元の空間から極値を検出するイメージになります。これを「スケール・スペース ( Scale Space )」といいます。
スケール・スペースを使って極値を検出することによって、任意の大きさのブロブを検出することが可能になります。例えば、元の画像が 2 倍に拡大されていても、同じ領域を検出することができるわけです。このように、スケールに対して不変な ( Scale invariant ) 検出器であるというのが、LoG フィルタを用いたスケール・スペースの大きな特徴になります。

ところが、以前説明した通り、LoG フィルタは処理が重いという欠点があります。そのため、SIFT では LoG フィルタとほぼ同等の結果が得られる「DoG フィルタ ( Difference of Gaussian Filter )」を採用しています。


SIFT では、初期値 σ から始めて定数 k 倍しながら画像をガウシアン・フィルタに掛けていきます。こうして σ, kσ, k2σ ... とフィルタを掛けた結果 G( x, y σ ), G( x, y, kσ ) ... の差分 ( G( x, y, kσ ) - G( x, y, σ ), G( x, y, k2σ ) - G( x, y, kσ ) ... ) から DoG フィルタの結果が得られます。

σ = ( σ12 + σ22 )1/2 でフィルタを掛ける場合、

f( x, y ) * g( x, y ; σ ) = f( x, y ) * g( x, y ; σ1 ) * g( x, y ; σ2 )

が成り立ちます ( 補足 1 )。つまり、原画像に σ, kσ ... と畳み込み積分する代わりに

σn = ( knσ2 - kn-1σ2 )1/2

を使って前回の処理結果を畳み込み積分しても同じ結果が得られます。よって、フィルタリングするために常に原画像を残しておく必要はなく、前回の結果さえあれば処理を継続することができます。

さらに処理を高速化するため一工夫します。2σ のガウシアン・フィルタで処理された画像を 1 / 2 に縮小した画像を、1 / 2 に縮小してから σ のガウシアン・フィルタで処理された画像とみなします。こうすることによって、σ の値をあまり大きくすることなくガウシアン・フィルタで処理することができます。σ が大きくなると、フィルタのサイズも大きくしなければならないので、処理は重くなっていきます。σ を大きくする代わりに画像を縮小することで、処理時間は大幅に短縮することができます。

このため通常は k = 21/s として σ, 21/sσ, 22/sσ ... と順に処理を行い、2σ になったら 1 / 2 に縮小して、再度 21/sσ, 22/sσ ... と順番に処理します ( もちろん、実際は前の処理結果を 2n/sσ - 2(n-1)/sσ でフィルタリングします )。画像を縮小するまでの一連の処理をオクターブ ( Octave ) といいます。2σ までフィルタリングをするためには、オクターブあたり s 回の処理を行えばいいのですが、後述する極値抽出には 3 枚の DoG フィルタ結果が必要になるため、合計で s + 2 枚の DoG フィルタが必要になり、さらに DoG フィルタには 2 枚のガウシアン・フィルタ処理結果が必要なので、実際には s + 3 枚のガウシアン・フィルタ結果を生成することになります。なお、s はオクターブあたりのスケール数を表します。

図 1-2. オクターブ ( s = 2, k = √2 の場合 )
オクターブ

オクターブの回数は、画像の大きさで判定します。縮小した画像が一定の大きさよりも小さくなったら、処理を終了します。


スケール・スペースを生成するためのサンプル・プログラムを以下に示します。

/*
  CreateGaussianFilter : ガウシアン・フィルタの生成

  フィルタのサイズは、配列 filter のサイズで決定する。

  sigma : σの値
  filter : 生成したフィルタを保持する配列へのポインタ
*/
void CreateGaussianFilter( double sigma, vector< double >* filter )
{
  assert( sigma > 0 );
  assert( filter != 0 );

  // フィルタのサイズ(片側のサイズとする)
  int sz = filter->size();
  if ( sz < 1 ) return;
  sz = ( sz - 1 ) / 2;

  for ( int i = -sz ; i <= sz ; ++i ) {
    double norm = std::pow( static_cast< double >( i ), 2.0 );
    double s2_2 = 2 * std::pow( sigma, 2.0 );
    (*filter)[i + sz] = std::exp( -norm / s2_2 ) / ( std::sqrt( M_PI * s2_2 ) );
  }
}

/*
  Convolusion : フィルタによる畳込み積分

  src : フィルタを掛ける対象のデータ
  c : フィルタを掛ける対象の位置
  filter : フィルタ
  isHorizontal : 水平方向の畳み込み積分なら True

  戻り値 : 畳み込み積分の結果
*/
double Convolusion( const Matrix< double >& src, const Coord< int >& c,
                    const vector< double >& filter, bool isHorizontal )
{
  double ans = 0;

  int sz = filter.size();
  if ( sz < 1 ) return( 0 );
  sz = ( sz - 1 ) / 2;

  for ( int i = -sz ; i <= sz ; ++i ) {
    if ( isHorizontal ) {
      ans += src[c.y][Move( c.x, i, src.cols() )] * filter[sz - i];
    } else {
      ans += src[Move( c.y, i, src.rows() )][c.x] * filter[sz - i];
    }
  }

  return( ans );
}

/*
  Convolusion_2D : 2 次元方向の畳み込み積分

  src : 対象のデータ
  dest : 処理結果を保持するデータへのポインタ
  filterX,filterY X/Y 方向のフィルタ
*/
void Convolusion_2D( const Matrix< double >& src, Matrix< double >* dest,
                     const vector< double >& filterX, const vector< double >& filterY )
{
  assert( dest != 0 );

  // 水平方向のガウシアン・フィルタ結果を保持するバッファ
  Matrix< double > buffer( src.rows(), src.cols() );

  dest->assign( src.rows(), src.cols() );

  Coord< int > c;

  // 水平方向のフィルタ
  for ( c.y = 0 ; c.y < static_cast< int >( src.rows() ) ; ++( c.y ) )
    for ( c.x = 0 ; c.x < static_cast< int >( src.cols() ) ; ++( c.x ) )
      buffer[c.y][c.x] = Convolusion( src, c, filterX, true );

  // 垂直方向のフィルタ
  for ( c.y = 0 ; c.y < static_cast< int >( src.rows() ) ; ++( c.y ) ) {
    for ( c.x = 0 ; c.x < static_cast< int >( src.cols() ) ; ++( c.x ) ) {
      (*dest)[c.y][c.x] = Convolusion( buffer, c, filterY, false );
    }
  }
}

/*
  Diff : 二つのデータの差分を計算する

  data1, data2 : 二つのデータ
  res : 結果を保持するデータへのポインタ
*/
void Diff( const Matrix< double >& data1, const Matrix< double >& data2, Matrix< double >* res )
{
  assert( res != 0 );
  assert( data1.rows() == data2.rows() && data1.cols() == data2.cols() );

  res->assign( data1.rows(), data1.cols() );
  for ( Matrix< double >::size_type r = 0 ; r < res->rows() ; ++r ) {
    for ( Matrix< double >::size_type c = 0 ; c < res->cols() ; ++c ) {
      (*res)[r][c] = data1[r][c] - data2[r][c];
    }
  }
}

/*
  DownSampling : データの縦・横を半分にダウンサンプリングする

  src : サンプリング対象データ
  dest : 結果を保持する変数へのポインタ
*/
void DownSampling( const Matrix< double >& src, Matrix< double >* dest )
{
  assert( dest != 0 );

  dest->assign( src.rows() / 2, src.cols() / 2 );
  for ( Matrix< double >::size_type r = 0 ; r < dest->rows() ; ++r )
    for ( Matrix< double >::size_type c = 0 ; c < dest->cols() ; ++c )
      // 4 ピクセルの平均を採用する
      (*dest)[r][c] = ( src[r * 2][c * 2] + src[r * 2][c * 2 + 1] +
                        src[r * 2 + 1][c * 2] + src[r * 2 + 1][c * 2 + 1] ) / 4;
}

/*
  CreateDoG : DoG 画像の作成

  src : 対象の画像データ
  dog : 作成した DoG 画像を保持する配列へのポインタ
  gauss : 作成した平滑化画像を保持する配列へのポインタ
  s0 : ガウシアンフィルタの標準偏差の初期値
  minSize : 処理を行う画像データの最小サイズ
*/
void CreateDoG( Matrix< double > src, vector< Matrix< double > >* dog, vector< Matrix< double > >* gauss,
                double s0, size_t minSize )
{
  const size_t GAUSS_COUNT = KeyPointValue::scalePerOctave() + 3; // ガウシアンフィルタを行う回数
  const size_t DOG_COUNT = KeyPointValue::scalePerOctave() + 2;   // DoG画像の枚数

  // ガウシアンフィルターの作成
  // 再生性を利用して前のレイヤにσの差分でフィルタを掛ける
  vector< double > filter[GAUSS_COUNT];
  double s1 = s0; // 次のσは s0 で初期化
  s0 = 0;         // 最初のσはゼロとする
  for ( size_t i = 0 ; i < GAUSS_COUNT ; ++i ) {
    double sigma = std::sqrt( s1 * s1 - s0 * s0 );
    filter[i].resize( static_cast< size_t >( std::round( sigma * 6 ) ) | 1 ); // フィルタのサイズは 6σ とする
    CreateGaussianFilter( sigma, &filter[i] );
    s0 = s1;
    s1 *= std::pow( 2.0, 1.0 / KeyPointValue::scalePerOctave() );
  }

  dog->clear();
  gauss->clear();
  gauss->push_back( Matrix< double >() );
  Convolusion_2D( src, &( gauss->back() ), filter[0], filter[0] );

  while ( ( gauss->back() ).rows() >= minSize && ( gauss->back() ).cols() >= minSize ) {
    // ガウシアンフィルターによる畳み込み積分
    for ( size_t k = 1 ; k < GAUSS_COUNT - 1 ; ++k ) {
      gauss->push_back( Matrix< double >() );
      Convolusion_2D( *( gauss->rbegin() + 1 ), &( gauss->back() ), filter[k], filter[k] );
    }
    Matrix< double > gaussBack;
    Convolusion_2D( gauss->back(), &gaussBack, filter[GAUSS_COUNT - 1], filter[GAUSS_COUNT - 1] );
    // DoG 画像の作成
    for ( size_t k = gauss->size() - DOG_COUNT + 1 ; k < gauss->size() ; ++k ) {
      dog->push_back( Matrix< double >() );
      Diff( (*gauss)[k], (*gauss)[k - 1], &( dog->back() ) );
    }
    dog->push_back( Matrix< double >() );
    Diff( gaussBack, gauss->back(), &( dog->back() ) );
    // ダウンサンプリング
    gauss->push_back( Matrix< double >() );
    DownSampling( (*gauss)[gauss->size() - DOG_COUNT + KeyPointValue::scalePerOctave() - 1], &( gauss->back() ) );
  }
  gauss->pop_back();
}

サンプル・プログラムの中で使用されている Matrix クラスは二次元配列と同等であり、vector< vector< double > > などと置き換えて考えれば理解しやすいと思います。また、KeyPointValue というクラスの静的メンバ関数 scalePerOctave が使われていますが、これはオクターブあたりのスケール数 ( 極値計算の結果の枚数 ) を返す関数です。KeyPointValue については後ほど紹介します。

最後に示した CreateDoG が、DoG 画像を生成するプログラムのメイン・ルーチンになります。最初に、ガウシアン・フィルタに使う配列を生成していますが、この生成に利用する関数が CreateGaussianFilter です。前の処理結果を利用してフィルタリングするために、σ に差分を使っているのに注意して下さい。
次に、原画像 src を σ でガウシアン・フィルタに掛けます。ここで利用する関数が Convolusion_2D で、二次元方向に畳み込み積分を行うための関数です。この関数は、一次元方向用の畳み込み積分関数 Convolusion を縦・横方向に 2 回呼び出しているだけの単純なものです。通常ならば、二次元方向に生成したカーネルを使って畳み込み積分を行う必要がありますが、ガウシアン・フィルタの場合

G( x, y ; σ ) = [ 1 / ( 2πσ2 )1/2 ] exp( -x2 / 2σ2 )・[ 1 / ( 2πσ2 )1/2 ] exp( -y2 ) / 2σ2 ) = G( x ; σ )G( y ; σ )

より縦・横を別々に処理しても同じ結果が得られます。

DoG フィルタ、ガウシアン・フィルタのいずれの結果も後で利用するため、引数として渡した配列に保持するようになっています。DoG フィルタの結果は極値計算のためにすべて必要になりますが、ガウシアン・フィルタの結果は不要なものもあります。このプログラムでは、σ が一致するように両方とも同数の結果を保持するようにして、添字の管理を簡略化しています。しかし、生成される枚数はガウシアン・フィルタの結果の方が一枚多いので、その一枚分は gaussBack という変数で保持するようにしています。必要な枚数分のガウシアン・フィルタ処理結果が得られたら、あとは DoG を生成するために関数 Diff を使って差分を計算し、最後に DownSampling で 2σ でガウシアン・フィルタ処理した結果を半分のサイズに縮小します。これを、サイズが一定量より小さくなるまで繰り返します。


2) キーポイントの生成と絞り込み

DoG フィルタの結果が得られたら、次に画像の特徴を求める箇所 ( キーポイント ) を抽出します。最初に、DoG フィルタ処理結果を使って極値を検索します。極値は連続した 3 枚の DoG フィルタ処理結果を使い、中央の DoG フィルタ処理結果の一点に着目したとき、同一スケールの縦・横・斜め方向にある近隣ピクセル 8 点と、上下のスケールの近隣ピクセル 9 x 2 点の計 26 点と比較して、最大または最小となっていたら、それをキーポイントの候補点とします。

図 2-1. 極値の検出
極値の検出

DoG フィルタの結果は、オクターブあたり s + 2 枚あります。両端の極値は抽出できないので、2 枚目から s + 1 枚目までの s 枚の結果を使って極値を抽出することになります。この処理を全オクターブについて行い、キーポイントの候補点とします。なお、1 オクターブ目の極値の抽出は 21/sσ から始まり、2σ で終わります。2 オクターブ目の 1 枚目は 2σ の処理結果を表すので、極値の抽出は 2(s+1)/sσ から始まることになります。このように、各オクターブについてスケールは連続していることに注意してください。


極値の抽出を行うためのサンプル・プログラムを以下に示します。

/**
   KeyPointValue : キーポイント用クラス
**/
class KeyPointValue
{
  static size_t scalePerOctave_; // オクターブごとのスケール数

  size_t dogIndex_;    // DoG 画像の添字
  Coord< int > dogP_; // DoG 画像に合わせた座標値

  Coord< double > subP_; // サブピクセルの座標
  double subScale_;      // サブピクセルのスケール

public:

  /*
    DoG 画像上の座標とDoG 画像への添字を指定して構築

    dogP : DoG 画像上の座標
    dogIndex : DoG 画像への添字
  */
  KeyPointValue( const Coord< int >& dogP, size_t dogIndex = 0 )
    : dogIndex_( dogIndex ), dogP_( dogP ), subP_( 0, 0 ), subScale_( 0 ) {}

  /*
    setScalePerOctave : オクターブあたりのスケール数をセットする

    scalePerOctave : オクターブあたりのスケール数
  */
  static void setScalePerOctave( size_t scalePerOctave )
  { scalePerOctave_ = scalePerOctave; }

  /*
    scalePerOctave : オクターブあたりのスケール数を返す
  */
  static size_t scalePerOctave() { return( scalePerOctave_ ); }

  /*
    octave : オクターブ番号を返す
  */
  int octave() const { return( dogIndex_ / ( scalePerOctave_ + 2 ) ); }

  /*
    scale : スケール番号を返す
  */
  int scale() const { return( octave() * scalePerOctave_ + ( dogIndex_ % ( scalePerOctave_ + 2 ) ) ); }

  /*
    ratio : 原画像に対する DoG 画像の大きさの比率(逆数)を返す
  */
  int ratio() const { return( std::pow( 2, octave() ) ); }

  /*
    dogIndex : DoG 画像への添字を返す
  */
  size_t dogIndex() const { return( dogIndex_ ); }

  /*
    p : 原画像上に合わせた座標値を返す
  */
  Coord< int > p() const { return( Coord< int >( dogP_.x * ratio() + ratio() / 2, dogP_.y * ratio() + ratio() / 2 ) ); }

  /*
    dogP : DoG 画像上の座標値を返す
  */
  Coord< int > dogP() const { return( dogP_ ); }

  /*
    setDogP : DoG 画像上の座標値をセットする

    dogP : DoG 画像上の座標値
  */
  void setDogP( const Coord< int >& dogP ) { dogP_ = dogP; }

  /*
    setDogIndex : DoG 画像への添字をセットする

    dogIndex : DoG 画像への添字
  */
  void setDogIndex( size_t dogIndex ) { dogIndex_ = dogIndex; }

  /*
    subP : サブピクセルの座標値を返す
  */
  Coord< double > subP() const { return( subP_ ); }

  /*
    subScale : サブピクセルのスケール値を返す
  */
  double subScale() const { return( subScale_ ); }

  /*
    setSubPixel : サブピクセルの値をセットする

    p : サブピクセルの座標値
    scale : サブピクセルのスケール値
  */
  void setSubPixel( const Coord< double >& p, double scale )
  {
    subP_ = p;
    subScale_ = scale;
  }
};

/*
  KeyPointValue 同士の比較

  kp1 < kp2 なら true を返す
*/
inline bool operator<( const KeyPointValue& kp1, const KeyPointValue& kp2 )
{
  if ( kp1.p().x == kp2.p().x )
    return( kp1.p().y < kp2.p().y );
  else
    return( kp1.p().x < kp2.p().x );
}

/**
   KeyPointValue 同士の等号チェック用関数オブジェクト
**/
struct EqKeyPointValue
{
  const KeyPointValue& kp_;

  /*
    比較対象の KeyPointValue を指定して構築

    kp : 比較対象の KeyPointValue
  */
  EqKeyPointValue( const KeyPointValue& kp ) : kp_( kp ) {}

  /*
    両者の KeyPointValue が等しければ true を返す
  */
  bool operator()( const KeyPointValue& kp )
  {
    return( kp_.p().x == kp.p().x && kp_.p().y == kp.p().y && kp_.dogIndex() == kp.dogIndex() );
  }
};

// KeyPointValue の配列の型 
using KeyPoint = std::set< KeyPointValue >;

/*
  DetectExtreme : 極値の検出

  dog : DoG 画像の配列
  k : 注目画素のある DoG 画像の添え字
  dest : 極値を登録する配列へのポインタ
*/
void DetectExtreme( const vector< Matrix< double > >& dog, vector< Matrix< double > >::size_type k,
                    KeyPoint* dest )
{
  for ( Matrix< double >::size_type y = 1 ; y + 1 < dog[k - 1].rows() ; ++y ) {
    for ( Matrix< double >::size_type x = 1 ; x + 1 < dog[k - 1].cols() ; ++x ) {
      double cp = dog[k][y][x]; // 注目画素
      bool isMax = true;        // 注目画素は極大値か?
      bool isMin = true;        // 注目画素は極小値か?
      for ( Matrix< double >::size_type dy = y - 1 ; dy <= y + 1 ; ++dy ) {
        for ( Matrix< double >::size_type dx = x - 1 ; dx <= x + 1 ; ++dx ) {
          for ( auto l = k - 1 ; l <= k + 1 ; ++l ) {
            if ( l == k && dx == x && dy == y ) continue;
            if ( cp < dog[l][dy][dx] ) {
              isMax = false;
              if ( ! isMin ) break;
            } else if ( cp > dog[l][dy][dx] ) {
              isMin = false;
              if ( ! isMax ) break;
            }
          }
          if ( ! ( isMax || isMin ) ) break;
        }
        if ( ! ( isMax || isMin ) ) break;
      }
      // 変化ない場合は isMax && isMin == true
      if ( ( isMax || isMin ) && ( ! ( isMax && isMin ) ) ) {
        // DoG 画像のサイズに対する原画像の倍率
        size_t ratio = std::pow( 2, k / ( KeyPointValue::scalePerOctave() + 2 ) );
        // 原画像にサイズ変換したときのピクセルの範囲内に、すでに登録されたピクセルが存在するかをチェック
        bool notMatched = true;
        for ( Coord< int > c( 0, y * ratio ) ; c.y < static_cast< int >( ( y + 1 ) * ratio ) ; ++( c.y ) ) {
          for ( c.x = x * ratio ; c.x < static_cast< int >( ( x + 1 ) * ratio ) ; ++( c.x ) ) {
            if ( dest->find( c ) != dest->end() ) {
              notMatched = false;
              break;
            }
          }
          if ( ! notMatched ) break;
        }
        if ( notMatched ) {
          KeyPointValue kp( Coord< int >( x, y ), k );
          dest->insert( kp );
        }
      }
    }
  }
}

/*
  DetectExtreme : DoG 画像からの極値検出

  dog : DoG 画像の配列
  dest : 検出した極値を登録する配列へのポインタ
*/
void DetectExtreme( const vector< Matrix< double > >& dog, KeyPoint* dest )
{
  const size_t DOG_COUNT = KeyPointValue::scalePerOctave() + 2; // オクターブあたりの DoG 画像の枚数

  for ( vector< Matrix< double > >::size_type o = 0 ; o + DOG_COUNT <= dog.size() ; o += DOG_COUNT ) {
    for ( vector< Matrix< double > >::size_type k = o + 1 ; k < o + DOG_COUNT - 1 ; ++k ) {
      DetectExtreme( dog, k, dest );
    }
  }
}

クラス KeyPointValue は、キーポイントを保持するために使います。メンバ変数の scalePerOctave_ はオクターブごとのスケール数を保持するためのものですが、これは全キーポイントに共通なため静的メンバ変数としています。dogIndex_ は極値を抽出した DoG 画像の添字、dogP_ はその座標を保持する変数です。dogP_ には、原画像ではなく DoG 画像での座標を代入します。つまり、画像が縮小された場合、その大きさでの座標値がそのまま入るということになります。そのため、原画像に合わせた座標に変換するメンバ関数 p が用意されています。
その他に、サブピクセルの座標とスケールを保持するためのメンバ変数 subP_, subScale_ がありますが、これらについては後述します。

極値を求めるための関数は上側の DetectExtreme です。注目画素の近傍にある計 26 点の画素と値を比較して、全ての点より大きければ極大、小さければ極小と判定します。注意点として、異なるスケールの同じ箇所で極値が見つかった場合、スケールの小さい側を採用します。そのため、極値と判定された後で、同じ箇所にすでに登録済みのキーポイントがないかチェックする必要があります。


キーポイントが得られたら、次に極値とその位置を正確に求める作業を行います。このとき、コントラストの小さい ( 極値の値が低い ) 箇所についてはノイズに反応しやすいため削除します。

DoG フィルタの結果を D( x, y, σ ) としたとき、そのテーラー - マクローリン展開による 2 次までの展開式は ( *2-1 )

D( x + Δx, y + Δy, σ + Δσ )=D( x, y, σ ) + [ Δx( ∂ / ∂x ) + Δy( ∂ / ∂y ) + Δσ( ∂ / ∂σ ) ]D( x, y, σ )
 + [ Δx( ∂ / ∂x ) + Δy( ∂ / ∂y ) + Δσ( ∂ / ∂σ ) ]2D( x, y, σ ) / 2

で表されます。x + Δx → x, y + Δy → y, σ + Δσ → σ と置き直して、α = x - Δx, β = y - Δy, γ = σ - Δσ とすると、

D( x, y, σ )=D( α, β, γ ) + [ ( ∂ / ∂x )( x - α ) + ( ∂ / ∂y )( y - β ) + ( ∂ / ∂σ )( σ - γ ) ]D( α, β, γ )
 + [ ( ∂ / ∂x )( x - α ) + ( ∂ / ∂y )( y - β ) + ( ∂ / ∂σ )( σ - γ ) ]2D( α, β, γ ) / 2

という形で表すこともできます。この式を、x, y, σ で偏微分すると、

( ∂ / ∂x )D( x, y, σ )=( ∂ / ∂x )D( α, β, γ ) + ( ∂2 / ∂x2 )D( α, β, γ )( x - α )
 + ( ∂2 / ∂x∂y )D( α, β, γ )( y - β ) + ( ∂2 / ∂σ∂x )D( α, β, γ )( σ - γ )
( ∂ / ∂y )D( x, y, σ )=( ∂ / ∂y )D( α, β, γ ) + ( ∂2 / ∂y2 )D( α, β, γ )( y - β )
 + ( ∂2 / ∂x∂y )D( α, β, γ )( x - α ) + ( ∂2 / ∂y∂σ )D( α, β, γ )( σ - γ )
( ∂ / ∂σ )D( x, y, σ )=( ∂ / ∂σ )D( α, β, γ ) + ( ∂2 / ∂σ2 )D( α, β, γ )( σ - γ )
 + ( ∂2 / ∂y∂σ )D( α, β, γ )( y - β ) + ( ∂2 / ∂σ∂x )D( α, β, γ )( x - α )

なので、D( α, β, γ ) を簡潔に D で表せば、( ( ∂ / ∂x )D( x, y, σ ), ( ∂ / ∂y )D( x, y, σ ), ( ∂ / ∂σ )D( x, y, σ ) ) = 0 のとき(つまり極値において)

-∂D / ∂x=( ∂2D / ∂x2 )( x - α ) + ( ∂2D / ∂x∂y )( y - β ) + ( ∂2D / ∂σ∂x )( σ - γ )
-∂D / ∂y=( ∂2D / ∂x∂y )( x - α ) + ( ∂2D / ∂y2 )( y - β ) + ( ∂2D / ∂y∂σ )( σ - γ )
-∂D / ∂σ=( ∂2D / ∂σ∂x )( x - α ) + ( ∂2D / ∂y∂σ )( y - β ) + ( ∂2D / ∂σ2 )( σ - γ )

と連立方程式になります。( α, β, γ ) はキーポイントの座標値、( x, y, σ ) が実際に極値をとる座標とすれば、連立方程式で求められる ( x - α, y - β, σ - γ ) はその差分であり、文献ではサブピクセルと呼ばれています。この連立方程式は、行列で表すと

|2D / ∂x2,2D / ∂x∂y,2D / ∂σ∂x||x - α|= -|∂D / ∂x|
|2D / ∂x∂y,2D / ∂y2,2D / ∂y∂σ||y - β||∂D / ∂y|
|2D / ∂σ∂x,2D / ∂y∂σ,2D / ∂σ2||σ - γ||∂D / ∂σ|

となるので、

|x - α|= -|2D / ∂x2,2D / ∂x∂y,2D / ∂σ∂x|-1|∂D / ∂x|
|y - β||2D / ∂x∂y,2D / ∂y2,2D / ∂y∂σ||∂D / ∂y|
|σ - γ||2D / ∂σ∂x,2D / ∂y∂σ,2D / ∂σ2||∂D / ∂σ|

を計算すれば、連立方程式の解法を使うことなく解を得ることができます。

ところで、テーラー - マクローリンの展開式は、( x, y, σ )T = x, ( α, β, γ )T = x^ で表し、

2D / ∂x2=|2D / ∂x2,2D / ∂x∂y,2D / ∂σ∂x|
|2D / ∂x∂y,2D / ∂y2,2D / ∂y∂σ|
|2D / ∂σ∂x,2D / ∂y∂σ,2D / ∂σ2|
∂D / ∂x = ( ∂D / ∂x, ∂D / ∂y, ∂D / ∂σ )T

とすることで

D(x) = D(x^) + ( ∂D(x^) / ∂x )T( x - x^ ) + ( x - x^ )T( ∂2D(x^) / ∂x2 )( x - x^ ) / 2

と表せます。連立方程式は

x - x^ = -( ∂2D(x^) / ∂x2 )-1( ∂D(x^) / ∂x )

なので、これを一部だけに代入すれば ( ∂2D(x^) / ∂x2 が対称行列であることに注意して )

D(x)=D(x^) + ( ∂D(x^) / ∂x )T( x - x^ ) - ( ∂D(x^) / ∂x )T( ∂2D(x^) / ∂x2 )-1( ∂2D(x^) / ∂x2 )( x - x^ ) / 2
=D(x^) + ( ∂D(x^) / ∂x )T( x - x^ ) - ( ∂D(x^) / ∂x )T( x - x^ ) / 2
=D(x^) + ( ∂D(x^) / ∂x )T( x - x^ ) / 2

という結果が得られます。D(x^) はキーポイントにおける DoG の値、∂D(x^) / ∂x は連立方程式の計算で求められ、x - x^ はその解なので、これらを使って極値における DoG の値を得ることができます。

こうして得られた極値での DoG の値がしきい値よりも小さければ、コントラストが小さいとみなして、そのキーポイントは削除します。文献では、DoG の値が 0 から 1 までの範囲だとした場合、0.03 をしきい値としています。


サブピクセルを求め、コントラストの低いキーポイントを削除するためのサンプル・プログラムを以下に示します。

/*
  CalcHessian : へッセ行列の計算 ( 3 x 3 )

  dog : DoG 画像の配列
  kp : 対象のKeyPointValue
  s0 : ガウシアンフィルタの標準偏差の初期値
  ds : スケール成分の増分
  dxx, dxy, dyy, dxs, dys, dss : 計算したへッセ行列の要素を保持する変数へのポインタ
*/
void CalcHessian( const vector< Matrix< double > >& dog, const KeyPointValue& kp,
                  double s0, double ds,
                  double* dxx, double* dxy, double* dyy, double* dxs, double* dys, double* dss )
{
  Coord< int > c = kp.dogP(); // DoG 画像上の座標
  size_t k = kp.dogIndex();   // DoG 画像の添字
  size_t lc = kp.scale();     // 実際のスケール

  CalcHessian( dog[k], c, dxx, dxy, dyy );

  *dxs = ( dog[k + 1][c.y][c.x + 1] - dog[k + 1][c.y][c.x - 1]
         - dog[k - 1][c.y][c.x + 1] + dog[k - 1][c.y][c.x - 1] ) /
         ( 2 * s0 * ( std::sqrt( std::pow( ds, lc + 1 ) ) - std::sqrt( std::pow( ds, lc - 1 ) ) ) );
  *dys = ( dog[k + 1][c.y + 1][c.x] - dog[k + 1][c.y - 1][c.x]
         - dog[k - 1][c.y + 1][c.x] + dog[k - 1][c.y - 1][c.x] ) /
         ( 2 * s0 * ( std::sqrt( std::pow( ds, lc + 1 ) ) - std::sqrt( std::pow( ds, lc - 1 ) ) ) );
  *dss = ( ( dog[k + 1][c.y][c.x] - dog[k][c.y][c.x] ) /
         ( s0 * ( std::sqrt( std::pow( ds, lc + 1 ) ) - std::sqrt( std::pow( ds, lc ) ) ) )
         - ( dog[k][c.y][c.x] - dog[k - 1][c.y][c.x] ) /
         ( s0 * ( std::sqrt( std::pow( ds, lc ) ) - std::sqrt( std::pow( ds, lc - 1 ) ) ) ) ) /
         ( s0 * ( std::sqrt( std::pow( ds, lc ) ) - std::sqrt( std::pow( ds, lc - 1 ) ) ) );
}

/*
  CalcDiff : 微分値の計算

  dog :DoG 画像の配列
  kp : 対象の KeyPointValue
  s0 : ガウシアンフィルタの標準偏差の初期値
  d : スケール成分の増分
  dx, dy, ds : 計算した微分値を保持する変数へのポインタ
*/
void CalcDiff( const vector< Matrix< double > >& dog, const KeyPointValue& kp,
               double s0, double d, double* dx, double* dy, double* ds )
{
  Coord< int > c = kp.dogP(); // DoG 画像上の座標
  size_t k = kp.dogIndex();   // DoG 画像の添字
  size_t lc = kp.scale();     // 実際のスケール

  *dx = ( dog[k][c.y][c.x + 1] - dog[k][c.y][c.x - 1] ) / 2;
  *dy = ( dog[k][c.y + 1][c.x] - dog[k][c.y - 1][c.x] ) / 2;
  *ds = ( dog[k + 1][c.y][c.x] - dog[k - 1][c.y][c.x] ) /
        ( s0 * ( std::sqrt( std::pow( d, lc + 1 ) ) - std::sqrt( std::pow( d, lc - 1 ) ) ) );
}

/*
  GetSubPixel : サブピクセルの計算

  dog : DoG 画像の配列
  kp : 対象の KeyPointValue
  s0 : ガウシアンフィルタの標準偏差の初期値
  d : スケールの増分
  x, y, s : 取得するサブピクセル
  dx, dy, ds : 勾配ベクトルの成分
*/
bool GetSubPixel( const vector< Matrix< double > >& dog, const KeyPointValue& kp,
                  double s0, double d,
                  double* x, double* y, double* s, double* dx, double* dy, double* ds )
{
  double dxx, dxy, dyy, dxs, dys, dss; // ヘッセ行列の成分

  CalcHessian( dog, kp, s0, d, &dxx, &dxy, &dyy, &dxs, &dys, &dss );
  CalcDiff( dog, kp, s0, d, dx, dy, ds );

  // 行列式
  double det = dxx * dyy * dss + dxy * dys * dxs + dxs * dxy * dys
             - dxs * dyy * dxs - dxy * dxy * dss - dxx * dys * dys;
  // 逆行列の成分を計算
  double a00 = -( dyy * dss - dys * dys ) / det;
  double a10 = ( dxy * dss - dys * dxs ) / det;
  double a20 = -( dxy * dys - dyy * dxs ) / det;
  double a11 = -( dxx * dss - dxs * dxs ) / det;
  double a21 = ( dxx * dys - dxy * dxs ) / det;
  double a22 = -( dxx * dyy - dxy * dxy ) / det;

  // サブピクセルの座標値を求める
  *x = a00 * *dx + a10 * *dy + a20 * *ds;
  *y = a10 * *dx + a11 * *dy + a21 * *ds;
  *s = a20 * *dx + a21 * *dy + a22 * *ds;

  return( std::abs( *x ) < 0.5 && std::abs( *y ) < 0.5 && std::abs( *s ) < 0.5 );
}

/*
  EstimateCoord : 正確なサブピクセル位置の推定

  dog :DoG 画像の配列
  src : 推定対象のキーポイント
  s0 : ガウシアンフィルタの標準偏差の初期値
  dth : DoG のしきい値
*/
void EstimateCoord( const vector< Matrix< double > >& dog, KeyPoint* src, double s0, double dth )
{
  const int DOG_COUNT = KeyPointValue::scalePerOctave() + 2;

  double dx, dy, ds; // 勾配ベクトルの成分
  double x, y, s;    // 求めるサブピクセル
  double d = std::pow( 2.0, 1.0 / KeyPointValue::scalePerOctave() ); // スケールの増分

  KeyPoint cand( *src );
  src->clear();
  for ( auto i = cand.begin() ; i != cand.end() ; ++i ) {
    KeyPointValue kp = *i;
    bool isDeleted = false;
    vector< KeyPointValue > buff;
    while ( ! GetSubPixel( dog, kp, s0, d, &x, &y, &s, &dx, &dy, &ds ) ) {
      kp.setDogP( Coord< int >( kp.dogP().x + std::round( x ), kp.dogP().y + std::round( y ) ) );
      kp.setSubPixel( Coord< double >( x, y ), s );
      int mod = kp.dogIndex() % DOG_COUNT;
      int is = std::round( s );
      if ( mod + is < 1 ) { // DoG 画像の比率が変わる場合
        // 現 DoG の先頭と前の DoG の末尾をスキップ
        for ( ; is < 0 ; ++is ) {
          if ( --mod == 0 ) {
            if ( kp.dogIndex() < 2 ) {
              isDeleted = true;
              break;
            }
            kp.setDogIndex( kp.dogIndex() - 2 );
            kp.setDogP( Coord< int >( kp.dogP().x * 2, kp.dogP().y * 2 ) );
            mod = KeyPointValue::scalePerOctave();
          }
          if ( kp.dogIndex() < 1 ) {
            isDeleted = true;
            break;
          }
          kp.setDogIndex( kp.dogIndex() - 1 );
        }
      } else if ( mod + is > static_cast< int >( KeyPointValue::scalePerOctave() ) ) { // DoG 画像の比率が変わる場合
        // 現 DoG の先頭と前の DoG の末尾をスキップ
        for ( ; is > 0 ; --is ) {
          if ( ++mod > static_cast< int >( KeyPointValue::scalePerOctave() ) ) {
            kp.setDogIndex( kp.dogIndex() + 2 );
            kp.setDogP( Coord< int >( kp.dogP().x / 2, kp.dogP().y / 2 ) );
            mod = 1;
          }
          kp.setDogIndex( kp.dogIndex() + 1 );
          if ( kp.dogIndex() + 1 >= dog.size() ) {
            isDeleted = true;
            break;
          }
        }
      } else {
        kp.setDogIndex( kp.dogIndex() + is );
      }
      if ( isDeleted ) break;
      // 発散のチェック
      if ( kp.dogP().x < 1 || kp.dogP().x + 1 >= static_cast< int >( dog[kp.dogIndex()].cols() ) ||
           kp.dogP().y < 1 || kp.dogP().y + 1 >= static_cast< int >( dog[kp.dogIndex()].rows() ) ) {
        isDeleted = true;
        break;
      }
      // 循環していないかチェック
      // 循環していた場合は、最後の(最初に重複が見つかった)キーポイントを採用する
      auto buffIdx = std::find_if( buff.begin(), buff.end(), EqKeyPointValue( kp ) );
      if ( buffIdx != buff.end() )
        break;
      buff.push_back( kp );
    }

    if ( ! isDeleted ) {
      // サブピクセル上の強度
      double subD = dog[kp.dogIndex()][kp.dogP().y][kp.dogP().x] + ( dx * x + dy * y + ds * s ) / 2;
      if ( std::abs( subD ) >= dth ) {
        auto dup = src->find( kp );
        if ( dup != src->end() ) {
          if ( dup->dogIndex() > kp.dogIndex() ) {
            src->erase( dup );
            src->insert( kp );
          }
        } else {
          src->insert( kp );
        }
      }
    }
  }
}

連立方程式の係数部分は CalcHessian、右辺は CalcDiff で計算しています。係数行列は 3 次のヘッセ行列であり、2 次のヘッセ行列は後で紹介するエッジ判定のサンプル・プログラムで実装されるので、残りの σ を含む要素の係数を計算します。必要な値が得られたら、ヘッセ行列の逆行列の成分を求め、右辺の値と掛け合わせることでサブピクセルの座標値を求めます。この処理は GetSubPixel が行います。
GetSubPixel は、求めたサブピクセルのうちいずれか一つでも 0.5 以上になったら false を返すようになっています。サブピクセルが 0.5 を超えたということは、キーポイントは別の位置にあることを示しています。そのため、キーポイントを変更して再度補完処理を行う必要があります。

EstimateCoordは、この処理のメインルーチンです。キーポイントを配列から一つずつ取り出し、サブピクセルの座標値を求めて極値を計算して、しきい値以上なら採用するという処理を繰り返すのですが、前述の通り、キーポイントを変更する必要が生じた場合の処理を含んでいるため少し複雑な実装になっています。
座標については単純に入れ替えるだけで済みますが、スケール成分は画像の比率が変化する箇所をまたぐ場合、重複した画像をスキップする必要がある分、余計な処理が加わります。
気をつけなければならない点は他にもあって、極点が発散したり循環していないかチェックする必要もあります。発散のチェックは、極点がスケール・スペースの外側に出てしまったことで判定できます。循環は、前に求めた極点の位置を記録しておいて、新たに求めた点と同じ要素がないかチェックすることで判断できます。循環した場合、最初に重複が見つかった点を採用するようにしています。


DoG フィルタはエッジ抽出にも使われます。そのため、キーポイントの中にはエッジも多く含まれていることになります。しかし、特徴としてはブロブだけを利用したいので、次にエッジを削除することを検討します。以前「ヘシアンを利用したコーナー検出」で説明したように、ヘッセ行列

|2I/∂x2,2I/∂x∂y|
H = |2I/∂x∂y,2I/∂y2|

の行列式はガウス曲率の分子部分であり、同時に H の固有値の積を表しています。また、H のトレースは固有値の和に等しいので ( *2-2 )、固有値を α, β ( 但し α ≥ β ) とすると

αβ = ( ∂2I/∂x2 )( ∂2I/∂y2 ) - ( ∂2I/∂x∂y )2
α + β = ∂2I/∂x2 + ∂2I/∂y2

で得られます。ここで α = γβ ( γ ≥ 1 ) とおくと

( α + β )2 / αβ = ( γβ + β )2 / γβ2 = ( γ + 1 )2 / γ

となって、( α + β )2 / αβ が γ だけの式で表せます。αβ がガウス曲率に比例していることから、γ が大きいほど曲率の最大値と最小値の差が大きくなると考えれば、γ の値をしきい値としてエッジ判定ができるようになります。なぜなら、曲率は大雑把に言えば曲がり具合を表す値なので、曲率の最大と最小の差が大きい場合は方向によって曲がり方が大きくなったり小さくなったりしているのに対し、差が小さい場合はどの方向も曲率が変わらないため、コーナーのような箇所が想定されるからです。

γ は固有値を求めることなく行列式とトレースから求めることが可能です。従って、γ のしきい値を決めて、γ がそれ以上になったらエッジとみなして削除します。


エッジ判定のサンプル・プログラムを以下に示します。

/*
  CalcHessian : へッセ行列の計算 ( 2 x 2 )

  dog : DoG 画像
  c : 計算対象の座標
  h00, h01, h11 : 計算したへッセ行列の要素を保持する変数へのポインタ
*/
void CalcHessian( const Matrix< double >& dog, const Coord< int >& c,
                  double* h00, double* h01, double* h11 )
{
  *h00 = dog[c.y][c.x + 1] - 2 * dog[c.y][c.x] + dog[c.y][c.x - 1];
  *h01 = ( dog[c.y + 1][c.x + 1] - dog[c.y + 1][c.x - 1] - dog[c.y - 1][c.x + 1] + dog[c.y - 1][c.x - 1] ) / 4;
  *h11 = dog[c.y + 1][c.x] - 2 * dog[c.y][c.x] + dog[c.y - 1][c.x];
}

/*
  EliminateKeyPoint : 不要なキーポイントの削除

  dog : DoG 画像の配列
  key : キーポイントの入った配列
  rth : しきい値
*/
void EliminateKeyPoint( const vector< Matrix< double > >& dog, KeyPoint* key, double rth )
{
  double h00, h01, h11;
  double r = ( rth + 1 ) * ( rth + 1 ) / rth; // tr^2 / det に対するしきい値

  vector< KeyPointValue > deleted; // 削除対象の座標
  for ( auto k = key->begin() ; k != key->end() ; ++k ) {
    CalcHessian( dog[k->dogIndex()], k->dogP(), &h00, &h01, &h11 );
    double tr = h00 + h11;              // ヘッセ行列のトレース
    double det = h00 * h11 - h01 * h01; // ヘッセ行列の行列式
    if ( tr * tr / det >= r ) {
      deleted.push_back( *k );
    }
  }
  for ( auto k = deleted.begin() ; k != deleted.end() ; ++k )
    key->erase( *k );
}

キーポイントの削除は EliminateKeyPoint 関数で行っています。キーポイントごとにヘッセ行列を計算し、行列式とトレースを求めてしきい値と比較する処理を繰り返します。ヘッセ行列の計算では二階微分が必要になりますが、ここでは差分法を使って次のように求めています。

∂I / ∂x(x,y) = I( x + 1, y ) - I( x, y )
∂I / ∂x(x-1,y) = I( x, y ) - I( x - 1, y ) より

2I / ∂x2=∂I / ∂x(x,y) - ∂I / ∂x(x-1,y)
=I( x + 1, y ) - 2I( x, y ) + I( x - 1, y )

∂I / ∂x(x,y+1) = [ I( x + 1, y + 1 ) - I( x - 1, y + 1 ) ] / 2
∂I / ∂y(x,y-1) = [ I( x + 1, y - 1 ) - I( x - 1, y - 1 ) ] / 2 より

2I / ∂x∂y=[ ∂I / ∂x(x,y+1) - ∂I / ∂x(x,y-1) ] / 2
=[ I( x + 1, y + 1 ) - I( x - 1, y + 1 ) - I( x + 1, y - 1 ) + I( x - 1, y - 1 ) ] / 4

2I / ∂x2 は差分法の前進差分、∂2I / ∂x∂y は中心差分を利用しています。

ヘッセ行列が得られたら、あとは求めた要素からトレース tr と行列式 det を計算し、tr * tr / det をしきい値 rth から計算した r = ( rth + 1 ) * ( rth + 1 ) / rth と比較して、r 以上ならエッジとみなして除外します。なお、文献では rth = 10 を最適値としているので、そのとき r = 12.1 になります。

*2-1)確率・統計 (18) 一般化線形モデル ( Generalized Linear Model )」の「補足 2) 多変量のテイラー - マクローリン展開」を参照

*2-2)画像処理 (5) コーナー検出」の「補足 3) 行列式・トレースと固有値の関係」を参照


3) キーポイントの方向の算出

得られたキーポイントの情報を使って特徴量を求める前に、その方向の算出を行います。特徴量を求める際に、その方向を使って軸を揃えれば、画像の回転に影響しない特徴量を得ることができます。

方向の算出にはガウシアン・フィルタで処理した画像データを使います。キーポイントの持つスケールの画像データを L( x, y ) で表したとき、キーポイントの周辺領域における勾配強度 m( x, y ) と勾配方向 θ( x, y ) は下式で求めることができます。

m( x, y ) = ( fx2 + fy2 )1/2
θ( x, y ) = tan-1( fy / fx )

但し fx = L( x + 1, y ) - L( x - 1, y )
fy = L( x, y + 1 ) - L( x, y - 1 )

キーポイント周辺の矩形領域に対して m と θ を計算し、θ について 36 分割した上で m の和を求めてヒストグラムを生成します。但し、m はキーポイントを中心としたガウシアン・フィルタで重み付けをするものとします。これは、中心から遠ざかるほど影響が小さくなるようにするためです。この時のガウシアン・フィルタの標準偏差は、文献ではキーポイントのスケールの 1.5 倍としています。

生成されたヒストグラムの中でピークを検出し、その中の最大値を得たら、その 80% 以上の値を持つピークを全てそのキーポイントの方向として採用します。従って、一つのキーポイントに対して複数の方向を持つ場合があります。

図 3-1. キーポイントの方向の算出
キーポイントの方向の算出

キーポイントの方向を算出するためのサンプル・プログラムを以下に示します。

const int BLOCK_SIZE = 4; // ブロックのサイズ
const int BLOCK_COUNT = 4; // 縦横のブロック数
const int IMAGE_SIZE = BLOCK_SIZE * BLOCK_COUNT; // 局所画像のサイズ
const int THETA_COUNT1 = 36; // ピーク検出用オリエンテーションの量子化数
const double THETA_PER_BIN = 2 * M_PI / THETA_COUNT1; // ピーク検出用オリエンテーションの単位角度

/**
   FeatureValue : 特徴量を保持する構造体
**/
struct FeatureValue
{
  KeyPointValue kp; // キーポイント
  double ori;       // ピークのある方向
  std::vector< double > hist; // 特徴量

  /*
    キーポイントとピークのある方向を指定して構築

    kp_ : キーポイント
    ori_ : ピークのある方向
  */
  FeatureValue( const KeyPointValue& kp_, double ori_ )
   : kp( kp_ ), ori( ori_ ) {}
};

/*
  CreateGaussianFilter : ガウシアン・フィルタの生成(二次元)

  フィルタのサイズは、行列 filter のサイズで決定する。

  sigma : σの値
  filter : 生成したフィルタを保持する行列へのポインタ
*/
void CreateGaussianFilter( double sigma, Matrix< double >* filter )
{
  assert( sigma > 0 );
  assert( filter != 0 );

  // フィルタのサイズ(片側のサイズとする)
  Coord< int > sz( filter->cols(), filter->rows() );
  if ( sz.x < 1 || sz.y < 1 ) return;
  sz.x = ( sz.x - 1 ) / 2;
  sz.y = ( sz.y - 1 ) / 2;

  for ( int y = -sz.y ; y <= sz.y ; ++y ) {
    for ( int x = -sz.x ; x <= sz.x ; ++x ) {
      double norm = std::pow( static_cast< double >( x ), 2.0 ) + std::pow( static_cast< double >( y ), 2.0 );
      double s2_2 = 2 * std::pow( sigma, 2.0 );
      (*filter)[y + sz.y][x + sz.x] = exp( -norm / s2_2 ) / ( M_PI * s2_2 );
    }
  }
}

/*
  GetPeakOfParabora : ピークを正確に求める

  y0, y1, y2 : 連なった 3 点の値
*/
double GetPeakOfParabora( double y0, double y1, double y2 )
{
  double a = ( y0 - 2 * y1 + y2 ) / ( 2 * std::pow( THETA_PER_BIN, 2 ) );
  double b = ( y0 - y2 ) / ( 2 * THETA_PER_BIN );

  return( b / ( 2 * a ) );
}

/*
  GetOrientation : ピークを持つ方向の探索

  gauss : ガウシアン・フィルタを行った画像データ
  kp : 対象のキーポイント
  s0 : スケールの初期値
  fv : 特徴量を登録する配列へのポインタ
*/
void GetOrientation( const vector< Matrix< double > >& gauss, const KeyPointValue& kp,
                     const double s0, vector< FeatureValue >* fv )
{
  const double THRESHOLD_RATIO = 0.8; // ピークの THRESHOLD_RATIO 分の方向は採用する
  const double SIGMA_MAG = 1.5;       // ガウス窓のσに掛ける定数

  vector< double > bin( THETA_COUNT1, 0 );
  Matrix< double > filter( IMAGE_SIZE, IMAGE_SIZE );

  const Matrix< double >& gp = gauss[kp.dogIndex()];
  Coord< int > sz( gp.cols(), gp.rows() );
  CreateGaussianFilter( s0 * std::pow( 2.0, static_cast< double >( kp.scale() ) / KeyPointValue::scalePerOctave() ) * SIGMA_MAG, &filter );
  Coord< int > p( kp.dogP() );

  // ヒストグラムの生成
  for ( Coord< int > c( 0, -IMAGE_SIZE / 2 ) ; c.y < IMAGE_SIZE / 2 ; ++( c.y ) ) {
    for ( c.x = -IMAGE_SIZE / 2 ; c.x < IMAGE_SIZE / 2 ; ++( c.x ) ) {
      double dx = gp[Mirror( p.y + c.y, sz.y )][Mirror( p.x + c.x + 1, sz.x )] -
                  gp[Mirror( p.y + c.y, sz.y )][Mirror( p.x + c.x - 1, sz.x )];
      double dy = gp[Mirror( p.y + c.y + 1, sz.y )][Mirror( p.x + c.x, sz.x )] -
                  gp[Mirror( p.y + c.y - 1, sz.y )][Mirror( p.x + c.x, sz.x )];
      double mag = std::sqrt( std::pow( dx, 2 ) + std::pow( dy, 2 ) );
      double theta = std::atan2( dy, dx );
      if ( theta < 0 ) theta += 2 * M_PI;
      size_t binNum = theta / THETA_PER_BIN;
      bin.at( binNum ) += filter[c.y + IMAGE_SIZE / 2][c.x + IMAGE_SIZE / 2] * mag;
    }
  }

  // ピークを検索する
  map< size_t, double > peak;
  double maxBin = 0;
  for ( size_t i = 0 ; i < THETA_COUNT1 ; ++i ) {
    if ( bin[i] > bin[( i + THETA_COUNT1 - 1 ) % THETA_COUNT1] &&
         bin[i] > bin[( i + 1 ) % THETA_COUNT1] ) {
      peak.insert( std::make_pair( i, bin[i] ) );
      if ( maxBin < bin[i] ) maxBin = bin[i];
    }
  }
  maxBin *= THRESHOLD_RATIO;

  for ( auto i = peak.begin() ; i != peak.end() ; ++i )
    if ( i->second >= maxBin )
      fv->push_back
        ( FeatureValue( kp, i->first * THETA_PER_BIN +
                        GetPeakOfParabora( bin[( i->first + THETA_COUNT1 - 1 ) % THETA_COUNT1],
                                           bin[i->first],
                                           bin[( i->first + 1 ) % THETA_COUNT1]
                                           ) ) );
}

FeatureValue は特徴量を保持するために使う構造体です。この中には、キーポイントの情報とピークのある方向、そして特徴量を持たせることができます。

キーポイントの周辺領域として、4 x 4 ( BLOCK_SIZE ) のサイズのブロックをさらに 4 x 4 ( BLOCK_COUNT ) 個並べた 16 x 16 ( IMAGE_SIZE ) の矩形を利用します。これは、文献の中での実験で実際に利用されたサイズになります。ブロックについては次の節で詳しく説明します。

GetOrientaion が、ピークを持つ方向を求めるための関数になりますが、二点ほど注意すべき箇所があります。まず、最後に FeatureValue を構築する箇所で GetPeakOfParabora を呼び出していますが、これはピークの位置をその近隣の値から補間する関数です。ピークをはさんだ左右を含めたヒストグラムの値が y0, y1, y2 の順番で並んでいたとします。x 軸は中央を原点として左右が -x, +x であるとしたとき、各点を通る曲線が二次関数 y = ax2 - bx + c で近似できるとすれば、

ax2 + bx + c = y0
c = y1
ax2 - bx + c = y2

という簡単な連立方程式になります。これを解けば

a = ( y0 - 2y1 + y2 ) / 2x2
b = ( y0 - y2 ) / 2x

という結果が得られ、ピークでの x の値 b / 2a を求めることができます。

もう一つの注意点は atan2 の使い方です。atan2( y, x ) で tan-1( y / x ) の値を求めることができますが、その値域は [ -π, π ] です。従って、[ 0, 2π ] の範囲にしたい場合は、負数だった時に 2π を加える処理が必要になります。標準関数には atan というものもあり、こちらは atan( t ) で tan-1( t ) の値を返しますが、その値域は [ -π / 2, π / 2 ] です。この違いは、atan の場合、t の値が y / x と -y / -x のどちらなのか区別できないことから生じます。


サンプル・プログラムを使って極値の抽出からキーポイントの方向の算出までを行った結果を以下に示します。利用した画像サンプルは、南カリフォルニア大学の The USC-SIPI Image Database にある "Splash" です。処理に際して、DoG 画像生成時に行うガウシアン・フィルタの標準偏差の初期値は 1.6、オクターブ数は 2、スケールを小さくする最小の画像サイズは 10、コントラストに対するしきい値は文献よりも小さく 0.01 x 256、エッジ判定のしきい値は文献の通り 10 としました。

図 3-2. Splash (実サイズ 512 x 512)
Splash

 

図 3-3. キーポイント抽出結果
極値抽出結果正確なキーポイントの抽出・低コントラスト除外
極値抽出正確なキーポイントの抽出・低コントラスト除外
エッジ除外・方向の算出
エッジ除外・方向の算出

画像の中で、キーポイントを示す円の大きさはスケール成分 σ の大きさを表しています。また、方向の算出結果は円の中心から円周への線分で表しています。
極値を抽出したときには画面全体にキーポイントが得られていますが、低コントラストなものを除外することで王冠の部分以外のキーポイントは消え去っています。さらにエッジを除外することで、画像の両端に残ったエッジ起因のキーポイントも消去されていることが分かります。

4) 特徴量の計算

最後に、特徴量の計算を行います。特徴量の計算に使うのは、前述の方向算出に用いたものと同じ領域になりますが、大きく異なるのは、その方向ベクトルに座標軸を合わせるという点です。こうすることで、画像が回転していたとしても同じ箇所が等しい特徴量を持つことが期待できます。なお、矩形の回転によって、領域に含まれる画素は四方で変化します。すると、その変化に特徴量が影響を受けるため、それを避ける目的で、領域の幅の半分の大きさを標準偏差とするガウシアン・フィルタで重み付けを行います。これにより、領域の四方の影響を小さくすることができます。

図 4-1. 特徴量の算出
特徴量の算出

勾配強度 m( x, y ) と勾配方向 θ( x, y ) の求め方は方向算出時と同様です。しかし、ここでは 4 x 4 のブロック毎にヒストグラムを作成します。また、ヒストグラムの分割数は 8 になります。
こうして得られるビンの数は 4 x 4 x 8 = 128 個になります。この 128 次元の値が特徴量となります。


特徴量を求めるためのサンプル・プログラムを以下に示します。

/*
  GetFeatureValue : 特徴量の抽出

  gauss : ガウシアン・フィルタを行った画像データ
  fv : 抽出した特徴量を登録する変数へのポインタ
  s0 : スケールの初期値
*/
void GetFeatureValue( const vector< Matrix< double > >& gauss, FeatureValue* fv, double s0 )
{
  double rot = fv->ori;
  Matrix< double > filter( IMAGE_SIZE, IMAGE_SIZE );

  const Matrix< double >& gp = gauss[( fv->kp ).dogIndex()];
  Coord< int > sz( gp.cols(), gp.rows() );
  CreateGaussianFilter( IMAGE_SIZE / 2, &filter );

  // DoG 画像に合わせて座標変換
  Coord< int > o( ( fv->kp ).dogP() );

  for ( Coord< int > b( 0, 0 ) ; b.y < BLOCK_COUNT ; ++( b.y ) ) {
    for ( b.x = 0 ; b.x < BLOCK_COUNT ; ++( b.x ) ) {
      vector< double > bin( THETA_COUNT2, 0 );
      for ( Coord< int > p( 0, 0 ) ; p.y < BLOCK_SIZE ; ++( p.y ) ) {
        for ( p.x = 0 ; p.x < BLOCK_SIZE ; ++( p.x ) ) {
          Coord< int > c( b.x * BLOCK_SIZE + p.x - IMAGE_SIZE / 2, b.y * BLOCK_SIZE + p.y - IMAGE_SIZE / 2 );
          Coord< int > rR( round( ( c.x + 1 ) * cos( rot ) - c.y * sin( rot ) ),
                           round( ( c.x + 1 ) * sin( rot ) + c.y * cos( rot ) ) );
          Coord< int > rL( round( ( c.x - 1 ) * cos( rot ) - c.y * sin( rot ) ),
                           round( ( c.x - 1 ) * sin( rot ) + c.y * cos( rot ) ) );
          Coord< int > rU( round( c.x * cos( rot ) - ( c.y - 1 ) * sin( rot ) ),
                           round( c.x * sin( rot ) + ( c.y - 1 ) * cos( rot ) ) );
          Coord< int > rD( round( c.x * cos( rot ) - ( c.y + 1 ) * sin( rot ) ),
                           round( c.x * sin( rot ) + ( c.y + 1 ) * cos( rot ) ) );
          double dx = gp[Mirror( o.y + rR.y, sz.y )][Mirror( o.x + rR.x, sz.x )] -
                      gp[Mirror( o.y + rL.y, sz.y )][Mirror( o.x + rL.x, sz.x )];
          double dy = gp[Mirror( o.y + rD.y, sz.y )][Mirror( o.x + rD.x, sz.x )] -
                      gp[Mirror( o.y + rU.y, sz.y )][Mirror( o.x + rU.x, sz.x )];
          double mag = std::sqrt( std::pow( dx, 2 ) + std::pow( dy, 2 ) );
          double theta = std::atan2( dy, dx );
          if ( theta < 0 ) theta += 2 * M_PI;
          size_t binNum = theta / THETA_PER_ORI;
          bin.at( binNum ) += filter[c.y + IMAGE_SIZE / 2][c.x + IMAGE_SIZE / 2] * mag;
        }
      }
      ( fv->hist ).insert( ( fv->hist ).end(), bin.begin(), bin.end() );
    }
  }
}

前節で示した関数 GetOrientaion によく似た構成ですが、大きな違いはピクセルを取得する範囲がピーク方向に回転している点です。rR, rL, rU, rD は、矩形領域を回転した時の各ピクセルの相対座標を表していますが、これらはキーポイントに対して毎回同じ値になるので、あらかじめテーブル化しておけば計算量を少なくすることができるでしょう。

こうして、スケールや回転に不変な特徴量が得られたのですが、画像の明るさに対してはまだ不変となっていません。画像の明るさの変化の影響を小さくするため、特徴量をノルムの総和で割って正規化します。画像全体が定数倍されたときの変化量は、やはり全体が定数倍された状態になるので、正規化することで影響を取り除くことができます。
しかし、局所的な明るさの変化にはまだ未対応の状態です。これは、勾配の大きさには強い変化をもたらす可能性がありますが、勾配の方向についてはそれほど影響を与えないことが期待できます。そこで、正規化した値をしきい値でカットして、再度正規化処理を行います。つまり、大きな勾配での値を一致させることはもはや重要ではなく、方向の分布をより重視することを意味します。しきい値は、文献の中での実験で 0.2 となっています。


正規化処理のためのサンプル・プログラムを以下に示します。

/*
  L2Norm : L2 ノルム用関数オブジェクト
*/
struct L2Norm
{
  /*
    operator() : 二乗値を加算する

    x : 被加算数
    y : 加算数 ( 二乗値の計算 )

    戻り値 : 計算結果
  */
  double operator()( double x, double y ) const
  { return( x + y * y ); }

  /*
    cleanUp : 計算後の後処理 ( 平方根の計算 )
  */
  double cleanUp( double d ) const
  { return( std::sqrt( d ) ); }
};

/*
  L1Norm : L1 ノルム用関数オブジェクト
*/
struct L1Norm
{
  /*
    operator() : 絶対値を加算する

    x : 被加算数
    y : 加算数 ( 絶対値の計算 )

    戻り値 : 計算結果
  */
  double operator()( double x, double y ) const
  { return( x + std::abs( y ) ); }

  /*
    cleanUp : 計算後の後処理 ( 何もしない )
  */
  double cleanUp( double d ) const
  { return( d ); }
};

/*
  Normalize : 正規化用関数

  ノルムを求めて各要素を割る

  fv : 正規化対象の配列
  op : ノルム計算用関数オブジェクト
*/
template< class Op >
void Normalize( vector< double >* fv, Op op )
{
  double norm = op.cleanUp( std::accumulate( fv->begin(), fv->end(), 0.0, op ) );

  std::transform( fv->begin(), fv->end(), fv->begin(), std::bind( std::divides< double >(), _1, norm ) );
}

/*
  CutOff : しきい値以上の要素をしきい値にカットする
*/
class CutOff
{
  double threshold_; // しきい値

  public:

  /*
    しきい値を指定して構築

    th : しきい値
  */
  CutOff( double th )
  : threshold_( th ) {}

  /*
    operator() : しきい値以上の値をしきい値にして返す
  */
  double operator()( double value ) const
  { return( std::min( threshold_, value ) ); }
};

/*
  NormalizeFeatureValue : 特徴量の正規化

  fv : 特徴量
  op : ノルム計算用関数オブジェクト
*/
template< class Op >
void NormalizeFeatureValue( vector< double >* fv, Op op )
{
  const double THRESHOLD = 0.2;

  CutOff cutOff( THRESHOLD );

  Normalize( fv, op );
  std::transform( fv->begin(), fv->end(), fv->begin(), cutOff );
  Normalize( fv, op );
}

ノルムの計算には、L2NormL1Norm の 2 種類を用意しています。L2Norm は二乗和を、L1Norm は絶対値の和を計算するよう operator() を定義して、STL ( Standard Template Library ) の accumulate 関数で処理することで結果が得られる仕組みです。L2Norm では最後に平方根を求める必要があるので、そのためのメンバ関数 cleanUp が用意されています。なお、L1Norm では後処理不要なので、cleanUp は引数の値をそのまま返すようになっています。

正規化処理は Normalize が行います。処理内容は非常に簡単で、先述の通り L2NormL1Norm を利用して accumulate 関数でノルムを計算し、transoform 関数を使ってノルムで各要素を除算するだけです。L2NormL1Norm はテンプレート引数 op で渡すようになっているので、呼び出し時に切り替えることができます。

CutOff は引数の値がしきい値以上だったらしきい値にして返す簡単な関数オブジェクトです。最後の NormalizeFeatureValue が特徴量を正規化処理するメインルーチンで、正規化処理の後にしきい値以上の値をしきい値にカットし、もう一度正規化処理を行います。


5) 性能評価

サンプル・プログラムを使って、画像の変形に対する特徴量の不変性を評価してみました。利用したサンプルは、南カリフォルニア大学の The USC-SIPI Image Database にある 12 枚の画像です。処理に際して、DoG 画像生成時に行うガウシアン・フィルタの標準偏差の初期値は 1.6、オクターブ数は 2、スケールを小さくする最小の画像サイズは 10、コントラストに対するしきい値は 0.01 x 256、エッジ判定のしきい値は 10 としています。

例として、"Splash" ( 4.2.01 ) に対して 30 度時計回りに回転した画像と元の画像の特徴量を求めたとします。得られたキーポイントの数は、元の画像が 367、回転した画像が 430 でした。次に元の画像の各キーポイントについて、特徴量とのノルムが最小のキーポイントと、2 番めに小さいキーポイントを回転した画像から抽出します。そして、ノルムの最小値が、2 番めに小さいノルムの 0.75 倍より小さければ、同じキーポイントだとみなしてペアを形成します。この処理によって、206 個のペアが形成されました。最後に、それぞれのキーポイントの座標が一致しているかをチェックします。回転した画像のキーポイントは、反時計回りに 30 度回転させれば元の画像のキーポイントと一致するはずですが、処理時の誤差を加味して、ズレ量が x, y それぞれ 3 ピクセル以内なら一致したとみなします。その結果、座標の一致したキーポイントは 168 個得られました。

以上の結果から、二つの画像のキーポイントの合計 ( 但し、一致したキーポイントは一つとカウントする ) に対する一致したキーポイントの比率をマッチ率、特徴量から得られたペアの中で、実際に一致したキーポイントの比率を正解率として求めます。この例では

マッチ率 = 168 / ( 367 + 430 - 168 ) = 26.7%
正解率 = 168 / 206 = 81.6%

という結果になります。

表 5-1. SIFT 特徴量の一致率
4.1.014.1.024.1.034.1.044.1.054.1.064.1.074.1.084.2.014.2.034.2.054.2.06平均
30度回転マッチ率38.134.127.230.320.032.736.445.126.729.935.835.333.1
正解率93.387.978.885.967.788.588.289.881.691.887.890.888.4
1.5倍拡大マッチ率9.05.79.09.87.96.78.58.35.41.73.52.84.3
正解率63.871.170.062.768.466.766.063.043.034.350.345.053.0
30度回転+
1.5倍拡大
マッチ率4.83.53.83.83.43.14.77.02.21.11.41.42.3
正解率55.855.350.053.152.050.054.547.627.230.823.627.837.1
180度回転マッチ率41.236.947.343.843.632.246.739.134.639.243.138.639.5
正解率90.189.489.892.985.083.785.787.183.592.390.889.889.3
2倍拡大マッチ率17.413.921.122.611.113.821.330.917.29.013.412.913.5
正解率91.090.995.195.077.189.088.589.670.080.978.882.882.8
20%せん断マッチ率35.233.230.025.223.832.634.933.527.026.531.431.730.1
正解率85.783.076.584.975.084.988.084.279.282.384.086.983.8
100%せん断マッチ率0.81.93.72.61.01.60.92.11.80.30.70.20.9
正解率8.619.420.819.46.714.59.523.515.83.87.52.39.8
水平反転マッチ率2.73.06.14.75.24.45.45.74.42.21.51.92.8
正解率25.627.837.546.238.540.734.338.128.025.017.723.227.9
垂直反転マッチ率3.24.26.36.07.32.83.05.54.92.51.81.62.9
正解率26.535.937.542.945.226.925.938.729.830.419.619.428.0
ランダムノイズマッチ率13.48.529.919.118.931.025.733.714.417.117.121.319.3
正解率72.949.280.078.571.988.882.488.165.075.676.780.077.2
明度1.2倍マッチ率92.486.280.278.967.277.395.894.955.682.471.377.577.1
正解率100.098.8100.098.494.495.398.399.193.298.796.898.597.7
明度2倍マッチ率58.571.624.029.110.423.223.616.717.815.46.720.519.8
正解率96.597.465.279.157.177.878.865.672.471.839.373.273.0

 

図 5-1. マッチ率・正解率比較グラフ
マッチ率・正解率比較グラフ

表の縦軸 ( グラフの横軸 ) には画像に対して行った加工内容が示してあります。せん断は画像を斜めに傾けたような加工を表し、20% のせん断で画像の幅の 20%、100% では画像の幅と同じだけ上辺が右側に移動した状態になります。また、回転・拡大・せん断のいずれも補間処理はしていません ( 正確には最近傍法 )。ランダムノイズは、画素数の 10% に対してランダムな場所にランダムな色で点を描画した画像、明度は YUV 変換した上で Y 成分を 1.2 倍または 2 倍にしてから RGB に戻してプロットする操作を行っています。

マッチ率は、得られた特徴量の中でどれだけ共通のものがあるかを示す指標になります。この値は明度 1.2 倍を除き、50% を超えるものはありませんでした。極端に低いのは、拡大・100% せん断・水平反転・垂直反転で、拡大処理以外は予想通りの結果となっています。正解率は、特徴量の近いキーポイントが同じ位置を示している割合を表しており、マッチ率の高いものは正解率も高い傾向にあります。

予想に反して、拡大処理でのマッチ率や正解率があまり高くないという結果になりました。特に 1.5 倍という中途半端な比率で拡大した場合は正解率が 5 割程度で回転処理などに比べて低いです。2 倍に拡大した場合は高い正解率が得られましたが、マッチ率はやはりあまり高くありません。せん断処理ではマッチ率や正解率は低いだろうと予想していましたが、20% 程度ならマッチさせることは可能なようです。面白いのは垂直反転と 180 度回転の比較で、似たような加工であるにもかかわらず垂直反転のほうが極端に低いという結果になっています。垂直反転の場合、キーポイントの方向は元画像のそれと裏返しの関係になるので、特徴量が一致しなくなるわけですね。


SIFT は、画像のマッチングの他にも、物体認識や画像分類などに応用されてきました。現在では畳み込みニューラル・ネットワークの誕生で影が薄くなったように感じますが、学習が不要ですぐに利用できるという強みが SIFT にはあります。昔は特許が取られていたので利用しにくい面もありました。しかし、2020/03/06 に特許が切れているので、今は自由に利用可能な状態です。ぜひ、この記事を参考に一度試してみて下さい。


補足 1) ガウシアン・フィルタの性質

x( t ) を f( t ) で畳み込み積分した結果を w( t ) とし、さらに w( t ) を g( t ) で畳み込み積分した結果を y( t ) とします。このとき、

w( t ) = ∫{-∞→∞} f( s )x( t - s ) ds

y( t ) = ∫{-∞→∞} g( u )w( t - u ) du

下側の式に上側の式を代入して

y( t ) = ∫{-∞→∞} g( u ) ∫{-∞→∞} f( s )x( t - u - s ) ds du

v = u + s とすると

y( t )=∫{-∞→∞} g( u ) ∫{-∞→∞} f( v - u )x( t - v ) dv du
=∫{-∞→∞} [ ∫{-∞→∞} g( u )f( v - u ) du ] x( t - v ) dv

よって、

h( v ) = ∫{-∞→∞} g( u )f( v - u ) du

とすれば

y( t ) = ∫{-∞→∞} h( v )x( t - v ) dv

となって、直列に接続したフィルタによる畳み込み積分は、2 つのフィルタの畳み込み積分によって得られたフィルタによるそれと等しくなることが分かります。

特に 2 つのフィルタがガウシアン・フィルタの場合、

f( t ) = [ 1 / ( 2π )1/2σ1 ]exp( -t2 / 2σ12 )

g( t ) = [ 1 / ( 2π )1/2σ2 ]exp( -t2 / 2σ22 )

として

h( t )=∫{-∞→∞} g( u )f( t - u ) du
=[ 1 / 2πσ1σ2 ]∫{-∞→∞} exp( -u2 / 2σ12 )exp( -( t - u )2 / 2σ22 ) du
=[ 1 / 2πσ1σ2 ]∫{-∞→∞} exp( -u2 / 2σ12 - ( t - u )2 / 2σ22 ) du
=[ 1 / 2πσ1σ2 ]∫{-∞→∞} exp( -( 1 / 2σ12 + 1 / 2σ22 )u2 + tu / σ22 - t2 / 2σ22 ) du
=[ 1 / 2πσ1σ2 ]∫{-∞→∞} exp( -( 1 / 2σ12 + 1 / 2σ22 ){ u - [ σ12 / ( σ12 + σ22 ) ]t }2 + ( 1 / 2σ12 + 1 / 2σ22 )[ σ12 / ( σ12 + σ22 ) ]2t2 - t2 / 2σ22 ) du
=[ 1 / 2πσ1σ2 ]∫{-∞→∞} exp( -( 1 / 2σ12 + 1 / 2σ22 ){ u - [ σ12 / ( σ12 + σ22 ) ]t }2 - [ 1 / 2( σ12 + σ22 ) ]t2 ) du
=[ 1 / 2πσ1σ2 ]exp( - [ 1 / 2( σ12 + σ22 ) ]t2 )∫{-∞→∞} exp( -( 1 / 2σ12 + 1 / 2σ22 ){ u - [ σ12 / ( σ12 + σ22 ) ]t }2 ) du

積分は、

∫{-∞→∞} exp( -( t - b )2 / 2c2 ) dt = c( 2π )1/2

より

∫{-∞→∞} exp( -( 1 / 2σ12 + 1 / 2σ22 ){ u - [ σ12 / ( σ12 + σ22 ) ]t }2 ) du = [ σ12σ22 / ( σ12 + σ22 ) ]1/2( 2π )1/2

なので

h( t ) = [ 1 / ( 2π )1/2( σ1 + σ2 )1/2 ]exp( - [ 1 / 2( σ12 + σ22 ) ]t2 )

となって、( σ1 + σ2 )1/2 を新たな標準偏差とするガウシアン・フィルタと等しくなります。

以上の結果を整理すると、ガウシアン・フィルタで 2 回畳み込み積分をする代わりに、それぞれのフィルタの標準偏差の二乗和の平方根を新たな標準偏差としたガウシアン・フィルタで 1 回畳み込み積分をしても結果は同じであることになります。SIFT では、複数の標準偏差 ( σ, kσ, k2σ ... ) で何度もガウシアン・フィルタによる畳み込み積分を行うわけですが、例えば原画像に対して σ で処理した後、同じ原画像を kσ で処理する代わりに

σ' = ( ( kσ )2 - σ2 )1/2

を使って σ で処理した結果を畳み込み積分しても結果は変わらないことになります。従って、前回の処理結果を使って繰り返し処理を行うことができるわけです。



<参考文献>
  1. 「画像認識」 原田達也著 ( 講談社 )
  2. 機械知覚 & ロボティクスグループ 中部大学」 - 「Gradient ベースの特徴抽出 - SIFT と HOG -
  3. The University of British Columbia」 - 「David Lowe」 - 「Distinctive Image Features from Scale-Invariant Keypoints

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system