画像処理

(5) コーナー検出

画像認識では主に、画像の持つ特徴を表す代表点の抽出が最初に行われます。代表点として使われるものには、前章で紹介したエッジの他に、物体の角を表すコーナーがあります。今回は、コーナー検出の方法について紹介したいと思います。

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

1) ヘシアンを利用したコーナー検出

画素を平面 ( x, y ) の強度 I による曲面で表したとき、「ガウス曲率 ( Gaussian Curvature )」K が次式で定義できます ( 補足 1 )。

K = ( IxxIyy - Ixy2 ) / ( 1 + Ix2 + Iy2 )2

但し、Ixx = ∂2I / ∂x2、Iyy = ∂2I / ∂y2、Ixy = ∂2I / ∂x∂y をそれぞれ表します。コーナーではガウス曲率は大きくなります。そこで、ガウス曲率の分子部分を使い、この値が極大となる点をコーナーとして検出する手法を「ヘシアン・コーナー検出 ( Hessian Corner Detector )」といいます。ガウス曲率の分子部分は、ちょうど I( x, y ) の x, y による二階微分からなる下記のような行列「ヘッセ行列 ( Hessian Matrix )」の行列式なので、この名が付けられています。

|Ixx, Ixy|
|Ixy, Iyy|

早速、サンプル・プログラムを示したいと思います。

// 検出したコーナーを保持するためのコンテナ
typedef std::multimap< double, Coord< int >, std::greater< double > > Corner;

/*
  GetExtremeValue : 極値の補正

  中央 iC、左側 iL、右側 iR、上側 iU、下側 iD を通る楕円体の極大値を計算する
*/
double GetExtremeValue( double iC, double iL, double iR, double iU, double iD )
{
  double a = ( iL + iR - iC * 2 ) / 2;
  double b = ( iU + iD - iC * 2 ) / 2;

  return( iC
          - ( ( a == 0 ) ? 0 : a * std::pow( ( ( iR - iC ) / a - 1 ) / 2, 2 ) )
          - ( ( b == 0 ) ? 0 : b * std::pow( ( ( iD - iC ) / b - 1 ) / 2, 2 ) )
          );
}

/*
  Corner_Hessian : ヘシアンを利用したコーナー検出

  data : 検出対象の画素データ
  corner : 検出したコーナーの座標値を保持する multimap ( キー : 極値 ; 値 : 座標 )
  threshold : 極値抽出時のしきい値 ( threshold 以上の差異がなければ極値とみなさない )
*/
void Corner_Hessian( const Matrix< double >& data,
                     Corner* corner, double threshold )
{
  corner->clear();

  // サイズが 5 より小さい画像に対しては検出しない
  if ( data.rows() < 5 || data.cols() < 5 ) return;

  // ヘシアン値
  Matrix< double > det( data.rows() - 2, data.cols() - 2 );

  // ヘシアンの計算
  Coord< int > c;
  for ( c.y = 1 ; c.y < static_cast< int >( data.rows() - 1 ) ; ++( c.y ) ) {
    for ( c.x = 1 ; c.x < static_cast< int >( data.cols() - 1 ) ; ++( c.x ) ) {
      double ixx = data[c.y][c.x - 1] + data[c.y][c.x + 1] - data[c.y][c.x] * 2;
      double iyy = data[c.y - 1][c.x] + data[c.y + 1][c.x] - data[c.y][c.x] * 2;
      double ixy = ixx + iyy;
      det[c.y - 1][c.x - 1] = ixx * iyy - ixy * ixy;
    }
  }

  // 極大値を求める
  for ( c.y = 1 ; c.y < static_cast< int >( det.rows() - 1 ) ; ++( c.y ) ) {
    for ( c.x = 1 ; c.x < static_cast< int >( det.cols() - 1 ) ; ++( c.x ) ) {
      if ( det[c.y][c.x - 1] + threshold >= det[c.y][c.x] || det[c.y][c.x] < det[c.y][c.x + 1] + threshold )
        continue;
      if ( det[c.y - 1][c.x] + threshold < det[c.y][c.x] && det[c.y][c.x] >= det[c.y + 1][c.x] + threshold ) {
        corner->insert( std::pair< double, Coord< int > >
                        ( GetExtremeValue( det[c.y][c.x], det[c.y][c.x - 1], det[c.y][c.x + 1], det[c.y - 1][c.x], det[c.y + 1][c.x] ),
                          Coord< int >( c.x + 1, c.y + 1 ) ) );
      }
    }
  }
}

メインとなる関数は Corner_Hessian です。引数として二次元配列 Matrix を型とする data と、検出したコーナーの座標値を取り込むための corner、そして極値判定のためのしきい値 threshold の三つがあります。corner の型は、キーを double ( これは極大値の値を表します )、値を Coord< int > ( 極大値を持つ座標を表します ) とする multimap で、比較用の叙述関数は greater なので、極大値の値が大きいものほど先頭側に並ぶことになります。なお、Matrix クラスは列数を表す cols と行数を表す rows の 2 つのメンバ関数を持ちます。また、Coord 型は 2 つのメンバ変数 xy を持つ構造体です。
処理の内容はいたって単純で、最初に各画素のヘシアンを計算し、その値が周囲よりも大きくなる位置 ( このとき、周囲との差が threshold 以上でないものは極大値とは見なさないようにしています ) を corner に抽出するだけです。

注意点として、ヘシアンの計算に周囲 2 ピクセル分の画素が必要になるため、画像の端 2 ピクセル分に対しては計算を行っていません。また極大値は、周囲の値を用いて補間してから抽出しています。この計算は関数 GetExtremeValue で行っていますが、具体的には中央・右・左・上・下 5 ピクセルを通る楕円体の極大値を計算して返しています ( 補足 2 )。


コーナー点に対して、検出された位置が近隣にたくさん固まっているような場合があります。その場合、周囲の検出点の中で最もコーナーらしい ( 上記処理の場合、ヘシアン値の最も大きい ) 点だけを検出することで、不要な点を削除することができます。そのような処理のサンプル・プログラムは次のようになります。

/*
  LessCoord : 座標値の大小判定をするための関数オブジェクト
*/
struct LessCoord
{
  /*
    operator() : 座標値の大小判定をする

    c0, c1 : 対象の座標
    戻り値 : c0 < c1 なら true
  */
  bool operator()( const Coord< int >& c0, const Coord< int >& c1 ) const
  {
    if ( c0.x == c1.x )
      return( c0.y < c1.y );
    else
      return( c0.x < c1.x );
  }
};

/*
  CornerSieve : 取得したコーナーをふるいにかける

  corner : 検出したコーナーの座標値を保持する multimap ( キー : 極値 ; 値 : 座標 )
  size : 画像のサイズ
  r : 極値判定をするときの周辺画素の範囲
*/
void CornerSieve( Corner* corner, const Coord< int >& size, int r )
{
  // 座標をキーにして登録
  typedef std::map< Coord< int >, double, LessCoord > Buffer;
  Buffer buffer;
  for ( Corner::const_iterator i = corner->begin() ; i != corner->end() ; ++i )
    buffer[i->second] = i->first;

  corner->clear();
  for ( Buffer::const_iterator i = buffer.begin() ; i != buffer.end() ; ++i ) {
    const Coord< int >& c = i->first;
    bool isPeak = true;
    for ( int y = -r ; y < r ; ++y ) {
      for ( int x = -r ; x < r ; ++x ) {
        if ( x * x + y * y > r * r ) continue;
        Buffer::const_iterator j = buffer.find( Coord< int >( c.x + x, c.y + y ) );
        if ( j == buffer.end() ) continue;
        if ( i->second < j->second ) {
          isPeak = false;
          break;
        }
      }
      if ( ! isPeak ) break;
    }
    if ( isPeak )
      corner->insert( std::make_pair( i->second, i->first ) );
  }
}

処理の内容は非常に単純で、corner のキーと座標を反転して buffer に登録し、半径 r 以内に検出点があったら値を比較します。もし、周囲の検出点のほうが値が大きければ対象の点は削除し、そのような点が一つもない場合だけ corner に再登録します。

