Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Merge pull request #1 from jdagilliland/master

Repo cleanup
  • Loading branch information...
commit 4b0ecfeeaa12c8b5bab71377444987e00df6a4d5 2 parents f41d89d + 7c45690
@jesolem jesolem authored
Showing with 3,487 additions and 3,458 deletions.
  1. +1 −0  .gitignore
  2. BIN  PCV/__init__.pyc
  3. BIN  PCV/classifiers/__init__.pyc
  4. +51 −51 PCV/classifiers/bayes.py
  5. +27 −27 PCV/classifiers/knn.py
  6. BIN  PCV/clustering/__init__.pyc
  7. +121 −121 PCV/clustering/hcluster.py
  8. BIN  PCV/geometry/__init__.pyc
  9. +61 −61 PCV/geometry/camera.py
  10. +143 −143 PCV/geometry/homography.py
  11. +191 −191 PCV/geometry/sfm.py
  12. +129 −129 PCV/geometry/warp.py
  13. BIN  PCV/imagesearch/__init__.pyc
  14. +167 −167 PCV/imagesearch/imagesearch.py
  15. BIN  PCV/imagesearch/imagesearch.pyc
  16. +51 −51 PCV/imagesearch/vocabulary.py
  17. BIN  PCV/localdescriptors/__init__.pyc
  18. +28 −28 PCV/localdescriptors/dsift.py
  19. +131 −131 PCV/localdescriptors/harris.py
  20. +98 −98 PCV/localdescriptors/sift.py
  21. BIN  PCV/tools/__init__.pyc
  22. +91 −91 PCV/tools/graphcut.py
  23. +84 −84 PCV/tools/imregistration.py
  24. +66 −66 PCV/tools/imtools.py
  25. +50 −50 PCV/tools/ncut.py
  26. +42 −42 PCV/tools/pca.py
  27. +42 −42 PCV/tools/rof.py
  28. BIN  PCV/tools/rof.pyc
  29. +96 −96 examples/lktrack.py
  30. +60 −60 examples/searchdemo.py
  31. +28 −0 expand_tabs.py
  32. +51 −51 pcv_book/bayes.py
  33. +61 −61 pcv_book/camera.py
  34. +28 −28 pcv_book/dsift.py
  35. +91 −91 pcv_book/graphcut.py
  36. +131 −131 pcv_book/harris.py
  37. +121 −121 pcv_book/hcluster.py
  38. +143 −143 pcv_book/homography.py
  39. +167 −167 pcv_book/imagesearch.py
  40. +84 −84 pcv_book/imregistration.py
  41. +66 −66 pcv_book/imtools.py
  42. +27 −27 pcv_book/knn.py
  43. +96 −96 pcv_book/lktrack.py
  44. +50 −50 pcv_book/ncut.py
  45. +42 −42 pcv_book/pca.py
  46. +42 −42 pcv_book/rof.py
  47. +60 −60 pcv_book/searchdemo.py
  48. +191 −191 pcv_book/sfm.py
  49. +98 −98 pcv_book/sift.py
  50. +51 −51 pcv_book/vocabulary.py
  51. +129 −129 pcv_book/warp.py
