Opencv两图像融合

class WaveTransform
{
public:
	Mat WDT(const Mat &_src, const string _wname, const int _level);//小波分解
	Mat IWDT(const Mat &_src, const string _wname, const int _level);//小波重构
	void wavelet_D(const string _wname, Mat &_lowFilter, Mat &_highFilter);//分解包
	void wavelet_R(const string _wname, Mat &_lowFilter, Mat &_highFilter);//重构包
	Mat waveletDecompose(const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter);
	Mat waveletReconstruct(const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter);
};
Mat WaveTransform::WDT(const Mat &_src, const string _wname, const int _level)
{
	Mat_ src = Mat_(_src);
	Mat dst = Mat::zeros(src.rows, src.cols, src.type());
	int row = src.rows;
	int col = src.cols;
	//高通低通滤波器
	Mat lowFilter;
	Mat highFilter;
	wavelet_D(_wname, lowFilter, highFilter);
	//小波变换
	int t = 1;

	while (t <= _level)
	{
		//先进行 行小波变换
		//#pragma omp parallel for
		for (int i = 0; i < row; i++)
		{
			//取出src中要处理的数据的一行
			Mat oneRow = Mat::zeros(1, col, src.type());

			for (int j = 0; j < col; j++)
			{
				oneRow.at(0, j) = src.at(i, j);
			}

			oneRow = waveletDecompose(oneRow, lowFilter, highFilter);
			for (int j = 0; j < col; j++)
			{
				dst.at(i, j) = oneRow.at(0, j);
			}
		}

		//小波列变换
		//#pragma omp parallel for
		for (int j = 0; j < col; j++)
		{
			Mat oneCol = Mat::zeros(row, 1, src.type());

			for (int i = 0; i < row; i++)
			{
				oneCol.at(i, 0) = dst.at(i, j);//dst,not src
			}
			oneCol = (waveletDecompose(oneCol.t(), lowFilter, highFilter)).t();

			for (int i = 0; i < row; i++)
			{
				dst.at(i, j) = oneCol.at(i, 0);
			}
		}
		//更新 
		row /= 2;
		col /= 2;
		t++;
		src = dst;
	}
	return dst;
}
Mat WaveTransform::IWDT(const Mat &_src, const string _wname, const int _level)
{

	Mat src = Mat_(_src);
	Mat dst;
	src.copyTo(dst);
	int N = src.rows;
	int D = src.cols;

	//高低通滤波器
	Mat lowFilter;
	Mat highFilter;
	wavelet_R(_wname, lowFilter, highFilter);

	//小波变换
	int t = 1;
	int row = N / std::pow(2., _level - 1);
	int col = D / std::pow(2., _level - 1);

	while (row <= N && col <= D)
		//while(t<=_level)
	{
		//列逆变换
		for (int j = 0; j < col; j++)
		{
			Mat oneCol = Mat::zeros(row, 1, src.type());

			for (int i = 0; i < row; i++)
			{
				oneCol.at(i, 0) = src.at(i, j);
			}
			oneCol = (waveletReconstruct(oneCol.t(), lowFilter, highFilter)).t();

			for (int i = 0; i < row; i++)
			{
				dst.at(i, j) = oneCol.at(i, 0);
			}

		}

		//行逆变换
		for (int i = 0; i < row; i++)
		{
			Mat oneRow = Mat::zeros(1, col, src.type());
			for (int j = 0; j < col; j++)
			{
				oneRow.at(0, j) = dst.at(i, j);
			}
			oneRow = waveletReconstruct(oneRow, lowFilter, highFilter);
			for (int j = 0; j < col; j++)
			{
				dst.at(i, j) = oneRow.at(0, j);
			}
		}

		row *= 2;
		col *= 2;
		t++;
		src = dst;
	}

	return dst;
}
void WaveTransform::wavelet_D(const string _wname, Mat &_lowFilter, Mat &_highFilter)
{
	if (_wname == "haar" || _wname == "db1")
	{
		int N = 2;
		_lowFilter = Mat::zeros(1, N, CV_32F);
		_highFilter = Mat::zeros(1, N, CV_32F);

		_lowFilter.at(0, 0) = 1 / sqrtf(N);
		_lowFilter.at(0, 1) = 1 / sqrtf(N);

		_highFilter.at(0, 0) = -1 / sqrtf(N);
		_highFilter.at(0, 1) = 1 / sqrtf(N);
	}
	else if (_wname == "sym2")
	{
		int N = 4;
		float h[] = { -0.4830, 0.8365, -0.2241, -0.1294 };
		float l[] = { -0.1294, 0.2241,  0.8365, 0.4830 };

		_lowFilter = Mat::zeros(1, N, CV_32F);
		_highFilter = Mat::zeros(1, N, CV_32F);

		for (int i = 0; i < N; i++)
		{
			_lowFilter.at(0, i) = l[i];
			_highFilter.at(0, i) = h[i];
		}
	}
}
void WaveTransform::wavelet_R(const string _wname, Mat &_lowFilter, Mat &_highFilter)
{
	if (_wname == "haar" || _wname == "db1")
	{
		int N = 2;
		_lowFilter = Mat::zeros(1, N, CV_32F);
		_highFilter = Mat::zeros(1, N, CV_32F);


		_lowFilter.at(0, 0) = 1 / sqrtf(N);
		_lowFilter.at(0, 1) = 1 / sqrtf(N);

		_highFilter.at(0, 0) = 1 / sqrtf(N);
		_highFilter.at(0, 1) = -1 / sqrtf(N);
	}
	else if (_wname == "sym2")
	{
		int N = 4;
		float h[] = { -0.1294,-0.2241,0.8365,-0.4830 };
		float l[] = { 0.4830, 0.8365, 0.2241, -0.1294 };

		_lowFilter = Mat::zeros(1, N, CV_32F);
		_highFilter = Mat::zeros(1, N, CV_32F);

		for (int i = 0; i < N; i++)
		{
			_lowFilter.at(0, i) = l[i];
			_highFilter.at(0, i) = h[i];
		}
	}
}
Mat WaveTransform::waveletDecompose(const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter)
{
	assert(_src.rows == 1 && _lowFilter.rows == 1 && _highFilter.rows == 1);
	assert(_src.cols >= _lowFilter.cols && _src.cols >= _highFilter.cols);
	Mat &src = Mat_(_src);

	int D = src.cols;

	Mat &lowFilter = Mat_(_lowFilter);
	Mat &highFilter = Mat_(_highFilter);

	//频域滤波或时域卷积；ifft( fft(x) * fft(filter)) = cov(x,filter) 
	Mat dst1 = Mat::zeros(1, D, src.type());
	Mat dst2 = Mat::zeros(1, D, src.type());

	filter2D(src, dst1, -1, lowFilter);
	filter2D(src, dst2, -1, highFilter);

	//下采样
	//数据拼接
	for (int i = 0, j = 1; i < D / 2; i++, j += 2)
	{
		src.at(0, i) = dst1.at(0, j);//lowFilter
		src.at(0, i + D / 2) = dst2.at(0, j);//highFilter
	}
	return src;
}
Mat WaveTransform::waveletReconstruct(const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter)
{
	assert(_src.rows == 1 && _lowFilter.rows == 1 && _highFilter.rows == 1);
	assert(_src.cols >= _lowFilter.cols && _src.cols >= _highFilter.cols);
	Mat &src = Mat_(_src);

	int D = src.cols;

	Mat &lowFilter = Mat_(_lowFilter);
	Mat &highFilter = Mat_(_highFilter);


	/// 插值;
	Mat Up1 = Mat::zeros(1, D, src.type());
	Mat Up2 = Mat::zeros(1, D, src.type());


	for (int i = 0, cnt = 0; i < D / 2; i++, cnt += 2)
	{
		Up1.at(0, cnt) = src.at(0, i);     ///< 前一半
		Up2.at(0, cnt) = src.at(0, i + D / 2); ///< 后一半
	}

	/// 前一半低通，后一半高通
	Mat dst1 = Mat::zeros(1, D, src.type());
	Mat dst2 = Mat::zeros(1, D, src.type());
	filter2D(Up1, dst1, -1, lowFilter);
	filter2D(Up2, dst2, -1, highFilter);

	/// 结果相加
	dst1 = dst1 + dst2;
	return dst1;
}
void RGB2HSI(Mat src, Mat &dst) {
	Mat HSI(src.rows, src.cols, CV_32FC3);
	float r, g, b, H, S, I, num, den, theta, sum, min_RGB;
	for (int i = 0; i < src.rows; i++)
	{
		for (int j = 0; j < src.cols; j++)
		{
			b = src.at(i, j)[0];
			g = src.at(i, j)[1];
			r = src.at(i, j)[2];

			num = 0.5 * ((r - g) + (r - b));
			den = sqrt((r - g)*(r - g) + (r - b)*(g - b));

			if (den == 0) {
				H = 0; // 分母不能为0
			}
			else {
				theta = acos(num / den);
				if (b <= g) {
					H = theta;
				}
				else {
					H = (2 * PI - theta);
				}
			}

			min_RGB = min(min(b, g), r); // min(R,G,B)
			sum = b + g + r;
			if (sum == 0)
			{
				S = 0;
			}
			else {
				S = 1 - 3 * min_RGB / sum;
			}

			I = sum / 3.0;

			HSI.at(i, j)[0] = H;
			HSI.at(i, j)[1] = S;
			HSI.at(i, j)[2] = I;
		}
	}
	dst = HSI;
	return;
}
void HSI2RGB(Mat src, Mat &dst) {
	Mat RGB(src.size(), CV_32FC3);
	for (int i = 0; i < src.rows; i++) {
		for (int j = 0; j < src.cols; j++) {
			float DH = src.at(i, j)[0];
			float DS = src.at(i, j)[1];
			float DI = src.at(i, j)[2];
			//分扇区显示
			float R, G, B;
			if (DH < (2 * PI / 3) && DH >= 0) {
				B = DI * (1 - DS);
				R = DI * (1 + (DS * cos(DH)) / cos(PI / 3 - DH));
				G = (3 * DI - (R + B));
			}
			else if (DH < (4 * PI / 3) && DH >= (2 * PI / 3)) {
				DH = DH - (2 * PI / 3);
				R = DI * (1 - DS);
				G = DI * (1 + (DS * cos(DH)) / cos(PI / 3 - DH));
				B = (3 * DI - (G + R));
			}
			else {
				DH = DH - (4 * PI / 3);
				G = DI * (1 - DS);
				B = DI * (1 + (DS * cos(DH)) / cos(PI / 3 - DH));
				R = (3 * DI - (G + B));
			}
			RGB.at(i, j)[0] = B;
			RGB.at(i, j)[1] = G;
			RGB.at(i, j)[2] = R;
		}
	}
	dst = RGB;
	return;
}
void harr_fusion(Mat src1, Mat src2, Mat &dst) {
	assert(src1.rows == src2.rows&&src1.cols == src2.cols);
	int row = src1.rows;
	int col = src1.cols;
	Mat src1_gray, src2_gray;
	normalize(src1, src1_gray, 0, 255, CV_MINMAX);
	cvtColor(src2, src2_gray, CV_RGB2GRAY);
	normalize(src2_gray, src2_gray, 0, 255, CV_MINMAX);
	src1_gray.convertTo(src1_gray, CV_32F);
	src2_gray.convertTo(src2_gray, CV_32F);
	WaveTransform m_waveTransform;
	const int level = 2;
	Mat src1_dec = m_waveTransform.WDT(src1_gray, "haar", level);
	Mat src2_dec = m_waveTransform.WDT(src2_gray, "haar", level);
	Mat dec = Mat(row, col, CV_32FC1);
	//融合规则：高频部分采用加权平均的方法，低频部分采用模值取大的方法
	int halfRow = row / (2 * level);
	int halfCol = col / (2 * level);
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++) {
			if (i < halfRow&&j < halfCol) {
				dec.at(i, j) = (src1_dec.at(i, j) + src2_dec.at(i, j)) / 2;
			}
			else {
				float p = abs(src1_dec.at(i, j));
				float q = abs(src2_dec.at(i, j));
				if (p > q) {
					dec.at(i, j) = src1_dec.at(i, j);
				}
				else {
					dec.at(i, j) = src2_dec.at(i, j);
				}

			}
		}
	}
	dst = m_waveTransform.IWDT(dec, "haar", level);
}
void fusion(Mat Visible, Mat Infrared, Mat &dst) {
	Mat HSI(Visible.size(), CV_32FC3);
	Mat Visible_I(Visible.size(), CV_32FC1);
	RGB2HSI(Visible, HSI);
	for (int i = 0; i < Visible.rows; i++) {
		for (int j = 0; j < Visible.cols; j++) {
			Visible_I.at(i, j) = HSI.at(i, j)[2];
		}
	}
	Mat fusion_I;
	harr_fusion(Visible_I, Infrared, fusion_I);
	Mat fusion_dst = Mat::zeros(Visible.size(), CV_32FC3);
	for (int i = 0; i < Visible.rows; i++) {
		for (int j = 0; j < Visible.cols; j++) {
			fusion_dst.at(i, j)[2] = fusion_I.at(i, j);
			fusion_dst.at(i, j)[0] = HSI.at(i, j)[0];
			fusion_dst.at(i, j)[1] = HSI.at(i, j)[1];
		}
	}
	HSI2RGB(fusion_dst, dst);
	dst.convertTo(dst, CV_8UC3);
	normalize(dst, dst, 0, 255, CV_MINMAX);
}

int _test_0605()
{
	Mat Visible = imread("C:\Users\HP\Desktop\test\0605\1.bmp");    //右图
	Mat Infrared = imread("C:\Users\HP\Desktop\test\0605\3.bmp");    //左图
	imshow("左上", Visible);
	imshow("左下", Infrared);

	Mat dst;
	fusion(Visible, Infrared, dst);
	Mat dst1;
	imshow("fusion", dst);
	normalize(part4_r, part4_r, 0, 1, NORM_MINMAX);//归一化处理
	normalize(part4_r, part4_r, 0, 255, NORM_MINMAX);
	imshow("fusion1", part4_r);

	imwrite("C:\Users\HP\Desktop\test\0605\test3.bmp", part4_r);
	waitKey();
	destroyAllWindows();
	return 0;
}








int main()
{
	_test_0605();
	return 0;
}