Coord< int > 型は大小比較ができないため、比較用の関数オブジェクトとして LessCoord 構造体を定義しています。したがって、buffer の型は比較用の関数オブジェクトもテンプレート引数の中に含めてあることに注意して下さい。


2) ハリス・コーナー検出器 ( Harris Corner Detector )

画素のある領域を ( u, v ) だけずらしたときの強度の変化を、下式のように表現します。

E( u, v ) = Σx,yw( x, y )[ I( x + u, y + v ) - I( x, y ) ]2

w( x, y ) は窓関数で、次式のような矩形窓

w( x, y )=1( 0 ≤ x ≤ 1 ; 0 ≤ y ≤ 1 )
=0( それ以外 )

や、ガウス窓

w( x, y ) = exp( -( x2 + y2 ) / 2σ2 )

などが用いられ、画像のある有限な領域のみを演算の対象とする役割があります。

ここで、ニ変数でのテーラー-マクローリン展開 ( *2-1 )

f( x + u, y + v ) = Σn{0→∞}( [ u( ∂ / ∂x ) + v( ∂ / ∂y ) ]nf( x, y ) / n! )

の一次までの近似式

f( x + u, y + v ) ≅ f( x, y ) + u( ∂f / ∂x ) + v( ∂f / ∂y )

を使って E( u, v ) を計算すると、

E( u, v )Σx,y( w( x, y )[ I( x, y ) + uIx + vIy - I( x, y ) ]2 )
=Σx,y( w( x, y )( uIx + vIy )2 )
=Σx,y( w( x, y )( u2Ix2 + 2uvIxIy + v2Iy2 ) )
=u2Σx,y( w( x, y )Ix2 ) + 2uvΣx,y( w( x, y )IxIy ) + v2Σx,y( w( x, y )Iy2 )
=| u, v ||Σx,y( w( x, y )Ix2 ), Σx,y( w( x, y )IxIy ) || u |
|Σx,y( w( x, y )IxIy ), Σx,y( w( x, y )Iy2 ) || v |
aTMa

となります。但し、

a = ( u, v )T

M=|Σx,y( w( x, y )Ix2 ), Σx,y( w( x, y )IxIy )|
|Σx,y( w( x, y )IxIy ), Σx,y( w( x, y )Iy2 )|

です。

窓関数 w( x, y ) が矩形窓であるならば、行列 M の要素はある有限な領域について和をとった値になります。下図のような画像に対して、a から d までの 4 つの矩形領域で和を計算した結果を考えてみます。

図 2-1. 行列 M とコーナー検出の関係
行列Mとコーナー検出の関係

a の領域では、Ix2 と Iy2 に関する和はゼロより十分大きいある値になることが期待できますが、IxIy については Ix > 0 のとき Iy = 0、Iy > 0 のとき Ix = 0 となるので、その和はゼロです。

よって、

M=|Σx,y( Ix2 ),0|
|0,Σx,y( Iy2 )|

という結果が得られます。b の領域の場合、Ix = 0 ですから、

M=|0,0|
|0,Σx,y( Iy2 )|

c の領域の場合、Iy = 0 より

M=|Σx,y( Ix2 ),0|
|0,0|

となります。d の領域では Ix = Iy = 0 なので、M = 0 になります。

したがって、M の対角成分が大きな値であれば、コーナーであると判断することができることになります。

ところが、コーナー部分は常に矩形に対して整列しているとは限らず、斜め方向に傾いている場合もあるでしょう。そのため、行列 M を対角化 ( つまり回転 ) する必要があります ( *2-2 )。

対角化によって、M は固有値を対角成分に持つ対角行列 D と固有ベクトルからなる行列 U を使って

M = UDUT=U|λ1,0| UT
|0,λ2|

となります。この固有値 λ1 と λ2 を使ってコーナー検出をすればいいことになります。

しかし、小領域ごとに毎回固有値を計算するのはあまり効率的ではないため、固有値を求めずにコーナーを判定する手法があります。それを「ハリス・コーナー検出器 ( Harris Corner Detector )」といいます。

ハリス・コーナー検出器では、次の「コーナネス ( Cornerness )」と呼ばれる値 C を計算します。

C = det( M ) - k[ tr( M ) ]2

但し、det( M ) は行列 M の行列式、tr( M ) は行列 M のトレース、k は定数 ( 通常は 0.04 ~ 0.06 ) を表します。固有値 λ1, λ2 に対して

det( M ) = λ1λ2

tr( M ) = λ1 + λ2

なので ( 補足 3 )、C の値が大きければコーナーであると判定することができます。

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

*2-2)固有値問題 (1) 対称行列の固有値」を参照


微分値を利用したコーナー検出のサンプル・プログラムを以下に示します。

/*
  Mirror : 座標 c0 が境界 sz 外だった場合、鏡像 ( 折り返し ) の座標を返す

  例)

  c0 = -3 のとき、3 を返す。但し、sz <= 3 ならば再び境界外となるので、そのときは sz - 1 を返す。
  sz = 3, c0 = 4 のとき、( sz - 1 ) x 2 - c0 = 0 を返す。c0 = 5 のときは計算結果が負数になるので、この場合はゼロを返す。
*/
inline int Mirror( int c0, int sz )
{
  if ( c0 < 0 ) {
    c0 = -c0;
    if ( c0 >= sz ) c0 = sz - 1;
  } else if ( c0 >= sz ) {
    c0 = ( sz - 1 ) * 2 - c0;
    if ( c0 < 0 ) c0 = 0;
  }

  return( c0 );
}

/*
  Move : 座標 c0 から d 進んだ座標値を返す

  結果が領域外だった場合は Mirror を使って鏡像 ( 折り返し ) の座標を返す。
*/
inline int Move( int c0, int d, int sz )
{
  c0 += d;

  return( Mirror( c0, sz ) );
}

/*
  CalcMatrix : 対角成分を Σix^2, Σiy^2、非対角成分を Σix * iy とする 2 x 2 の対称行列を求める

  和の成分は窓関数 w によって重み付けされる。和を求める範囲は w の大きさで決まる。
  w の大きさは偶数でも算出できるが、中央の座標はずれた状態になる。

  ix, iy : 算出に使うデータ
  c : 算出対象の座標値
  w : 窓関数
  m00, m01, m11 : 求める行列の成分を返す変数へのポインタ ( m00, m11 : 対角成分、m01 ; 非対角成分 )
*/
void CalcMatrix( const Matrix< double >& ix, const Matrix< double >& iy, Coord< int > c,
                 const Matrix< double >& w, double* m00, double* m01, double* m11 )
{
  *m00 = *m01 = *m11 = 0;

  if ( w.rows() < 1 || w.cols() < 1 ) return;

  // 左上隅の相対座標
  int x0 = -( ( w.cols() - 1 ) / 2 );
  int y0 = -( ( w.rows() - 1 ) / 2 );

  Coord< int > dc;
  for ( dc.y = y0 ; dc.y < static_cast< int >( w.rows() ) + y0 ; ++( dc.y ) ) {
    for ( dc.x = x0 ; dc.x < static_cast< int >( w.cols() ) + x0 ; ++( dc.x ) ) {
      std::pair< double, double > i( ix[Move( c.y, dc.y, ix.rows() )][Move( c.x, dc.x, ix.cols() )],
                                     iy[Move( c.y, dc.y, ix.rows() )][Move( c.x, dc.x, ix.cols() )] );
      *m00 += std::pow( i.first, 2 ) * w[dc.y - y0][dc.x - x0];
      *m01 += i.first * i.second * w[dc.y - y0][dc.x - x0];
      *m11 += std::pow( i.second, 2 ) * w[dc.y - y0][dc.x - x0];
    }
  }
}

/*
  DecisionByEigenvalue : 固有値を計算して判定する関数オブジェクト
*/
class DecisionByEigenvalue
{
  double threshold_; // 判定用のしきい値

public:

  /*
    判定用のしきい値を指定して構築
  */
  DecisionByEigenvalue( double threshold )
    : threshold_( threshold ) {}

