opencv逆滤波程序的问题

[longbozhang0918]

楼主小白一枚，综合网上找的一些现有程序写了一个逆滤波的程序，但是复原结果一团黑，找不着问题在哪，求大神帮助。测试图像是用matlab加的模糊尺度20，角度45的lena图像，代码如下。

#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <stdio.h>
#include<math.h>
#include <cv.h> 
using namespace cv;
using namespace std;
const double EPS = 1e-12;
const double PI=3.141592653;
void generatePSF(Mat &psf, double len,double angle)
{
	double half=len/2;											
	double alpha = (angle-floor(angle/ 180) *180) /180* PI;
	double cosalpha = cos(alpha);
	double sinalpha = sin(alpha);
	int xsign;
	if (cosalpha < 0)
	{
		xsign = -1;
	}
	else
	{
		if (angle == 90)
		{
			xsign = 0;
		}
		else
		{
			xsign = 1;
		}
	}
	int psfwdt = 1;
	int sx = (int)fabs(half*cosalpha + psfwdt*xsign - len*EPS);
	int sy = (int)fabs(half*sinalpha + psfwdt - len*EPS);
	Mat_<double> psf1(sy, sx, CV_64F);
	Mat_<double> psf2(sy * 2, sx * 2, CV_64F);
	int row = 2 * sy;
	int col = 2 * sx;
	/*为减小运算量，先计算一半大小的PSF*/
	for (int i = 0; i < sy; i++)
	{
		double* pvalue = psf1.ptr<double>(i);
		for (int j = 0; j < sx; j++)
		{
			pvalue[j] = i*fabs(cosalpha) - j*sinalpha;

			double rad = sqrt((double)(i*i + j*j));
			if (rad >= half && fabs(pvalue[j]) <= psfwdt)
			{ 
				double temp = half - fabs((j + pvalue[j] * sinalpha) / cosalpha);
				pvalue[j] = sqrt(pvalue[j] * pvalue[j] + temp*temp);
			}
			pvalue[j] = psfwdt + EPS - fabs(pvalue[j]);
			if (pvalue[j] < 0)
			{
				pvalue[j] = 0;
			}
		}
	}
	/*将模糊核矩阵扩展至实际大小*/
	for (int i = 0; i < sy; i++)
	{
		double* pvalue1 = psf1.ptr<double>(i);
		double* pvalue2 = psf2.ptr<double>(i);
		for (int j = 0; j < sx; j++)
		{
			pvalue2[j] = pvalue1[j];
		}
	}
	
	
	//保持图像总能量不变，归一化矩阵
	double sum = 0;
	for (int i = 0; i < row; i++)
	{
		for (int j = 0; j < col; j++)
		{
			sum+= psf2[i][j];
		}
	}
	psf2 = psf2 / sum;
	if (cosalpha>0)
	{
		flip(psf2, psf2, 0);
	}

	//cout << "psf2=" << psf2 << endl;
	psf = psf2;


}
int dft2(Mat &complexImg, Mat input)
{

    Mat img = input;
	//	imread("C:/Users/JohnSmith/Desktop/lena.bmp", CV_LOAD_IMAGE_GRAYSCALE);
    if( img.empty() )
    {
        printf("Cannot read image file: %s\n", "C:/Users/JohnSmith/Desktop/lena.bmp");
        return -1;
    }
//这里获取符合快速傅里叶变换（FFT）的大小，m和n都可以分解为2、3、5相乘的组合，参见 注1
    int M = getOptimalDFTSize( img.rows );
    int N = getOptimalDFTSize( img.cols );
    Mat padded;
//将原图像的大小变为m*n的大小，补充的位置填0，
    copyMakeBorder(img, padded, (M - img.rows)/2, (M - img.rows)/2, (N - img.cols)/2, (N - img.cols)/2, BORDER_CONSTANT, Scalar::all(0));
//这里是获取了两个mat，一个用于存放dft变换的实部，一个用于存放虚部，初始的时候，实部就是图像本身，虚部全为0
    Mat planes[] = {Mat_<double>(padded), Mat::zeros(padded.size(), CV_64F)};
//将几个单通道的mat融合成一个多通道的mat，这里融合的complexImg即有实部，又有虚部
    merge(planes, 2, complexImg);
//dft变换，因为complexImg本身就是两个通道的mat，所以dft变换的结果也可以保存在其中
    dft(complexImg, complexImg);
	
}
int main(int argc, const char ** argv)
{
	Mat blur=imread("lena2045.bmp", CV_LOAD_IMAGE_GRAYSCALE);
	Mat motion,mpadded;
	Mat recover;
    Mat h,g,f;
	
	double a,b,c,d;
	int i,j;
    int isizex=getOptimalDFTSize(blur.rows);
	int isizey=getOptimalDFTSize(blur.cols);
	generatePSF(motion,20,45);
	
	copyMakeBorder(motion, mpadded, (isizex - motion.rows)/2, (isizex - motion.rows)/2, (isizey - motion.cols)/2, (isizey - motion.cols)/2, BORDER_CONSTANT, Scalar::all(0));	
	
	dft2(h,mpadded);
	dft2(g,blur);
	
	Mat hplanes[] = {Mat_<double>(mpadded), Mat::zeros(mpadded.size(), CV_64F)};
	Mat gplanes[] = {Mat_<double>(mpadded), Mat::zeros(mpadded.size(), CV_64F)};
	Mat fplanes[] = {Mat_<double>(mpadded), Mat::zeros(mpadded.size(), CV_64F)};
	
	split(h, hplanes);
	split(g, gplanes);
	

	for(i=0;i<isizex;i++)
	{
	  for(j=0;j<isizey;j++)
	  {
		a=gplanes[0].at<double>(i,j);
		b=gplanes[1].at<double>(i,j);
        c=hplanes[0].at<double>(i,j);
		d=hplanes[0].at<double>(i,j);
	    fplanes[0].at<double>(i,j)=(a*c+b*d)/(c*c+d*d);
		fplanes[1].at<double>(i,j)=(b*c-a*d)/(c*c+d*d);
	  }
	}
	
	merge(fplanes, 2, f);
	
idft(f,recover,DFT_SCALE|DFT_REAL_OUTPUT);
	recover.convertTo(recover,CV_8U);

	imshow("模糊图像",blur);
	imshow("复原结果",recover);
	 waitKey();
    return 0;
}

