Attachment 'image.py'

Download

   1 #/usr/bin/python
   2 # -*- coding: utf8 -*-
   3 """
   4 image.py
   5 
   6 
   7 
   8 
   9 A set of image processing sub-routines
  10 
  11 See http://incm.cnrs-mrs.fr/LaurentPerrinet/Publications/Perrinet08spie
  12 
  13 $Id: image.py 2526 2009-01-22 08:32:16Z lup $
  14 """
  15 import pylab, numpy, scipy, progressbar
  16 import scipy.ndimage as nd
  17 from scipy.fftpack import fft2, fftshift, ifft2, ifftshift
  18 from scipy import real,absolute
  19 #### init / database handling
  20 def load_in_database(url_database):
  21     """
  22     Loads a random image from directory url_database
  23 
  24     """
  25     from numpy import ceil, random
  26     from pylab import imread
  27     import os
  28     filelist = os.listdir(url_database)
  29     filelist.remove('.svn')
  30     i = random.randint(0,len(filelist))
  31     #name=['../database/images/', str(i), '.png']
  32     name = url_database + '/' + filelist[i]
  33     #print name
  34     # a = numpy.frombuffer(Image.open('myimage.jpg'))
  35     # Image.frombuffer(b).save('filtered.png')
  36     image = imread(name)[:,:,0]
  37 
  38     #image = normalize(image)
  39     return image
  40 
  41 def patch(image, size, threshold = 0.2):
  42     """
  43     takes a subimage of size s (a tuple)
  44 
  45     does not accept if energy is relatively below a threshold (flat image)
  46 
  47     """
  48     from numpy import std, ceil, random
  49 
  50     image_size_h,image_size_v = image.shape
  51     if size[0] > image_size_h or size[1] > image_size_v:
  52         raise('Patch too big')
  53 
  54 
  55     energy = std(image[:])
  56     energy_ = 0;
  57 
  58     while energy_ < threshold*energy:
  59         #if energy_ > 0: print 'dropped patch'
  60         x_rand = ceil((image_size_h-size[0])*random.rand())
  61         y_rand = ceil((image_size_v-size[1])*random.rand())
  62         image_ = image[(x_rand):(x_rand+size[0]),(y_rand):(y_rand+size[1])]
  63         energy_ = std(image_[:])
  64 
  65 
  66     return image_
  67 
  68 def filter_grid(shape = (16, 16)):
  69     """
  70 
  71     N_X, N_Y = 32, 32 # size of the filter
  72 
  73     returns a meshgrid
  74     /!\ values between -1 and 1
  75 
  76     """
  77     from scipy import mgrid
  78 
  79     return mgrid[-1:1:1j*shape[0],-1:1:1j*shape[1]]
  80 
  81 
  82 def filter_list(shape = (16, 16), n_orientation = 6 ):
  83     """
  84 
  85     N_X, N_Y = 32, 32 # size of the filter
  86 
  87     returns a list of gabor filters
  88 
  89     """
  90     from scipy import mgrid
  91     list = []
  92     N_X,N_Y = shape
  93     [fX,fY] =  filter_grid(N_X,N_Y)
  94 
  95 
  96     for i_orientation in range(n_orientation):
  97         orientation = 2*numpy.pi * i_orientation / n_orientation
  98         filter = gabor(fX,fY, orientation)
  99         list.append(filter)
 100 
 101     return list
 102 
 103 #### low-level core processing
 104 
 105 def coco(image, filter, normalize = True):
 106     """
 107     Returns the correlation coefficient
 108 
 109     """
 110     from scipy.signal import correlate2d
 111     from numpy import sqrt
 112     coco = correlate2d(image, filter, mode = 'same')
 113     if normalize:
 114         coco /= sqrt(energy(image)*energy(filter))
 115 
 116     return coco
 117 
 118 def energy(image):
 119     return  numpy.mean(image.ravel()**2)
 120 
 121 def normalize(image):
 122     image -= numpy.mean(image.ravel())
 123     image /= energy(image)**.5
 124 
 125     return image
 126 
 127 
 128 def FTfilter(image,FTfilter):
 129     from scipy.fftpack import fft2, fftshift, ifft2, ifftshift
 130     from scipy import real
 131     FTimage = fftshift(fft2(image)) * FTfilter
 132     return real(ifft2(ifftshift(FTimage)))
 133     #return real(ifft2(fft2(image)* FTfilter))
 134 
 135 
 136 #### whitening
 137 def olshausen_whitening_filt(size, f_0 = .78, alpha = 4., N = 0.01):
 138     """
 139     Returns the whitening filter used by (Olshausen, 98)
 140 
 141     f_0 = 200 / 512
 142 
 143     /!\ you will have some problems at dewhitening without a low-pass
 144 
 145     """
 146     from scipy import mgrid, absolute
 147     fx, fy = mgrid[-1:1:1j*size[0],-1:1:1j*size[1]]
 148     rho = numpy.sqrt(fx**2+fy**2)
 149     K_ols = (N**2 + rho**2)**.5 * low_pass(size, f_0 = f_0, alpha = alpha)
 150     K_ols /= numpy.max(K_ols)
 151 
 152     return  K_ols
 153 
 154 def low_pass(size, f_0, alpha):
 155     """
 156     Returns the low_pass filter used by (Olshausen, 98)
 157 
 158     parameters from Atick (p.240)
 159     f_0 = 22 c/deg in primates: the full image is approx 45 deg
 160     alpha makes the aspect change (1=diamond on the vert and hor, 2 = anisotropic)
 161 
 162     """
 163 
 164     from scipy import mgrid, absolute
 165     fx, fy = mgrid[-1:1:1j*size[0],-1:1:1j*size[1]]
 166     rho = numpy.sqrt(fx**2+fy**2)
 167     low_pass = numpy.exp(-(rho/f_0)**alpha)
 168 
 169     return  low_pass
 170 
 171 
 172 def whitening_filt(size = (512,512),
 173                          url_database = 'database/Yelmo',
 174                          n_learning = 400,
 175                          N = 0.1,
 176                          f_0 = .8, alpha = 1.4,
 177                          N_0 = .5,
 178                          recompute = False):
 179     """
 180     Returns the correlation filter in FT space.
 181 
 182     Computes the average power spectrum = FT of cross-correlation, the mean decorrelation
 183     is given for instance by (Attick, 92).
 184 
 185     """
 186     import shelve
 187 
 188     results = shelve.open('white'+ str(size[0]) + '-' + str(size[1]) + '.mat')
 189     try:
 190         if recompute:
 191             print('Recomputing the whitening filter')
 192             boing
 193         K = results['K']
 194 
 195     except:
 196         print ' Learning the whitening filter'
 197         pbar=progressbar.ProgressBar(widgets=["calculating", " ", progressbar.Percentage(), ' ',
 198             progressbar.Bar(), ' ', progressbar.ETA()], maxval=n_learning)
 199         power_spectrum = 0. # power spectrum
 200         for i_learning in range(n_learning):
 201             image = patch(load_in_database(url_database = url_database), size)
 202             image = normalize(image) #TODO : is this fine?
 203             power_spectrum += numpy.abs(fft2(image))**2
 204             pbar.update(i_learning)
 205 
 206         power_spectrum = fftshift(power_spectrum)
 207 
 208 
 209         power_spectrum /= numpy.mean(power_spectrum)
 210         # formula from Attick:
 211         M = numpy.sqrt( power_spectrum / (N**2 + power_spectrum)) * low_pass(size, f_0 = f_0, alpha = alpha)
 212         K = M / numpy.sqrt( M**2 * (N**2 + power_spectrum) + N_0**2)
 213         K /= numpy.max(K) # normalize energy :  DC is one <=> xcorr(0) = 1
 214 
 215         pbar.finish()
 216 
 217         results['K'] = K
 218         results['power_spectrum'] = power_spectrum
 219     results.close()
 220 
 221     return K
 222 
 223 def whitening(image):
 224     """
 225     Returns the whitened image
 226     """
 227     K = whitening_filt(size = image.shape, recompute = False)
 228 
 229     if not(K.shape == image.shape):
 230         K = whitening_filt(size = image.shape, recompute = True)
 231 
 232     return FTfilter(image,K)
 233 
 234 def dewhitening(white):
 235     """
 236     Returns the dewhitened image
 237 
 238     """
 239     K = whitening_filt(white.shape)
 240 
 241     return FTfilter(white,1/K)
 242 
 243 
 244 def retina(image):
 245     """
 246     A dummy retina processsing
 247 
 248     """
 249     white = whitening(image)
 250     white = normalize(white) # mean = 0, std = 1
 251     return white
 252 
 253 #### filter definition
 254 def dog(shape, orientation , sigma = .3, K = 3.):
 255     """
 256     returns a dog shape
 257 
 258     /!\ X, Y values between -1 and 1
 259 
 260     """
 261     [X,Y] =  filter_grid(shape)
 262     R2 = X**2 + Y**2
 263     dog = numpy.exp( - R2 / 2 / sigma**2)- 1/K**2*numpy.exp( - R2 / 2 / sigma**2/K**2)
 264     print sum(dog.ravel())/energy(dog)**.5
 265     dog /= sum(dog.ravel()**2)**.5
 266 
 267     return dog
 268 
 269 
 270 def trans(shape, (x_pos,y_pos)):
 271     """
 272     returns a translation filter
 273 
 274     /!\ X, Y values between -1 and 1
 275 
 276 
 277 
 278     """
 279     raise("not implemented yet!")
 280 
 281     [fX,fY] =  filter_grid(shape)
 282     R2 = X**2 + Y**2
 283     trans = numpy.exp( - R2 / 2 / sigma**2)- 1/K**2*numpy.exp( - R2 / 2 / sigma**2/K**2)
 284     print sum(dog.ravel())/energy(dog)**.5
 285     trans /= sum(trans.ravel()**2)**.5
 286 
 287     return trans
 288 
 289 def gabor(width, orientation , sigma = .3, f_s = 1.6, phase = 0.):
 290     """
 291     returns a gabor shape
 292     /!\ X, Y values between -1 and 1
 293 
 294     """
 295     [X,Y] =  filter_grid(shape=(width,width))
 296     R2 = X**2 + Y**2
 297     S = numpy.cos(orientation)* X + numpy.sin(orientation)* Y
 298     gabor = numpy.exp( - R2 / 2 / sigma**2) * numpy.sin(2* numpy.pi * f_s * S + phase)
 299 
 300     return gabor
 301 
 302 def ring(shape, f_0 = .5, df = numpy.sqrt(2./3.)):
 303     """
 304     ring filter : corresponds to a log-dog
 305 
 306     """
 307     eps = 10e-8
 308     [fX,fY] =  filter_grid(shape)
 309     logf = numpy.log2(numpy.sqrt(fX**2 + fY**2)+eps) # logf<0 bacause f<1
 310     ring = numpy.exp( - (logf-numpy.log(f_0+eps))**2 / 2 / df**2)
 311     ring /= numpy.max(ring.ravel())
 312 
 313     return ring
 314 
 315     
 316 #### multi-scale / pyramid
 317 
 318 def resize(image,n_x_d,n_y_d):
 319     #from scipy.interpolate import interpolate
 320     assert image.ndim == 2
 321 
 322     n_x, n_y = image.shape
 323     
 324 ##    index_x = [int(idx) for idx in numpy.linspace(0,n_x-1,n_x_d)]
 325 ##    image_r = image[index_x,:]
 326 ##    index_y = [int(idx) for idx in numpy.linspace(0,n_y-1,n_y_d)]
 327 ##    image_r = image_r[:,index_y]
 328 #    else:
 329     image_r = nd.zoom(image,(n_x_d/float(n_x),n_y_d/float(n_y)),  order=1, mode='mirror')
 330 ##    if n_x < n_x_d:
 331 ##        ratio = float(n_x)/float(n_x_d)
 332 ##        image_r = nd.gaussian_filter(image_r,sigma = 1./ratio )
 333         
 334 ##    if n_x > n_x_d:
 335 ##        ratio = float(n_x)/float(n_x_d)
 336 ##        image_r = nd.gaussian_filter(image_r, sigma = ratio )
 337 
 338 
 339     return image_r
 340 
 341 
 342 def fibo(n):
 343     """ Fibonacci series starting at rank 0 and 1 with (1,1)
 344 
 345     >> fibo(1)
 346     1
 347 
 348     >> fibo(12)
 349     144
 350 
 351     """
 352     if n < 0:
 353         raise('n>0')
 354     elif n<=2:
 355         return 1
 356     else:
 357         return fibo(n-1) + fibo(n-2)
 358 
 359 
 360 def n_fibo(f):
 361     """
 362     Returns the minimum n which Fibonacci number is superior to f
 363 
 364     >> n_fibo(143)
 365     11
 366 
 367     >> n_fibo(144)
 368     12
 369 
 370     """
 371     n = 0
 372     while fibo(n)<=f:
 373         n+=1
 374 
 375     return n-1
 376 
 377 def rot_trans(square,stab):
 378     rot_trans_square =  square.copy()
 379     # rotate by -pi/2 around 0,0
 380     # xmin is ymin
 381     rot_trans_square[0] = square[2]
 382     # xmax is ymax
 383     rot_trans_square[1] = square[3]
 384     # ymin is -xmax
 385     rot_trans_square[2] = - square[1]
 386     # ymax is -xmin
 387     rot_trans_square[3] = - square[0]
 388 
 389     # and now translate by stab
 390     rot_trans_square[0] += stab # = [pos += fibonacci(n_level) for pos in _square]
 391     rot_trans_square[1] += stab
 392     rot_trans_square[2] += stab
 393     rot_trans_square[3] += stab
 394     return rot_trans_square
 395 
 396 
 397 def golden_grid(n_level, discrete = True):
 398     """ define the grid of a golden pyramid
 399 
 400     n_level  is the number of scales
 401 
 402     outputs a list of [xmin, xmax, ymin, ymax] squares of the pyramid
 403 
 404     """
 405     PHI= (1 + numpy.sqrt(5))/2.
 406 
 407     if discrete:
 408         stab = fibonacci(n_level)
 409     else:
 410         stab = PHI**n_level
 411 
 412     grid = [numpy.array([0,stab,0,stab])]
 413 
 414     if n_level>0:
 415         for _square in golden_grid(n_level-1):
 416             grid.append(rot_trans(_square,stab))
 417 
 418     return grid
 419 
 420 def square(edge):
 421     """ transforms a square [xmin, xmax, ymin, ymax] into a path (a list) starting from
 422     the bottom left
 423 
 424     """
 425     #print 'edge', edge
 426     return [edge[0],edge[1],edge[1],edge[0],edge[0]],[edge[2],edge[2],edge[3],edge[3],edge[2]]
 427 
 428 ###########################################################################
 429 if __name__ == '__main__':
 430 
 431 
 432 
 433     ##url_database = '/data/archive/DyVA/NatImages/tuebingen/Png/'
 434     ##url_database = 'database/icabench' # olshausen98'#
 435     ##
 436     ##image = patch(load_in_database(url_database),  (128,128))
 437 
 438     #image = pylab.imread('database/gris512.png')[:,:,0]
 439     image = pylab.imread('database/gris128.png')[:,:,0]
 440     #image = normalize(pylab.imread('database/square_arnold.png')[:,:,0])
 441         #pylab.show()
 442     ##list = filter_list()
 443     ##
 444     ##for filter in list:
 445     ##    I_coco = coco(image, filter)
 446     ##    pylab.subplot(121)
 447     ##    pylab.imshow(I_coco)
 448     ##    pylab.subplot(122)
 449     ##    pylab.imshow(filter)
 450     ##    #pylab.show()
 451     #n_x, n_y = image.shape
 452     #factor = (numpy.sqrt(5) +1. )/2.
 453     #n_x_d, n_y_d = numpy.ceil(n_x/factor), numpy.ceil(n_y/factor)
 454     #
 455     #pylab.subplot(121)
 456     #image_r = congrid(image,(n_x_d, n_y_d),method ='spline')
 457     #pylab.imshow(image_r)
 458     #image_u = congrid(image_r,(n_x, n_y),method ='spline')
 459     #pylab.subplot(122)
 460     #pylab.imshow(image_u)
 461     #pylab.show()
 462     [X,Y] = filter_grid(shape = (512,512))
 463     image = .5 + .5*numpy.cos(2* 12 * numpy.arctan2(X,Y))
 464     pylab.imshow(image)
 465     from scipy.misc import imsave
 466     imsave('database/test.png', image)
 467 
 468 # make movies: http://www.psychopy.org/reference/psychopy.makeMovies-pysrc.html

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