```
import numpy as np
import matplotlib.pyplot as plt
from skimage.color import rgb2gray
from scipy.signal import convolve2d
```

In the early 1960s, David Hubel and Torsten Wiesel performed a set of seminal experiments in visual neuroscience. They systematically examined the response properties of neurons in primary visual cortex (V1) in cats and monkeys, and discovered that neurons are tuned to particular spatial locations and orientations across a range of spatial scales. Their discoveries on information processing in the visual system were awarded with a Nobel Prize in 1981. The response functions reported by Hubel and Wiesel can be described well by Gabor filters. Gabor filters are spatially localized, oriented, and bandpass. They form a basis set for sparse coding of natural images, as shown by Olshausen and Field in 1996.

A Gabor filter is defined by the equation:

\begin{equation} (x,y;\lambda , \theta , \psi , \sigma , \gamma) = \exp (-\frac{{x}'^{2}+\gamma^{2}{y}'^{2}}{2\sigma^{2}})\cos (2\pi\frac{{x}'}{\lambda} + \psi) \tag{1}\label{gabor} \end{equation}where

\begin{equation} {x}' = x \cos \theta + y \sin \theta \tag{2}\label{xprime} \end{equation}and

\begin{equation} {y}' = -x \sin \theta + y \cos \theta \tag{3}\label{yprime} \end{equation}where \( \lambda\) is the wavelength of the sinusoid, \( \theta\) is the orientation of the filter, \( \psi\) is the phase offset of the sinusoid, \( \sigma\) is the standard deviation of the Gaussian kernel, and \( \gamma\) is the spatial aspect ratio of the Gaussian kernel. This equation shows that a Gabor filter is a Gaussian kernel modulated by a sinusoidal wave.

**Step 1: Build our own Gabor filters**

```
# build gabor filter
def gabor_filter(size, wavelength, orientation, monitor):
# size = filter size in pixels (w = h)
# wavelength = wavelength of sinusoid relative to half the size of the filter
# orientation = desired orientation of filter
# set parameters (based on neural data from V1)
lambda_ = size * 2. / wavelength # wavelength of sinusoidal plane wave (how many waves 'fit' in the filter)
sigma = lambda_ * 0.8 # standard deviation of gaussian kernel
gamma = 0.3 # spatial aspect ratio of gaussian kernel
theta = np.deg2rad(orientation + 90) # orientation
# gaussian
x, y = np.mgrid[:size, :size] - (size // 2)
rotx = x * np.cos(theta) + y * np.sin(theta)
roty = -x * np.sin(theta) + y * np.cos(theta)
gauss = np.exp(-(rotx**2 + gamma**2 * roty**2) / (2 * sigma ** 2))
# sinusoid
sinusoid = np.cos(2 * np.pi * rotx / lambda_)
# gabor
filt = gauss * sinusoid
filt[np.sqrt(x**2 + y**2) > (size / 2)] = 0
filt = filt - np.mean(filt)
filt = filt / np.sqrt(np.sum(filt ** 2))
# show
if monitor:
plt.figure(); plt.imshow(gauss, cmap=plt.gray()); plt.title('2D Gaussian kernel')
plt.figure(); plt.imshow(sinusoid, cmap=plt.gray()); plt.title('Sinusoidal plane wave')
plt.figure(); plt.imshow(filt, cmap=plt.gray()); plt.title('Gabor filter')
return filt
# generate a filter with the following parameters:
# size = 23, wavelength = 3.6, orientation = -45 degrees
filt = gabor_filter(23, 3.6, -45, 1)
```

**Step 2: Create a bank of Gabor filters**

```
# generate a filter bank
sizes = np.arange(7, 39, 2) # 16 sizes
wavelengths = np.arange(4, 3.2, -0.05) # 16 associated wavelengths
orientations = np.arange(-45, 135, 45) # 4 orientations
params = []
i = 0;
for s in sizes:
i = i + 1
for o in orientations:
w = wavelengths[i-1]
params.append((s,w,o))
filterBank = []
gaborParams = []
for (size, wavelength, orientation) in params:
gaborParam = {'size':size, 'wavelength':wavelength, 'orientation':orientation, 'monitor':0}
filt = gabor_filter(**gaborParam)
filterBank.append(filt)
gaborParams.append(gaborParam)
# show filter bank
plt.figure()
n = len(filterBank)
for i in range(n):
plt.subplot(16,4,i+1)
plt.axis('off'); plt.imshow(filterBank[i])
```

```
# read in images
face = rgb2gray(plt.imread('images/image19_66440000_rs.jpg'))
zebra = rgb2gray(plt.imread('images/image32_41110000_rs.jpg'))
# show image and filter
plt.figure(); plt.imshow(zebra, cmap=plt.gray()); plt.title('image')
filt = filterBank[30]
plt.figure(); plt.imshow(filt, cmap=plt.gray()); plt.title('filter')
# convolve filter with image (convolution = dot product)
%time output = convolve2d(zebra, filt, mode='valid')
plt.figure(); plt.imshow(output, cmap=plt.gray()); plt.title('response')
```

Play around with different images and filters.

Gabor filters model early stages of visual processing. To extend our model beyond V1, we need to find ways of creating selectivity for increasingly complex shapes as well as invariance to position and scale. The HMAX model implements this extension. Importantly, the filters in the early stages of HMAX are manually engineered to simulate response properties of neurons in V1 and V2 (not learned through experience). Today, we will run HMAX on our images.

Python code for HMAX can be downloaded here. Furthermore, we need to install pytorch to run the code. To reduce processing time, select 5 images from our image set (images 23, 27, 32, 62, and 84) and save them in a subfolder of the folder 'example_images' that comes with the HMAX code. Move the original example images to a different folder.

```
# run HMAX code
cd hmax_path # replace with local path
run example.py
# if you want to run the code for all 96 images, change the batch size in example.py
dataloader = DataLoader(example_images, batch_size=96)
```

The code ouputs a pickle file with the image activations, which we will examine in this week's problem set.