基于陆地移动距离（EMD）的彩色图像直方图距离计算

直方图匹配策略

直方图是用来比较两幅图像相似度的重要途径，常用的距离标准包括：相关、卡方、直方图相交以及Bhattacharyya距离。

光线变化能使得图像颜色值的飘逸，尽管颜色直方图形状没有改变（整体位移），这样的变化使得以上匹配策略失效【我的理解为，此时匹配仍可以得到正确结论，即具有越相似直方图统计特性，会在四种评判中得到更好的结果，但是可能已经丧失其原本的意义了，数值大小不能反应图像之间的差异，即不再具有判断图像相似度的能力】。

举个例子图像A,B,C，我们想从B和C中找到与A更相似的图片（视觉上），当这两者根据某一标准的结果都很差时，即使相对较好的结果我们也不能说它和图像A更相似（虽然在统计上是这样的）。所以需要提出一种比较靠谱的度量方法。

陆地移动距离

更靠谱的距离度量方法就是陆地移动距离（EMD），下面是The Earth Mover’s Distance as a Metric for Image Retrieval（EMD）这篇文献摘要的翻译，看了可能对EMD有更好的认识:

我们研究了一种基于图像内容检索方法的直方图分布特性，陆地移动距离（EMD）。 EMD的意义是将一种直方图分布转化为另一种分布所要付出的最小成本。对于图像检索来说，我们将EMD和矢量量化结合一起，该方法相比以往方法在视觉相似性表现出了更好的方法。EMD相比其他直方图匹配算法，具有更佳的鲁棒性。当用来比较相同质量分布时，EMD度量更加真实。在本文中，我们专注于应用于颜色和纹理，比较EMD了的检索性能与其他距离的检索性能。

因此，EMD距离模型应该还是挺有用途的。

PS：

OPENCV中cvCalcEMD2（）函数能计算的最大直方图签名的行数为15440。

（1）对于1维直方图，区间数目不能超过15440；

（2）对于2维直方图，2个区间数目的乘积不能超过15440，区间相同的话每个区间不能超过124；

（3）对于3维直方图，3个区间数目的乘积不能超过15440，区间相同的话每个区间不能超过24.9；

（4）实际使用的时候，因为同一个程序中的其他部分还要占用内存，所以会比计算的值更少。

具体参见：http://www.cnblogs.com/xrwang/archive/2010/02/02/cvCalcEMD2Exception.html

实验：card1、card2分别是光照较强时的图像，card2相对card1有旋转，card3是光照较弱时的图像（由于我是手机拍的，可能不太清楚）

计算了相互间直方图相交以及EMD度量。

#include<cv.h>
#include<highgui.h>

//输入bgr彩色图像，进行分割，返回计算的直方图
//由于彩色图像，维数默认3，size为bin尺寸，range为bin范围
//如果 norm =1，返回归一化的直方图
CvHistogram* cvCalcbgrHist(IplImage *src,int* size,float** ranges,int norm=1)
{
if (src->nChannels==3)
{
//分割单通道
IplImage *b_plane = cvCreateImage(cvGetSize(src),8,1);
IplImage *g_plane = cvCreateImage(cvGetSize(src),8,1);
IplImage *r_plane = cvCreateImage(cvGetSize(src),8,1);
IplImage *plane[] = {b_plane,g_plane,r_plane};
cvCvtPixToPlane(src,b_plane,g_plane,r_plane,0);

CvHistogram *hist;
hist = cvCreateHist(3,size,CV_HIST_ARRAY,ranges,1);
cvCalcHist(plane,hist,0,0);
if(norm ==1)
cvNormalizeHist(hist,1.0);
return hist;
}
else
{
printf("it must be color image");
}
}

//计算输入两个彩色图像直方图的陆地移动距离（EMD）
//输入两个直方图，bins_size相同
float cvCalbgrEMD(CvHistogram* hist1,CvHistogram* hist2,int *size)
{
//首先创建signature，用于cvCalcEMD2的输入
CvMat *sig1;
CvMat *sig2;
int b=0,g=0,r=0;
int b_bins = size[0];
int g_bins = size[1];
int r_bins = size[2];
int numrows = b_bins*g_bins*r_bins;

sig1 = cvCreateMat(numrows,4,CV_32FC1);  //1count + 3 coords = 4
sig2 = cvCreateMat(numrows,4,CV_32FC1);

for( b=0;b<b_bins;b++)
{
for( g=0;g<g_bins;g++)
{
for( r=0;r<r_bins;r++)
{
float bin_val = cvQueryHistValue_3D(hist1,b,g,r);
int pos = b*g_bins*r_bins+g*r_bins+r;
cvSet2D(sig1,pos,0,cvScalar(bin_val));
cvSet2D(sig1,pos,1,cvScalar(b));
cvSet2D(sig1,pos,2,cvScalar(g));
cvSet2D(sig1,pos,3,cvScalar(r));
bin_val = cvQueryHistValue_3D(hist2,b,g,r);
cvSet2D(sig2,pos,0,cvScalar(bin_val));
cvSet2D(sig2,pos,1,cvScalar(b));
cvSet2D(sig2,pos,2,cvScalar(g));
cvSet2D(sig2,pos,3,cvScalar(r));
}
}
}
float emd = cvCalcEMD2(sig1,sig2,CV_DIST_L2);
return emd;
}

int main()
{
IplImage *card1;
IplImage *card2;
IplImage *card3;
card1 = cvLoadImage("card1.jpg");
card2 = cvLoadImage("card2.jpg");
card3 = cvLoadImage("card3.jpg");

//定义基本参数
int b_bins = 20, g_bins =20, r_bins =20;
int hist_size[]  = {b_bins,g_bins,r_bins};
float b_range[] = {0,255};
float g_range[] = {0,255};
float r_range[] = {0,255};
float *range[] = {b_range,g_range,r_range};

//计算直方图
CvHistogram *hist_card1=cvCalcbgrHist(card1,hist_size,range,1);
CvHistogram *hist_card2=cvCalcbgrHist(card2,hist_size,range,1);
CvHistogram *hist_card3=cvCalcbgrHist(card3,hist_size,range,1);

//直方图相交法INTERSECT
double value1 = cvCompareHist(hist_card1,hist_card2,CV_COMP_INTERSECT);
double value2 = cvCompareHist(hist_card1,hist_card3,CV_COMP_INTERSECT);

//EMD
double value3 = cvCalbgrEMD(hist_card1,hist_card2,hist_size);
double value4 = cvCalbgrEMD(hist_card1,hist_card3,hist_size);

cvNamedWindow("card1",1);
cvShowImage("card1",card1);
cvNamedWindow("card2",1);
cvShowImage("card2",card2);
cvNamedWindow("card3",1);
cvShowImage("card3",card3);

printf("INTERSECT Distance : card1&card2 %f\n",value1);
printf("INTERSECT Distance : card1&card3 %f\n",value2);
printf("EMD Distance : card1&card2 %f\n",value3);
printf("EMD Distance : card1&card3 %f\n",value4);

cvWaitKey(0);

return 0;
}