pyvision.point.GaborJets

1 # PyVision License 2 # 3 # Copyright (c) 2009 David S. Bolme 4 # All rights reserved. 5 # 6 # Redistribution and use in source and binary forms, with or without 7 # modification, are permitted provided that the following conditions 8 # are met: 9 # 10 # 1. Redistributions of source code must retain the above copyright 11 # notice, this list of conditions and the following disclaimer. 12 # 13 # 2. Redistributions in binary form must reproduce the above copyright 14 # notice, this list of conditions and the following disclaimer in the 15 # documentation and/or other materials provided with the distribution. 16 # 17 # 3. Neither name of copyright holders nor the names of its contributors 18 # may be used to endorse or promote products derived from this software 19 # without specific prior written permission. 20 # 21 # 22 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 23 # ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 24 # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 25 # A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR 26 # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 27 # EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 28 # PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 29 # PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 30 # LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 31 # NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 32 # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 33 34 ''' 35 This is implemented based on Wiskott et.al. "Elastic Bunch Graph Matching" PAMI 1997 36 ''' 37 38 39 import numpy as np 40 import unittest 41 import os.path 42 import pyvision as pv 43 from pyvision.analysis.face import EyesFile 44 #import time 45

46 -class GaborWavelet:

47 - def __init__(self, freq, oreint, sigma):

48 ''' 49 From Wiskott 96 50 ''' 51 self.freq = freq 52 self.oreint = oreint 53 self.sigma = sigma 54 55 self.k = np.array([freq*np.cos(oreint),freq*np.sin(oreint)],'f') 56 self.k2 = np.dot(self.k,self.k) 57 self.sigma2 = sigma*sigma

58

59 - def mask(self,size):

60 w,h = size 61 62 #m = np.zeros(size,np.complex64) 63 x = np.arange(-w/2,w/2).reshape(w,1)*np.ones(size) 64 x = np.concatenate((x[w/2:,:],x[:w/2,:]),axis=0) 65 y = np.arange(-h/2,h/2).reshape(1,h)*np.ones(size) 66 y = np.concatenate((y[:,h/2:],y[:,:h/2]),axis=1) 67 68 m2 = self(x,y) 69 return m2

70 71

72 - def __call__(self,x,y):

73 if isinstance(x,int) and isinstance(y,int): 74 #print "ints" 75 x = np.array([x,y],'d') 76 77 x2 = np.dot(x,x) 78 k = self.k 79 k2 = self.k2 80 sigma2 = self.sigma2 81 82 dot_kx = np.dot(k,x) 83 tmp = (k2/sigma2)*np.exp(-k2*x2/(2.0*sigma2))*(np.exp(1j*dot_kx) - np.exp(-sigma2/2.0)) 84 return tmp 85 else: 86 x = np.array([x,y],'d') 87 88 x2 = np.array(x*x).sum(axis=0) 89 k = self.k 90 k2 = self.k2 91 sigma2 = self.sigma2 92 93 dot_kx = np.array(k.reshape(2,1,1)*x).sum(axis=0) 94 tmp = (k2/sigma2)*np.exp(-k2*x2/(2.0*sigma2))*(np.exp(1j*dot_kx) - np.exp(-sigma2/2.0)) 95 return tmp

96 97

98 -class FilterBank:

99 ''' 100 This class uses the FFT to quickly compute corelations and convolutions 101 for an image. The algorithm precomputes the filters in frequency space 102 and uses only one FFT to to transform the image into frequency space. 103 '''

104 - def __init__(self,tile_size=(128,128),window=None,preprocess=None):

105 self.tile_size = tile_size 106 self.filters = [] 107 self.window = None 108 self.preprocess = preprocess 109 if window != None: 110 self.window = window(tile_size)

111 112

113 - def addFilter(self,f,recenter=False):

