Gabor的OpenCV代码

    唯一持续维护地址：http://52coding.com/opencv-gabor

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

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

//GaborFR.h

#pragma once#include "opencv2\opencv.hpp"#include <vector>using namespace std;using namespace cv;class GaborFR{public:GaborFR();static MatgetImagGaborKernel(Size ksize, double sigma, double theta, double nu,double gamma=1, int ktype= CV_32F);static MatgetRealGaborKernel( Size ksize, double sigma, double theta, double nu,double gamma=1, int ktype= CV_32F);static MatgetPhase(Mat &real,Mat &imag);static MatgetMagnitude(Mat &real,Mat &imag);static void getFilterRealImagPart(Mat& src,Mat& real,Mat& imag,Mat &outReal,Mat &outImag);static MatgetFilterRealPart(Mat& src,Mat& real);static MatgetFilterImagPart(Mat& src,Mat& imag);voidInit(Size ksize=Size(19,19), double sigma=2*CV_PI,double gamma=1, int ktype=CV_32FC1);private:vector<Mat> gaborRealKernels;vector<Mat> gaborImagKernels;bool isInited;};

#include "stdafx.h"#include "GaborFR.h"GaborFR::GaborFR(){isInited = false;}void GaborFR::Init(Size ksize, double sigma,double gamma, int ktype){gaborRealKernels.clear();gaborImagKernels.clear();double mu[8]={0,1,2,3,4,5,6,7};double nu[5]={0,1,2,3,4};int i,j;for(i=0;i<5;i++){for(j=0;j<8;j++){gaborRealKernels.push_back(getRealGaborKernel(ksize,sigma,mu[j]*CV_PI/8,nu[i],gamma,ktype));gaborImagKernels.push_back(getImagGaborKernel(ksize,sigma,mu[j]*CV_PI/8,nu[i],gamma,ktype));}}isInited = true;}Mat GaborFR::getImagGaborKernel(Size ksize, double sigma, double theta, double nu,double gamma, int ktype){doublesigma_x= sigma;doublesigma_y= sigma/gamma;intnstds= 3;doublekmax= CV_PI/2;doublef= cv::sqrt(2.0);int xmin, xmax, ymin, ymax;double c = cos(theta), s = sin(theta);if( ksize.width > 0 ){xmax = ksize.width/2;}else//这个和matlab中的结果一样，默认都是19 !{xmax = cvRound(std::max(fabs(nstds*sigma_x*c), fabs(nstds*sigma_y*s)));}if( ksize.height > 0 ){ymax = ksize.height/2;}else{ymax = cvRound(std::max(fabs(nstds*sigma_x*s), fabs(nstds*sigma_y*c)));}xmin = -xmax;ymin = -ymax;CV_Assert( ktype == CV_32F || ktype == CV_64F );float*pFloat;double*pDouble;Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);double k=kmax/pow(f,nu);double scaleReal=k*k/sigma_x/sigma_y;for( int y = ymin; y <= ymax; y++ ){if( ktype == CV_32F ){pFloat = kernel.ptr<float>(ymax-y);}else{pDouble = kernel.ptr<double>(ymax-y);}for( int x = xmin; x <= xmax; x++ ){double xr = x*c + y*s;double v = scaleReal*exp(-(x*x+y*y)*scaleReal/2);double temp=sin(k*xr);v=  temp*v;if( ktype == CV_32F ){pFloat[xmax - x]= (float)v;}else{pDouble[xmax - x] = v;}}}return kernel;}//sigma一般为2*piMat GaborFR::getRealGaborKernel( Size ksize, double sigma, double theta, double nu,double gamma, int ktype){doublesigma_x= sigma;doublesigma_y= sigma/gamma;intnstds= 3;doublekmax= CV_PI/2;doublef= cv::sqrt(2.0);int xmin, xmax, ymin, ymax;double c = cos(theta), s = sin(theta);if( ksize.width > 0 ){xmax = ksize.width/2;}else//这个和matlab中的结果一样，默认都是19 !{xmax = cvRound(std::max(fabs(nstds*sigma_x*c), fabs(nstds*sigma_y*s)));}if( ksize.height > 0 )ymax = ksize.height/2;elseymax = cvRound(std::max(fabs(nstds*sigma_x*s), fabs(nstds*sigma_y*c)));xmin = -xmax;ymin = -ymax;CV_Assert( ktype == CV_32F || ktype == CV_64F );float*pFloat;double*pDouble;Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);double k=kmax/pow(f,nu);double exy=sigma_x*sigma_y/2;double scaleReal=k*k/sigma_x/sigma_y;int   x,y;for( y = ymin; y <= ymax; y++ ){if( ktype == CV_32F ){pFloat = kernel.ptr<float>(ymax-y);}else{pDouble = kernel.ptr<double>(ymax-y);}for( x = xmin; x <= xmax; x++ ){double xr = x*c + y*s;double v = scaleReal*exp(-(x*x+y*y)*scaleReal/2);double temp=cos(k*xr) - exp(-exy);v=temp*v;if( ktype == CV_32F ){pFloat[xmax - x]= (float)v;}else{pDouble[xmax - x] = v;}}}return kernel;}Mat GaborFR::getMagnitude(Mat &real,Mat &imag){CV_Assert(real.type()==imag.type());CV_Assert(real.size()==imag.size());int ktype=real.type();int row = real.rows,col = real.cols;int i,j;float*pFloat,*pFloatR,*pFloatI;double*pDouble,*pDoubleR,*pDoubleI;Matkernel(row, col, real.type());for(i=0;i<row;i++){if( ktype == CV_32FC1 ){pFloat = kernel.ptr<float>(i);pFloatR= real.ptr<float>(i);pFloatI= imag.ptr<float>(i);}else{pDouble = kernel.ptr<double>(i);pDoubleR= real.ptr<double>(i);pDoubleI= imag.ptr<double>(i);}for(j=0;j<col;j++){if( ktype == CV_32FC1 ){pFloat[j]= sqrt(pFloatI[j]*pFloatI[j]+pFloatR[j]*pFloatR[j]);}else{pDouble[j] = sqrt(pDoubleI[j]*pDoubleI[j]+pDoubleR[j]*pDoubleR[j]);}}}return kernel;}Mat GaborFR::getPhase(Mat &real,Mat &imag){CV_Assert(real.type()==imag.type());CV_Assert(real.size()==imag.size());int ktype=real.type();int row = real.rows,col = real.cols;int i,j;float*pFloat,*pFloatR,*pFloatI;double*pDouble,*pDoubleR,*pDoubleI;Matkernel(row, col, real.type());for(i=0;i<row;i++){if( ktype == CV_32FC1 ){pFloat = kernel.ptr<float>(i);pFloatR= real.ptr<float>(i);pFloatI= imag.ptr<float>(i);}else{pDouble = kernel.ptr<double>(i);pDoubleR= real.ptr<double>(i);pDoubleI= imag.ptr<double>(i);}for(j=0;j<col;j++){if( ktype == CV_32FC1 ){// if(pFloatI[j]/(pFloatR[j]+pFloatI[j]) > 0.99)// {// pFloat[j]=CV_PI/2;// }// else// {//pFloat[j] = atan(pFloatI[j]/pFloatR[j]);pFloat[j] = asin(pFloatI[j]/sqrt(pFloatR[j]*pFloatR[j]+pFloatI[j]*pFloatI[j]));/*}*///pFloat[j] = atan2(pFloatI[j],pFloatR[j]);}//CV_32Felse{if(pDoubleI[j]/(pDoubleR[j]+pDoubleI[j]) > 0.99){pDouble[j]=CV_PI/2;}else{pDouble[j] = atan(pDoubleI[j]/pDoubleR[j]);}//pDouble[j]=atan2(pDoubleI[j],pDoubleR[j]);}//CV_64F}}return kernel;}Mat GaborFR::getFilterRealPart(Mat& src,Mat& real){CV_Assert(real.type()==src.type());Mat dst;Mat kernel;flip(real,kernel,-1);//中心镜面//filter2D(src,dst,CV_32F,kernel,Point(-1,-1),0,BORDER_CONSTANT);filter2D(src,dst,CV_32F,kernel,Point(-1,-1),0,BORDER_REPLICATE);return dst;}Mat GaborFR::getFilterImagPart(Mat& src,Mat& imag){CV_Assert(imag.type()==src.type());Mat dst;Mat kernel;flip(imag,kernel,-1);//中心镜面//filter2D(src,dst,CV_32F,kernel,Point(-1,-1),0,BORDER_CONSTANT);filter2D(src,dst,CV_32F,kernel,Point(-1,-1),0,BORDER_REPLICATE);return dst;}void GaborFR::getFilterRealImagPart(Mat& src,Mat& real,Mat& imag,Mat &outReal,Mat &outImag){outReal=getFilterRealPart(src,real);outImag=getFilterImagPart(src,imag);}

