python 拉普拉斯增强算法 图像增强拉普拉斯变换

拉普拉斯算子进行图像增强，以及算法优化

环境：vs2017 + OpenCV3.4.1

步骤：

（1）新建工程LapFilter

（2）确定项目阶段

（3）FFT变换部分

w = getOptimalDFTSize(gray_image.cols);//将输入图像延展到最佳尺寸，用0填充
	h = getOptimalDFTSize(gray_image.rows);//将输入图像延展到最佳尺寸，用0填充
	Mat padded;
	//将添加的像素用0填充
	copyMakeBorder(gray_image, padded, 0, h - gray_image.rows, 0, w - gray_image.cols, BORDER_CONSTANT, Scalar::all(0));
	padded.convertTo(padded, CV_32FC1);//转化为单通道浮点数
	//imshow("padded", padded);
	//将图像移到中心
	for (int i = 0; i < padded.rows; i++)
	{
		float *ptr = padded.ptr<float>(i);
		for (int j = 0; j < padded.cols; j++)	ptr[j] *= pow(-1, i + j);
	}
	//为傅里叶变换结果分配存储空间
	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
	Mat complexImg;
	//将plane数组合并成一个多通道数组complexImg
	merge(plane, 2, complexImg);
	//进行傅里叶变换
	dft(complexImg, complexImg);
	//laplace算子分配存储空间
	Mat laplace(padded.size(), CV_32FC2);

	for (int i = 0; i < padded.rows; i++)
	{
		float*k = laplace.ptr<float>(i);
		for (int j = 0; j < padded.cols; j++)
		{
			//laplace算子-4*PI*PI*（i-u/2）^2+(j-v/2)^2
			float d = pow(i - padded.rows / 2, 2) + pow(j - padded.cols / 2, 2);
			k[2 * j] = -4 * PI * PI * d;
			k[2 * j + 1] = -4 * PI * PI * d;
		}
	}
	//频域乘积H(u,v)*F(u,v)
	multiply(complexImg, laplace, laplace);
	//将多维数组拆分
	split(complexImg, plane);
	magnitude(plane[0], plane[1], plane[0]);//求幅值
	plane[0] += Scalar::all(1);
	//log(plane[0], plane[0]);
	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);//归一化，形成能够显示的结构
	imshow("dft", plane[0]);//dft图像

结果：

（4）laplace算子部分

//将H(u,v)*F(u,v)多维数组拆分
	split(laplace, plane);
	magnitude(plane[0], plane[1], plane[0]);//求幅值
	plane[0] += Scalar::all(1);
	//log(plane[0], plane[0]);
	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);//归一化，形成能够显示的结构
	imshow("laplace", plane[0]);//显示H(u,v)*F(u,v)

//将H(u,v)*F(u,v)做逆FFT
	idft(laplace, laplace);
	//将变换后的数组拆分
	split(laplace, plane);
	magnitude(plane[0], plane[1], plane[0]);//求幅值
	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);//归一化，形成能够显示的结构
	temp = plane[0];
imshow("idft-laplace", temp );

结果：

（5）输出图像

temp.convertTo(out, CV_8UC1);
	copyMakeBorder(gray_image, temp1, 0, h - gray_image.rows, 0, w - gray_image.cols, BORDER_CONSTANT, Scalar::all(0));
	output = temp1 - out;//g(x,y)= f(x,y)-(idft-laplace)
	imshow("output", output);
	//计算峰值信噪比
	double snr;
	snr = getPSNR(temp1, output);
	cout << "laplace_PSNR:" << snr << endl

结果：

（6）通过幂次变换进行优化

Mat matDst2;
//幂次变换
matDst2 = PowerTranseform(output, 1.4, 1);
matDst2.convertTo(matDst2, CV_8UC1);
double snr;
//计算峰值信噪比
snr = getPSNR(temp1, matDst2);
cout << "mici_after_PSNR:" << snr << endl;
imshow("幂次变换 y=1.4", matDst2);

优化结果：

结论：经过幂次变换，峰值信噪比降为原laplace算子的三分之一左右。

附程序：

