# Speed distributions

In [1]:
%matplotlib inline
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
#!rm -fr ../files/speed*

In [2]:
import MotionClouds as mc
name = 'noisy-speed'
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
print(mc.envelope_speed.__doc__)

    Returns the speed envelope:
selects the plane corresponding to the speed (V_X, V_Y) with some bandwidth B_V.

* (V_X, V_Y) = (0,1) is downward and  (V_X, V_Y) = (1, 0) is rightward in the movie.
* A speed of V_X=1 corresponds to an average displacement of 1/N_X per frame.
To achieve one spatial period in one temporal period, you should scale by
V_scale = N_X/float(N_frame)
If N_X=N_Y=N_frame and V=1, then it is one spatial period in one temporal
period. It can be seen along the diagonal in the fx-ft face of the MC cube.

A special case is used when B_V=0, where the fx-ft plane is used as
the speed plane: in that case it is desirable to set (V_X, V_Y) to (0, 0)
to avoid aliasing problems.

Run the 'test_speed' notebook to explore the speed parameters, see
http://motionclouds.invibe.net/posts/testing-speed.html


In [3]:
# explore parameters
for B_V in [0.0, 0.01, 0.1, 1.0, 10.0]:
name_ = name + '-B_V-' + str(B_V).replace('.', '_')
z = mc.envelope_gabor(fx, fy, ft, V_X=0, B_V=B_V)
mc.figures(z, name_)
mc.in_show_video(name_)


# More is not always better

## MotionClouds¶

MotionClouds are random dynamic stimuli optimized to study motion perception.

Notably, this method was used in the following paper:

• Claudio Simoncini, Laurent U. Perrinet, Anna Montagnini, Pascal Mamassian, Guillaume S. Masson. More is not always better: dissociation between perception and action explained by adaptive gain control. Nature Neuroscience, 2012 URL

In this notebook, we describe the scripts used to generate such stimuli.

# Colored Motion Clouds

## Exploring colored Motion Clouds¶

By construction, Motion Clouds are grayscale:

In [1]:
import os
import numpy as np
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
name = 'color'

In [2]:
env = mc.envelope_gabor(fx, fy, ft, V_X=0., V_Y=0.)
mc.figures(env, name + '_gray', do_figs=False)
mc.in_show_video(name + '_gray')


However it is not hard to imagine extending them to the color space. A first option is to create a diferent motion cloud for each channel:

In [3]:
colored_MC = np.zeros((env.shape[0], env.shape[1], 3, env.shape[2]))

for i in range(3):
colored_MC[:, :, i, :] = mc.rectif(mc.random_cloud(mc.envelope_gabor(fx, fy, ft, V_X=0., V_Y=0.)))

mc.anim_save(colored_MC, os.path.join(mc.figpath, name + '_color'))
mc.in_show_video(name + '_color')


Note that the average luminance is also a random cloud:

In [4]:
mc.anim_save(mc.rectif(colored_MC.sum(axis=2)), os.path.join(mc.figpath, name + '_gray2'))
mc.in_show_video(name + '_gray2')


We may create a strictly isoluminant cloud:

In [5]:
luminance = colored_MC.sum(axis=2)[:, :, np.newaxis, :]
mc.anim_save(colored_MC/luminance, os.path.join(mc.figpath, name + '_isocolor'))
mc.in_show_video(name + '_isocolor')


There are now many more possibilities, such as

• weighting the different R, G and B channels to obtain a better tuning to psychophysics,
• let only the predominant channels be activated (like the one corresponding to the red and green channels which correspond to the maximal responses of cones).

# Testing the utility functions of Motion Clouds

## Motion Clouds utilities¶

Here, we test some of the utilities that are delivered with the MotionClouds package.

In [1]:
%load_ext autoreload

In [2]:
import os
import MotionClouds as mc

In [3]:
mc.N_X, mc.N_Y, mc.N_frame = 30, 40, 50
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)


### generating figures¶

As they are visual stimuli, the main outcome of the scripts are figures. Utilities allow to plot all figures, usually marked by a name:

