• simulations avec différents MCs
• ring avec représentation de la vitesse angulaire
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt

import numpy as np
import MotionClouds as mc

fig_width_pt = 646.79  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inches
fig_width = fig_width_pt*inches_per_pt  # width in inches

In [2]:
sf_0 = 0.15
B_sf = 0.05
theta_0 = np.pi/2
B_theta = 0.15
loggabor = True

def get_grids(N_X, N_Y):
fx, fy = np.mgrid[(-N_X//2):((N_X-1)//2 + 1), (-N_Y//2):((N_Y-1)//2 + 1)]
fx, fy = fx*1./N_X, fy*1./N_Y
return fx, fy

N_X, N_Y = fx.shape[0], fy.shape[1]
R2 = fx**2 + fy**2
R2[N_X//2 , N_Y//2] = np.inf
return np.sqrt(R2)

def envelope_orientation(fx, fy, theta_0=theta_0, B_theta=B_theta, norm = True):
theta = np.arctan2(fx, fy)
env =  np.exp(np.cos(2*(theta-theta_0))/B_theta**2)
#return (1/(2*np.pi*iv(0, (1/((B_theta)**2)))))*env
if norm: env /= np.sqrt((env**2).sum())
return env

def envelope_radial(fx, fy, sf_0=sf_0, B_sf=B_sf, loggabor=loggabor, norm = True):
if sf_0 == 0.: return 1.
if loggabor:
env = 1./fr*np.exp(-.5*(np.log(fr/sf_0)**2)/(np.log((sf_0+B_sf)/sf_0)**2))
if norm: env /= np.sqrt((env**2).sum())
return env
else:

fx, fy = get_grids(mc.N_X, mc.N_Y)

In [27]:
env = envelope_radial(fx, fy)

fig, ax = plt.subplots(1, 4, figsize=(14, 9))
for i, (f, label) in enumerate(zip([fx, fy, fr, env], ['x', 'y', 'f_r', 'env'])):
ax[i].matshow(f)
ax[i].set_title(label)
ax[i].set_xticks([])
ax[i].set_yticks([])


In [28]:
env = envelope_radial(fx, fy) * envelope_orientation(fx, fy)

fig, ax = plt.subplots(1, 4, figsize=(14, 9))
for i, (f, label) in enumerate(zip([fx, fy, fr, env], ['x', 'y', 'f_r', 'env'])):
ax[i].matshow(f)
ax[i].set_title(label)
ax[i].set_xticks([])
ax[i].set_yticks([])



## Figure¶

In [5]:
#  réalisation du motion clouds

theta0 = np.pi/2
Btheta = 0.15
theta_0 = [0, np.pi/4, np.pi/2]
B_theta = [0.1, 0.5, 1.]

def texture(env):
return np.fft.fft2(np.fft.ifftshift(env * np.exp(1j * 2 * np.pi * np.random.rand(mc.N_X, mc.N_Y)))).real

fig, ax = plt.subplots(1, 3, figsize=(fig_width, fig_width*6/18))
for i, (theta0_, label) in enumerate(zip(theta_0, [r'$\theta = \pi$', r'$\theta = \pi/4$', r'$\theta = \pi/2$']) ) :
env = envelope_radial(fx, fy) * envelope_orientation(fx, fy, theta_0=theta0_, B_theta=Btheta)
I = texture(env)
ax[i].matshow(I, cmap=plt.cm.gray)
ax[i].set_title(label)
#    ax[i].set_xticks(np.linspace(0, 120, 5))
#    ax[i].set_yticks(np.linspace(0, 120, 5))
ax[i].set_xticks([])
ax[i].set_yticks([])
#plt.title(u'Réalisation du Motion Clouds')
plt.tight_layout()
fig.savefig('figures/realisation_MC_theta0(week4).pdf')

fig, ax = plt.subplots(1, 3, figsize=(fig_width, fig_width*6/18))
for i, (Btheta_, label) in enumerate(zip(B_theta, [r'$B_\theta = 0.1$', r'$B_\theta = 0.5$', r'$B_\theta = 1.0$']) ) :
env = envelope_radial(fx, fy) * envelope_orientation(fx, fy, theta_0=theta0, B_theta=Btheta_)
I = texture(env)
ax[i].matshow(I, cmap=plt.cm.gray)
ax[i].set_title(label)
ax[i].set_xticks([])
ax[i].set_yticks([])
#plt.title(u'Réalisation du Motion Clouds')
plt.tight_layout()
fig.savefig('figures/realisation_MC_Btheta(week4).pdf')


### Réponse théorique de l'énergie du filtre linéaire¶

In [6]:
# prédiction de la réponse spikante maximale

env = envelope_radial(fx, fy) * envelope_orientation(fx, fy)

reponse = env**2

fig, ax = plt.subplots(figsize=(14, 9))
ax.matshow(reponse)
plt.title(u'Réponse spikante')

Out[6]:
<matplotlib.text.Text at 0x110738f90>
In [7]:
# montrons que la convolution d'un MC avec un gabor donne un von - Mises
sf_0 = 0.15
B_sf = 0.05
theta_0 = np.pi/2
B_theta = 0.15
loggabor = True
# testons d'abord la convolution avec un angle arbitraire
seed = None
np.random.seed(seed=seed)
def texture(env):
I= np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(env * np.exp(1j * 2 * np.pi * np.random.rand(mc.N_X, mc.N_Y)))).real)
#I /= np.sqrt((env**2).sum())
I /= env.sum()
return I

def impulse(env, phi=2 * np.pi):
I = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(env * np.exp(1j * phi))).real)
#I /= np.sqrt((env**2).sum())
I /= env.sum()
return I

