灰度图L0最小梯度平滑opencv+C++代码

Matlab代码：​​​​​​http://t.csdn.cn/vX6Hj​​​​​​

对上述文章中的“代码一”进行C++转换

本文C++代码只针对灰度图像，且只进行了暴力转换，尚未对算法速度进行优化，还望各位大佬手下留情，如有需要可留言

代码如下：

#include<opencv2/opencv.hpp>
#include<iostream>
using namespace cv;
using namespace std;

//Matlab中psf2otf()函数
Mat psf2otf(const cv::Mat& psf, const cv::Size& outSize)
{
	Mat otf;

	if (countNonZero(psf) && !psf.empty())
	{
		Size psfSize = psf.size();

		Size paddSize = outSize - psfSize;

		copyMakeBorder(psf, otf, 0, paddSize.height,
			0, paddSize.width,
			BORDER_CONSTANT, Scalar(0));
		dft(otf, otf, DFT_COMPLEX_OUTPUT);
		
		//计算psf的元素个数
		int nElem = psf.cols * psf.rows;

		double nOps = 0;

		int nfft = nElem / psfSize.width;
		nOps += nfft + psfSize.width * log2(psfSize.width) * nfft;
		nfft = nElem / psfSize.height;
		nOps += nfft + psfSize.height * log2(psfSize.height) * nfft;
		//将复数otf分解成实部real和虚部imagin
		Mat planes[2];
		split(otf, planes);

		double max1, min1;
		cv::minMaxIdx(abs(planes[1]), &min1, &max1);
		//double imagin_max = mymax(abs(planes[1]));
		double imagin_max = max1;

		double max2, min2;
		cv::minMaxIdx(abs(otf), &min2, &max2);
		//double mag_max = mymax(myMagnitude(otf));
		double mag_max = max2;

		double eps = nOps * pow(2, -52);
		if ((imagin_max / mag_max) < eps)
			return planes[0];
	}
	return otf;
}
//Matlab中diff()函数
Mat diff(Mat& src, int a, int b)
{
	if (b == 2)//列数（宽）减少
	{
		Mat dst = Mat::zeros(Size(src.cols - 1, src.rows), src.type());
		for (int i = 0; i < dst.rows; i++)
		{
			double* p = src.ptr<double>(i);
			double* d = dst.ptr<double>(i);
			for (int j = 0; j < dst.cols; j++)
			{
				d[j] = (p[j + 1] - p[j]);
			}
		}
		return dst;	
	}
	else
	{
		Mat dst = Mat::zeros(Size(src.cols, src.rows - 1), src.type());
		for (int i = 0; i < dst.rows; i++)
		{
			double* p = src.ptr<double>(i);
			double* p_down = src.ptr<double>(i + 1);
			double* d = dst.ptr<double>(i);
			for (int j = 0; j < dst.cols; j++)
			{
				d[j] = (p_down[j] - p[j]);
			}
		}
		return dst;
	}
}
//L0梯度最小平滑
void L0(Mat& t_Im, float t_lambda, float t_kappa)
{
	double lambda = (double)t_lambda;
	double kappa = (double)t_kappa;

	//S = im2double(Im);
	Mat S;
	t_Im.convertTo(S, CV_64FC1, 1.0 / 255.0);

	//betamax = 1e5;
	double betamax = 1e5;
	/*
	* fx = [-1, 1];
	* fy = [-1; 1];
	* [N,M,D] = size(Im);
	* sizeI2D = [N,M];
	* otfFx = psf2otf(fx,sizeI2D);
	* otfFy = psf2otf(fy,sizeI2D);
	*/
	Mat fx = (Mat_<double>(1, 2) << 1, -1);
	Mat fy = (Mat_<double>(2, 1) << 1, -1);
	Mat otfFx = psf2otf(fx, t_Im.size());
	Mat otfFy = psf2otf(fy, t_Im.size());

	//Normin1 = fft2(S);
	Mat Normin1;
	dft(S, Normin1, DFT_COMPLEX_OUTPUT);

	//Denormin2 = abs(otfFx).^2 + abs(otfFy ).^2;
	Mat Denormin2(t_Im.rows, t_Im.cols, CV_64FC1);
	for (int i = 0; i < t_Im.rows; i++)
	{
		double* d = Denormin2.ptr<double>(i);
		for (int j = 0; j < t_Im.cols; j++)
		{
			Vec2d& c1 = otfFx.at<Vec2d>(i, j);
			Vec2d& c2 = otfFy.at<Vec2d>(i, j);
			// 0: Real, 1: Image
			d[j] = pow(c1[0], 2) + pow(c1[1], 2) + pow(c2[0], 2) + pow(c2[1], 2);
		}
	}

	//beta = 2*lambda;
	double beta = 2.0 * lambda;

	while (beta < betamax)
	{
		// Denormin   = 1 + beta*Denormin2;
		Mat Denormin = 1.0 + beta * Denormin2;

		// h-v subproblem
		
		//diff(S,1,2)
		Mat diffS12 = diff(S, 1, 2);
		/*
		* matlab中S(:, 1, : ), S中的第一列数据
		* matlab中S(:, end, : ), S中的最后一列数据
		*/
		Mat S01 = S.colRange(0, 1);
		Mat S0end = S.colRange(S.cols - 1, S.cols);
		/*
		* h = [diff(S,1,2), S(:,1,:) - S(:,end,:)];
		*/
		Mat Scol = S01 - S0end;
		Mat h = Mat::zeros(S.size(), CV_64FC1);
		for (int i = 0; i < S.rows; i++)
		{
			double* p = diffS12.ptr<double>(i);
			double* sp = Scol.ptr<double>(i);
			double* d = h.ptr<double>(i);
			for (int j = 0; j < S.cols; j++)
			{
				if (j < diffS12.cols)
				{
					d[j] = p[j];
				}
				else
				{
					d[j] = sp[0];
				}
			}
		}

		//diff(S,1,1)
		Mat diffS11 = diff(S, 1, 1);
		/*
		* matlab中S(1, :, : ), S中的第一行数据
		* matlab中S(end, :, : ), S中的最后一行数据
		*/
		Mat S10 = S.rowRange(0, 1);
		Mat Send0 = S.rowRange(S.rows - 1, S.rows);
		/*
		* v = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
		*/
		Mat Srow = S10 - Send0;
		Mat v = Mat::zeros(S.size(), CV_64FC1);
		for (int i = 0; i < S.rows; i++)
		{
			if (i < diffS11.rows)
			{
				double* p = diffS11.ptr<double>(i);
				double* d = v.ptr<double>(i);
				for (int j = 0; j < S.cols; j++)
				{
					d[j] = p[j];
				}
			}
			else
			{
				double* sp = Srow.ptr<double>(0);
				double* d = v.ptr<double>(i);
				for (int j = 0; j < S.cols; j++)
				{
					d[j] = sp[j];
				}
			}
		}
		/*
		* matlab中：
	    *   if D==1
				t = (h.^2+v.^2)<lambda/beta;
			else
				t = sum((h.^2+v.^2),3)<lambda/beta;
				t = repmat(t,[1,1,D]);
			end
			h(t)=0; v(t)=0;
		*/
		for (int i = 0; i < S.rows; i++)
		{
			double* ht = h.ptr<double>(i);
			double* vt = v.ptr<double>(i);
			for (int j = 0; j < S.cols; j++)
			{
				// (∂S/∂x)^2 + (∂S/∂y)^2			
				double val = pow(ht[j], 2) + pow(vt[j], 2);
				// (∂S/∂x)^2 + (∂S/∂y)^2  < λ/β
				if (val < lambda / beta) {
					ht[j] = 0.0;
					vt[j] = 0.0;
				}
			}
		}
		// S subproblem

		//Normin2 = [h(:,end,:) - h(:, 1,:), -diff(h,1,2)];
		Mat diffh12 = diff(h, 1, 2);
		/*
		* matlab中h(:, 1,:), h中的第一列数据
		* matlab中h(:,end,:), h中的最后一列数据
		*/
		Mat h0end = h.colRange(h.cols - 1, h.cols);
		Mat h01 = h.colRange(0, 1);
		Mat hcol = h0end - h01;
		Mat Normin21 = Mat::zeros(S.size(), CV_64FC1);
		for (int i = 0; i < S.rows; i++)
		{
			double* p = diffh12.ptr<double>(i);
			double* sp = hcol.ptr<double>(i);
			double* d = Normin21.ptr<double>(i);
			for (int j = 0; j < S.cols; j++)
			{
				if (j == 0)
				{
					d[j] = sp[0];
				}
				else
				{
					d[j] = -p[j - 1];
				}
			}
		}
		//Normin2 = Normin2 + [v(end,:,:) - v(1, :,:); -diff(v,1,1)];
		/*
		* matlab中v(1, :,:), v中的第一行数据
		* matlab中v(end,:,:), v中的最后一行数据
		*/
		Mat diffv11 = diff(v, 1, 1);

		Mat v10 = v.rowRange(0, 1);
		Mat vend0 = v.rowRange(v.rows - 1, v.rows);
		Mat vrow = vend0 - v10;
		Mat Normin22 = Mat::zeros(S.size(), CV_64FC1);
		for (int i = 0; i < S.rows; i++)
		{
			if (i == 0)
			{
				double* sp = vrow.ptr<double>(0);
				double* d = Normin22.ptr<double>(i);
				for (int j = 0; j < S.cols; j++)
				{
					d[j] = sp[j];
				}
			}
			else
			{
				double* p = diffv11.ptr<double>(i - 1);
				double* d = Normin22.ptr<double>(i);
				for (int j = 0; j < S.cols; j++)
				{
					d[j] = -p[j];
				}
			}
		}
		Mat Normin2 = Normin21 + Normin22;

		//fft2(Normin2)
		Mat fft2Normin2;
		dft(Normin2, fft2Normin2, DFT_COMPLEX_OUTPUT);

		//FS = (Normin1 + beta*fft2(Normin2))./Denormin;
		Mat FS = Normin1 + beta * fft2Normin2;
		for (int i = 0; i < S.rows; i++)
		{
			double* p = Denormin.ptr<double>(i);
			for (int j = 0; j < S.cols; j++)
			{
				FS.at<Vec2d>(i, j)[0] /= p[j];
		     	FS.at<Vec2d>(i, j)[1] /= p[j];
			}
		}
		//ifft2(FS)
		Mat ifft;
		idft(FS, ifft, cv::DFT_SCALE | cv::DFT_COMPLEX_OUTPUT);
		for (int i = 0; i < S.rows; i++)
		{
			double* p = S.ptr<double>(i);
			for (int j = 0; j < S.cols; j++)
			{
				//S = real(ifft2(FS));
				p[j] = ifft.at<cv::Vec2d>(i, j)[0];
			}
		}
		//beta = beta*kappa;
		beta = beta * kappa;
	}
	//显示图片
	normalize(S, S, 0, 255, NORM_MINMAX);
	S.convertTo(S, CV_8UC1);
	resize(S, S, Size(S.cols / 1, S.rows / 1));
	imshow("S", S);
}
int main()
{
    Mat src = imread("C:\\Users\\DELL\\Desktop\\Lena.jpg", 0);
    L0(src,0.02,2.0);
	waitKey(0);
	return 0;
}

选取Lena图像，原图、Matlab与C++结果图分别如下：

原图：

 Matlab的L0梯度最小平滑结果图：

  C++的L0梯度最小平滑结果图：