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

不幸的是，我没有许可证，无法运行代码。此外，imgaborfilt的官方 Matlab 文档没有准确指定函数的作用。

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

我发现了这个问题，其中在 C++ 中使用 OpenCV 来解决问题。我假设这些函数的工作方式非常相似（另请注意 官方 C++ 文档）。但是，它们有许多附加参数。我怎样才能找出 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

MATLABimgaborfilt和 OpenCV的文档getGaborKernel几乎足够完整，能够进行 1:1 翻译。只需进行一点实验即可弄清楚如何将 MATLAB 的“ SpatialFrequencyBandwidth”转换为高斯包络的西格玛。

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

Gabor 过滤器，例如可以从同一维基百科页面了解到了解到）是一个复值过滤器。因此，将其应用于图像的结果也是复值的。MATLAB 正确返回复数结果的幅度和相位，而不是复数值图像本身，因为它主要是我们感兴趣的幅度。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

在使用 OpenCV 的 Python 中，需要执行以下操作：

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)