114 ''' 115 f can be a function f(x,y) is defined over x = (-w/2, w/2] and 116 y = (-h/2,h/2] and should be centered on the coord 0,0. 117 118 TODO: At some point this function should be expanded to take filters 119 represented by arrays. 120 ''' 121 if recenter == True: 122 raise NotImplementedError 123 if isinstance(f,GaborWavelet): 124 filt = np.fft.fft2(f.mask(self.tile_size)) 125 self.filters.append(filt) 126 else: 127 w,h = self.tile_size 128 m = np.zeros((w,h),np.complex64) 129 for x in range(-w/2,w/2): 130 for y in range(-h/2,h/2): 131 m[x,y] = f(x,y) 132 filt = np.fft.fft2(m) 133 self.filters.append(filt.conj())

134 135

136 - def convolve(self,im,ilog=None):

137 if isinstance(im,pv.Image): 138 im = im.asMatrix2D() 139 140 w,h = self.tile_size 141 assert im.shape[0] == w 142 assert im.shape[1] == h 143 144 if self.preprocess != None: 145 im = self.preprocess(im) 146 147 if self.window != None: 148 im = self.window*im 149 150 if ilog != None: 151 ilog.log(pv.Image(im)) 152 153 c = len(self.filters) 154 155 result = np.zeros((w,h,c),np.complex64) 156 157 fft_image = np.fft.fft2(im) 158 159 for i in range(c): 160 product = self.filters[i]*fft_image 161 result[:,:,i] = np.fft.ifft2(product) 162 163 return result

164 165

166 -def createGaborKernels():

167 '''Create gabor wavelets from Wiskott 1999''' 168 kernels = [] 169 sigma = 2*np.pi 170 for freq in [np.pi*2.**(-(i+2.)/2.) for i in range(5)]: 171 for orient in [i*np.pi/8. for i in range(8)]: 172 #print "Freq: %8.5f Orient: %8.5f Sigma: %8.5f"%(freq,orient,sigma) 173 wavelet = GaborWavelet(freq,orient,sigma) 174 kernels.append(wavelet) 175 176 return kernels

177 178

179 -class GaborFilters:

180

181 - def __init__(self, kernels=createGaborKernels(), tile_size=(128,128),window=None,preprocess=None):

182 self.kernels = kernels 183 self.bank = FilterBank(tile_size=tile_size,window=window,preprocess=preprocess) 184 self.k = np.zeros((len(kernels),2),dtype=np.float64) 185 for i in range(len(kernels)): 186 wavelet = kernels[i] 187 self.bank.addFilter(wavelet) 188 self.k[i,0] = wavelet.k[0] 189 self.k[i,1] = wavelet.k[1]

190 #print "K",self.k 191 192

193 - def convolve(self,im, ilog=None):

194 data = self.bank.convolve(im,ilog=ilog) 195 return GaborImage(data,self.kernels,self.k)

196 197

198 -class GaborImage:

199 - def __init__(self,data,kernels,k):

200 self.data = data 201 self.kernels = kernels 202 self.k = k

203

204 - def extractJet(self,pt,subpixel=True):

205 x = int(np.round(pt.X())) 206 y = int(np.round(pt.Y())) 207 x = max(min(x,self.data.shape[0]-1),0) 208 y = max(min(y,self.data.shape[0]-1),0) 209 jet = self.data[x,y,:] 210 211 return GaborJet(jet,self.kernels,self.k,x,y,pt,subpixel)

212

213 - def locatePoint(self,jet,start_pt=None,method="Simple"):