w = getOptimalDFTSize(gray_image.cols);//将输入图像延展到最佳尺寸，用0填充
	h = getOptimalDFTSize(gray_image.rows);//将输入图像延展到最佳尺寸，用0填充
	Mat padded;
	//将添加的像素用0填充
	copyMakeBorder(gray_image, padded, 0, h - gray_image.rows, 0, w - gray_image.cols, BORDER_CONSTANT, Scalar::all(0));
	padded.convertTo(padded, CV_32FC1);//转化为单通道浮点数
	//imshow("padded", padded);
	//将图像移到中心
	for (int i = 0; i < padded.rows; i++)
	{
		float *ptr = padded.ptr<float>(i);
		for (int j = 0; j < padded.cols; j++)	ptr[j] *= pow(-1, i + j);
	}
	//为傅里叶变换结果分配存储空间
	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
	Mat complexImg;
	//将plane数组合并成一个多通道数组complexImg
	merge(plane, 2, complexImg);
	//进行傅里叶变换
	dft(complexImg, complexImg);
	//laplace算子分配存储空间
	Mat laplace(padded.size(), CV_32FC2);

	for (int i = 0; i < padded.rows; i++)
	{
		float*k = laplace.ptr<float>(i);
		for (int j = 0; j < padded.cols; j++)
		{
			//laplace算子-4*PI*PI*（i-u/2）^2+(j-v/2)^2
			float d = pow(i - padded.rows / 2, 2) + pow(j - padded.cols / 2, 2);
			k[2 * j] = -4 * PI * PI * d;
			k[2 * j + 1] = -4 * PI * PI * d;
		}
	}
	//频域乘积H(u,v)*F(u,v)
	multiply(complexImg, laplace, laplace);
	//将多维数组拆分
	split(complexImg, plane);
	magnitude(plane[0], plane[1], plane[0]);//求幅值
	plane[0] += Scalar::all(1);
	//log(plane[0], plane[0]);
	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);//归一化，形成能够显示的结构
	imshow("dft", plane[0]);//dft图像
	//将H(u,v)*F(u,v)多维数组拆分
	split(laplace, plane);
	magnitude(plane[0], plane[1], plane[0]);//求幅值
	plane[0] += Scalar::all(1);
	//log(plane[0], plane[0]);
	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);//归一化，形成能够显示的结构
	imshow("laplace", plane[0]);//显示H(u,v)*F(u,v)
	//将H(u,v)*F(u,v)做逆FFT
	idft(laplace, laplace);
	//将变换后的数组拆分
	split(laplace, plane);
	magnitude(plane[0], plane[1], plane[0]);//求幅值
	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);//归一化，形成能够显示的结构
	temp = plane[0];
imshow("idft-laplace", temp );
	temp1;
	temp.convertTo(out, CV_8UC1);
	copyMakeBorder(gray_image, temp1, 0, h - gray_image.rows, 0, w - gray_image.cols, BORDER_CONSTANT, Scalar::all(0));
	output = temp1 - out;//g(x,y)= f(x,y)-(idft-laplace)
	imshow("output", output);
	//计算峰值信噪比
	double snr;
	snr = getPSNR(temp1, output);
	cout << "laplace_PSNR:" << snr << endl;
Mat matDst2;
	//幂次变换
	matDst2 = PowerTranseform(output, 1.4, 1);
	matDst2.convertTo(matDst2, CV_8UC1);
	double snr;
	//计算峰值信噪比
	snr = getPSNR(temp1, matDst2);
	cout << "mici_after_PSNR:" << snr << endl;
	imshow("幂次变换 y=1.4", matDst2);
}
幂次变换函数：PowerTranseform
Mat CLapFilterDlg::PowerTranseform(Mat &src, double gamma, int parameter)
{ 
	Mat dst;
	dst.create(src.size(), src.type());
	vector<double> value;
	for (int r = 0; r < src.rows; r++)
	{
		uchar* srcRowData = src.ptr<uchar>(r);
		for (int c = 0; c < src.cols; c++)
		{
			//幂次变换的公式为 s = c * r ^ v	r为输入图像像素
			value.push_back(parameter * pow((double)srcRowData[c], gamma));
		}
	}
	double max = 0.0;
	double min = 0.0;
	for (int i = 0; i < value.size(); i++)
	{
		if (value[i] > max)
			max = value[i];
		if (value[i] < min)
			min = value[i];
	}
	int index = 0;
	for (int r = 0; r < dst.rows; r++)
	{
		uchar* dstRowData = dst.ptr<uchar>(r);
		for (int c = 0; c < dst.cols; c++)
		{
			dstRowData[c] = (uchar)(255 * ((value[index++] - min) * 1.0 / (max - min)));
		}
	}
	return dst;
}