Gabor滤波器

发布日期: 2020-11-02

Gabor滤波器

Gabor Filter

在图像处理中，Gabor函数是一个用于边缘提取的线性滤波器。Gabor滤波器的频率和方向表达同人类视觉系统类似。研究发现，Gabor滤波器十分适合纹理表达和分离。在空间域中，一个二维Gabor滤波器是一个由正弦平面波调制的高斯核函数。还有，生物学实验发现，Gabor滤波器可以很好地近似单细胞的感受野函数（光强刺激下的传递函数），什么视皮层内的超柱，bla…bla，总之是这方面仿生的数学模型。另外，网上有一种说法，gabor分为实部和虚部，用实部进行滤波后图像会平滑；虚部滤波后用来检测边缘。

Gabor滤波器的脉冲响应，可以定义为一个正弦波（对于二维Gabor滤波器是正弦平面波）乘以高斯函数。由于乘法卷积性质，Gabor滤波器的脉冲响应的傅立叶变换是其调和函数的傅立叶变换和高斯函数傅立叶变换的卷积。该滤波器由实部和虚部组成，二者相互正交。一组不同频率不同方向的Gabor函数数组对于图像特征提取非常有用。

一维Gabor滤波器

Gabor滤波器是处理一维信号（比如音频）最佳的带通滤波器¹，一个复杂的Gabor滤波器是一个高斯核函数²乘以一个复杂的sin函数：

$$
Gabor(t)=ke^{i\theta}w(at)s(t) \tag{1}
$$
其中：
$$
\begin{cases}
w(t)=e^{- \pi t^2} \\
s(t)=e^{i(2\pi f_0 t)}
\end{cases} \tag{2}
$$

其中，$f_0$是复数波$s(t)$的频率。
将复数波$s(t)=e^{i(2\pi f_0 t)}$带入（1）式中，得到：
$$
\begin{aligned}
Gabor(t) & = k \omega(at) e^{i(2\pi f_0 t + \theta)} \\ & = k \omega(at) \left[ \cos(2\pi f_0 t+\theta) + i\sin(2\pi f_0 t+\theta) \right]
\end{aligned} \tag{3}
$$
上面最后一步得到了 Gabor 核的复数表示，我们就可以按实部和虚部将其拆分为实核和虚核，在很多应用中，我们只需要应用 Gabor核的实数部分即可：
$$
\begin{cases} Gabor_{real}(t) = \omega(at)\cos(2\pi f_0 t + \theta) \\ Gabor_{imag}(t) = \omega(at)\sin(2\pi f_0 t + \theta) \end{cases} \tag{4}
$$

二维Gabor滤波器

二维Gabor滤波器，也叫空间Gabor滤波器（The Spatial (2-D) Gabor Filter）。类似一维 Gabor 核，我们将二维高斯函数与二维复数波相乘，就得到了二维的Gabor核：
$$
\begin{aligned}
Gabor(x_0, y_0, \theta, \sigma_x, \sigma_y, u_0, v_0) & = s(x,y) \omega_r(x,y) \\
& = K \exp\left(-\pi \left( (x-x_0)_r^2/\sigma_x^2 + (y-y_0)_r^2/\sigma_y^2 \right) \right) \exp\left(i 2\pi (u_0 x + v_0 y) \right) \
\end{aligned} \tag{5}
$$
它的各个参数含义如下：

$(x_0, y_0)$: 高斯核的中心点
$\theta$: 高斯核的旋转方向（顺时针）
$(\sigma_x, \sigma_y)$: 高斯核两个方向上的尺度
$(u_0, v_0)$: 频域坐标
$K$: 高斯核的幅度（magnitude）的比例

Gabor之Python实现

Gabor变换

其实，skimage包中就有封装好的gabor变换函数，可以直接调用。因此，这里介绍一种简单偷懒的方式。从skimage中导入filters（from skimage import filters)后,便可调用gabor函数了，其函数说明如下：
skimage.filters.gabor(image, frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None, n_stds=3, offset=0, mode=‘reflect’, cval=0)