214 ''' 215 If start_pt == None perform a grid search with a spacing of one half 216 the longest Gabor wavelength. Otherwize start at start_pt and follow 217 the Jacobian to the local minimum. 218 219 @param jet: the an example jet. 220 @param start_pt The point to start the search from. 221 ''' 222 if start_pt == None: 223 # Compute the gate to use in the search 224 kx = self.k[:,0] 225 ky = self.k[:,1] 226 kt = np.sqrt(kx*kx + ky*ky).min() 227 gate = int(round(0.5*np.pi/kt)) 228 229 # search for best similarity 230 231 best_sim = -1.0 232 best_pt = pv.Point(0,0) 233 w,h,_ = self.data.shape[:] 234 for y in range(gate,h,gate): 235 for x in range(gate,w,gate): 236 pt = pv.Point(x,y) 237 novel = self.extractJet(pt,subpixel=False) 238 sim = novel.simDisplace(jet) 239 if sim > best_sim: 240 best_sim = sim 241 best_pt = pt 242 243 start_pt = best_pt 244 245 pt = start_pt 246 247 #print pt 248 249 novel = self.extractJet(pt,subpixel=False) 250 d = novel.displace(jet,method=method) 251 pt = pv.Point(novel.x + d[0], novel.y+d[1]) 252 #print pt 253 254 novel = self.extractJet(pt,subpixel=False) 255 d = novel.displace(jet,method=method) 256 pt = pv.Point(novel.x + d[0], novel.y+d[1]) 257 #print pt 258 259 novel = self.extractJet(pt,subpixel=False) 260 d = novel.displace(jet,method=method) 261 pt = pv.Point(novel.x + d[0], novel.y+d[1]) 262 #print pt 263 264 novel = self.extractJet(pt,subpixel=False) 265 sim = novel.simPhase(jet) 266 #pt = pv.Point(novel.x + d[0], novel.y+d[1]) 267 #print pt 268 269 return pt,sim,novel

270

271 - def show(self,*args,**kwargs):

272 print self.data.shape 273 tiles = [] 274 w,h,n = self.data.shape 275 for i in range(n): 276 mat = self.data[:,:,i] 277 tiles.append(pv.Image(mat.real)) 278 tiles.append(pv.Image(mat.imag)) 279 mont = pv.ImageMontage(tiles,layout=(8,10),tileSize=(w,h)) 280 mont.show(*args,**kwargs)

281 282

283 -class GaborJet:

284 - def __init__(self,jet,kernels,k,x,y,pt,subpixel):

285 self.jet = jet 286 self.kernels = kernels 287 self.k = k 288 self.x = x 289 self.y = y 290 291 re = np.real(jet) 292 im = np.imag(jet) 293 294 self.mag = np.sqrt(re*re + im*im) 295 self.phase = np.arctan2(re,im) 296 297 if subpixel: 298 d = np.array([[pt.X()-x],[pt.Y()-y]]) 299 comp = np.dot(self.k,d) 300 self.phase -= comp.flatten() 301 self.jet = self.mag*np.exp(1.0j*self.phase)

302 303

304 - def displace(self,gj,method="Blocked",**kwargs):

305 ''' 306 @param method: can be one of "Blocked", "Iter", "Simple" 307 ''' 308 if method=="Blocked": 309 return self.displaceBlocked(gj,**kwargs) 310 elif method=="Iter": 311 return self.displaceIter(gj,**kwargs) 312 elif method=="Simple": 313 return self.displaceSimple(gj,**kwargs) 314 else: 315 raise ValueError("Unknown displacement estimation method: %s"%method)

316

317 - def displaceSimple(self,gj):

318 m1 = self.mag 319 m2 = gj.mag 320 p1 = self.phase 321 p2 = gj.phase 322 kx = self.k[:,0] 323 ky = self.k[:,1] 324 325 p = p1 - p2 326 327 mask = p > np.pi 328 p -= 2.0*np.pi * mask 329 330 mask = p < -np.pi 331 p += 2.0*np.pi * mask 332 333 334 phi_x = (m1*m2*kx*p).sum() 335 phi_y = (m1*m2*ky*p).sum() 336 gam_xx = (m1*m2*kx*kx).sum() 337 gam_xy = (m1*m2*kx*ky).sum() 338 gam_yy = (m1*m2*ky*ky).sum() 339 340 denom = gam_xx*gam_yy - gam_xy*gam_xy 341 342 if denom == 0: 343 print "Warning: divide by zero error in gabor displacement. returning (0.0,0.0)" 344 return 0.,0. 345 346 else: 347 tmp1 = np.array([[gam_yy, -gam_xy],[-gam_xy,gam_xx]]) 348 tmp2 = np.array([[phi_x],[phi_y]]) 349 tmp = np.dot(tmp1,tmp2)/denom 350 351 return tmp[0,0],tmp[1,0]

352 353

