Statistics of time-varying images

In this notebook, we are interested in replicating the results from Dong and Attick (1995).

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
figsize=(21, 13) # fib
opts = {'fontsize':24}

In [2]:
import numpy as np


Takes a MPEG movie and produces a numpy array:

TODO:¶

We here provide a large data set of eye movement data on high-resolution natural movies, Hollywood trailers, and static images. Our analysis of this data and details on data collection can be found in

Michael Dorr, Thomas Martinetz, Karl Gegenfurtner, and Erhardt Barth. Variability of eye movements when viewing dynamic natural scenes. Journal of Vision, 10(10):1-17, 2010.

http://webmail.inb.uni-luebeck.de/inb-toolsdemos/FILES/movies-mpg.zip

In [3]:
try:
except:
# Let's take whatever movie that is supported by ffmpeg
MOVIE = '/Volumes/archives/11-07-14_RTC/montypython.mpg'
fs = 25. # framerate of the movie (in Hz)

N_frame = 512 # 8192 # 4096 # 16348 # 44197
start_frame = 2000

start, duration = start_frame / fs, N_frame / fs

import os
if not(os.path.isfile('/var/tmp/foo-001.png')):
os.system('ffmpeg -i ' + MOVIE + ' -ss ' + str(start) + ' -t ' + str(duration) +' -f image2 /var/tmp/foo-%03d.png')

# convert them to numpy
from numpy import zeros, save, mean, flipud
#from scipy import
from scipy.ndimage import zoom

border_x, border_y, actual_size, expected_size = 20, 80, 190, 128.
zfactor = expected_size/actual_size
xmin, xmax, ymin, ymax = border_x, border_x+actual_size , border_y, border_y+actual_size
n_x, n_y = zoom(mean(imread('/var/tmp/foo-001.png'), axis=2)[xmin:xmax, ymin:ymax], zfactor).shape

movie = zeros((n_y, n_x, N_frame))
for i_frame in range(N_frame):
name = '/var/tmp/foo-%03d.png' % (i_frame +1)
movie[:, :, i_frame] = zoom(image, zfactor)

#movie -= movie.mean()
##    print mov.shape
os.system('rm -f /var/tmp/foo-*.png')

save('/var/tmp/movie.npy', movie)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-3-02422c48ca21> in <module>()
1 try:
3 except:

/usr/local/lib/python3.5/site-packages/numpy/lib/npyio.py in load(file, mmap_mode, allow_pickle, fix_imports, encoding)
361     if isinstance(file, basestring):
--> 362         fid = open(file, "rb")
363         own_fid = True

FileNotFoundError: [Errno 2] No such file or directory: '/var/tmp/movie.npy'

During handling of the above exception, another exception occurred:

FileNotFoundError                         Traceback (most recent call last)
<ipython-input-3-02422c48ca21> in <module>()
25     zfactor = expected_size/actual_size
26     xmin, xmax, ymin, ymax = border_x, border_x+actual_size , border_y, border_y+actual_size
---> 27     n_x, n_y = zoom(mean(imread('/var/tmp/foo-001.png'), axis=2)[xmin:xmax, ymin:ymax], zfactor).shape
28
29     movie = zeros((n_y, n_x, N_frame))

2291
2292

1321             return handler(fd)
1322         else:
-> 1323             with open(fname, 'rb') as fd:
1324                 return handler(fd)
1325     else:

FileNotFoundError: [Errno 2] No such file or directory: '/var/tmp/foo-001.png'
In [ ]:
N_X, N_Y, N_frame = movie.shape
print 'The movie consists of ', N_frame, ' frames @ size (', N_X, ' x ', N_Y, ')'

In [ ]:
movie.mean(), movie.min(), movie.max()


Sample images:

In [ ]:
fig, axs = plt.subplots(5, 8, figsize=figsize) # fib
np.random.seed(42)
for axs_ in axs:
for ax in axs_:
_ = ax.imshow(movie[:, :, int(N_frame*np.random.rand())].T, cmap=plt.gray(), origin='lower')
ax.set_xticks([])
ax.set_yticks([])


