MATLAB中imfilter函数用Opencv 实现遇到的若干问题

文章出处：http://blog.csdn.net/wx123456123/article/details/37658377

最近在做科研的过程中用到图像的滤波，由于图像数据量比较大，就考虑空域转化频域相乘来提高效率，测试的的时候出现MATLAB 与OPENCV的结果不相同，以前也知道MATLAB和OPENCV的滤波机理不相同，但是具体的过程也是不是太清楚，于是用程序测试filter2,conv2,imfilter2,cvFilter2D函数的异同。

conv2用来计算源图像A和滤波模板H的卷积，计算过程为将滤波模板旋转180度后，进行掩膜计算，计算的过程中要对边界补0计算。

[html]  view plain  copy

1. a=[1 2;3 4];  
2. b=[-1 1;-2 2];  
3. c=conv2(a,b,'full')  

输出为

conv2的计算可以使用频域的乘积计算，但是要对矩阵的边界补零

[plain]  view plain  copy

1. a1=[1 2 0;3 4 0;0 0 0];  
2. b1=[-1 1 0;-2 2 0;0 0 0];  
3. a1f=fft2(a1);  
4. b1f=fft2(b1);  
5. abf1=a1f.*b1f;  
6. c1=ifft2(abf1)  

输出为

filter2计算过程为对滤波模板H进行180度旋转后，再调用conv2函数。由于conv2函数又对模板旋转180度，等于直接对滤波模板H直接进行掩膜计算。

[html]  view plain  copy

1. a=[1 2;3 4];  
2. b=[-1 1;-2 2];  
3. c2=filter2(a,b,'full')  
4. c3=filter2(b,a,'full')  

输出结果：

可以看出filter2函数已经不满足卷积的可交换律了

imfilter则是对函数filter2和conv2函数的综合，当选择参数conv时则，进行卷积计算与函数conv2函数相同， 选择corr参数时，进行相关运算和filter2过程相同

 opencv中的滤波函数  cvFilter2D( const CvArr* src, CvArr* dst, const CvMat* kernel,  CvPoint anchor CV_DEFAULT(cvPoint(-1,-1)))可以实现图像的卷及滤波，但是cvFilter2D处理时不对滤波核进行180度旋转，边界处理时实行镜像处理

[plain]  view plain  copy

1. #include "iostream"  
2. #include "cv.h"  
3. #include "highgui.h"  
4. using namespace std;  
5.   
6. void speedy_convolution(const CvMat *A,const CvMat *B,CvMat *C);  
7. void ShowMat(CvMat *m);  
8.   
9. int main()  
10. {    
11.     float A1[]={1,2,3,4};float B1[]={-1,1,-2,2};  
12.     //float A[]={1,2,0,3,4,0,0,0,0};  
13.     //float B[]={-1,1,0,-2,2,0,0,0,0};  
14.     CvMat Ma1=cvMat(2,2,CV_32FC1,&A1);  
15.     CvMat Mb1=cvMat(2,2,CV_32FC1,&B1);  
16.       
17.     CvMat *C1=cvCreateMat(3,3,CV_32FC1);  
18.     CvMat *C=cvCreateMat(2,2,CV_32FC1);  
19.     double t=cvGetTickCount();  
20.     speedy_convolution(&Ma1,&Mb1,C1);//  
21.     t=cvGetTickCount()-t;  
22.     cout<<"fft时间："<<t/((double)cvGetTickFrequency()*1000.)<<endl;  
23.     cout<<"FFT卷积计算结果:"<<endl;  
24.     ShowMat(C1);  
25.   
26.     t=cvGetTickCount();  
27.     cvFilter2D(&Mb1,C,&Ma1,cvPoint(1,1));  
28.     t=cvGetTickCount()-t;  
29.     cout<<endl<<"cvFilter2D时间："<<t/((double)cvGetTickFrequency()*1000.)<<endl;  
30.     cout<<"cvFilter2D函数测试结果:"<<endl;  
31.     ShowMat(C);   
32.     cout<<endl;  
33.       
34. }    
35. void ShowMat(CvMat *m)  
36. {  
37.     int i,j;  
38.     for (i=0;i<m->rows;i++)  
39.     {  
40.         for (j=0;j<m->cols;j++)  
41.         {  
42.             cout<<cvGet2D(m,i,j).val[0]<<"  ";  
43.         }  
44.         cout<<endl;  
45.     }  
46.   
47. }  
48. void speedy_convolution(  
49.     const CvMat *A,//size:M1 * N1  
50.     const CvMat *B, //size : M2 * N2  
51.     CvMat *C    //size:(M1+M2-1) * (N1+N2-1)  
52.     )  
53. {  
54.     int dft_M = cvGetOptimalDFTSize(A->rows + B->rows -1);  
55.     int dft_N = cvGetOptimalDFTSize(A->cols + B->cols -1);  
56.   
57.     CvMat *dft_A = cvCreateMat(dft_M , dft_N , A->type);  
58.     CvMat *dft_B = cvCreateMat(dft_M, dft_N, B->type);  
59.     CvMat tmp;  
60.   
61.     //copy A TO dft_A and pad dft_A with zeros  
62.     cvGetSubRect(dft_A,&tmp , cvRect(0,0,A->cols,A->rows));  
63.     cvCopy(A,&tmp);  
64.     cvGetSubRect(dft_A,&tmp,cvRect(A->cols,0,dft_A->cols - A->cols , A->rows));  
65.     cvZero(&tmp);  
66.   
67.     //no need to pad bottom part of dft_A with zeros because of   
68.     //use nonzero_rows parameter in cvDFT() call below  
69.     cvDFT(dft_A,dft_A,CV_DXT_FORWARD,A->rows);  
70.   
71.     //repeat the same with the second array  
72.     cvGetSubRect(dft_B,&tmp,cvRect(0,0,B->cols,B->rows));  
73.     cvCopy(B,&tmp);  
74.     cvGetSubRect(dft_B,&tmp,cvRect(B->cols,0,dft_B->cols - B->cols , B->rows));  
75.     cvZero(&tmp);  
76.   
77.     //no need to pad bottom part of dft_B with zeros because of   
78.     //use nonzero-rows parameter in cvDFT() call below  
79.     cvDFT(dft_B,dft_B,CV_DXT_FORWARD,B->rows);  
80.   
81.     //or CV_DXT_MUL_CONJ to get correlation rather than convolution   
82.     cvMulSpectrums(dft_A,dft_B,dft_A,0);  
83.   
84.     //calculate only the top part  
85.     cvDFT(dft_A,dft_A,CV_DXT_INV_SCALE,C->rows);  
86.     cvGetSubRect(dft_A,&tmp,cvRect(0,0,C->cols,C->rows));  
87.     cvCopy(&tmp,C);  
88.     cvReleaseMat(&dft_A);  
89.     cvReleaseMat(&dft_B);  
90. }  

上面代码：函数speedy_convolution代码参考http://blog.csdn.net/fdl19881/article/details/6716202

测试结果：
        

 cvFilter2D 处理过程对Mb1的边界镜像补数为：

 

 再以Ma1：的cvPoint（1,1）为滤波核为中心进行掩膜计算

一般图像滤波的模板都是关于中心对称的，因此在计算图像卷积滤波时，是否对卷积模板旋转180度关系不是太大。

由于矩阵比较小，fft的效果不是太明显，如果是图像的维数比较大的话，效率的提高就比较可观了

 

以上只是本人的个人感悟，不当之处希望和大家一起进步。