OpenCV-EMDアルゴリズムを勉強します.



 
//adapted by the example of leanring opencv by mikewolf
//2008-11-26
#include "stdafx.h"
#include <cv.h>
#include <highgui.h>

int _tmain(int argc, _TCHAR* argv[])
{
   IplImage* src, *src1;


   src=cvLoadImage("lena.jpg", 1);
   src1=cvLoadImage("lena3.jpg", 1);

   // Compute the HSV image and decompose it into separate planes.
   //
   IplImage* hsv = cvCreateImage( cvGetSize(src), 8, 3 );
   IplImage* hsv1 = cvCreateImage( cvGetSize(src), 8, 3 );
   cvCvtColor( src, hsv, CV_BGR2HSV );
   cvCvtColor( src1, hsv1, CV_BGR2HSV );
   IplImage* h_plane = cvCreateImage( cvGetSize(src), 8, 1 );
   IplImage* s_plane = cvCreateImage( cvGetSize(src), 8, 1 );
   IplImage* v_plane = cvCreateImage( cvGetSize(src), 8, 1 );
   IplImage* h_plane1 = cvCreateImage( cvGetSize(src), 8, 1 );
   IplImage* s_plane1 = cvCreateImage( cvGetSize(src), 8, 1 );
   IplImage* v_plane1 = cvCreateImage( cvGetSize(src), 8, 1 );
   IplImage* planes[] = { h_plane, s_plane,v_plane };
   IplImage* planes1[] = { h_plane1, s_plane1,v_plane1 };
   cvCvtPixToPlane( hsv, h_plane, s_plane, v_plane, 0 );
   cvCvtPixToPlane( hsv1, h_plane1, s_plane1, v_plane1, 0 );
   // Build the histogram and compute its contents.
   //     h_bins*s_bins*v_bins     ,   int h_bins = 32, s_bins = 30, v_bins = 8  ,        
   int h_bins = 8, s_bins = 5, v_bins = 5;  
   CvHistogram* hist, *hist1;
   {
      int hist_size[] = { h_bins, s_bins ,v_bins};
      float h_ranges[] = { 0, 180 }; // hue is [0,180]
      float s_ranges[] = { 0, 255 };
      float v_ranges[] = { 0, 255 };
      float* ranges[] = { h_ranges, s_ranges ,v_ranges};
      hist = cvCreateHist(
         3,
         hist_size,
         CV_HIST_ARRAY,
         ranges,
         1
         );
      hist1 = cvCreateHist(
         3,
         hist_size,
         CV_HIST_ARRAY,
         ranges,
         1
         );
   }
   cvCalcHist( planes, hist, 0, 0 ); //Compute histogram
   cvNormalizeHist( hist, 1.0 ); //Normalize it
   cvCalcHist( planes1, hist1, 0, 0 ); //Compute histogram
   cvNormalizeHist( hist1, 1.0 ); //Normalize it

   CvMat* sig1,*sig2;
   int numrows = h_bins*s_bins*v_bins;
   //Create matrices to store the signature in
   //
   sig1 = cvCreateMat(numrows, 4, CV_32FC1); 
   sig2 = cvCreateMat(numrows, 4, CV_32FC1); 
   int h,s,v;
   for( h = 0; h < h_bins; h++ ) 
   {
      for( s = 0; s < s_bins; s++ )
      {
         for( v=0; v < v_bins; v++)
         {
            float bin_val = cvQueryHistValue_3D( hist, h, s,v );
            
            cvSet2D(sig1,h*s_bins*v_bins + s*v_bins+v,0,cvScalar(bin_val)); //bin value
            cvSet2D(sig1,h*s_bins*v_bins + s*v_bins+v,1,cvScalar(h)); //Coord 1
            cvSet2D(sig1,h*s_bins*v_bins + s*v_bins+v,2,cvScalar(s)); //Coord 2
            cvSet2D(sig1,h*s_bins*v_bins + s*v_bins+v,3,cvScalar(v)); //Coord 3

            bin_val = cvQueryHistValue_3D( hist1, h, s,v );
            cvSet2D(sig2,h*s_bins*v_bins + s*v_bins+v,0,cvScalar(bin_val)); //bin value
            cvSet2D(sig2,h*s_bins*v_bins + s*v_bins+v,1,cvScalar(h)); //Coord 1
            cvSet2D(sig2,h*s_bins*v_bins + s*v_bins+v,2,cvScalar(s)); //Coord 2
            cvSet2D(sig2,h*s_bins*v_bins + s*v_bins+v,3,cvScalar(v)); //Coord 3

         }
      }
   }
    
   float emd = cvCalcEMD2(sig1,sig2,CV_DIST_L2);
   printf("    :%3.1f%%",(1-emd)*100);

   cvNamedWindow( "Source1", 1 );
   cvMoveWindow("Source1", 0, 0);
   cvShowImage( "Source1", src );
   cvNamedWindow( "Source2", 1 );
   cvMoveWindow("Source2", 550, 0);
   cvShowImage( "Source2", src1 );
   cvWaitKey(0);

   return 0;
}


