画像処理 (3) 顕著性マップ

人は、目から得られるたくさんの情報を常に処理し、それに対して判断を行っています。目に映るモノ全てに対してそれが何かをいちいち判断していたのでは、特に瞬時の判断が必要な場合は生死に関わることになる可能性もあります。車の運転中、前の車が急停車したという情報を処理するときは、目の前に近づいてきた車やそのブレーキランプを見て判断すれば充分なのに、周囲の景色や道路標識なども常に処理されていては判断が間に合わなくなります。そのため、本当に注目すべきモノだけに無意識に注意を引くような仕組みが人には備わっています。これは、人に限らず他の動物にもあり、天敵などの危険から逃げたり、逆に獲物を追う時に逃げる対象の動きを瞬時にとらえるなど、常に視覚情報から注意を引く部分を抽出して、それ以外の情報はカットしてしまう機能が必要になる機会は人よりも多いことが想像できます。

今回は、注意を引く対象を見つけるための仕組みとして考案された「顕著性マップ ( Saliency Map )」を紹介し、これをシーム・カービング処理に応用してみたいと思います。

1) ガウシアン・フィルタ

前章で紹介したガボール・フィルタは、「正弦波」と「ガウス関数」の積として定義されていました。この中から「正弦波」を取り除いて、「ガウス関数」だけを使ったフィルタが「ガウシアン・フィルタ ( Gaussian Filter )」です。

以下の式が成り立つので、積分の結果を 1 にするために係数 1 / sqrt( 2π )σ, 1 / 2πσ² が付けられています (*1-1)。

ガウス関数のフーリエ変換はやはりガウス関数になります。

g(x) は平均 0、分散 σ² の「正規分布 ( Normal Distribution )」を表し、0 を中心に釣鐘状の分布を持っています。また、分散が大きいほどデータのばらつきが大きいことになり、分布は広範囲に広がることになります。そのフーリエ変換である G(ω) は平均 0、分散 1/σ² の正規分布となっていることが上式からわかります。このフィルタを使って信号を畳み込み積分すると、低周波数成分だけを増幅させて高周波領域をカットする「低域フィルタ」として作用することになり、分散を大きくするほどそのフーリエ変換の幅は狭くなるため、低周波領域が極端に強調されるようになります。これは、畳み込み積分によって周囲の信号をいっしょに積算する度合いが増えることにより値が平滑化されることからも推測ができます。このことから、ガウシアン・フィルタを画像に作用されることによって "ぼかし" を掛けたような効果が得られます。これは「ガウシアンぼかし ( Gaussian Blur )」と呼ばれる画像処理です。

ガウシアンぼかしのサンプル・プログラムを以下に示します。

/*
  Check_RGBPrimary : p が RGB 各成分値の範囲内にあるかを判定し、そうでなければ範囲内に収めて返す
*/
template< class T > unsigned char Check_RGBPrimary( T p )
{
  if ( p < 0 ) p = 0;
  if ( p > std::numeric_limits< unsigned char >::max() )
    p = std::numeric_limits< unsigned char >::max();

  return( static_cast< unsigned char >( p ) );
}

/******************************************
 GaussianFilter : Gaussian Filter Kernel
******************************************/
class GaussianFilter
{
  std::vector< double > gaussKernel_; // カーネル
  unsigned int wSize_;                // カーネルサイズ
  int s3_;                            // 3σ

 public:

  // ガウス関数の標準偏差 sigma を指定して構築
  GaussianFilter( double sigma = 1.0 );

  // 座標 c のガウス関数の値を返す
  double operator[]( Coord< int > c ) const
  {
    int i = ( c.y + s3_ ) * wSize_ + c.x + s3_;

    return( ( i < 0 || (unsigned int)i >= gaussKernel_.size() ) ? 0 : gaussKernel_[i] );
  }

  // ガウス関数のパラメータ 3σ を返す
  int s3() const { return( s3_ ); }
};

/*
  GaussianFilter : コンストラクタ

  sigma : 標準偏差

  カーネルサイズは ± 3σ までとする。
*/
GaussianFilter::GaussianFilter( double sigma )
{
  if ( sigma <= 0 ) sigma = 1; // σ が負数なら 1 にする

  s3_ = round( 3 * sigma );    // 3σ
  if ( s3_ == 0 ) s3_ = 1;
  wSize_ = static_cast< unsigned int >( s3_ ) * 2 + 1; // カーネルサイズ

  gaussKernel_.resize( wSize_ * wSize_ );

  // カーネルの計算
  for ( int y = -s3_ ; y <= s3_ ; ++y ) {
    for ( int x = -s3_ ; x <= s3_ ; ++x ) {
      gaussKernel_[static_cast< size_t >( y + s3_ ) * wSize_ + static_cast< size_t >( x + s3_ )] =
        exp( - ( pow( x, 2 ) + pow( y, 2 ) ) / ( 2 * pow( sigma, 2 ) ) ) / ( 2 * M_PI * pow( sigma, 2 ) );
    }
  }
}

/*
  GaussianConvolution : ガウス関数による畳み込み積分処理

  draw : 描画対象オブジェクト
  gc : 画像の対象位置
  filter : ガウシアン・フィルタ
  rgb : 積分結果
*/
void GaussianConvolution( const DrawingArea_IF& draw, Coord< int > gc, const GaussianFilter& filter, RGB* rgb )
{
  RGB p;
  double r = 0, g = 0, b = 0;
  Coord< int > c;
  int s3 = filter.s3();

  for ( c.y = -s3 ; c.y <= s3 ; ++( c.y ) ) {
    for ( c.x = -s3 ; c.x <= s3 ; ++( c.x ) ) {
      Coord< int > gc2( ( gc.x - c.x < 0 ) ? 0 : ( ( gc.x - c.x >= draw.size().x ) ? draw.size().x - 1 : gc.x - c.x ),
                        ( gc.y - c.y < 0 ) ? 0 : ( ( gc.y - c.y >= draw.size().y ) ? draw.size().y - 1 : gc.y - c.y ) );
      draw.point( gc2, &p );
      r += p.r * filter[c];
      g += p.g * filter[c];
      b += p.b * filter[c];
    }
  }

  *rgb = RGB( Check_RGBPrimary( r ), Check_RGBPrimary( g ), Check_RGBPrimary( b ) );
}

/*
  GaussianBlur : 描画対象 draw にガウシアンぼかしを行う

  sigma : ガウシアン・フィルタの標準偏差
*/
void GaussianBlur( DrawingArea_IF* draw, double sigma )
{
  GaussianFilter filter( sigma );
  vector< RGB > result( ( draw->size() ).x * ( draw->size() ).y );
  Coord< int > c;

  // 畳み込み積分処理
  for ( c.y = 0 ; c.y < ( draw->size() ).y ; ++( c.y ) ) {
    for ( c.x = 0 ; c.x < ( draw->size() ).x ; ++( c.x ) ) {
      GaussianConvolution( *draw, c, filter, &result[c.y * ( draw->size() ).x + c.x] );
    }
  }

  // 積分結果を画像に反映
  for ( c.y = 0 ; c.y < ( draw->size() ).y ; ++( c.y ) )
    for ( c.x = 0 ; c.x < ( draw->size() ).x ; ++( c.x ) )
      draw->pset( c, result[c.y * ( draw->size() ).x + c.x] );
}

