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，又说到。上代码和链接地址！下载地址~

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

//GaborFR.h

#pragma once
#include "opencv2opencv.hpp"
#include <vector>
using namespace std;
using namespace cv;
class GaborFR
{
public:
    GaborFR();
    static Mat  getImagGaborKernel(Size ksize, double sigma, double theta,
                                    double nu,double gamma=1, int ktype= CV_32F);
    static Mat  getRealGaborKernel( Size ksize, double sigma, double theta,
                                    double nu,double gamma=1, int ktype= CV_32F);
    static Mat  getPhase(Mat &real,Mat &imag);
    static Mat  getMagnitude(Mat &real,Mat &imag);
    static void getFilterRealImagPart(Mat& src,Mat& real,Mat& imag,Mat &outReal,Mat &outImag);
    static Mat  getFilterRealPart(Mat& src,Mat& real);
    static Mat  getFilterImagPart(Mat& src,Mat& imag);

    void        Init(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)
{
    double  sigma_x     = sigma;
    double  sigma_y     = sigma/gamma;
    int     nstds       = 3;
    double  kmax        = CV_PI/2;
    double  f           = 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*pi
Mat GaborFR::getRealGaborKernel( Size ksize, double sigma, double theta,
    double nu,double gamma, int ktype)
{
    double  sigma_x     = sigma;
    double  sigma_y     = sigma/gamma;
    int     nstds       = 3;
    double  kmax        = CV_PI/2;
    double  f           = 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 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;
    Mat     kernel(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;
    Mat     kernel(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_32F
            else
            {
                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 "opencv2opencv.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左右，但没看到他们描述过卷积核的大小。