python - 在 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 ，其中C++中使用OpenCV来解决问题。我假设这些函数的工作方式非常相似(还要注意 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 和 OpenCV 的 getGaborKernel 几乎可以完成 1:1 的翻译。只需要一点点实验就可以弄清楚如何将 MATLAB 的“SpatialFrequencyBandwidth”转换为高斯包络的 sigma。

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

Gabor 滤波器，例如可以从相同的 Wikipedia page 中学习OpenCV 的文档链接到的是一个复数值过滤器。因此，将其应用于图像的结果也是复值的。 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)