Gabor的OpenCV代码

转自http://blog.csdn.net/guoming0000/article/details/7839917

最近弄人脸识别，用到Gabor卷积核，但网上的代码似乎没有和我心意的，于是参考了自己写了下！参考了Zhou Mian以及matlab的Gabor实现代码的代码。虽然OpenCV的imporc下面有个gabor.cpp，但那个是一般形式的公式，不是用来做人脸识别的，可以参考文献A review on Gabor wavelets for face recognition，又说到。上代码和链接地址！下载地址~

   目前代码未经过更多的测试，不少功能为加入，但可以满足许多人的使用和参考了吧，很多人肯定非常非常需要，先开源下，欢迎指出错误之处。

[cpp] view plaincopyprint?
  1. //GaborFR.h  
[cpp] view plaincopyprint?
  1. #pragma once  
  2. #include "opencv2\opencv.hpp"  
  3. #include <vector>  
  4. using namespace std;  
  5. using namespace cv;  
  6. class GaborFR  
  7. {  
  8. public:  
  9.     GaborFR();  
  10.     static Mat  getImagGaborKernel(Size ksize, double sigma, double theta,   
  11.                                     double nu,double gamma=1, int ktype= CV_32F);  
  12.     static Mat  getRealGaborKernel( Size ksize, double sigma, double theta,   
  13.                                     double nu,double gamma=1, int ktype= CV_32F);  
  14.     static Mat  getPhase(Mat &real,Mat &imag);  
  15.     static Mat  getMagnitude(Mat &real,Mat &imag);  
  16.     static void getFilterRealImagPart(Mat& src,Mat& real,Mat& imag,Mat &outReal,Mat &outImag);  
  17.     static Mat  getFilterRealPart(Mat& src,Mat& real);  
  18.     static Mat  getFilterImagPart(Mat& src,Mat& imag);  
  19.   
  20.     void        Init(Size ksize=Size(19,19), double sigma=2*CV_PI,  
  21.                     double gamma=1, int ktype=CV_32FC1);  
  22. private:  
  23.     vector<Mat> gaborRealKernels;  
  24.     vector<Mat> gaborImagKernels;  
  25.     bool isInited;  
  26. };  

