1 #/usr/bin/python
   2 # -*- coding: utf8 -*-
   3 """
   4 motion_plant.py
   5 
   6 Fonctions de base pour la detection de mouvement.
   7 
   8 Voir http://incm.cnrs-mrs.fr/LaurentPerrinet/ProjetBim
   9 
  10 """
  11 __author__ = "Laurent Perrinet INCM - CNRS"
  12 __revision__ = "$Id: MotionParticles.py 2818 2009-03-12 17:27:34Z lup $"
  13 
  14 import pylab
  15 from numpy import *
  16 import matplotlib.cm as cm
  17 
  18 def generate(N_X = 128., N_Y = 128., width = 2., noise = 0.02,
  19                  frequence = 6., V_X= .1):
  20     """
  21     Crée 2 images successives d'une sinusoide à la vitesse V_X
  22 
  23     """
  24 
  25     from scipy import mgrid
  26     r = width / 2.
  27     x, y = mgrid[r*(-1+.5/(N_X+1)):r*(1-.5/(N_X+1)):1j*N_X,
  28                  r*(-1+.5/(N_Y+1)):r*(1-.5/(N_X+1)):1j*N_Y]
  29 
  30     I1 = sin( 2 * pi * frequence *  y/width  )
  31     I1 += noise * random.randn(N_X, N_Y)
  32     I2 = sin( 2 * pi * (frequence *  y/width - V_X) ) 
  33     I2 += noise * random.randn(N_X, N_Y)
  34     return I1, I2
  35 
  36 
  37 def show_images(I1,I2):
  38     """
  39     créée une figure à 2 panneaux avec les 2 images
  40     
  41     TODO: rajouter des titres
  42     """
  43     fig = pylab.figure(figsize=(12,6))
  44     ax1 = pylab.subplot(121)
  45     ax1.imshow(I1, cmap = cm.gray)
  46     ax2 = pylab.subplot(122)
  47     ax2.imshow(I2, cmap = cm.gray)
  48     pylab.draw()
  49 
  50 
  51 def proba(V, image1, image2, contrast = 1.):
  52     """" 
  53     Calcule la proba sur les vitesses V pour les 2 images.
  54     
  55     Il faut définir V avec la fonction velocity_grid, par exemple:
  56     
  57     >>> V= velocity_grid(v_max = 25.)
  58     
  59     (v_max est important il donne la valeur max des proba à tester)
  60     
  61     il faut aussi définir les images. par exemple:
  62     
  63     >>> I1, I2 = generate(V_X= 1.)
  64 
  65     pour invoquer:
  66     
  67     >>> P= proba(V, I1, I2, contrast = 1.)
  68 
  69     """
  70     E = MotionEnergy(image1, image2, V, width = 2.0)
  71 
  72     return MotionProba(E, V, contrast = contrast)
  73 
  74 def show_global(V, proba):
  75     """
  76     Montre la densité de probabilité.
  77     
  78     TODO: rajouter le nom des axes
  79     
  80     """
  81     fig = pylab.figure(figsize=(6,6))
  82     a = fig.add_axes((0.15,0.1,.8,.8))
  83     pylab.scatter(V[0,:], V[1,:], c =proba) # X is going down, Y going to the right
  84     pylab.axis('tight')
  85 ##    pylab.draw()
  86     return fig
  87 
  88 
  89 def ArgMaxProba(V, proba):
  90     """
  91     Retourne la vitesse correpond au maximum de proba
  92     
  93     
  94     """
  95     
  96     address = proba.argmax()
  97 
  98     return V[:,address]
  99 
 100 
 101 
 102 #### fonctions de bas niveau
 103 
 104 def prior(V):
 105     """
 106     
 107     Returns a prior PDF of
 108     http://en.wikipedia.org/wiki/V._A._Epanechnikov
 109 
 110     """
 111     N_V = V.shape[1]
 112     prior = ones(N_V)
 113     v_max = sqrt(max( V[0,:]**2 + V[1,:]**2 ))
 114     
 115     for i_V in range(N_V):
 116         V_r = sqrt(V[0,i_V]**2 + V[1,i_V]**2)
 117         prior[i_V] = (v_max - V_r)**2
 118     return prior
 119 
 120 
 121 def MotionProba(energy, V, contrast = 1.):
 122     """" Returns the PDF of motion energy on RFs.
 123 
 124     :param contrast measures
 125 
 126     """
 127 
 128     N_V = energy.shape
 129     P = zeros(N_V)
 130 
 131     E = energy.copy()
 132     E -= amin(E)
 133     sigma_V = median(sqrt(E.ravel())) / contrast
 134     P = exp( - E / 2 / sigma_V**2)
 135     #P *= prior(V)
 136     P /= sum(P)
 137 
 138     return P
 139 
 140 def MotionEnergy(image1, image2, V, width = 2.0):
 141     """" Returns the motion energy on V.
 142 
 143     Input
 144     -----
 145     image : a 3D input sequence,
 146     V: a velocity grid, (2,N_V) array,
 147     RF: a list of center and widths, (4,N_RF) array.ArrayType
 148 
 149     Output
 150     ------
 151     E: the energy at velocity V, (N_V,) array
 152 
 153     """
 154     N_V = V.shape[1]
 155     energy = zeros(N_V)
 156 
 157     for i_V in range(N_V):
 158         I_residual = (image2 - translate(image1,V[:,i_V]))**2
 159         I_energy = image2**2
 160         I_energy_past = translate(image1,V[:,i_V])**2
 161         energy[i_V] = sum( I_residual )
 162         energy[i_V] /= sqrt(sum( I_energy ))
 163         energy[i_V] /= sqrt(sum( I_energy_past ))
 164 
 165     return energy
 166 
 167 
 168 def translate(image,vec):
 169     from scipy import mgrid
 170     f_x, f_y = mgrid[-1:1:1j*image.shape[0],-1:1:1j*image.shape[1]]
 171     trans = exp(-1j*pi*(vec[1]*f_x + vec[0]*f_y))
 172     return FTfilter(image, trans)
 173 
 174 def FTfilter(image,FTfilter):
 175     from scipy.fftpack import fftn, fftshift, ifftn, ifftshift
 176     from scipy import real
 177     FTimage = fftshift(fftn(image)) * FTfilter
 178     return real(ifftn(ifftshift(FTimage)))
 179 
 180 
 181 def velocity_grid(N_v = 12, N_theta=24, v_max= 25., hex = True):
 182     """
 183     Creates a set of velocity samples on a radial grid
 184 
 185     The velocity is in pixels/frame.
 186 
 187     """
 188     V = zeros((2,N_v*N_theta +1))
 189     N_V = N_v*N_theta
 190 
 191     v_rho = v_max - linspace(0 ,v_max, N_v, endpoint=False)
 192     if hex:
 193         parity = arange(N_v) %2
 194     else:
 195         parity = 0
 196     v_theta = linspace(0,2*pi, N_theta, endpoint=False)
 197     for i_theta, theta in enumerate(v_theta):
 198         V[0,i_theta*N_v:(i_theta+1)*N_v] = cos(theta + parity * pi / N_theta) * v_rho
 199         V[1,i_theta*N_v:(i_theta+1)*N_v] = sin(theta + parity * pi / N_theta) * v_rho
 200     V[0,-1] = 0.
 201     V[1,-1] = 0.
 202 
 203     return V