354 - def displaceBlocked(self,gj,block_size=8):

355 m1 = self.mag 356 m2 = gj.mag 357 p1 = self.phase 358 p2 = gj.phase 359 kx = self.k[:,0] 360 ky = self.k[:,1] 361 k = self.k 362 d = np.array([[0.],[0.]]) 363 364 s = len(m1) 365 bs = block_size 366 nb = s/bs 367 assert len(m1)%bs == 0 368 369 for b in range(nb): 370 correction = np.dot(k,d).flatten() 371 low,high = s-(b+1)*bs,s 372 373 p = p1 - p2 - correction 374 375 mask = p > np.pi 376 while mask.max() == True: 377 p -= 2.0*np.pi * mask 378 mask = p > np.pi 379 380 mask = p < -np.pi 381 while mask.max() == True: 382 p += 2.0*np.pi * mask 383 mask = p < -np.pi 384 385 phi_x = (m1*m2*kx*p)[low:high].sum() 386 phi_y = (m1*m2*ky*p)[low:high].sum() 387 gam_xx = (m1*m2*kx*kx)[low:high].sum() 388 gam_xy = (m1*m2*kx*ky)[low:high].sum() 389 gam_yy = (m1*m2*ky*ky)[low:high].sum() 390 391 denom = gam_xx*gam_yy - gam_xy*gam_xy 392 393 if denom == 0: 394 print "Warning: divide by zero error in gabor displacement. returning (0.0,0.0)" 395 return 0.,0. 396 397 else: 398 tmp1 = np.array([[gam_yy, -gam_xy],[-gam_xy,gam_xx]]) 399 tmp2 = np.array([[phi_x],[phi_y]]) 400 tmp = np.dot(tmp1,tmp2)/denom 401 d += tmp 402 403 return d[0,0],d[1,0]

404 405

406 - def displaceIter(self,gj,N=8):

407 m1 = self.mag 408 m2 = gj.mag 409 p1 = self.phase 410 p2 = gj.phase 411 kx = self.k[:,0] 412 ky = self.k[:,1] 413 k = self.k 414 d = np.array([[0.],[0.]]) 415 416 for _ in range(N): 417 correction = np.dot(k,d).flatten() 418 419 p = p1 - p2 - correction 420 421 mask = p > np.pi 422 while mask.max() == True: 423 p -= 2.0*np.pi * mask 424 mask = p > np.pi 425 426 mask = p < -np.pi 427 while mask.max() == True: 428 p += 2.0*np.pi * mask 429 mask = p < -np.pi 430 431 phi_x = (m1*m2*kx*p).sum() 432 phi_y = (m1*m2*ky*p).sum() 433 gam_xx = (m1*m2*kx*kx).sum() 434 gam_xy = (m1*m2*kx*ky).sum() 435 gam_yy = (m1*m2*ky*ky).sum() 436 437 denom = gam_xx*gam_yy - gam_xy*gam_xy 438 439 if denom == 0: 440 print "Warning: divide by zero error in gabor displacement. returning (0.0,0.0)" 441 return 0.,0. 442 443 else: 444 tmp1 = np.array([[gam_yy, -gam_xy],[-gam_xy,gam_xx]]) 445 tmp2 = np.array([[phi_x],[phi_y]]) 446 tmp = np.dot(tmp1,tmp2)/denom 447 d += tmp 448 449 return d[0,0],d[1,0]

450 451

452 - def simMag(self,gj):

453 ''' 454 Magnitude similarity measure. 455 ''' 456 m1 = self.mag 457 m2 = gj.mag 458 459 tmp1 = (m1*m2).sum() 460 tmp2 = (m1*m1).sum() 461 tmp3 = (m2*m2).sum() 462 463 return tmp1/np.sqrt(tmp2*tmp3)

464 465

466 - def simPhase(self,gj):

467 ''' 468 Magnitude similarity measure. 469 ''' 470 m1 = self.mag 471 m2 = gj.mag 472 p1 = self.phase 473 p2 = gj.phase 474 475 tmp1 = (m1*m2*np.cos(p1 - p2)).sum() 476 tmp2 = (m1*m1).sum() 477 tmp3 = (m2*m2).sum() 478 479 return tmp1/np.sqrt(tmp2*tmp3)

