Skip to content

Instantly share code, notes, and snippets.

@thouis
Forked from amueller/filterbank.py
Created July 23, 2012 11:39
Show Gist options
  • Save thouis/3163187 to your computer and use it in GitHub Desktop.
Save thouis/3163187 to your computer and use it in GitHub Desktop.
Filterbank responses for low level vision
##########################################################################
# Maximum Response filterbank from
# http://www.robots.ox.ac.uk/~vgg/research/texclass/filters.html
# based on several edge and bar filters.
# Adapted to Python by Andreas Mueller amueller@ais.uni-bonn.de
# Share and enjoy
#
import numpy as np
import matplotlib.pyplot as plt
from itertools import product, chain
from scipy.misc import lena
#from sklearn.externals.joblib import Parallel, delayed
def makeRFSfilters(radius=24, sigmas=[1, 2, 4], n_orientations=6):
""" Generates filters for RFS filterbank.
Parameters
----------
radius : int, default 28
radius of all filters. Size will be 2 * radius + 1
sigmas : list of floats, default [1, 2, 4]
define scales on which the filters will be computed
n_orientations : int
number of fractions the half-angle will be divided in
Returns
-------
edge : ndarray (len(sigmas), n_orientations, 2*radius+1, 2*radius+1)
Contains edge filters on different scales and orientations
bar : ndarray (len(sigmas), n_orientations, 2*radius+1, 2*radius+1)
Contains bar filters on different scales and orientations
rot : ndarray (2, 2*radius+1, 2*radius+1)
contains two rotation invariant filters, Gaussian and Laplacian of
Gaussian
"""
def make_gaussian_filter(x, sigma, order=0):
if order > 2:
raise ValueError("Only orders up to 2 are supported")
# compute unnormalized Gaussian response
response = np.exp(-x ** 2 / (2. * sigma ** 2))
if order == 1:
response = -response * x
elif order == 2:
response = response * (x ** 2 - sigma ** 2)
# normalize
response /= np.abs(response).sum()
return response
def makefilter(scale, phasey, pts, sup):
gx = make_gaussian_filter(pts[0, :], sigma=3 * scale)
gy = make_gaussian_filter(pts[1, :], sigma=scale, order=phasey)
f = (gx * gy).reshape(sup, sup)
# normalize
f /= np.abs(f).sum()
return f
support = 2 * radius + 1
x, y = np.mgrid[-radius:radius + 1, radius:-radius - 1:-1]
orgpts = np.vstack([x.ravel(), y.ravel()])
rot, edge, bar = [], [], []
for sigma in sigmas:
for orient in xrange(n_orientations):
# Not 2pi as filters have symmetry
angle = np.pi * orient / n_orientations
c, s = np.cos(angle), np.sin(angle)
rotpts = np.dot(np.array([[c, -s], [s, c]]), orgpts)
edge.append(makefilter(sigma, 1, rotpts, support))
bar.append(makefilter(sigma, 2, rotpts, support))
length = np.sqrt(x ** 2 + y ** 2)
rot.append(make_gaussian_filter(length, sigma=10))
rot.append(make_gaussian_filter(length, sigma=10, order=2))
# reshape rot and edge
edge = np.asarray(edge)
edge = edge.reshape(len(sigmas), n_orientations, support, support)
bar = np.asarray(bar).reshape(edge.shape)
rot = np.asarray(rot)[:, np.newaxis, :, :]
return edge, bar, rot
def apply_filterbank(img, filterbank):
from scipy.ndimage import convolve
result = []
for battery in filterbank:
response = [convolve(img, filt) for filt in battery]
#response = Parallel(n_jobs=5)(
#delayed(convolve)(img, filt) for filt in battery)
max_response = np.max(response, axis=0)
result.append(max_response)
print("battery finished")
return result
if __name__ == "__main__":
sigmas = [1, 2, 4]
n_sigmas = len(sigmas)
n_orientations = 6
edge, bar, rot = makeRFSfilters(sigmas=sigmas,
n_orientations=n_orientations)
n = n_sigmas * n_orientations
# plot filters
# 2 is for bar / edge, + 1 for rot
fig, ax = plt.subplots(n_sigmas * 2 + 1, n_orientations)
for k, filters in enumerate([bar, edge]):
for i, j in product(xrange(n_sigmas), xrange(n_orientations)):
row = i + k * n_sigmas
ax[row, j].imshow(filters[i, j, :, :], cmap=plt.cm.gray)
ax[row, j].set_xticks(())
ax[row, j].set_yticks(())
ax[-1, 0].imshow(rot[0, 0], cmap=plt.cm.gray)
ax[-1, 0].set_xticks(())
ax[-1, 0].set_yticks(())
ax[-1, 1].imshow(rot[1, 0], cmap=plt.cm.gray)
ax[-1, 1].set_xticks(())
ax[-1, 1].set_yticks(())
for i in xrange(2, n_orientations):
ax[-1, i].set_visible(False)
# apply filters to lena
img = lena().astype(np.float)
filterbank = chain(edge, bar, rot)
n_filters = len(edge) + len(bar) + len(rot)
response = apply_filterbank(img, filterbank)
# plot responses
fig2, ax2 = plt.subplots(3, 3)
for axes, res in zip(ax2.ravel(), response):
axes.imshow(res, cmap=plt.cm.gray)
axes.set_xticks(())
axes.set_yticks(())
ax2[-1, -1].set_visible(False)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment