Introduction to Neural Networks

Programming Lab 8: Fundamental models of the primate visual system


Getting started

As before, we start by importing the libraries we will use:
import numpy as np 
import matplotlib.pyplot as plt
from skimage.color import rgb2gray
from scipy.signal import convolve2d

Gabor filters

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])    
	

Applying Gabor filters to images

Next, we are going to apply our Gabor filters to images. Download the images from the course website and unzip them.
# 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.

Fig 1 Example output for a particular image and filter. Image size is 250 x 250 pixels, filter size is 21 x 21 pixels. To create the output, the filter is convolved with the image, i.e. we move the filter across the image and compute the dot product of the image intensities and the filter at each location.

The HMAX model: beyond V1

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.


Fig 2
Schematic of the HMAX model. During the programming lab, we will focus on the first four layers of HMAX, which roughly capture visual processing from V1 to V4 in the primate visual system. Note that the first layer (S1) is identical to the bank of Gabor filters we created today. Figure adapted from Serre et al. 2007.

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.