Python OpenCV实现Log Gabor滤波器（由LGHD描述符扩展）

引言

笔者在研究红外图像与可见光图像配准时接触到了很多描述符，其中关于LGHD描述符的Log Gabor滤波器很有意思，与大家分享

LGHD（Log-Gabor Histogram Descriptor）

描述符的思想是基于高频分量分布的描述符对于不同的非线性强度变化具有鲁棒性。可以理解为虽然非线性强度差异可以影响图像，但是场景中所包含对象的形状的整体外观仍趋于保持恒定。Log-Gabor 滤波器是根据Gabor滤波器衍生而来，Gabor变换能达到时频局部化的目的，它能够在整体上提供信号的全部信息又能提供在任意局部时间内信号变化剧烈程度的信息。

Log-Gabor滤波器

Log-Gabor 滤波器的表达是从Gabor衍生而来的，相比于Gabor滤波器，它的传递函数是对数频率尺度下的高斯函数，始终没有直流分量，在图像处理时可以不受亮度的影响。
并经研究Log-Gabor函数更符合哺乳动物的视觉观察系统，有更优的纹理提取、目标检测效果。
Log-Gabor的传递函数为：
G(ω)=e−log⁡(ω/ω0)2/(2(log⁡(k/ω0)2))G(\omega)=e^{-\log(\omega/\omega_0)^2/(2(\log(k/\omega_0)^2))} G(ω)=e−log(ω/ω0)2/(2(log(k/ω0)2))
Log-Gabor函数并不能在空间域中得到表达式，滤波器的构造须在频域中。对图像来说，在空域的卷积等价于在频域的乘积，所以我们直接创建与图像大小相同的Log-Gabor滤波器，在频域中完成图像处理，之后再回到空域中进行统计。一个二维的Log-Gabor滤波器可以分解为径向滤波器和角度滤波器两部分，对应极坐标公式为：

完整的Log-Gabor滤波器由这两部分相乘得到：

其中r是经向坐标，θ\thetaθ为角度坐标，f0,θ0f_0,\theta_0f0,θ0分别为滤波器中心频率和方向角度，参数可σr,σθ\sigma_r,\sigma_\thetaσr,σθ分别用于决定滤波器的径向带宽和角度带宽。对提供波长的滤波器，频率由滤波器的波长设置，在不同尺度下，其转换公式为：

在这里我将6个方向4个尺度的Log-Gabor 滤波器做展示，滤波器与图像的大小相同，一般情况下将图像进行FFT变换后与滤波器对应位相乘，再转换至空域，就是滤波后的图像效果。

构造滤波器

import cv2
import numpy as np
import matplotlib.pyplot as pltn_scales=4
n_angles=6
img = cv2.imread("img_path")
H,W,_ = ir_img.shapedef lowpassfilter(H, W, cutoff, n):if cutoff < 0 or cutoff > 0.5:raise ValueError('the cutoff frequency needs to be between 0 and 0.5')if not n == int(n) or n < 1.0:raise ValueError('n must be an integer >= 1')xrange = np.linspace(-0.5, 0.5, W)yrange = np.linspace(-0.5, 0.5, H)x, y = np.meshgrid(xrange, yrange)radius = np.sqrt(x ** 2 + y ** 2)radius = np.fft.ifftshift(radius)return 1.0 / (1.0 + (radius / cutoff) ** (2 * n))xrange = np.linspace(-0.5, 0.5, W)
yrange = np.linspace(-0.5, 0.5, H)
x, y = np.meshgrid(xrange, yrange)
radius = np.sqrt(x ** 2 + y ** 2)
theta = np.arctan2(-y, x)# numpy.fft模块中的fftshift函数可以将FFT输出中的直流分量移动到频谱的中央。ifftshift函数则是其逆操作
radius = np.fft.ifftshift(radius)
theta = np.fft.ifftshift(theta)sintheta = np.sin(theta)
costheta = np.cos(theta)lp_filter = lowpassfilter(H, W, 0.45, 15)
# lp_filter = np.fft.ifftshift(lp_filter)
log_gabor_filters = np.zeros((n_scales, H, W))#  不同尺度
for sc in range(n_scales):wavelength = min_wavelength * multiplier ** sclog_gabor_filters[sc] = np.exp((-(np.log(radius * wavelength + 1e-5)) ** 2) / (2 * np.log(sigma_onf + 1e-5) ** 2)) * lp_filterspreads = np.zeros((n_angles, H, W))
#  不同方向
for o in range(n_angles):angle = o * np.pi / n_anglesds = sintheta * np.cos(angle) - costheta * np.sin(angle)dc = costheta * np.cos(angle) + sintheta * np.sin(angle)dtheta = abs(np.arctan2(ds, dc))dtheta = np.minimum(dtheta * n_angles * 0.5, np.pi)spreads[o] = (np.cos(dtheta) + 1) / 2
#  构造集合的filter
filter_bank = np.zeros((n_scales * n_angles, H, W))
for sc in range(n_scales):for o in range(n_angles):filter_bank[sc * n_angles + o] = log_gabor_filters[sc] * spreads[o]#  可视化
for sc in range(n_scales):plt.figure(sc,figsize=(30,120))first = filter_bank[sc * n_angles] - np.min(filter_bank[sc * n_angles])first /= np.max(first)first *= 255shuiping = firstfor o in range(n_angles-1):out = filter_bank[sc * n_angles + o+1] - np.min(filter_bank[sc * n_angles + o+1])out /= np.max(out)out *= 255shuiping = np.hstack((shuiping,out))plt.imshow(shuiping,cmap="gray")

滤波器可视化

进行图像测试

原图

vi_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
# 图像的傅里叶变换
image_fft = cv2.dft(np.float32(vi_gray), flags=cv2.DFT_COMPLEX_OUTPUT)
eo = np.zeros((filter_bank.shape[0], filter_bank.shape[1], filter_bank.shape[2], 2))
for i, filter in enumerate(filter_bank):eo[i] = cv2.idft(np.multiply(np.expand_dims(filter, -1), image_fft))for sc in range(n_scales):plt.figure(sc,figsize=(30,120))first = eo_magnitude[n_scales]first /= np.max(first)first *= 255shuiping = firstfor o in range(n_angles-1):out = eo_magnitude[sc * n_angles + o+1] - np.min(eo_magnitude[sc * n_angles + o+1])out /= np.max(out)out *= 255shuiping = np.hstack((shuiping,out))plt.imshow(shuiping,cmap="gray")

滤波后