Attachment 'GoldenPyramid.py'

Download

   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 
  72     def add(self,pyr_1, pyr_2):
  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
 340         ranked_adresses = numpy.argsort(-numpy.abs(flat_pyramid))
 341         # see Eq. 10: estimate the z-value from the rank (first is 1, last is 0)
 342         flat_pyramid_q[ranked_adresses] = numpy.linspace(1.,0.,len(flat_pyramid))
 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
 439     image = pylab.imread('database/gris512.png')[:,:,0]
 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 
 468         pyramid.show(pyramid.add(im_pyr,pyramid.scalar(pyr_q,-1.)))
 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

You are not allowed to attach a file to this page.

welcome: please sign in