env_in = envelope_radial(fx, fy, sf_0=sf_0, B_sf=B_sf) * envelope_orientation(fx, fy, theta_0=theta_0, B_theta=B_theta)
env_V1 = envelope_radial(fx, fy, sf_0=sf_0, B_sf=B_sf) * envelope_orientation(fx, fy, theta_0=np.random.rand()*np.pi, B_theta=B_theta)

fig, ax = plt.subplots(1, 2, figsize=(14, 9))
for i, (f, label) in enumerate(zip([texture(env_in), impulse(env_V1)], [u'Texture', u'Gabor'])):
ax[i].matshow(f, cmap=plt.cm.gray)
ax[i].set_title(label)
ax[i].set_xticks([])
ax[i].set_yticks([])


Pour justifier la convolution circulaire que nous allons utiliser voir :

http://stackoverflow.com/questions/18172653/convolving-a-periodic-image-with-python

In [8]:
#from scipy.signal import fftconvolve as convolve
# convolution circulaire
def convolve(image_in, image_V1):
env_in = np.fft.fft2(image_in)
env_V1 = np.fft.fft2(image_V1)
return np.fft.fftshift(np.fft.ifft2((env_in*env_V1)).real)

R = convolve(texture(env_in), impulse(env_V1))
print R.shape
fig, ax = plt.subplots(figsize=(14, 9))
ax.matshow(R, cmap=plt.cm.gray)
plt.title(u"Convolution d'un MC avec un gabor")
print (R**2).sum() / mc.N_X**4 # à normaliser

(128, 128)
6.5244606908e-46

In [9]:
# images des convolutions avec différents angles
N_theta=360
theta0 = np.pi/2
theta_0 = np.linspace(0., np.pi, 10)

for i, theta0_ in enumerate(theta_0) :
env_in = envelope_radial(fx, fy, sf_0=sf_0, B_sf=B_sf) * envelope_orientation(fx, fy, theta_0=theta0, B_theta=B_theta)
env_V1 = envelope_radial(fx, fy, sf_0=sf_0, B_sf=B_sf) * envelope_orientation(fx, fy, theta_0=theta0_, B_theta=B_theta)
R = convolve(texture(env_in), impulse(env_V1))
fig, ax = plt.subplots(1, 3, figsize=(14, 9))
for i, (f, label) in enumerate(zip([texture(env_in), impulse(env_V1), R], [u'Texture', u'Gabor', u'convolution'])):
ax[i].matshow(f, cmap=plt.cm.gray)
ax[i].set_title(label)
ax[i].set_xticks([])
ax[i].set_yticks([])