特に難しい部分はなく、各ピクセル位置におけるガウス関数の値をコンストラクタで作成して、画像内の各画素について、その周囲のピクセルといっしょに畳み込み積分を行っているだけです。GaussianBlur がガウシアンぼかしを行うためのメイン・ルーチンで、その中で畳み込み積分を行うため GaussianConvolution を呼び出しています。また、GaussianFilter はガウシアン・フィルタを作成するためのクラスになります。積分を行う範囲は、対象画素から ±3σ の範囲にしてあります。離散値を使った和の形で畳み込み積分を行っているため、3σがピクセル幅の 1 より小さくなった場合は正しく計算することができません。ここでは、3σが 1 より小さい値になった場合は強制的に 1 に補正しています。

ガウシアンぼかしを使ったテスト結果は次のようになります。テストに利用した画像は The USC-SIPI Image Database にある Peppers という画像です(処理に利用した実際の画像は 512 x 512 の大きさがあります)。

下図は、処理結果の中央部分を 256 x 256 の大きさ分だけ切り出した画像を示しています。

σを大きくするほど低周波領域の成分が強調されるため、結果として画像がぼやけた感じになる様子を上の結果から見ることができます。

ガウシアンぼかしを行いながら、画素のサンプリングを縦横ともに一つおきに行うと、結果として得られる画像の大きさは縦横それぞれ半分となります。サンプリングされた画素には、その周囲の画素の色成分がガウス関数の分布に従って加味されるので、ある程度周囲の画素の影響を受けることになります。この性質を利用して、通常のサンプル補間によって発生するエイリアシングを抑えることができます (*1-2)。この縮小処理を繰り返して作られた画像を重ねていくと、ちょうどピラミッドのようになるので、これを「ガウシアン・ピラミッド(Gaussian Pyramid)」といいます。

ガウシアン・ピラミッドを作成するには処理する位置を一つおきにするのを繰り返すだけで実現できます。ガウシアンぼかしの機能を利用してガウシアン・ピラミッドを作成するサンプル・プログラムを以下に示します。

/*
  CreateGaussianPyramid : ガウシアン・ピラミッドの作成

  一つおきに畳み込み積分処理を行う。

  src : 元画像オブジェクト
  dest : 処理結果を格納する描画オブジェクト
  sigma : ガウシアン・フィルタの標準偏差
*/
void CreateGaussianPyramid( const DrawingArea_IF& src, DrawingArea_IF* dest, double sigma )
{
  GaussianFilter filter( sigma );
  Coord< int > d( ( src.size().x + 1 ) / 2, ( src.size().y + 1 ) / 2 );
  vector< RGB > result( d.x * d.y );
  Coord< int > c;

  dest->resize( d );

  // 畳み込み積分処理
  for ( c.y = 0 ; c.y < d.y ; ++( c.y ) ) {
    for ( c.x = 0 ; c.x < d.x ; ++( c.x ) ) {
      GaussianConvolution( src, Coord< int >( c.x * 2, c.y * 2 ), filter, &result[c.y * d.x + c.x] );
    }
  }

  // 積分結果を画像に反映
  for ( c.y = 0 ; c.y < d.y ; ++( c.y ) )
    for ( c.x = 0 ; c.x < d.x ; ++( c.x ) )
      dest->pset( c, result[c.y * d.x + c.x] );
}

σ = 1 として、上記プログラムによる処理を繰り返した結果を以下に示します。下図は、原画像の中央部分 256 x 256 の領域に対してガウシアン・ピラミッドの作成を 5 回繰り返した結果です。

下図は、各処理結果を原画像と同じ大きさに拡大したものです。

処理後に原画像と同じ大きさに拡大処理して同位置にある画素を比較すると、処理後の画素はその周囲にある画素の影響を受けています。例えば、一点だけ白で周囲は全て黒であるような場合、周囲にある黒色によって処理後の画素はかなり黒に近い色になるので、処理前後の値の差は大きくなります。逆に、周囲の色とほとんど差のないような場合は、処理前後の値の差は小さくなり、極端な例として全てが同色であるようなときは、差はゼロになります。このように、ガウシアン・ピラミッドによって作成された画像で差を取ってみると、エッジなどのような周囲と差のある部分の値が大きくなる傾向が見られることになります。

2) 網膜神経節細胞の働き

目の網膜には、網膜神経節細胞 ( Retinal Ganglion Cell ; 以下 RGC と略記 ) と呼ばれる神経細胞があり、人の場合は 120 万～ 150 万個程度存在すると言われています。さらに、RGC の中には受容野 ( Receptive Field ) と呼ばれる領域があって、ここで光による刺激を受けて、その情報が脳へと伝達されます。受容野は中央にある円形の部分 ( Center ) とその周辺にある領域 ( Surround ) の二つで構成されており、それぞれの光の受け方によって挙動が変化します。
受容野には「オン中心型 ( On Center Cell )」と「オフ中心型 ( Off Center Cell )」の二種類が存在し、前者は Center のみ、後者は Surround のみに光による刺激があったときに強く反応します。光による刺激がなければ当然反応はしませんが、Center と Surround の両方に光による刺激があってもその反応は抑制され、弱い信号しか脳へ伝達されません。このような仕組みによって、光の強度の差 ( コントラスト ; Contrast ) やエッジの検出が可能になります。

ガウシアン・ピラミッドにおいて、スケールの異なる画像を比較すると、スケールの大きい ( 画素の細かい ) 画像上の画素に対して、スケールの小さい ( 画素の粗い ) 画像上の画素はその周辺の値を表しています。前述の通り、それぞれの画像の差分を取ってみると、周囲に比べて差異のある画素が顕著に現れるようになります。この仕組みは、ちょうど RGC にある受容野の働きによく似ています。

受容野は様々な大きさのものが用意されており、小さな受容野は画素の細かい画像 ( 高周波空間 )、大きな受容野は画素の粗い画像 (低周波空間 ) を使った処理に等しくなります。よって、ガウシアン・ピラミッドによって得られたスケールの異なるペアをいくつか使って差異を求め、それらを重ね合わせることで、脳へ伝わる信号が強くなる箇所 ( = 注意を引く場所 ) を特定できると考えられます。この手法は「Center-Surround Operation」と呼ばれ、「Laurent Itti」「Christof Koch」「 Ernst Niebur」によって 1998 年に考案されました。また、この手法によって得られた、画像上の注意を引く箇所を示したマップは「顕著性マップ ( Saliency Map )」と名付けられました。次の章で、Itti らが考案した顕著性マップの作成方法を文献にしたがって紹介したいと思います。

3) 顕著性マップ

顕著性マップは、画像のガウシアン・ピラミッドを作成するところからスタートします。段階的に縮小された画像は、原画像から 1 / 256 ( = 1 / 2⁸ ) 縮尺まで計 9 枚用意されます。これらの画像に対して連番 σ = 0 ( 原画像 ) ～ 8 ( 1 / 256 縮尺 ) を付加しておきます。次に、各画像それぞれに対して、輝度 ( Intensity )・色相成分 ( Color Channel )・方向成分 ( Orientation ) を計算していきます。以下、この画像を「スケール画像」と呼ぶことにします。

σ 番目のスケール画像の RGB 成分を、赤成分 r(σ)・緑成分 g(σ)・青成分 b(σ) としたとき、σ 番目のスケール画像に対する各ピクセルの輝度 I(σ) は次のように計算されます。

輝度を求めたら、色彩に対する輝度の影響を除外するため RGB 各成分を輝度で正規化します。例えば、二つのピクセルにおいて、片側は ( r, g, b ) = ( 100, 0, 0 )、もう片側は ( r, g, b ) = ( 10, 0, 0 ) であった場合、輝度は前者の方が大きくなりますが、どちらも赤成分しか持っていないので、彩度としては同程度と見なすことになります。但し、輝度の非常に低い領域に対しては色彩を知覚することができないので、輝度の最大値の 1/10 より小さい輝度であるピクセルに対しては RGB 成分をゼロとします。
正規化した RGB 成分を使って、色相成分として赤 R(σ)・緑 G(σ)・青 B(σ)・黄 Y(σ) 各値を以下の式で求めます。

但し、計算結果が負数になった場合はゼロとします。
他の色相成分の値も大きければ、減算により上記計算の結果は小さくなるので、ある色相成分のみ大きな値を持っている場合のみ計算結果が大きくなります。

最後に、方向成分を求めます。ここでは前回紹介したガボール・フィルタを利用して、θ = 0°, 45°, 90°, 135° の四方向に対して畳み込み積分を行います。また、畳み込みに利用する値は輝度 I(σ) になります。ガボール・フィルタを ψ(θ) としたとき、方向成分 O( σ, θ ) は次のように計算されます。

上式において、I( σ, m, n ) は σ 番目のスケール画像における ( m, n ) 座標のピクセルの輝度、ψ( θ, x - m, y - n ) はガボール・フィルタ ψ(θ) の ( x - m, y - n ) 座標の成分を表しています。また、h_m と h_n はそれぞれガボール・フィルタの x, y 方向の窓サイズ ( フィルタ中央からの距離 ) を表しています。ガボール・フィルタは原点からある程度離れるとゼロに収束するので、途中で計算処理を打ち切ってもある程度の精度を保つことができます。

以上の処理により、I(σ), R(σ), G(σ), B(σ), Y(σ) がそれぞれ 9 枚ずつ、O( σ, θ ) は四方向分あるので計 36 枚、全て合わせて 81 枚の画像が作られます。

ガウシアン・ピラミッドを作成するためのサンプル・プログラムを以下に示します。

/******************************************
 GaussianPyramid : ガウシアン・ピラミッド
******************************************/
class GaussianPyramid
{
  // 各スケールの画像用パラメータ
  struct ScaleImage
  {
    static const size_t O_SIZE = 4; // 方向成分の数

    struct ColorComponent
    {
      double r_; // 赤成分
      double g_; // 緑成分
      double b_; // 青成分
      double y_; // 黄成分
    };

    Matrix< double > i_;          // 輝度
    Matrix< ColorComponent > cc_; // 色成分
    Matrix< double > o_[O_SIZE];  // 方向成分(0,45,90,135)

    // 画像データ draw とガボール・フィルタ群 gabor から構築
    ScaleImage( const DrawingArea_IF& draw, const GaborFilterGroup& gabor );
  };

 public:

  /* 公開メンバ定数 */
  static const size_t IMAGE_COUNT = 7; // ScaleImageの数(σ = 2～8)

 private:

  std::vector< ScaleImage > scaleImage_; // 各スケールの画像(σ = 2～8)

  // スケール画像の番号 i と座標値 c からパラメータの有効性を評価
  bool isValid( size_t i, Coord< int >* c ) const;

 public:

  /*
    画像データとフィルタのパラメータを指定して構築

    draw : 対象画像データ
    sigma : ガウシアン・フィルタの標準偏差
    lambda : 方向成分で利用するガボール・フィルタの波長
  */
  GaussianPyramid( const DrawingArea_IF& draw, double sigma = 1.0, double lambda = GaborFilter::DEFAULT_LAMBDA );

  /* 各成分の値を返す */

  // 画像番号 i・位置 c の輝度を返す
  double lumi( size_t i, Coord< int > c ) const
  { return( ( ! isValid( i, &c ) ) ? 0 :
            scaleImage_[i].i_[c.y][c.x] ); }

  // 画像番号 i・位置 c の色彩(赤)を返す
  double red( size_t i, Coord< int > c ) const
  { return( ( ! isValid( i, &c ) ) ? 0 :
            scaleImage_[i].cc_[c.y][c.x].r_ ); }

  // 画像番号 i・位置 c の色彩(緑)を返す
  double green( size_t i, Coord< int > c ) const
  { return( ( ! isValid( i, &c ) ) ? 0 :
            scaleImage_[i].cc_[c.y][c.x].g_ ); }

  // 画像番号 i・位置 c の色彩(青)を返す
  double blue( size_t i, Coord< int > c ) const
  { return( ( ! isValid( i, &c ) ) ? 0 :
            scaleImage_[i].cc_[c.y][c.x].b_ ); }

  // 画像番号 i・位置 c の色彩(黄)を返す
  double yellow( size_t i, Coord< int > c ) const
  { return( ( ! isValid( i, &c ) ) ? 0 :
            scaleImage_[i].cc_[c.y][c.x].y_ ); }

  // 画像番号 i・位置 c・方向 theta の方向成分を返す
  //
  // theta : 方向成分( 0 = 0° ; 1 = 45° ; 2 = 90° ; 3 = 135° )
  double ori( size_t i, size_t theta, Coord< int > c ) const;

  // 画像番号 i の大きさを返す
  Coord< int > imageSize( size_t i ) const;
};

/*
  GaussianPyramid::ScaleImage コンストラクタ

  draw : 描画対象オブジェクト
  gabor : 方向成分計算用のガボール・フィルタ群
*/
GaussianPyramid::ScaleImage::ScaleImage( const DrawingArea_IF& draw, const GaborFilterGroup& gabor )
{
  assert( gabor.size() == O_SIZE );

  Coord< int > size = draw.size(); // 画像のサイズ
  double iMax = std::numeric_limits< double >::min(); // 輝度の最大値

  // 要素数分だけ領域を確保
  i_.assign( size.y, size.x );
  cc_.assign( size.y, size.x );

  RGB rgb; // drawから取得する色コード

  // 輝度と色成分のセット
  Matrix< double >::iterator iIt = i_.begin();
  Matrix< ColorComponent >::iterator cIt = cc_.begin();
  for ( Coord< int > c( 0, 0 ) ; c.x < size.x ; ++( c.x ) ) {
    for ( c.y = 0 ; c.y < size.y ; ++( c.y ) ) {
      draw.point( c, &rgb );

      double red = rgb.r;
      double green = rgb.g;
      double blue = rgb.b;

      *iIt = ( red + green + blue ) / 3; // I = ( r + g + b ) / 3
      iMax = std::max( iMax, *iIt );
      cIt->r_ = std::max( ( red - ( green + blue ) / 2 ) / *iIt, 0.0 ); // R = r - ( g + b ) / 2
      cIt->g_ = std::max( ( green - ( red + blue ) / 2 ) / *iIt, 0.0 ); // G = g - ( r + b ) / 2
      cIt->b_ = std::max( ( blue - ( red + green ) / 2 ) / *iIt, 0.0 ); // B = b - ( r + g ) / 2
      cIt->y_ = std::max( ( ( red + green ) / 2 - std::abs( red - green ) / 2 - blue ) / *iIt, 0.0 ); // Y = ( r + g ) / 2 - | r - g | / 2 - b
      ++iIt; ++cIt;
    }
  }

  // 色成分を正規化
  for ( iIt = i_.begin(), cIt = cc_.begin() ; iIt != i_.end() ; ++iIt, ++cIt )
    if ( *iIt < iMax / 10.0 )
      cIt->r_ = cIt->g_ = cIt->b_ = cIt->y_ = 0;

  // 方向成分のセット
  size_t k = 0; // 各要素をセットする位置
  for ( GaborFilterGroup::const_iterator g = gabor.begin() ; g != gabor.end() ; ++g ) {
    o_[k].assign( size.y, size.x );
    Matrix< double >::iterator oIt = o_[k].begin();
    for ( Coord< int > c( 0, 0 ) ; c.x < size.x ; ++( c.x ) )
      for ( c.y = 0 ; c.y < size.y ; ++( c.y ) )
        *oIt++ = GaborFunction::Convolution( i_, c, *g );
    ++k;
  }
}

