在 Python 中重现 MATLAB 的 imgaborfilt

Question

我正在尝试在 Python 中重现以下 MATLAB 代码的行为：

% Matlab code
wavelength = 10
orientation = 45
image = imread('filename.tif') % grayscale image
[mag,phase] = imgaborfilt(image, wavelength, orientation)
gabor_im = mag .* sin(phase)

很遗憾，我没有许可证，无法运行代码。此外，official Matlab documentation of imgaborfilt 并未准确说明函数的作用。

由于缺乏明显的替代方案，我尝试在 Python 中使用 OpenCV（接受其他建议）。我没有使用 OpenCV 的经验。我正在尝试使用 cv2.getGaborKernel 和 cv2.filter2D。我也找不到这些函数行为的详细文档。 Afaik 没有关于 OpenCV 的 Python 包装器的官方文档。函数的文档字符串提供了一些信息，但它不完整且不精确。

我找到了 this question, where OpenCV is used in C++ to solve the problem. I assume the functions work in a very similar way (also note the official C++ documentation)。但是，它们有许多附加参数。 我如何找出 matlab 函数真正做了什么 来重现该行为？

# python 3.6
import numpy as np
import cv2

wavelength = 10
orientation = 45
shape = (500, 400)  # arbitrary values to get running example code...
sigma = 100  # what to put for Matlab behaviour?
gamma = 1  # what to put for Matlab behaviour?
gabor_filter = cv2.getGaborKernel(shape, sigma, orientation, wavelength, gamma)
print(gabor_filter.shape)  # =(401, 501). Why flipped?

image = np.random.random(shape)  # create some random data.
out_buffer = np.zeros(shape)

destination_depth = -1  # like dtype for filter2D. Apparantly, -1="same as input".
thing = cv2.filter2D(image, destination_depth, gabor_filter, out_buffer)
print(out_buffer.shape, out_buffer.dtype, out_buffer.max())  # =(500, 400) float64 65.2..
print(thing.shape, thing.dtype, thing.max())  # =(500, 400) float64 65.2..

编辑：

在收到 Cris Luengo 的精彩回答后，我用它制作了两个函数，分别使用 OpenCV 和 scikit-image，（希望）重现 MATLAB imgaborfit 函数行为。我把它们包括在这里。请注意，scikit 实现比 OpenCV 慢很多。

我对这些函数还有进一步的疑问：

精确到什么程度 OpenCV解法和MATLAB解法结果一致？
对于不想使用 OpenCV 的人，我还提供了一个 scikit-image 解决方案这里。我找到的参数，使得幅度几乎相等。但是，scikit-image 解决方案的阶段似乎与 OpenCV 解决方案不同。为什么是这样？

import numpy as np
import math
import cv2

def gaborfilt_OpenCV_likeMATLAB(image, wavelength, orientation, SpatialFrequencyBandwidth=1, SpatialAspectRatio=0.5):
    """Reproduces (to what accuracy in what MATLAB version??? todo TEST THIS!) the behaviour of MATLAB imgaborfilt function using OpenCV."""

    orientation = -orientation / 180 * math.pi # for OpenCV need radian, and runs in opposite direction
    sigma = 0.5 * wavelength * SpatialFrequencyBandwidth
    gamma = SpatialAspectRatio
    shape = 1 + 2 * math.ceil(4 * sigma)  # smaller cutoff is possible for speed
    shape = (shape, shape)
    gabor_filter_real = cv2.getGaborKernel(shape, sigma, orientation, wavelength, gamma, psi=0)
    gabor_filter_imag = cv2.getGaborKernel(shape, sigma, orientation, wavelength, gamma, psi=math.pi / 2)
    filtered_image = cv2.filter2D(image, -1, gabor_filter_real) + 1j * cv2.filter2D(image, -1, gabor_filter_imag)
    mag = np.abs(filtered_image)
    phase = np.angle(filtered_image)
    return mag, phase

import numpy as np
import math
from skimage.filters import gabor

def gaborfilt_skimage_likeMATLAB(image, wavelength, orientation, SpatialFrequencyBandwidth=1, SpatialAspectRatio=0.5):
    """TODO (does not quite) reproduce the behaviour of MATLAB imgaborfilt function using skimage."""
    sigma = 0.5 * wavelength * SpatialFrequencyBandwidth
    filtered_image_re, filtered_image_im = gabor(
        image, frequency=1 / wavelength, theta=-orientation / 180 * math.pi,
        sigma_x=sigma, sigma_y=sigma/SpatialAspectRatio, n_stds=5,
    )
    full_image = filtered_image_re + 1j * filtered_image_im
    mag = np.abs(full_image)
    phase = np.angle(full_image)
    return mag, phase

测试上述功能的代码：

