Instructions for Psychophysicists

making experiments using motionclouds

making psychophysics using psychopy

installation (preparation 3min, cooking / downloading 10 min)

  • it is easy to create a simple psychophysical experiment using a combination of PsychoPy and MotionClouds. Simply download the standalone psychopy binary for your computer and try to load the example script psychopy_competing.py.
  • for this, you may:
    1. clone / download the MC repository on github.
    2. download PsychoPy
  • then, open the psychopy_competing.py file in the editor (the so-called "PsychoPy Coder")

running

  • Hit the run button (the green round button with a runner on it), it will then:
  • tell you about some configuration parameters (be sure to enter the correct size of your screen)
  • pre-generate data
  • run the experiment, that is present a sequence of trials consisting of :
    1. a movie presentation,
    2. a blank screen where the program waits for a key press: UP arrow if you saw the movie going up, DOWN if down or Q / Esc to exit

results

  • the results are displayed at the end showing that :
  • the movies consisted of 2 MotionClouds in opposite directions (UP and DOWN) but added with different contrasts,
  • the percept was going up or down (on the y-axis) as a function of this contrast (x-axis)
  • it seems I have a slight bias for UP movements as the contrast of equal probability seems inferior to 1/2 (a proper analysis would necessitate a fit).

TODO: Frequently asked questions

Experimental protocol

These instructions were given on VSG system but should work on similar working platforms:

  1. Measure your experimental settings (screen size, frame rate, resolution, viewing distance) as these values will constraint your spatial and temporal sampling frequencies.

  2. Calibrate your screen (gamma / luminance).

  3. Experimental constants (with some default values):

    a. $D$ : Viewing distance distance = 570 mm

    b. $S$: screen size = 390 mm

    c. $dt$ : Max frame period 0,01 s

    d. $\text{VA}$ : Stimulus Width in degrees of visual angle at viewing distance $D$: 38,0958

    e. $f_{d}$ : Desired spatial frequency in cycles per degree - cpd

    f. $V_{drift}$: Drifting velocity in degrees per second - d/s

    g. $X$, $Y$ : sampling step in degrees.

    h. $F_X$, $F_Y$: spatial sampling rates in cpd

    i. $T$ : Stimulus duration in seconds.

    j. $\alpha$: mean orientation of the frequency component (called $\theta$ in the scripts)

    k. $\theta$: direction of image motion.

  4. Parameters for the stimulus:

    a. ($N_X$, $N_Y$) : Stimuli’s frame size in pixels as generated in Fourier space.

    b. ($n_X$, $n_Y$): Screen resolution = 640 x 480 px (approx 0.06 deg/px)

    c. ($N_frame$): Number of movie frames as generated in Fourier space.

    d. $f_pref$ : Central frequency normalized in Fourier space

    e. $V_X$: left - rightward motion

    f. $V_Y$: up-downward motion

    g. $W_X$, $W_Y$: spatial sampling frequencies in Fourier space

How do I compute a visual angle?

$$VA = 2*arctan\left( \frac{S}{2D} \right)$$

where $S$ is stimulus size on the screen (X or Y) and $D$ is the viewing distance.

How do I set the spatial sampling?

$X = \frac{\text{VA}}{n_{x}}$ is the spatial step in the x-axis and for the y-axis the sampling step, Y, has the same value the x-axis.

In practice, the spatial sampling step (X and/or Y) can not be larger than the semi-period of the finest spatial frequency ($F_X$ and or $F_Y$) that is represented in the image. The image will not contain those frequency components are higher than the Nyquist frequency along the two axes.

What is the unit for frequencies?

The frequency distribution is periodic with period $F_S$. The normalized frequency ranges from 0.0 to 1.0, which corresponds to a real frequency range of 0 to the sampling frequency $F_S$.

The normalized frequency also wraps around 1.0 so a normalized frequency of 1.1 is equivalent to 0.1.

When the actual frequency, $f_t$, has units of Hz, for example, the normalized frequencies, also denoted by $f_t$, have units of cycles per sample, and the periodicity of the normalized distribution is 1.

What is the maximal frequency ?

If a sampled signal is real-valued the periodicity of the frequency distribution is still $F_S$. But due to symmetry, it is completely defined by the content within a span of just $F_S/2$ (Nyquist frequency). Accordingly, some applications use that as the normalization reference (and the resulting units are half-cycles per sample).

Normalization produces a distribution that is independent of the sample-rate, and thus one plot is sufficient for all possible sample-rates.