main

// Win32TestPure.cpp : 定义控制台应用程序的入口点。#include "stdafx.h"#include <vector>#include <deque>#include <iomanip>#include <stdexcept>#include <string>#include <iostream>#include <fstream>#include <direct.h>//_mkdir()#include "opencv2\opencv.hpp"#include "GaborFR.h"using namespace std;using namespace cv;int main(){ //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);Mat saveM;//s8-4//s1-5//s1中年男人Mat I=imread("H:\\pic\\s1-5.bmp",-1);normalize(I,I,1,0,CV_MINMAX,CV_32F);Mat showM,showMM;Mat M,MatTemp1,MatTemp2;Mat line;int iSize=50;//如果数值比较大，比如50则接近论文中所述的情况了！估计大小和处理的源图像一样！for(int i=0;i<8;i++){showM.release();for(int j=0;j<5;j++){Mat M1= GaborFR::getRealGaborKernel(Size(iSize,iSize),2*CV_PI,i*CV_PI/8+CV_PI/2, j,1);Mat M2 = GaborFR::getImagGaborKernel(Size(iSize,iSize),2*CV_PI,i*CV_PI/8+CV_PI/2, j,1);//加了CV_PI/2才和大部分文献的图形一样，不知道为什么！Mat outR,outI;GaborFR::getFilterRealImagPart(I,M1,M2,outR,outI);//M=GaborFR::getPhase(M1,M2);//M=GaborFR::getMagnitude(M1,M2);//M=GaborFR::getPhase(outR,outI);//M=GaborFR::getMagnitude(outR,outI); //M=GaborFR::getMagnitude(outR,outI);// MatTemp2=GaborFR::getPhase(outR,outI);// M=outR; M=M1;// resize(M,M,Size(100,100));normalize(M,M,0,255,CV_MINMAX,CV_8U);showM.push_back(M);line=Mat::ones(4,M.cols,M.type())*255;showM.push_back(line);}showM=showM.t();line=Mat::ones(4,showM.cols,showM.type())*255;showMM.push_back(showM);showMM.push_back(line);}showMM=showMM.t();//bool flag=imwrite("H:\\out.bmp",showMM);imshow("saveMM",showMM);waitKey(0);return 0;}//endof   main()

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