OpenCV 实现
话不多说，附上C++代码
cv::Mat L0Smoothing(cv::Mat &im8uc3, double lambda = 2e-2, double kappa = 2.0) { // convert the image to double format int 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: Image Denormin2.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/∂y cv::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)^2 double 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 subproblem for (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; }
// 矩阵的循环平移 void circshift(cv::Mat &A, int shift_row, int shift_col) { int row = A.rows, col = A.cols; // 考虑负数和超出边界值的情况 shift_row = (row + (shift_row % row)) % row; shift_col = (col + (shift_col % col)) % col; printf("shift_row = %d, shift_col = %d\n", shift_row, shift_col); cv::Mat temp = A.clone(); if (shift_row) { // temp[row - shift_row, row][*] -> A[0, shift_row][*] temp.rowRange(row - shift_row, row).copyTo(A.rowRange(0, shift_row)); // 上移 row - shift_row 个像素距离（移动 shift_row 个像素） // temp[0, row - shift_row][*] -> A[shift_row, row][*] temp.rowRange(0, row - shift_row).copyTo(A.rowRange(shift_row, row)); // 将原来的上部（row - shift_row 个像素）下移 shift_row 个像素距离 } if (shift_col) { temp.colRange(col - shift_col, col).copyTo(A.colRange(0, shift_col)); // 左移 temp.colRange(0, col - shift_col).copyTo(A.colRange(shift_col, col)); // 将原来的左部右移 } return; }
/* % PSF2OTF Convert point-spread function to optical transfer function. % OTF = PSF2OTF(PSF) computes the Fast Fourier Transform (FFT) of the % point-spread function (PSF) array and creates the optical transfer % function (OTF) array that is not influenced by the PSF off-centering. By % default, the OTF array is the same size as the PSF array. % % PSF2OTF 将点扩展函数转换为光传输函数。OTF = PSF2OTF(PSF) 计算点扩展函数 (PSF) 数组的快速傅里叶变换 (FFT)， % 创建不受 PSF 偏离中心影响的光传输函数(OTF)阵列。默认情况下，OTF 数组的大小与 PSF 数组相同。 % OTF = PSF2OTF(PSF,OUTSIZE) converts the PSF array into an OTF array of % specified size OUTSIZE. The OUTSIZE cannot be smaller than the PSF % array size in any dimension. % % OTF = PSF2OTF(PSF,OUTSIZE )将 PSF 数组转换为指定 OUTSIZE 大小的 OTF 数组。在任何维度上， OUTSIZE 都不能小于 PSF 数组的大小。 % To ensure that the OTF is not altered due to PSF off-centering, PSF2OTF % post-pads the PSF array (down or to the right) with zeros to match % dimensions specified in OUTSIZE, then circularly shifts the values of % the PSF array up (or to the left) until the central pixel reaches (1,1) % position. % 确保 OTF 不会因为 PSF 的偏离中心所改变，PSF2OTF 会向 PSF 数组后向填充 0 (向下或向右)来匹配指定的 OUTSIZE 的大小。 % 然后循环移位 PSF 数组的值(向左或向上)直到*像素达到 (1,1) 的位置。 % Note that this function is used in image convolution/deconvolution % when the operations involve the FFT. % % 注意，当涉及 FFT 操作时，此函数用于图像卷积/去卷积。 */ cv::Mat psf2otf(const cv::Mat &psf, const cv::Size &outSize) { cv::Size psfSize = psf.size(); cv::Mat new_psf = cv::Mat(outSize, CV_64FC2);// outSize 指定的大小，而且是两个通道， psf 只有一个通道 new_psf.setTo(0); // 只赋值到第一个通道 for (int i = 0; i < psfSize.height; i++) { for (int j = 0; j < psfSize.width; j++) { new_psf.at<cv::Vec2d>(i, j)[0] = psf.at<double>(i, j); // (0, 0) (0, 1)/(1, 0) } } // std::cout << //-1*int(floor(psfSize.height*0.5)) << //", " << //-1*int(floor(psfSize.width*0.5)) << //std::endl; // 循环平移矩阵（将矩阵中心移动到左上角） circshift(new_psf, -1 * int(floor(psfSize.height * .5f)), -1 * int(floor(psfSize.width * .5f))); // (0, -1) (-1, 0) cv::Mat otf; // 离散傅里叶变换，指定输出复数格式（默认是CCS格式） cv::dft(new_psf, otf, cv::DFT_COMPLEX_OUTPUT); return otf; }
效果图：
原图大小是640x640
内存是ddr4 2133MHZ
注意：
1.circshift, psf2otf 函数
circshift 的作用是循环平移矩阵；
psf2otf 的作用是将卷积核的中心移动到左上角，然后再进行FFT/DFT变换；
详细的psf2otf、circshift的解释
2. 二次规划问题：即问题的目标函数是二次函数
延伸阅读： P和NP问题的通俗解释
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);
效果如下：
时间只有第一个方法的一半不到！
CUDA 实现
我自己实现了一个基于 CUDA 的版本（Github 上实现的是错误的），处理一张640*640的图片，仅需要 447 ms （迭代相同次数的情况下）