OpenCV每日函数 直方图相关(3)

一、EMD函数

1、函数原型

        该函数计算推土机距离和/或两个加权点配置之间距离的下边界。 应用之一是用于图像检索的多维直方图比较。 EMD 是一个使用单纯形算法的一些修改来解决的运输问题，因此在最坏的情况下复杂度是指数级的，但平均而言它要快得多。 在真实度量的情况下，可以更快地计算下边界（使用线性时间算法），它可以用来粗略地确定两个签名是否足够远，以至于它们不能与同一个对象相关。

float 	cv::EMD (InputArray signature1, InputArray signature2, int distType, InputArray cost=noArray(), float *lowerBound=0, OutputArray flow=noArray())

2、参数详解

signature1	第一个签名，一个 size1×dims+1 浮点矩阵。 每行存储点权重，后跟点坐标。 如果使用用户定义的成本矩阵，则允许该矩阵具有单列（仅权重）。 权重必须是非负的并且至少有一个非零值。
signature2	与 signature1 格式相同的第二个签名，尽管行数可能不同。 总重量可能不同。 在这种情况下，一个额外的“虚拟”点被添加到签名 1 或签名 2。 权重必须是非负的并且至少有一个非零值。
distType	使用的公制。 请参阅距离类型。
cost	用户定义的 size1×size2 成本矩阵。 此外，如果使用成本矩阵，则无法计算下边界 lowerBound 因为它需要一个度量函数。
lowerBound	可选输入/输出参数：两个签名之间距离的下边界，即质心之间的距离。 如果使用用户定义的成本矩阵，点配置的总权重不相等，或者签名仅由权重组成（签名矩阵只有一列），则可能无法计算下边界。 您必须**初始化 *lowerBound 。 如果计算的质心之间的距离大于或等于 *lowerBound（这意味着签名足够远），则该函数不计算 EMD。 在任何情况下，*lowerBound 都设置为返回时计算的质心之间的距离。 因此，如果要计算质心和 EMD 之间的距离，*lowerBound 应设置为 0。
flow	得到的 size1×size2 流矩阵： 是从 i-th 签名 1 的点到 j-th 签名 2 的点的流。

3、OpenCV源码

（1）源码路径

opencv\modules\imgproc\src\emd.cpp

（2）源码代码