In [4]:
name = 'testing_utilities'

In [5]:
help(mc.figures_MC)

Help on function figures_MC in module MotionClouds:

figures_MC(fx, fy, ft, name, V_X=1.0, V_Y=0.0, do_figs=True, do_movie=True, B_V=0.5, sf_0=0.125, B_sf=0.1, loggabor=True, recompute=False, theta=0.0, B_theta=0.19634954084936207, alpha=0.0, vext='.mp4', seed=None, impulse=False, do_amp=False, verbose=False, figpath='../files/', return_envelope=False, **kwargs)
Generates the figures corresponding to the Fourier spectra and the stimulus cubes and
movies directly from the parameters.

The figures names are automatically generated.


In [6]:
mc.figures_MC(fx, fy, ft, name, recompute=True)

In [7]:
help(mc.in_show_video)

Help on function in_show_video in module MotionClouds:

in_show_video(name, vext='.mp4', loop=True, autoplay=True, controls=True, embed=False, figpath='../files/', **kwargs)
Columns represent isometric projections of a cube. The left column displays
iso-surfaces of the spectral envelope by displaying enclosing volumes at 5
different energy values with respect to the peak amplitude of the Fourier spectrum.
The middle column shows an isometric view of the faces of the movie cube.
The first frame of the movie lies on the x-y plane, the x-t plane lies on the
top face and motion direction is seen as diagonal lines on this face (vertical
motion is similarly see in the y-t face). The third column displays the actual
movie as an animation.

Given a name, displays the figures corresponding to the Fourier spectra, the
stimulus cubes and movies within the notebook.


In [8]:
mc.in_show_video(name)


This function embeds the images and video within the notebook. Sometimes you want to avoid that:

In [9]:
mc.in_show_video(name, embed=False)


Sometimes, you may have already computed some envelope or just want to distort it, then you can use mc.figures:

In [10]:
env = mc.envelope_gabor(fx, fy, ft)

In [11]:
help(mc.figures)

Help on function figures in module MotionClouds:

figures(z=None, name='MC', vext='.mp4', do_movie=True, do_figs=True, recompute=False, seed=None, impulse=False, verbose=False, masking=False, do_amp=False, figpath='../files/', **kwargs)
Given an envelope, generates the figures corresponding to the Fourier spectra
and the stimulus cubes and movies.

The figures names are automatically generated.


In [12]:
import numpy as np
mc.figures(np.sqrt(env), name + '_0')

In [13]:
mc.in_show_video(name + '_0')


### low-level figures : 3D visualizations¶

In [14]:
help(mc.cube)

Help on function cube in module MotionClouds:

cube(im_in, azimuth=30.0, elevation=45.0, name=None, ext='.png', do_axis=True, show_label=True, cube_label={'x': 'x', 'y': 'y', 't': 't'}, colormap='gray', roll=-180.0, vmin=0.0, vmax=1.0, figsize=(800, 800), figpath='../files/', **kwargs)
Visualization of the stimulus as a cube


In [15]:
help (mc.visualize)

Help on function visualize in module MotionClouds:

visualize(z_in, azimuth=25.0, elevation=30.0, thresholds=[0.94, 0.89, 0.75, 0.5, 0.25, 0.1], opacities=[0.9, 0.8, 0.7, 0.5, 0.2, 0.1], fourier_label={'f_x': 'f_x', 'f_y': 'f_y', 'f_t': 'f_t'}, name=None, ext='.png', do_axis=True, do_grids=False, draw_projections=True, colorbar=False, f_N=2.0, f_tN=2.0, figsize=(800, 800), figpath='../files/', **kwargs)
Visualization of the Fourier spectrum by showing 3D contour plots at different thresholds

parameters
----------
z : envelope of the cloud



### Handling filenames¶

By default, the folder for generating figures or data is mc.figpath:

In [16]:
print(mc.figpath)

../files/

print(os.listdir(mc.figpath))

To generate figures, we assign file names, such as:

In [17]:
filename = os.path.join(mc.figpath, name)


It is then possible to check if that figures exist:

