## Attachment 'GoldenPyramid.py'

```   1 #/usr/bin/python
2 ## -*- coding: utf8 -*-
3 """
4 GoldenPyramid.py
5
6 See http://incm.cnrs-mrs.fr/LaurentPerrinet/Publications/Perrinet08spie
7
8 \$Id: GoldenPyramid.py 2526 2009-01-22 08:32:16Z lup \$
9 """
10 import pylab, numpy, progressbar
11 from scipy.fftpack import fft2, fftshift, ifft2, ifftshift
12 from  image import load_in_database, patch, retina, normalize, resize, whitening, dewhitening
13
14
15 class GoldenPyramid(object):
16     """
17     defines a GoldenPyramid
18
19     """
20     def __init__(self, (n_x,n_y), n_quant = 512):
21         """
22         initializes the GoldenPyramid structure
23
24         shape is an array with the size of each level in the pyramid
25         n_quant is the quantization of coefficients for the modulation function
26         """
27         self.phi = (numpy.sqrt(5) +1. )/2. # golden number
28         #self.n_level = n_fibo(size) # alternate way to get perfect pyramids
29         self.n_levels = int(numpy.max(numpy.ceil(numpy.log(n_x)/numpy.log(self.phi)),
30                             numpy.ceil(numpy.log(n_y)/numpy.log(self.phi))))
31
32         self.n_x = n_x
33         self.n_y = n_y
34
35         self.shape =[]
36         self.total_pixels = 0
37         #print self.shape
38         for i_level in range(self.n_levels):
39             n_x_d, n_y_d = int(numpy.ceil(n_x/(self.phi**i_level))), int(numpy.ceil(n_y/(self.phi**i_level)))
40             self.total_pixels += n_x_d * n_y_d
41             self.shape.append((n_x_d, n_y_d))
42
43         self.cumlevel_size = numpy.cumsum([numpy.prod(shape) for shape in self.shape])
44         self.n_quant = n_quant
45
46
47     ## LOW LEVEL OPERATIONS ON PYRAMIDS
48     def zeros(self, column = False):
49         """
50         Returns a zero pyramid
51
52         """
53         zero_pyr = []
54         for i_level in range(self.n_levels):
55             if column:
56                 zero_pyr.append(numpy.zeros((self.shape[i_level][0],self.shape[i_level][1],self.n_column)))
57             else:
58                 zero_pyr.append(numpy.zeros(self.shape[i_level]))
59         return zero_pyr
60
61
62     def max(self,im_pyr):
63         flat_pyr = self.flat(im_pyr, quantize = False)
64         return numpy.max(abs(flat_pyr))
65
66     def scalar(self,im_pyr, scalar):
67         pyr = []
68         for i_level in range(self.n_levels):
69             pyr.append( im_pyr[i_level] * scalar )
70         return pyr
71
73         pyr = []
74         for i_level in range(self.n_levels):
75             pyr.append(pyr_1[i_level] + pyr_2[i_level])
76         return pyr
77
78
79     def threshold(self, pyr, rank):
80         """
81         Thresholds a pyramid by zeroing out all coefficients inferior in norm to
82         the threshold computed with the mean Mod.
83
84         0<= rank <= 1 : relative rank
85
86         """
87
88         pyr_max = self.max(pyr)
89
90         im_pyr = []
91         for i_level in range(self.n_levels):
92             im = pyr[i_level]
93             Mod = - numpy.sort(-numpy.abs(im.ravel()))
94             index = numpy.int(numpy.prod(im.shape)*rank)
95             value =  Mod[index]
96             im_pyr.append( im * ( numpy.abs(im) > value ) )
97
98         return im_pyr
99
100
101     def show(self,pyramid, border = 1, text = None):
102         height = 8.  #inches
103         fig = pylab.figure(num=1,figsize=(height*self.phi,height))
104
105         xmin, ymin, size = 0, 0, 1.
106 #        for i_level in range(min(self.n_levels,6)):
107         for i_level in range( min(8,self.n_levels) ):
108             #print xmin, ymin, size
109             ax = fig.add_axes([xmin/self.phi,ymin, size/self.phi, size])
110
111             ax.axis('off')
112             if i_level ==0 and text is not(None):
113                 ax.text(0.01, .9, text,
114                 horizontalalignment='left',
115                 verticalalignment='bottom',
116                 transform = ax.transAxes,
117                 bbox=dict(facecolor='green', alpha=0.5),
118                 fontsize=36
119                 )
120
121
122             im = pyramid[i_level].copy()
123
124             pmax = numpy.max(numpy.abs(im).ravel())
125             if pmax == 0: pmax = 1
126
127 #            print numpy.min(col_im.ravel()), numpy.max(col_im.ravel())
128             if len(im.shape) == 3:
129                 im[0:border,:,:] = -pmax
130                 im[:,0:border,:] = -pmax
131                 im[-1,:,:] = -pmax
132                 im[:,-1,:] = -pmax
133                 rgba_im = numpy.zeros((im.shape[0],im.shape[1],4))
134                 rgba_im[:,:,0:3] = numpy.dot(numpy.abs(im),self.column_color)
135                 rgba_im[:,:,3] = numpy.sum(numpy.abs(im),axis = 2)
136
137                 ax.imshow(rgba_im  / pmax, interpolation = 'nearest')#,  vmin = -pmax, vmax = pmax)#, cmap=pylab.gray())
138             else:
139                 im[0:border,:] = -pmax
140                 im[:,0:border] = -pmax
141                 im[-1,:] = -pmax
142                 im[:,-1] = -pmax
143                 ax.imshow(im,  vmin = -pmax, vmax = pmax, interpolation = 'nearest', cmap=pylab.gray())
144
145             i_orientation = numpy.mod(i_level,4)
146             #print i_orientation
147             if i_orientation==0:
148                 xmin += size
149                 ymin += size /self.phi**2
150             elif i_orientation==1:
151                 xmin += size /self.phi**2
152                 ymin +=  -size /self.phi
153             elif i_orientation==2:
154                 xmin += -size/self.phi
155             elif i_orientation==3:
156                 ymin += size
157             size /= self.phi
158         pylab.draw()
159         return fig
160
161     def show_column(self, border = 2e-3):
162         height = 8.  #inches
163         fig = pylab.figure(num=3,figsize=(height,height))
164         N = numpy.floor(numpy.sqrt(self.n_column))
165         size = 1/N - 2*border
166         size_x, size_y = self.column[0].shape
167         for i_column in range(self.n_column):
168             ax = fig.add_axes([border/2 + i_column//N/N ,border/2 + i_column%N/N, size, size])
169             ax.axis('off')
170             ax.imshow(self.column[i_column], interpolation = 'nearest', cmap=pylab.gray())
171             pylab.text(.0, .85, str(i_column), fontsize=12, bbox=dict(facecolor='white', alpha=0.5))
172         pylab.draw()
173         return fig
174
175     ## QUANTIZATION
176
177     def get_P(self, im_pyr, normed = True):
178         """
179         Computes the cumulative probability for each level of the pyramid.
180
181         """
182         P =[]
183         if normed:
184             norm = self.max(im_pyr)
185             if norm >0: im_pyr = self.scalar(im_pyr, 1./norm)
186
187         ax = numpy.linspace(0.,1., self.n_quant)
188         #zax = numpy.linspace(0.,1., self.n_quant)
189         for i_level in range(self.n_levels):
190             vec = numpy.abs(im_pyr[i_level].ravel())
191             vec = numpy.sort(vec)
192             vec = numpy.concatenate(([0], vec, [1]))
193             vec = numpy.unique(vec)
194
195             P_ = numpy.interp(ax, vec , numpy.linspace(0.,1.,len(vec)))
196             P.append(P_)
197
198         return P
199
200     def invert(self, P):
201         """
202         Transforms the cumulative probability in the Modulation function.
203
204         """
205         Mod =[]
206         ax = numpy.linspace(0.,1., self.n_quant)
207         for i_level in range(self.n_levels):
208             Mod.append(numpy.interp(ax[::-1], P[i_level] + ax * 1e-12, ax))
209
210         return Mod
211
212     def get_Mod(self, im_pyr):
213         """
214         Computes the modulation function for each level of the pyramid.
215
216         """
217         P = self.get_P(im_pyr)
218         return self.invert(P)
219
220     def learn(self,url_database, n_learning, coding = 'laplacian', normed = True):
221         """
222         Returns the mean Mod over n_learning steps by choosing randomly image
223         patches in the database.
224
225         """
226         pbar=progressbar.ProgressBar(widgets=["calculating", " ", progressbar.Percentage(), ' ',
227             progressbar.Bar(), ' ', progressbar.ETA()], maxval=n_learning)
228         print 'Learning the modulation function'
229         for i_learning in range(n_learning):
230             image_ = retina(patch(load_in_database(url_database = url_database), (self.n_x,self.n_y)))
231
232             # gets the correct coding method
233             if coding == 'laplacian':
234                 pyramid_ = self.LaplacianPyramid(image_)
235             elif coding == 'ssc':
236                 Lpyramid_ = self.LaplacianPyramid(image_)
237                 col_pyr = self.ColumnPyramid(Lpyramid_)
238                 pyramid_ = self.mp(col_pyr, threshold = 1e-3, step_max = 1e8, verbose = False)
239             elif coding == 'column': # linear column coding
240                 pyramid_ = self.ColumnPyramid(self.LaplacianPyramid(image_))
241             else:
242                 raise('/!\ Coding method ', coding , ' is not known')
243
244             P_ = self.get_P(pyramid_, normed = normed)
245             # actually learns (by a shifting average)
246             eta= 1./(i_learning+1.)
247             if i_learning==0: # booting
248                 P = P_
249             else:
250                 P = self.add(self.scalar(P, 1. - eta),self.scalar(P_, eta))
251             pbar.update(i_learning)
252
253         self.P = P
254         # transform in a modulation function
255         pbar.finish()
256
257
258
259     def flat(self, im_pyr, quantize = True):#, normed = True):
260         """ Flattens the pyramid
261
262         Transforms the pyramid into a vector. Useful to use other numpy stuff.
263         if quantize, gives for every level the f_j(value) the signed z value in
264         [-1,1] (1 means high).
265
266         It is assumed the z value is normalized if quantized
267
268         """
269
270         flat_pyramid = numpy.array([])
271
272         if quantize :
273             rho_axis = numpy.linspace(0.,1.,self.n_quant)
274
275         for i_level in range(self.n_levels):
276             # the  signal
277             vec = im_pyr[i_level].ravel()
278
279             if quantize :
280                 # keeps the sign
281                 vec_sign = numpy.sign(vec)
282                 # the  absolute signal
283                 vec = numpy.abs(vec)
284                 vec = numpy.interp(vec, rho_axis, self.P[i_level])
285                 # restores the sign
286                 vec = vec * vec_sign
287
288             flat_pyramid=numpy.concatenate((flat_pyramid,vec))
289
290         return flat_pyramid
291
292     def deflat(self, flat_pyramid, dequantize = True):
293         """ Reconstructs the shape of the pyramid from the flat vector.
294
295         if dequantize is True, retrieves the value from the z-value (= 1-rank),
296         it then assumes that the pyramid is normalized.
297
298         """
299         q_pyramid = []
300         counter_left = 0
301         n_col = len(flat_pyramid) / self.total_pixels
302
303         rho_axis = numpy.linspace(0.,1.,self.n_quant)
304         for i_level in range(self.n_levels):
305             # takes a sub-set corresponding to one level
306             width = numpy.prod(self.shape[i_level])*n_col
307             # retrieves signed Z values between -1 and 1
308             vec = flat_pyramid[counter_left:(counter_left+width)]
309
310             if dequantize :
311                 # retrieves Z values between 0 and 1
312                 vec_sign = numpy.sign(vec)
313                 vec = numpy.abs(vec)
314                 vec = numpy.interp(vec, self.P[i_level], rho_axis)
315                 # this was fully coded
316                 vec = vec  * vec_sign
317
318             if n_col > 1:
319                 im = vec.reshape(self.shape[i_level][0],self.shape[i_level][1],n_col)
320             else:
321                 im = vec.reshape(self.shape[i_level])
322
323             q_pyramid.append(im)
324             counter_left += width
325
326         return q_pyramid
327
328     def quantize(self,im_pyr):#, normed = True):
329         """
330         Returns the pyramid with the rank-quantized values.
331
332         """
333
334         norm = self.max(im_pyr)
335         if norm >0: im_pyr = self.scalar(im_pyr, 1./norm)
336         flat_pyramid = self.flat(im_pyr, quantize = True)
337         flat_pyramid_sign = numpy.sign(flat_pyramid)
338         flat_pyramid_q = numpy.zeros(flat_pyramid.shape)
339         # from the z-values across all levels compute the real rank quantization
341         # see Eq. 10: estimate the z-value from the rank (first is 1, last is 0)
343         # in the coding channel just the addresses (and polarities) are passed
344         flat_pyramid_q = flat_pyramid_q * flat_pyramid_sign
345         # get back the scalar values from the estimated z values
346         q_pyramid = self.deflat(flat_pyramid_q, dequantize = True)
347         q_pyramid = self.scalar(q_pyramid,norm)
348
349         return q_pyramid
350
351
352     ## LINEAR LAPLACIAN PYRAMID
353     def LaplacianPyramid(self,image):
354         """
355         Returns a pyramid of details from the base until the lower level as a
356         list of images.
357
358         """
359         pyramid =[]
360
361         for i_level in range(self.n_levels-1):
362             n_x, n_y = self.shape[i_level]
363             n_x_, n_y_ = image.shape
364             if not(n_x == n_x_ and  n_y == n_y_):
365                 raise("Incompatible shape for image at level ", i_level)
366             n_x_d, n_y_d = self.shape[i_level+1]
367             image_downscaled = resize(image,n_x_d, n_y_d)
368             image_upscaled = resize(image_downscaled,n_x, n_y)
369             image_residual = image - image_upscaled
370             pyramid.append(image_residual)
371             image = image_downscaled
372
373         pyramid.append(image_downscaled) # last level
374
375         return pyramid
376
377
378     def LaplacianDecode(self,pyramid):
379         """
380         Returns the image reconstructed from the pyramid.
381
382         """
383         image = pyramid[self.n_levels-1]
384         for i_level in range(self.n_levels-1,0,-1):
385             n_x, n_y = self.shape[i_level-1]
386             image = resize(image,n_x, n_y)
387             image +=  pyramid[i_level-1]
388
389         return image
390
391     ## COLUMNLET TRANSFORM
392     def make_column(self, width = 13, # should be uneven to define a center
393                     n_orientation = 16, type = 'logGabor'):
394         """
395         Creates the column
396
397         """
398
399         self.column = []
400         self.n_column = n_orientation
401         self.column_color = numpy.zeros((self.n_column,3)) # convert to RGB to show
402
403         for i_column, orientation in enumerate(numpy.linspace(0, numpy.pi, self.n_column, endpoint = False)):
404             # define filters
405             filter = gabor(width = width, orientation = orientation, phase = numpy.pi / 2.)
406             #filter *= filter>0
407             filter /= numpy.sqrt(numpy.sum(filter.ravel()**2))
408             self.column.append(filter)
409
410             # define transform for the color in figures (n_columns -> RGB)
411             self.column_color[i_column,0] = 1-3.*(1./2.-numpy.abs(orientation/numpy.pi -.5))
412             self.column_color[i_column,1] = 1-3.*(1./2.-numpy.abs((numpy.mod(orientation/numpy.pi-1./3.,1)) -.5))
413             self.column_color[i_column,2] = 1-3.*(1./2.-numpy.abs((numpy.mod(orientation/numpy.pi-2./3.,1)) -.5))
414
415
416         self.column_color = self.column_color * (self.column_color>0)
417         from scipy.signal import correlate as corr
418         self.lateral = []
419         for i_column in range(self.n_column):
420             R = numpy.zeros((2*width-1,2*width-1,self.n_column))
421             for j_column in range(self.n_column):
422                 R[:,:,j_column] = corr(self.column[i_column],self.column[j_column])
423             self.lateral.append(R)
424
425
426 #def _test():
427 #    import doctest
428 #    doctest.testmod()
429 #####################################
430 #
431 #_test()
432 if __name__ == '__main__':
433     #### Main
434     """
435     Some examples of use for the class
436
437     """
438     # whitening
440     image = normalize(image)
441     white = whitening(image)
442     if True:
443         im_rec = dewhitening(white)
444         print 'whitening error=', (image-im_rec).std()/ image.std()
445
446     # (linear) Golden Laplacian Pyramid
447     pyramid =  GoldenPyramid(image.shape)
448     #pyramid.multiscale()
449     #im_pyr = []
450     #for i_level in range(pyramid.n_levels):
451     #    im_pyr.append(pyramid.filter[i_level])
452     im_pyr = pyramid.LaplacianPyramid(white)
453     if True:
454         fig = pyramid.show(im_pyr)
455         #pylab.show()
456         im_rec = pyramid.LaplacianDecode(im_pyr)
457         print 'pyramid error=', (white-im_rec).std() / white.std()
458         pylab.matshow(dewhitening(im_rec))
459
460     if True:
461         # quantification of pyramid
462         pyramid.P = pyramid.get_P(im_pyr)
463         pyr_q = pyramid.quantize(im_pyr)
464         im_rec = pyramid.LaplacianDecode(pyr_q)
465         print 'quantification error (exact Mod)=', (white-im_rec).std() / white.std()
466         pylab.matshow(dewhitening(im_rec))
467
469
470         pyramid.learn('database/Yelmo', 20)
471         im_rec = pyramid.LaplacianDecode(pyramid.quantize(im_pyr))
472         print 'quantification error (learn)=', (white-im_rec).std() / white.std()
473         pylab.matshow(dewhitening(im_rec))
474         #pylab.scatter( pyramid.flat(im_pyr, quantize = False), pyramid.flat(pyr_q, quantize = False))
475
476     #for level in range(pyramid.n_levels):
477     #    numpy.sum(scipy.real(ifft2(ifftshift(pyramid.filter[level]))).ravel()**2)/ pyramid.max( pyramid.pyramid(pyramid.get_filter((x_pos, y_pos, level)) ))
```

## Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
• [get | view] (2014-10-25 14:02:40, 16.6 KB) [[attachment:GoldenPyramid.py]]
• [get | view] (2014-10-25 14:02:40, 3.3 KB) [[attachment:Makefile]]
• [get | view] (2014-10-25 14:02:42, 8.4 KB) [[attachment:SpikeCoding.tex]]
• [get | view] (2012-03-17 18:25:02, 80198.6 KB) [[attachment:database.zip]]
• [get | view] (2014-10-25 14:14:14, 22.4 KB) [[attachment:experiment_SpikeCoding.py]]
• [get | view] (2014-10-25 14:14:14, 6.9 KB) [[attachment:experiment_whitening.py]]
• [get | view] (2014-10-25 14:14:15, 11.7 KB) [[attachment:image.py]]
• [get | view] (2012-03-17 18:27:17, 918.3 KB) [[attachment:perrinet08spie.pdf]]
• [get | view] (2014-10-25 14:14:25, 54.8 KB) [[attachment:perrinet08spie.tex]]
• [get | view] (2012-03-17 18:27:07, 14240.3 KB) [[attachment:perrinet08spie_talk.pdf]]
All files | Selected Files: delete move to page