【Gabor滤波器】提取图像纹理（Python、C++两种实现）

主要参考：（后续链接重复）

• 分析opencv实现的gabor + 不同参数的gabor核（图）：https://blog.csdn.net/lhanchao/article/details/55006663
• 分析wiki百科的gabor实现：https://zhuanlan.zhihu.com/p/33311267
• 原理、opencv+c++实现：https://blog.csdn.net/ibelievesunshine/article/details/105113458
• 本文代码下载：https://gitee.com/zengxy2020/csdn/tree/master/gabor_filter

一、简介

什么是Gabor滤波器？

From：https://blog.csdn.net/ibelievesunshine/article/details/105113458

效果图

二、如何实现

【代码下载】https://gitee.com/zengxy2020/csdn/tree/master/gabor_filter

2.1 Gabor滤波器核函数实部的公式

一般而言，网上gabor滤波器核函数的公式为复数形式，opencv等实际使用仅使用实数部分公式：

参数说明：

• Sigma : σ 表示高斯函数的标准差
• Gamma： γ 表示长宽比，决定这Gabor核函数图像的椭圆率
• Lambda：λ 表示滤波的波长
• Psi ： ψ 相位偏移量，取值范围是-180~180
• Theta： θ 表示Gabor核函数图像的倾斜角度

x，y为gabor滤波器核函数，在opencv实现中，x,y的具体值由gabor核的大小决定，目测为经验值：

    xmax = ksize.width/2;
    ymax = ksize.height/2;
   
   when ksize=9时 ： xmin, ymin,xmax, ymax -4 -4 4 4
    [[-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]] 
    x:
    [[-4 -4 -4 -4 -4 -4 -4 -4 -4]
    [-3 -3 -3 -3 -3 -3 -3 -3 -3]
    [-2 -2 -2 -2 -2 -2 -2 -2 -2]
    [-1 -1 -1 -1 -1 -1 -1 -1 -1]
    [ 0  0  0  0  0  0  0  0  0]
    [ 1  1  1  1  1  1  1  1  1]
    [ 2  2  2  2  2  2  2  2  2]
    [ 3  3  3  3  3  3  3  3  3]
    [ 4  4  4  4  4  4  4  4  4]]

2.2 Python实现：来源于wiki百科

参考：https://zhuanlan.zhihu.com/p/33311267

import cv2
import numpy as np


def gabor_kernel(ksize, sigma, gamma, lamda, alpha, psi):
    '''
    reference
      https://en.wikipedia.org/wiki/Gabor_filter
    '''

    sigma_x = sigma
    sigma_y = sigma / gamma

    ymax = xmax = ksize // 2  # 9//2
    xmin, ymin = -xmax, -ymax
    # print("xmin, ymin,xmin, ymin",xmin, ymin,ymax ,xmax)
    # X(第一个参数，横轴)的每一列一样，  Y（第二个参数，纵轴）的每一行都一样
    (y, x) = np.meshgrid(np.arange(ymin, ymax + 1), np.arange(xmin, xmax + 1))  # 生成网格点坐标矩阵
    # print("y\n",y)
    # print("x\n",x)

    x_alpha = x * np.cos(alpha) + y * np.sin(alpha)
    y_alpha = -x * np.sin(alpha) + y * np.cos(alpha)
    print("x_alpha[0][0]", x_alpha[0][0], y_alpha[0][0])
    exponent = np.exp(-.5 * (x_alpha ** 2 / sigma_x ** 2 +
                             y_alpha ** 2 / sigma_y ** 2))
    # print(exponent[0][0])
    # print(x[0],y[0])
    kernel = exponent * np.cos(2 * np.pi / lamda * x_alpha + psi)
    print(kernel)
    # print(kernel[0][0])
    return kernel


def gabor_filter(gray_img, ksize, sigma, gamma, lamda, psi):
    filters = []
    for alpha in np.arange(0, np.pi, np.pi / 4):
        print("alpha", alpha)
        kern = gabor_kernel(ksize=ksize, sigma=sigma, gamma=gamma,
                            lamda=lamda, alpha=alpha, psi=psi)
        filters.append(kern)

    gabor_img = np.zeros(gray_img.shape, dtype=np.uint8)

    i = 0
    for kern in filters:
        fimg = cv2.filter2D(gray_img, ddepth=cv2.CV_8U, kernel=kern)
        gabor_img = cv2.max(gabor_img, fimg)
        cv2.imwrite("2." + str(i) + "gabor.jpg", gabor_img)
        i += 1
    p = 1.25
    gabor_img = (gabor_img - np.min(gabor_img, axis=None))**p
    _max = np.max(gabor_img, axis=None)
    gabor_img = gabor_img / _max
    print(gabor_img)
    gabor_img = gabor_img * 255
    return gabor_img.astype(dtype=np.uint8)


def main():
    src = cv2.imread("1.jpg")
    src_gray = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)
    gabor_img = gabor_filter(src_gray,
                             ksize=9,
                             sigma=1,
                             gamma=0.5,
                             lamda=5,
                             psi=-np.pi / 2)
    cv2.imwrite("gabor.jpg", gabor_img)


main()

2.3 C++实现（wiki和opencv源码）

// #include "opencv2/core.hpp"
#include <opencv2/opencv.hpp>
#include <opencv2/ximgproc.hpp>
#include <iostream>
#include <time.h>

using namespace cv;
using namespace std;