CalcEMD 2つの重み付けポイントセットの間で最小動作距離float cvCalcEMD 2を計算します.                                  Cv Distance Function distanc=NULL、const CvArr*costit_maturix=NULL、                                  CvArr*flow=NULL、float*lowerband=NULL、void*userdata=NULL);                                  typedef float(*Cv Distance Function)(const float*f 1,const float*f 2,void*userdata)
signature 1最初の署名のサイズはsize 1です.×(dims+1)の浮動小数マトリクスは、各ラインにポイントの重みとポイントの座標を順次記憶する.行列は、一列のみ(すなわち、重みのみ)を許可します.ユーザー定義の代価行列を使用する場合は.signature 2番目の署名は、signature 1のフォーマットと同じです.×(dims+1)行数は異なるが(列数は同じ).追加の仮想ポイントが、signature 1またはsignature 2に参加する場合、重みも異なることがある.distancetype使用の基準、CV_DIST_L 1,CV_DIST_L 2とCV_DIST_Cはそれぞれ標準的な基準である.CV_DIST_USERとは、ユーザーカスタム関数distance(u)を使用することを意味します.funcまたは事前に計算した代価行列costt_matixdistancefuncユーザがカスタマイズした距離関数.2点の座標で2点間の距離を計算します.コストカmatrixカスタムサイズはsize 1です.×size 2の対価行列コストカmatixとdistance_funcの両方は少なくとも一つはNULLでなければなりません.そして、もし対価関数を使うなら、下の境界は計算できません.基準関数が必要です.flowのサイズはsize 1です.×size 2フロー行列(flow matrix):flow i jは、signature 1のi番目の点からsignature 2のj番目の点までの流れです.lower.boundオプションの入出力パラメータ:2つの署名の間の距離の下で境界は、2つのセントロイドの間の距離です.ユーザー定義の代価行列を使用すると、ポイントセットの所有権が異なるか、または署名がある場合には重みだけが含まれます.すなわち、署名行列は単独の列だけです.ユーザーは*lowerを初期化しなければなりません.bound.セントロイド間の距離が取得イコール*lower_bound(これは署名間が十分遠いという意味です.関数はEMDを計算しません.いずれの場合も関数が戻るときには*lowerbandは計算されたセントロイド距離に設定されます.したがって、ユーザーがセントロイド距離とT EMDを同時に計算したい場合、*lowerbandは0に設定されます.userdataがカスタム距離関数に送信するオプションのデータポインタ関数cvCalcEMD 2は、2つの重み付きポイントを計算します.セット間の移動距離や距離は、下界にあります.「RubnerSept 98」に記載されているアプリケーションの一つは、画像抽出の多次元ヒストグラム比較です.EMDは、ある種類の単純形アルゴリズム(simplex algorithm)を使用しています.これによって解決される交通問題です.複雑さは最悪の場合は指数形式ですが、平均的にはかなり速いです.実の基準に対して、下境界の計算はより速く(線形時間アルゴリズムを使用して)できます.2つの点セットが十分に遠いかどうかを大まかに判断するために使用できます.同じ目標に接続できないほどです.