计算Affinity矩阵，根据能量函数计算output能量... (C++)

目的

全局匹配方法往往会给到一个能量函数，思想不管如何，目的都是希望能够使这个能量最小，本文
根据快速双边求解器（简称fbs）的思想（见下式），直接求解图像的能量函数，可有效帮助自己理解全局匹配方法，也可从特定的层面，检验fbs输出的正确与否。

本文提供了自己写的计算图像能量函数的C++代码，以做记录，可应用于三通道的reference图像,单通道的target和confidence图像。

1. vector_idx2row_col_idx ： 可以将1D的下标转换为2D的行列下标。

---

1. 参数说明
vector_idx： 黑框中的数字即为vector_idx, 即为2D矩阵一行一行取出来，归置为1D向量的下标。

• cols ： 以上图为例，cols为5，为2D向量的宽度。

	// 将横向的vector序标，转换为2D的行列号
	pair<int, int> vector_idx2row_col_idx(int vector_idx, int cols) {
		int row_idx = floor(vector_idx / cols);
		int col_idx = vector_idx - row_idx*cols;
		pair<int, int> p1(row_idx, col_idx);
		return p1;
	}

2. calSimilarity（三通道） ： 计算两点之间的相似度

---

• 计算公式：
  

/*
	calSimilarity
	--- 作用：根据公式计算两点之间的相似度
	exp(-(x1-x2)^2 / 2 * sigam_spatial^2) * ...
	@ yuv: 三通道的yuv图像
	@ pair_base: base点的行列号
	@ pair_search: search点的行列号
	@ sigma_spatial, sigma_luma, sigma_chroma: 控制空间域,色彩域的参数
	*/
	double calSimilarity(Mat& yuv, pair<int, int> pair_base, pair<int, int> pair_search, int sigma_spatial, int sigma_luma, int sigma_chroma) {
		int delta_x = pair_base.first - pair_search.first;
		double part_x = exp(-(delta_x * delta_x) * 0.5 / (sigma_spatial * sigma_spatial));
		int delta_y = pair_base.second - pair_search.second;
		double part_y = exp(-(delta_y * delta_y) * 0.5 / (sigma_spatial * sigma_spatial));

		// 分割yuv通道
		cv::Mat imageY(yuv.rows, yuv.cols, 1);
		cv::Mat imageU(yuv.rows, yuv.cols, 1);
		cv::Mat imageV(yuv.rows, yuv.cols, 1);
		std::vector<Mat> mv;
		split(yuv, (vector<Mat>&)mv);
		imageY = mv[0].clone();
		imageU = mv[1].clone();
		imageV = mv[2].clone();

		double delta_luma = imageY.at<uchar>(pair_base.first, pair_base.second) - imageY.at<uchar>(pair_search.first, pair_search.second);
		double part_luma = exp(-(delta_luma*delta_luma) * 0.5 / (sigma_luma * sigma_luma));
		double delta_chroma1 = imageU.at<uchar>(pair_base.first, pair_base.second) - imageU.at<uchar>(pair_search.first, pair_search.second);
		double part_chroma1 = exp(-(delta_chroma1*delta_chroma1) * 0.5 / (sigma_chroma * sigma_chroma));
		double delta_chroma2 = imageV.at<uchar>(pair_base.first, pair_base.second) - imageV.at<uchar>(pair_search.first, pair_search.second);
		double part_chroma2 = exp(-(delta_chroma2 * delta_chroma2) * 0.5 / (sigma_chroma * sigma_chroma));
		double res = part_x * part_y * part_luma * part_chroma1 * part_chroma2;

		return res;
	}

3. 计算Affinity Matrix

• 填充图像中每两个点之间的相似度（计算公式见标题2）

	void calAffinityMatrix(const Mat &reference, double sigma_spatial, double sigma_luma, double sigma_chroma, Mat &Affinity_matrix) {
		cv::Mat reference_yuv;
		cv::cvtColor(reference, reference_yuv, COLOR_BGR2YCrCb);
		int cols = reference_yuv.cols;
		int rows = reference_yuv.rows;
		//const unsigned char* pref = (const unsigned char*)reference_yuv.data;

		int npixels = rows*cols;

		for (int i = 0; i<npixels; ++i) {
			pair<int, int> base_ind = vector_idx2row_col_idx(i, cols);
			float* ptr_Affinity = Affinity_matrix.ptr<float>(i);
			for (int j = 0; j<npixels; ++j) {
				pair<int, int> search_ind = vector_idx2row_col_idx(j, cols);
				ptr_Affinity[j] = calSimilarity(reference_yuv, base_ind, search_ind, sigma_spatial, sigma_luma, sigma_chroma);
			}
		}
	}

4. 计算图像能量

---

1. 计算公式

double calEnergy(const Mat& Affinty_matrix, const Mat& output, const Mat& confidence, const Mat& target, double lambda) {

	double energy;
	double smooth_part = 0;
	double keepReal_part = 0;
	
	for (int i = 0;i < Affinty_matrix.rows;++i) {
		const float* ptr_Affinity = Affinty_matrix.ptr<float>(i);
		for (int j = 0; j < Affinty_matrix.cols;++j) {
			pair<int, int> base_idx = vector_idx2row_col_idx(i, output.cols);
			pair<int, int> search_idx = vector_idx2row_col_idx(j, output.cols);
			float delta_x = (output.ptr<float>(base_idx.first)[base_idx.second] - output.ptr<float>(search_idx.first)[search_idx.second]);
			double temp_smooth = ptr_Affinity[j] * delta_x * delta_x;
			smooth_part += temp_smooth;

			double ci = confidence.ptr<float>(base_idx.first)[base_idx.second];
			double delta_xi_minus_ti = output.ptr<float>(base_idx.first)[base_idx.second] - target.ptr<float>(base_idx.first)[base_idx.second];
			double temp_keepReal = ci * delta_xi_minus_ti * delta_xi_minus_ti;
			keepReal_part += temp_keepReal;
		}
	}
	smooth_part = 0.5 * lambda * smooth_part;
	energy = keepReal_part + smooth_part;
	return energy;
}