480 481

482 - def simDisplace(self,gj,d=None):

483 ''' 484 Displacement similarity measure. 485 ''' 486 m1 = self.mag 487 m2 = gj.mag 488 p1 = self.phase 489 p2 = gj.phase 490 k = self.k 491 492 if d: 493 pass 494 else: 495 d = np.array(self.displace(gj)).reshape(2,1) 496 497 correction = np.dot(k,d).flatten() 498 499 tmp1 = (m1*m2*np.cos(p1 - p2 - correction)).sum() 500 tmp2 = (m1*m1).sum() 501 tmp3 = (m2*m2).sum() 502 503 return tmp1/np.sqrt(tmp2*tmp3)

504 505

506 -class _FastFilterTest(unittest.TestCase):

507 - def setUp(self):

508 SCRAPS_FACE_DATA = os.path.join(pv.__path__[0],"data","csuScrapShots") 509 self.test_images = [] 510 self.eyes = EyesFile(os.path.join(SCRAPS_FACE_DATA,"coords.txt")) 511 for filename in self.eyes.files()[0:10]: 512 im = pv.Image(os.path.join(SCRAPS_FACE_DATA, filename + ".pgm")) 513 eyes = self.eyes.getEyes(filename) 514 #print eyes 515 affine = pv.AffineFromPoints(eyes[0][0],eyes[0][1],pv.Point(40,40),pv.Point(88,40),(128,128)) 516 im = affine.transformImage(im) 517 518 self.test_images.append(im)

519 520

521 - def test_gabor1(self):

522 ilog = None # pv.ImageLog(name="GaborTest1") 523 524 bank = FilterBank(tile_size=(128,128)) 525 kernels = createGaborKernels() 526 for wavelet in kernels: 527 bank.addFilter(wavelet) 528 529 for i in range(len(bank.filters)): 530 corr_filter = np.fft.ifft2(bank.filters[i]) 531 if ilog: 532 ilog.log(pv.Image(np.fft.fftshift(corr_filter.real)),label="Filter_RE_%d"%i) 533 ilog.log(pv.Image(np.fft.fftshift(corr_filter.imag)),label="Filter_IM_%d"%i) 534 535 536 for im in self.test_images[:1]: 537 if ilog: 538 ilog.log(im,label="ORIG") 539 results = bank.convolve(im) 540 #print "RShape",results.shape[2] 541 if ilog: 542 for i in range(results.shape[2]): 543 ilog.log(pv.Image(results[:,:,i].real),label="CONV_RE_%d"%i) 544 ilog.log(pv.Image(results[:,:,i].imag),label="CONV_IM_%d"%i) 545 if ilog: 546 ilog.show()

547 548

549 - def test_gabor2(self):