函数返回：
- Gabor变换后的实部和虚部real，imag；
- 大小与输入图像尺寸相同。
函数参数：
- 图像（image）：二维数组输入图像（灰度图像）；
- 频率（frequency）：浮点谐波函数的空间频率，控制尺度。$k_v=k_{max}/f^v$,$k_{max}=\pi/2$为最大的采用频率。在图像特征提取领域通常取5个不同尺度,v = [0,1,2,3,4]。
- 方向（theta）：float，可选的弧度方向，如果为0，则谐波处于x方向。$\varphi_u=\pi u/8$,在图像特征提取领域通常取8个不同方向u=[0,1,2,3,4,5,6,7]。
- 带宽（bandwidth)：浮点，可选由过滤器捕获的带宽。对于固定带宽，sigma_x和sigma_y将随着频率的增加而降低，如果用户设置了sigma_x和sigma_y，则该值将被忽略；
- 标准偏差(sigma_x，sigma_y)：float，x和y方向上的可选标准偏差；
- 内核的线性大小(n_stds)：标量，可选内核的线性大小为n_stds（默认为3）标准偏差；
- 偏移量(offset)：浮点数，可选项以弧度表示的谐波函数的相位偏移；
- 模式(mode）：{‘constant’，‘near’，‘reflect’，‘mirror’，‘wrap’}，可选用于将图像与内核进行卷积的模式，传递给ndi.convolve；
- cval：标量，可选值如果卷积模式为’不变’，该参数被传递给ndi.convolve。

源码：

import matplotlib.pyplot as plt
from skimage import filters,io,color
import numpy as np

filename='D:/lena.jpg'
img = io.imread(filename)#读取图像
img_gray = color.rgb2gray(img)#RGB转灰度
frequency=0.6
#调用gabor函数
real, imag = filters.gabor(img_gray, frequency=0.6,theta=45,n_stds=5)
#取模图像
img_mod=np.sqrt(real.astype(float)**2+imag.astype(float)**2)
#图像显示
plt.figure()
plt.subplot(2,2,1)
plt.imshow(img_gray,cmap='gray')
plt.subplot(2,2,2)
plt.imshow(img_mod,cmap='gray')
plt.subplot(2,2,3)
plt.imshow(real,cmap='gray')
plt.subplot(2,2,4)
plt.imshow(imag,cmap='gray')
plt.show()

注：

通过theta参数调不同方向，这里输入是弧度，不是角度。
通过frequency参数调不同尺度变化。

Gabor卷积核

利用skimage生成gobor卷积核：gabor_kernel
skimage.filters.gabor_kernel(frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None, n_stds=3, offset=0)

函数返回：返回2D Gabor滤波器内核，包含实部与虚部。
参数与skimage.filters.gabor（）函数相同。

源码：

import matplotlib.pyplot as plt
from skimage import filters
gk = filters.gabor_kernel(frequency=0.1,theta=np.pi*30/180.0,n_stds=5)
mod=np.sqrt(gk.real.astype(float) ** 2 + gk.imag.astype(float) ** 2)
plt.figure()
plt.subplot(1,3,1)
plt.imshow(gk.real*255,cmap='gray')
plt.subplot(1,3,2)
plt.imshow(gk.imag*255,cmap='gray')
plt.subplot(1,3,3)
plt.imshow(mod*255,cmap='gray')
plt.show()

注：

参数n_stds=3并不是说滤波核大小为（3,3)，滤波核大小由参数frequency、theta、n_stds三个参数共同决定。
gabor_kernel函数返回是带有实部和虚部的复矩阵。
theta是弧度，不是角度

Gabor特征

Gabor滤波器可以在频域上不同尺度、不同方向上提取相关的特征。另外，Gabor函数与人眼的作用相仿，所以经常用作纹理识别上，并取得了较好的效果。Gabor特征提取一般包括取模、特征降维、特征归一化和特征选取等操作：

取模：skimage.filters.gabor（）函数返回的是图像变换后的实部和虚部，在图像识别领域一般使用其模作为图像特征。
$$mod=\sqrt{real^2+imag^2} \tag{6}$$
特征降维：skimage.filters.gabor（）函数返回值大小和原图像一样；另外，图像识别领域一般使用8个方向和5个尺度的gabor滤波器，得到40幅变换图像；导致提取的图像特征维度（40幅图像特征串联）比较高。因此，需要对取模图像进行降维处理。最简单的方法就是对每个模图像进行下采样，如使用Opencv中的resize函数cv2.resize(img,(0,0),fx=1/4,fy=1/4)。
特征归一化：在机器学习领域中，不同评价指标（即特征向量中的不同特征就是所述的不同评价指标）往往具有不同的量纲和量纲单位，这样的情况会影响到数据分析的结果，为了消除指标之间的量纲影响，需要进行数据标准化处理，以解决数据指标之间的可比性。原始数据经过数据标准化处理后，各指标处于同一数量级，适合进行综合对比评价。其中，最典型的就是数据的归一化处理。简而言之，归一化的目的就是使得预处理的数据被限定在一定的范围内。一种Z-score标准化方法：
$$x=\frac{x-\mu}{\sigma} \tag{7}$$
- 其中μ为所有样本数据的均值，σ为所有样本数据的标准差。
- 本方法要求原始数据的分布可以近似为高斯分布，否则归一化的效果会变得很糟糕；
- 应用场景：在分类、聚类算法中，需要使用距离来度量相似性的时候、或者使用PCA技术进行降维的时候，Z-score standardization表现更好。
特征选取：经过以上三步后，可以将各模图像特征串联起来作为输入图像的特征向量。但此时得到的特征向量存在大量冗余，最好再进行一次特征降维（特征选取）操作。可以采用PCA或Fisher等线性子空间分析方法。

源码：

import cv2
import numpy as np
from skimage import filters
filename='D:/lena.jpg'
img = cv2.imread(filename)#读图像
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)#转灰度
frequency=0.6
#gabor变换
real, imag = filters.gabor(img_gray, frequency=0.6,theta=45,n_stds=5)
#取模
img_mod=np.sqrt(real.astype(float)**2+imag.astype(float)**2)
#图像缩放（下采样）
newimg = cv2.resize(img_mod,(0,0),fx=1/4,fy=1/4,interpolation=cv2.INTER_AREA)
tempfea = newimg.flatten()#矩阵展平
tmean = np.mean(tempfea)#求均值
tstd = np.std(tempfea)#求方差
newfea = (tempfea - tmean)/tstd#数值归一化
print(newfea)