CV_IMPL float cvCalcEMD2( const CvArr* signature_arr1,
            const CvArr* signature_arr2,
            int dist_type,
            CvDistanceFunction dist_func,
            const CvArr* cost_matrix,
            CvArr* flow_matrix,
            float *lower_bound,
            void *user_param )
{
    cv::AutoBuffer<char> local_buf;
    CvEMDState state;
    float emd = 0;

    memset( &state, 0, sizeof(state));

    double total_cost = 0;
    int result = 0;
    float eps, min_delta;
    CvNode2D *xp = 0;
    CvMat sign_stub1, *signature1 = (CvMat*)signature_arr1;
    CvMat sign_stub2, *signature2 = (CvMat*)signature_arr2;
    CvMat cost_stub, *cost = &cost_stub;
    CvMat flow_stub, *flow = (CvMat*)flow_matrix;
    int dims, size1, size2;

    signature1 = cvGetMat( signature1, &sign_stub1 );
    signature2 = cvGetMat( signature2, &sign_stub2 );

    if( signature1->cols != signature2->cols )
        CV_Error( CV_StsUnmatchedSizes, "The arrays must have equal number of columns (which is number of dimensions but 1)" );

    dims = signature1->cols - 1;
    size1 = signature1->rows;
    size2 = signature2->rows;

    if( !CV_ARE_TYPES_EQ( signature1, signature2 ))
        CV_Error( CV_StsUnmatchedFormats, "The array must have equal types" );

    if( CV_MAT_TYPE( signature1->type ) != CV_32FC1 )
        CV_Error( CV_StsUnsupportedFormat, "The signatures must be 32fC1" );

    if( flow )
    {
        flow = cvGetMat( flow, &flow_stub );

        if( flow->rows != size1 || flow->cols != size2 )
            CV_Error( CV_StsUnmatchedSizes,
            "The flow matrix size does not match to the signatures' sizes" );

        if( CV_MAT_TYPE( flow->type ) != CV_32FC1 )
            CV_Error( CV_StsUnsupportedFormat, "The flow matrix must be 32fC1" );
    }

    cost->data.fl = 0;
    cost->step = 0;

    if( dist_type < 0 )
    {
        if( cost_matrix )
        {
            if( dist_func )
                CV_Error( CV_StsBadArg,
                "Only one of cost matrix or distance function should be non-NULL in case of user-defined distance" );

            if( lower_bound )
                CV_Error( CV_StsBadArg,
                "The lower boundary can not be calculated if the cost matrix is used" );

            cost = cvGetMat( cost_matrix, &cost_stub );
            if( cost->rows != size1 || cost->cols != size2 )
                CV_Error( CV_StsUnmatchedSizes,
                "The cost matrix size does not match to the signatures' sizes" );

            if( CV_MAT_TYPE( cost->type ) != CV_32FC1 )
                CV_Error( CV_StsUnsupportedFormat, "The cost matrix must be 32fC1" );
        }
        else if( !dist_func )
            CV_Error( CV_StsNullPtr, "In case of user-defined distance Distance function is undefined" );
    }
    else
    {
        if( dims == 0 )
            CV_Error( CV_StsBadSize,
            "Number of dimensions can be 0 only if a user-defined metric is used" );
        user_param = (void *) (size_t)dims;
        switch (dist_type)
        {
        case CV_DIST_L1:
            dist_func = icvDistL1;
            break;
        case CV_DIST_L2:
            dist_func = icvDistL2;
            break;
        case CV_DIST_C:
            dist_func = icvDistC;
            break;
        default:
            CV_Error( CV_StsBadFlag, "Bad or unsupported metric type" );
        }
    }

    result = icvInitEMD( signature1->data.fl, size1,
                        signature2->data.fl, size2,
                        dims, dist_func, user_param,
                        cost->data.fl, cost->step,
                        &state, lower_bound, local_buf );

    if( result > 0 && lower_bound )
    {
        emd = *lower_bound;
        return emd;
    }

    eps = CV_EMD_EPS * state.max_cost;

    /* if ssize = 1 or dsize = 1 then we are done, else ... */
    if( state.ssize > 1 && state.dsize > 1 )
    {
        int itr;

        for( itr = 1; itr < MAX_ITERATIONS; itr++ )
        {
            /* find basic variables */
            result = icvFindBasicVariables( state.cost, state.is_x,
                                            state.u, state.v, state.ssize, state.dsize );
            if( result < 0 )
                break;

            /* check for optimality */
            min_delta = icvIsOptimal( state.cost, state.is_x,
                                      state.u, state.v,
                                      state.ssize, state.dsize, state.enter_x );

            if( min_delta == CV_EMD_INF )
                CV_Error( CV_StsNoConv, "" );

            /* if no negative deltamin, we found the optimal solution */
            if( min_delta >= -eps )
                break;

            /* improve solution */
            if(!icvNewSolution( &state ))
                CV_Error( CV_StsNoConv, "" );
        }
    }

    /* compute the total flow */
    for( xp = state._x; xp < state.end_x; xp++ )
    {
        float val = xp->val;
        int i = xp->i;
        int j = xp->j;

        if( xp == state.enter_x )
          continue;

        int ci = state.idx1[i];
        int cj = state.idx2[j];

        if( ci >= 0 && cj >= 0 )
        {
            total_cost += (double)val * state.cost[i][j];
            if( flow )
                ((float*)(flow->data.ptr + flow->step*ci))[cj] = val;
        }
    }

    emd = (float) (total_cost / state.weight);
    return emd;
}

float cv::EMD( InputArray _signature1, InputArray _signature2,
               int distType, InputArray _cost,
               float* lowerBound, OutputArray _flow )
{
    CV_INSTRUMENT_REGION();

    Mat signature1 = _signature1.getMat(), signature2 = _signature2.getMat();
    Mat cost = _cost.getMat(), flow;

    CvMat _csignature1 = cvMat(signature1);
    CvMat _csignature2 = cvMat(signature2);
    CvMat _ccost = cvMat(cost), _cflow;
    if( _flow.needed() )
    {
        _flow.create(signature1.rows, signature2.rows, CV_32F);
        flow = _flow.getMat();
        flow = Scalar::all(0);
        _cflow = cvMat(flow);
    }

    return cvCalcEMD2( &_csignature1, &_csignature2, distType, 0, cost.empty() ? 0 : &_ccost,
                       _flow.needed() ? &_cflow : 0, lowerBound, 0 );
}

float cv::wrapperEMD(InputArray _signature1, InputArray _signature2,
               int distType, InputArray _cost,
               Ptr<float> lowerBound, OutputArray _flow)
{
    return EMD(_signature1, _signature2, distType, _cost, lowerBound.get(), _flow);
}

4、参考代码

        EMD（earth mover distance）方法是比较图像相似度的一种很好的方法。但是处理时间很慢。为了使用 EMD 比较，我们应该制作签名值。EMD 方法比较两个签名值。

        首先，我们准备 2 张图像的直方图。并将直方图的值转换为签名。

#include <iostream>
#include <vector>

#include <stdio.h>      
#include <opencv2\opencv.hpp>    


#ifdef _DEBUG           
#pragma comment(lib, "opencv_core249d.lib")   
#pragma comment(lib, "opencv_imgproc249d.lib")   //MAT processing   
#pragma comment(lib, "opencv_highgui249d.lib")   
#else   
#pragma comment(lib, "opencv_core249.lib")   
#pragma comment(lib, "opencv_imgproc249.lib")      
#pragma comment(lib, "opencv_highgui249.lib")   
#endif   


using namespace cv;   
using namespace std;   
  
int main()   
{   

 //read 2 images for histogram comparing   
 ///   
 Mat imgA, imgB;   
 imgA = imread(".\\image1.jpg");   
 imgB = imread(".\\image2.jpg");   


 imshow("img1", imgA);
 imshow("img2", imgB);


 //variables preparing   
 ///   
 int hbins = 30, sbins = 32;    
 int channels[] = {0,  1};   
 int histSize[] = {hbins, sbins};   
 float hranges[] = { 0, 180 };   
 float sranges[] = { 0, 255 };   
 const float* ranges[] = { hranges, sranges};    

 Mat patch_HSV;   
 MatND HistA, HistB;   

 //cal histogram & normalization   
 ///   
 cvtColor(imgA, patch_HSV, CV_BGR2HSV);   
 calcHist( &patch_HSV, 1, channels,  Mat(), // do not use mask   
  HistA, 2, histSize, ranges,   
  true, // the histogram is uniform   
  false );   
 normalize(HistA, HistA,  0, 1, CV_MINMAX);   


 cvtColor(imgB, patch_HSV, CV_BGR2HSV);   
 calcHist( &patch_HSV, 1, channels,  Mat(),// do not use mask   
  HistB, 2, histSize, ranges,   
  true, // the histogram is uniform   
  false );   
 normalize(HistB, HistB, 0, 1, CV_MINMAX);   

 //compare histogram   
 ///   
 int numrows = hbins * sbins;

 //make signature
 Mat sig1(numrows, 3, CV_32FC1);
 Mat sig2(numrows, 3, CV_32FC1);

 //fill value into signature
 for(int h=0; h< hbins; h++)
 {
  for(int s=0; s< sbins; ++s)
  {
   float binval = HistA.at< float>(h,s);
   sig1.at<float>( h*sbins + s, 0) = binval;
   sig1.at<float>( h*sbins + s, 1) = h;
   sig1.at<float>( h*sbins + s, 2) = s;

   binval = HistB.at< float>(h,s);
   sig2.at<float>( h*sbins + s, 0) = binval;
   sig2.at<float>( h*sbins + s, 1) = h;
   sig2.at<float>( h*sbins + s, 2) = s;
  }
 }

 //compare similarity of 2images using emd.
 float emd = cv::EMD(sig1, sig2, CV_DIST_L2); //emd 0 is best matching. 
 printf("similarity %5.5f %%\n", (1-emd)*100 );
 
 waitKey(0);   

 return 0;   
}