;

[赵4老师]

代码功能归根结底不是别人帮自己看或讲解或注释出来的；而是被自己静下心来花足够长的时间和精力亲自动手单步或设断点或对执行到某步获得的中间结果显示或写到日志文件中一步一步分析出来的。 提醒：再牛×的老师也无法代替学生自己领悟和上厕所！ 单步调试和设断点调试（VS IDE中编译连接通过以后，按F10或F11键单步执行，按Shift+F11退出当前函数；在某行按F9设断点后按F5执行停在该断点处。）是程序员必须掌握的技能之一。

[赵4老师]

请检查每个函数调用的返回值。

[Mr_Zhouzl]

单步加断点调试检查每个接口的数值

[cattpon]

learning~

[hugh_z]

66666666666666666

[shiter]

能否简短的解释一下逆滤波？你期望的效果是啥

[nettman]

[qq_36769369]

66666666666666666666666666666666

[cattpon]

learning~

[WowMusic]

好久没来CSDN，现在都在推图片处理的帖子啊

[hugh_z]

6666666666666666

[hugh_z]

666666666666666

[赵4老师]

引用 11 楼 WowMusic 的回复:

好久没来CSDN，现在都在推图片处理的帖子啊

没准再过几个月，都推AI的帖子了。

[meadow]

楼主说的是反卷积吗， 你看看最后图像像素值的尺度范围 现在是毕业季

[nettman]

[luqiang6q]

666666666666

[qq_34650785]

666666666666