实施 log Gabor 滤波器组
Implementing log Gabor filter bank
我正在阅读这篇论文"Self-Invertible 2D Log-Gabor Wavelets"它定义了 2D log gabor 滤波器:
论文还指出滤波器仅覆盖频率的一侧 space 并在这张图片中显示
在我尝试实施过滤器时,我得到的结果与论文中所说的不符。让我从我的实现开始,然后我将说明问题。
实施:
我创建了一个包含滤波器的二维数组并转换了每个索引,以便频域的原点位于数组的中心,正 x-axis 向右和正 y-axis 上升。
number_scales = 5 # scale resolution
number_orientations = 9 # orientation resolution
N = constantDim # image dimensions
def getLogGaborKernal(scale, angle, logfun=math.log2, norm = True):
# setup up filter configuration
center_scale = logfun(N) - scale
center_angle = ((np.pi/number_orientations) * angle) if (scale % 2) \
else ((np.pi/number_orientations) * (angle+0.5))
scale_bandwidth = 0.996 * math.sqrt(2/3)
angle_bandwidth = 0.996 * (1/math.sqrt(2)) * (np.pi/number_orientations)
# 2d array that will hold the filter
kernel = np.zeros((N, N))
# get the center of the 2d array so we can shift origin
middle = math.ceil((N/2)+0.1)-1
# calculate the filter
for x in range(0,constantDim):
for y in range(0,constantDim):
# get the transformed x and y where origin is at center
# and positive x-axis goes right while positive y-axis goes up
x_t, y_t = (x-middle),-(y-middle)
# calculate the filter value at given index
kernel[y,x] = logGaborValue(x_t,y_t,center_scale,center_angle,
scale_bandwidth, angle_bandwidth,logfun)
# normalize the filter energy
if norm:
Kernel = kernel / np.sum(kernel**2)
return kernel
为了计算每个索引处的过滤器值,在我们转到 log-polar space
的地方进行另一个转换
def logGaborValue(x,y,center_scale,center_angle,scale_bandwidth,
angle_bandwidth, logfun):
# transform to polar coordinates
raw, theta = getPolar(x,y)
# if we are at the center, return 0 as in the log space
# zero is not defined
if raw == 0:
return 0
# go to log polar coordinates
raw = logfun(raw)
# calculate (theta-center_theta), we calculate cos(theta-center_theta)
# and sin(theta-center_theta) then use atan to get the required value,
# this way we can eliminate the angular distance wrap around problem
costheta, sintheta = math.cos(theta), math.sin(theta)
ds = sintheta * math.cos(center_angle) - costheta * math.sin(center_angle)
dc = costheta * math.cos(center_angle) + sintheta * math.sin(center_angle)
dtheta = math.atan2(ds,dc)
# final value, multiply the radial component by the angular one
return math.exp(-0.5 * ((raw-center_scale) / scale_bandwidth)**2) * \
math.exp(-0.5 * (dtheta/angle_bandwidth)**2)
问题:
角度:论文指出,从 1->8 索引角度会产生良好的方向覆盖,但在我的实现中角度从 1 ->n 除了半方向外不覆盖。甚至垂直方向也没有被正确覆盖。这可以在此图中显示,其中包含比例为 3 的过滤器集和方向范围为 1->8:
覆盖范围: 从上面的过滤器可以清楚地看出,过滤器覆盖了 space 的两侧,这不是论文所说的。这可以通过使用范围从 -4 -> 4 的 9 个方向变得更加明确。下图包含一张图像中的所有过滤器,以显示它如何覆盖光谱的两侧(该图像是通过在每个位置取最大值创建的来自所有过滤器):
Middle Column (orientation $\pi / 2$):第一张图中orientation from 3 -> 8可以看出filter在 $ \pi / 2$ 方向消失。这是正常的吗?当我将所有过滤器(所有 5 个尺度和 9 个方向)组合在一张图像中时,也可以看到这一点:
更新:
在空间域中添加滤波器的脉冲响应,可以看到 -4 和 4 方向有明显的失真:
经过大量的代码分析,我发现我的实现是正确的,但是 getPolar
功能被搞砸了,所以上面的代码应该可以正常工作。这是没有 getPolar
功能的新代码,如果有人正在寻找的话:
number_scales = 5 # scale resolution
number_orientations = 8 # orientation resolution
N = 128 # image dimensions
def getFilter(f_0, theta_0):
# filter configuration
scale_bandwidth = 0.996 * math.sqrt(2/3)
angle_bandwidth = 0.996 * (1/math.sqrt(2)) * (np.pi/number_orientations)
# x,y grid
extent = np.arange(-N/2, N/2 + N%2)
x, y = np.meshgrid(extent,extent)
mid = int(N/2)
## orientation component ##
theta = np.arctan2(y,x)
center_angle = ((np.pi/number_orientations) * theta_0) if (f_0 % 2) \
else ((np.pi/number_orientations) * (theta_0+0.5))
# calculate (theta-center_theta), we calculate cos(theta-center_theta)
# and sin(theta-center_theta) then use atan to get the required value,
# this way we can eliminate the angular distance wrap around problem
costheta = np.cos(theta)
sintheta = np.sin(theta)
ds = sintheta * math.cos(center_angle) - costheta * math.sin(center_angle)
dc = costheta * math.cos(center_angle) + sintheta * math.sin(center_angle)
dtheta = np.arctan2(ds,dc)
orientation_component = np.exp(-0.5 * (dtheta/angle_bandwidth)**2)
## frequency componenet ##
# go to polar space
raw = np.sqrt(x**2+y**2)
# set origin to 1 as in the log space zero is not defined
raw[mid,mid] = 1
# go to log space
raw = np.log2(raw)
center_scale = math.log2(N) - f_0
draw = raw-center_scale
frequency_component = np.exp(-0.5 * (draw/ scale_bandwidth)**2)
# reset origin to zero (not needed as it is already 0?)
frequency_component[mid,mid] = 0
return frequency_component * orientation_component
我正在阅读这篇论文"Self-Invertible 2D Log-Gabor Wavelets"它定义了 2D log gabor 滤波器:
论文还指出滤波器仅覆盖频率的一侧 space 并在这张图片中显示
在我尝试实施过滤器时,我得到的结果与论文中所说的不符。让我从我的实现开始,然后我将说明问题。
实施:
我创建了一个包含滤波器的二维数组并转换了每个索引,以便频域的原点位于数组的中心,正 x-axis 向右和正 y-axis 上升。
number_scales = 5 # scale resolution number_orientations = 9 # orientation resolution N = constantDim # image dimensions def getLogGaborKernal(scale, angle, logfun=math.log2, norm = True): # setup up filter configuration center_scale = logfun(N) - scale center_angle = ((np.pi/number_orientations) * angle) if (scale % 2) \ else ((np.pi/number_orientations) * (angle+0.5)) scale_bandwidth = 0.996 * math.sqrt(2/3) angle_bandwidth = 0.996 * (1/math.sqrt(2)) * (np.pi/number_orientations) # 2d array that will hold the filter kernel = np.zeros((N, N)) # get the center of the 2d array so we can shift origin middle = math.ceil((N/2)+0.1)-1 # calculate the filter for x in range(0,constantDim): for y in range(0,constantDim): # get the transformed x and y where origin is at center # and positive x-axis goes right while positive y-axis goes up x_t, y_t = (x-middle),-(y-middle) # calculate the filter value at given index kernel[y,x] = logGaborValue(x_t,y_t,center_scale,center_angle, scale_bandwidth, angle_bandwidth,logfun) # normalize the filter energy if norm: Kernel = kernel / np.sum(kernel**2) return kernel
为了计算每个索引处的过滤器值,在我们转到 log-polar space
的地方进行另一个转换def logGaborValue(x,y,center_scale,center_angle,scale_bandwidth, angle_bandwidth, logfun): # transform to polar coordinates raw, theta = getPolar(x,y) # if we are at the center, return 0 as in the log space # zero is not defined if raw == 0: return 0 # go to log polar coordinates raw = logfun(raw) # calculate (theta-center_theta), we calculate cos(theta-center_theta) # and sin(theta-center_theta) then use atan to get the required value, # this way we can eliminate the angular distance wrap around problem costheta, sintheta = math.cos(theta), math.sin(theta) ds = sintheta * math.cos(center_angle) - costheta * math.sin(center_angle) dc = costheta * math.cos(center_angle) + sintheta * math.sin(center_angle) dtheta = math.atan2(ds,dc) # final value, multiply the radial component by the angular one return math.exp(-0.5 * ((raw-center_scale) / scale_bandwidth)**2) * \ math.exp(-0.5 * (dtheta/angle_bandwidth)**2)
问题:
角度:论文指出,从 1->8 索引角度会产生良好的方向覆盖,但在我的实现中角度从 1 ->n 除了半方向外不覆盖。甚至垂直方向也没有被正确覆盖。这可以在此图中显示,其中包含比例为 3 的过滤器集和方向范围为 1->8:
覆盖范围: 从上面的过滤器可以清楚地看出,过滤器覆盖了 space 的两侧,这不是论文所说的。这可以通过使用范围从 -4 -> 4 的 9 个方向变得更加明确。下图包含一张图像中的所有过滤器,以显示它如何覆盖光谱的两侧(该图像是通过在每个位置取最大值创建的来自所有过滤器):
Middle Column (orientation $\pi / 2$):第一张图中orientation from 3 -> 8可以看出filter在 $ \pi / 2$ 方向消失。这是正常的吗?当我将所有过滤器(所有 5 个尺度和 9 个方向)组合在一张图像中时,也可以看到这一点:
更新: 在空间域中添加滤波器的脉冲响应,可以看到 -4 和 4 方向有明显的失真:
经过大量的代码分析,我发现我的实现是正确的,但是 getPolar
功能被搞砸了,所以上面的代码应该可以正常工作。这是没有 getPolar
功能的新代码,如果有人正在寻找的话:
number_scales = 5 # scale resolution
number_orientations = 8 # orientation resolution
N = 128 # image dimensions
def getFilter(f_0, theta_0):
# filter configuration
scale_bandwidth = 0.996 * math.sqrt(2/3)
angle_bandwidth = 0.996 * (1/math.sqrt(2)) * (np.pi/number_orientations)
# x,y grid
extent = np.arange(-N/2, N/2 + N%2)
x, y = np.meshgrid(extent,extent)
mid = int(N/2)
## orientation component ##
theta = np.arctan2(y,x)
center_angle = ((np.pi/number_orientations) * theta_0) if (f_0 % 2) \
else ((np.pi/number_orientations) * (theta_0+0.5))
# calculate (theta-center_theta), we calculate cos(theta-center_theta)
# and sin(theta-center_theta) then use atan to get the required value,
# this way we can eliminate the angular distance wrap around problem
costheta = np.cos(theta)
sintheta = np.sin(theta)
ds = sintheta * math.cos(center_angle) - costheta * math.sin(center_angle)
dc = costheta * math.cos(center_angle) + sintheta * math.sin(center_angle)
dtheta = np.arctan2(ds,dc)
orientation_component = np.exp(-0.5 * (dtheta/angle_bandwidth)**2)
## frequency componenet ##
# go to polar space
raw = np.sqrt(x**2+y**2)
# set origin to 1 as in the log space zero is not defined
raw[mid,mid] = 1
# go to log space
raw = np.log2(raw)
center_scale = math.log2(N) - f_0
draw = raw-center_scale
frequency_component = np.exp(-0.5 * (draw/ scale_bandwidth)**2)
# reset origin to zero (not needed as it is already 0?)
frequency_component[mid,mid] = 0
return frequency_component * orientation_component