In [18]:
print('filename=', filename, ', exists? : ', mc.check_if_anim_exist(filename))

filename= ../files/testing_utilities , exists? :  False


Note that the file won't be recomputed if it exists:

In [19]:
mc.figures(env, name)


This behavior can be overriden using the recompute option

In [20]:
mc.figures(env, name, recompute=True)


Warning: be sure that when you display a given file, it corresponds to the parameters you have set for your stimulus.

### low-level figures : exporting to various formats¶

It is possible to export motion clouds to many different formats. Here are some examples:

In [21]:
!rm -fr ../files/export

In [22]:
name = 'export'
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
z = mc.rectif(mc.random_cloud(mc.envelope_gabor(fx, fy, ft)))
mc.PROGRESS = False
for vext in mc.SUPPORTED_FORMATS:
print ('Exporting to format: ', vext)
mc.anim_save(z, os.path.join(mc.figpath, name), display=False, vext=vext, verbose=False)

Exporting to format:  .h5
Exporting to format:  .mpg
Exporting to format:  .mp4
Exporting to format:  .gif
Exporting to format:  .webm
Exporting to format:  .zip
Exporting to format:  .mat
Exporting to format:  .png


### showing a video¶

To show a video in a notebook, issue:

In [23]:
mc.notebook = True # True by default
mc.in_show_video('export')


### Rectifying the contrast¶

The mc.rectif function allows to rectify the amplitude of luminance values within the whole generated texture between $0$ and $1$:

In [24]:
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
envelope = mc.envelope_gabor(fx, fy, ft)
image = mc.random_cloud(envelope)
print('Min :', image.min(), ', mean: ', image.mean(), ', max: ', image.max())

Min : -2.21666797481 , mean:  -1.70234197109e-19 , max:  2.09408478901

In [25]:
image = mc.rectif(image)
print('Min :', image.min(), ', mean: ', image.mean(), ', max: ', image.max())

Min : 0.0 , mean:  0.5 , max:  0.972349673655

In [26]:
import pylab
import numpy as np
import matplotlib.pyplot as plt
import math
%matplotlib inline
#%config InlineBackend.figure_format='retina' # high-def PNGs, quite bad when using file versioning
%config InlineBackend.figure_format='svg'

In [27]:
name = 'contrast_methods-'
#initialize
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
ext = '.zip'
contrast = 0.25
B_sf = 0.3

for method in ['Michelson', 'Energy']:
z = mc.envelope_gabor(fx, fy, ft, B_sf=B_sf)
im = np.ravel(mc.random_cloud(z, seed =1234))
im_norm = mc.rectif(mc.random_cloud(z), contrast, method=method, verbose=True)
plt.figure()
plt.subplot(111)
plt.title(method + ' Histogram Ctr: ' + str(contrast))
plt.ylabel('pixel counts')
plt.xlabel('grayscale')
bins = int((np.max(im_norm[:])-np.min(im_norm[:])) * 256)
plt.xlim([0, 1])
plt.hist(np.ravel(im_norm), bins=bins, normed=False, facecolor='blue', alpha=0.75)
#plt.savefig(name_)

def image_entropy(img):
"""calculate the entropy of an image"""
histogram = img.histogram()
histogram_length = np.sum(histogram)
samples_probability = [float(h) / histogram_length for h in histogram]
return -np.sum([p * math.log(p, 2) for p in samples_probability if p != 0])

Before Rectification of the frames
Mean= 1.06581410364e-18 , std= 0.84527878442 , Min= -3.07067785094 , Max= 3.47591337022  Abs(Max)= 3.47591337022
After Rectification of the frames
Mean= 0.5 , std= 0.0303977219219 , Min= 0.389572986872 , Max= 0.625
percentage pixels clipped= 0.0
Before Rectification of the frames
Mean= -1.18423789293e-19 , std= 0.817871986507 , Min= -3.42450006883 , Max= 3.40014935729  Abs(Max)= 3.42450006883
After Rectification of the frames
Mean= 0.5 , std= 0.125 , Min= -0.0233857078689 , Max= 1.01966405094
percentage pixels clipped= 0.005