  /*
    operator() : 判定用ルーチン

    m00, m01, m11 : 行列の成分 ( m00, m11 : 対角成分、m01 ; 非対角成分 )
    cornerness : 固有値の積を返す変数へのポインタ ( 判定結果が false の場合は代入されない )
  */
  bool operator()( double m00, double m01, double m11, double* cornerness )
  {
    double d = std::sqrt( std::pow( m00 + m11, 2 ) - 4 * ( m00 * m11 - std::pow( m01, 2 ) ) );
    double ev1 = ( m00 + m11 + d ) / 2;
    double ev2 = ( m00 + m11 - d ) / 2;

    if ( ev1 < threshold_ || ev2 < threshold_ ) return( false );

    *cornerness = ev1 * ev2;

    return( true );
  }
};

/*
  DecisionByCornerness : コーナネスを計算して判定する関数オブジェクト
*/
class DecisionByCornerness
{
  double threshold_; // 判定用のしきい値
  double k_;         // トレースに掛ける定数

public:

  /*
    判定用のしきい値とトレースに掛ける定数を指定して構築
  */
  DecisionByCornerness( double threshold, double k )
    : threshold_( threshold ), k_( k ) {}

  /*
    operator() : 判定用ルーチン

    m00, m01, m11 : 行列の成分 ( m00, m11 : 対角成分、m01 ; 非対角成分 )
    cornerness : コーナネスを返す変数へのポインタ
  */
  bool operator()( double m00, double m01, double m11, double* cornerness )
  {
    *cornerness = m00 * m11 - m01 * m01 - k_ * std::pow( m00 + m11, 2 );

    return( *cornerness >= threshold_ );
  }
};

/*
  Corner_Defferential : 微分値を利用したコーナ判定

  data : 画素データ
  w : 窓関数
  corner : 検出したコーナーの座標値を保持する multimap ( キー : 極値 ; 値 : 座標 )
  op : 判定用関数オブジェクト
*/
template< class Op >
void Corner_Defferential( const Matrix< double >& data, const Matrix< double >& w,
                          Corner* corner, Op op )
{
  corner->clear();

  // サイズが 3 より小さい画像に対しては検出しない
  if ( data.rows() < 3 || data.cols() < 3 ) return;

  // 微分値
  Matrix< double > ix( data.rows() - 1, data.cols() - 1 );
  Matrix< double > iy( data.rows() - 1, data.cols() - 1 );

  // 微分値の計算
  Coord< int > c;
  for ( c.y = 1 ; c.y < static_cast< int >( data.rows() ) ; ++( c.y ) ) {
    for ( c.x = 1 ; c.x < static_cast< int >( data.cols() ) ; ++( c.x ) ) {
      ix[c.y - 1][c.x - 1] = data[c.y][c.x] - data[c.y][c.x - 1];
      iy[c.y - 1][c.x - 1] = data[c.y][c.x] - data[c.y - 1][c.x];
    }
  }

  for ( c.y = 0 ; c.y < static_cast< int >( ix.rows() ) ; ++( c.y ) ) {
    for ( c.x = 0 ; c.x < static_cast< int >( ix.cols() ) ; ++( c.x ) ) {
      double m00, m01, m11, cornerness;
      CalcMatrix( ix, iy, c, w, &m00, &m01, &m11 );
      if ( op( m00, m01, m11, &cornerness ) )
        corner->insert( std::pair< double, Coord< int > >( cornerness, Coord< int >( c.x + 1, c.y + 1 ) ) );
    }
  }
}

/*
  Corner_Eigenvalue : 固有値を利用したコーナ判定

  data : 画素データ
  w : 窓関数
  corner : 検出したコーナーの座標値を保持する multimap ( キー : 極値 ; 値 : 座標 )
  threshold : 判定用のしきい値
*/
void Corner_Eigenvalue( const Matrix< double >& data, const Matrix< double >& w,
                        Corner* corner, double threshold )
{
  DecisionByEigenvalue op( threshold );

  Corner_Defferential( data, w, corner, op );
}

/*
  Corner_Harris : ハリスコーナ検出器

  data : 画素データ
  w : 窓関数
  corner : 検出したコーナーの座標値を保持する multimap ( キー : 極値 ; 値 : 座標 )
  threshold : 判定用のしきい値
  k :トレースに掛ける定数
*/
void Corner_Harris( const Matrix< double >& data, const Matrix< double >& w,
                    Corner* corner, double threshold, double k )
{
  DecisionByCornerness op( threshold, k );

  Corner_Defferential( data, w, corner, op );
}

CalcMatrix は、あらかじめ計算した微分値 ixiy を使い、指定した座標 c での行列 M の成分を求める関数です。窓関数 w はあらかじめ計算して矩形 ( Matrix 型 ) の形で渡されるようにしてあります。処理内容は、窓関数の領域にある画素値の和を求めるだけの単純なものです。
この関数の中で、インライン関数として Move が呼び出されていますが、これは第一引数の値を第二引数分だけ進めた値を返す関数です。Move の中で、さらに Mirror というインライン関数が呼び出されていますが、これは座標値が範囲外になった場合に鏡像、つまり範囲外の画像がちょうど鏡に映したように反転した状態にあるとして値を返すような仕組みになっています。

コーナーかどうかを判定する関数オブジェクトが二つ用意されています。一つは DecisionByEigenvalue で、固有値を直接求めて判定する形をとります。構築時にしきい値を渡し、二つの固有値がどちらもしきい値以上であれば、true を返す叙述関数 operator() が実装されています。また、引数の cornerness には二つの固有値の積が返されます。もう一つは DecisionByCornerness で、ハリス・コーナー検出器で使われるコーナネスを使って判定する形をとります。構築時にはしきい値の他に定数 k を渡し、叙述関数 operator() でコーナネスを計算してしきい値以上なら true を返します。

Corner_Defferential がコーナー判定のためのメイン・ルーチンであり、テンプレート引数の OpDecisionByEigenvalueDecisionByCornerness のいずれかを表すことになります。最初に微分値を計算し、行列 M の要素を CalcMatrix で計算したら、DecisionByEigenvalueDecisionByCornerness のいずれかを使ってコーナー判定を行い、コーナーだと判定された場合は corner に登録を行います。


3) セグメント・テスト ( Segment Test )

コーナーかどうかを判定する画素の周囲の状態を調べてコーナーを検出する方法を検討してみます。具体的には、対象の画素を中心とする円形の領域について、円周上にある画素が対象の画素と近い値を持つかどうかを調べます。

図 3-1. 周囲の画素による判定方法
周囲の画素による判定方法

上図において、a, b, c それぞれの位置で円周上の画素を調べたとき、赤い部分は中心の画素と異なり、青い部分は中心の画素と一致していることがわかります。a の場合、赤い部分は全体の中の 3 / 4 程度を占めており、しかも連続しています。b の場合も赤い部分は連続していますが、全体の半分程度となっています。c の場合は赤い部分はなく、全てが青い部分で占められています。
このことから、コーナーを検出するためには、赤い部分が連続して 3 / 4 以上を占めている場所を見つければよいことがわかります。

そこで、ある画素 c を中心として決められた半径 r 上にある画素について、c よりしきい値以上大きい ( または小さい ) 数をカウントして、それが定数 n 個以上であれば c をコーナーと判定します。この手法を「セグメント・テスト ( Segment Test )」といいます。例えば下図にあるような、r = 3 の場合を考えてみましょう。

図 3-2. セグメント・テスト ( r = 3 )
セグメント・テスト(r=3)

中心の画素 c の周囲には 16 個の画素があり、それぞれに 1 から 16 までの番号が割り当てられています。あるしきい値 t ( > 0 ) を使い、これらの画素が連続して 12 個以上、c より t 以上大きいか、または t 以上小さくなっていた場合、c をコーナーとみなします。その判定は、次のようにすることで高速に行うことができます。

  1. 画素 1 と画素 9 を調べる。両方とも、画素 c との差の絶対値がしきい値 t 未満なら c をコーナーと見なさず処理を終了する。
  2. 画素 5 と画素 13 を調べる。画素 1, 5, 9, 13 のうち、少なくとも 3 個の画素が t 以上大きいか、t 以上小さい場合は処理を続け、そうでなければ c をコーナーと見なさず処理を終了する。
  3. 残り 12 個の画素を調べ、連続して 12 個以上の画素が、c より t 以上大きいか、または t 以上小さくなっていた場合、c をコーナーとみなす。

