基于L0测度的图像平滑

视频网址：https://www.bilibili.com/video/av35811232

完整代码：https://download.youkuaiyun.com/download/lykymy/10780311

关键代码

Mat L0Smoothing(Mat & image8UC3, double lambda = 2e-2, double kappa = 2.0)
{
	//将输入的图片转换为三通道的Double类型,按原来的像素值进行1：255的缩放
	Mat image64D3;
	image8UC3.convertTo(image64D3, CV_64FC3, 1.0 / 255.0);

	//定义roberts算子在x，y方向上的卷积核
	Mat fx(1, 2, CV_64FC1);
	Mat fy(2, 1, CV_64FC1);
	fx.at<double>(0) = 1; fx.at<double>(1) = -1;
	fy.at<double>(0) = 1; fy.at<double>(1) = -1;


	//把一个空间点扩散函数转换为频谱面的光学传递函数,即对图像进行FFT变换
	Mat otfFx = psf2otf(fx, image8UC3.size());
	Mat otfFy = psf2otf(fy, image8UC3.size());
	//定义DFT以后的输出图像，FNormal为两个通道
	Mat FNormin1[3];
	Mat single_channel[3];
	//分离为三个通道
	split(image64D3, single_channel);
	for (int k = 0; k < 3; k++) {
		// 离散傅里叶变换，指定输出复数格式（默认是CCS格式）
		dft(single_channel[k], FNormin1[k], DFT_COMPLEX_OUTPUT);
	}
	
	//F(∂x)∗F(∂x)+F(∂y)∗F(∂y);
	Mat Denormin2(image8UC3.rows, image8UC3.cols, CV_64FC1);
	for (int i = 0; i < image8UC3.rows; i++) {
		for (int j = 0; j < image8UC3.cols; j++) {
			Vec2f &c1 = otfFx.at<Vec2f>(i, j);
			Vec2f &c2 = otfFy.at<Vec2f>(i, j);
			// 0: Real, 1: Image
			Denormin2.at<double>(i, j) = pow(c1[0], 2) + pow(c1[1], 2) + pow(c2[0], 2) + pow(c2[1], 2);
		}
	}
	double beta = 2.0 * lambda;
	double betamax = 1e5;

	while (beta < betamax) {
		// F(1)+β(F(∂x)∗F(∂x)+F(∂y)∗F(∂y))
		Mat Denormin = 1.0 + beta * Denormin2;

		// h-v subproblem
		// 三个通道的 ∂S/∂x ∂S/∂y
		Mat dx[3], dy[3];
		for (int k = 0; k < 3; k++) {
			//利用Roberts算子进行求解
			//求解X方向上的梯度
			Mat shifted_x = single_channel[k].clone();
			circshift(shifted_x, 0, -1);
			dx[k] = shifted_x - single_channel[k];
			//求解Y方向上的梯度
			Mat shifted_y = single_channel[k].clone();
			circshift(shifted_y, -1, 0);
			dy[k] = shifted_y - single_channel[k];

		}
		for (int i = 0; i < image8UC3.rows; i++) {
			for (int j = 0; j < image8UC3.cols; j++) {
				// (∂S/∂x)^2 + (∂S/∂y)^2
				double val = pow(dx[0].at<double>(i, j), 2) + pow(dy[0].at<double>(i, j), 2) +
					         pow(dx[1].at<double>(i, j), 2) + pow(dy[1].at<double>(i, j), 2) +
					         pow(dx[2].at<double>(i, j), 2) + pow(dy[2].at<double>(i, j), 2);
				// (∂S/∂x)^2 + (∂S/∂y)^2  < λ/β
				if (val < lambda / beta) {
					dx[0].at<double>(i, j) = dx[1].at<double>(i, j) = dx[2].at<double>(i, j) = 0.0;
					dy[0].at<double>(i, j) = dy[1].at<double>(i, j) = dy[2].at<double>(i, j) = 0.0;
				}
			}
		}

		// S subproblem
		for (int k = 0; k < 3; k++) {
			//利用Roberts算子求解二阶导数
			//求解X方向上的梯度
			Mat shift_dx = dx[k].clone();
			circshift(shift_dx, 0, 1);
			Mat ddx = shift_dx - dx[k];
			//求解Y方向上的梯度
			Mat shift_dy = dy[k].clone();
			circshift(shift_dy, 1, 0);
			Mat ddy = shift_dy - dy[k];

			Mat Normin2 = ddx + ddy;
			Mat FNormin2;
			// 离散傅里叶变换，指定输出复数格式（默认是CCS格式）
			dft(Normin2, FNormin2, cv::DFT_COMPLEX_OUTPUT);
			// F(g)+β(F(∂x)∗F(h)+F(∂y)∗F(v))
			Mat FS = FNormin1[k] + beta*FNormin2;

			// 论文的公式(8)：F^-1括号中的内容
			for (int i = 0; i < image8UC3.rows; i++) {
				for (int j = 0; j < image8UC3.cols; j++) {
					FS.at<Vec2d>(i, j)[0] /= Denormin.at<double>(i, j);
					FS.at<Vec2d>(i, j)[1] /= Denormin.at<double>(i, j);
				}
			}
			// 论文的公式(8)：傅里叶逆变换
			Mat ifft;
			idft(FS, ifft, cv::DFT_SCALE | cv::DFT_COMPLEX_OUTPUT);
			for (int i = 0; i < image8UC3.rows; i++) {
				for (int j = 0; j < image8UC3.cols; j++) {
					single_channel[k].at<double>(i, j) = ifft.at<cv::Vec2d>(i, j)[0];
				}
			}
		}
		beta *= kappa;
	}

	//合并三个通道
	merge(single_channel, 3, image64D3);
	//创建显示图片的窗口
	namedWindow("L0 Smoothing Image", CV_WINDOW_AUTOSIZE);
	//显示原始图片
	imshow("L0 Smoothing Image", image64D3);
	return image64D3;
}

效果

原始图片

L0平滑后的图片