【短道速滑十一】标准的Gabor滤波器及Log_Gabor滤波器的实现、解析、速度优化及其和Halcon中gen_gabor的比较。

/*
 Gabor filters and such. To be greatly extended to have full texture analysis.
 For the formulas and the explanation of the parameters see:
 http://en.wikipedia.org/wiki/Gabor_filter
*/
cv::Mat cv::getGaborKernel( Size ksize, double sigma, double theta,
                            double lambd, double gamma, double psi, int ktype )
{
    double sigma_x = sigma;
    double sigma_y = sigma/gamma;
    int nstds = 3;
    int xmin, xmax, ymin, ymax;
    double c = cos(theta), s = sin(theta);
    if( ksize.width > 0 )
        xmax = ksize.width/2;
    else
        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 );
    Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);
    double scale = 1;
    double ex = -0.5/(sigma_x*sigma_x);
    double ey = -0.5/(sigma_y*sigma_y);
    double cscale = CV_PI*2/lambd;
    for( int y = ymin; y <= ymax; y++ )
        for( int x = xmin; x <= xmax; x++ )
        {
            double xr = x*c + y*s;
            double yr = -x*s + y*c;

            double v = scale*std::exp(ex*xr*xr + ey*yr*yr)*cos(cscale*xr + psi);
            if( ktype == CV_32F )
                kernel.at<float>(ymax - y, xmax - x) = (float)v;
            else
                kernel.at<double>(ymax - y, xmax - x) = v;
        }
    return kernel;
}

static std::vector<cv::Mat> build_filters()
{
    std::vector<cv::Mat> filters;
    const int ksize = 31;
    const double sigma = 4.0;
    const double lambd = 10.0;
    const double gamma = 0.5;
    const double psi = 0;
    // 此处创建16个滤波器， 只有 getGaborkernel 的第三个参数 theta 不同
    for (int i = 0; i < 16; i++)
    {
        double theta = CV_PI * i / 16;
        cv::Mat kernel = cv::getGaborKernel(cv::Size(ksize, ksize), sigma, theta, lambd, gamma, psi, CV_32F);
        kernel /= 1.5 * cv::sum(kernel)[0];
        filters.emplace_back(kernel);
    }
    return filters;
}
cv::Mat process(const cv::Mat& src, std::vector<cv::Mat>& filters)
{
    cv::Mat accum = cv::Mat::zeros(src.size(), src.type());
    for (cv::Mat kernel: filters)
    {
        cv::Mat fimg;
        AutoTimer timer("filter2D");
        cv::filter2D(src, fimg, CV_8UC3, kernel); // 这里是耗时的瓶颈
        AutoTimer timer("getmax");
        accum = cv::max(accum, fimg);
    }
    return accum;
}
int main()
{
    cv::Mat src = cv::imread(image_path);
    std::vector<cv::Mat> filters = build_filters();
    cv::Mat res = process(src, filters);
}

WaveLength = 10;
SigmaR = 0.4;
cols = 500, rows=500;
[x,y] = meshgrid( [-cols/2:(cols/2-1)]/cols,[-rows/2:(rows/2-1)]/rows);
radius = sqrt(x.^2 + y.^2);       
Frequency = 1.0 / WaveLength;                  % 频率等于波长的导数
logGabor = exp((-(log(radius / Frequency)).^2) / (2 * SigmaR * SigmaR));  % log gabor函数的传递函数表达式
imshow(logGabor,[])

Angle = 45 / 180 *3.1415926;
SigmaA = 0.4;
cols = 500, rows=500;
[x,y] = meshgrid( [-cols/2:(cols/2-1)]/cols,[-rows/2:(rows/2-1)]/rows);
theta = atan2(-y,x);              
sintheta = sin(theta);
costheta = cos(theta);
ds = sintheta * cos(Angle) - costheta * sin(Angle);    
dc = costheta * cos(Angle) + sintheta * sin(Angle);     
dtheta = atan2(abs(ds),abs(dc));                           
spread = exp((-dtheta.^2) / (2 * SigmaA * SigmaA));      
imshow(spread,[])

WaveLength = 10;
SigmaR = 0.4;
Angle = 45 / 180 *3.1415926;
SigmaA = 0.4;
cols = 500, rows=500;
[x,y] = meshgrid( [-cols/2:(cols/2-1)]/cols,[-rows/2:(rows/2-1)]/rows);
radius = sqrt(x.^2 + y.^2);  
Frequency = 1.0 / WaveLength;                  % 频率等于波长的导数
logGabor = exp((-(log(radius / Frequency)).^2) / (2 * SigmaR * SigmaR));  % log gabor函数的传递函数表达式
theta = atan2(-y,x);              
sintheta = sin(theta);
costheta = cos(theta);
ds = sintheta * cos(Angle) - costheta * sin(Angle);    
dc = costheta * cos(Angle) + sintheta * sin(Angle);     
dtheta = atan2(abs(ds),abs(dc));                           
spread = exp((-dtheta.^2) / (2 * SigmaA * SigmaA));      
imshow(logGabor .* spread,[])

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　WaveLength = 10，SigmaR = 0.4，Angle = 45 / 180 *3.1415926， SigmaA = 0.4;

read_image (Image, 'fabrik')
get_image_size (Image, Width, Height)
gen_gabor(Filter,10,0.1,50,1.57,'n','dc_center',Width, Height)
fft_generic(Image,ImageFFT,'to_freq',-1,'none','dc_center','complex')
convol_gabor(ImageFFT,Filter,Gabor,Hilbert)
fft_generic(Gabor,GaborInv,'from_freq',1,'none','dc_center','byte')