图像相似性度量—— 归一化互信息实现

//************************************************************************
// 函数名称:    	GetMutualInfo
// 访问权限:    	public 
// 创建日期:		2016/12/09
// 创 建 人:		
// 函数说明:		计算输入图像A和输入图像B的互信息
// 函数参数: 	cv::Mat & img_a	输入图像A
// 函数参数: 	cv::Mat & img_b	输入图像B
// 返 回 值:   	double
//************************************************************************
double CCalcMutualInfo::GetMutualInfo(cv::Mat& img_a, cv::Mat& img_b)
{
	if (!img_a.data || !img_b.data)
	{
		cout << "no input img" << endl;
		return 0.0;
	}
	if (img_a.rows!=img_a.rows || img_a.cols!=img_b.cols)
	{
		cout << "input img's row and cosl not eqaul" << endl;
		return 0.0;
	}

	int rows(img_a.rows);
	int cols(img_b.cols);
	int total_pixel(rows*cols);
	unsigned char* data_a = nullptr;
	unsigned char* data_b = nullptr;
	double value(0.0);						//计算得到的互信息的结果
	double H_a(0.0), H_b(0.0), H_ab(0.0);	//图像A,B或是AB的信息熵
	double* count_array_a = new double[256];
	memset(count_array_a, 0, sizeof(double)*256);
	double* count_array_b = new double[256];
	memset(count_array_b, 0, sizeof(double)* 256);
	double* count_array_ab = new double[256*256];
	memset(count_array_ab, 0, sizeof(double)* 256 * 256);

	//计算H_a, H_b
	for (int i=0; i<rows; i++)
	{
		data_a = img_a.ptr<unsigned char>(i);
		data_b = img_b.ptr<unsigned char>(i);
		for (int j=0; j<cols; j++)
		{
			count_array_a[data_a[j]]++;
			count_array_b[data_b[j]]++;
		}
	}
	for (int i=0; i<255; i++)
	{
		if (0.0 != count_array_a[i])
		{
			double p(count_array_a[i] / (double)total_pixel);
			H_a = H_a + (-1.0*p*(std::log(p) / std::log(2)));
		}
		
		if (0.0 != count_array_b[i])
		{
			double p(count_array_b[i] / (double)total_pixel);
			H_b = H_b + (-1.0*p*(std::log(p) / std::log(2)));
		}
	}

	//计算H_ab
	for (int m = 0; m < 256; m++)	//8位的灰度级
	{
		for (int n = 0; n < 256; n++)
		{
			for (int i = 0; i < rows; i++)
			{
				data_a = img_a.ptr<unsigned char>(i);
				data_b = img_b.ptr<unsigned char>(i);
				for (int j = 0; j < cols; j++)
				{
					//if ((std::abs(m-data_a[j])<20) && (std::abs(n-data_b[j])<20))
					if ((m == data_a[j]) && (n == data_b[j]))	//由联合概率密度的实际物理意义，统计点数
					{
						count_array_ab[m*256 + n]++;
					}
				}
			}
		}
	}
	for (int m = 0; m < 256; m++)
	{
		for (int n = 0; n < 256; n++)
		{
			if (0.0 != count_array_ab[m*256 + n])
			{
				double p(count_array_ab[m*256 + n] / (double)total_pixel);
				H_ab = H_ab + (-1.0*p*(std::log(p) / std::log(2)));
			}
		}
	}
	
	value = (H_a + H_b) / H_ab;	//得出归一化互信息
	delete[] count_array_a;
	delete[] count_array_b;
	delete[] count_array_ab;
	count_array_ab = nullptr;
	count_array_a = nullptr;
	count_array_b = nullptr;
	return value;
}