我试图用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包装器的官方文档。函数的docstring提供了一些信息,但它是不完整和不精确的
我发现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慢得多
关于这些功能,我还有进一步的问题:
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()
MATLAB的^{} 和OpenCV的^{} 的文档几乎都足够完整,可以进行1:1的翻译。只需要一点点实验就可以找出如何将MATLAB的“
SpatialFrequencyBandwidth
”转换为高斯包络的西格玛我在这里注意到的一件事是OpenCV对Gabor滤波的实现似乎表明Gabor滤波器还没有被很好地理解。Google的一个快速练习表明OpenCV中最流行的Gabor过滤教程没有正确理解Gabor过滤器
Gabor滤波器是一个复数滤波器,例如可以从OpenCV文档链接到的同一个Wikipedia page中学习到。因此,将其应用于图像的结果也是复杂的。MATLAB正确地返回复数结果的幅度和相位,而不是复数图像本身,因为它主要是感兴趣的幅度。Gabor滤波器的大小指示图像的哪些部分具有给定波长和方向的频率
例如,可以对该图像(左)应用Gabor滤波器以产生该结果(右)(这是复数输出的幅度):
然而,OpenCV的过滤似乎是严格实值的。可以构建具有任意相位的Gabor滤波器核的实值分量。Gabor滤波器的实部为0相位,虚部为π/2相位(即实部为偶数,虚部为奇数)。结合偶数和奇数滤波器可以分析具有任意相位的信号,不需要创建具有其他相位的滤波器
要复制以下MATLAB代码:
在使用OpenCV的Python中,需要执行以下操作:
请注意,输入图像必须为浮点类型,否则计算结果将转换为无法表示表示Gabor滤波器结果所需的所有值的类型
OP中的最后一行代码是
对我来说,这是非常奇怪的,我想知道这个代码是用来做什么的。它完成的是获得Gabor滤波器的虚部的结果:
相关问题 更多 >
编程相关推荐