2014-10-21 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: 
    movie = np.load('/var/tmp/movie.npy')
except:
    # Let's take whatever movie that is supported by ffmpeg
    MOVIE = '/Volumes/archives/11-07-14_RTC/montypython.mpg'
    #MOVIE = '/Users/lolo/Downloads/movies_mpg/beach.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 pylab import imread
    #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)
        image = flipud(mean(imread(name)[xmin:xmax, ymin:ymax], axis=2)).T
        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:
----> 2     movie = np.load('/var/tmp/movie.npy')
      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))

/usr/local/lib/python3.5/site-packages/matplotlib/pyplot.py in imread(*args, **kwargs)
   2288 @docstring.copy_dedent(_imread)
   2289 def imread(*args, **kwargs):
-> 2290     return _imread(*args, **kwargs)
   2291 
   2292 

/usr/local/lib/python3.5/site-packages/matplotlib/image.py in imread(fname, format)
   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>
            </table></center>""" % ('data:image/png;base64,' + b64encode(image_file.read()))
        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_xlabel('radial frequency', **opts)
_ = 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)
grad =  np.abs(ft_bins)
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
grad =  np.abs(np.linspace(-1, 1., N_v))
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.set_xlabel('radial frequency')
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
grad =  np.abs(np.linspace(-1, 1., N_v))
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 [ ]:
 

Comments

Comments powered by Disqus