If we normalise the histogram then the entropy base on gray levels is going to be the almost the same.
TODO: Review the idea of entropy between narrowband and broadband stimuli.

# Testing basic functions of Motion Clouds

## Motion Clouds: raw principles¶

Motion Clouds are synthesized textures which aim at having similar characteristics as natural images but with controlled parameters. There are many ways to achieve these results and this notebook aims at showing that different procedures from different communities (neurioscience, modelling, computer vision, ...) may produce similar results.

In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline


### Using Fourier ("official Motion Clouds")¶

In [2]:
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)


### Using mixtures of images¶

In [3]:
from scipy.misc import face
lena = face()
print(lena.shape)
lena = lena[:, (1024-768):, :].mean(axis=-1)
print(lena.shape)
lena -= lena.mean()
lena /= lena.std()
print(lena.shape)

(768, 1024, 3)
(768, 768)
(768, 768)

In [4]:
plt.imshow(lena, cmap=plt.cm.gray)

Out[4]:
<matplotlib.image.AxesImage at 0x103200d68>
In [5]:
lena.shape

Out[5]:
(768, 768)
In [6]:
lena[0, :]

Out[6]:
array([ 0.723,  0.824,  1.025,  1.132,  0.895,  0.45 ,  0.159,  0.1  ,
-0.185, -0.292, -0.392, -0.339, -0.143,  0.124,  0.343,  0.45 ,
0.907,  1.405,  1.529,  1.138,  0.444, -0.446, -0.908, -0.736,
0.011,  0.367,  0.402, -0.114, -0.47 , -0.238,  0.26 ,  0.58 ,
0.883,  0.491,  0.064, -0.096,  0.011,  0.171,  0.171,  0.1  ,
0.474,  0.349,  0.082, -0.292, -0.665, -0.932, -1.075, -1.11 ,
-0.499, -0.214,  0.106,  0.207,  0.064, -0.131, -0.22 , -0.22 ,
-0.315, -0.327, -0.339, -0.297, -0.22 , -0.161, -0.137, -0.161,
0.07 , -0.048, -0.137, -0.226, -0.155,  0.201,  0.438,  0.367,
0.248,  0.296,  0.402,  0.539,  0.646,  0.663,  0.58 ,  0.527,
0.527,  0.402,  0.136, -0.025,  0.136,  0.385,  0.26 , -0.06 ,
-0.007, -0.042,  0.171,  0.325, -0.066, -0.742, -0.956, -0.725,
-0.558, -0.416, -0.452, -0.63 , -0.665, -0.505, -0.416, -0.517,
-0.713, -0.73 , -0.819, -0.926, -0.897, -0.879, -1.039, -1.252,
-1.235, -0.505,  0.189,  0.029, -0.576, -0.713, -0.375, -0.143,
-0.197, -0.398, -0.487, -0.452, -0.398, -0.22 , -0.203, -0.381,
-0.096, -0.078, -0.042, -0.203, -0.47 , -0.487, -0.452, -0.558,
-0.416, -0.381, -0.416, -0.363, -0.007,  0.313,  0.064, -0.434,
-0.529, -0.458, -0.386, -0.28 , -0.173, -0.155, -0.297, -0.475,
-0.262, -0.369, -0.44 , -0.386, -0.262, -0.084,  0.183,  0.45 ,
0.242, -0.09 , -0.321, -0.321, -0.339, -0.464, -0.446, -0.274,
-0.114,  0.242,  0.58 ,  0.628,  0.432,  0.219,  0.076,  0.005,
-0.072, -0.41 , -0.819, -0.98 , -0.819, -0.535, -0.41 , -0.41 ,
-0.161,  0.319,  0.462, -0.09 , -0.784, -0.908, -0.428,  0.07 ,
0.165,  0.076, -0.037, -0.037, -0.066, -0.226, -0.042,  0.438,
0.367,  0.171, -0.137, -0.155,  0.219,  0.598,  0.918,  1.239,
1.156,  0.978,  0.711,  0.497,  0.444,  0.462,  0.402,  0.26 ,
0.242,  0.1  , -0.149, -0.523, -0.914, -0.986, -0.487,  0.118,
0.141,  0.497,  0.681,  0.752,  0.764,  0.414,  0.112,  0.195,
0.195,  0.064,  0.005, -0.114, -0.274, -0.173, -0.001,  0.017,
0.094,  0.165,  0.219,  0.13 , -0.066, -0.155, -0.013,  0.219,
0.521,  0.29 ,  0.094,  0.183,  0.503,  0.752,  0.752,  0.628,
-0.019, -0.001,  0.248,  0.782,  1.227,  1.227,  0.764,  0.296,
-0.096,  0.242,  0.847,  1.452,  1.595,  1.221,  0.669,  0.313,
0.646,  1.108,  1.28 ,  0.96 ,  0.569,  0.396,  0.628,  1.037,
1.405,  1.28 ,  0.77 ,  0.278,  0.207,  0.248,  0.201,  0.183,
-0.084, -0.031,  0.379,  0.788,  0.966,  1.126,  0.913,  0.308,
0.361,  0.895,  1.268,  1.322,  1.268,  0.984,  0.61 ,  0.432,
0.527,  0.491,  0.343,  0.213,  0.302,  0.592,  0.776,  0.829,
0.604,  0.337,  0.159,  0.373,  0.776,  0.936,  0.776,  0.527,
0.426,  0.551,  0.853,  1.162,  1.179,  0.841,  0.48 ,  0.284,
0.622,  0.408,  0.296,  0.509,  0.889,  1.12 ,  1.031,  0.853,
0.515,  0.497,  0.337, -0.001, -0.297, -0.286, -0.037,  0.224,
0.883,  0.966,  1.031,  1.049,  0.996,  0.705,  0.13 , -0.422,
-0.019,  0.017, -0.001, -0.072, -0.019,  0.195,  0.355,  0.408,
0.462,  0.764,  0.996,  0.907,  0.675,  0.515,  0.462,  0.444,
0.254,  0.094, -0.084, -0.143, -0.072, -0.06 , -0.185, -0.363,
0.183,  0.521,  0.835,  0.907,  0.764,  0.515,  0.195, -0.037,
-0.096, -0.143, -0.131,  0.064,  0.396,  0.669,  0.628,  0.468,
0.355,  0.355,  0.408,  0.491,  0.438,  0.26 ,  0.094,  0.023,
0.011,  0.1  ,  0.1  , -0.001,  0.035,  0.189,  0.29 ,  0.254,
0.195,  0.213,  0.236,  0.254,  0.337,  0.408,  0.373,  0.313,
0.118, -0.007, -0.078, -0.025,  0.023,  0.041,  0.094,  0.195,
0.23 ,  0.213,  0.195,  0.195,  0.177,  0.118,  0.029, -0.06 ,
-0.078, -0.096, -0.078,  0.047,  0.159,  0.106, -0.09 , -0.274,
-0.452, -0.167,  0.136,  0.224,  0.171,  0.183,  0.219,  0.236,
-0.108, -0.09 , -0.037,  0.052,  0.141,  0.177,  0.177,  0.159,
-0.037, -0.179, -0.25 , -0.161, -0.09 , -0.125, -0.25 , -0.321,
-0.428, -0.357, -0.143,  0.047, -0.078, -0.363, -0.511, -0.493,
-0.381, -0.47 , -0.517, -0.428, -0.292, -0.238, -0.327, -0.452,
-0.452, -0.434, -0.416, -0.434, -0.434, -0.434, -0.398, -0.363,
-0.363, -0.381, -0.398, -0.47 , -0.541, -0.576, -0.612, -0.612,
-0.547, -0.529, -0.511, -0.493, -0.493, -0.511, -0.529, -0.529,
-0.653, -0.636, -0.618, -0.618, -0.636, -0.618, -0.582, -0.547,
-0.594, -0.63 , -0.689, -0.725, -0.742, -0.76 , -0.796, -0.814,
-0.814, -0.814, -0.814, -0.831, -0.861, -0.879, -0.897, -0.914,
-0.903, -0.92 , -0.938, -0.956, -0.938, -0.903, -0.849, -0.814,
-0.938, -0.938, -0.956, -0.974, -0.991, -1.009, -1.027, -1.045,
-1.009, -0.956, -0.831, -0.725, -0.653, -0.677, -0.748, -0.819,
-0.849, -0.849, -0.855, -0.837, -0.819, -0.766, -0.736, -0.695,
-0.855, -0.873, -0.873, -0.891, -0.867, -0.843, -0.796, -0.766,
-0.742, -0.73 , -0.683, -0.677, -0.671, -0.671, -0.647, -0.63 ,
-0.606, -0.564, -0.523, -0.446, -0.381, -0.345, -0.327, -0.309,
-0.333, -0.404, -0.517, -0.612, -0.683, -0.677, -0.636, -0.606,
-0.736, -0.689, -0.606, -0.529, -0.464, -0.41 , -0.327, -0.286,
-0.268, -0.179, -0.072,  0.029,  0.082,  0.165,  0.254,  0.349,
0.396,  0.361,  0.325,  0.325,  0.313,  0.284,  0.213,  0.147,
0.13 ,  0.141,  0.124,  0.106,  0.082,  0.082,  0.124,  0.195,
0.106, -0.535, -1.229, -1.478, -1.324, -1.027, -0.908, -0.956,
-0.855, -0.962, -1.152, -1.312, -1.407, -1.472, -1.49 , -1.508,
-1.407, -1.508, -1.596, -1.543, -1.43 , -1.359, -1.401, -1.478,
-1.223, -1.68 , -1.834, -1.531, -0.481,  0.468,  0.764,  0.646,
0.379, -0.297, -1.158, -1.561, -1.436, -1.098, -1.015, -1.14 ,
-1.057, -1.051, -1.051, -1.051, -1.063, -1.08 , -1.098, -1.098,
-1.039, -1.069, -1.069, -1.051, -1.027, -1.021, -1.027, -1.045,
-0.956, -0.956, -0.962, -0.944, -0.938, -0.932, -0.95 , -0.938,
-1.069, -1.08 , -1.08 , -1.092, -1.098, -1.098, -1.086, -1.086,
-1.158, -1.134, -1.128, -1.116, -1.104, -1.14 , -1.146, -1.163,
-1.152, -1.169, -1.205, -1.223, -1.223, -1.169, -1.134, -1.11 ,
-1.027, -0.974, -0.903, -0.867, -0.867, -0.903, -0.92 , -0.92 ,
-1.039, -0.997, -0.997, -1.075, -1.033, -0.873, -0.671, -0.576,
-0.309, -0.256, -0.292, -0.392, -0.475, -0.47 , -0.47 , -0.523,
-0.475, -0.458, -0.44 , -0.404, -0.327, -0.173, -0.007,  0.106,
0.171, -0.025, -0.191, -0.179,  0.035,  0.367,  0.752,  1.037,
0.954,  0.99 ,  0.972,  0.883,  0.901,  0.942,  0.913,  0.806,
0.752,  0.58 ,  0.302,  0.052, -0.054, -0.001,  0.177,  0.319])
In [7]:
def noise(image=lena):
for axis in [0, 1]:
image = np.roll(image, np.random.randint(image.shape[axis]), axis=axis)
return image

In [8]:
plt.imshow(noise(), cmap=plt.cm.gray)

Out[8]:
<matplotlib.image.AxesImage at 0x109532cc0>
In [9]:
plt.imshow(noise(), cmap=plt.cm.gray)

Out[9]:
<matplotlib.image.AxesImage at 0x109f67b38>

### Using ARMA processes¶

Now, we define the ARMA process as an averaging process with a certain time constant $\tau=30.$ (in frames).

In [10]:
def ARMA(image, tau=30.):
image = (1 - 1/tau)* image + 1/tau * noise()
return image


initializing

In [11]:
image = ARMA(lena)
plt.imshow(image, cmap=plt.cm.gray)

Out[11]:
<matplotlib.image.AxesImage at 0x10cea6390>