fx, fy, ft = np.mgrid[(-N\_X//2):((N\_X-1)//2 + 1),
(-N\_Y//2):((N\_Y-1)//2 + 1),(-N\_frame//2):((N\_frame-1)//2 + 1)]

fx, fy, ft = fx\*1./N\_X, fy\*1./N\_Y, ft\*1./N\_frame

The range is always $(-0.5, 0.5)$ which corresponds to the maximum frequency we can represent.

Example: N_X = N_Y = 512 and N_frames = 128 then,

Nyquist limit is $2^{- 1}W\ cyc/width$, W is the width of the stimulus.

The highest frequency band considered will be centered at $2^{- 2}W\ cyc/width$

The blobs change spatial frequency in octave steps, the scale-k band will have a spatial frequency pass band centered at: $2^{- (k + 2)}W\ cyc/width$

easy: it's just a design choice, we could have used B_tf, but defining B_sf then slicing in the envelope with a B_V envelope was more intuitive.

Convention used in DFT

From numpy.fft.ifftn(a, s=None, axes=None)

From http://docs.scipy.org/doc/numpy/reference/routines.fft.html:

The values in the result follow so-called “standard” order: If A = fft(a,n), then A[0] contains the zero-frequency term (the mean of the signal), which is always purely real for real inputs. Then A[1:n/2] contains the positive-frequency terms, and A[n/2+1:] contains the negative-frequency terms, in order of decreasingly negative frequency.

For an even number of input points, A[n/2] represents both positive and negative Nyquist frequency (+/-0.5), and is also purely real for real input.

For an odd number of input points, A[(n-1)/2] contains the largest positive frequency, while A[(n+1)/2] contains the largest negative frequency. The routine np.fft.fftfreq(A) returns an array giving the frequencies of corresponding elements in the output.

The routine np.fft.fftshift(A) shifts transforms and their frequencies to put the zero-frequency components in the middle, and np.fft.ifftshift(A) undoes that shift. (...)

When the input is purely real, its transform is Hermitian, i.e., the component at frequency is the complex conjugate of the component at frequency , which means that for real inputs there is no information in the negative frequency components that is not already available from the positive frequency components. The family of rfft functions is designed to operate on real inputs, and exploits this symmetry by computing only the positive frequency components, up to and including the Nyquist frequency. Thus, n input points produce n/2+1 complex output points. The inverses of this family assumes the same symmetry of its input, and for an output of n points uses n/2+1 input points.

Theoretical background: Nyquist–Shannon sampling theorem

Our function is band limited to [${- F}_{B},F_{B}{}_{}$] where B stands for bandwidth.

The Nyquist sampling theorem provides a prescription for the nominal sampling interval required to avoid aliasing and to perfectly reconstruct the signal:

The sampling frequency should be at least twice the highest frequency contained in the signal.

$F_{S} \geq {2F}_{B}$

$F_{N} = \ \frac{F_{S}}{2}$ is called Nyquist frequency .

Why should we use polar coordinates for frequencies?

For the $f_x$ and $f_y$ axes units are cycles/px and for ft axis is

Example: $f_x= 0,25 cpd$ (vertical grating)

$f_{\text{s\ }} =$\

fd is the magnitude of the spatial frequency in polar coordinates. The direction \theta is $$ \theta=\arctan (fx/fy) $$ (w.r.t cartesian coordinates) It can be seen as the rotation of the gaussian envelope.

$f_{x_{\text{norm}}\ } = \ f_{s}\text{cos\ }\theta$

$f_{y_{\text{norm}}\ } = \ f_{s}\text{sin\ }\theta$

${f_{\text{spatial}{}_{cyc/pix}} = {(f}_{\text{pre}f_{\text{cpd}}}*\ VA)/N_{x}}_{}$

$V_{X} = \ \frac{V_{\text{drift}}*T}{N_{f}}$

$ft\ cyc/frame = fspatial\ /Vx$

It’s not a deformation, projections (marginalization I would say) of gaussians are still gaussians

B_tf = B_sf*V

$V_{\max} = \ \frac{VA*T}{n_{f}}$

Vx=(Vdrift * T) / VAx

VAx or N_X for horizontal dispalcement

Vmax= Width /Length = (deg/sec)

Ex: vmax= 30deg/64frames*0,02sec/frame = 23deg/sec

Vx= 16 deg/sec

Vx_norm = 0,68

f_t = Vx_norm * f_s norm (case of drifting grating, no rotation in x-y axes)

f_t = 0,68 * 0,03152 = 0,02

Log-Gabor filter construction process in the Fourier domain (Hanssen and Hess, 2006, Journal of Vision)

an example physiology experiment using VSDI

In [1]:
%%writefile ../files/experiment_VSDI.py
#!/usr/bin/env python
"""
experiment_VSDI.py
Experiment designed for optical imaging showing a conversion of size from Fourier 
to screen coordinates and an export to a zipped file conaining BMP files (as the
video card has limited memory)

(c) Paula Sanz Leon - INT/CNRS


"""
import os
import scipy
import MotionClouds as mc
import numpy as np

# uncomment to preview movies
# vext, display = None, True

#------------------------- Zipped grating ------------------------------- #

name = 'zipped_grating'

#initialize - 
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
sf_0 = 0.2
B_sf = sf_0 / 4.
alpha = 1.0
V_X = 0.5
B_V = V_X / 10.

name_ = mc.figpath + name
for seed in range(424242, 424242+8):
    name_ = mc.figpath + name + '-seed-' + str(seed)
    mc.figures_MC(fx, fy, ft, name_, B_V=B_V, sf_0=sf_0, B_sf=B_sf, V_X=V_X, theta=np.pi/4., alpha=alpha, seed=seed, vext='.zip')

#-------------------- Narrowband vs Braodband experiment ---------------- #
vext = '.mpg'
#vext = '.mat'
#vext = '.zip'
#display = False

# Initialize frequency cube
N_X = 640.
N_Y = 480.
N_frame = 30. # a full period in time frames
fx, fy, ft = mc.get_grids(N_X, N_Y, N_frame)

# Experimental constants 
contrast = 0.5
seeds = 1
VA = 38.0958       # VSG constants for a viewing ditance 570 mm. 
framerate = 50.    # Refreshing rate in [Hz]
T = 0.6            # Stimulus duration [s] 
f_d = 0.5         # Desired spatial frequency [cpd]

# Clouds parameters      
B_V = 0.2     # BW temporal frequency (speed plane thickness)
B_sf = 0.15   # BW spatial frequency
theta = 0.0   # Central orientation
B_theta = np.pi/12 # BW orientation
verbose = False
alpha = 1.0

# Get normalised units
sf_0=0.1
V_X=0.5

def physicalUnits2discreteUnits(sf0_cpd, Bsf_cpd, v_dps, B_dps,pixelPitch,viewingDistance,frameRate):
    """%   % laptop monitor
    %   pixelPitch      = 0.22/10; % in cm
    %   viewingDistance = 50;
    %   sf0_cpd         = 4 ;
    %   Bsf_cpd         = 1 ;
    %   v_dps           = 4 ;
    %   B_dps           = 1 ;
    %   frameRate       = 20 ;
    % convert to machine units
    """
    cmPerDegree    = 2*viewingDistance*tand(1/2)
    pxPerDegree    = cmPerDegree/pixelPitch
    sf0            = sf0_cpd/pxPerDegree
    Bsf            = Bsf_cpd/pxPerDegree

    v              = v_dps/frameRate*pxPerDegree
    Bv             = B_dps/frameRate*pxPerDegree
    return sf0, Bsf, v, Bv


# take monitor parameters
# planar display
monitorRefresh = 60 ;
pixelPitch     = 0.2865/10 ;% pixelsize in cm

# experiments parameter
viewingDistance = 70 ;

frameRefresh   = monitorRefresh/3 ;
stimModRate    = 1 ;
numberOfReps   = 12 ;
framesPerMod   = frameRefresh/stimModRate/2 ; 


# Masks
# gaussian mask
sigma_mask_x = 0.15
sigma_mask_y = 0.2
x, y, t = mc.get_grids(N_X, N_Y, N_frame)
n_x, n_y = N_X, N_Y
gauss = np.exp(-(((x-172./n_x)**2/(2*sigma_mask_x**2)) + (((y-108./n_y)**2)/(2*sigma_mask_y**2))))


def tukey(n, r=0.5):
    '''The Tukey window, also known as the tapered cosine window, can be regarded as a 
    cosine lobe of width r * N / 2 that is convolved with a rectangle window of width 
    (1 - r / 2). At r = 1 it becomes rectangular, and at r = 0 it becomes a Hann window.
    http://www.mathworks.com/access/helpdesk/help/toolbox/signal/tukeywin.html
    '''
    # Special cases
    if r <= 0:
        return np.ones(n.shape) #rectangular window
    elif r >= 1:
        return np.hanning(n.shape)

    # Normal case
    x = np.linspace(0, 1, n)
    w = np.ones(x.shape)

    # first condition 0 <= x < r/2
    first_condition = x<r/2
    w[first_condition] = 0.5 * (1 + np.cos(2*np.pi/r * (x[first_condition] - r/2) ))

    # second condition already taken care of

    # third condition 1 - r / 2 <= x <= 1
    third_condition = x>=(1 - r/2)
    w[third_condition] = 0.5 * (1 + np.cos(2*np.pi/r * (x[third_condition] - 1 + r/2))) 

    return w

# Tukey mask - fading effect
tw_x = tukey(n=n_x, r=0.15)
tw_y = tukey(n=n_y, r=0.15)
w = np.tile(((np.outer(tw_y,tw_x))), (N_frame,1,1))
tukey_mask = w.T


# Get Random Clouds
name_ = mc.figpath + name

for seed in [123456 + step for step in range(seeds)]:
    name__ = mc.figpath + name + '-seed-' + str(seed) + '-sf0-' + str(sf_0).replace('.', '_') + '-V_X-' + str(V_X).replace('.', '_')
    # broadband 
    z = mc.envelope_gabor(fx, fy, ft, name_, B_sf=Bsf, sf_0=sf_0, theta=theta, B_V=B_V, B_theta = B_theta, alpha=alpha)
    movie = mc.figures(z, name=None, vext=vext, seed=seed, masking=True)    
    for label, mask in zip(['_mask', '_tukey_mask'], [gauss, tukey_mask]):
        name_ = name__ + '-cloud-' + label
        if anim_exist(name_): 
            movie = mc.rectif(movie*mask)
            mc.anim_save(movie, name_, display=False, vext=vext)

   # narrowband 
    z = mc.envelope(fx, fy, ft, name_, B_sf=B_sf/10., sf_0=sf_0, theta=theta, B_V=B_V, B_theta=B_theta, alpha=alpha)
    movie = mc.figures(z, name=None, vext=vext, seed=seed, masking=True)
    for label, mask in zip(['_mask', 'tukey_mask'], [gauss, tukey_mask]):
        name_ = name__ + '-blob-' + label
        if anim_exist(name_):
            movie = mc.rectif(movie*mask)
            mc.anim_save(movie, name_, display=False, vext=vext)
Overwriting ../files/experiment_VSDI.py

Testing Color

Testing the Color

In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
print(mc.envelope_color.__doc__)
    Returns the color envelope.

    In some circonstances, it is desirable to create a texture with a different "color"
    than that of natural images (that is where the envelope is in 1/f).

    Run 'test_color' notebook to see the effect of alpha
    alpha = 0 white
    alpha = 1 pink
    alpha = 2 red/brownian
    (see http://en.wikipedia.org/wiki/1/f_noise )
    
In [3]:
import os
name = 'color'
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
z = mc.envelope_color(fx, fy, ft, alpha=1.0, ft_0=1.)
mc.figures(z, name)
mc.in_show_video(name)
In [4]:
# explore parameters
for alpha in [0.001, 0.5, 1.0, 1.5, 2.0]:
    # resp. white(0), pink(1), red(2) or brownian noise (see http://en.wikipedia.org/wiki/1/f_noise
    name_ = name + '-alpha-' + str(alpha).replace('.', '_')
    z = mc.envelope_color(fx, fy, ft, alpha)
    mc.figures(z, name_)
    mc.in_show_video(name_)
In [5]:
for ft_0 in [0.125, 0.25, 0.5, 1., 2., 4., np.inf]:# time space scaling
    name_ = name + '-ft_0-' + str(ft_0).replace('.', '_')
    z = mc.envelope_color(fx, fy, ft, alpha=1., ft_0=ft_0)
    mc.figures(z, name_)
    mc.in_show_video(name_)
In [6]:
for size in range(5, 7):
    N_X, N_Y, N_frame = 2**size, 2**size, 2**size
    fx, fy, ft = mc.get_grids(N_X, N_Y, N_frame)
    ft_0 = N_X/float(N_frame)
    name_ = name + '-size-' + str(size).replace('.', '_')
    z = mc.envelope_color(fx, fy, ft, alpha=1., ft_0=ft_0)
    mc.figures(z, name_)
    mc.in_show_video(name_)

Testing Radial

Testing spatial frequency component

In [1]:
import numpy as np
np.set_printoptions(precision=5, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
#!rm -fr ../files/radial*
In [2]:
import MotionClouds as mc
import os
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
help(mc.envelope_radial)
Help on function envelope_radial in module MotionClouds:

envelope_radial(fx, fy, ft, sf_0=0.125, B_sf=0.1, ft_0=inf, loggabor=True)
    Returns the radial frequency envelope:
    
    Selects a preferred spatial frequency ``sf_0`` and a bandwidth ``B_sf``.
    
    Run the 'test_radial' notebook to see the explore the effect of ``sf_0`` and ``B_sf``, see
    http://motionclouds.invibe.net/posts/testing-radial.html

In [3]:
name = 'radial'
#initialize
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)

mc.figures_MC(fx, fy, ft, name)
verbose = False
mc.in_show_video(name)
In [4]:
N = 5
np.logspace(-5, 0., 7, base=2)
Out[4]:
array([ 0.03125,  0.05568,  0.09921,  0.17678,  0.31498,  0.56123,  1.     ])

exploring different frequency bandwidths

In [5]:
for B_sf in np.logspace(-4, 0., N, base=2):
    name_ = name + '-B_sf-' + str(B_sf).replace('.', '_')
    mc.figures_MC(fx, fy, ft, name_, B_sf=B_sf, verbose=verbose)
    mc.in_show_video(name_)

exploring different (median) central frequencies

In [6]:
for sf_0 in [0.01, 0.05, 0.1, 0.2, 0.4, 0.8]:
    name_ = name + '-sf_0-' + str(sf_0).replace('.', '_')
    mc.figures_MC(fx, fy, ft, name_, sf_0=sf_0, verbose=verbose)
    mc.in_show_video(name_)

Note that some information is lost in the last case as $f_0$ exceeds the Nyquist frequency.

exploring different speed bandwidths

In [7]:
for B_V in [0, 0.001, 0.01, 0.05, 0.1, 0.5, 1.0, 10.0]:
    name_ = name + '-B_V-' + str(B_V).replace('.', '_')
    mc.figures_MC(fx, fy, ft, name_, B_V=B_V, verbose=verbose)
    mc.in_show_video(name_)

Testing without using a log-normal distribution

In [8]:
for B_sf in np.logspace(-4, 0., N, base=2):
    name_ = name + '-B_sf_nologgabor-' + str(B_sf).replace('.', '_')
    mc.figures_MC(fx, fy, ft, name_, B_sf=B_sf, loggabor=False, verbose=verbose)
    mc.in_show_video(name_)    

checking that sf_0 is consistant with the number of cycles per period

By design, MotionClouds are created to be characterized by a given frequency range. Let's check on the raw image that this is correct.

In [9]:
downscale = 4
fx, fy, ft = mc.get_grids(mc.N_X/downscale, mc.N_Y/downscale, mc.N_frame/downscale)
N_X, N_Y, N_frame = fx.shape

def spatial_xcorr(im):
    import scipy.ndimage as nd
    N_X, N_Y, N_frame = im.shape
    Xcorr = np.zeros((N_X, N_Y))
    for t in range(N_frame):
        Xcorr += nd.correlate(im[:, :, t], im[:, :, t], mode='wrap')
    return Xcorr

mc_i = mc.envelope_gabor(fx, fy, ft, 
                         V_X=0., V_Y=0., 
                         sf_0=0.15, B_sf=0.03)
im = mc.random_cloud(mc_i)
Xcorr = spatial_xcorr(im)
plt.imshow(Xcorr)
Out[9]:
<matplotlib.image.AxesImage at 0x115ce2128>
In [10]:
xcorr = Xcorr.sum(axis=1)
xcorr /= xcorr.max()
plt.plot(np.arange(-N_X/2, N_X/2), xcorr)
print('First maximum reached at index ' + str(np.argmax(xcorr)) + ' - that is normal as we plot the autocorrelation')
print('Max = ', xcorr[np.argmax(xcorr)])
First maximum reached at index 32 - that is normal as we plot the autocorrelation
Max =  1.0

To find the period, intuitively, one could just find the zero of the gradient:

In [11]:
half_xcorr = xcorr[N_X//2:]
half_xcorr_gradient = np.gradient(xcorr)[N_X//2:]
plt.plot(half_xcorr)
plt.plot(half_xcorr_gradient, 'g')
idx = np.argmax(half_xcorr_gradient>0)
print('First zero gradient reached at index ' + str(idx) + ' - at the minimum of the xcorr')
print (half_xcorr_gradient[idx-1], half_xcorr_gradient[idx])
First zero gradient reached at index 4 - at the minimum of the xcorr
-0.210252703916 0.435683947279

We thus get in pixels half of the width of a period.

More precisely, let's fit with a sinusoid (that's yet another FFT...):

In [12]:
xcorr_MC = np.absolute(mc_i**2).sum(axis=-1).sum(axis=-1)
half_xcorr_MC = xcorr_MC[N_X//2:]
half_xcorr_MC /= half_xcorr_MC.max()
plt.plot(half_xcorr_MC, 'r--')

half_xcorr_FT = np.absolute(np.fft.fft(xcorr))[:N_X//2]
half_xcorr_FT /= half_xcorr_FT.max()
plt.plot(half_xcorr_FT)
idx = np.argmax(half_xcorr_FT)
print('Max auto correlation reached at frequency ' + str(idx))
print('Max auto correlation reached with a period ' + str(1.*N_X/idx))
Max auto correlation reached at frequency 9
Max auto correlation reached with a period 7.11111111111

Wrapping things up:

In [13]:
sf_0_ = [0.01, 0.05, 0.1, 0.2, 0.4, 0.8]
sf_0_ = np.linspace(0.01, .5, 12)
downscale = 4
fx, fy, ft = mc.get_grids(mc.N_X//downscale, mc.N_Y//downscale, mc.N_frame//downscale)

def period_px(mc_i):
    N_X, N_Y, N_frame = im.shape
    xcorr_MC = np.absolute(mc_i**2).sum(axis=-1).sum(axis=-1)
    half_xcorr_MC = xcorr_MC[N_X//2:]
    idx = np.argmax(half_xcorr_MC)
    return np.fft.fftfreq(N_X)[idx]

ppx = np.zeros_like(sf_0_)
for i, sf_0 in enumerate(sf_0_):
    mc_i = mc.envelope_gabor(fx, fy, ft, 
                             V_X=0., V_Y=0., 
                             sf_0=sf_0, B_sf=0.03)
    ppx[i] = period_px(mc_i)
    
plt.plot(sf_0_, ppx)
plt.xlabel('sf_0 [a. u.]')
plt.ylabel('observed freq [a. u.]')
Out[13]:
<matplotlib.text.Text at 0x11607a320>

Note that the range of frequencies is accessible using:

In [14]:
print(np.fft.fftfreq(fx.shape[0]))
[ 0.       0.01562  0.03125  0.04688  0.0625   0.07812  0.09375  0.10938
  0.125    0.14062  0.15625  0.17188  0.1875   0.20312  0.21875  0.23438
  0.25     0.26562  0.28125  0.29688  0.3125   0.32812  0.34375  0.35938
  0.375    0.39062  0.40625  0.42188  0.4375   0.45312  0.46875  0.48438
 -0.5     -0.48438 -0.46875 -0.45312 -0.4375  -0.42188 -0.40625 -0.39062
 -0.375   -0.35938 -0.34375 -0.32812 -0.3125  -0.29688 -0.28125 -0.26562
 -0.25    -0.23438 -0.21875 -0.20312 -0.1875  -0.17188 -0.15625 -0.14062
 -0.125   -0.10938 -0.09375 -0.07812 -0.0625  -0.04688 -0.03125 -0.01562]

Testing Grating

testing the orientation component

On this page, we concerntrate on the orientation variable and test the effect of changing the average orientation as well as the bandwidth.

In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
%load_ext autoreload
%autoreload 2
In [2]:
import MotionClouds as mc
import os
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
help(mc.envelope_orientation)
Help on function envelope_orientation in module MotionClouds:

envelope_orientation(fx, fy, ft, theta=0.0, B_theta=0.19634954084936207)
    Returns the orientation envelope:
    We use a von-Mises distribution on the orientation:
    - mean orientation is ``theta`` (in radians),
    - ``B_theta`` is the bandwidth (in radians). It is equal to the standard deviation of the Gaussian
    envelope which approximate the distribution for low bandwidths. The Half-Width at Half Height is
    given by approximately np.sqrt(2*B_theta_**2*np.log(2)).
    
    Run 'testing-grating.py' notebook to see the effect of changing theta and B_theta, see
    http://motionclouds.invibe.net/posts/testing-grating.html

In [3]:
name = 'grating'
#initialize
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)

z = mc.envelope_gabor(fx, fy, ft)
mc.figures(z, name, vext='.gif')
mc.figures(z, name)
mc.in_show_video(name)
In [4]:
# explore parameters
for sigma_div in [1, 2, 3, 5, 8, 13 ]:
    name_ = name + '-largeband-B_theta-pi-over-' + str(sigma_div).replace('.', '_')
    z = mc.envelope_gabor(fx, fy, ft, B_theta=np.pi/sigma_div)
    mc.figures(z, name_)
    mc.in_show_video(name_)
In [5]:
for div in [1, 2, 4, 3, 5, 8, 13, 20, 30]:
    name_ = name + '-theta-pi-over-' + str(div).replace('.', '_')
    z = mc.envelope_gabor(fx, fy, ft, theta=np.pi/div)
    mc.figures(z, name_)
    mc.in_show_video(name_)
In [6]:
V_X = 1.0
for sigma_div in [1, 2, 3, 5, 8, 13 ]:
    name_ = name + '-B_theta-pi-over-' + str(sigma_div).replace('.', '_') + '-V_X-' + str(V_X).replace('.', '_')
    z = mc.envelope_gabor(fx, fy, ft, V_X=V_X, B_theta=np.pi/sigma_div)
    mc.figures(z, name_)
    mc.in_show_video(name_)
In [7]:
for B_sf in [0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.8]:
    name_ = name + '-B_sf-' + str(B_sf).replace('.', '_')
    z = mc.envelope_gabor(fx, fy, ft, B_sf=B_sf)
    mc.figures(z, name_)
    mc.in_show_video(name_)

tuning the bandwidth

In the MotionClouds script, we define the orientation envelope based on a VonMises distribution:

In [8]:
def vonmises(th, theta, kappa, norm=True):
    if kappa==0:
        env = np.ones_like(th) 
    elif kappa==np.inf:
        env = np.zeros_like(th)
        env[np.argmin(th < theta)] = 1.
    else:
        env = np.exp(kappa*(np.cos(2*(th-theta))-1))
    if norm:
        return env/env.sum()
    else:
        return env

This definition handles also the extreme cases when the bandwidth is infinite (in which case we return a flat distribution) of null (so we get a dirac). Note also that the von Mises distribution is on a continuous variable (orientation), while we use a discrete fourier transform. In particular, some care should be taken for very narrow bandwidths.

In [9]:
N_theta = 16
theta, B_sf, B_V = np.pi/2, .1, .3
bins = 60
th = np.linspace(0, np.pi, bins, endpoint=False)
fig, ax = plt.subplots(1, 1, figsize=(13, 8))
B_theta_ = [np.pi/640, np.pi/32, np.pi/16, np.pi/8, np.pi/4, np.pi/2, np.inf]
for B_theta, color in zip(B_theta_, [plt.cm.hsv(h) for h in np.linspace(0, 1, len(B_theta_))]):
    kappa = 1./4/B_theta**2
    ax.plot(th, vonmises(th, theta, kappa), alpha=.6, color=color, lw=3)
    ax.fill_between(th, 0, vonmises(th, theta, B_theta), alpha=.1, color=color)
    ax.set_xlabel('orientation (radians)')
_ = ax.set_xlim([0, np.pi])

Getting everything peaking at one:

In [10]:
N_theta = 16
theta, B_sf, B_V = np.pi/2, .1, .3
bins = 360
th = np.linspace(0, np.pi, bins, endpoint=False)
fig, ax = plt.subplots(1, 1, figsize=(13, 8))
for B_theta, color in zip(B_theta_, [plt.cm.hsv(h) for h in np.linspace(0, 1, len(B_theta_))]):
    kappa = 1./4/B_theta**2
    ax.plot(th, vonmises(th, theta, kappa, norm=False), alpha=.6, color=color, lw=3)
    ax.fill_between(th, 0, vonmises(th, theta, kappa, norm=False), alpha=.1, color=color)
    ax.set_xlabel('orientation (radians)')
_ = ax.set_xlim([0, np.pi])

Let's check that the HWHH (half-width at half-height) is given by the bandwidth parameter B_theta. Note that for high values of B_theta, the minimum of the distribution may be higher than 1/2 and that this decriptive value may be ill-defined (see for instance Swindale, N. V. (1998). Orientation tuning curves: empirical description and estimation of parameters. Biological Cybernetics, 78(1):45-56.). We will use the broader definition by using the fully normalized von mises that peaks at one and with a minimum at 0 : $$M(\theta) = \frac {\exp( \kappa \cdot cos(2\theta)) - \exp( -\kappa) }{\exp(\kappa) - \exp( -\kappa) }$$ The different expression for the HWHH are for the Gaussian approximation: $$\theta_{HWHH} = \sqrt{2\ln(2)} \cdot \kappa $$ for the non-normalized von-Mises pdf (as in Swindale): $$\theta_{HWHH} = \frac 1 2 \cdot \arccos (1 + \frac{\kappa - \ln(2)}{\kappa})$$for the normalized (exact) solution: $$\theta_{HWHH} = \frac 1 2 \cdot \arccos (1 + \frac{\ln(\frac{1+e^{-2\kappa})}{2}}{\kappa})$$

In [11]:
N_theta = 160
theta, B_sf, B_V = 0, .1, .3
B_theta_ = np.linspace(0, np.pi/2, N_theta, endpoint=False)
bins = 3600
th = np.linspace(0, np.pi, bins, endpoint=False)
HWHH = np.zeros(N_theta)
for _, B_theta in enumerate(B_theta_):
    kappa = 1./4/B_theta**2
    dist = vonmises(th, theta, kappa)
    dist = (dist - dist.min())/(dist.max() - dist.min())
    ind_HWHH = np.argmax(dist<.5)
    HWHH[_] = th[ind_HWHH]
fig, ax = plt.subplots(1, 1, figsize=(13, 8))
ax.plot(B_theta_, HWHH, 'k+', label='observed')
ax.plot(B_theta_, B_theta_*np.sqrt(2*np.log(2)), 'r', label='gaussian')
k= 1./4/B_theta_**2
ax.plot(B_theta_, .5*np.arccos((-np.log(2)+k)/k), 'b', label='swindale')
ax.plot(B_theta_, .5*np.arccos(1+ np.log((1+np.exp(-2*k))/2)/k), 'g', label='exact')
ax.set_xlabel('B_theta (radians)')
ax.set_ylabel('HWHH (radians)')
ax.legend()
_ = ax.set_xlim([0, np.pi/2])
/usr/local/lib/python3.6/site-packages/ipykernel/__main__.py:8: RuntimeWarning: divide by zero encountered in double_scalars
/usr/local/lib/python3.6/site-packages/ipykernel/__main__.py:16: RuntimeWarning: divide by zero encountered in true_divide
/usr/local/lib/python3.6/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in true_divide
/usr/local/lib/python3.6/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in arccos

This shows that in the range of relevant values of B_theta (that is inferior to approximately 60 degrees) we can control the bandwidth as for a Gaussian or a von-Mises curve, while for higher values, the distribution can be considered to be flat.

Testing Speed

Testing the Speed

A moving pattern (stimulus/image) can be thought as a set of pixels and that each pixel intensities are translated from one frame to the next frame (the intensities, not the pixel:

I(x,t)=I(x+u,t+1); space is a vector with coordinates (x1,x2)’ and time is t. u=(u1,u2) is the 2D velocity. Motion Clouds represent 2D + t images. x(t), each pixel follows a 2D path as a function of time. So we have a space-time function f(x,t) by translating a 2D signal f0(x) with velocity u, i.e, f(x,t)=f0(x-ut).

The Fourier transform of such a function is

F(wx, wy, wt)=F0(wx,wy) \delta(u1wx + u2wy + wt)

F0 is the 2D transform of f0 and \delta is the dirac delta. The energy spectrum is nonzero on a plane whose orientation gives the velocity.

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 = 'speed'
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
help(mc.envelope_speed)
Help on function envelope_speed in module MotionClouds:

envelope_speed(fx, fy, ft, V_X=1.0, V_Y=0.0, B_V=0.5)
    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 V_X in [-1.0, -0.5, 0.0, 0.1, 0.5, 1.0, 4.0]:
    name_ = name + '-V_X-' + str(V_X).replace('.', '_')
    z = mc.envelope_gabor(fx, fy, ft, V_X=V_X)
    mc.figures(z, name_)
    mc.in_show_video(name_)
In [4]:
for V_Y in [-1.0, -0.5, 0.5, 1.0, 2.0]:
    name_ = name + '-V_Y-' + str(V_Y).replace('.', '_')
    z = mc.envelope_gabor(fx, fy, ft, V_X=0.0, V_Y=V_Y)
    mc.figures(z, name_)
    mc.in_show_video(name_)
In [5]:
for B_V in [0, 0.001, 0.01, 0.05, 0.1, 0.5, 1.0, 10.0]:
    name_ = name + '-B_V-' + str(B_V).replace('.', '_')
    z = mc.envelope_gabor(fx, fy, ft, B_V=B_V)
    mc.figures(z, name_)
    mc.in_show_video(name_)

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
%autoreload 2
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>