r = 3、n = 12 の場合、上記の方法によって比較的高速に処理を行うことができますが、r や n の条件を変えたい場合はこの手法は利用できません。また、最初に調べる画素を 1, 9、次に調べる画素を 5, 13 としていますが、棄却するために調べる画素がこれで最適かどうかはわからないという問題もあります。

セグメント・テストや後述する FAST コーナー検出器では、真のコーナーに近い画素についてもコーナーと判定される可能性があるため、それらの画素の中で最も指標の高い点をコーナーと点として抽出する必要があります。しかし、セグメント・テストや FAST コーナー検出器の場合、指標となる値がありません。そのため、次式のような式を使って指標 V を計算して極大値を求めるという手法が用いられます。

V = max( Σi{i∈Cb}( | Ii - Ic | - t ), Σi{i∈Cd}( | Ic - Ii | - t ) )

セグメント・テストを使ったコーナー検出のサンプル・プログラムを以下に示します。

/*
  IsDarker : 差分 d がしきい値 t より小さい(暗い)かどうかを判定する関数オブジェクト
*/
struct IsDarker
{
  bool operator()( double d, double t )
  { return( d < -t ); }
};

/*
  IsBrighter : 差分 d がしきい値 t より大きい(明るい)かどうかを判定する関数オブジェクト
*/
struct IsBrighter
{
  bool operator()( double d, double t )
  { return( d > t ); }
};

/*
  DecisionSegment : コーナ点の連鎖があるかどうかを判定する

  data : 画素データ
  c : コーナかどうかを判定する座標
  arc : c に対する相対位置を保持する配列の先頭
  szArc : c に対する相対位置を保持する配列の長さ
  op : 中央より明るい(暗い)を判定する関数オブジェクト
  threshold : op による判定時のしきい値

  戻り値 : 連鎖があれば true
*/
template< class Op >
bool DecisionSegment( const Matrix< double >& data, Coord< int > c,
                      const Coord< int >* arc, size_t szArc,
                      Op op, double threshold )
{
  double p = data[c.y][c.x]; // 中央の画素データ
  const size_t cornerLen = szArc * 3 / 4; // 連鎖の最小の長さ

  // 先頭から順にコーナ点でない場所を探索
  size_t s = 0;
  for ( ; s < szArc ; ++s )
    if ( ! op( p - data[c.y + arc[s].y][c.x + arc[s].x], threshold ) )
      break;
  if ( s >= szArc ) return( true );

  // 末尾から順にコーナ点でない場所を探索
  size_t e = szArc;
  for ( ; e > 0 ; --e )
    if ( ! op( p - data[c.y + arc[e - 1].y][c.x + arc[e - 1].x], threshold ) )
      break;

  // 先頭(末尾)をはさんでコーナと判定できるかチェック
  if ( s + ( szArc - e ) >= cornerLen ) return( true );

  // 残りの部分についてチェック
  while ( e > s + cornerLen + 2 ) {
    // コーナ点の開始位置を検出
    size_t i = s + 1;
    for ( ; i < e - 1 ; ++i )
      if ( op( p - data[c.y + arc[i].y][c.x + arc[i].x], threshold ) )
        break;
    if ( i == e - 1 ) return( false );
    // コーナ点の連鎖の末尾の次の位置まで探索
    for ( s = i + 1 ; s < e - 1 ; ++s )
      if ( ! op( p - data[c.y + arc[s].y][c.x + arc[s].x], threshold ) )
        break;
    if ( s - i >= cornerLen ) return( true );
  }

  return( false );
}

/*
  DiffSum : 中央より暗い(明るい)点の差分値 - しきい値の和を求める

  data : 画素データ
  c : コーナかどうかを判定する座標
  arc : c に対する相対位置を保持する配列の先頭
  szArc : c に対する相対位置を保持する配列の長さ
  op : 中心より明るい(暗い)を判定する関数オブジェクト
  threshold : op による判定時のしきい値
*/
template< class Op >
double DiffSum( const Matrix< double >& data, Coord< int > c,
                const Coord< int >* arc, size_t szArc,
                Op op, double threshold )
{
  double p = data[c.y][c.x]; // 中央の値

  // 中央との差分としきい値との差を計算して和を求める
  double sum = 0;
  for ( const Coord< int >* a = arc ; a < arc + szArc ; ++a ) {
    double diff = p - data[c.y + a->y][c.x + a->x];
    if ( ! op( diff, threshold ) )
      continue;
    sum += std::abs( diff ) - threshold;
  }

  return( sum );
}

// 周囲のピクセルの相対位置
const Coord< int > SEGMENT[] = {
  Coord< int >(  0, -3 ), Coord< int >(  1, -3 ), Coord< int >(  2, -2 ), Coord< int >(  3, -1 ),
  Coord< int >(  3,  0 ), Coord< int >(  3,  1 ), Coord< int >(  2,  2 ), Coord< int >(  1,  3 ),
  Coord< int >(  0,  3 ), Coord< int >( -1,  3 ), Coord< int >( -2,  2 ), Coord< int >( -3,  1 ),
  Coord< int >( -3,  0 ), Coord< int >( -3, -1 ), Coord< int >( -2, -2 ), Coord< int >( -1, -3 ),
};
const size_t SIZE_OF_SEGMENT = sizeof( SEGMENT ) / sizeof( SEGMENT[0] );

/*
  SegmentTest : セグメント・テストによるコーナ判定

  data : 画素データ
  c : コーナかどうかを判定する座標
  threshold : コーナ点判定時のしきい値
*/
bool SegmentTest( const Matrix< double >& data, Coord< int > c, double threshold )
{
  double d[4]; // 周囲(上下左右)のピクセル

  double p = data[c.y][c.x]; // 中央の値
  // 上下のピクセルがどちらも暗くも明るくもなければスキップ
  d[0] = p - data[c.y + SEGMENT[0].y][c.x + SEGMENT[0].x];
  d[2] = p - data[c.y + SEGMENT[8].y][c.x + SEGMENT[8].x];
  if ( std::abs( d[0] ) < threshold && std::abs( d[2] ) < threshold ) return( false );
  // 左右のピクセルを含めて暗い・明るいピクセルの数をカウント
  d[1] = p - data[c.y + SEGMENT[4].y][c.x + SEGMENT[4].x];
  d[3] = p - data[c.y + SEGMENT[12].y][c.x + SEGMENT[12].x];
  unsigned int darker = 0;
  unsigned int brighter = 0;
  for ( size_t i = 0 ; i < 4 ; ++i ) {
    if ( d[i] < -threshold ) {
      ++darker;
    } else if ( d[i] > threshold ) {
      ++brighter;
    }
  }
  // 少なくとも 3 つ以上同じピクセルがなければ連鎖は生じない
  if ( darker < 3 && brighter < 3 ) return( false );
  // 暗い・明るい場合は全ての周囲ピクセルをチェック
  if ( darker >= 3 ) {
    return( DecisionSegment( data, c, &SEGMENT[0], SIZE_OF_SEGMENT, IsDarker(), threshold ) );
  } else {
    return( DecisionSegment( data, c, &SEGMENT[0], SIZE_OF_SEGMENT, IsBrighter(), threshold ) );
  }
}

/*
  Corner_SegmentTest : セグメント・テストによるコーナ検出

  data : 画素データ
  corner : 検出したコーナーの座標値を保持する multimap ( キー : 極値 ; 値 : 座標 )
  threshold : 判定用のしきい値
*/
void Corner_SegmentTest( const Matrix< double >& data, Corner* corner, double threshold )
{
  corner->clear();
  if ( data.rows() < 7 || data.cols() < 7 ) return;

  Coord< int > c;
  for ( c.y = 3 ; c.y < static_cast< int >( data.rows() - 3 ) ; ++( c.y )  ) {
    for ( c.x = 3 ; c.x < static_cast< int >( data.cols() - 3 ) ; ++( c.x ) ) {
      if ( SegmentTest( data, c, threshold ) )
        corner->insert( std::pair< double, Coord< int > >
                        ( std::max( DiffSum( data, c, &SEGMENT[0], SIZE_OF_SEGMENT, IsDarker(), threshold ),
                                    DiffSum( data, c, &SEGMENT[0], SIZE_OF_SEGMENT, IsBrighter(), threshold ) ),
                          c ) );
    }
  }
}