[cpp] view plaincopyprint?
  1. #include "stdafx.h"  
  2. #include "GaborFR.h"  
  3. GaborFR::GaborFR()  
  4. {  
  5.     isInited = false;  
  6. }  
  7. void GaborFR::Init(Size ksize, double sigma,double gamma, int ktype)  
  8. {  
  9.     gaborRealKernels.clear();  
  10.     gaborImagKernels.clear();  
  11.     double mu[8]={0,1,2,3,4,5,6,7};  
  12.     double nu[5]={0,1,2,3,4};  
  13.     int i,j;  
  14.     for(i=0;i<5;i++)  
  15.     {  
  16.         for(j=0;j<8;j++)  
  17.         {  
  18.             gaborRealKernels.push_back(getRealGaborKernel(ksize,sigma,mu[j]*CV_PI/8,nu[i],gamma,ktype));  
  19.             gaborImagKernels.push_back(getImagGaborKernel(ksize,sigma,mu[j]*CV_PI/8,nu[i],gamma,ktype));  
  20.         }  
  21.     }  
  22.     isInited = true;  
  23. }  
  24. Mat GaborFR::getImagGaborKernel(Size ksize, double sigma, double theta, double nu,double gamma, int ktype)  
  25. {  
  26.     double  sigma_x     = sigma;  
  27.     double  sigma_y     = sigma/gamma;  
  28.     int     nstds       = 3;  
  29.     double  kmax        = CV_PI/2;  
  30.     double  f           = cv::sqrt(2.0);  
  31.     int xmin, xmax, ymin, ymax;  
  32.     double c = cos(theta), s = sin(theta);  
  33.     if( ksize.width > 0 )  
  34.     {  
  35.         xmax = ksize.width/2;  
  36.     }  
  37.     else//这个和matlab中的结果一样，默认都是19 !  
  38.     {  
  39.         xmax = cvRound(std::max(fabs(nstds*sigma_x*c), fabs(nstds*sigma_y*s)));  
  40.     }  
  41.     if( ksize.height > 0 )  
  42.     {  
  43.         ymax = ksize.height/2;  
  44.     }  
  45.     else  
  46.     {  
  47.         ymax = cvRound(std::max(fabs(nstds*sigma_x*s), fabs(nstds*sigma_y*c)));  
  48.     }  
  49.     xmin = -xmax;  
  50.     ymin = -ymax;  
  51.     CV_Assert( ktype == CV_32F || ktype == CV_64F );  
  52.     float*  pFloat;  
  53.     double* pDouble;  
  54.     Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);  
  55.     double k        =   kmax/pow(f,nu);  
  56.     double scaleReal=   k*k/sigma_x/sigma_y;  
  57.     for( int y = ymin; y <= ymax; y++ )  
  58.     {  
  59.         if( ktype == CV_32F )  
  60.         {  
  61.             pFloat = kernel.ptr<float>(ymax-y);  
  62.         }  
  63.         else  
  64.         {  
  65.             pDouble = kernel.ptr<double>(ymax-y);  
  66.         }  
  67.         for( int x = xmin; x <= xmax; x++ )  
  68.         {  
  69.             double xr = x*c + y*s;  
  70.             double v = scaleReal*exp(-(x*x+y*y)*scaleReal/2);  
  71.             double temp=sin(k*xr);  
  72.             v   =  temp*v;  
  73.             if( ktype == CV_32F )  
  74.             {  
  75.                 pFloat[xmax - x]= (float)v;  
  76.             }  
  77.             else  
  78.             {  
  79.                 pDouble[xmax - x] = v;  
  80.             }  
  81.         }  
  82.     }  
  83.     return kernel;  
  84. }  
  85. //sigma一般为2*pi  
  86. Mat GaborFR::getRealGaborKernel( Size ksize, double sigma, double theta,   
  87.     double nu,double gamma, int ktype)  
  88. {  
  89.     double  sigma_x     = sigma;  
  90.     double  sigma_y     = sigma/gamma;  
  91.     int     nstds       = 3;  
  92.     double  kmax        = CV_PI/2;  
  93.     double  f           = cv::sqrt(2.0);  
  94.     int xmin, xmax, ymin, ymax;  
  95.     double c = cos(theta), s = sin(theta);  
  96.     if( ksize.width > 0 )  
  97.     {  
  98.         xmax = ksize.width/2;  
  99.     }  
  100.     else//这个和matlab中的结果一样，默认都是19 !  
  101.     {  
  102.         xmax = cvRound(std::max(fabs(nstds*sigma_x*c), fabs(nstds*sigma_y*s)));  
  103.     }  
  104.   
  105.     if( ksize.height > 0 )  
  106.         ymax = ksize.height/2;  
  107.     else  
  108.         ymax = cvRound(std::max(fabs(nstds*sigma_x*s), fabs(nstds*sigma_y*c)));  
  109.     xmin = -xmax;  
  110.     ymin = -ymax;  
  111.     CV_Assert( ktype == CV_32F || ktype == CV_64F );  
  112.     float*  pFloat;  
  113.     double* pDouble;  
  114.     Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);  
  115.     double k        =   kmax/pow(f,nu);  
  116.     double exy      =   sigma_x*sigma_y/2;  
  117.     double scaleReal=   k*k/sigma_x/sigma_y;  
  118.     int    x,y;  
  119.     for( y = ymin; y <= ymax; y++ )  
  120.     {  
  121.         if( ktype == CV_32F )  
  122.         {  
  123.             pFloat = kernel.ptr<float>(ymax-y);  
  124.         }  
  125.         else  
  126.         {  
  127.             pDouble = kernel.ptr<double>(ymax-y);  
  128.         }  
  129.         for( x = xmin; x <= xmax; x++ )  
  130.         {  
  131.             double xr = x*c + y*s;  
  132.             double v = scaleReal*exp(-(x*x+y*y)*scaleReal/2);  
  133.             double temp=cos(k*xr) - exp(-exy);  
  134.             v   =   temp*v;  
  135.             if( ktype == CV_32F )  
  136.             {  
  137.                 pFloat[xmax - x]= (float)v;  
  138.             }  
  139.             else  
  140.             {  
  141.                 pDouble[xmax - x] = v;  
  142.             }  
  143.         }  
  144.     }  
  145.     return kernel;  
  146. }  
  147. Mat GaborFR::getMagnitude(Mat &real,Mat &imag)  
  148. {  
  149.     CV_Assert(real.type()==imag.type());  
  150.     CV_Assert(real.size()==imag.size());  
  151.     int ktype=real.type();  
  152.     int row = real.rows,col = real.cols;  
  153.     int i,j;  
  154.     float*  pFloat,*pFloatR,*pFloatI;  
  155.     double* pDouble,*pDoubleR,*pDoubleI;  
  156.     Mat     kernel(row, col, real.type());  
  157.     for(i=0;i<row;i++)  
  158.     {  
  159.         if( ktype == CV_32FC1 )  
  160.         {  
  161.             pFloat = kernel.ptr<float>(i);  
  162.             pFloatR= real.ptr<float>(i);  
  163.             pFloatI= imag.ptr<float>(i);  
  164.         }  
  165.         else  
  166.         {  
  167.             pDouble = kernel.ptr<double>(i);  
  168.             pDoubleR= real.ptr<double>(i);  
  169.             pDoubleI= imag.ptr<double>(i);  
  170.         }  
  171.         for(j=0;j<col;j++)  
  172.         {  
  173.             if( ktype == CV_32FC1 )  
  174.             {  
  175.                 pFloat[j]= sqrt(pFloatI[j]*pFloatI[j]+pFloatR[j]*pFloatR[j]);  
  176.             }  
  177.             else  
  178.             {  
  179.                 pDouble[j] = sqrt(pDoubleI[j]*pDoubleI[j]+pDoubleR[j]*pDoubleR[j]);  
  180.             }  
  181.         }  
  182.     }  
  183.     return kernel;  
  184. }  
  185. Mat GaborFR::getPhase(Mat &real,Mat &imag)  
  186. {  
  187.     CV_Assert(real.type()==imag.type());  
  188.     CV_Assert(real.size()==imag.size());  
  189.     int ktype=real.type();  
  190.     int row = real.rows,col = real.cols;  
  191.     int i,j;  
  192.     float*  pFloat,*pFloatR,*pFloatI;  
  193.     double* pDouble,*pDoubleR,*pDoubleI;  
  194.     Mat     kernel(row, col, real.type());  
  195.     for(i=0;i<row;i++)  
  196.     {  
  197.         if( ktype == CV_32FC1 )  
  198.         {  
  199.             pFloat = kernel.ptr<float>(i);  
  200.             pFloatR= real.ptr<float>(i);  
  201.             pFloatI= imag.ptr<float>(i);  
  202.         }  
  203.         else  
  204.         {  
  205.             pDouble = kernel.ptr<double>(i);  
  206.             pDoubleR= real.ptr<double>(i);  
  207.             pDoubleI= imag.ptr<double>(i);  
  208.         }  
  209.         for(j=0;j<col;j++)  
  210.         {  
  211.             if( ktype == CV_32FC1 )  
  212.             {  
  213. //              if(pFloatI[j]/(pFloatR[j]+pFloatI[j]) > 0.99)  
  214. //              {  
  215. //                  pFloat[j]=CV_PI/2;  
  216. //              }  
  217. //              else  
  218. //              {  
  219. //                  pFloat[j] = atan(pFloatI[j]/pFloatR[j]);  
  220.                 pFloat[j] = asin(pFloatI[j]/sqrt(pFloatR[j]*pFloatR[j]+pFloatI[j]*pFloatI[j]));  
  221. /*              }*/  
  222. //              pFloat[j] = atan2(pFloatI[j],pFloatR[j]);  
  223.             }//CV_32F  
  224.             else  
  225.             {  
  226.                 if(pDoubleI[j]/(pDoubleR[j]+pDoubleI[j]) > 0.99)  
  227.                 {  
  228.                     pDouble[j]=CV_PI/2;  
  229.                 }  
  230.                 else  
  231.                 {  
  232.                     pDouble[j] = atan(pDoubleI[j]/pDoubleR[j]);  
  233.                 }  
  234. //              pDouble[j]=atan2(pDoubleI[j],pDoubleR[j]);  
  235.             }//CV_64F  
  236.         }  
  237.     }  
  238.     return kernel;  
  239. }  
  240. Mat GaborFR::getFilterRealPart(Mat& src,Mat& real)  
  241. {  
  242.     CV_Assert(real.type()==src.type());  
  243.     Mat dst;  
  244.     Mat kernel;  
  245.     flip(real,kernel,-1);//中心镜面  
  246. //  filter2D(src,dst,CV_32F,kernel,Point(-1,-1),0,BORDER_CONSTANT);  
  247.     filter2D(src,dst,CV_32F,kernel,Point(-1,-1),0,BORDER_REPLICATE);  
  248.     return dst;  
  249. }  
  250. Mat GaborFR::getFilterImagPart(Mat& src,Mat& imag)  
  251. {  
  252.     CV_Assert(imag.type()==src.type());  
  253.     Mat dst;  
  254.     Mat kernel;  
  255.     flip(imag,kernel,-1);//中心镜面  
  256. //  filter2D(src,dst,CV_32F,kernel,Point(-1,-1),0,BORDER_CONSTANT);  
  257.     filter2D(src,dst,CV_32F,kernel,Point(-1,-1),0,BORDER_REPLICATE);  
  258.     return dst;  
  259. }  
  260. void GaborFR::getFilterRealImagPart(Mat& src,Mat& real,Mat& imag,Mat &outReal,Mat &outImag)  
  261. {  
  262.     outReal=getFilterRealPart(src,real);  
  263.     outImag=getFilterImagPart(src,imag);  
  264. }  

