Statistics of the natural input to a ring model

What does the input to a population of neurons in the primary visual cortex look like? In this post, we will try to have a feeling of the structure and statistics of the natural input to such a "ring" model.

This notebook explores this question using a retina-like temporal filtering and oriented Gabor-like filters. It produces this polar plot of the instantaneous energy in the different orientations for a natural movie :



One observes different striking features in the structure of this input to populations of V1 neurons:

  • input is sparse: often, a few orientations dominate - the shape of these components (bandwidth) seem to be similar,
  • there are many "switches": at some moments, the representations flips to another. This is due to cuts in the movie (changes from one scene to the other for instance). In a more realistic setting where we would add eye movements, these switches should also happen during saccades (but is there any knowledge of the occurence of the switch by the visual system?),
  • between switches, there is some coherence in amplitude (a component will slowly change its energy) but also in time (a component is more likely to have a ghradually changing oriantation, for instance when the scene rotates).

This structure is specific to the structure of natural images and to the way they transform (translations, rotations, zooms due to the motion and deformation of visual objects). This is certainly incorporated as a "prior" information in the structure of the visual cortex. As to know how and where this is implemented is an open scientific question.

This is joint work with Hugo Ladret.

Let's first initialize the notebook:

In [1]:
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
FIG_WIDTH = 8
FIGSIZE=(FIG_WIDTH, .618 * FIG_WIDTH)

loading images

First, I used a sequence of the "Life of Brian" movie which were converted to a numpy array using the imageio library.

In [2]:
video = np.load('/tmp/montypython.npy')
video = np.rot90(video)
#video = np.load('LifeOfBrian.npy')
#video = np.swapaxes(video, 0, -1)
#video = np.swapaxes(video, 0, 1)

N_frame = 8
N_frame = 1024
#video = video[:, 10:(10+36), :] # squaring
video = video[:, :, :N_frame]
N_X, N_Y, N_frame = video.shape
print('video.shape=', video.shape)
video.shape= (54, 54, 1024)
In [3]:
one_image = N_frame//2
image = video[:, :, one_image]
fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_WIDTH))
ax.imshow(image, cmap=plt.gray())
Out[3]:
<matplotlib.image.AxesImage at 0x1184a8908>

AR(1) filtering

Let's define different strategies for the temporal filtering, first with a moving average, aka AR(1):

In [4]:
alpha = -.61803
In [5]:
%%time
video_ar = np.zeros_like(video)
for i in range(N_frame):
    video_ar[:, :, i] = video[:, :, i] + alpha *video[:, :, i-1]
    
CPU times: user 109 ms, sys: 5.29 ms, total: 114 ms
Wall time: 113 ms
In [6]:
%%time
video_ar = video + alpha * np.roll(video, 1, axis=-1)
CPU times: user 12.5 ms, sys: 5.34 ms, total: 17.8 ms
Wall time: 16.6 ms
In [7]:
image = video_ar[:, :, one_image]
fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_WIDTH))
ax.imshow(image, cmap=plt.gray())
Out[7]:
<matplotlib.image.AxesImage at 0x101ed7fd0>

AR(2) filtering

What if we use a second-order linear filtering?

In [8]:
alpha, beta = -.61803, .61803-1. +.5
In [9]:
%%time
video_ar = video + alpha * np.roll(video, 1, axis=-1) + beta * np.roll(video, 2, axis=-1)
CPU times: user 20.6 ms, sys: 5.37 ms, total: 26 ms
Wall time: 24.8 ms
In [10]:
image = video_ar[:, :, one_image]
fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_WIDTH))
ax.imshow(image, cmap=plt.gray())
Out[10]:
<matplotlib.image.AxesImage at 0x10bc37978>

ARMA(2, 1) filtering

An even more realistic combination, yet still linear:

In [11]:
alpha, beta, gamma = -.61803, .61803-1., .5
In [12]:
%%time
video_ar = video + alpha * np.roll(video, 1, axis=-1) + beta * np.roll(video, 2, axis=-1)
for t in range(5):
    video_ar += gamma**t * np.roll(video_ar, t+1, axis=-1)