メイン・ルーチンは Corner_SegmentTest ですが、判定処理は SegmentTest が行っています。まず、画素 1, 9 ( 配列は 0 からスタートするので、プログラム内では 0, 8 になります ) の中心画素との差の絶対値がしきい値 threshold より小さければ、false を返して処理を終了します。次に画素 5, 13 の中心画素との差を求め、-threshold より小さい個数と threshold より大きい個数を数えます。その個数がどちらも 3 より小さければ、false を返して終了します。これらのテストを通過した場合、DecisionSegment を呼び出して全ての画素について判定を行います。このとき、画素 1, 5, 9, 13 の中の 3 つ以上の画素が中心画素より大きい場合 ( brighter ≥ 3 ) と小さい場合 ( darker ≥ 3 ) で、判定用の関数オブジェクトを切り替えていることに注意して下さい。

関数 DecisionSegment では、先頭からと末尾からのそれぞれで、中央の画素より値が大きくない ( または小さくない ) 場所を探索します。実際には先頭と末尾はつながっているので、その部分をはさんでコーナーと判定できるかをまずチェックします。もし、コーナーでなければ、あとは残りの部分について連鎖がないかを確認します。

関数 DiffSum は、コーナーとして検出された点において、どの程度コーナーらしいかを判断する指標 V を計算するために使われます。


4) FAST コーナー検出器 ( FAST Corner Detector )

セグメント・テストでは、中心画素より値がしきい値以上大きい ( または小さい ) 連鎖を見つける処理が必要でした。この処理が、処理を重くする要因となっています。連鎖を全てチェックするのではなく、一部の限られた画素だけをチェックしてコーナーかどうかを判定することができれば、処理時間を短くすることができます。この手法を「FAST コーナー検出器 ( FAST Corner Detector )」といい、コーナー判定に「決定木 ( Decision Tree )」を用いるという特徴があります。

まず、すでに決定木があるとして、コーナー判定をどのように行うかを説明します。

中心画素 c の周囲の画素の番号を i ( 1 ≤ i ≤ 16 ) として、それぞれの画素の値を Ic、Ii とします。次に、周囲の画素の状態 Si を、Ic と Ii、またしきい値 t を用いて、三つの状態に分類します。

Si=d( Ii ≤ Ic - t )
=s( Ic - t < Ii < Ic + t )
=b( Ii ≥ Ic + t )

つまり、周囲の画素が中心画素より t 以上小さければ d、中心画素より t 以上大きければ b、いずれでもなければ s とするということになります。

下図には、決定木と、周囲の画素 1 から 16 までを横に並べた例を示しています。黒色が d、白色が b、灰色が s の状態であることを表しています。

図 4-1. 決定木を使ったコーナー検出
決定木を使ったコーナー検出

決定木が上図のような状態だった場合、まず最初に決定木の最上位にある 8 番目の画素をチェックします。それは d なので、決定木を d の方向に進むと、次は 2 番目の画素になります。それも d なので、その方向に進みます。今度は 7 番目の画素になります。これも d なので、その方向に進むと、それ以上分岐しない葉の節に到達し、この例ではコーナーであると判定されます。もし、7 番目の画素が b だったとしたら、非コーナーと判定されることになります。

決定木の学習には、あらかじめコーナー点と非コーナー点に分類された画像の集合を使って「ID3 アルゴリズム」によって行われます。

まず、学習に使う全画像の持つ中心画素の集合 C に対し、番号 i の画素について、Si が d, s, b のいずれかであるかによって集合 Cd, Cs, Cb に分割します。集合 C について、「エントロピー ( Entropy )」H[C] を次の式で計算します。

H[C] = -[ c1 / ( c1 + c2 ) ]log[ c1 / ( c1 + c2 ) ] - [ c2 / ( c1 + c2 ) ]log[ c2 / ( c1 + c2 ) ]

ここで c1 は集合 C の中のコーナーと判定された画素の数、c2 は非コーナーと判定された画素の数をそれぞれ表します。

Cd, Cs, Cb についても同様にエントロピー H[Cd], H[Cs], H[Cb] を計算し、最後に画素 i に対する「情報利得 ( Information Gain )」Hi を下式で求めます。

Hi = H[C] - ( |Cb| / |C| )H[Cb] - ( |Cs| / |C| )H[Cs] - ( |Cd| / |C| )H[Cd]

各番号 i の情報利得 Hi を求めたら、その中で Hi が最大になる i を選択し、C をその番号で Cd, Cs, Cb に分割して子ノードとします。あとは、子ノードに対して同様の処理を繰り返していきます。


FAST コーナー検出器を使ったコーナー検出のサンプル・プログラムを以下に示します。

/*
  PixelState : 周囲ピクセルの状態(中央ピクセルより暗い・明るい・同じ)
*/
enum PixelState {
  DARK,   // 暗い
  SAME,   // 同じ
  BRIGHT, // 明るい
};

/*
   FASTアルゴリズムによるコーナー検出用クラス
*/
class FAST
{
  ID3_Algorithm< PixelState, bool > id3_; // ID3 アルゴリズム(決定木)

 public:

  /*
    FASTコンストラクタ : 学習パターン(セグメント・テストの結果) x と教師信号(コーナーであるか) y を指定して構築
  */
  FAST( const Matrix< PixelState >& x, const std::vector< bool >& y )
    : id3_( x, y ) {}

  /*
    FAST::decision : data の座標 c がコーナーかどうかを判定する

    コーナーなら 1、非コーナーなら 0 を返す。セグメントが決定木になかった場合は - 1 を返す。
  */
  int decision( const Matrix< double >& data, Coord< int > c, double threshold ) const;
};

/*
  FAST::decision : 画素データ data の c の位置がコーナーかどうかを判定する

  data : 画素データ
  c : コーナかどうかを判定する座標
  threshold : コーナ点判定時のしきい値
*/
int FAST::decision( const Matrix< double >& data, Coord< int > c, double threshold ) const
{
  PixelState state[SIZE_OF_SEGMENT]; // セグメントの状態を保持する配列

  double p = data[c.y][c.x];
  for ( size_t i = 0 ; i < SIZE_OF_SEGMENT ; ++i ) {
    double diff = p - data[c.y + SEGMENT[i].y][c.x + SEGMENT[i].x];
    if ( diff < -threshold )
      state[i] = DARK;
    else if ( diff > threshold )
      state[i] = BRIGHT;
    else
      state[i] = SAME;
  }

  std::pair< bool, bool > ans = id3_.predict( &state[0] );

  // state と同じパターンが決定木に存在しなければ負数を返す
  if ( ! ans.second )
    return( -1 );

  return( ( ans.first ) ? 1 : 0 );
}

/*
  FAST_Accumulate : FASTアルゴリズムを利用するための教師データ x, y を画素データ data から取得する

  data : 画素データ
  x : 学習パターン(周囲ピクセルの状態)を保持する変数へのポインタ
  y : 教師ベクトル(コーナ・非コーナの判定)を保持する変数へのポインタ
  threshold : 判定用のしきい値
*/
void FAST_Accumulate( const Matrix< double >& data, Matrix< PixelState >* x, std::vector< bool >* y, double threshold )
{
  assert( x->cols() == SIZE_OF_SEGMENT || x->rows() == 0 );
  if ( x->rows() > 0 ) {
    ErrLib::CheckLinearModel( *x, y->begin(), y->end() );
  } else {
    x->assign( 0, SIZE_OF_SEGMENT );
    y->clear();
  }

  if ( data.rows() < 7 || data.cols() < 7 ) return;
  size_t s = x->rows(); // 書き込む行の開始位置
  x->resize( s + ( data.rows() - 6 ) * ( data.cols() - 6 ), x->cols() );
  y->resize( s + ( data.rows() - 6 ) * ( data.cols() - 6 ) );
  Coord< int > c;
  for ( c.y = 3 ; c.y < static_cast< int >( data.rows() - 3 ) ; ++( c.y )  ) {
    for ( c.x = 3 ; c.x < static_cast< int >( data.cols() - 3 ) ; ++( c.x ) ) {
      double p = data[c.y][c.x];
      for ( size_t i = 0 ; i < SIZE_OF_SEGMENT ; ++i ) {
        double diff = p - data[c.y + SEGMENT[i].y][c.x + SEGMENT[i].x];
        if ( diff < -threshold )
          (*x)[s][i] = DARK;
        else if ( diff > threshold )
          (*x)[s][i] = BRIGHT;
        else
          (*x)[s][i] = SAME;
      }
      (*y)[s] = SegmentTest( data, c, threshold );
      ++s;
    }
  }
}

/*
  FAST_Detection : FAST アルゴリズムによって画素データ data のコーナーを検出する

  data : 画素データ
  fast : FAST アルゴリズム
  corner : 検出したコーナーの座標値を保持する multimap ( キー : 極値 ; 値 : 座標 )
  threshold : 判定用のしきい値
*/
void FAST_Detection( const Matrix< double >& data, const FAST& fast, Corner* corner, double threshold )
{
  corner->clear();
  if ( data.rows() < 7 || data.cols() < 7 ) return;

  Coord< int > c;
  for ( c.y = 3 ; c.y < static_cast< int >( data.rows() - 3 ) ; ++( c.y ) ) {
    for ( c.x = 3 ; c.x < static_cast< int >( data.cols() - 3 ) ; ++( c.x ) ) {
      if ( fast.decision( data, c, threshold ) )
        corner->insert( std::pair< double, Coord< int > >
                        ( std::max( DiffSum( data, c, &SEGMENT[0], SIZE_OF_SEGMENT, IsDarker(), threshold ),
                                    DiffSum( data, c, &SEGMENT[0], SIZE_OF_SEGMENT, IsBrighter(), threshold ) ),
                          c ) );
    }
  }
}

まず、学習用の教師データを画像から作成する必要がありますが、その処理を行うための関数が FAST_Accumulate です。最初に呼び出されている ErrLib::CheckLinearModel は、Matrix 型の引数 x の行数と vector 型の引数 y のサイズが一致しているかをチェックするための関数ですが、この中では実装していません。処理の内容は単純で、画素 data の周囲の画素が threshold 以上に大きい、または小さいかを判定して DARK、BRIGHT、SAME の 3 つの状態に振り分けて x に代入し、SegmentTest 関数を使ってコーナーか非コーナーかを判定してその結果を y に登録します。

教師データが得られたら、FAST オブジェクトを構築します。FAST は決定木の「ID3 アルゴリズム」である id3_ を内部に保持しており、教師データ x, y を引数として渡せば決定木が得られるようになっています。

決定木が得られたら、あとは FAST_Detection を呼び出してコーナーを検出するだけです。実際のコーナー検出は FAST クラスのメンバ関数 decision が行います。まず、周囲の画素の状態を調べて DARK、BRIGHT、SAME の 3 つの状態に振り分けた後、決定木からコーナー・非コーナーを判定します。


5) 性能評価

最後に、サンプル・プログラムを使い、テスト画像を使ってコーナー検出を実行した結果を紹介します。どのテストも、画像の輝度を抽出してコーナーを検出しています。

図 5-1. テスト画像 (1)
テスト画像

下図は、上に示した「テスト画像 (1)」のような、図形を並べた画像に対してコーナー検出をした結果で、円で示された部分がコーナーと判定された箇所になります。ヘシアンを利用した検出、固有値を利用した検出、ハリス・コーナー検出器については、σ = 1 のガウシアン・フィルタを行ってからコーナー検出処理を行っています。また、固有値を利用した検出とハリス・コーナー検出器では、窓関数として 5 x 5 の矩形窓を利用しています。さらに、抽出したコーナーに対して、CornerSieve 関数を使い、r = 5 の範囲から極大値だけを抽出する処理を行っています。

図 5-2. コーナー検出結果
ヘシアンを利用したコーナー検出 ( threshold=1000 )
Hessian Corner Detector
固有値を利用したコーナー検出 ( threshold=1000 )
Corner Detector by Eigenvalue
ハリス・コーナー検出器 ( k = 0.04, threshold=1000 )
Harris Corner Detector
セグメント・テスト ( threshold=10 )
Segment Test
FASTコーナー検出器 ( threshold=10 )
FAST Corner Detector

ヘシアンを利用したコーナー検出は、コーナーが検出できていないだけでなく、エッジ部分が誤判定されるなど、正常に検出できていないことがわかります。固有値を利用したコーナー検出の場合、コーナーは検出されているものの、エッジ部分を誤判定している箇所はまだ多くあります。
5 つのコーナー検出器の中で最も正しく検出できたのがハリス・コーナー検出器です。ラインを誤判定したり、円に対してコーナーと検出しているところもありますが、固有値を利用したコーナー検出よりも誤判定が少ないという結果になりました。
セグメント・テストの場合、エッジやラインに対する誤判定はないですが、コーナーとして判定されない箇所がいくつかあります。また、FAST コーナー検出器はセグメント・テストと全く同じ結果となっています。

図 5-3. テスト画像 (2)
テスト画像

次は、自然画に近いグレースケール画像に対してコーナー検出を行った結果を示します。テスト条件は、先程の画像の場合と同一です。

図 5-2. コーナー検出結果
ヘシアンを利用したコーナー検出 ( threshold=50 )
Hessian Corner Detector
固有値を利用したコーナー検出 ( threshold=500 ) ハリス・コーナー検出器 ( k = 0.04, threshold=100000 )
Corner Detector by Eigenvalue Harris Corner Detector
セグメント・テスト ( threshold=14 ) FASTコーナー検出器 ( threshold=14 )
Segment Test FAST Corner Detector

このテストでも、固有値を利用したコーナー検出とハリス・コーナー検出器の結果が最も良好だという結果になりました。


今回は、コーナー検出に関するいろいろなアルゴリズムについて紹介しました。画像認識の分野では、抽出したコーナーを使って特徴量を計算し、画像の中にある物体を検出したり、同じ特徴の画像を選択するような処理が行われるため、コーナー検出の精度が非常に重要なものになります。これからも、より精度の高いコーナー検出法が誕生するかもしれませんね。


補足 1) ガウス曲率 ( Gaussian Curvature )

まずは、平面上の曲線に対して「曲率 ( Curvature )」を定義します。

平面上の曲線は、パラメータ t を使って次のように表すことができます。

p(t) = ( x(t), y(t) )

例えば、

p(t) = ( x(t), y(t) ) = ( a cos t, a sin t )

は、x2 + y2 = a2 となることから、半径 a の円を表します。

0 から t までの曲線の長さ s(t) は

s(t) = ∫{0→t} |p'(u)| du

で求めることができるので、

s'(t) = ds / dt = |p'(t)| = [ ( dx / dt )2 + ( dy / dt )2 ]1/2

となって、t を s の関数として表すことができます。そこで、

p(s) = ( x(s), y(s) )

p が s で表されているとします。すると、

| p'(s) |=( x'(s)2 + y'(s)2 )1/2
=[ ( dx / ds )2 + ( dy / ds )2 ]1/2
={ [ ( dx / dt )( dt / ds ) ]2 + [ ( dy / dt )( dt / ds ) ]2 }1/2
=[ ( dx / dt )2 + ( dy / dt )2 ]1/2( dt / ds )
=( ds / dt )( dt / ds ) = 1

より、p(s) の速さは一定で 1 であることがわかります。

e1 = p'

とすると、e1(s0) は曲線 p(s) の p(s0) における接線方向の単位ベクトルを表します。

ここでもう一つ、単位ベクトル e2 を定義します。e2e1 に直交する単位ベクトルです。そのようなベクトルは 2 つありますが、e2 は、e1 に対して左の方へ向いている側を選ぶものとします。

このとき

( e1, e1 ) = 1

が成り立つので、両辺を微分して

( e'1, e1 ) + ( e1, e'1 ) = 0

より ( e'1, e1 ) = 0 すなわち、e'1e1 に直交します。したがって、ある変数 κ(s) を使って

e'1(s) = κ(s)e2(s)

と表すことができます。この κ(s) を、曲線 p(s) の「曲率 ( Curvature )」といいます。


例として半径 a の円を考えてみましょう。

x(t) = a cos t, y(t) = a sin t

より

s=∫{0→t} [ x'(t)2 + y'(t)2 ]1/2 dt
=∫{0→t} [ ( -a sin t )2 + ( a cos t )2 ]1/2 dt
=∫{0→t} a dt = at

なので、s で表せば

x(s) = a cos ( s / a ), y(s) = a sin ( s / a )

s で微分すると

x'(s) = -sin ( s / a ), y(s) = cos ( s / a )

となって

e1(s) = ( -sin ( s / a ), cos ( s / a ) )

もう一度 s で微分すると

e1'(s) = ( -( 1 / a )cos ( s / a ), -( 1 / a )sin ( s / a ) )

となりますが、e2(s) は e1(s) に直交する単位ベクトルなので

e2(s) = ( -cos ( s / a ), -sin ( s / a ) )

より

e1'(s) = ( 1 / a )e2(s)

となって、κ(s) = 1 / a という結果が得られます。


対象を平面から空間にしても、同様の考え方で曲率を定義することができます。曲線 p(s) によって表される運動の速さが一定で 1 になるようにパラメータ s がとられているとします。すなわち、

e1(s) = p'(s) = ( x'(s), y'(s), z'(s) )

に対して

|e1(s)|2 = x'(s)2 + y'(s)2 + z'(s)2 = 1

であるとします。このとき、加速度ベクトル e1'(s) は

2( e1'(s), e1(s) ) = 0

となるので、e1'(s) は e1(s) に直交していることになります。よって、単位ベクトル e2(s) を

e1'(s) = κ(s)e2(s)

によって定めると、e2(s) は e1(s) に直交する単位ベクトルとなり、| p''(s) | が曲線 p(s) の曲率であるということになります。


曲面上の点 p は、パラメーター u, v を使って

p( u, v ) = ( x( u, v ), y( u, v ), z( u, v ) )

と表すことができます。v = b と固定すると、u が変化することで空間に曲線が描かれます。その曲線が、v を少しずつ変化させることで何本も描かれ、曲面ができると考えれば理解しやすいと思います。

v = b と固定した時に描かれる曲線 p( u, b ) の、点 ( a, b ) における速度ベクトルは

pu( a, b ) ( 但し pu = ∂p / ∂u )

で与えられ、同様に u = a と固定した時に描かれる曲線 p( a, v ) の、点 ( a, b ) における速度ベクトルは

pv( a, b ) ( 但し pv = ∂p / ∂v )

で与えられます。下図は、曲面上の速度ベクトル pu( a, b ), pv( a, b ) の様子を表したものです。

図 n1-1. 曲面上の速度ベクトル
曲面上の速度ベクトル

曲面上の点 ( a, b ) の近傍はほぼ平面と見ることができます。u, v がパラメータ t の関数であるとして、t が極少量 Δt だけ変化したとき、u(t) が u(t) + Δu に、v(t) が v(t) + Δv に変化したとすると、

p( t + Δt ) - p(t) ≅ pu( u(t), v(t) )Δu + pv( u(t), v(t) )Δv

と表すことができます。

図 n1-2. 極少量の変化
極少量の変化

p(t) が p( t + Δt ) に動いたときの距離の 2 乗は、上式の右辺の自分自身との内積をとればいいので

( puΔu + pvΔv, puΔu + pvΔv )=( pu, pu )ΔuΔu + 2( pu, pv )ΔuΔv + ( pv, pv )ΔvΔv
=EΔuΔu + 2FΔuΔv + GΔvΔv

となります。但し、

E = ( pu, pu )
F = ( pu, pv )
G = ( pv, pv )

とします。この式の Δt→0 での極限を

I = Edudu + 2Fdudv + Gdvdv

として、曲面 p( u, v ) の「第一基本形式 ( First Fundamental Form )」と呼びます。

Δp = p( t + Δt ) - p(t)

として Δt→0 での極限を dp と考えれば

I = ( dp, dp )

と表すことができます。

次に、pupv の外積 pu x pv に平行な単位ベクトルを e とします。すなわち

e = pu x pv / | pu x pv |

です。epupv の外積と同じ方向を向いているので、

( pu, e ) = 0
( pv, e ) = 0

が成り立ちます。この式をさらに u, v について偏微分すると

( puu, e ) + ( pu, eu ) = 0
( puv, e ) + ( pu, ev ) = 0
( pvu, e ) + ( pv, eu ) = 0
( pvv, e ) + ( pv, ev ) = 0

となります。puv = pvu なので、

L = ( puu, e ) = -( pu, eu )
M = ( puv, e ) = -( pu, ev )
M = ( pvu, e ) = -( pv, eu )
N = ( pvv, e ) = -( pv, ev )

とすると、

II-( dp, de )
=-( pudu + pvdv, eudu + evdv )
=-( pu, eu )dudu - ( pu, ev )dudv - ( pv, eu )dvdu - ( pv, ev )dvdv
=Ldudu + 2Mdudv + Ndvdv

という結果が得られます。これを曲面 p( u, v ) の「第二基本形式 ( Second Fundamental Form )」といいます。

u, v がパラメータ s の関数で表せるとします。p(s) = p( u(s), v(s) ) として、|dp / ds| = 1 が成り立っているとします。

p''(s) の長さは曲線 p(s) の曲率となるのでした。p'(s) は曲面 p( u, v ) に接していますが、p''(s) は一般には曲面に接してはいません。そこで、

p''(s) = kg + kn

として、曲面に接するベクトル kg と、垂直なベクトル kn に分解します。kg は、曲線 p(s) の「測地的曲率ベクトル ( Geodesic Curvature Vector )」、kn は「法曲率ベクトル ( Normal Curvature Vector )」と呼ばれます。

kn は曲面に垂直なので、e を使って

kn = κne

と表すことができます。κn を曲面上の曲線の「法曲率 ( Normal Curvature )」といいます。

κn=κn( e, e )
=( kn, e )
=( p'' - kg, e )
=( p'', e )

となりますが、

( p', e ) = 0

より、両辺を微分すれば

( p'', e ) + ( p', e' ) = 0

なので

κn=-( p', e' )
=-( pu( du / ds ) + pv( dv / ds ), eu( du / ds ) + ev( dv / ds ) )
=L( du / ds )( du / ds ) + 2M( du / ds )( dv / ds ) + N( dv / ds )( dv / ds )

で求めることができます。|dp| = ds であり、|dp|2 = Edudu + 2Fdudv + Gdvdv だったので、上式は

κn = ( Ldudu + 2Mdudv + Ndvdv ) / ( Edudu + 2Fdudv + Gdvdv )

と変形することができます。

e を含む平面は、e を中心に一周回せば様々に取ることができます。曲面 p と平面との交線である曲線も様々で、法曲率 κn の値も変化するので、その中で必ず最大値と最小値が存在します。

ξ = du
η = dv

として、

2 + 2Mξη + Nη2 - κn( Eξ2 + 2Fξη + Gη2 ) = 0

より ξ, η で偏微分して ∂κn / ∂ξ = ∂κn / ∂η = 0 とすると

2Lξ + 2Mη - κn( 2Eξ + 2Fη ) = 0
2Mξ + 2Nη - κn( 2Fξ + 2Gη ) = 0

となるので、ξ と η について整理して

( L - κnE )ξ + ( M - κnF )η = 0
( M - κnF )ξ + ( N - κnG )η = 0

という連立方程式が得られます。( ξ, η ) ≠ ( 0, 0 ) である解を得るためには、行列式

|L - κnE,M - κnF|
det|M - κnF,N - κnG|

がゼロである必要があります。よって、

( L - κnE )( N - κnG ) - ( M - κnF )( M - κnF ) = 0

を κn について整理して

( EG - F2n2 - ( EN + GL - 2FM )κn + ( LN - M2 ) = 0

となります。解と係数の関係から、解を κ1, κ2 とすると

K ≡ κ1κ2 = ( LN - M2 ) / ( EG - F2 )
H ≡ ( κ1 + κ2 ) / 2 = ( EN + GL - 2FM ) / 2( EG - F2 )

という結果が得られます。K は「ガウス曲率 ( Gaussian Curvature )」、H は「平均曲率 ( Mean Curvature )」といいます。また κ1, κ2 は κn の極値なので、最大値と最小値を表し「主曲率 ( Principal Curvatures )」と呼ばれます。


曲面が z = f( x, y ) の形に表されるとします。このとき、曲面の点 p

p = ( x, y, f( x, y ) )

と表せるので、x, y 方向の偏微分は

px = ( 1, 0, ∂f / ∂x )
py = ( 0, 1, ∂f / ∂y )

であり、∂f / ∂x = fx、∂f / ∂y = fy とすれば

E = ( ( 1, 0, fx ), ( 1, 0, fx ) ) = 1 + fx2
F = ( ( 1, 0, fx ), ( 0, 1, fy ) ) = fxfy
G = ( ( 0, 1, fy ), ( 0, 1, fy ) ) = 1 + fy2

と求められます。また、

px x py=( 0・fy - fx・1, fx・0 - 1・fy, 1・1 - 0・0 )
=( -fx, -fy, 1 )

なので、∂2f / ∂x2 = fxx、∂2f / ∂x∂y = fxy、∂2f / ∂y2 = fyy とすれば

ex = ( -fxx, -fxy, 0 ) / ( fx2 + fy2 + 1 )1/2
ey = ( -fxy, -fyy, 0 ) / ( fx2 + fy2 + 1 )1/2

より

L = -( ( 1, 0, fx ), ( -fxx, -fxy, 0 ) / ( fx2 + fy2 + 1 )1/2 ) = fxx / ( fx2 + fy2 + 1 )1/2
M = -( ( 1, 0, fx ), ( -fxy, -fyy, 0 ) / ( fx2 + fy2 + 1 )1/2 ) = fxy / ( fx2 + fy2 + 1 )1/2
N = -( ( 0, 1, fy ), ( -fxy, -fyy, 0 ) / ( fx2 + fy2 + 1 )1/2 ) = fyy / ( fx2 + fy2 + 1 )1/2

となります。よって、ガウス曲率 K と平均曲率 H は

K={ [ fxx / ( fx2 + fy2 + 1 )1/2 ][ fyy / ( fx2 + fy2 + 1 )1/2 ] - [ fxy / ( fx2 + fy2 + 1 )1/2 ]2 } / [ ( 1 + fx2 )( 1 + fy2 ) - ( fxfy )2 ]
=( fxxfyy - fxy2 ) / ( fx2 + fy2 + 1 )2
H={ ( 1 + fx2 )[ fyy / ( fx2 + fy2 + 1 )1/2 ] + ( 1 + fy2 )[ fxx / ( fx2 + fy2 + 1 )1/2 ] - 2fxfy[ fxy / ( fx2 + fy2 + 1 )1/2 ] } /
 2[ ( 1 + fx2 )( 1 + fy2 ) - ( fxfy )2 ]
=[ ( 1 + fx2 )fyy + ( 1 + fy2 )fxx - 2fxfyfxy ] / 2( fx2 + fy2 + 1 )3/2

で計算することができます。


補足 2) 楕円体の極大値

楕円体を a( x - x0 )2 + b( y - y0 )2 + I0 とします。( 0, 0 ), ( 1, 0 ), ( -1, 0 ), ( 0, 1 ), ( 0, -1 ) 上の値がそれぞれ既知の値 IC, IR, IL, IU, ID だとしたとき、

IC = ax02 + by02 + I0
IR = a( 1 - x0 )2 + by02 + I0
IL = a( -1 - x0 )2 + by02 + I0
IU = ax02 + b( 1 - y0 )2 + I0
ID = ax02 + b( -1 - y0 )2 + I0

が成り立ちます。よって、

IR - IC=a( 1 - x0 )2 - ax02
=a - 2ax0
IC - IL=ax02 - a( -1 - x0 )2
=-a - 2ax0

より

a = [ ( IR - IC ) - ( IC - IL ) ] / 2

また、

IU - IC=b( 1 - y0 )2 - by02
=b - 2by0
IC - ID=by02 - b( -1 - y0 )2
=-b - 2by0

より

b = [ ( IU - IC ) - ( IC - ID ) ] / 2

となります。また、

x0 = [ 1 - ( IR - IC ) / a ] / 2
y0 = [ 1 - ( IU - IC ) / b ] / 2

なので、

I0 = IC - a[ 1 - ( IR - IC ) / a ]2 / 4 - b[ 1 - ( IU - IC ) / b ]2 / 4

で極大値の I0 を求めることができます。


補足 3) 行列式・トレースと固有値の関係

任意の正方行列 A に対して、行列式 det( A ) は固有値の積に等しく、トレース tr( A ) は固有値の和に等しくなります。すなわち、A の固有値を λi としたとき、

det( A ) = Πi( λi )
tr( A ) = Σi( λi )

となります。

行列 A の「固有方程式 ( Characteristic Polynomial )」は次のような式になるのでした ( n3-1 )。

det( tE - A ) = 0

但し、E は単位行列を表します。行列 A のサイズが n であるとし、第 j 行 i 列の成分を aji としたとき、固有方程式は

|t - a11,-a12,...-a1n|
|-a21,t - a22,...-a2n|
|::...:|
det|-an1,-an2,...t - ann| = 0

と表すことができます。行列式は、各列から相異なる行の要素を選んで ( または各行から相異なる列の要素を選んで ) 積を計算し、奇置換ならば符号を反転した上で、その和を求めた結果になります ( n3-2 )。したがって、上式の tn の項は対角成分の積からしか得られず、その係数は明らかに 1 です。また、tn-1 の項も、対角成分の積以外は多くても tn-2 の項までしか得られないので、やはり対角成分の積にしか存在せず、その係数は -tr( A ) = -Σi{1→n}( aii ) に等しくなります。なぜならば、n = 1 のときは明らかであり、n - 1 まで tn-2 の項の係数が -Σi{1→n-1}( aii ) ならば、

[ tn-1 - Σi{1→n-1}( aii )tn-2 + o(tn-3) ]( t - ann )
=tn - anntn-1 - Σi{1→n-1}( aii )tn-1 + o(tn-2)
=tn - Σi{1→n}( aii )tn-1 + o(tn-2)

と計算できます。但し、o(tn) は n 次以下の項を表します。よって、帰納法により成り立つことがわかります。

n 個の固有値を λi ( i = 1, 2, ... n ) としたとき、固有方程式は

Πi{1→n}( t - λi ) = 0

となります。このとき、左辺の tn-1 項の係数は、先に証明したように -Σi{1→n}( λi ) に等しくなります。よって、係数を見比べることで

tr( A ) = Σi{1→n}( λi )

という結果が得られます。


次に、固有方程式の定数項が det( A ) に等しいことを証明します。

det( tE - A ) の定数項は、t = 0 のときの値を求めればいいので、det( -A ) になります。行列式の定義から明らかなように、det( -A ) = ( -1 )ndet( A ) です。

一方、固有値 λi を使って固有方程式を表すと

Πi{1→n}( t - λi ) = 0

だったので、この左辺の定数項も t = 0 のときの値を考えればよく、( -1 )nΠi{1→n}( λi ) となります。

二つの式を見比べれば

det( A ) = Πi{1→n}( λi )

ということになります。

*n3-1)数値演算法 (8) 連立方程式を解く -2-」の「5) 固有値と固有ベクトル」を参照

*n3-2)確率・統計 (5) 正規分布」の「補足4) 行列の積の行列式」を参照



<参考文献>
  1. 「画像認識」 原田達也著 ( 講談社 )
  2. 「曲線と曲面の微分幾何」 小林昭七著 ( 裳華房 )
  3. 近藤弘一様のページ」 - 「線形代数学 II」 - 「5.9 固有多項式とトレース
  4. Wikipedia

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