« back to all changes in this revision
Viewing changes to lib/sfop.py
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: tranx
- Date: 2010-10-01 16:56:14 UTC
- Revision ID: tranx@havane-20101001165614-u938mdd1y1fgd0o5
initial commit
- files added:
- demos
- lib
- lib/gb
- nodes
- wrappers
added
removed
lib/sfop.py
1
# -*- coding: utf-8 -*-
2
"""
3
SFOP :
4
5
6
Algorithm authors:
7
Wolfgang Förstner, Timo Dickscheid and Falko Schindler
8
Related Publication
9
Detecting Interpretable and Accurate Scale-Invariant Keypoints, International Conference on Computer Vision (ICCV'09), Kyoto, 200
10
11
Matlab Code Author:
12
Falko Schindler (falko.schindler@uni-bonn.de)
13
Department of Photogrammetry
14
Institute of Geodesy and Geoinformation
15
University of Bonn
16
Bonn, Germany
17
Copyright 2009-2009
18
19
Ported to Python from matlab by B.Nouvel(bertrand.nouvel@gmail.com)
20
21
22
23
24
25
Licence:
26
For research and academic use only.
27
No commercial application.
28
29
Warranty:
30
No warranty for validity of this implementation.
31
32
"""
33
34
import numpy, scipy,os
35
import scipy.stats
36
from scipy import ndimage
37
import pycvfext.keypoints.lib
38
39
##
40
## MATLAB Equivalents + Utility functions
41
##
42
43
44
def ind2sub(shp,ind):
45
return [ [ (i/p[1])%p[0] for i in ind] for p in zip(shp,reduce( lambda x,y:x+[x[-1]*y] ,shp[:-1],[1],) ) ]
46
47
def imregionalmax(I):
48
return (scipy.ndimage.maximum_filter(I,size=(3,3,3))==I) #*(1- ((scipy.ndimage.minimum_filter(I,size=(3,3,3))==I)))
49
50
def chi2inv(a,b):
51
return scipy.stats.gamma.cdf(a/2.0,b/2.0)
52
53
def rcond(A):
54
return 1/(numpy.linalg.norm(A,1)*numpy.linalg.norm(numpy.linalg.inv(A), 1))
55
56
def pad(img,p,v=1):
57
a=numpy.ones(map(lambda x,y:x+2*y,img.shape,p))*v
58
a[p[0]:-p[0],p[1]:-p[1]]=img
59
return a
60
61
def unpad(img,p):
62
if p[0]>0 or p[1]>0:
63
return img[p[0]:-p[0],p[1]:-p[1]]
64
else:
65
return img
66
67
def shiftdim(M,n):
68
strides=M.strides
69
shape=M.shape
70
strides=strides[n:]+strides[:n]
71
shape=shape[n:]+shape[:n]
72
return numpy.ndarray(buffer=M.copy('C').data,strides=strides,shape=shape,dtype=M.dtype)
73
74
75
76
LUT_time = numpy.load(os.path.join(pycvfext.keypoints.lib.__path__[0],'fftLUT_time.npy'));
77
78
79
80
81
gauss_kernels={
82
'G':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: numpy.outer(Hy, Hx)),
83
'Gx':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: numpy.outer(Hy, (Hx*dx))),
84
'Gxx':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: numpy.outer(Hy, (Hx*dxx))),
85
'Gy':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: numpy.outer((Hy*dy), Hx)),
86
'Gyy':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: numpy.outer((Hy*dyy), Hx)),
87
'y2G':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: numpy.outer((Hy*dyy), Hx)*sigma[0]**4+numpy.outer(Hy, Hx)*sigma[1]**2),
88
'x2G':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: numpy.outer(Hy, (Hx*dxx))*sigma[1]**4+numpy.outer(Hy, Hx)*sigma[1]**2),
89
'xyG':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: numpy.outer((Hy*dy), (Hx*dx))*sigma[0]*sigma[1]**2)
90
}
91
92
93
##
94
## These are separated versions for CUDA (later)
95
##
96
97
gauss_kernels_sep={
98
'G':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: (Hy, Hx)),
99
'Gx':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: (Hy, (Hx*dx))),
100
'Gxx':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: (Hy, (Hx*dxx))),
101
'Gy':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: ((Hy*dy), Hx)),
102
'Gyy':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: ((Hy*dyy), Hx)),
103
'y2G':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: ((Hy*dyy), Hx)*sigma[0]**4+numpy.outer(Hy, Hx)*sigma[1]**2),
104
'x2G':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: (Hy, (Hx*dxx))*sigma[1]**4+numpy.outer(Hy, Hx)*sigma[1]**2),
105
'xyG':(lambda Hy,Hx,dx,dy,dxx,dyy,sigma: ((Hy*dy), (Hx*dx))*sigma[0]*sigma[1]**2)
106
}
107
108
109
110
##
111
## CONSTANTS
112
##
113
114
D0=numpy.array([[16., 8., 16.], [8., 4., 8.],[ 16., 8., 16.]])**-1
115
D1=numpy.array([[16., 8., 16.], [numpy.inf, numpy.inf, numpy.inf],[ -16, -8, -16]])**-1
116
D33 = numpy.dstack([D0, -2 * D0, D0])
117
D11 = shiftdim(D33, 2);
118
D22 = shiftdim(D11, 2);
119
D13 = numpy.dstack([ D1, numpy.zeros((3,3)), -D1])
120
D12 = shiftdim(D13, 2);
121
D23 = shiftdim(D12, 2);
122
DA=numpy.array([ D11, D12, D13, D12, D22, D23, D13, D23, D33])
123
124
125
GD3 = numpy.dstack([D0, numpy.zeros((3,3)), -D0])
126
GD1 = shiftdim(GD3, 2)
127
GD2 = shiftdim(GD1, 2)
128
DGD=numpy.array([GD1,GD2,GD3])
129
130
##
131
## SPECIFIC IMPLEMENTATION UTILITY FUNC
132
##
133
134
135
def hessian(C):
136
dx=numpy.inner(DA.reshape(9,27) , C.reshape(1,27))
137
H = dx.reshape(( 3, 3))
138
return H
139
140
141
def gradient(C):
142
return numpy.inner(DGD.reshape((3,27)) , C.reshape((1,27)));
143
144
145
146
147
def cubicInterp3(ar,x,y,z):
148
if (x.shape[0]==0):
149
return numpy.array([])
150
#print ar.shape,ar.dtype
151
#print x.shape,x.min(),x.max()
152
#print y.shape,y.min(),y.max()
153
#print z.shape,z.min(),z.max()
154
coordinates=numpy.hstack([x.reshape((-1,1)),y.reshape((-1,1)),z.reshape((-1,1))])
155
#print coordinates.shape
156
#print numpy.isnan(ar).sum()
157
#assert(not numpy.isnan(ar).any())
158
#print coordinates,coordinates.dtype
159
r=ndimage.map_coordinates(ar.astype(numpy.float32),coordinates.T.astype(numpy.float32),mode='nearest')
160
#print r.shape
161
#print numpy.isnan(r).sum()
162
#assert(not numpy.isnan(r).any())
163
return r
164
165
166
def gaussFFT(img, sigma,ffttime, *kernels):
167
"""
168
% GAUSSFFT convolves an image IMG with Gaussian kernels via FFT.
169
%
170
% Kernels with SIGMA smaller than 0.4 pixel will be neglected since the
171
% algorithm gets very slow and thus not suited for small scales.
172
%
173
% Input:
174
% img image to be convolved
175
% sigma kernel parameter of length 1 or 2
176
% varargin first optional argument can be a look-up table for fft2
177
% computational costs in order to predict optimal image
178
% sizes (vector format),
179
% all other arguments are kernel names, possible names are:
180
% 'G',
181
% 'Gx', 'Gy',
182
% 'Gxx', 'Gyy', 'Gxy',
183
% 'x2G', 'y2G', 'xyG'
184
%
185
% Output:
186
% varargout convolved images, corresponding to one kernel each
187
%
188
% Example:
189
% img = zeros(512);
190
% img(150, 100) = 1;
191
% fftResult = gaussFFT(img, 20, 'G');
192
% imshow(fftResult, []);
193
%% use standard convolution
194
% varargout = cell(size(varargin));
195
% for i = 2 : numel(varargin)
196
% varargout{i - 1} = conv2(img, getFilter(sigma, varargin{i}), 'same');
197
% end
198
% return
199
"""
200
201
sigma = numpy.array([sigma, sigma])
202
sigma = sigma * (sigma >= 0.4); # < neglect small values
203
204
205
# pad and transform image
206
# image padding in spatial domain, depending on sigma
207
padding = (4,)*len(sigma)
208
# default size of padding
209
paddingv = map(lambda x,y:round(x*y), sigma , padding)
210
if paddingv[0]>0 or paddingv[1]>0:
211
f = pad(img,paddingv,1)
212
else:
213
f=img
214
215
216
# Fast Fourier Transformation of the image
217
F = numpy.fft.fft2(f);
218
N = F.shape
219
220
#% gaussians Hc for rows and columns separately
221
# transition function: gaussians
222
Hc=[None,None]
223
# build them for rows and columns separately
224
for dim in range(2):
225
#circular repetitions in frequency domain
226
if sigma[dim] > 0:
227
nShah = numpy.ceil(6 / sigma[dim])
228
else:
229
nShah = 0
230
231
# multiple repetitions of each pixel
232
j, n = numpy.mgrid[-nShah : (nShah+1), 0 : N[dim] ]
233
234
# gaussian value at each position
235
Hc[dim] = numpy.exp(-2 * (sigma[dim] * numpy.pi * (j + n / N[dim]))**2);
236
# sum over all repetitions
237
Hx_temp =Hc[1].sum(axis=0);
238
Hy_temp = Hc[0].sum(axis=0);
239
240
# normalization (divide by N, so it won't have to be done when using fft2 instead of ifft2)
241
Hy = Hy_temp / N[0]
242
Hx = Hx_temp / N[1]
243
#% construct kernels
244
# see handwritten notes
245
dy = - 1j * numpy.sin(2 * numpy.pi * numpy.arange(0,1,1./N[0]))
246
dx = 1j * numpy.sin(2 * numpy.pi * numpy.arange(0,1,1./N[1]))
247
dyy = - 2 * (1 - numpy.cos(2 * numpy.pi * numpy.arange(0,1,1./N[0])))
248
dxx = - 2 * (1 - numpy.cos(2 * numpy.pi * numpy.arange(0,1,1./N[1])))
249
250
#return map(lambda k:unpad(numpy.real(numpy.fft.ifft2((F * gauss_kernels[k](Hy,Hx,dx,dy,dxx,dyy,sigma)))),paddingv),kernels)
251
return map(lambda k:numpy.flipud(unpad(numpy.real(numpy.fft.fft2((F * gauss_kernels[k](Hy,Hx,dx,dy,dxx,dyy,sigma)).T).T),paddingv)),kernels)
252
253
254
255
256
def omega(alphaDeg, g, sig, M, fftTime):
257
"""
258
OMEGA computes the model error for a certain ALPHA on a gradient image G.
259
"""
260
261
# convert angle
262
alpha = alphaDeg / 180. * numpy.pi;
263
264
# rotate gradients
265
gAlpha_r = numpy.cos(alpha) * g[0] + numpy.sin(alpha) * g[1];
266
gAlpha_c = - numpy.sin(alpha) * g[0] + numpy.cos(alpha) * g[1];
267
268
# compute error messure
269
return M * ( gaussFFT(gAlpha_r * gAlpha_r, sig, fftTime, 'y2G')[0]
270
+ 2 * gaussFFT(gAlpha_r * gAlpha_c, sig, fftTime, 'xyG')[0]
271
+ gaussFFT(gAlpha_c * gAlpha_c, sig, fftTime, 'x2G')[0])
272
273
274
275
def buildScaleSpace(img, params):
276
# process each scale
277
L2, P = numpy.zeros((img.shape)+(params.nS,)), numpy.zeros(img.shape+(params.nS,))
278
for s in range(params.nS):
279
##
280
## FIRST We COMPUTE L2
281
##
282
283
# differentiation
284
gTau_r, gTau_c = gaussFFT(img, params.tau(s), LUT_time, 'Gy', 'Gx')
285
M = 12 * params.sig(s)**2 + 1
286
287
##
288
## smallest eigenvalue of the structure tensor
289
##
290
291
N0_rr = gaussFFT(gTau_r * gTau_r, params.sig(s), LUT_time, 'G')[0]*M
292
N0_cc = gaussFFT(gTau_c * gTau_c, params.sig(s), LUT_time, 'G')[0]*M
293
N0_rc = gaussFFT(gTau_r * gTau_c, params.sig(s), LUT_time, 'G')[0]*M
294
295
tr = N0_rr + N0_cc
296
det = N0_rr * N0_cc - N0_rc**2
297
L2[:, :, s] = numpy.real(0.5 * (tr - numpy.sqrt(tr**2 - 4 * det)))
298
299
300
##
301
## NOW WE COMPUTE P
302
##
303
# precision
304
if params._type== 'min':
305
S0 = omega( 0, (gTau_r,gTau_c), params.sig(s), M, LUT_time)
306
S1 = omega( 60, (gTau_r,gTau_c), params.sig(s), M, LUT_time)
307
S2 = omega(120, (gTau_r,gTau_c), params.sig(s), M, LUT_time)
308
# model: S = a - b * cos(2 * alpha - 2 * alpha0)
309
a = (S0 + S1 + S2) / 3;
310
b = 2 / 3 * numpy.sqrt( S0**2 + S1**2 + S2**2 - S0 * S1 - S1 * S2 - S2 * S0);
311
# alpha0 = 1 / 2 * atan2(sqrt(3) * (S2 - S1), S1 + S2 - 2 * S0);
312
S = a - b
313
R = M - 3
314
P[:, :, s] = L2[:, :, s] * R / S
315
elif type(params._type) in [float,int]:
316
S = omega(params._type, (gTau_r,gTau_c), params.sig(s), M, LUT_time)
317
R = M - 2
318
P[:, :, s] = L2[:, :, s] * R / S
319
else:
320
raise ValueError,('Invalid scale \"type\". Must be \"min\" or \"numeric\".');
321
322
# remove imaginary and undefined precisions
323
P = numpy.real(P);
324
P[numpy.isnan(P)] = -numpy.inf;
325
return P,L2
326
327
328
329
330
331
#######################################################################################################################
332
### PART II KEYPOINT DETECTION
333
########################################################################################################################
334
335
336
def nonMax(P, r, c, s, params):
337
"""
338
NONMAX does a non-maximum suppression at (R, C, S) w. r. t. precision P.
339
"""
340
341
table = numpy.hstack([r.reshape(-1,1), c.reshape(-1,1), s.reshape(-1,1), cubicInterp3(P, r, c, s).reshape(-1,1)]) #, (len(r), 4));
342
if table.shape[0]==0:
343
return table.T
344
idx = numpy.argsort(table[:,3])[::-1]
345
table = table[idx]
346
347
for line1 in range(table.shape[0]):
348
# compare to following points (with lower precision)
349
line2=numpy.arange(line1 + 1, table.shape[0])
350
351
# extract sigmas
352
sig1 = params.sig(table[line1, 2])
353
sig2 = params.sig(table[line2, 2])
354
355
#print sig1,sig2
356
357
# Mahalanobis distance
358
d2 = (((table[line1, 0] - table[line2, 0])**2 + (table[line1, 1] - table[line2, 1])**2) / (sig1**2 + sig2**2)
359
+ (table[line1, 2] - table[line2, 2])**2 / ((params.nonmaxOctave * params.layersPerOctave)**2))
360
361
# remove all points closer than threshold
362
table[line2[d2 < params.nonmaxTd2], :] = None
363
364
return table[:,0],table[:,1],table[:,2],table[:,3]
365
366
367
368
369
def detectPoints(P, L2, params,eps=1e-27):
370
"""
371
DETECTPOINTS extracts interest points from the scale space.
372
"""
373
374
#% local maxima "
375
376
377
r, c, s = imregionalmax(P).nonzero()
378
print 'Found %d local maxima, '%(len(r),)
379
380
#% remove points close to the boundary
381
## slice assigning wourl be better hre ?
382
border = 2.5 * params.sig(s)
383
border[border<2]=2
384
valid = (r - 1 >= border)
385
valid *= ((P.shape[0] - r) > border)
386
valid *=( c - 1 >= border)
387
valid *= ((P.shape[1] - c) > border )
388
valid *= (s - 1 >= 1)
389
valid *= ((P.shape[2] - s) > 1)
390
391
#print "valid=",valid,valid.sum(),valid.shape
392
393
r, c, s = r[valid], c[valid], s[valid];
394
395
print('%d off boundary, '% len(r))
396
397
#% threshold on p and on lambda
398
h = params.noise**2 / 16 / numpy.pi / params.tau(s)**4
399
400
401
402
Tl2 = h * params.Tlambda2 * chi2inv(0.999, 24 * params.sig(s)**2 + 2);
403
404
ciP=P[(r.reshape((-1,1)),c.reshape((-1,1)),s.reshape((-1,1)))] #P[numpy.hstack([r.reshape((-1,1)),c.reshape((-1,1)),s.reshape((-1,1))]) ]
405
ciL=L2[(r.reshape((-1,1)),c.reshape((-1,1)),s.reshape((-1,1)))]
406
print valid.shape, ciP.shape, ciL.shape,params.Tp,Tl2.shape
407
valid = numpy.asarray(ciP > params.Tp).ravel() * numpy.asarray(ciL > Tl2).ravel()
408
print valid.shape, valid.dtype
409
print ('%d above thresholds, '%( valid.sum(),))
410
411
#% optimize each point
412
R, C, S = numpy.mgrid[-1:2,-1:2,-1:2]
413
414
r=r.reshape(-1,1)
415
c=c.reshape(-1,1)
416
s=s.reshape(-1,1)
417
valid=valid.reshape(-1,1)
418
for p in valid.nonzero()[0]:
419
# iterate multiple times
420
#print p
421
for it in range( params.koetheMaxIter + 1):
422
# gradient and Hessian in current 27-neighborhood
423
#if it == 0:
424
# P27 = P(r[p] + R, c[p] + C, s[p] + S);
425
#else:
426
#print R.shape, r[p].shape,r.shape,p
427
#P27 = cubicInterp3(P, R+r[p], C+c[p] , S+s[p]);
428
try:
429
P27=P[( R+r[p], C+c[p] , S+s[p])]
430
except:
431
valid[p]=False
432
break
433
434
H = hessian(P27)
435
436
if (H.shape==0):
437
valid[p] = False;
438
break # update position, but check for image boundary
439
#print H
440
if numpy.isnan(H).any():
441
valid[p]=False
442
break
443
if (numpy.linalg.eig(H)[0] > 0).all():
444
#print "eig rej"
445
valid[p]=False
446
break
447
if (rcond(H) < eps):
448
print "rcond rej"
449
valid[p] = False;
450
break # update position, but check for image boundary
451
452
##
453
## adding bias to positions
454
##
455
#print H.shape
456
G=gradient(P27)
457
#print G.shape,G
458
bias = (numpy.asmatrix(H) * G**-1)#numpy.linalg.inv(numpy.asmatrix(G)))
459
#print bias
460
r[p] = r[p] + bias[0]
461
c[p] = c[p] + bias[1]
462
s[p] = s[p] + bias[2]
463
464
if min(r[p], c[p], s[p], P.shape[0] - r[p], P.shape[1]-c[p], P.shape[2]-s[p]) < 2:
465
valid[p] = False
466
break
467
468
# check convergence
469
if numpy.linalg.norm(bias) < params.koetheEpsilon:
470
break
471
472
if it == params.koetheMaxIter:
473
valid[p] = False
474
break;
475
476
print('%d converged'%( valid.sum()));
477
478
print('removing non maxima')
479
#% supress non-maxima
480
x = nonMax(P, r[valid], c[valid], s[valid], params);
481
#print x.shape
482
valid=(1-numpy.isnan(x[0])).nonzero()
483
#print valid,x[0].shape,valid.shape
484
r,c,s=x[0],x[1],x[2]
485
return r[valid],c[valid],s[valid]
486
487
488
489
490
##
491
## Ok we just create a class as in the matlab implementation to store the parameters
492
##
493
494
495
496
class Parameters(object):
497
layersPerOctave = 4
498
numberOfOctaves = 3
499
tau = staticmethod(lambda sigma: sigma / 3)
500
Tp = -numpy.inf
501
Tlambda2 = 2
502
noise = 0.02
503
_type = 'min'
504
#_type = 30
505
nonmaxTd2 = 1
506
nonmaxOctave = 0.5
507
sizeFactor = 1
508
koetheMaxIter = 5
509
koetheEpsilon = 0.2
510
511
def __init__(self):
512
self.sigmaFrom = 2**(1 + 1 / self.layersPerOctave)
513
514
def update(self):
515
sigStep = 2**(1. / self.layersPerOctave);
516
sigFrom = self.sigmaFrom / sigStep**2;
517
sigTo = self.sigmaFrom * sigStep**2 * 2**self.numberOfOctaves;
518
self.nS = int( numpy.floor(numpy.log(sigTo / sigFrom) / numpy.log(sigStep)));
519
self.sig = lambda s: sigFrom * sigStep**(s - 1);
520
self.otau=self.tau
521
self.tau = lambda s: self.otau(self.sig(s))
522
523
524
525
def sfop(img,**kwargs):
526
"""
527
Please see license usage restrictions.
528
529
kwargs may be used to change one of these default parameter
530
layersPerOctave = 4
531
numberOfOctaves = 3
532
tau = staticmethod(lambda sigma: sigma / 3)
533
Tp = -numpy.inf
534
Tlambda2 = 2
535
noise = 0.02
536
_type = 'min'
537
#_type = 30
538
nonmaxTd2 = 1
539
nonmaxOctave = 0.5
540
sizeFactor = 1
541
koetheMaxIter = 5
542
koetheEpsilon = 0.2
543
"""
544
if (img.ndim==3):
545
img = img.mean(axis=2)
546
img=img.astype(numpy.float32)
547
img/=255.
548
549
params=Parameters()
550
for k in kwargs.items():
551
setattr(params,k[0],k[1])
552
params.update()
553
554
555
556
P, L2 = buildScaleSpace(img, params);
557
r, c, s = detectPoints(P, L2, params);
558
559
560
return numpy.hstack([r.reshape(-1,1),c.reshape(-1,1),s.reshape(-1,1)])
561
562
563
564
565
if __name__=="__main__":
566
#print sfop(scipy.lena()[:256,:512])
567
print sfop(scipy.lena()[:512,:512])
568
#print "done"
Loggerhead is a web-based interface for Breezy
Version: 2.0.1