CPU times: user 60.5 ms, sys: 912 µs, total: 61.4 ms
Wall time: 60.1 ms
In [13]:
image = video_ar[:, :, one_image]
fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_WIDTH))
ax.imshow(image, cmap=plt.gray())
Out[13]:
<matplotlib.image.AxesImage at 0x11d683e80>
In [14]:
image = video_ar[:, :, one_image]
video = video_ar 

designing filters for the modeling of the feed-forward processing of V1 filters

Such a population of neurons may be rendered in an abstract way as a ring where the different orientations from $0$ to $\pi$ are arranged on the circle.

In [15]:
#from SLIP import Image, imread
from LogGabor import LogGabor
parameterfile = 'https://raw.githubusercontent.com/bicv/LogGabor/master/default_param.py'
lg = LogGabor(parameterfile)
lg.set_size((N_X, N_Y))

B_theta = 15 /180 * np.pi
params= {'sf_0':.12, 'B_sf': .5, 'B_theta': B_theta}
print('lg.pe.B_sf', lg.pe.B_sf)

N_theta, N_phase = 180, 2
thetas = np.linspace(0, np.pi, N_theta, endpoint=False)

phi = np.zeros((N_theta, N_phase, N_X, N_Y))

for i_theta, theta in enumerate(thetas):
    x = N_X/2
    y = N_Y/2
    params.update(theta=theta)
    for i_phase in range(N_phase):
        phase = i_phase * np.pi/2
        kernel = lg.invert(lg.loggabor(x, y, **params)*np.exp(-1j*phase))
        phi[i_theta, i_phase, :] = lg.normalize(kernel)

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_WIDTH))

ax.imshow(phi[i_theta, 0, :], cmap=plt.gray())

for i_theta in range(N_theta):
    env = np.sqrt(phi[i_theta, 0, :]**2 + phi[i_theta, 1, :]**2)
    ax.contour(env, levels=[env.max()/2], lw=1, colors=[plt.cm.viridis(i_theta/(N_theta))])
#fig.suptitle('Tiling of visual space using phi')
ax.set_xlabel(r'$Y$')
ax.set_ylabel(r'$X$')
ax.axis('equal')
plt.tight_layout()
#plt.savefig('logPolar_shape.pdf')
lg.pe.B_sf 0.4
In [16]:
print('precision in orientation (as defined by the bandwidth)=', '%.3f' % (lg.pe.B_theta*180/np.pi), 'deg')
precision in orientation (as defined by the bandwidth)= 10.000 deg

filtering images

Let's now apply these spatial filters to any temporally filtered input image:

