魔幻滤镜Gabor Filter 原理与实现（python & C++）

Posted 2023-02-01 冰不语

Gabor Filter

以下定义翻译引用自维基百科：

在图像处理领域，以Dennis Gabor命名的Gabor滤波器，是一种用于纹理分析的线性滤波器，即它主要分析的是，图像在某一特定区域的特定方向上是否有特定的频率内容。当代许多视觉科学家认为，Gabor滤波器的频率和方向的表达与人类的视觉系统很相似，尽管并没有实验性证据和函数原理能证明这一观点。研究发现，Gabor滤波器特别适合于纹理表示和辨别。在空间域中，2D Gabor滤波器是由正弦平面波调制的高斯核函数。

其脉冲响应定义为定义为一个正弦波（对于二维Gabor滤波器是平面波）乘以高斯函数。因为乘法卷积性质（卷积定理），由于乘法卷积性质，Gabor滤波器的脉冲响应的傅立叶变换是调和函数的傅立叶变换和高斯函数傅立叶变换的卷积。该滤波器具有表示正交方向的实部和虚部。两个分量可以构成复数或单独使用。

• 复数形式
  

• 实部
  

• 虚部
  

其中

如果你了解傅立叶变换的基础，那么你应该知道，图像可以视为各个方向的一系列不同频率的正弦波叠加而成。变换中的“像素”表示波的强度，“像素”的位置表示波的频率和方向。 在实践中，人们只想选择特定频率和特定方向的特定波。Gabor变换是许多所谓的带通滤波器之一，它允许你“切割”傅立叶变换，只隔离特定的信息。另一个重要信息是，每个傅立叶“像素”是一个复数值（包含实部和虚部）。

二维Gabor核函数由一个高斯包络函数和一个余弦函数相乘得出，其中θ,ϕ,γ,λ,σ为参数。

A complex Gabor filter is defined as the product of a Gaussian kernel times a complex sinusoid.

Gaussian envelope mentiond in this method determines the width of the Gaussian kernel.

Paper— Fingerprint Image Enhancement: Algorithm and Performance Evaluation

需要注意的是，（x,y）是原始坐标，（x’,y’）是旋转后的坐标。在（x,y）坐标系中，Gabor只能沿x轴或y轴拉伸，但不能对角拉伸。但是上述最后两个公式允许Gabor能够在（x’,y’）坐标系中沿任何方向被拉伸（由θ定义）。

Python实现

维基百科给出了Gabor函数的Python实现。放在下面供参考。

import numpy as np

def gabor_fn(sigma, theta, Lambda, psi, gamma):
    sigma_x = sigma
    sigma_y = float(sigma) / gamma

    # ------这部分内容是为了确定卷积核的大小------
    # Bounding box
    nstds = 3 # Number of standard deviation sigma
    xmax = max(abs(nstds * sigma_x * np.cos(theta)), abs(nstds * sigma_y * np.sin(theta)))
    xmax = np.ceil(max(1, xmax))
    ymax = max(abs(nstds * sigma_x * np.sin(theta)), abs(nstds * sigma_y * np.cos(theta)))
    ymax = np.ceil(max(1, ymax))
    xmin = -xmax
    ymin = -ymax
    (y, x) = np.meshgrid(np.arange(ymin, ymax + 1), np.arange(xmin, xmax + 1))
    # ------这部分内容是为了确定卷积核的大小------

    # Rotation 
    x_theta = x * np.cos(theta) + y * np.sin(theta)
    y_theta = -x * np.sin(theta) + y * np.cos(theta)

    # ------这部分正是上面的公式（1）------
    gb = np.exp(-.5 * (x_theta ** 2 / sigma_x ** 2 + y_theta ** 2 / sigma_y ** 2)) * np.cos(2 * np.pi / Lambda * x_theta + psi)
    return gb

Gabor函数的OpenCV实现

跟其他滤波器（高斯，Sobel）等一样，需要先确定一个卷积核，然后用次卷积核通过滑窗法对输入图像进行卷积运算，得到输出图像。Gabor函数的OpenCV实现getGaborKernel函数，正是为了得到这个卷积核的。为在OpenCV中的getGaborKernel函数里需要传入的参数除了上述5个外，还需要传入卷积核的大小。

cv::Mat getGaborKernel(Size ksize, double sigma, double theta, double lambd, double gamma, double psi=CV_PI*0.5, int ktype=CV_64F );

cv2.getGaborKernel(ksize, sigma, theta, lambda, gamma, psi, ktype)

各参数含义

ksize：返回的滤波器的大小，即Gabor核的大小。如果ksize=(a，b)，那么Gabor核的大小即为axb像素。与许多其他的卷积核一样，这个尺寸一般取奇数而且为正方形。

sigma：σ是Gabor过滤器中实用的高斯函数的标准差。

theta：θ是Gabor函数的平行条纹的法线方向。有效值为从0~360度的实数。

lambda：λ是上述方程中正弦因子的波长。有效值应大于2，以像素为单位。

gamma：γ是空间宽高比，表示Gabor滤波器的椭圆度

psi：ψ是相位偏移。有效值为-180度~180度，0度和180度对应的方程与原点对称，-90度和90度的方程分别于原点成中心对称，可两项直角坐标系中的余弦函数。

ktype：像素类型。

OpenCV的Gabor源码

源码在opencv\\sources\\modules\\imgproc\\src\\gabor.cpp。仔细看上下两种实现的话，可以明显看出OpenCV的实现参考了维基百科的Python实现。只不过在加了一个确定卷积和大小的参数而已。而且在这个参数为（0,0）的时候，采用的正是维基百科中的方法。

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);

    // ------这部分内容是为了确定卷积核的大小------
    // ------可以看到，如果ksize不为0，则取ksize的一半------
    if( ksize.width > 0 )
        xmax = ksize.width/2;
    // ------如果ksize为0，则按照维基百科的方式来取值------
    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;

            // ------此处正是上面的公式（1）------
            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;

参考

1. https://en.wikipedia.org/wiki/Gabor_filter
2. https://dsp.stackexchange.com/questions/14714/understanding-the-gabor-filter-function
3. https://math.stackexchange.com/questions/259877/value-of-x-y-in-computing-gabor-filter-function
4. https://cvtuts.wordpress.com/2014/04/27/gabor-filters-a-practical-overview/
5. http://blog.csdn.net/xufangpang/article/details/6426244
6. https://math.stackexchange.com/questions/177977/how-to-intuitively-interpret-gabor-lambda-param

---

公众号CVPy，分享不仅限于OpenCV和Python的有意思的内容。