python - 在 Python 中重现 MATLAB 的 imgaborfilt
问题描述
我正在尝试在 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
and 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()
解决方案
MATLABimgaborfilt
和 OpenCV 的文档getGaborKernel
几乎都足够完整,可以进行 1:1 的翻译。只需要一点点实验就可以弄清楚如何将 MATLAB 的 " SpatialFrequencyBandwidth
" 转换为高斯包络的 sigma。
我在这里注意到的一件事是 OpenCV 的 Gabor 过滤器实现似乎表明 Gabor 过滤器还没有被很好地理解。一个快速的 Google 练习表明,OpenCV 中最流行的 Gabor 过滤教程没有正确理解 Gabor 过滤器。
Gabor 过滤器,例如可以从OpenCV 的文档链接到的同一个Wikipedia 页面中了解到,是一个复值过滤器。因此,将其应用于图像的结果也是复值。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)
推荐阅读
- git - 从功能分支推送推送不同的分支?
- pandas - 如何使用 Pandas 重新排序多索引列?
- python - 烧瓶应用的模型模块
- javascript - 结合python和JS(mapbox)进行交互渲染
- selenium-webdriver - protractor/webDriver 中的 isDisplayed() 是如何工作的?
- performance - 需要一种方法来创建多个(例如每个 250 个)google 和 Outlook 测试帐户,以对 oauth 相关的 api 进行负载测试。(例如:gmail api)
- neo4j - 类型不匹配:已使用冲突类型 Any(预期节点)定义了 connectedcompany
- regex - 正则表达式如果第二个参数是注册
- python - 实时检查 2 个 QLineEdits 文本
- python - 访问熊猫数据框中的单元格