1
2
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)
84 pylab.axis('tight')
85
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
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
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