/*
  GaussianPyramid コンストラクタ

  draw : 描画対象オブジェクト
  sigma : ガウシアン・フィルタにセットする標準偏差 σ
  lambda : ガボール・フィルタにセットする波長 λ
*/
GaussianPyramid::GaussianPyramid( const DrawingArea_IF& draw, double sigma, double lambda )
{
  using GaussianFunction::CreateGaussianPyramid;

  GaborFilterGroup gabor; // ガボール・フィルタ
  // 0,45,90,135° でフィルタを作成
  GaborFunction::CreateSymmetric( &gabor, lambda, ScaleImage::O_SIZE );

  // ガウシアン・ピラミッド
  GPattern p0( Coord< int >( 0, 0 ) );
  GPattern p1( Coord< int >( 0, 0 ) );

  // ガウシアン・ピラミッドを作成
  // 画像は先に 1/4 の大きさに処理してからパラメータをセットする
  CreateGaussianPyramid( draw, &p0, sigma ); // p0 → 1/2 scale
  CreateGaussianPyramid( p0, &p1, sigma );   // p1 → 1/4 scale
  for ( size_t i = 0 ; i < IMAGE_COUNT ; ++i ) {
    if ( ( i & 1 ) == 0 ) {
      scaleImage_.push_back( ScaleImage( p1, gabor ) );
      CreateGaussianPyramid( p1, &p0, sigma );
    } else {
      scaleImage_.push_back( ScaleImage( p0, gabor ) );
      CreateGaussianPyramid( p0, &p1, sigma );
    }
  }
}

/*
  GaussianPyramid::isValid : スケール画像の番号 i と座標値 c からパラメータの有効性を評価

  座標値 c が画像範囲外にある場合は補正を行う

  戻り値 番号 i が有効なら true を返す
*/
bool GaussianPyramid::isValid( size_t i, Coord< int >* c ) const
{
  if ( i >= IMAGE_COUNT ) return( false );
  if ( scaleImage_[i].i_.cols() == 0 || scaleImage_[i].i_.rows() == 0 )
    return( false );

  if ( c->x < 0 ) c->x = 0;
  if ( c->y < 0 ) c->y = 0;
  if ( static_cast< size_t >( c->x ) >= scaleImage_[i].i_.cols() )
    c->x = scaleImage_[i].i_.cols() - 1;
  if ( static_cast< size_t >( c->y ) >= scaleImage_[i].i_.rows() )
    c->y = scaleImage_[i].i_.rows() - 1;

  return( true );
}

/*
  GaussianPyramid::ori : 方向成分の取得

  i : スケール画像の番号
  theta : 方向番号( 0° = 0 , 45° = 1 , 90° = 2 , 135° = 3 )
  c : 座標値

  戻り値 : 方向成分(指定パラメータが無効ならゼロを返す)
*/
double GaussianPyramid::ori( size_t i, size_t theta, Coord< int > c ) const
{
  if ( ! isValid( i, &c ) ) return( 0 );
  if ( theta >= ScaleImage::O_SIZE ) return( 0 );

  return( ( scaleImage_[i].o_ )[theta][c.y][c.x] / 8 );
}

/*
  GaussianPyramid::imageSize : スケール画像の番号 i の大きさを返す

  戻り値 : マップの大きさ ( 指定パラメータが無効なら Coord< int >( 0, 0 ) を返す )
*/
Coord< int > GaussianPyramid::imageSize( size_t i ) const
{
  if ( i >= IMAGE_COUNT ) return( Coord< int >( 0, 0 ) );

  return( Coord< int >( scaleImage_[i].i_.cols(), scaleImage_[i].i_.rows() ) );
}

各スケール画像は GaussianPyramid オブジェクトの中で作成されます。引数として画像オブジェクトを渡すと、ガウシアン・フィルタでスケールを縮小させながら、輝度・色相・方向の各成分を順番に処理していきます。各スケール画像毎に ScaleImage オブジェクトが作られ、各成分の計算はこの ScaleImage のインスタンス化の時に行われます。輝度による色成分の正規化は、前述の順番と異なり、先に R(σ), G(σ), B(σ), Y(σ) を計算してから行っています。ここでは単純に輝度 I(σ) で徐算するだけとしているので、順番を変えても結果は変わりません。最後に方向成分を計算するためガボール・フィルタによる畳み込み積分処理を行います。ガボール・フィルタは各スケール画像共通で利用できるので、GaussianPyramid コンストラクタ側で作成したものを ScaleImage コンストラクタの引数として渡しています。
計算した各成分の値を返すため、lumi, red, green, blue, yellow, ori がメンバ関数として用意され、それぞれ輝度、色相成分 ( 赤・緑・青・黄 )、方向成分を返します。引数としてスケール画像の番号と対象ピクセルの座標値、また方向成分では角度を表す連番を渡す仕様になっています。

この後、各スケール画像毎の差分を計算していくわけですが、後述するように差分計算で利用するのは σ ≥ 2 のスケール画像であり、原画像と σ = 1 のスケール画像は利用しないため、GaussianPyramid オブジェクトにはそれらのデータは含まれていません。よって、スケール画像の数を表す imageCount は 9 ではなく 7 となっています。また、スケール画像の番号を指定するとき、0 が σ = 2 の画像のインデックスを表すことになります。

続いて、各スケール画像間の差分を求めます。各スケール画像の大きさはバラバラなので、サイズの小さい側を大きい側に合わせて伸張した上でピクセル毎に差分を計算します。以下、この操作を演算子 "(-)" で表すこととします。サイズの大きい側は中心 ( Center ) の画素、サイズの小さい側は周辺 ( Surround ) の画素を含む平均をそれぞれ表すことになるので、それぞれの番号を c, s ( c < s ) とします。まず、輝度に対する差異 I( c, s ) は

と単純な式で計算できます。
色相成分は少し処理が複雑です。受容野の中央部分において、ある色相によってニューロンが活性化し、逆に別の色相では抑制される仕組みがあります。また、その周辺部においては逆に、中央で抑制される色相で活性化し、中央で活性化する色相では抑制されます。この仕組みを文献では「Color Double-Opponent System」と呼んでいます。人の第一視覚野では、そのような色相の組み合わせとして赤と緑、青と黄があります。従って、赤と緑 ( 青と黄 ) の差異について中央部と周辺部の差を計算することで、色相成分での顕著性を得ます。式としては次のようになります。

RG( c, s ) が赤と緑の、BY( c, s ) が青と黄の色相の差を同時に表すことになります。実際に値を代入して計算した例が下表に示してあります。

表 3-1. **Color Double-Opponent System** による色相成分の顕著性計算例
Center	Surround	RG(c,s)
R(c)	G(c)	R(s)	G(s)
0	0	0	0	0
100	0	100
0	100	100
100	100	0
100	0	0	0	100
100	0	200
0	100	0
100	100	100
0	100	0	0	100
100	0	0
0	100	200
100	100	100
100	100	0	0	0
100	0	100
0	100	100
100	100	0

中央部と周辺部が赤と緑で反転しているような場合、受容野はどちらも活性化するか、もしくはどちらも抑制された形になり、その場合は RGC の働きのところで説明したように信号は弱くなります。上記の式は、その様子をうまく表現しています。

最後の方向成分は、各方向毎に Center - Surround 間で差分を計算します。

文献では、c ∈ { 2, 3, 4 } 、s ∈ { c + δ | δ ∈ { 3, 4 } } としています。従って、差分を計算した結果は計 6 パターン ( X(2)-X(5), X(2)-X(6), X(3)-X(6), X(3)-X(7), X(4)-X(7), X(4)-X(8) ) になります。また、輝度が一種類、色相成分が二種類、方向成分は四種類あるので、合計すると 42 枚の画像が得られます。これらのデータは文献では「特性マップ ( Feature Map )」と呼んでいます。

最後に特性マップを重ね合わせて顕著性マップが完成するのですが、ここで問題になるのが各成分における値の範囲 ( Dynamic Range ) の違いです。それぞれの成分は異なる計算方法を取っているので、例えば輝度と 0 度の方向成分に対して値の強弱を判断することができません。また、42 枚もの画像を重ね合わせることになるので、それぞれの画像が持っているノイズが積算されることで、重要な顕著性が隠れてしまう恐れもあります。

そこで、各特性マップに対して正規化処理を行います。文献では次のような手法を用いて正規化処理を行っています。

マップ上のある一点だけ他よりも大きな値を持っていれば、最大値は M、極大値はゼロになるので、正規化によってその一点だけ M³ となって、その他は全てゼロになります。逆に、最大値に近い値を持った極大値がたくさんあれば、極大値の平均値 m は M に非常に近くなり、最大値を含む全てのピクセルの値は非常に小さな値となります。

特性マップは個々に正規化処理を行います。重ね合わせ処理は輝度、色相成分、方向成分それぞれに対して実施して、さらにその結果を正規化処理した上で全てのマップを重ね合わせます。方向成分の重ね合わせ処理は他と異なり、各方向 (0°, 45°, 90°, 135°) 単位で重ね合わせた結果を正規化処理した上でそれらを重ね合わせる処理を行っています。重ね合わせ処理を行う各マップの大きさは異なるため、差分の時と同様に小さい側を拡大処理してピクセル毎に加算処理を行います。これを和の計算と同じような書式で [+]_k{k₁→k₂}<・> と表し、また正規化処理を N(・) で表すと、輝度 I、色相成分 C、方向成分 O の重ね合わせ処理は下式で表されます。

最後にもう一度各成分に対して正規化処理を行い、重ね合わせ処理を行って顕著性マップ S が完成します。

K_I, K_C, K_O は各成分の重み付けを表しており、文献では全て等しい値 ( 1 / 3 ) にしています。

顕著性マップを求めるためのサンプル・プログラムを以下に示します。

/****************************************************
 CenterSurroundDiff : Center-Surroundスケール間差分
****************************************************/
class CenterSurroundDiff
{
  typedef Matrix< double >::size_type size_type; // サイズの型

  static const size_t O_SIZE = 4; // 方向成分の数

  Coord< int > sz_; // Gauusian Pyramidの一番目の画像(σ = 2)の大きさ

  // スケール間差分の和
  Matrix< double > iSum_;             // 輝度成分
  Matrix< double > cSum_;             // 色相成分
  Matrix< double > oSum_[O_SIZE + 1]; // 方向成分( 0 = 0° ; 1 = 45° ; 2 = 90° ; 3 = 135°)
  Matrix< double > oSumComb_;         // 方向成分(合成)

  // 顕著性
  Matrix< double > saliency_;

  // Center-Surround差分を求める
  void createCSMap( const GaussianPyramid& pyramid );

  // 座標が有効であるか判定する
  bool isValid( Coord< int > c ) const
  { return( c.x >= 0 && c.x < sz_.x &&
            c.y >= 0 && c.y < sz_.y ); }

 public:

  // コンストラクタ
  CenterSurroundDiff( const GaussianPyramid& pyramid, double ki = 1.0/3.0, double kc = 1.0/3.0, double ko = 1.0/3.0 );

  // Gauusian Pyramidの一番目の画像(σ = 2)の大きさを返す
  Coord< int > size() const { return( sz_ ); }

  /*
    処理結果を返す

    c : 座標値
    theta : 方向成分( 0 = 0° ; 1 = 45° ; 2 = 90° ; 3 = 135° )
  */

  // 顕著性
  double operator[]( Coord< int > c ) const
  { return( ( ! isValid( c ) ) ? 0 :
            saliency_[c.y][c.x] ); }
  // 輝度成分
  double lumi( Coord< int > c ) const
  { return( ( ! isValid( c ) ) ? 0 :
            iSum_[c.y][c.x] ); }
  // 色相成分
  double color( Coord< int > c ) const
  { return( ( ! isValid( c ) ) ? 0 :
            cSum_[c.y][c.x] ); }
  // 方向成分
  double ori( size_t theta, Coord< int > c ) const
  { return( ( ( ! isValid( c ) ) || ( theta >= O_SIZE ) ) ? 0 :
            oSum_[theta][c.y][c.x] ); }
  // 方向成分(各方向の和)
  double ori( Coord< int > c ) const
  { return( ( ! isValid( c ) ) ? 0 :
            oSumComb_[c.y][c.x] ); }
};

/*
  Combine : マップ data の要素の和を計算して res に返す
*/
void Combine( const Matrix< double >& data, Matrix< double >* res )
{
  for ( Matrix< double >::size_type ry = 0 ; ry < res->rows() ; ++ry ) {
    for ( Matrix< double >::size_type rx = 0 ; rx < res->cols() ; ++rx ) {
      Matrix< double >::size_type dx = round( rx * static_cast< double >( data.cols() ) / res->cols() );
      Matrix< double >::size_type dy = round( ry * static_cast< double >( data.rows() ) / res->rows() );
      if ( dx < data.cols() && dy < data.rows() )
        (*res)[ry][rx] += data[dy][dx];
    }
  }
}

/*
  CalcGlobalMaxMin : data から最大値 glbMax・最小値 glbMin を抽出する
*/
void CalcGlobalMaxMin( const Matrix< double >& data, double* glbMax, double* glbMin )
{
  *glbMax = std::numeric_limits< double >::min();
  *glbMin = std::numeric_limits< double >::max();

  for ( Matrix< double >::const_iterator cit = data.begin() ; cit != data.end() ; ++cit ) {
    double d = *cit;
    if ( *glbMin > d ) *glbMin = d;
    if ( *glbMax < d ) *glbMax = d;
  }
}

/*
  CalcLocalMax : 極大値の平均を求める

  data : マップの要素
  locMaxAve : 求める極大値
  glbMax : 求める最大値
*/
void CalcLocalMax( const Matrix< double >& data, double* locMaxAve, double* glbMax )
{
  typedef Matrix< Coord< double > >::size_type size_type;

  Matrix< Coord< double > > d1( data.rows(), data.cols() ); // 差分
  Matrix< Coord< double > > d2( data.rows(), data.cols() ); // 二階差分
  Coord< size_type > globalMaxCoord; // 最大値を取る座標
  vector< double > localMax;         // 抽出した極大値(局所最大値)

  *glbMax = std::numeric_limits< double >::min();

  // 差分の計算と最大値の抽出
  for ( size_type y1 = 0 ; y1 < data.rows() ; ++y1 ) {
    size_type y0 = ( y1 == 0 ) ? y1 : y1 - 1;
    size_type y2 = ( y1 == data.rows() - 1 ) ? y1 : y1 + 1;
    for ( size_type x1 = 0 ; x1 < data.cols() ; ++x1 ) {
      if ( *glbMax < data[y1][x1] ) {
        *glbMax = data[y1][x1];
        globalMaxCoord = Coord< size_type >( x1, y1 );
      }
      size_type x0 = ( x1 == 0 ) ? x1 : x1 - 1;
      size_type x2 = ( x1 == data.cols() - 1 ) ? x1 : x1 + 1;
      d1[y1][x1].x = data[y1][x2] - data[y1][x0];
      d1[y1][x1].y = data[y2][x1] - data[y0][x1];
    }
  }

  // 二階差分の計算
  for ( size_type y1 = 0 ; y1 < data.rows() ; ++y1 ) {
    size_type y0 = ( y1 == 0 ) ? y1 : y1 - 1;
    size_type y2 = ( y1 == data.rows() - 1 ) ? y1 : y1 + 1;
    for ( size_type x1 = 0 ; x1 < data.cols() ; ++x1 ) {
      size_type x0 = ( x1 == 0 ) ? x1 : x1 - 1;
      size_type x2 = ( x1 == data.cols() - 1 ) ? x1 : x1 + 1;
      d2[y1][x1].x = d1[y1][x2].x - d1[y1][x0].x;
      d2[y1][x1].y = d1[y2][x1].y - d1[y0][x1].y;
    }
  }

  // 極大値の抽出
  //   差分が0.01より小さく、二階差分が負ならば極大値
  for ( size_type y = 0 ; y < data.rows() ; ++y ) {
    for ( size_type x = 0 ; x < data.cols() ; ++x ) {
      // 最大値は対象外とする
      if ( globalMaxCoord.x == x && globalMaxCoord.y == y ) continue;
      if ( std::abs( d1[y][x].x ) < 0.01 && std::abs( d1[y][x].y ) < 0.01 ) {
        if ( d2[y][x].x < 0 && d2[y][x].y < 0 )
          localMax.push_back( data[y][x] );
      }
    }
  }

  // 極大値の平均値を算出
  *locMaxAve = std::accumulate( localMax.begin(), localMax.end(), double() );
  if ( localMax.size() > 0 )
    *locMaxAve /= localMax.size();
}

/*
  Normalize : data の正規化処理
*/
void Normalize( Matrix< double >* data )
{
  if ( data->cols() < 1 || data->rows() < 1 ) return;

  double glbMax = std::numeric_limits< double >::min(); // 最大値
  double glbMin = std::numeric_limits< double >::max(); // 最小値
  double locMaxAve; // 極大値の平均

  // 最大・最小値の抽出
  CalcGlobalMaxMin( *data, &glbMax, &glbMin );

  // 正規化処理(0～1の範囲の値にする)
  double diff = ( glbMax - glbMin > 0 ) ? glbMax - glbMin : 1.0; // 差分(全て等しいデータなら 1 とする)
  for ( Matrix< double >::iterator it = data->begin() ; it != data->end() ; ++it )
    *it = ( *it - glbMin ) / diff;

  // 極大値の平均を求める
  CalcLocalMax( *data, &locMaxAve, &glbMax );

  // (最大値 - 極大値平均)の二乗を各値に掛ける
  double sqDiff = std::pow( glbMax - locMaxAve, 2.0 );
  for ( Matrix< double >::iterator it = data->begin() ; it != data->end() ; ++it )
    *it *= sqDiff;
}

/*
  CleanUp : 最後に行う data の正規化処理

  値を 0 - 255 の範囲に再配分する。
*/
void CleanUp( Matrix< double >* data )
{
  double max = std::numeric_limits< double >::min(); // 最大値
  double min = std::numeric_limits< double >::max(); // 最小値

  // 最大値・最小値の取得
  CalcGlobalMaxMin( *data, &max, &min );

  // 正規化処理
  double diff = ( max - min > 0 ) ? max - min : 1.0; // 差分(全て等しいデータなら 1 とする)
  for ( Matrix< double >::iterator it = data->begin() ; it != data->end() ; ++it )
    *it = ( *it - min ) * UCHAR_MAX / diff;
}

/*
  CenterSurroundDiff : コンストラクタ

  pyramid : 処理に使うガウシアン・ピラミッド
  ki, kc, ko : 強度・色成分・方向成分の線形結合時の係数(重み)
*/
CenterSurroundDiff::CenterSurroundDiff( const GaussianPyramid& pyramid, double ki, double kc, double ko )
{
  typedef Matrix< double >::const_iterator const_iterator;
  typedef Matrix< double >::iterator iterator;

  // 差分の加算結果用配列などを初期化
  sz_ = pyramid.imageSize( 0 );
  iSum_.assign( sz_.y, sz_.x, 0.0 );
  cSum_.assign( sz_.y, sz_.x, 0.0 );
  for ( size_t o = 0 ; o < O_SIZE ; ++o )
    oSum_[o].assign( sz_.y, sz_.x, 0.0 );
  oSumComb_.assign( sz_.y, sz_.x );
  saliency_.assign( sz_.y, sz_.x );

  // Center-Surround Map の作成
  createCSMap( pyramid );

  // 和に対して正規化処理
  Normalize( &iSum_ );
  Normalize( &cSum_ );
  for ( size_t o = 0 ; o < O_SIZE ; ++o ) {
    Normalize( &oSum_[o] );
    Combine( oSum_[o], &oSumComb_ );
  }
  Normalize( &oSumComb_ );

  // 各成分を線形結合
  iterator s = saliency_.begin();
  const_iterator i = iSum_.begin();
  const_iterator c = cSum_.begin();
  const_iterator o = oSumComb_.begin();
  for ( ; s != saliency_.end() ; ++s, ++i, ++c, ++o )
    *s = *i * ki + *c * kc + *o * ko;

  CleanUp( &saliency_ );
}

/*
  CenterSurroundDiff::createCSMap : Center-Surround Map の作成

  pyramid : 処理に使うガウシアン・ピラミッド
*/
void CenterSurroundDiff::createCSMap( const GaussianPyramid& pyramid )
{
  // 差分計算用配列
  Matrix< double > lumi;        // 強度 I(c,s)
  Matrix< double > rg;          // red - green = RG(c,s)
  Matrix< double > by;          // blue - yellow = BY(c,s)
  Matrix< double > ori[O_SIZE]; // 方向成分 = O(c,s,θ)

  // Centerマップ
  // c = { 0, 1, 2 } (実際は 1/4,1/8,1/16 のマップ)
  for ( size_t c = 0 ; c < 3 ; ++c ) {

    // 差分計算用配列をマップの大きさに初期化
    Coord< int > sz = pyramid.imageSize( c );
    lumi.assign( sz.y, sz.x );
    rg.assign( sz.y, sz.x );
    by.assign( sz.y, sz.x );
    for ( size_t o = 0 ; o < O_SIZE ; ++o )
      ori[o].assign( sz.y, sz.x );

    // Surroundマップ
    // s = { c + 3, c + 4 } (実際は 1/8c, 1/16c のマップ)
    for ( size_t s = c + 3 ; s < c + 5 ; ++s ) {
      double scale = std::pow( 2.0, s - c ); // Center-Surround の大きさの比率

      // 差分の計算
      for ( Coord< int > cc( 0, 0 ) ; cc.y < sz.y ; ++( cc.y ) ) {
        for ( cc.x = 0 ; cc.x < sz.x ; ++( cc.x ) ) {

          Coord< int > sc( round( cc.x / scale ), round( cc.y / scale ) ); // Surroundマップ上の座標

          // I(c,s) = | I(c) (-) I(s) |
          lumi[cc.y][cc.x] =
            std::abs( pyramid.lumi( c, cc ) - pyramid.lumi( s, sc ) );
          // RG(c,s) = | ( R(c) - G(c) ) (-) ( G(s) - R(s) ) |
          rg[cc.y][cc.x] =
            std::abs( ( pyramid.red( c, cc ) - pyramid.green( c, cc ) ) -
                      ( pyramid.green( s, sc ) - pyramid.red( s, sc ) ) );
          // BY(c,s) = | ( B(c) - Y(c) ) (-) ( Y(s) - B(s) ) |
          by[cc.y][cc.x] =
            std::abs( ( pyramid.blue( c, cc ) - pyramid.yellow( c, cc ) ) -
                      ( pyramid.yellow( s, sc ) - pyramid.blue( s, sc ) ) );
          // O(c,s,θ) = | O(c,θ) (-) O(s,θ) |
          for ( unsigned int o = 0 ; o < O_SIZE ; ++o )
            ori[o][cc.y][cc.x] =
              std::abs( pyramid.ori( c, o, cc ) - pyramid.ori( s, o, sc ) );
        }
      }

      // 正規化処理と和の計算
      Normalize( &lumi );
      Combine( lumi, &iSum_ );

      Normalize( &rg );
      Normalize( &by );
      Combine( rg, &cSum_ );
      Combine( by, &cSum_ );

      for ( unsigned int o = 0 ; o < O_SIZE ; ++o ) {
        Normalize( &ori[o] );
        Combine( ori[o], &oSum_[o] );
      }
    }
  }
}

各スケール画像間の差分計算など、主な処理はほとんど全てメンバ関数 createCSMap が行っています。引数として GaussianPyramid インスタンスがそのまま渡され、各スケール画像の成分を利用して処理を行います。
処理の内容は今まで説明したものがそのままコーディングされているだけですが、正規化処理を行う関数 Normalize について少し補足しておきます。Normalize では最初にマップ上から最大値と最小値を抽出しています。この二つの値を使って、各値が 0 から 1 までの値を取るように正規化処理を行います。次に CalcLocalMax を使って極大値 ( 局所的な最大値 ) の平均を抽出し、最大値と極大値平均の平方値を各値に掛けて処理が完了します。
最大・最小値抽出後の正規化処理によって最大値が変化するため、CalcLocalMax では極小値平均だけでなく最大値も抽出しています。最大値は通常 1 になるため、この処理は本来不要なのですが、極小値を求める時に最大値の座標も抽出しているので、ついでに処理を追加しています。
CalcLocalMax で極大値平均を求めるところでは、各ピクセルの持つ値の差分を使っています。差分については「画像処理 (1) シーム・カービング」の中でも簡単に説明してありますが、サンプル・プログラムでは「中央差分 ( Central Difference )」として以下の式で差分を求めています。

ここで、d0( X, Y ) は ( X, Y ) におけるマップ上の値、d1_X( X, Y ) と d1_Y( X, Y ) はそれぞれ x, y 方向における ( X, Y ) 上の差分、すなわち値の変化量を表しています。正しくは、座標値の差である 2 で徐算する必要がありますが、値そのものを利用するのではないため省略しています。

こうして求めた差分を使って、さらに同様の処理により二階差分 d2 を求めます。

差分 d1 は、微分係数の近似値と見なすことができます。同様に、d2 は二階微分係数の近似値になります。微分係数がゼロの点は極値点か変曲点のいずれかになり、極大値を取る点では二階微分係数は負数になるので、d1 がゼロで d2 が負数となる点を極大値として抽出します。d1 が完全にゼロになることはあり得ないので、サンプル・プログラムではかなり大雑把に、0.01 より小さい場合はゼロと見なしています。

正規化処理を行う関数としてもう一つ CleanUp があります。こちらは顕著性マップを画面に出力するときのために値を 0 ～ 255 の範囲に正規化する処理を行っているだけで、顕著性マップを求める処理そのものにはあまり関係していないものです。

サンプル・プログラムを使って顕著性マップを出力した結果を以下に示します。使用した画像は Kodak Lossless True Color Image Suite から拝借した以下のテスト用サンプルです。なお、実際の画像のサイズは 768 x 512 です。

まずは、スケール画像を出力した結果です。ガウシアン・フィルタにおける σ の値は 1.0、ガボール・フィルタの各パラーメータは、波長 λ = 8、アスペクト比 γ = 1、位相のオフセット φ = 0、帯域幅 b = 1、畳み込み積分を行う窓サイズ -8 ～ 8 としています。以下には、スケール σ = 2 ( 1 / 4 縮尺 ) の場合と σ = 5 ( 1 / 32 縮尺 ) の場合の二種類を示しています。なお、スケールが異なれば画像のサイズもバラバラになるのですが、見やすさを優先して全ての画像サイズを合わせています。

色相成分 ( 特に赤 ) において、ドアや、窓にある開き戸の部分がかなり強調されていることが分かります。

各成分毎にスケール画像間の差分を求め、それらを重ね合わせた結果が下図になります。

やはり色相成分において、ドアや開き戸の部分が強調されています。これらをさらに重ね合わせると、顕著性マップを得ることができます。

このテスト画像は確かに、ドアや開き戸の部分に目が移りやすいようにみえるので、割と正しく表現できているのではないかと思います。

顕著性マップに関する情報は University of Southern California の iLab Homepage にある Bottom-Up Visual Attention Home Page から得ることができます。この中にあるサンプル画像を使って顕著性マップを出力してみたいと思います。テスト画像は下図にあるようなもので、緑色のパターンの中に一つだけ赤色のパターンがあって、この赤色のパターンが最も注意を引く箇所になります。

処理条件は先と全く同じです。スケール画像と、その差分の重ね合わせ処理結果を以下に示します。

輝度の方が色相成分よりも強く出力されてしまい、顕著性マップは下図に示すように全パターンに高く現れる結果となりました。シーム・カービング処理に応用するだけなのであればこれでも充分なのかも知れませんが、注目度の最も高い箇所を求める必要がある場合はまだ感度が低いと言わざるを得ません。スケール画像を見ると、画素の細かい画像と粗い画像で色相成分の分布 ( 特に赤・緑の差分 ) に差が見られないので、これがスケール画像間の差分において顕著性が現れない原因の一つとなっていることが考えられます。

ガウシアン・フィルタの σ を 4.0 にして、各画素がより広範囲のピクセルの影響を受けるようにしてみます。

今度は、色相成分においてスケール画像間での差が生じることで、その差分で赤色パターンの箇所に強い顕著性が発生し、下図のように、ある程度改善された顕著性マップを得ることができました。しかし、輝度成分は全体的にもう少し低くなるべきであり、ここは正規化処理の改善が必要だと考えられます。