View
1  .gitignore
@@ -0,0 +1 @@
+*.pyc
View
BIN  PCV/__init__.pyc
Binary file not shown
View
BIN  PCV/classifiers/__init__.pyc
Binary file not shown
View
102 PCV/classifiers/bayes.py
@@ -3,61 +3,61 @@
class BayesClassifier(object):
- def __init__(self):
- """ Initialize classifier with training data. """
-
- self.labels = [] # class labels
- self.mean = [] # class mean
- self.var = [] # class variances
- self.n = 0 # nbr of classes
-
- def train(self,data,labels=None):
- """ Train on data (list of arrays n*dim).
- Labels are optional, default is 0...n-1. """
-
- if labels==None:
- labels = range(len(data))
- self.labels = labels
- self.n = len(labels)
-
- for c in data:
- self.mean.append(mean(c,axis=0))
- self.var.append(var(c,axis=0))
-
- def classify(self,points):
- """ Classify the points by computing probabilities
- for each class and return most probable label. """
+ def __init__(self):
+ """ Initialize classifier with training data. """
- # compute probabilities for each class
- est_prob = array([gauss(m,v,points) for m,v in zip(self.mean,self.var)])
-
- print 'est prob',est_prob.shape,self.labels
- # get index of highest probability, this gives class label
- ndx = est_prob.argmax(axis=0)
-
- est_labels = array([self.labels[n] for n in ndx])
-
- return est_labels, est_prob
+ self.labels = [] # class labels
+ self.mean = [] # class mean
+ self.var = [] # class variances
+ self.n = 0 # nbr of classes
+
+ def train(self,data,labels=None):
+ """ Train on data (list of arrays n*dim).
+ Labels are optional, default is 0...n-1. """
+
+ if labels==None:
+ labels = range(len(data))
+ self.labels = labels
+ self.n = len(labels)
+
+ for c in data:
+ self.mean.append(mean(c,axis=0))
+ self.var.append(var(c,axis=0))
+
+ def classify(self,points):
+ """ Classify the points by computing probabilities
+ for each class and return most probable label. """
+
+ # compute probabilities for each class
+ est_prob = array([gauss(m,v,points) for m,v in zip(self.mean,self.var)])
+
+ print 'est prob',est_prob.shape,self.labels
+ # get index of highest probability, this gives class label
+ ndx = est_prob.argmax(axis=0)
+
+ est_labels = array([self.labels[n] for n in ndx])
+
+ return est_labels, est_prob
def gauss(m,v,x):
- """ Evaluate Gaussian in d-dimensions with independent
- mean m and variance v at the points in (the rows of) x.
- http://en.wikipedia.org/wiki/Multivariate_normal_distribution """
-
- if len(x.shape)==1:
- n,d = 1,x.shape[0]
- else:
- n,d = x.shape
-
- # covariance matrix, subtract mean
- S = diag(1/v)
- x = x-m
- # product of probabilities
- y = exp(-0.5*diag(dot(x,dot(S,x.T))))
-
- # normalize and return
- return y * (2*pi)**(-d/2.0) / ( sqrt(prod(v)) + 1e-6)
+ """ Evaluate Gaussian in d-dimensions with independent
+ mean m and variance v at the points in (the rows of) x.
+ http://en.wikipedia.org/wiki/Multivariate_normal_distribution """
+
+ if len(x.shape)==1:
+ n,d = 1,x.shape[0]
+ else:
+ n,d = x.shape
+
+ # covariance matrix, subtract mean
+ S = diag(1/v)
+ x = x-m
+ # product of probabilities
+ y = exp(-0.5*diag(dot(x,dot(S,x.T))))
+
+ # normalize and return
+ return y * (2*pi)**(-d/2.0) / ( sqrt(prod(v)) + 1e-6)
View
54 PCV/classifiers/knn.py
@@ -1,35 +1,35 @@
from numpy import *
class KnnClassifier(object):
-
- def __init__(self,labels,samples):
- """ Initialize classifier with training data. """
-
- self.labels = labels
- self.samples = samples
-
- def classify(self,point,k=3):
- """ Classify a point against k nearest
- in the training data, return label. """
-
- # compute distance to all training points
- dist = array([L2dist(point,s) for s in self.samples])
-
- # sort them
- ndx = dist.argsort()
-
- # use dictionary to store the k nearest
- votes = {}
- for i in range(k):
- label = self.labels[ndx[i]]
- votes.setdefault(label,0)
- votes[label] += 1
-
- return max(votes)
+
+ def __init__(self,labels,samples):
+ """ Initialize classifier with training data. """
+
+ self.labels = labels
+ self.samples = samples
+
+ def classify(self,point,k=3):
+ """ Classify a point against k nearest
+ in the training data, return label. """
+
+ # compute distance to all training points
+ dist = array([L2dist(point,s) for s in self.samples])
+
+ # sort them
+ ndx = dist.argsort()
+
+ # use dictionary to store the k nearest
+ votes = {}
+ for i in range(k):
+ label = self.labels[ndx[i]]
+ votes.setdefault(label,0)
+ votes[label] += 1
+
+ return max(votes)
def L2dist(p1,p2):
- return sqrt( sum( (p1-p2)**2) )
+ return sqrt( sum( (p1-p2)**2) )
def L1dist(v1,v2):
- return sum(abs(v1-v2))
View
BIN  PCV/clustering/__init__.pyc
Binary file not shown
View
242 PCV/clustering/hcluster.py
@@ -3,144 +3,144 @@
class ClusterNode(object):
- def __init__(self,vec,left,right,distance=0.0,count=1):
- self.left = left
- self.right = right
- self.vec = vec
- self.distance = distance
- self.count = count # only used for weighted average
-
- def extract_clusters(self,dist):
- """ Extract list of sub-tree clusters from
- hcluster tree with distance<dist. """
- if self.distance < dist:
- return [self]
- return self.left.extract_clusters(dist) + self.right.extract_clusters(dist)
-
- def get_cluster_elements(self):
- """ Return ids for elements in a cluster sub-tree. """
- return self.left.get_cluster_elements() + self.right.get_cluster_elements()
-
- def get_height(self):
- """ Return the height of a node,
- height is sum of each branch. """
- return self.left.get_height() + self.right.get_height()
-
- def get_depth(self):
- """ Return the depth of a node, depth is
- max of each child plus own distance. """
- return max(self.left.get_depth(), self.right.get_depth()) + self.distance
-
- def draw(self,draw,x,y,s,imlist,im):
- """ Draw nodes recursively with image
- thumbnails for leaf nodes. """
-
- h1 = int(self.left.get_height()*20 / 2)
- h2 = int(self.right.get_height()*20 /2)
- top = y-(h1+h2)
- bottom = y+(h1+h2)
-
- # vertical line to children
- draw.line((x,top+h1,x,bottom-h2),fill=(0,0,0))
-
- # horizontal lines
- ll = self.distance*s
- draw.line((x,top+h1,x+ll,top+h1),fill=(0,0,0))
- draw.line((x,bottom-h2,x+ll,bottom-h2),fill=(0,0,0))
-
- # draw left and right child nodes recursively
- self.left.draw(draw,x+ll,top+h1,s,imlist,im)
- self.right.draw(draw,x+ll,bottom-h2,s,imlist,im)
-
+ def __init__(self,vec,left,right,distance=0.0,count=1):
+ self.left = left
+ self.right = right
+ self.vec = vec
+ self.distance = distance
+ self.count = count # only used for weighted average
+
+ def extract_clusters(self,dist):
+ """ Extract list of sub-tree clusters from
+ hcluster tree with distance<dist. """
+ if self.distance < dist:
+ return [self]
+ return self.left.extract_clusters(dist) + self.right.extract_clusters(dist)
+
+ def get_cluster_elements(self):
+ """ Return ids for elements in a cluster sub-tree. """
+ return self.left.get_cluster_elements() + self.right.get_cluster_elements()
+
+ def get_height(self):
+ """ Return the height of a node,
+ height is sum of each branch. """
+ return self.left.get_height() + self.right.get_height()
+
+ def get_depth(self):
+ """ Return the depth of a node, depth is
+ max of each child plus own distance. """
+ return max(self.left.get_depth(), self.right.get_depth()) + self.distance
+
+ def draw(self,draw,x,y,s,imlist,im):
+ """ Draw nodes recursively with image
+ thumbnails for leaf nodes. """
+
+ h1 = int(self.left.get_height()*20 / 2)
+ h2 = int(self.right.get_height()*20 /2)
+ top = y-(h1+h2)
+ bottom = y+(h1+h2)
+
+ # vertical line to children
+ draw.line((x,top+h1,x,bottom-h2),fill=(0,0,0))
+
+ # horizontal lines
+ ll = self.distance*s
+ draw.line((x,top+h1,x+ll,top+h1),fill=(0,0,0))
+ draw.line((x,bottom-h2,x+ll,bottom-h2),fill=(0,0,0))
+
+ # draw left and right child nodes recursively
+ self.left.draw(draw,x+ll,top+h1,s,imlist,im)
+ self.right.draw(draw,x+ll,bottom-h2,s,imlist,im)
+
class ClusterLeafNode(object):
- def __init__(self,vec,id):
- self.vec = vec
- self.id = id
+ def __init__(self,vec,id):
+ self.vec = vec
+ self.id = id
- def extract_clusters(self,dist):
- return [self]
+ def extract_clusters(self,dist):
+ return [self]
- def get_cluster_elements(self):
- return [self.id]
+ def get_cluster_elements(self):
+ return [self.id]
- def get_height(self):
- return 1
+ def get_height(self):
+ return 1
- def get_depth(self):
- return 0
-
- def draw(self,draw,x,y,s,imlist,im):
- nodeim = Image.open(imlist[self.id])
- nodeim.thumbnail([20,20])
- ns = nodeim.size
- im.paste(nodeim,[int(x),int(y-ns[1]//2),int(x+ns[0]),int(y+ns[1]-ns[1]//2)])
+ def get_depth(self):
+ return 0
+
+ def draw(self,draw,x,y,s,imlist,im):
+ nodeim = Image.open(imlist[self.id])
+ nodeim.thumbnail([20,20])
+ ns = nodeim.size
+ im.paste(nodeim,[int(x),int(y-ns[1]//2),int(x+ns[0]),int(y+ns[1]-ns[1]//2)])
def L2dist(v1,v2):
- return sqrt(sum((v1-v2)**2))
+ return sqrt(sum((v1-v2)**2))
def L1dist(v1,v2):
- return sum(abs(v1-v2))
+ return sum(abs(v1-v2))
def hcluster(features,distfcn=L2dist):
- """ Cluster the rows of features using
- hierarchical clustering. """
-
- # cache of distance calculations
- distances = {}
-
- # initialize with each row as a cluster
- node = [ClusterLeafNode(array(f),id=i) for i,f in enumerate(features)]
-
- while len(node)>1:
- closest = float('Inf')
-
- # loop through every pair looking for the smallest distance
- for ni,nj in combinations(node,2):
- if (ni,nj) not in distances:
- distances[ni,nj] = distfcn(ni.vec,nj.vec)
-
- d = distances[ni,nj]
- if d<closest:
- closest = d
- lowestpair = (ni,nj)
- ni,nj = lowestpair
-
- # average the two clusters
- new_vec = (ni.vec + nj.vec) / 2.0
-
- # create new node
- new_node = ClusterNode(new_vec,left=ni,right=nj,distance=closest)
- node.remove(ni)
- node.remove(nj)
- node.append(new_node)
-
- return node[0]
+ """ Cluster the rows of features using
+ hierarchical clustering. """
+
+ # cache of distance calculations
+ distances = {}
+
+ # initialize with each row as a cluster
+ node = [ClusterLeafNode(array(f),id=i) for i,f in enumerate(features)]
+
+ while len(node)>1:
+ closest = float('Inf')
+
+ # loop through every pair looking for the smallest distance
+ for ni,nj in combinations(node,2):
+ if (ni,nj) not in distances:
+ distances[ni,nj] = distfcn(ni.vec,nj.vec)
+
+ d = distances[ni,nj]
+ if d<closest:
+ closest = d
+ lowestpair = (ni,nj)
+ ni,nj = lowestpair
+
+ # average the two clusters
+ new_vec = (ni.vec + nj.vec) / 2.0
+
+ # create new node
+ new_node = ClusterNode(new_vec,left=ni,right=nj,distance=closest)
+ node.remove(ni)
+ node.remove(nj)
+ node.append(new_node)
+
+ return node[0]
from PIL import Image,ImageDraw
def draw_dendrogram(node,imlist,filename='clusters.jpg'):
- """ Draw a cluster dendrogram and save to a file. """
-
- # height and width
- rows = node.get_height()*20
- cols = 1200
-
- # scale factor for distances to fit image width
- s = float(cols-150)/node.get_depth()
-
- # create image and draw object
- im = Image.new('RGB',(cols,rows),(255,255,255))
- draw = ImageDraw.Draw(im)
-
- # initial line for start of tree
- draw.line((0,rows/2,20,rows/2),fill=(0,0,0))
-
- # draw the nodes recursively
- node.draw(draw,20,(rows/2),s,imlist,im)
- im.save(filename)
- im.show()
+ """ Draw a cluster dendrogram and save to a file. """
+
+ # height and width
+ rows = node.get_height()*20
+ cols = 1200
+
+ # scale factor for distances to fit image width
+ s = float(cols-150)/node.get_depth()
+
+ # create image and draw object
+ im = Image.new('RGB',(cols,rows),(255,255,255))
+ draw = ImageDraw.Draw(im)
+
+ # initial line for start of tree
+ draw.line((0,rows/2,20,rows/2),fill=(0,0,0))
+
+ # draw the nodes recursively
+ node.draw(draw,20,(rows/2),s,imlist,im)
+ im.save(filename)
+ im.show()
View
BIN  PCV/geometry/__init__.pyc
Binary file not shown
View
122 PCV/geometry/camera.py
@@ -3,72 +3,72 @@
class Camera(object):
- """ Class for representing pin-hole cameras. """
-
- def __init__(self,P):
- """ Initialize P = K[R|t] camera model. """
- self.P = P
- self.K = None # calibration matrix
- self.R = None # rotation
- self.t = None # translation
- self.c = None # camera center
-
-
- def project(self,X):
- """ Project points in X (4*n array) and normalize coordinates. """
-
- x = dot(self.P,X)
- for i in range(3):
- x[i] /= x[2]
- return x
-
-
- def factor(self):
- """ Factorize the camera matrix into K,R,t as P = K[R|t]. """
-
- # factor first 3*3 part
- K,R = linalg.rq(self.P[:,:3])
-
- # make diagonal of K positive
- T = diag(sign(diag(K)))
- if linalg.det(T) < 0:
- T[1,1] *= -1
-
- self.K = dot(K,T)
- self.R = dot(T,R) # T is its own inverse
- self.t = dot(linalg.inv(self.K),self.P[:,3])
-
- return self.K, self.R, self.t
-
-
- def center(self):
- """ Compute and return the camera center. """
+ """ Class for representing pin-hole cameras. """
- if self.c is not None:
- return self.c
- else:
- # compute c by factoring
- self.factor()
- self.c = -dot(self.R.T,self.t)
- return self.c
+ def __init__(self,P):
+ """ Initialize P = K[R|t] camera model. """
+ self.P = P
+ self.K = None # calibration matrix
+ self.R = None # rotation
+ self.t = None # translation
+ self.c = None # camera center
+
+
+ def project(self,X):
+ """ Project points in X (4*n array) and normalize coordinates. """
+
+ x = dot(self.P,X)
+ for i in range(3):
+ x[i] /= x[2]
+ return x
+
+
+ def factor(self):
+ """ Factorize the camera matrix into K,R,t as P = K[R|t]. """
+
+ # factor first 3*3 part
+ K,R = linalg.rq(self.P[:,:3])
+
+ # make diagonal of K positive
+ T = diag(sign(diag(K)))
+ if linalg.det(T) < 0:
+ T[1,1] *= -1
+
+ self.K = dot(K,T)
+ self.R = dot(T,R) # T is its own inverse
+ self.t = dot(linalg.inv(self.K),self.P[:,3])
+
+ return self.K, self.R, self.t
+
+
+ def center(self):
+ """ Compute and return the camera center. """
+
+ if self.c is not None:
+ return self.c
+ else:
+ # compute c by factoring
+ self.factor()
+ self.c = -dot(self.R.T,self.t)
+ return self.c
-# helper functions
+# helper functions
def rotation_matrix(a):
- """ Creates a 3D rotation matrix for rotation
- around the axis of the vector a. """
- R = eye(4)
- R[:3,:3] = linalg.expm([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])
- return R
-
+ """ Creates a 3D rotation matrix for rotation
+ around the axis of the vector a. """
+ R = eye(4)
+ R[:3,:3] = linalg.expm([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])
+ return R
+
def rq(A):
- from scipy.linalg import qr
-
- Q,R = qr(flipud(A).T)
- R = flipud(R.T)
- Q = Q.T
-
- return R[:,::-1],Q[::-1,:]
+ from scipy.linalg import qr
+
+ Q,R = qr(flipud(A).T)
+ R = flipud(R.T)
+ Q = Q.T
+
+ return R[:,::-1],Q[::-1,:]
View
286 PCV/geometry/homography.py
@@ -3,161 +3,161 @@
class RansacModel(object):
- """ Class for testing homography fit with ransac.py from
- http://www.scipy.org/Cookbook/RANSAC"""
-
- def __init__(self,debug=False):
- self.debug = debug
-
- def fit(self, data):
- """ Fit homography to four selected correspondences. """
-
- # transpose to fit H_from_points()
- data = data.T
-
- # from points
- fp = data[:3,:4]
- # target points
- tp = data[3:,:4]
-
- # fit homography and return
- return H_from_points(fp,tp)
-
- def get_error( self, data, H):
- """ Apply homography to all correspondences,
- return error for each transformed point. """
-
- data = data.T
-
- # from points
- fp = data[:3]
- # target points
- tp = data[3:]
-
- # transform fp
- fp_transformed = dot(H,fp)
-
- # normalize hom. coordinates
- fp_transformed = normalize(fp_transformed)
-
- # return error per point
- return sqrt( sum((tp-fp_transformed)**2,axis=0) )
-
+ """ Class for testing homography fit with ransac.py from
+ http://www.scipy.org/Cookbook/RANSAC"""
+
+ def __init__(self,debug=False):
+ self.debug = debug
+
+ def fit(self, data):
+ """ Fit homography to four selected correspondences. """
+
+ # transpose to fit H_from_points()
+ data = data.T
+
+ # from points
+ fp = data[:3,:4]
+ # target points
+ tp = data[3:,:4]
+
+ # fit homography and return
+ return H_from_points(fp,tp)
+
+ def get_error( self, data, H):
+ """ Apply homography to all correspondences,
+ return error for each transformed point. """
+
+ data = data.T
+
+ # from points
+ fp = data[:3]
+ # target points
+ tp = data[3:]
+
+ # transform fp
+ fp_transformed = dot(H,fp)
+
+ # normalize hom. coordinates
+ fp_transformed = normalize(fp_transformed)
+
+ # return error per point
+ return sqrt( sum((tp-fp_transformed)**2,axis=0) )
+
def H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):
- """ Robust estimation of homography H from point
- correspondences using RANSAC (ransac.py from
- http://www.scipy.org/Cookbook/RANSAC).
-
- input: fp,tp (3*n arrays) points in hom. coordinates. """
-
- import ransac
-
- # group corresponding points
- data = vstack((fp,tp))
-
- # compute H and return
- H,ransac_data = ransac.ransac(data.T,model,4,maxiter,match_theshold,10,return_all=True)
- return H,ransac_data['inliers']
+ """ Robust estimation of homography H from point
+ correspondences using RANSAC (ransac.py from
+ http://www.scipy.org/Cookbook/RANSAC).
+
+ input: fp,tp (3*n arrays) points in hom. coordinates. """
+
+ import ransac
+
+ # group corresponding points
+ data = vstack((fp,tp))
+
+ # compute H and return
+ H,ransac_data = ransac.ransac(data.T,model,4,maxiter,match_theshold,10,return_all=True)
+ return H,ransac_data['inliers']
def H_from_points(fp,tp):
- """ Find homography H, such that fp is mapped to tp
- using the linear DLT method. Points are conditioned
- automatically. """
-
- if fp.shape != tp.shape:
- raise RuntimeError('number of points do not match')
-
- # condition points (important for numerical reasons)
- # --from points--
- m = mean(fp[:2], axis=1)
- maxstd = max(std(fp[:2], axis=1)) + 1e-9
- C1 = diag([1/maxstd, 1/maxstd, 1])
- C1[0][2] = -m[0]/maxstd
- C1[1][2] = -m[1]/maxstd
- fp = dot(C1,fp)
-
- # --to points--
- m = mean(tp[:2], axis=1)
- maxstd = max(std(tp[:2], axis=1)) + 1e-9
- C2 = diag([1/maxstd, 1/maxstd, 1])
- C2[0][2] = -m[0]/maxstd
- C2[1][2] = -m[1]/maxstd
- tp = dot(C2,tp)
-
- # create matrix for linear method, 2 rows for each correspondence pair
- nbr_correspondences = fp.shape[1]
- A = zeros((2*nbr_correspondences,9))
- for i in range(nbr_correspondences):
- A[2*i] = [-fp[0][i],-fp[1][i],-1,0,0,0,
- tp[0][i]*fp[0][i],tp[0][i]*fp[1][i],tp[0][i]]
- A[2*i+1] = [0,0,0,-fp[0][i],-fp[1][i],-1,
- tp[1][i]*fp[0][i],tp[1][i]*fp[1][i],tp[1][i]]
-
- U,S,V = linalg.svd(A)
- H = V[8].reshape((3,3))
-
- # decondition
- H = dot(linalg.inv(C2),dot(H,C1))
-
- # normalize and return
- return H / H[2,2]
+ """ Find homography H, such that fp is mapped to tp
+ using the linear DLT method. Points are conditioned
+ automatically. """
+
+ if fp.shape != tp.shape:
+ raise RuntimeError('number of points do not match')
+
+ # condition points (important for numerical reasons)
+ # --from points--
+ m = mean(fp[:2], axis=1)
+ maxstd = max(std(fp[:2], axis=1)) + 1e-9
+ C1 = diag([1/maxstd, 1/maxstd, 1])
+ C1[0][2] = -m[0]/maxstd
+ C1[1][2] = -m[1]/maxstd
+ fp = dot(C1,fp)
+
+ # --to points--
+ m = mean(tp[:2], axis=1)
+ maxstd = max(std(tp[:2], axis=1)) + 1e-9
+ C2 = diag([1/maxstd, 1/maxstd, 1])
+ C2[0][2] = -m[0]/maxstd
+ C2[1][2] = -m[1]/maxstd
+ tp = dot(C2,tp)
+
+ # create matrix for linear method, 2 rows for each correspondence pair
+ nbr_correspondences = fp.shape[1]
+ A = zeros((2*nbr_correspondences,9))
+ for i in range(nbr_correspondences):
+ A[2*i] = [-fp[0][i],-fp[1][i],-1,0,0,0,
+ tp[0][i]*fp[0][i],tp[0][i]*fp[1][i],tp[0][i]]
+ A[2*i+1] = [0,0,0,-fp[0][i],-fp[1][i],-1,
+ tp[1][i]*fp[0][i],tp[1][i]*fp[1][i],tp[1][i]]
+
+ U,S,V = linalg.svd(A)
+ H = V[8].reshape((3,3))
+
+ # decondition
+ H = dot(linalg.inv(C2),dot(H,C1))
+
+ # normalize and return
+ return H / H[2,2]
def Haffine_from_points(fp,tp):
- """ Find H, affine transformation, such that
- tp is affine transf of fp. """
-
- if fp.shape != tp.shape:
- raise RuntimeError('number of points do not match')
-
- # condition points
- # --from points--
- m = mean(fp[:2], axis=1)
- maxstd = max(std(fp[:2], axis=1)) + 1e-9
- C1 = diag([1/maxstd, 1/maxstd, 1])
- C1[0][2] = -m[0]/maxstd
- C1[1][2] = -m[1]/maxstd
- fp_cond = dot(C1,fp)
-
- # --to points--
- m = mean(tp[:2], axis=1)
- C2 = C1.copy() #must use same scaling for both point sets
- C2[0][2] = -m[0]/maxstd
- C2[1][2] = -m[1]/maxstd
- tp_cond = dot(C2,tp)
-
- # conditioned points have mean zero, so translation is zero
- A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
- U,S,V = linalg.svd(A.T)
-
- # create B and C matrices as Hartley-Zisserman (2:nd ed) p 130.
- tmp = V[:2].T
- B = tmp[:2]
- C = tmp[2:4]
-
- tmp2 = concatenate((dot(C,linalg.pinv(B)),zeros((2,1))), axis=1)
- H = vstack((tmp2,[0,0,1]))
-
- # decondition
- H = dot(linalg.inv(C2),dot(H,C1))
-
- return H / H[2,2]
+ """ Find H, affine transformation, such that
+ tp is affine transf of fp. """
+
+ if fp.shape != tp.shape:
+ raise RuntimeError('number of points do not match')
+
+ # condition points
+ # --from points--
+ m = mean(fp[:2], axis=1)
+ maxstd = max(std(fp[:2], axis=1)) + 1e-9
+ C1 = diag([1/maxstd, 1/maxstd, 1])
+ C1[0][2] = -m[0]/maxstd
+ C1[1][2] = -m[1]/maxstd
+ fp_cond = dot(C1,fp)
+
+ # --to points--
+ m = mean(tp[:2], axis=1)
+ C2 = C1.copy() #must use same scaling for both point sets
+ C2[0][2] = -m[0]/maxstd
+ C2[1][2] = -m[1]/maxstd
+ tp_cond = dot(C2,tp)
+
+ # conditioned points have mean zero, so translation is zero
+ A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
+ U,S,V = linalg.svd(A.T)
+
+ # create B and C matrices as Hartley-Zisserman (2:nd ed) p 130.
+ tmp = V[:2].T
+ B = tmp[:2]
+ C = tmp[2:4]
+
+ tmp2 = concatenate((dot(C,linalg.pinv(B)),zeros((2,1))), axis=1)
+ H = vstack((tmp2,[0,0,1]))
+
+ # decondition
+ H = dot(linalg.inv(C2),dot(H,C1))
+
+ return H / H[2,2]
def normalize(points):
- """ Normalize a collection of points in
- homogeneous coordinates so that last row = 1. """
+ """ Normalize a collection of points in
+ homogeneous coordinates so that last row = 1. """
- for row in points:
- row /= points[-1]
- return points
-
+ for row in points:
+ row /= points[-1]
+ return points
+
def make_homog(points):
- """ Convert a set of points (dim*n array) to
- homogeneous coordinates. """
-
- return vstack((points,ones((1,points.shape[1]))))
+ """ Convert a set of points (dim*n array) to
+ homogeneous coordinates. """
+
+ return vstack((points,ones((1,points.shape[1]))))
View
382 PCV/geometry/sfm.py
@@ -4,236 +4,236 @@
def compute_P(x,X):
- """ Compute camera matrix from pairs of
- 2D-3D correspondences (in homog. coordinates). """
-
- n = x.shape[1]
- if X.shape[1] != n:
- raise ValueError("Number of points don't match.")
-
- # create matrix for DLT solution
- M = zeros((3*n,12+n))
- for i in range(n):
- M[3*i,0:4] = X[:,i]
- M[3*i+1,4:8] = X[:,i]
- M[3*i+2,8:12] = X[:,i]
- M[3*i:3*i+3,i+12] = -x[:,i]
-
- U,S,V = linalg.svd(M)
-
- return V[-1,:12].reshape((3,4))
+ """ Compute camera matrix from pairs of
+ 2D-3D correspondences (in homog. coordinates). """
+
+ n = x.shape[1]
+ if X.shape[1] != n:
+ raise ValueError("Number of points don't match.")
+
+ # create matrix for DLT solution
+ M = zeros((3*n,12+n))
+ for i in range(n):
+ M[3*i,0:4] = X[:,i]
+ M[3*i+1,4:8] = X[:,i]
+ M[3*i+2,8:12] = X[:,i]
+ M[3*i:3*i+3,i+12] = -x[:,i]
+
+ U,S,V = linalg.svd(M)
+
+ return V[-1,:12].reshape((3,4))
def triangulate_point(x1,x2,P1,P2):
- """ Point pair triangulation from
- least squares solution. """
-
- M = zeros((6,6))
- M[:3,:4] = P1
- M[3:,:4] = P2
- M[:3,4] = -x1
- M[3:,5] = -x2
+ """ Point pair triangulation from
+ least squares solution. """
+
+ M = zeros((6,6))
+ M[:3,:4] = P1
+ M[3:,:4] = P2
+ M[:3,4] = -x1
+ M[3:,5] = -x2
- U,S,V = linalg.svd(M)
- X = V[-1,:4]
+ U,S,V = linalg.svd(M)
+ X = V[-1,:4]
- return X / X[3]
+ return X / X[3]
def triangulate(x1,x2,P1,P2):
- """ Two-view triangulation of points in
- x1,x2 (3*n homog. coordinates). """
-
- n = x1.shape[1]
- if x2.shape[1] != n:
- raise ValueError("Number of points don't match.")
+ """ Two-view triangulation of points in
+ x1,x2 (3*n homog. coordinates). """
+
+ n = x1.shape[1]
+ if x2.shape[1] != n:
+ raise ValueError("Number of points don't match.")
- X = [ triangulate_point(x1[:,i],x2[:,i],P1,P2) for i in range(n)]
- return array(X).T
+ X = [ triangulate_point(x1[:,i],x2[:,i],P1,P2) for i in range(n)]
+ return array(X).T
def compute_fundamental(x1,x2):
- """ Computes the fundamental matrix from corresponding points
- (x1,x2 3*n arrays) using the 8 point algorithm.
- Each row in the A matrix below is constructed as
- [x'*x, x'*y, x', y'*x, y'*y, y', x, y, 1] """
-
- n = x1.shape[1]
- if x2.shape[1] != n:
- raise ValueError("Number of points don't match.")
-
- # build matrix for equations
- A = zeros((n,9))
- for i in range(n):
- A[i] = [x1[0,i]*x2[0,i], x1[0,i]*x2[1,i], x1[0,i]*x2[2,i],
- x1[1,i]*x2[0,i], x1[1,i]*x2[1,i], x1[1,i]*x2[2,i],
- x1[2,i]*x2[0,i], x1[2,i]*x2[1,i], x1[2,i]*x2[2,i] ]
-
- # compute linear least square solution
- U,S,V = linalg.svd(A)
- F = V[-1].reshape(3,3)
-
- # constrain F
- # make rank 2 by zeroing out last singular value
- U,S,V = linalg.svd(F)
- S[2] = 0
- F = dot(U,dot(diag(S),V))
-
- return F/F[2,2]
+ """ Computes the fundamental matrix from corresponding points
+ (x1,x2 3*n arrays) using the 8 point algorithm.
+ Each row in the A matrix below is constructed as
+ [x'*x, x'*y, x', y'*x, y'*y, y', x, y, 1] """
+
+ n = x1.shape[1]
+ if x2.shape[1] != n:
+ raise ValueError("Number of points don't match.")
+
+ # build matrix for equations
+ A = zeros((n,9))
+ for i in range(n):
+ A[i] = [x1[0,i]*x2[0,i], x1[0,i]*x2[1,i], x1[0,i]*x2[2,i],
+ x1[1,i]*x2[0,i], x1[1,i]*x2[1,i], x1[1,i]*x2[2,i],
+ x1[2,i]*x2[0,i], x1[2,i]*x2[1,i], x1[2,i]*x2[2,i] ]
+
+ # compute linear least square solution
+ U,S,V = linalg.svd(A)
+ F = V[-1].reshape(3,3)
+
+ # constrain F
+ # make rank 2 by zeroing out last singular value
+ U,S,V = linalg.svd(F)
+ S[2] = 0
+ F = dot(U,dot(diag(S),V))
+
+ return F/F[2,2]
def compute_epipole(F):
- """ Computes the (right) epipole from a
- fundamental matrix F.
- (Use with F.T for left epipole.) """
-
- # return null space of F (Fx=0)
- U,S,V = linalg.svd(F)
- e = V[-1]
- return e/e[2]
-
-
+ """ Computes the (right) epipole from a
+ fundamental matrix F.
+ (Use with F.T for left epipole.) """
+
+ # return null space of F (Fx=0)
+ U,S,V = linalg.svd(F)
+ e = V[-1]
+ return e/e[2]
+
+
def plot_epipolar_line(im,F,x,epipole=None,show_epipole=True):
- """ Plot the epipole and epipolar line F*x=0
- in an image. F is the fundamental matrix
- and x a point in the other image."""
-
- m,n = im.shape[:2]
- line = dot(F,x)
-
- # epipolar line parameter and values
- t = linspace(0,n,100)
- lt = array([(line[2]+line[0]*tt)/(-line[1]) for tt in t])
-
- # take only line points inside the image
- ndx = (lt>=0) & (lt<m)
- plot(t[ndx],lt[ndx],linewidth=2)
-
- if show_epipole:
- if epipole is None:
- epipole = compute_epipole(F)
- plot(epipole[0]/epipole[2],epipole[1]/epipole[2],'r*')
-
+ """ Plot the epipole and epipolar line F*x=0
+ in an image. F is the fundamental matrix
+ and x a point in the other image."""
+
+ m,n = im.shape[:2]
+ line = dot(F,x)
+
+ # epipolar line parameter and values
+ t = linspace(0,n,100)
+ lt = array([(line[2]+line[0]*tt)/(-line[1]) for tt in t])
+
+ # take only line points inside the image
+ ndx = (lt>=0) & (lt<m)
+ plot(t[ndx],lt[ndx],linewidth=2)
+
+ if show_epipole:
+ if epipole is None:
+ epipole = compute_epipole(F)
+ plot(epipole[0]/epipole[2],epipole[1]/epipole[2],'r*')
+
def skew(a):
- """ Skew matrix A such that a x v = Av for any v. """
+ """ Skew matrix A such that a x v = Av for any v. """
- return array([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])
+ return array([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])
def compute_P_from_fundamental(F):
- """ Computes the second camera matrix (assuming P1 = [I 0])
- from a fundamental matrix. """
-
- e = compute_epipole(F.T) # left epipole
- Te = skew(e)
- return vstack((dot(Te,F.T).T,e)).T
+ """ Computes the second camera matrix (assuming P1 = [I 0])
+ from a fundamental matrix. """
+
+ e = compute_epipole(F.T) # left epipole
+ Te = skew(e)
+ return vstack((dot(Te,F.T).T,e)).T
def compute_P_from_essential(E):
- """ Computes the second camera matrix (assuming P1 = [I 0])
- from an essential matrix. Output is a list of four
- possible camera matrices. """
-
- # make sure E is rank 2
- U,S,V = svd(E)
- if det(dot(U,V))<0:
- V = -V
- E = dot(U,dot(diag([1,1,0]),V))
-
- # create matrices (Hartley p 258)
- Z = skew([0,0,-1])
- W = array([[0,-1,0],[1,0,0],[0,0,1]])
-
- # return all four solutions
- P2 = [vstack((dot(U,dot(W,V)).T,U[:,2])).T,
- vstack((dot(U,dot(W,V)).T,-U[:,2])).T,
- vstack((dot(U,dot(W.T,V)).T,U[:,2])).T,
- vstack((dot(U,dot(W.T,V)).T,-U[:,2])).T]
-
- return P2
+ """ Computes the second camera matrix (assuming P1 = [I 0])
+ from an essential matrix. Output is a list of four
+ possible camera matrices. """
+
+ # make sure E is rank 2
+ U,S,V = svd(E)
+ if det(dot(U,V))<0:
+ V = -V
+ E = dot(U,dot(diag([1,1,0]),V))
+
+ # create matrices (Hartley p 258)
+ Z = skew([0,0,-1])
+ W = array([[0,-1,0],[1,0,0],[0,0,1]])
+
+ # return all four solutions
+ P2 = [vstack((dot(U,dot(W,V)).T,U[:,2])).T,
+ vstack((dot(U,dot(W,V)).T,-U[:,2])).T,
+ vstack((dot(U,dot(W.T,V)).T,U[:,2])).T,
+ vstack((dot(U,dot(W.T,V)).T,-U[:,2])).T]
+
+ return P2
def compute_fundamental_normalized(x1,x2):
- """ Computes the fundamental matrix from corresponding points
- (x1,x2 3*n arrays) using the normalized 8 point algorithm. """
+ """ Computes the fundamental matrix from corresponding points
+ (x1,x2 3*n arrays) using the normalized 8 point algorithm. """
- n = x1.shape[1]
- if x2.shape[1] != n:
- raise ValueError("Number of points don't match.")
+ n = x1.shape[1]
+ if x2.shape[1] != n:
+ raise ValueError("Number of points don't match.")
- # normalize image coordinates
- x1 = x1 / x1[2]
- mean_1 = mean(x1[:2],axis=1)
- S1 = sqrt(2) / std(x1[:2])
- T1 = array([[S1,0,-S1*mean_1[0]],[0,S1,-S1*mean_1[1]],[0,0,1]])
- x1 = dot(T1,x1)
-
- x2 = x2 / x2[2]
- mean_2 = mean(x2[:2],axis=1)
- S2 = sqrt(2) / std(x2[:2])
- T2 = array([[S2,0,-S2*mean_2[0]],[0,S2,-S2*mean_2[1]],[0,0,1]])
- x2 = dot(T2,x2)
+ # normalize image coordinates
+ x1 = x1 / x1[2]
+ mean_1 = mean(x1[:2],axis=1)
+ S1 = sqrt(2) / std(x1[:2])
+ T1 = array([[S1,0,-S1*mean_1[0]],[0,S1,-S1*mean_1[1]],[0,0,1]])
+ x1 = dot(T1,x1)
+
+ x2 = x2 / x2[2]
+ mean_2 = mean(x2[:2],axis=1)
+ S2 = sqrt(2) / std(x2[:2])
+ T2 = array([[S2,0,-S2*mean_2[0]],[0,S2,-S2*mean_2[1]],[0,0,1]])
+ x2 = dot(T2,x2)
- # compute F with the normalized coordinates
- F = compute_fundamental(x1,x2)
+ # compute F with the normalized coordinates
+ F = compute_fundamental(x1,x2)
- # reverse normalization
- F = dot(T1.T,dot(F,T2))
+ # reverse normalization
+ F = dot(T1.T,dot(F,T2))
- return F/F[2,2]
+ return F/F[2,2]
class RansacModel(object):
- """ Class for fundmental matrix fit with ransac.py from
- http://www.scipy.org/Cookbook/RANSAC"""
-
- def __init__(self,debug=False):
- self.debug = debug
-
- def fit(self,data):
- """ Estimate fundamental matrix using eight
- selected correspondences. """
-
- # transpose and split data into the two point sets
- data = data.T
- x1 = data[:3,:8]
- x2 = data[3:,:8]
-
- # estimate fundamental matrix and return
- F = compute_fundamental_normalized(x1,x2)
- return F
-
- def get_error(self,data,F):
- """ Compute x^T F x for all correspondences,
- return error for each transformed point. """
-
- # transpose and split data into the two point
- data = data.T
- x1 = data[:3]
- x2 = data[3:]
-
- # Sampson distance as error measure
- Fx1 = dot(F,x1)
- Fx2 = dot(F,x2)
- denom = Fx1[0]**2 + Fx1[1]**2 + Fx2[0]**2 + Fx2[1]**2
- err = ( diag(dot(x1.T,dot(F,x2))) )**2 / denom
-
- # return error per point
- return err
+ """ Class for fundmental matrix fit with ransac.py from
+ http://www.scipy.org/Cookbook/RANSAC"""
+
+ def __init__(self,debug=False):
+ self.debug = debug
+
+ def fit(self,data):
+ """ Estimate fundamental matrix using eight
+ selected correspondences. """
+
+ # transpose and split data into the two point sets
+ data = data.T
+ x1 = data[:3,:8]
+ x2 = data[3:,:8]
+
+ # estimate fundamental matrix and return
+ F = compute_fundamental_normalized(x1,x2)
+ return F
+
+ def get_error(self,data,F):
+ """ Compute x^T F x for all correspondences,
+ return error for each transformed point. """
+
+ # transpose and split data into the two point
+ data = data.T
+ x1 = data[:3]
+ x2 = data[3:]
+
+ # Sampson distance as error measure
+ Fx1 = dot(F,x1)
+ Fx2 = dot(F,x2)
+ denom = Fx1[0]**2 + Fx1[1]**2 + Fx2[0]**2 + Fx2[1]**2
+ err = ( diag(dot(x1.T,dot(F,x2))) )**2 / denom
+
+ # return error per point
+ return err
def F_from_ransac(x1,x2,model,maxiter=5000,match_theshold=1e-6):
- """ Robust estimation of a fundamental matrix F from point
- correspondences using RANSAC (ransac.py from
- http://www.scipy.org/Cookbook/RANSAC).
+ """ Robust estimation of a fundamental matrix F from point
+ correspondences using RANSAC (ransac.py from
+ http://www.scipy.org/Cookbook/RANSAC).
- input: x1,x2 (3*n arrays) points in hom. coordinates. """
+ input: x1,x2 (3*n arrays) points in hom. coordinates. """
- import ransac
+ import ransac
- data = vstack((x1,x2))
-
- # compute F and return with inlier index
- F,ransac_data = ransac.ransac(data.T,model,8,maxiter,match_theshold,20,return_all=True)
- return F, ransac_data['inliers']
+ data = vstack((x1,x2))
+
+ # compute F and return with inlier index
+ F,ransac_data = ransac.ransac(data.T,model,8,maxiter,match_theshold,20,return_all=True)
+ return F, ransac_data['inliers']
View
258 PCV/geometry/warp.py
@@ -7,151 +7,151 @@
def image_in_image(im1,im2,tp):
- """ Put im1 in im2 with an affine transformation
- such that corners are as close to tp as possible.
- tp are homogeneous and counter-clockwise from top left. """
-
- # points to warp from
- m,n = im1.shape[:2]
- fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
-
- # compute affine transform and apply
- H = homography.Haffine_from_points(tp,fp)
- im1_t = ndimage.affine_transform(im1,H[:2,:2],
- (H[0,2],H[1,2]),im2.shape[:2])
- alpha = (im1_t > 0)
-
- return (1-alpha)*im2 + alpha*im1_t
+ """ Put im1 in im2 with an affine transformation
+ such that corners are as close to tp as possible.
+ tp are homogeneous and counter-clockwise from top left. """
+
+ # points to warp from
+ m,n = im1.shape[:2]
+ fp = array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
+
+ # compute affine transform and apply
+ H = homography.Haffine_from_points(tp,fp)
+ im1_t = ndimage.affine_transform(im1,H[:2,:2],
+ (H[0,2],H[1,2]),im2.shape[:2])
+ alpha = (im1_t > 0)
+
+ return (1-alpha)*im2 + alpha*im1_t
def combine_images(im1,im2,alpha):
- """ Blend two images with weights as in alpha. """
- return (1-alpha)*im1 + alpha*im2
+ """ Blend two images with weights as in alpha. """
+ return (1-alpha)*im1 + alpha*im2
def alpha_for_triangle(points,m,n):
- """ Creates alpha map of size (m,n)
- for a triangle with corners defined by points
- (given in normalized homogeneous coordinates). """
-
- alpha = zeros((m,n))
- for i in range(min(points[0]),max(points[0])):
- for j in range(min(points[1]),max(points[1])):
- x = linalg.solve(points,[i,j,1])
- if min(x) > 0: #all coefficients positive
- alpha[i,j] = 1
- return alpha
+ """ Creates alpha map of size (m,n)
+ for a triangle with corners defined by points
+ (given in normalized homogeneous coordinates). """
+
+ alpha = zeros((m,n))
+ for i in range(min(points[0]),max(points[0])):
+ for j in range(min(points[1]),max(points[1])):
+ x = linalg.solve(points,[i,j,1])
+ if min(x) > 0: #all coefficients positive
+ alpha[i,j] = 1
+ return alpha
def triangulate_points(x,y):
- """ Delaunay triangulation of 2D points. """
-
- centers,edges,tri,neighbors = md.delaunay(x,y)
- return tri
+ """ Delaunay triangulation of 2D points. """
+
+ centers,edges,tri,neighbors = md.delaunay(x,y)
+ return tri
def plot_mesh(x,y,tri):
- """ Plot triangles. """
-
- for t in tri:
- t_ext = [t[0], t[1], t[2], t[0]] # add first point to end
- plot(x[t_ext],y[t_ext],'r')
+ """ Plot triangles. """
+
+ for t in tri:
+ t_ext = [t[0], t[1], t[2], t[0]] # add first point to end
+ plot(x[t_ext],y[t_ext],'r')
def pw_affine(fromim,toim,fp,tp,tri):
- """ Warp triangular patches from an image.
- fromim = image to warp
- toim = destination image
- fp = from points in hom. coordinates
- tp = to points in hom. coordinates
- tri = triangulation. """
-
- im = toim.copy()
-
- # check if image is grayscale or color
- is_color = len(fromim.shape) == 3
-
- # create image to warp to (needed if iterate colors)
- im_t = zeros(im.shape, 'uint8')
-
- for t in tri:
- # compute affine transformation
- H = homography.Haffine_from_points(tp[:,t],fp[:,t])
-
- if is_color:
- for col in range(fromim.shape[2]):
- im_t[:,:,col] = ndimage.affine_transform(
- fromim[:,:,col],H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
- else:
- im_t = ndimage.affine_transform(
- fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
-
- # alpha for triangle
- alpha = alpha_for_triangle(tp[:,t],im.shape[0],im.shape[1])
-
- # add triangle to image
- im[alpha>0] = im_t[alpha>0]
-
- return im
+ """ Warp triangular patches from an image.
+ fromim = image to warp
+ toim = destination image
+ fp = from points in hom. coordinates
+ tp = to points in hom. coordinates
+ tri = triangulation. """
+
+ im = toim.copy()
+
+ # check if image is grayscale or color
+ is_color = len(fromim.shape) == 3
+
+ # create image to warp to (needed if iterate colors)
+ im_t = zeros(im.shape, 'uint8')
+
+ for t in tri:
+ # compute affine transformation
+ H = homography.Haffine_from_points(tp[:,t],fp[:,t])
+
+ if is_color:
+ for col in range(fromim.shape[2]):
+ im_t[:,:,col] = ndimage.affine_transform(
+ fromim[:,:,col],H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
+ else:
+ im_t = ndimage.affine_transform(
+ fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
+
+ # alpha for triangle
+ alpha = alpha_for_triangle(tp[:,t],im.shape[0],im.shape[1])
+
+ # add triangle to image
+ im[alpha>0] = im_t[alpha>0]
+
+ return im
def panorama(H,fromim,toim,padding=2400,delta=2400):
- """ Create horizontal panorama by blending two images
- using a homography H (preferably estimated using RANSAC).
- The result is an image with the same height as toim. 'padding'
- specifies number of fill pixels and 'delta' additional translation. """
-
- # check if images are grayscale or color
- is_color = len(fromim.shape) == 3
-
- # homography transformation for geometric_transform()
- def transf(p):
- p2 = dot(H,[p[0],p[1],1])
- return (p2[0]/p2[2],p2[1]/p2[2])
-
- if H[1,2]<0: # fromim is to the right
- print 'warp - right'
- # transform fromim
- if is_color:
- # pad the destination image with zeros to the right
- toim_t = hstack((toim,zeros((toim.shape[0],padding,3))))
- fromim_t = zeros((toim.shape[0],toim.shape[1]+padding,toim.shape[2]))
- for col in range(3):
- fromim_t[:,:,col] = ndimage.geometric_transform(fromim[:,:,col],
- transf,(toim.shape[0],toim.shape[1]+padding))
- else:
- # pad the destination image with zeros to the right
- toim_t = hstack((toim,zeros((toim.shape[0],padding))))
- fromim_t = ndimage.geometric_transform(fromim,transf,
- (toim.shape[0],toim.shape[1]+padding))
- else:
- print 'warp - left'
- # add translation to compensate for padding to the left
- H_delta = array([[1,0,0],[0,1,-delta],[0,0,1]])
- H = dot(H,H_delta)
- # transform fromim
- if is_color:
- # pad the destination image with zeros to the left
- toim_t = hstack((zeros((toim.shape[0],padding,3)),toim))
- fromim_t = zeros((toim.shape[0],toim.shape[1]+padding,toim.shape[2]))
- for col in range(3):
- fromim_t[:,:,col] = ndimage.geometric_transform(fromim[:,:,col],
- transf,(toim.shape[0],toim.shape[1]+padding))
- else:
- # pad the destination image with zeros to the left
- toim_t = hstack((zeros((toim.shape[0],padding)),toim))
- fromim_t = ndimage.geometric_transform(fromim,
- transf,(toim.shape[0],toim.shape[1]+padding))
-
- # blend and return (put fromim above toim)
- if is_color:
- # all non black pixels
- alpha = ((fromim_t[:,:,0] * fromim_t[:,:,1] * fromim_t[:,:,2] ) > 0)
- for col in range(3):
- toim_t[:,:,col] = fromim_t[:,:,col]*alpha + toim_t[:,:,col]*(1-alpha)
- else:
- alpha = (fromim_t > 0)
- toim_t = fromim_t*alpha + toim_t*(1-alpha)
-
- return toim_t
+ """ Create horizontal panorama by blending two images
+ using a homography H (preferably estimated using RANSAC).
+ The result is an image with the same height as toim. 'padding'
+ specifies number of fill pixels and 'delta' additional translation. """
+
+ # check if images are grayscale or color
+ is_color = len(fromim.shape) == 3
+
+ # homography transformation for geometric_transform()
+ def transf(p):
+ p2 = dot(H,[p[0],p[1],1])
+ return (p2[0]/p2[2],p2[1]/p2[2])
+
+ if H[1,2]<0: # fromim is to the right
+ print 'warp - right'
+ # transform fromim
+ if is_color:
+ # pad the destination image with zeros to the right
+ toim_t = hstack((toim,zeros((toim.shape[0],padding,3))))
+ fromim_t = zeros((toim.shape[0],toim.shape[1]+padding,toim.shape[2]))
+ for col in range(3):
+ fromim_t[:,:,col] = ndimage.geometric_transform(fromim[:,:,col],
+ transf,(toim.shape[0],toim.shape[1]+padding))
+ else:
+ # pad the destination image with zeros to the right
+ toim_t = hstack((toim,zeros((toim.shape[0],padding))))
+ fromim_t = ndimage.geometric_transform(fromim,transf,
+ (toim.shape[0],toim.shape[1]+padding))
+ else:
+ print 'warp - left'
+ # add translation to compensate for padding to the left
+ H_delta = array([[1,0,0],[0,1,-delta],[0,0,1]])
+ H = dot(H,H_delta)
+ # transform fromim
+ if is_color:
+ # pad the destination image with zeros to the left
+ toim_t = hstack((zeros((toim.shape[0],padding,3)),toim))
+ fromim_t = zeros((toim.shape[0],toim.shape[1]+padding,toim.shape[2]))
+ for col in range(3):
+ fromim_t[:,:,col] = ndimage.geometric_transform(fromim[:,:,col],
+ transf,(toim.shape[0],toim.shape[1]+padding))
+ else:
+ # pad the destination image with zeros to the left
+ toim_t = hstack((zeros((toim.shape[0],padding)),toim))
+ fromim_t = ndimage.geometric_transform(fromim,
+ transf,(toim.shape[0],toim.shape[1]+padding))
+
+ # blend and return (put fromim above toim)
+ if is_color:
+ # all non black pixels
+ alpha = ((fromim_t[:,:,0] * fromim_t[:,:,1] * fromim_t[:,:,2] ) > 0)
+ for col in range(3):
+ toim_t[:,:,col] = fromim_t[:,:,col]*alpha + toim_t[:,:,col]*(1-alpha)
+ else:
+ alpha = (fromim_t > 0)
+ toim_t = fromim_t*alpha + toim_t*(1-alpha)
+
+ return toim_t
View
BIN  PCV/imagesearch/__init__.pyc
Binary file not shown
View
334 PCV/imagesearch/imagesearch.py
@@ -4,172 +4,172 @@
class Indexer(object):
-
- def __init__(self,db,voc):
- """ Initialize with the name of the database
- and a vocabulary object. """
-
- self.con = sqlite.connect(db)
- self.voc = voc
-
- def __del__(self):
- self.con.close()
-
- def db_commit(self):
- self.con.commit()
-
- def get_id(self,imname):
- """ Get an entry id and add if not present. """
-
- cur = self.con.execute(
- "select rowid from imlist where filename='%s'" % imname)
- res=cur.fetchone()
- if res==None:
- cur = self.con.execute(
- "insert into imlist(filename) values ('%s')" % imname)
- return cur.lastrowid
- else:
- return res[0]
-
- def is_indexed(self,imname):
- """ Returns True if imname has been indexed. """
-
- im = self.con.execute("select rowid from imlist where filename='%s'" % imname).fetchone()
- return im != None
-
- def add_to_index(self,imname,descr):
- """ Take an image with feature descriptors,
- project on vocabulary and add to database. """
-
- if self.is_indexed(imname): return
- print 'indexing', imname
-
- # get the imid
- imid = self.get_id(imname)
-
- # get the words
- imwords = self.voc.project(descr)
- nbr_words = imwords.shape[0]
-
- # link each word to image
- for i in range(nbr_words):
- word = imwords[i]
- # wordid is the word number itself
- self.con.execute("insert into imwords(imid,wordid,vocname) values (?,?,?)", (imid,word,self.voc.name))
-
- # store word histogram for image
- # use pickle to encode NumPy arrays as strings
- self.con.execute("insert into imhistograms(imid,histogram,vocname) values (?,?,?)", (imid,pickle.dumps(imwords),self.voc.name))
-
- def create_tables(self):
- """ Create the database tables. """
-
- self.con.execute('create table imlist(filename)')
- self.con.execute('create table imwords(imid,wordid,vocname)')
- self.con.execute('create table imhistograms(imid,histogram,vocname)')
- self.con.execute('create index im_idx on imlist(filename)')
- self.con.execute('create index wordid_idx on imwords(wordid)')
- self.con.execute('create index imid_idx on imwords(imid)')
- self.con.execute('create index imidhist_idx on imhistograms(imid)')
- self.db_commit()
+
+ def __init__(self,db,voc):
+ """ Initialize with the name of the database
+ and a vocabulary object. """
+
+ self.con = sqlite.connect(db)
+ self.voc = voc
+
+ def __del__(self):
+ self.con.close()
+
+ def db_commit(self):
+ self.con.commit()
+
+ def get_id(self,imname):
+ """ Get an entry id and add if not present. """
+
+ cur = self.con.execute(
+ "select rowid from imlist where filename='%s'" % imname)
+ res=cur.fetchone()
+ if res==None:
+ cur = self.con.execute(
+ "insert into imlist(filename) values ('%s')" % imname)
+ return cur.lastrowid
+ else:
+ return res[0]
+
+ def is_indexed(self,imname):
+ """ Returns True if imname has been indexed. """
+
+ im = self.con.execute("select rowid from imlist where filename='%s'" % imname).fetchone()
+ return im != None
+
+ def add_to_index(self,imname,descr):
+ """ Take an image with feature descriptors,
+ project on vocabulary and add to database. """
+
+ if self.is_indexed(imname): return
+ print 'indexing', imname
+
+ # get the imid
+ imid = self.get_id(imname)
+
+ # get the words
+ imwords = self.voc.project(descr)
+ nbr_words = imwords.shape[0]
+
+ # link each word to image
+ for i in range(nbr_words):
+ word = imwords[i]
+ # wordid is the word number itself
+ self.con.execute("insert into imwords(imid,wordid,vocname) values (?,?,?)", (imid,word,self.voc.name))
+
+ # store word histogram for image
+ # use pickle to encode NumPy arrays as strings
+ self.con.execute("insert into imhistograms(imid,histogram,vocname) values (?,?,?)", (imid,pickle.dumps(imwords),self.voc.name))
+
+ def create_tables(self):
+ """ Create the database tables. """
+
+ self.con.execute('create table imlist(filename)')
+ self.con.execute('create table imwords(imid,wordid,vocname)')
+ self.con.execute('create table imhistograms(imid,histogram,vocname)')
+ self.con.execute('create index im_idx on imlist(filename)')
+ self.con.execute('create index wordid_idx on imwords(wordid)')
+ self.con.execute('create index imid_idx on imwords(imid)')
+ self.con.execute('create index imidhist_idx on imhistograms(imid)')
+ self.db_commit()
class Searcher(object):
-
- def __init__(self,db,voc):
- """ Initialize with the name of the database. """
- self.con = sqlite.connect(db)
- self.voc = voc
-
- def __del__(self):
- self.con.close()
-
- def get_imhistogram(self,imname):
- """ Return the word histogram for an image. """
-
- im_id = self.con.execute(
- "select rowid from imlist where filename='%s'" % imname).fetchone()
- s = self.con.execute(
- "select histogram from imhistograms where rowid='%d'" % im_id).fetchone()
-
- # use pickle to decode NumPy arrays from string
- return pickle.loads(str(s[0]))
-
- def candidates_from_word(self,imword):
- """ Get list of images containing imword. """
-
- im_ids = self.con.execute(
- "select distinct imid from imwords where wordid=%d" % imword).fetchall()
- return [i[0] for i in im_ids]
-
- def candidates_from_histogram(self,imwords):
- """ Get list of images with similar words. """
-
- # get the word ids
- words = imwords.nonzero()[0]
-
- # find candidates
- candidates = []
- for word in words:
- c = self.candidates_from_word(word)
- candidates+=c
-
- # take all unique words and reverse sort on occurrence
- tmp = [(w,candidates.count(w)) for w in set(candidates)]
- tmp.sort(cmp=lambda x,y:cmp(x[1],y[1]))
- tmp.reverse()
-
- # return sorted list, best matches first
- return [w[0] for w in tmp]
-
- def query(self,imname):
- """ Find a list of matching images for imname. """
-
- h = self.get_imhistogram(imname)
- candidates = self.candidates_from_histogram(h)
-
- matchscores = []
- for imid in candidates:
- # get the name
- cand_name = self.con.execute(
- "select filename from imlist where rowid=%d" % imid).fetchone()
- cand_h = self.get_imhistogram(cand_name)
- cand_dist = sqrt( sum( self.voc.idf*(h-cand_h)**2 ) )
- matchscores.append( (cand_dist,imid) )
-
- # return a sorted list of distances and database ids
- matchscores.sort()
- return matchscores
-
- def get_filename(self,imid):
- """ Return the filename for an image id. """
-
- s = self.con.execute(
- "select filename from imlist where rowid='%d'" % imid).fetchone()
- return s[0]
+
+ def __init__(self,db,voc):
+ """ Initialize with the name of the database. """
+ self.con = sqlite.connect(db)
+ self.voc = voc
+
+ def __del__(self):
+ self.con.close()
+
+ def get_imhistogram(self,imname):
+ """ Return the word histogram for an image. """
+
+ im_id = self.con.execute(
+ "select rowid from imlist where filename='%s'" % imname).fetchone()
+ s = self.con.execute(
+ "select histogram from imhistograms where rowid='%d'" % im_id).fetchone()
+
+ # use pickle to decode NumPy arrays from string
+ return pickle.loads(str(s[0]))
+
+ def candidates_from_word(self,imword):
+ """ Get list of images containing imword. """
+
+ im_ids = self.con.execute(
+ "select distinct imid from imwords where wordid=%d" % imword).fetchall()
+ return [i[0] for i in im_ids]
+
+ def candidates_from_histogram(self,imwords):
+ """ Get list of images with similar words. """
+
+ # get the word ids
+ words = imwords.nonzero()[0]
+
+ # find candidates
+ candidates = []
+ for word in words:
+ c = self.candidates_from_word(word)
+ candidates+=c
+
+ # take all unique words and reverse sort on occurrence
+ tmp = [(w,candidates.count(w)) for w in set(candidates)]
+ tmp.sort(cmp=lambda x,y:cmp(x[1],y[1]))
+ tmp.reverse()
+
+ # return sorted list, best matches first
+ return [w[0] for w in tmp]
+
+ def query(self,imname):
+ """ Find a list of matching images for imname. """
+
+ h = self.get_imhistogram(imname)
+ candidates = self.candidates_from_histogram(h)
+
+ matchscores = []
+ for imid in candidates:
+ # get the name
+ cand_name = self.con.execute(
+ "select filename from imlist where rowid=%d" % imid).fetchone()
+ cand_h = self.get_imhistogram(cand_name)
+ cand_dist = sqrt( sum( self.voc.idf*(h-cand_h)**2 ) )
+ matchscores.append( (cand_dist,imid) )
+
+ # return a sorted list of distances and database ids
+ matchscores.sort()
+ return matchscores
+
+ def get_filename(self,imid):
+ """ Return the filename for an image id. """
+
+ s = self.con.execute(
+ "select filename from imlist where rowid='%d'" % imid).fetchone()
+ return s[0]
def tf_idf_dist(voc,v1,v2):
-
- v1 /= sum(v1)
- v2 /= sum(v2)
-
- return sqrt( sum( voc.idf*(v1-v2)**2 ) )
+
+ v1 /= sum(v1)
+ v2 /= sum(v2)
+
+ return sqrt( sum( voc.idf*(v1-v2)**2 ) )
def compute_ukbench_score(src,imlist):
- """ Returns the average number of correct
- images on the top four results of queries. """
-
- nbr_images = len(imlist)
- pos = zeros((nbr_images,4))
- # get first four results for each image
- for i in range(nbr_images):
- pos[i] = [w[1]-1 for w in src.query(imlist[i])[:4]]
-
- # compute score and return average
- score = array([ (pos[i]//4)==(i//4) for i in range(nbr_images)])*1.0
- return sum(score) / (nbr_images)
+ """ Returns the average number of correct
+ images on the top four results of queries. """
+
+ nbr_images = len(imlist)
+ pos = zeros((nbr_images,4))
+ # get first four results for each image
+ for i in range(nbr_images):
+ pos[i] = [w[1]-1 for w in src.query(imlist[i])[:4]]
+
+ # compute score and return average
+ score = array([ (pos[i]//4)==(i//4) for i in range(nbr_images)])*1.0
+ return sum(score) / (nbr_images)
# import PIL and pylab for plotting
@@ -177,13 +177,13 @@ def compute_ukbench_score(src,imlist):
from pylab import *
def plot_results(src,res):
- """ Show images in result list 'res'. """
-
- figure()
- nbr_results = len(res)
- for i in range(nbr_results):
- imname = src.get_filename(res[i])
- subplot(1,nbr_results,i+1)
- imshow(array(Image.open(imname)))
- axis('off')
- show()
+
+ figure()
+ nbr_results = len(res)
+ for i in range(nbr_results):
+ imname = src.get_filename(res[i])
+ subplot(1,nbr_results,i+1)
+ imshow(array(Image.open(imname)))
+ axis('off')
+ show()
View
BIN  PCV/imagesearch/imagesearch.pyc
Binary file not shown
View
102 PCV/imagesearch/vocabulary.py
@@ -4,54 +4,54 @@
class Vocabulary(object):
-
- def __init__(self,name):
- self.name = name
- self.voc = []
- self.idf = []
- self.trainingdata = []
- self.nbr_words = 0
-
- def train(self,featurefiles,k=100,subsampling=10):
- """ Train a vocabulary from features in files listed
- in featurefiles using k-means with k number of words.
- Subsampling of training data can be used for speedup. """
-
- nbr_images = len(featurefiles)
- # read the features from file
- descr = []
- descr.append(sift.read_features_from_file(featurefiles[0])[1])
- descriptors = descr[0] #stack all features for k-means
- for i in arange(1,nbr_images):
- descr.append(sift.read_features_from_file(featurefiles[i])[1])
- descriptors = vstack((descriptors,descr[i]))
-
- # k-means: last number determines number of runs
- self.voc,distortion = kmeans(descriptors[::subsampling,:],k,1)
- self.nbr_words = self.voc.shape[0]
-
- # go through all training images and project on vocabulary
- imwords = zeros((nbr_images,self.nbr_words))
- for i in range( nbr_images ):
- imwords[i] = self.project(descr[i])
-
- nbr_occurences = sum( (imwords > 0)*1 ,axis=0)
-
- self.idf = log( (1.0*nbr_images) / (1.0*nbr_occurences+1) )
- self.trainingdata = featurefiles
-
- def project(self,descriptors):
- """ Project descriptors on the vocabulary
- to create a histogram of words. """
-
- # histogram of image words
- imhist = zeros((self.nbr_words))
- words,distance = vq(descriptors,self.voc)
- for w in words:
- imhist[w] += 1
-
- return imhist
-
- def get_words(self,descriptors):
- """ Convert descriptors to words. """
- return vq(descriptors,self.voc)[0]
+
+ def __init__(self,name):
+ self.name = name
+ self.voc = []
+ self.idf = []
+ self.trainingdata = []
+ self.nbr_words = 0
+
+ def train(self,featurefiles,k=100,subsampling=10):
+ """ Train a vocabulary from features in files listed
+ in featurefiles using k-means with k number of words.
+ Subsampling of training data can be used for speedup. """
+
+ nbr_images = len(featurefiles)
+ # read the features from file
+ descr = []
+ descr.append(sift.read_features_from_file(featurefiles[0])[1])
+ descriptors = descr[0] #stack all features for k-means
+ for i in arange(1,nbr_images):
+ descr.append(sift.read_features_from_file(featurefiles[i])[1])
+ descriptors = vstack((descriptors,descr[i]))
+
+ # k-means: last number determines number of runs
+ self.voc,distortion = kmeans(descriptors[::subsampling,:],k,1)
+ self.nbr_words = self.voc.shape[0]
+
+ # go through all training images and project on vocabulary
+ imwords = zeros((nbr_images,self.nbr_words))
+ for i in range( nbr_images ):
+ imwords[i] = self.project(descr[i])
+
+ nbr_occurences = sum( (imwords > 0)*1 ,axis=0)
+
+ self.idf = log( (1.0*nbr_images) / (1.0*nbr_occurences+1) )
+ self.trainingdata = featurefiles
+
+ def project(self,descriptors):
+ """ Project descriptors on the vocabulary
+ to create a histogram of words. """
+
+ # histogram of image words
+ imhist = zeros((self.nbr_words))
+ words,distance = vq(descriptors,self.voc)
+ for w in words:
+ imhist[w] += 1
+
+ return imhist
+
+ def get_words(self,descriptors):
+ """ Convert descriptors to words. """
+ return vq(descriptors,self.voc)[0]
View
BIN  PCV/localdescriptors/__init__.pyc
Binary file not shown
View
56 PCV/localdescriptors/dsift.py
@@ -5,34 +5,34 @@
def process_image_dsift(imagename,resultname,size=20,steps=10,force_orientation=False,resize=None):
- """ Process an image with densely sampled SIFT descriptors
- and save the results in a file. Optional input: size of features,
- steps between locations, forcing computation of descriptor orientation
- (False means all are oriented upwards), tuple for resizing the image."""
+ """ Process an image with densely sampled SIFT descriptors
+ and save the results in a file. Optional input: size of features,
+ steps between locations, forcing computation of descriptor orientation
+ (False means all are oriented upwards), tuple for resizing the image."""
- im = Image.open(imagename).convert('L')
- if resize!=None:
- im = im.resize(resize)
- m,n = im.size
-
- if imagename[-3:] != 'pgm':
- #create a pgm file
- im.save('tmp.pgm')
- imagename = 'tmp.pgm'
+ im = Image.open(imagename).convert('L')
+ if resize!=None:
+ im = im.resize(resize)
+ m,n = im.size
+
+ if imagename[-3:] != 'pgm':
+ #create a pgm file
+ im.save('tmp.pgm')
+ imagename = 'tmp.pgm'
- # create frames and save to temporary file
- scale = size/3.0
- x,y = meshgrid(range(steps,m,steps),range(steps,n,steps))
- xx,yy = x.flatten(),y.flatten()
- frame = array([xx,yy,scale*ones(xx.shape[0]),zeros(xx.shape[0])])
- savetxt('tmp.frame',frame.T,fmt='%03.3f')
-
- if force_orientation:
- cmmd = str("sift "+imagename+" --output="+resultname+
- " --read-frames=tmp.frame --orientations")
- else:
- cmmd = str("sift "+imagename+" --output="+resultname+
- " --read-frames=tmp.frame")
- os.system(cmmd)
- print 'processed', imagename, 'to', resultname
+ # create frames and save to temporary file
+ scale = size/3.0
+ x,y = meshgrid(range(steps,m,steps),range(steps,n,steps))
+ xx,yy = x.flatten(),y.flatten()
+ frame = array([xx,yy,scale*ones(xx.shape[0]),zeros(xx.shape[0])])
+ savetxt('tmp.frame',frame.T,fmt='%03.3f')
+
+ if force_orientation:
+ cmmd = str("sift "+imagename+" --output="+resultname+
+ " --read-frames=tmp.frame --orientations")
+ else:
+ cmmd = str("sift "+imagename+" --output="+resultname+
+ " --read-frames=tmp.frame")
+ os.system(cmmd)
+ print 'processed', imagename, 'to', resultname
View
262 PCV/localdescriptors/harris.py
@@ -4,155 +4,155 @@
def compute_harris_response(im,sigma=3):
- """ Compute the Harris corner detector response function
- for each pixel in a graylevel image. """
-
- # derivatives
- imx = zeros(im.shape)
- filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
- imy = zeros(im.shape)
- filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
-
- # compute components of the Harris matrix
- Wxx = filters.gaussian_filter(imx*imx,sigma)
- Wxy = filters.gaussian_filter(imx*imy,sigma)
- Wyy = filters.gaussian_filter(imy*imy,sigma)
-
- # determinant and trace
- Wdet = Wxx*Wyy - Wxy**2
- Wtr = Wxx + Wyy
-
- return Wdet / Wtr
+ """ Compute the Harris corner detector response function
+ for each pixel in a graylevel image. """
+
+ # derivatives