[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([])
    
No description has been provided for this image
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([])
    
No description has been provided for this image

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')
No description has been provided for this image
No description has been provided for this image

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>
No description has been provided for this image
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([])    
No description has been provided for this image

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
No description has been provided for this image
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([])
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image