void gabor_filter(bool if_opencv_kernal,Mat gray_img, Mat& gabor_img, Mat gabor_tmp,int k=9, float sigma=1.0,float gamma=0.5,float lambda=5.0,float psi=-CV_PI / 2);
Mat gabor_kernal_wiki(Size ksize,double theta,double sigma=1.0,
                            double lambd=5.0, double gamma=0.5, double psi=-CV_PI / 2, int ktype=CV_32F);

int main()
{   
   
    clock_t total_start=clock();
    Mat src=imread("1.jpg");
    Mat src_gray,src_gray_gabor,gabor_tmp;
    cvtColor(src,src_gray,COLOR_BGR2GRAY);

    bool if_opencv_kernal=true;//false;
    gabor_filter(if_opencv_kernal,src_gray, src_gray_gabor, gabor_tmp);
    
    char save_name[100];
    sprintf(save_name, "gabor_filter_img.jpg"); 
    imwrite(save_name,src_gray_gabor);
    
    clock_t total_end=clock();
    cout<<"Total Time:"<<(double)(total_end - total_start)/ CLOCKS_PER_SEC*1000 <<endl;
    return 0;

}


/*Gabor 滤波处理相关
Sigma : σ 表示高斯函数的标准差
Gamma： γ 表示长宽比，决定这Gabor核函数图像的椭圆率
Lambda：λ 表示滤波的波长
Psi  ： ψ 相位偏移量，取值范围是-180~180
Theta： θ  表示Gabor核函数图像的倾斜角度 
*/
void  gabor_filter(bool if_opencv_kernal,Mat gray_img, Mat& gabor_img, Mat gabor_tmp,int k, float sigma,float gamma,float lambda, float psi) 
{   
    int ddepth=CV_8U;  
    //表示Gabor核函数图像的倾斜角度 // 0  45 90 135度 的核
    double theta[4]={
        0.0,
        CV_PI/4,
        CV_PI/2,
        CV_PI/4*3,
    }; 

    Size ksize=Size(k,k); //卷积核的大小
    Mat gabor_kernel;
    char save_name[100];
    gabor_img=Mat::zeros(gray_img.size(),CV_8UC1);
    // 分别取 0  45 90 135 四个方向gabor核进行计算
    for  (int i=0;i<4;i++)
    {       
        cout<<"theta θ= "<<i*45<<"="<<theta[i]<<endl;
        if(if_opencv_kernal) 
            gabor_kernel= getGaborKernel(ksize, sigma, theta[i], lambda, gamma, psi, CV_32F);
        else 
            gabor_kernel=gabor_kernal_wiki(ksize,theta[i]);
        // cout<<gabor_kernel<<endl;
        filter2D(gray_img,gabor_tmp,ddepth,gabor_kernel);

        /*合成各个角度的图*/
        //cv::max(src1,src2,dst);
        cv::max(gabor_tmp,gabor_img,gabor_img);
        sprintf(save_name, "im_resize_gray_gabor_%d.jpg",i*45); 
        imwrite(save_name,gabor_img);
    }



}


Mat gabor_kernal_wiki(Size ksize, double theta,double sigma, 
                            double lambd, double gamma, double psi, int ktype)
{
    double sigma_x = sigma;   //  σ 表示高斯函数的标准差
    double sigma_y = sigma/gamma; // γ 表示长宽比，决定这Gabor核函数图像的椭圆率

    int xmin, xmax, ymin, ymax;
    // double c = cos(theta), s = sin(theta);

    xmax = ksize.width/2;
    ymax = ksize.height/2;
    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);
    /*
    y
    xmin, ymin,xmax, ymax -4 -4 4 4
    [[-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]] 
    x:
    [[-4 -4 -4 -4 -4 -4 -4 -4 -4]
    [-3 -3 -3 -3 -3 -3 -3 -3 -3]
    [-2 -2 -2 -2 -2 -2 -2 -2 -2]
    [-1 -1 -1 -1 -1 -1 -1 -1 -1]
    [ 0  0  0  0  0  0  0  0  0]
    [ 1  1  1  1  1  1  1  1  1]
    [ 2  2  2  2  2  2  2  2  2]
    [ 3  3  3  3  3  3  3  3  3]
    [ 4  4  4  4  4  4  4  4  4]]

    */
   for( int x = xmin; x <= xmax; x++ )
   {
    // cout<<x<<endl;
    for( int y = ymin; y <= ymax; y++ )
        {   
            // cout<<y<<" ";
            double x_alpha=x*cos(theta)+y*sin(theta);
            double y_alpha=-x*sin(theta)+y*cos(theta);
            // double xr = x*c + y*s;
            // double yr = -x*s + y*c;
            // double exponent=std::exp(-0.5*(x_alpha*x_alpha/sigma_x*sigma_x+
            double exponent=exp(-0.5*(x_alpha*x_alpha/pow(sigma_x,2)+
                                           y_alpha*y_alpha/pow(sigma_y,2)
                                     )
                                );
            double v = exponent*cos(2*CV_PI/lambd*x_alpha + psi);
            if( ktype == CV_32F )
                {
                // printf("%d %d %d %d %d %d ",ymax,y, x,xmax, x+xmax,y+ymax);
                // printf("%d %d %d %d ", ymax - y, xmax - x,x+xmax,y+ymax);   
                kernel.at<float>(y+ymax,x+xmax) = (float)v;  // 从
                // kernel.at<float>(ymax - y, xmax - x) = (float)v;  // opencv的赋值
                }
            else
                kernel.at<double>(y+ymax,x+xmax) = v;
        }
    // cout<<endl;
   }

    cout<<kernel<<endl;
    return kernel;
}