2014-06-18 [Stage M1] Chloé Pasturel : week 4

  • 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

def frequency_radius(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:
        fr = frequency_radius(fx, fy)
        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:
        return np.exp(-.5*(frequency_radius(fx, fy) - sf_0)**2/B_sf**2)

fx, fy = get_grids(mc.N_X, mc.N_Y)
fr = frequency_radius(fx, fy)
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([])