视频网址:https://www.bilibili.com/video/av35811232
完整代码:https://download.csdn.net/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;
}