We wish to compute the motion energy of this movie, first globally.

In [ ]:
import MotionClouds as mc

In [ ]:
import os, tempfile
from base64 import b64encode

def show_spectrum(F_in):
F_ = F_in.copy()
from IPython.core.display import display, Image, HTML
name = os.path.join(tempfile.mkdtemp(), 'spectrum')
mc.visualize(F_, name=name)           # Visualize the Fourier Spectrum
with open(name + mc.ext, "r") as image_file:
s = """<center><table border=none width=100%% height=100%%>
<tr><td width=33%%><center><img src="%s" width=100%%/></td>        </tr>
display(HTML(s))


Let's try first on a first chunk of 64 frames:

In [ ]:
chunk = 128
fx, fy, ft = mc.get_grids(N_X, N_Y, chunk)
fr = np.sqrt(fx**2 + fy**2)
#print fr.min(), fr.max(), F.shape[-1]

In [ ]:
F = np.fft.fftshift(np.absolute(np.fft.fftn(movie[:, :, :chunk]))**2)
print 'F is maximum at :', np.unravel_index(np.argmax(F), dims=F.shape), ' (we expected (', chunk/2, ', ', chunk/2, ', ', chunk/2, ')'
show_spectrum(np.sqrt(F))

In [ ]:
N_f = chunk/2 - 1 # making an histogram with N_f bins
f_bins = np.linspace(0., 0.5, N_f+1)

In [ ]:
from lmfit.models import PowerLawModel
mod = PowerLawModel()
pars = mod.guess((F/fr**2)[(chunk/2+1):, chunk/2, chunk/2], x=(f_bins[:-1]+f_bins[1:])/2)
out  = mod.fit((F/fr**2)[(chunk/2+1):, chunk/2, chunk/2], pars, x=(f_bins[:-1]+f_bins[1:])/2)
print(out.fit_report(min_correl=0.25))

In [ ]:
show_spectrum(np.sqrt(F)*fr)#*np.abs(ft)**.5)


Averaging over the movie:

In [ ]:
F = np.zeros_like(F)
for i_frame in np.arange(0, N_frame, chunk):
F += np.fft.fftshift(np.absolute(np.fft.fftn(movie[:, :, i_frame:(i_frame+chunk)]))**2)
show_spectrum(np.sqrt(F))

In [ ]:
show_spectrum(np.sqrt(F)*fr)

In [ ]:
from lmfit.models import PowerLawModel
mod = PowerLawModel()
pars = mod.guess(F[(chunk/2+1):, chunk/2, chunk/2], x=(f_bins[:-1]+f_bins[1:])/2)
out  = mod.fit(F[(chunk/2+1):, chunk/2, chunk/2], pars, x=(f_bins[:-1]+f_bins[1:])/2)
print(out.fit_report(min_correl=0.25))

In [ ]:
print 'F is maximum at :', np.unravel_index(np.argmax(F), dims=F.shape), ' (we expected ', [f/2 for f in F.shape], ')'


Averaging circularly

In [ ]:
N_f = chunk/2 - 1 # making an histogram with N_f bins
f_bins = np.linspace(0., 0.5, N_f+1)
F_rot = np.zeros((N_f, F.shape[-1]))
for i_bin in range(N_f):
f_slice = (f_bins[i_bin] < fr) *  ( fr < f_bins[i_bin+1])
F_rot[i_bin, :] = (f_slice * F).sum(axis=(0, 1))
#F_rot[i_bin, :] /= (f_bins[i_bin] + f_bins[i_bin+1])/2 # normalize by the integration area (theory)
F_rot[i_bin, :] /= f_slice.sum() # normalize by the integration area (numeric)
F_rot /= F_rot.max()

In [ ]:
fig, axs = plt.subplots(1, 1, figsize=figsize) # fib
axs.contourf(f_bins[:-1], np.linspace(-.5, .5, chunk), np.log(F_rot).T, cmap=plt.hot())
_ = axs.set_ylabel('temporal frequency', **opts)


Let's reproduce the results from Dong and Attick (1995, 2001). Basically, this spectrum is not separable (see page 342 of http://www.naturalimagestatistics.net/nis_preprintFeb2009.pdf):

In [ ]:
fig, axs = plt.subplots(1, 2, figsize=figsize) # fib
ft_bins = np.linspace(-.5, .5, chunk, endpoint=False)
axs[0].set_color_cycle(np.array([[1., 0., 0.]]) * grad[:, np.newaxis])
axs[1].set_color_cycle(np.array([[0., 1., 0.]]) * f_bins[:, np.newaxis])
for ft_bin in range(chunk): axs[0].loglog(f_bins[:-1], F_rot[:, ft_bin], alpha=.7)
for f_bin in range(N_f): axs[1].loglog(ft_bins, F_rot[f_bin, :], alpha=.7)
axs[0].loglog(f_bins[:-1], f_bins[:-1]**-1*2.e-8, 'k', lw=2, alpha=.7)
axs[0].loglog(f_bins[:-1], f_bins[:-1]**-2*5.e-5, 'k', lw=2, alpha=.7)
axs[1].loglog(ft_bins[chunk/2:-1], ft_bins[chunk/2:-1]**-1*2.e-8, 'k', lw=2, alpha=.7)
axs[1].loglog(ft_bins[chunk/2:-1], ft_bins[chunk/2:-1]**-2*5.e-5, 'k', lw=2, alpha=.7)
axs[0].set_xlabel('spatial frequency', **opts)
axs[1].set_xlabel('temporal frequency', **opts)
axs[0].axis('tight')
axs[1].axis('tight')
_ = axs[0].set_ylabel('Power', **opts)

In [ ]:
fig, axs = plt.subplots(1, 2, figsize=figsize) # fib
fig.suptitle('Fitting a power-law', **opts)
from lmfit.models import PowerLawModel
mod = PowerLawModel()

f_exp, ft_exp = np.zeros(chunk), np.zeros(N_f)
for i_ft in range(chunk):
pars = mod.guess(F_rot[:, i_ft], x=(f_bins[:-1]+f_bins[1:])/2)
out  = mod.fit(F_rot[:, i_ft], pars, x=(f_bins[:-1]+f_bins[1:])/2)
if i_ft == chunk/2: print(out.fit_report(min_correl=0.25))
f_exp[i_ft] = -out.params.get('exponent').value
for i_f in range(N_f):
v_dyn = -(np.arange(chunk)==chunk/2)
pars = mod.guess(F_rot[i_f, v_dyn], x=np.abs(ft_bins[v_dyn]))
out  = mod.fit(F_rot[i_f, v_dyn], pars, x=np.abs(ft_bins[v_dyn]))
if i_f == 0: print(out.fit_report(min_correl=0.25))
ft_exp[i_f] = -out.params.get('exponent').value

axs[0].plot(f_bins[:-1], ft_exp)
axs[0].set_xlabel('spatial frequency', **opts)
axs[0].axis('tight', **opts)
axs[1].plot(ft_bins, f_exp)
axs[1].set_xlabel('temporal frequency', **opts)
axs[1].axis('tight')

_ = axs[0].set_ylabel('exponent of the best-fit of a power-law', **opts)


but will be if one looks at different slices corresponding to different speeds.

Let's first define a slice in the speed plane (as we handle real signals, it's fine to consider just one side of the plane):

In [ ]:
N_f, N_v, N_dir = 15, 24, 24 # making an histogram with N_v - 1 speeds
v_max = 3.
speed = np.linspace(0., v_max, N_v+1)
direction = np.linspace(0, 2*np.pi, N_dir+1)
f_bins = np.linspace(0., 0.5, N_f+1)
f_bins = np.logspace(-.5, 0, N_f+1, base=10)*0.5

In [ ]:
v_bin, dir_bin = -2, 11
#v_slice = (speed[v_bin]*np.cos(direction[dir_bin])*fx + speed[v_bin]*np.sin(direction[dir_bin])*fy + ft > 0. )
#v_slice *= (speed[v_bin+1]*np.cos(direction[dir_bin])*fx + speed[v_bin+1]*np.sin(direction[dir_bin])*fy + ft < 0. ) # speed plane
B_V = .8
#fr[fr<.001]=1.e6
v_slice = np.exp(-(speed[v_bin]*np.cos(direction[dir_bin])*fx + speed[v_bin]*np.sin(direction[dir_bin])*fy + ft )**2/2/(B_V*fr)**2)
show_spectrum(v_slice*1.)
#print v_slice.min(), v_slice.max()

In [ ]:
F_v = np.zeros((N_f, N_v, N_dir))
for f_bin in range(N_f):
for dir_bin in range(N_dir):
for v_bin in range(N_v):
#v_slice = (speed[v_bin]*np.cos(direction[dir_bin])*fx + speed[v_bin]*np.sin(direction[dir_bin])*fy + ft > 0. )
#v_slice *= (speed[v_bin+1]*np.cos(direction[dir_bin])*fx + speed[v_bin+1]*np.sin(direction[dir_bin])*fy + ft < 0. ) # speed plane
v_slice = np.exp(-(speed[v_bin]*np.cos(direction[dir_bin])*fx + speed[v_bin]*np.sin(direction[dir_bin])*fy + ft )**2/2/(B_V*fr)**2)
f_slice = (f_bins[f_bin] < fr) *  ( fr < f_bins[f_bin+1])
F_v[f_bin, v_bin, dir_bin] = (v_slice * f_slice * F).sum()
#print f_bin, v_bin, dir_bin, F_v[f_bin, v_bin, dir_bin], (v_slice * f_slice).sum(), ( f_slice).sum(), (v_slice ).sum()
F_v[f_bin, v_bin, dir_bin] /= (v_slice * f_slice).sum() # (f_bins[i_bin] + f_bins[i_bin+1])/2 # normalize by the integration area

In [ ]:
fig, axs = plt.subplots(1, 1, figsize=figsize) # fib
axs.set_color_cycle(np.array([[0., 1., 0.]]) * grad[:, np.newaxis])
for v_bin in range(N_v-1):
axs.loglog(f_bins[:-1], F_v[:, v_bin], 'o')
axs.loglog(f_bins[:-1], np.mean(F_v[:, v_bin], axis=-1), '-' )
axs.axis('tight')
_ = axs.set_ylabel('Amplitude')


If we multiply this spectrum by $f_r^\alpha$, with $\alpha$ being an appropriate scaling, we get:

In [ ]:
N_test = 100
error, alphas = np.zeros(N_test), np.linspace(1.5, 4.5, N_test)
for i, alpha in enumerate(alphas):
spectrum = ((f_bins[1:]+f_bins[:-1])[:, np.newaxis, np.newaxis]/2)**(alpha+1) * F_v[:, :, :]
spectrum /= spectrum.mean()
error[i] = np.std(spectrum, axis = (0,2)).mean()
fig, axs = plt.subplots(1, 1, figsize=figsize) # fib
axs.plot(alphas, error)
axs.text(2.5, 0.2, r'Error is minimal at $\alpha$=' + '%0.2f' % alphas[error.argmin()], **opts)
axs.set_xlabel('Alpha', **opts)
axs.axis('tight')
axs.set_ylim([0, error.max()])
_ = axs.set_ylabel('Error', **opts)

In [ ]:
alpha = alphas[error.argmin()] # 2.3 for Dong - 2.7 for Hyvarinen
fig, axs = plt.subplots(1, 1, figsize=figsize) # fib
axs.set_color_cycle(np.array([[0., 1., 0.]]) * grad[:, np.newaxis])
for f_bin in range(1, N_f):
axs.semilogy(speed[:-1], ((f_bins[f_bin]+f_bins[f_bin+1])/2)**(alpha+1) * F_v[f_bin, :, :], 'o')
axs.semilogy(speed[:-1], np.mean(((f_bins[f_bin]+f_bins[f_bin+1])/2)**(alpha+1) * F_v[f_bin, :, :], axis = -1), '-')
#    axs.loglog(speed[:-1], F_v[f_bin, :, :])
axs.set_xlabel('velocity')
#axs.axis('tight')
_ = axs.set_ylabel('Amplitude')

In [ ]: