Statistics of time-varying images

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

import numpy as np

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)

N_X, N_Y, N_frame = movie.shape
print 'The movie consists of ', N_frame, ' frames @ size (', N_X, ' x ', N_Y, ')'

movie.mean(), movie.min(), movie.max()

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([])

import MotionClouds as mc

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)) 

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]

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))

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

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))

show_spectrum(np.sqrt(F)*fr)#*np.abs(ft)**.5)

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))

show_spectrum(np.sqrt(F)*fr)

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))

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

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()

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)

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)

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)

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

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()

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

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')

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)

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')