Incremental PCA の実装

Candid Covariance-free Incremental Principal Component Analysis を OpenCV で、実装してみました。結果はいかほどに。アルゴリズムはこんな感じ。

 分析する画像。ただし、
 主要ベクトル。
 忘却の値。論文中では２〜４が良いと書いてある。
そして、から画像を与えていき、以下を実行する。
1. 
2. 
2.1  のようにi番目の主要ベクトル初期化する
2.2 
2.2.1 
2.2.2 

ソースも、参考までに。

IplImage** ForeExtractUsingIPCA::m_ArrImgMain; //!<主成分画像の配列
//!<ただし、[0]には擬似平均画像が入る、[1]からが主成分。
IplImage** ForeExtractUsingIPCA::m_ArrImgTmp; //!<一時画像の配列

void ForeExtractUsingIPCA::UpdateImgMain(IplImage *imgIn) //!< GRB画像を入力
{
  ///// 初期化 uchar -> float 変換
  cvConvert( imgIn, m_ArrImgTmp[ 0 ] );

  int iterateMax = m_MainMax -1;
  if( m_MainMax > m_FrameNow ) iterateMax = m_FrameNow;

  IplImage* imgA = cvCreateImage( cvSize(m_Width, m_Height), IPL_DEPTH_32F, 3 );
  IplImage* imgB = cvCreateImage( cvSize(m_Width, m_Height), IPL_DEPTH_32F, 3 );
  IplImage* imgC = cvCreateImage( cvSize(m_Width, m_Height), IPL_DEPTH_32F, 3 );
  double scalerA, scalerB, scalerC;
  double nrmV;
  double l = m_Amnesic; //!<忘却の値、２ぐらいがよい

  for( int i=0; i<=iterateMax; i++ )
  {
    if( i == m_FrameNow )
    {
      cvConvert( m_ArrImgTmp[ m_FrameNow ], m_ArrImgMain[ m_FrameNow ] );
      continue;
    }

    ///// Vi(n) = [a= (n-1-l)/n * Vi(n-1)] + [b= (1+l)/n * Ui(n)T Vi(n-1)/|Vi(n-1)| * Ui(n) ]
    nrmV = cvNorm( m_ArrImgMain[i], NULL, CV_L2 );

    scalerA = (double)(m_FrameNow-1-l) / (double)m_FrameNow;
    cvScale( m_ArrImgMain[i], imgA, scalerA );

    double dotUV = cvDotProduct( m_ArrImgTmp[i], m_ArrImgMain[i] );
    scalerB = ((double)(1+l) * dotUV) / ((double)m_FrameNow * nrmV );
    cvScale( m_ArrImgTmp[i], imgB, scalerB );

    cvAdd( imgA, imgB, m_ArrImgMain[i] );

    ///// Ui+1(n) = Ui(n) - [c= Ui(n)T Vi(n)/|Vi(n)| * Vi(n)/|Vi(n)| ]
    if( i == iterateMax || i >= m_MainMax -1 ) continue;

    nrmV = cvNorm( m_ArrImgMain[i], NULL, CV_L2 );
    dotUV = cvDotProduct( m_ArrImgTmp[i], m_ArrImgMain[i] );
    scalerC = dotUV / (nrmV * nrmV);
    cvScale( m_ArrImgMain[i], imgC, scalerC );
    cvSub( m_ArrImgTmp[i], imgC, m_ArrImgTmp[i+1] );
  }

  cvReleaseImage( &imgA );
  cvReleaseImage( &imgB );
  cvReleaseImage( &imgC );

  m_FrameNow++;
}

PCA で得られた主要ベクトルとの距離の比較をすると以下のように。
なお、画像は10080枚与えました。

いまいち収束性がなぁ。