最後に、顕著性マップを使ったシーム・カービング処理を行ってみたいと思います。サンプル画像も、前回と同じもの(下図)を使います。

下図が処理した結果で、左側が輝度による顕著性、中央がガボール・フィルタによる顕著性、右側が顕著性マップを使った時の画像です。

図 3-16. 顕著性マップを使ったシーム・カービング削除処理の結果
Splash
輝度	ガボール・フィルタ	顕著性マップ

Sailboat on lake
輝度	ガボール・フィルタ	顕著性マップ

以前の方法と比較すると、Splash では左右の端にある黒い部分が残る傾向が見られます。確かに、この部分は注意を引きやすい部分となるかもしれません。しかし、Sailboat の方では木の幹の部分がかなり削られてしまったことが気になります。他の箇所と比べると目立たない箇所と判定されたということでしょうか。帆船の部分はほとんど残っており、顕著性マップが最もきれいに残っているという結果になりました。
下図は、テスト画像の顕著性マップです。このように、注意を引かない ( 目立たない ) ところが低いエネルギーであれば、顕著性マップを利用したシーム・カービング処理もかなり有効であるといえます。

顕著性マップを利用したシーム・カービング処理のためのサンプル・プログラムを以下に示します。

/******************************************
   Saliency Map を使った顕著性データ
******************************************/
class ImageSaliency_SaliencyMap
{
  GaussianPyramid gp_;     // ガウシアン・ピラミッド
  CenterSurroundDiff csd_; // スケール間差分

  Coord< int > drawSize_;  // 元画像の大きさ

 public:

  /*
    IS_SaliencyMap コンストラクタ

    draw : 対象の画像オブジェクト
    sigma : ガウシアン・フィルタの σ
    lambda : ガボール・フィルタの波長
    ki, kc, ko : 輝度・色相・方向成分の線形結合時の係数
  */
  ImageSaliency_SaliencyMap( const DrawingArea_IF& draw,
                             double sigma, double lambda, double ki, double kc, double ko )
    : gp_( draw, sigma, lambda ), csd_( gp_, ki, kc, ko ), drawSize_( draw.size() ) {}

  // c の位置の顕著性データを求める
  double operator()( Coord< int > c ) const;
};

/*
  ImageSaliency_saliencyMap::operator() : c の位置の顕著性データを求める
*/
double ImageSaliency_SaliencyMap::operator()( Coord< int > c ) const
{
  Coord< int > csdSize = csd_.size(); // Center-Surround差分(σ = 2)の画像サイズ

  // Center-Surround差分上の座標値
  Coord< double > p( static_cast< double >( c.x * csdSize.x ) / drawSize_.x,
                     static_cast< double >( c.y * csdSize.y ) / drawSize_.y ); // double型
  Coord< int > c2( p.x, p.y ); // int型(格子の左上隅を示す)

  // 各格子の顕著性を取得
  double q11 = csd_[c2];                               // 左上
  double q21 = csd_[Coord< int >( c2.x + 1, c2.y )];     // 右上
  double q12 = csd_[Coord< int >( c2.x, c2.y + 1 )];     // 左下
  double q22 = csd_[Coord< int >( c2.x + 1, c2.y + 1 )]; // 右下

  // 線形補間処理
  return( ( ( c2.x + 1 ) - p.x ) * ( ( c2.y + 1 ) - p.y ) * q11 +
          ( p.x - c2.x ) * ( ( c2.y + 1 ) - p.y ) * q21 +
          ( ( c2.x + 1 ) - p.x ) * ( p.y - c2.y ) * q12 +
          ( p.x - c2.x ) * ( p.y - c2.y ) * q22
          );
}

/*
  CreateBySaliencyMap : Saliency Map を使った顕著性データへの変換

  draw : 対象の画像
  saliency : 変換した顕著性データを登録する変数へのポインタ
  sigma : ガウシアン・フィルタの σ
  lambda : ガボール・フィルタの波長
  ki, kc, ko : 輝度・色相・方向成分の線形結合時の係数
*/
void CreateBySaliencyMap( const DrawingArea_IF& draw, Matrix< double >* saliency,
                          double sigma, double lambda, double ki, double kc, double ko )
{
  ImageSaliency_SaliencyMap isSMap( draw, sigma, lambda, ki, kc, ko );
  Create( saliency, draw.size(), isSMap );
}

ImageSaliency_SaliencyMap は顕著性マップを利用してシーム・カービングに必要な顕著性を求めるためのクラスで、コンストラクタで、先に説明したガウシアン・ピラミッドと顕著性マップ用のインスタンスを構築しています。シーム・カービング処理では各ピクセルの持つ値を顕著性として、その画素間の差分をエネルギーとして求めていました。この章で「顕著性マップ」と呼んでいる値はシーム・カービングでの顕著性と同じ意味になるので、CenterSurroundDiff オブジェクトで求められた成分がそのまま利用できます。しかし、原画像とその顕著性マップでは画像のサイズが異なり、後者の方が 1/4 ほど縮尺されているので、原画像のサイズに合わせて補間を行う必要があります。この時、最も近くにあるピクセルの値をそのまま使う、いわゆる「最近傍法 ( Nearest Neighbour )」を利用すると、エネルギーを求める時に、差分がゼロの箇所が大量に発生してしまい正しく処理することができないため、ここでは「線形補間法 ( Bilinear Interpolation )」を利用して補間処理を行っています。サンプル補間については「グラフィック・パターンの扱い (5) サンプル補間」で詳しく説明してあります。

前述の通り、顕著性マップは 1998 年に Itti らによって発表されたモデルで、人の視覚の機能を元に考案されたものです。人が視覚情報の中から何に注目するかについては、たくさんの要素が相互に作用しているのですが、それらを分類すると二つの要素、「ボトムアップ ( Bottom-Up )」と「トップダウン ( Top-Down )」になるとしています。前者は入力された視覚情報によって瞬時に判断されるもので、内部器官の状態に依存しない要素、後者は内部器官の状態に依存するもので、いわゆる記憶や経験などによって影響される要素です。今回の仕組みは、当然前者の「ボトムアップ」だけを取り上げています。
今回は、シーム・カービング処理のための一手法として利用してみましたが、文献ではその他の応用例として、道路標識などの自動認識や、注意を引くデザインの評価などを挙げています。また、他にもオブジェクトを抽出する処理やロボットの視覚情報処理への応用など、様々な場面で応用できそうな技術です。

∫{-∞ → ∞}∫{-∞ → ∞} exp( - ( x² + y² ) / 2σ² ) dxdy	=	∫{-∞ → ∞}∫{-∞ → ∞} exp( -x² / 2σ² ) exp( -y² / 2σ² ) dxdy
	=	∫{-∞ → ∞} exp( -x² / 2σ² ) dx ∫{-∞ → ∞} exp( -y² / 2σ² ) dy
	=	2πσ²

Luminance	Orientation

Red	Green	Red - Green

Blue	Yellow	Blue - Yellow

Luminance	Orientation

Red	Green	Red - Green

Blue	Yellow	Blue - Yellow

Luminance	Orientation

Red	Green	Red - Green

Blue	Yellow	Blue - Yellow

Luminance	Orientation

Red	Green	Red - Green

Blue	Yellow	Blue - Yellow

O( σ, θ )	=	I(σ) * ψ(θ)
	=	Σ_m{-h_m → h_m}Σ_n{-h_n → h_n}( I( σ, m, n ) ψ( θ, x - m, y - n ) )