Image Smoothing via L0 Gradient Minimization基于C++的代码实现

论文解析图像梯度L0范数最小化图像梯度最小化平滑---一维信号图像梯度最小化平滑---二维图像源代码//实现L0测度平滑Mat L0Smoothing(Mat & image8UC3, double lambda = 2e-2, double kappa = 2.0){//将输入的图片转换为三通道的Double类型,按原来的像素值进行1...

铿锵的玫瑰

964人浏览 · 2018-11-11 22:08:44

铿锵的玫瑰 · 2018-11-11 22:08:44 发布

论文解析

图像梯度L0范数最小化

图像梯度最小化平滑---一维信号

图像梯度最小化平滑---二维图像

源代码

//实现L0测度平滑
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;

	定义sobel算子在x、y方向上进行卷积的参数
	int scale = 1;
	int delta = 0;
	int ddepth = CV_64F;
	Mat fx(3, 3, CV_64FC1);
	Mat fy(3, 3, CV_64FC1);
	//构造在x方向上的卷积核
	fx.at<double>(0, 0) = -1; fx.at<double>(0, 1) = -2; fx.at<double>(0, 2) = -1;
	fx.at<double>(1, 0) = 0; fx.at<double>(1, 1) = 0; fx.at<double>(1, 2) = 0;
	fx.at<double>(2, 0) = 1; fx.at<double>(2, 1) = 2; fx.at<double>(2, 2) = 1;
	//构造在y方向上的卷积核
	fy.at<double>(0, 0) = -1; fy.at<double>(0, 1) = 0; fy.at<double>(0, 2) = 1;
	fy.at<double>(1, 0) = -2; fy.at<double>(1, 1) = 0; fy.at<double>(1, 2) = 2;
	fy.at<double>(2, 0) = -1; fy.at<double>(2, 1) = 0; fy.at<double>(2, 2) = 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];

			利用Sobel算子进行求解
			//求解X方向上的梯度
			Mat shifted_x = single_channel[k].clone();
			Sobel(shifted_x, dx[k], ddepth, 1, 0, 3, scale, delta, BORDER_DEFAULT);
			//求解X方向上的梯度
			Mat shifted_y = single_channel[k].clone();
			Sobel(shifted_y, dy[k], ddepth, 0, 1, 3, scale, delta, BORDER_DEFAULT);
		}
		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];

			//利用Sobel算子求解二阶导数
			//求解X方向上的梯度
			Mat shifted_dx = dx[k].clone();
			Mat ddx;
			Sobel(shifted_dx, ddx, ddepth, 1, 0, 3, scale, delta, BORDER_DEFAULT);
			//求解X方向上的梯度
			Mat shifted_dy = dy[k].clone();
			Mat ddy;
			Sobel(shifted_dy, ddy, ddepth, 0, 1, 3, scale, delta, BORDER_DEFAULT);


			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;
	}

	创建显示图片的窗口
	//namedWindow("single_channel[0] Image", CV_WINDOW_AUTOSIZE);
	显示原始图片
	//imshow("single_channel[0] Image", single_channel[0]);
	创建显示图片的窗口
	//namedWindow("single_channel[1] Image", CV_WINDOW_AUTOSIZE);
	显示原始图片
	//imshow("single_channel[1] Image", single_channel[1]);
	创建显示图片的窗口
	//namedWindow("single_channel[2] Image", CV_WINDOW_AUTOSIZE);
	显示原始图片
	//imshow("single_channel[2] Image", single_channel[2]);

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