550 #print 551 bank = FilterBank(tile_size=(128,128)) 552 #print "generating filters" 553 freq = 5 554 oreint = 3 555 test_values = [[0.39182228, 0.39265844, 5.4866541e-07, 7.4505806e-08, -0.00035403497, -0.049227916], [0.39180198, 0.3926785, -1.0430813e-07, 4.4703484e-08, 0.017195048, -0.045993667], [0.39180198, 0.39267856, -5.9604645e-08, 7.4505806e-09, 0.045639627, 0.017549083], [0.1959009, 0.19633935, -1.0430813e-07, -1.4901161e-08, -0.054863166, -0.010506729], [0.19590084, 0.19633923, -1.0337681e-07, -3.7252903e-08, -0.050385918, -0.024042567], [0.19590083, 0.19633923, -1.8626451e-07, -2.6077032e-08, 0.053237237, 0.014138865], [0.097950377, 0.098169744, 7.5995922e-07, 4.0978193e-08, -0.029739097, 0.029439665], [0.097950377, 0.098169573, 8.5681677e-08, 7.4505806e-09, -0.034587309, 0.023616293], [0.097950369, 0.098169573, 3.7252903e-09, 2.7939677e-08, 0.040645007, 0.0075459499], [0.048975188, 0.049084734, -1.0859221e-06, 1.4528632e-07, -0.0026100441, 0.025389811], [0.048975196, 0.049084827, -9.3132257e-09, 3.3527613e-08, -0.0058529032, 0.02486741], [0.048975192, 0.049084831, -7.6368451e-08, -4.703179e-08, 0.025110571, 0.0032778508], [0.024487557, 0.024542304, -0.00027022697, 1.1050142e-06, 0.0053004427, 0.013041494], [0.024487549, 0.024542348, -4.4152141e-05, 3.0389987e-05, 0.0040912526, 0.013478963], [0.024487551, 0.024542345, -4.4191256e-05, 3.0397903e-05, 0.013955923, 0.001284558]] 556 i = 0 557 for f in range(freq): 558 #f=f+2 559 for o in range(oreint): 560 #print "F:%10.5f O:%10.5f"%(f,o) 561 w1 = GaborWavelet(np.pi*2.0**(-(f+2.0)/2.0),o*np.pi/oreint,np.pi) 562 mask = w1.mask((64,64)) 563 ssr = (np.real(mask)**2).sum() 564 ssi = (np.imag(mask)**2).sum() 565 sr = np.real(mask).sum() 566 si = np.imag(mask).sum() 567 568 self.assertAlmostEqual(test_values[i][0],ssr,places=5) 569 self.assertAlmostEqual(test_values[i][1],ssi,places=5) 570 self.assertAlmostEqual(test_values[i][2],sr,places=5) 571 self.assertAlmostEqual(test_values[i][3],si,places=5) 572 self.assertAlmostEqual(test_values[i][4],np.real(mask)[3,2],places=5) 573 self.assertAlmostEqual(test_values[i][5],np.imag(mask)[3,2],places=5) 574 575 bank.addFilter(w1) 576 i+=1

577 #mask = w1.mask((128,128)) 578 #Image(real(mask)).show() 579 #Image(imag(mask)).show() 580 #print "Filter Time:",stop-start 581 582 583 584

585 - def test_GaborFilters(self):

586 587 #bank = FilterBank(tile_size=(128,128)) 588 kernels = createGaborKernels() 589 filters = GaborFilters(kernels) 590 591 gim = filters.convolve(self.test_images[0]) 592 593 template = gim.extractJet(pv.Point(64,64)) 594 595 table = pv.Table() 596 597 for i in range(0,128): 598 table[i,'disp'] = i - 64 599 novel = gim.extractJet(pv.Point(i,64)) 600 table[i,'simMag'] = template.simMag(novel) 601 table[i,'simPhase'] = template.simPhase(novel) 602 table[i,'simDisplace'] = template.simDisplace(novel) 603 dx,dy = template.displace(novel,method="Simple") 604 table[i,'displace_dx'] = dx 605 table[i,'displace_dy'] = dy 606 dx,dy = template.displace(novel,method="Blocked") 607 table[i,'blocked_dx'] = dx 608 table[i,'blocked_dy'] = dy 609 dx,dy = template.displace(novel,method="Iter") 610 table[i,'iter_dx'] = dx 611 table[i,'iter_dy'] = dy

612 613 #print 614 #print table 615 #table.save("../../gabor_plot.csv") 616

617 - def test_GaborImage(self):

618 619 kernels = createGaborKernels() 620 filters = GaborFilters(kernels) 621 622 gim = filters.convolve(self.test_images[0]) 623 624 test_point = pv.Point(62.6,64.8) 625 template = gim.extractJet(test_point) 626 627 new_point,_,_ = gim.locatePoint(template,pv.Point(60,70)) 628 self.assertAlmostEqual(new_point.l2(test_point),0.0) 629 630 new_point,_,_ = gim.locatePoint(template,pv.Point(30,49)) 631 self.assert_(new_point.l2(test_point) > 1.0) 632 633 new_point,_,_ = gim.locatePoint(template) 634 self.assertAlmostEqual(new_point.l2(test_point),0.0)

635

Source Code for Module pyvision.point.GaborJets