In [17]:
image_ = (image*lg.mask).reshape((N_X*N_Y))
print('image_.shape', image_.shape)
image_.shape (2916,)
In [18]:
phi_ = phi.reshape((N_theta, N_phase, N_X*N_Y))
print('phi_.shape', phi_.shape)
phi_.shape (180, 2, 2916)
In [19]:
code = phi_ @ image_
print('All scalar products in one shot = ', (code).shape )
All scalar products in one shot =  (180, 2)
In [20]:
theta_bins = np.linspace(-np.pi/N_theta/2, -3*np.pi/N_theta/2+np.pi, N_theta)
print('theta_bins', theta_bins)
#binedges_theta = np.hstack((theta_bin, theta_bin[0]+np.pi))
#print('binedges_theta', binedges_theta)
width = np.pi/N_theta #binedges_theta[1:] - binedges_theta[:-1]
print('width', width)
theta_bins [-0.00872665  0.00872665  0.02617994  0.04363323  0.06108652  0.07853982
  0.09599311  0.1134464   0.13089969  0.14835299  0.16580628  0.18325957
  0.20071286  0.21816616  0.23561945  0.25307274  0.27052603  0.28797933
  0.30543262  0.32288591  0.3403392   0.3577925   0.37524579  0.39269908
  0.41015237  0.42760567  0.44505896  0.46251225  0.47996554  0.49741884
  0.51487213  0.53232542  0.54977871  0.56723201  0.5846853   0.60213859
  0.61959188  0.63704518  0.65449847  0.67195176  0.68940505  0.70685835
  0.72431164  0.74176493  0.75921822  0.77667152  0.79412481  0.8115781
  0.82903139  0.84648469  0.86393798  0.88139127  0.89884456  0.91629786
  0.93375115  0.95120444  0.96865773  0.98611103  1.00356432  1.02101761
  1.0384709   1.0559242   1.07337749  1.09083078  1.10828408  1.12573737
  1.14319066  1.16064395  1.17809725  1.19555054  1.21300383  1.23045712
  1.24791042  1.26536371  1.282817    1.30027029  1.31772359  1.33517688
  1.35263017  1.37008346  1.38753676  1.40499005  1.42244334  1.43989663
  1.45734993  1.47480322  1.49225651  1.5097098   1.5271631   1.54461639
  1.56206968  1.57952297  1.59697627  1.61442956  1.63188285  1.64933614
  1.66678944  1.68424273  1.70169602  1.71914931  1.73660261  1.7540559
  1.77150919  1.78896248  1.80641578  1.82386907  1.84132236  1.85877565
  1.87622895  1.89368224  1.91113553  1.92858882  1.94604212  1.96349541
  1.9809487   1.99840199  2.01585529  2.03330858  2.05076187  2.06821516
  2.08566846  2.10312175  2.12057504  2.13802833  2.15548163  2.17293492
  2.19038821  2.2078415   2.2252948   2.24274809  2.26020138  2.27765467
  2.29510797  2.31256126  2.33001455  2.34746784  2.36492114  2.38237443
  2.39982772  2.41728101  2.43473431  2.4521876   2.46964089  2.48709418
  2.50454748  2.52200077  2.53945406  2.55690735  2.57436065  2.59181394
  2.60926723  2.62672052  2.64417382  2.66162711  2.6790804   2.69653369
  2.71398699  2.73144028  2.74889357  2.76634686  2.78380016  2.80125345
  2.81870674  2.83616003  2.85361333  2.87106662  2.88851991  2.9059732
  2.9234265   2.94087979  2.95833308  2.97578637  2.99323967  3.01069296
  3.02814625  3.04559954  3.06305284  3.08050613  3.09795942  3.11541271]
width 0.017453292519943295
In [21]:
fig, ax = plt.subplots(1, 1, figsize=(FIG_WIDTH, FIG_WIDTH), subplot_kw=dict(projection='polar'))
for i, c in enumerate(['r', 'b']): 
    #ax.plot(theta, np.abs(code[:, i])**.5, color=c)
    #ax.plot(theta+np.pi, np.abs(code[:, i])**.5, color=c)
    # in polar coordinates, probability should be made proportional to
    # the surface, i.e. bars height to the square root of probability
    # see http://blog.invibe.net/posts/2014-12-09-polar-bar-plots.html
    ax.bar(theta_bins, (np.abs(code[:, i]))**.5, width=width, color='none', edgecolor=c, lw=2, align='edge', alpha=.8)# edgecolor="none")
    ax.bar(theta_bins+np.pi, (np.abs(code[:, i]))**.5, width=width, color='none', edgecolor=c, align='edge', alpha=.6)
    #ax.plot(theta_bin, np.ones_like(theta_bin)*np.sqrt(v_hist.mean()), 'r--')
    #ax.plot(theta_bin+np.pi, np.ones_like(theta_bin)*np.sqrt(v_hist.mean()), 'r--')
plt.setp(ax, yticks=[])
Out[21]:
[]

This shows the polar histogram of the energy of activation for symetric (red) and anti-symmetric (blue) filters.

This process is fast:

In [22]:
%%timeit
code = phi_ @ image_
523 µs ± 7.59 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

