非局部均值滤波(磨皮)

note

根据掩膜区域与其他区域的相似程度来赋予掩膜区域像素权重

相似程度由均方差来衡量

code

/*
 \brief 矩阵求邻和
 \param type=1，列方向；type=2，行方向
*/
static void MyCumSum(Mat& src, Mat& res, int type) {
	if ((src.channels() > 1) || (res.channels() > 1)) {
		return;
	}
	src.copyTo(res);
	res.convertTo(res, CV_64FC1);
	if (type == 1) {
		for (int i = 0; i < src.cols; ++i) {
			for (int j = 1; j < src.rows; ++j) {
				res.at<double>(j,i) = res.at<double>(j-1,i) + res.at<double>(j,i);
			}
		}
	}
	else if (type == 2) {
		for (int i = 0; i < src.rows; ++i) {
			for (int j = 1; j < src.cols; ++j) {
				res.at<double>(i,j) = res.at<double>(i,j-1) + res.at<double>(i,j);
			}
		}
	}
}

void GetWeight(Mat& distance, Mat& weight, double h) {
	for (int i = 0; i < weight.rows; ++i) {
		for (int j = 0; j < weight.cols; ++j) {
			double w = distance.at<double>(i,j);
			weight.at<double>(i,j) = exp(-1*w/h/h);
		}
	}
}

void MyNonLocalMeanFilter(Mat& src, Mat& res, int ds = 2, int Ds = 5, double h = 10) {
	int srcType = src.type();
	if (srcType != CV_64FC1) {
		src.convertTo(src, CV_64FC1);
		res.convertTo(res, CV_64FC1);
	}

	int m = src.rows;
	int n = src.cols;
	int offset = ds+Ds;

	Mat paddedImg;
	copyMakeBorder(src, paddedImg, offset, offset, offset, offset, BORDER_REFLECT);	// 上下左右镜像复制矩阵边界

	Mat sumImage(m, n, CV_64FC1, Scalar(0));
	Mat sumWeight(m, n, CV_64FC1, Scalar(0));
	Mat maxWeight(m, n, CV_64FC1, Scalar(0));

	Rect rc;
	rc.x = Ds;
	rc.y = Ds;
	rc.width = m+ds;
	rc.height = n+ds;
	Mat image;	// 和搜索图比较使用的原图
	paddedImg(rc).copyTo(image);

	int M = image.rows;
	int N = image.cols;

	// 移动搜索框
	for (int r = -Ds; r <= Ds; r++) {
		for (int s = -Ds; s <= Ds; s++) {
			if (r == 0 && s == 0) {
				continue;
			}

			// 求差值积分图
			Rect recTmp;
			recTmp.x = Ds+s;
			recTmp.y = Ds+r;
			recTmp.width = n+ds;
			recTmp.height = m+ds;
			Mat wimage;	// 某个搜索框对应的搜索图
			paddedImg(recTmp).copyTo(wimage);

			Mat diff;	// 差值平方图
			subtract(image, wimage, diff);
			multiply(diff, diff, diff);

			Mat J;
			MyCumSum(diff, J, 1);	// 列方向求邻和
			MyCumSum(J, J, 2);	// 行方向求邻和

			recTmp.x = N-n;
			recTmp.y = M-m;
			recTmp.width = n;
			recTmp.height = m;
			Mat matTmp1;
			J(recTmp).copyTo(matTmp1);

			recTmp.x = 0;
			recTmp.y = 0;
			recTmp.width = n;
			recTmp.height = m;
			Mat matTmp2;
			J(recTmp).copyTo(matTmp2);

			recTmp.x = 0;
			recTmp.y = M-m;
			recTmp.width = n;
			recTmp.height = m;
			Mat matTmp3;
			J(recTmp).copyTo(matTmp3);

			recTmp.x = N-n;
			recTmp.y = 0;
			recTmp.width = n;
			recTmp.height = m;
			Mat matTmp4;
			J(recTmp).copyTo(matTmp4);

			// 计算距离
			Mat distance;
			add(matTmp1, matTmp2, distance);
			subtract(distance, matTmp3, distance);
			subtract(distance, matTmp4, distance);
			double var = (2*ds+1) * (2*ds+1);
			distance = distance / var;

			// 计算权重并获得单个偏移下的加权图像
			Mat weight(distance.rows, distance.cols, CV_64FC1, Scalar(0));
			GetWeight(distance, weight, h);

			recTmp.x = ds;
			recTmp.y = ds;
			recTmp.width = n;
			recTmp.height = m;
			Mat matTmp5;
			wimage(recTmp).copyTo(matTmp5);

			multiply(weight, matTmp5, matTmp5);
			add(sumImage, matTmp5, sumImage);

			add(sumWeight, weight, sumWeight);

			maxWeight = cv::max(maxWeight, weight);
		}
	}
	
	rc.x = ds;
	rc.y = ds;
	rc.width = n;
	rc.height = m;
	Mat matTmp;
	image(rc).copyTo(matTmp);
	multiply(maxWeight, matTmp, matTmp);
	add(sumImage, matTmp, sumImage);
	add(sumWeight, maxWeight, sumWeight);
	divide(sumImage, sumWeight, res);

	src.convertTo(src, srcType);
	res.convertTo(res, srcType);
}

test