from matplotlib import pyplot as plt
import numpy as np

def show(im, title=""):
    plt.figure()
    plt.imshow(im)
    plt.title(f"{title}: dtype={im.dtype}, shape={im.shape},\n max={im.max():.3e}, min= {im.min():.3e}")
    plt.colorbar()

image = np.zeros((400, 400))
image[200, 200] = 1  # a delta impulse image to visualize the filtering kernel
wavelength = 10
orientation = 33  # in degrees (for MATLAB)

mag_cv, phase_cv = gaborfilt_OpenCV_likeMATLAB(image, wavelength, orientation)
show(mag_cv, "mag")  # normalized by maximum, non-zero noise even outside filter window region
show(phase_cv, "phase")  # all over the place

mag_sk, phase_sk = gaborfilt_skimage_likeMATLAB(image, wavelength, orientation)
show(mag_sk, "mag skimage")  # small values, zero outside filter region
show(phase_sk, "phase skimage")  # and hence non-zero only inside filter window region

show(mag_cv - mag_sk/mag_sk.max(), "cv - normalized(sk)")  # approximately zero-image.
show(phase_sk - phase_cv, "phase_sk - phase_cv") # phases do not agree at all! Not even in the window region!
plt.show()

Answer 1

两个 MATLAB imgaborfilt and OpenCV's getGaborKernel 的文档几乎已经完成，足以进行 1:1 翻译。只需要一点点实验就可以弄清楚如何将 MATLAB 的“SpatialFrequencyBandwidth”转换为高斯包络的西格玛。

我在这里注意到的一件事是，OpenCV 对 Gabor 过滤的实现似乎表明 Gabor 过滤器没有得到很好的理解。快速 Google 练习表明 OpenCV 中最流行的 Gabor 过滤教程没有正确理解 Gabor 过滤器。

Gabor 过滤器，例如可以从 OpenCV 的文档链接到的相同 Wikipedia page 中了解到，它是一个复值过滤器。因此，将其应用于图像的结果也是复值的。 MATLAB 正确 returns 复数结果的幅度和相位，而不是复值图像本身，因为它主要是感兴趣的幅度。 Gabor 滤波器的大小表示图像的哪些部分具有给定波长和方向的频率。

例如，可以将 Gabor 滤波器应用于此图像（左）以产生此结果（右）（这是复值输出的幅度）：

然而，OpenCV 的过滤似乎是严格实值的。可以构建具有任意相位的 Gabor 滤波器内核的实值分量。 Gabor 滤波器有一个相位为 0 的实部和一个相位为 π/2 的虚部（即实部为偶数，虚部为奇数）。结合偶数和奇数滤波器可以分析任意相位的信号，无需创建其他相位的滤波器。

要复制以下 MATLAB 代码：

image = zeros(64,64); 
image(33,33) = 1;     % a delta impulse image to visualize the filtering kernel

wavelength = 10;
orientation = 30; # in degrees
[mag,phase] = imgaborfilt(image, wavelength, orientation);
% defaults: 'SpatialFrequencyBandwidth'=1; 'SpatialAspectRatio'=0.5

在 Python 中，使用 OpenCV 需要做的是：

import cv2
import numpy as np
import math

image = np.zeros((64, 64))
image[32, 32] = 1          # a delta impulse image to visualize the filtering kernel

wavelength = 10
orientation = -30 / 180 * math.pi    # in radian, and seems to run in opposite direction
sigma = 0.5 * wavelength * 1         # 1 == SpatialFrequencyBandwidth
gamma = 0.5                          # SpatialAspectRatio
shape = 1 + 2 * math.ceil(4 * sigma) # smaller cutoff is possible for speed
shape = (shape, shape)
gabor_filter_real = cv2.getGaborKernel(shape, sigma, orientation, wavelength, gamma, psi=0)
gabor_filter_imag = cv2.getGaborKernel(shape, sigma, orientation, wavelength, gamma, psi=math.pi/2)

gabor = cv2.filter2D(image, -1, gabor_filter_real) + 1j * cv2.filter2D(image, -1, gabor_filter_imag)
mag = np.abs(gabor)
phase = np.angle(gabor)

请注意，输入图像必须是浮点类型，否则计算结果将转换为无法表示表示 Gabor 滤波器结果所需的所有值的类型。

OP 中代码的最后一行是

gabor_im = mag .* sin(phase)

这对我来说很奇怪，我想知道这段代码的用途。它完成的是获取Gabor滤波器虚部的结果：

gabor_im = cv2.filter2D(image, -1, gabor_filter_imag)

在 Python 中重现 MATLAB 的 imgaborfilt

Reproduce MATLAB's imgaborfilt in Python

python

matlab

opencv

scikit-image

gabor-filter