filtering video

We can now apply it to the whole video.

In [23]:
video_ = video.reshape((N_X*N_Y, N_frame))
print('video_.shape', video_.shape)
video_.shape (2916, 1024)
In [24]:
phi_ = phi.reshape((N_theta, N_phase, N_X*N_Y))
print('phi_.shape', phi_.shape)
phi_.shape (180, 2, 2916)
In [25]:
code = phi_ @ video_
print('All scalar products in one shot = ', (code).shape )
All scalar products in one shot =  (180, 2, 1024)
In [26]:
fig, axs = plt.subplots(2, 1, figsize=(FIG_WIDTH, FIG_WIDTH))

for i in range(2): 
    axs[i].pcolormesh(np.abs(code[:, i, :]), cmap=plt.viridis())
    if i: axs[i].set_xlabel('frame #')
    axs[i].set_ylabel(r'orientation $\theta$')
axs[0].set_title('symmetric')
axs[1].set_title('anti-symmetric')
Out[26]:
Text(0.5, 1.0, 'anti-symmetric')
In [27]:
%%timeit
code = phi_ @ video_
641 ms ± 34.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

temporal evolution for the energy

Following the formula $\cos(\theta)^2 + \sin(\theta)^2 = 1$, one can compute the total energy as the sum of the square of each component. This mimics complex cells.

In [28]:
code2 = np.sqrt(np.sum(code**2, axis=1))
print('code2.shape', code2.shape)
code2.shape (180, 1024)
In [29]:
fig, ax = plt.subplots(1, 1, figsize=FIGSIZE)

ax.pcolormesh(code2, cmap=plt.viridis())
ax.set_xlabel('frame #')
ax.set_ylabel(r'orientation $\theta$')
Out[29]:
Text(0, 0.5, 'orientation $\\theta$')

temporal evolution

Let's now show this temporal evolution as a dynamical polar plot.

In [30]:
theta_mean = np.sum(code2, axis=1)
print('theta_mean.shape', theta_mean.shape)
theta_mean.shape (180,)
In [31]:
fig, ax = plt.subplots(1, 1, figsize=(FIG_WIDTH, FIG_WIDTH), subplot_kw=dict(projection='polar'))
ax.bar(theta_bins, (theta_mean)**.5, width=width, color=c, align='edge', alpha=.5)
ax.bar(theta_bins+np.pi, (theta_mean)**.5, width=width, color=c, align='edge', alpha=.4)
ax.set_xticks([])
ax.set_yticks([])
#ax.grid(False)
ax.axis('off')
Out[31]:
(0.0, 6.283185307179586, 0.0, 44.6569249089471)
In [32]:
code2_mean = np.mean(code2, axis=0)
print('code2_mean.shape', code2_mean.shape)
code2_mean.shape (1024,)
In [33]:
fig, ax = plt.subplots(1, 1, figsize=(FIG_WIDTH, FIG_WIDTH))
ax.plot(code2_mean)
ax.plot(np.zeros_like(code2_mean), 'g--')
ax.set_xlabel('frame #')
ax.set_ylabel(r'energy in orientation $\theta$')
Out[33]:
Text(0, 0.5, 'energy in orientation $\\theta$')
In [34]:
fig, ax = plt.subplots(1, 1, figsize=FIGSIZE)

ax.pcolormesh(code2/code2_mean[None, :], cmap=plt.viridis())
ax.set_xlabel('frame #')
ax.set_ylabel(r'orientation $\theta$')
Out[34]:
Text(0, 0.5, 'orientation $\\theta$')
In [35]:
fig, ax = plt.subplots(1, 1, figsize=FIGSIZE)

ax.pcolormesh(code2/code2_mean[None, :]/theta_mean[:, None], cmap=plt.viridis())
ax.set_xlabel('frame #')
ax.set_ylabel(r'orientation $\theta$')
Out[35]:
Text(0, 0.5, 'orientation $\\theta$')