main

[cpp] view plaincopyprint?
  1. // Win32TestPure.cpp : 定义控制台应用程序的入口点。  
  2. #include "stdafx.h"  
  3. #include <vector>  
  4. #include <deque>  
  5. #include <iomanip>  
  6. #include <stdexcept>  
  7. #include <string>  
  8. #include <iostream>  
  9. #include <fstream>  
  10. #include <direct.h>//_mkdir()  
  11. #include "opencv2\opencv.hpp"  
  12. #include "GaborFR.h"  
  13. using namespace std;  
  14. using namespace cv;  
  15. int main()  
  16. {   
  17.     //Mat M = getGaborKernel(Size(9,9),2*CV_PI,u*CV_PI/8, 2*CV_PI/pow(2,CV_PI*(v+2)/2),1,0);  
  18.     Mat saveM;  
  19.     //s8-4  
  20.     //s1-5  
  21.     //s1中年男人  
  22.     Mat I=imread("H:\\pic\\s1-5.bmp",-1);  
  23.     normalize(I,I,1,0,CV_MINMAX,CV_32F);  
  24.     Mat showM,showMM;Mat M,MatTemp1,MatTemp2;  
  25.     Mat line;  
  26.     int iSize=50;//如果数值比较大，比如50则接近论文中所述的情况了！估计大小和处理的源图像一样！  
  27.     for(int i=0;i<8;i++)  
  28.     {  
  29.         showM.release();  
  30.         for(int j=0;j<5;j++)  
  31.         {  
  32.             Mat M1= GaborFR::getRealGaborKernel(Size(iSize,iSize),2*CV_PI,i*CV_PI/8+CV_PI/2, j,1);  
  33.             Mat M2 = GaborFR::getImagGaborKernel(Size(iSize,iSize),2*CV_PI,i*CV_PI/8+CV_PI/2, j,1);  
  34.             //加了CV_PI/2才和大部分文献的图形一样，不知道为什么！  
  35.             Mat outR,outI;  
  36.             GaborFR::getFilterRealImagPart(I,M1,M2,outR,outI);  
  37. //          M=GaborFR::getPhase(M1,M2);  
  38. //          M=GaborFR::getMagnitude(M1,M2);  
  39. //          M=GaborFR::getPhase(outR,outI);  
  40. //          M=GaborFR::getMagnitude(outR,outI);  
  41.  //         M=GaborFR::getMagnitude(outR,outI);  
  42. //          MatTemp2=GaborFR::getPhase(outR,outI);  
  43. //          M=outR;  
  44.              M=M1;  
  45.             //      resize(M,M,Size(100,100));  
  46.             normalize(M,M,0,255,CV_MINMAX,CV_8U);  
  47.             showM.push_back(M);  
  48.             line=Mat::ones(4,M.cols,M.type())*255;  
  49.             showM.push_back(line);  
  50.         }  
  51.         showM=showM.t();  
  52.         line=Mat::ones(4,showM.cols,showM.type())*255;  
  53.         showMM.push_back(showM);  
  54.         showMM.push_back(line);  
  55.     }  
  56.     showMM=showMM.t();  
  57. //  bool flag=imwrite("H:\\out.bmp",showMM);  
  58.     imshow("saveMM",showMM);waitKey(0);  
  59.     return 0;  
  60. }//endof   main()  

一下图片可能和程序实际运行结果有点不同，图片只是示意图，代码暂时没问题。需要考虑的是iSize大小问题，首先iSize要用奇数，然后大部分文献iSize都比较大，好像是100左右，但没看到他们描述过卷积核的大小。




posted on 2015-06-22 23:17  moffis  阅读(306)  评论(0)  编辑  收藏  举报

导航