L0 范数图像平滑(L0 Smooth) 代码及详细注释 【OpenCV】

原理可以参考原作者的论文 以及 这位大神写的博客

话不多说，附上C++代码

cv::Mat L0Smoothing(cv::Mat &im8uc3, double lambda = 2e-2, double kappa = 2.0) {// convert the image to double formatint row = im8uc3.rows, col = im8uc3.cols;cv::Mat S;im8uc3.convertTo(S, CV_64FC3, 1. / 255.);        // 2*2 的卷积核cv::Mat fx(1, 2, CV_64FC1);cv::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;// 把一个空间点扩散函数转换为频谱面的光学传递函数cv::Size sizeI2D = im8uc3.size();cv::Mat otfFx = psf2otf(fx, sizeI2D);cv::Mat otfFy = psf2otf(fy, sizeI2D);// FNormin1 = fft2(S);cv::Mat FNormin1[3];// 注意：DFT以后，FNormal为两个通道cv::Mat single_channel[3];cv::split(S, single_channel);// 分裂为三个通道for (int k = 0; k < 3; k++) {// 离散傅里叶变换，指定输出复数格式（默认是CCS格式）cv::dft(single_channel[k], FNormin1[k], cv::DFT_COMPLEX_OUTPUT);}// F(∂x)∗F(∂x)+F(∂y)∗F(∂y);cv::Mat Denormin2(row, col, CV_64FC1);for (int i = 0; i < row; i++) {for (int j = 0; j < col; j++) {cv::Vec2d &c1 = otfFx.at<cv::Vec2d>(i, j);cv::Vec2d &c2 = otfFy.at<cv::Vec2d>(i, j);// 0: Real, 1: ImageDenormin2.at<double>(i, j) = sqr(c1[0]) + sqr(c1[1]) + sqr(c2[0]) + sqr(c2[1]);}}double beta = 2.0*lambda;double betamax = 1e5;while (beta < betamax) {// F(1)+β(F(∂x)∗F(∂x)+F(∂y)∗F(∂y))cv::Mat Denormin = 1.0 + beta*Denormin2;// h-v subproblem// 三个通道的 ∂S/∂x ∂S/∂ycv::Mat dx[3], dy[3];for (int k = 0; k < 3; k++) {cv::Mat shifted_x = single_channel[k].clone();circshift(shifted_x, 0, -1);dx[k] = shifted_x - single_channel[k];cv::Mat shifted_y = single_channel[k].clone();circshift(shifted_y, -1, 0);dy[k] = shifted_y - single_channel[k];}for (int i = 0; i < row; i++) {for (int j = 0; j < col; j++) {// (∂S/∂x)^2 + (∂S/∂y)^2double val =sqr(dx[0].at<double>(i, j)) + sqr(dy[0].at<double>(i, j)) +sqr(dx[1].at<double>(i, j)) + sqr(dy[1].at<double>(i, j)) +sqr(dx[2].at<double>(i, j)) + sqr(dy[2].at<double>(i, j));// (∂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 subproblemfor (int k = 0; k < 3; k++) {// 二阶导cv::Mat shift_dx = dx[k].clone();circshift(shift_dx, 0, 1);cv::Mat ddx = shift_dx - dx[k];cv::Mat shift_dy = dy[k].clone();circshift(shift_dy, 1, 0);cv::Mat ddy = shift_dy - dy[k];cv::Mat Normin2 = ddx + ddy;cv::Mat FNormin2;// 离散傅里叶变换，指定输出复数格式（默认是CCS格式）cv::dft(Normin2, FNormin2, cv::DFT_COMPLEX_OUTPUT);// F(g)+β(F(∂x)∗F(h)+F(∂y)∗F(v))cv::Mat FS = FNormin1[k] + beta*FNormin2;// 论文的公式(8)：F^-1括号中的内容for (int i = 0; i < row; i++) {for (int j = 0; j < col; j++) {FS.at<cv::Vec2d>(i, j)[0] /= Denormin.at<double>(i, j);FS.at<cv::Vec2d>(i, j)[1] /= Denormin.at<double>(i, j);//std::cout<< FS.at<cv::Vec2d>(i, j)[0] / Denormin.at<double>(i, j) <<std::endl;}}// 论文的公式(8)：傅里叶逆变换cv::Mat ifft;cv::idft(FS, ifft, cv::DFT_SCALE | cv::DFT_COMPLEX_OUTPUT);for (int i = 0; i < row; i++) {for (int j = 0; j < col; j++) {single_channel[k].at<double>(i, j) = ifft.at<cv::Vec2d>(i, j)[0];//std::cout<< ifft.at<cv::Vec2d>(i, j)[0] <<std::endl;}}}beta *= kappa;std::cout << '.';}cv::merge(single_channel, 3, S);return S;}

效果图：

原图大小是640x640

内存是ddr4 2133MHZ

注意：

1.

circshift 的作用是循环平移矩阵；

psf2otf 的作用是将卷积核的中心移动到左上角，然后再进行FFT/DFT变换；

详细的psf2otf、circshift的解释

2. 

二次规划问题：即问题的目标函数是二次函数

延伸阅读： P和NP问题的通俗解释

3.OpenCV 3.x.x的扩展模块（ximgproc. Extended Image Processing）提供了l0smooth的API

附：在Windows下编译扩展OpenCV 3.1.0 + opencv_contrib ，值得注意的是，在编译OpenCV310源码时，opencv_aruco 模块的源码有bug，需要将

calibrateCamera(processedObjectPoints, processedImagePoints, imageSize, _cameraMatrix,                           _distCoeffs, _rvecs, _tvecs, _stdDeviationsIntrinsics, _stdDeviationsExtrinsics,                           _perViewErrors, flags, criteria);

修改为

calibrateCamera(processedObjectPoints, processedImagePoints, imageSize, _cameraMatrix,                           _distCoeffs, _rvecs, _tvecs, /*_stdDeviationsIntrinsics, _stdDeviationsExtrinsics,                           _perViewErrors,*/ flags, criteria);

main函数代码：

cv::Mat src = cv::imread("D:/Pictures/beard.jpg", 1);cv::Mat dst;int64 begin = cvGetTickCount();cv::ximgproc::l0Smooth(src, dst, 0.01, 2.);int64 end = cvGetTickCount();float time = (end - begin) / (cvGetTickFrequency() * 1000.);printf("time = %fms\n", time);imshow("l0Smooth", dst);cv::waitKey(0);

效果如下：

时间只有第一个方法的一半不到！