再探基于L0测度的流场平滑

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

代码网址：https://download.csdn.net/download/lykymy/10833542

实现思想

本次实验主要是根据《Denoising point sets via L0 minimization》以及《Image Smoothing via L0 Graient Minimization》这两篇论文实现的。在这两篇论文中它们都涉及到同一个思想，那就是L0范式的优化问题。其实际上就是借助辅助变量，将L0优化问题转化为L2优化问题，从而较为方便地进行求解。在这，我们可以发现，针对不同的处理对象，其设计的函数也不相同。下面将具体进行介绍。

关于L0以及其他的相关范式的定义，可以参考我写的相关博客。

L0范式实际上就是求解一个向量中非0向量的个数，其具体的表达式如下所示：

其中符号|*|0就是代表求解非0的个数。

我们的目的便是为了实现L0能量公式的最小化，即如下所示：

由于该公式很难直接求解出最优问题，因此需要引入辅助变量，将其转化为基于L2的范式进行最优解的求取。因此，上述公式可以转化为下列公式：

下面就需要根据具体进行具体分析：

针对一维的信号，S是表示输入图像，D(S)是表示向量中非0的个数

针对二维的图像，S是表示图像中每个像素的颜色，D(S)是表示颜色的梯度

针对三维的图像，S是表示网格中点的位置，D(S)是表示基于Laplcaian算子的边缘的导数

为了更好的求解上述公式，我们需要将其分解为两个子问题：一个是求解S，另一个是求解辅助变量。

在求解辅助变量的时候，我们需要将S看成是一个常量，利用下述公式进行求解

通过归纳证明，可以得出

具体的证明过程可以参考这篇论文《Image Smoothing via L0 Graient Minimization》。

在求解S的时候，我们需要将辅助变量看成是一个常量，利用下述公式进行求解

我们将持续迭代上述公式，从而求解最优化问题。

通过上述的讲解，我们发现：只要根据不同的情况，定义好不同的S以及D(S)便可以较为方便地将L0范式转化为L2范式，从而快速地进行求解。

这里我们定义：

S为图像的流场

D(S)为任意像素与周围上、下、左、右四个像素的欧式距离

具体的实现算法如下所示：

流场的绘制

正0度与正90度的图像

正25度与负25度的图像

正45度与负45度的图像

正75度与负75度的图像

部分代码

//实现l0平滑图像
Mat l0_smooth(Mat & image, int nu, int mu, int beta_start, int beta_end)
{
    char energy_name[50] = "l0 smooth energy";
	//创建grad_x和grad_y矩阵
	Mat grad_x = Mat::zeros(image.size(), CV_16S);
	Mat grad_y = Mat::zeros(image.size(), CV_16S);
	//创建roberts算子处理显示图片
	Mat roberts_img = Mat::zeros(image.size(), CV_16S);
	for (int i = 0; i < image.rows - 1; i++)
	{
		for (int j = 0; j < image.cols - 1; j++)
		{
			//求取roberts算子水平方向的梯度
			grad_x.at<short>(i, j) = image.at<uchar>(i, j) - image.at<uchar>(i + 1, j + 1);
			//求取roberts算子垂直方向的梯度
			grad_y.at<short>(i, j) = image.at<uchar>(i, j + 1) - image.at<uchar>(i + 1, j);
		}
	}
	//构建流场
	for (int i = 0; i < image.rows; i++)
	{
		for (int j = 0; j < image.cols; j++)
		{
			if (grad_x.at<short>(i, j) != 0)
			{
				roberts_img.at<short>(i, j) = atan(grad_y.at<short>(i, j) / grad_x.at<short>(i, j)) * 180 / 3.1415;
			}
			else
			{
				roberts_img.at<short>(i, j) = 90;
			}
		}
	}
	//Mat zeta_value = sobel_img.clone();
	//long temp = 0, up_difference = 0, down_difference = 0, left_difference = 0, right_difference = 0;
	int beta_current = beta_start;
	int time = 0;
	MatrixXf second_derivative = MatrixXf::Identity(8, 8);
	second_derivative  = second_derivative * 2;
	second_derivative = second_derivative.inverse();
	while (beta_current <= beta_end)
	{
		//计算的D(s)
		/*for (int i = 1; i < image.rows - 1; i++)
		{
			for (int j = 1; j < image.cols - 1; j++)
			{
				up_difference = pow((image.at<float>(i, j) - image.at<float>(i - 1, j)), 2);
				down_difference = pow((image.at<float>(i, j) - image.at<float>(i + 1, j)), 2);
				left_difference = pow((image.at<float>(i, j) - image.at<float>(i, j - 1)), 2);
				right_difference = pow((image.at<float>(i, j) - image.at<float>(i, j + 1)), 2);
				temp = up_difference + down_difference + left_difference + right_difference;
				if (temp < nu / beta_current)
				{
					zeta_value.at<float>(i, j) = 0;
				}
				else
				{
					zeta_value.at<float>(i, j) = temp;
				}
			}
		}*/
		MatrixXf first_derivative(8, 1);
		for (int i = 1; i < roberts_img.rows - 1; i++)
		{
			for (int j = 1; j < roberts_img.cols - 1; j++)
			{
				first_derivative(0, 0) = 2 * (roberts_img.at<short>(i - 1, j - 1) - roberts_img.at<short>(i, j));
				first_derivative(1, 0) = 2 * (roberts_img.at<short>(i - 1, j) - roberts_img.at<short>(i, j));
				first_derivative(2, 0) = 2 * (roberts_img.at<short>(i - 1, j + 1) - roberts_img.at<short>(i, j));
				first_derivative(3, 0) = 2 * (roberts_img.at<short>(i, j - 1) - roberts_img.at<short>(i, j));
				first_derivative(4, 0) = 2 * (roberts_img.at<short>(i, j + 1) - roberts_img.at<short>(i, j));
				first_derivative(5, 0) = 2 * (roberts_img.at<short>(i + 1, j - 1) - roberts_img.at<short>(i, j));
				first_derivative(6, 0) = 2 * (roberts_img.at<short>(i + 1, j) - roberts_img.at<short>(i, j));
				first_derivative(7, 0) = 2 * (roberts_img.at<short>(i + 1, j + 1) - roberts_img.at<short>(i, j));
				
				MatrixXf direction = second_derivative * first_derivative;
				roberts_img.at<short>(i - 1, j - 1) = roberts_img.at<short>(i - 1, j - 1) - direction(0, 0);
				roberts_img.at<short>(i - 1, j) = roberts_img.at<short>(i - 1, j) - direction(1, 0);
				roberts_img.at<short>(i - 1, j + 1) = roberts_img.at<short>(i - 1, j + 1) - direction(2, 0);
				roberts_img.at<short>(i, j - 1) = roberts_img.at<short>(i, j - 1) - direction(3, 0);
				roberts_img.at<short>(i, j + 1) = roberts_img.at<short>(i, j + 1) - direction(4, 0);
				roberts_img.at<short>(i + 1, j - 1) = roberts_img.at<short>(i + 1, j - 1) - direction(5, 0);
				roberts_img.at<short>(i + 1, j) = roberts_img.at<short>(i + 1, j) - direction(6, 0);
				roberts_img.at<short>(i + 1, j + 1) = roberts_img.at<short>(i + 1, j + 1) - direction(7, 0);
				
			}
		}
		beta_current = beta_current * mu;
		time++;
	}
	Mat smooth_img;
	roberts_img.convertTo(smooth_img, CV_32FC1);
	return smooth_img;
}

 

展示效果

原始图像

 

 

 

原始的流场

平滑后的流场

原始图片

原始的流场

平滑后的流场