...
 
Commits (2)
......@@ -39,7 +39,9 @@ int main(int argc, char** argv)
printf("<g class=\"houghlines\">\n");
printf(" <!-- %d total lines -->\n", lines.size());
fprintf(stderr, "%d total lines\n", lines.size());
// fprintf(stderr, "%d total lines\n", lines.size());
// CSV headers
fprintf(stderr, "type,angle,x1,y1,x2,y2\n");
for( size_t i = 0; i < lines.size(); i++ )
{
double Angle = atan2(lines[i][3]- lines[i][1], lines[i][2]- lines[i][0]) * 180.0 / CV_PI;
......@@ -47,13 +49,13 @@ int main(int argc, char** argv)
double dv = abs(abs(Angle) - 90);
if (dh < 5) {
printf(" <line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" class=\"horizontal\" />\n", lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
fprintf(stderr, "H angle:%0.1f line:%d,%d to %d,%d\n", Angle, lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
fprintf(stderr, "horizontal,%0.1f,%d,%d,%d,%d\n", Angle, lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
} else if (dv < 5) {
printf(" <line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" class=\"vertical\" />\n", lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
fprintf(stderr, "V angle:%0.1f line:%d,%d to %d,%d\n", Angle, lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
fprintf(stderr, "vertical,%0.1f,%d,%d,%d,%d\n", Angle, lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
} else {
printf(" <line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" />\n", lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
fprintf(stderr, "l angle:%0.1f line:%d,%d to %d,%d\n", Angle, lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
fprintf(stderr, "line,%0.1f,%d,%d,%d,%d\n", Angle, lines[i][0], lines[i][1], lines[i][2], lines[i][3]);
}
// if(abs(Angle)>=0 && abs(Angle)<=5) {
// could also consider:
......
......@@ -72,7 +72,8 @@ int main(int argc, char** argv)
fprintf(stderr, "SYM:'%s', conf: %.2f; BoundingBox: %d,%d,%d,%d; psize: %d, font id: %d\n", word, conf, x1, y1, x2, y2, psize, fid);
if(conf>50) {
// printf(" <rect x=\"%d\" y=\"%d\" width=\"%d\" height=\"%d\"/>\n", x1, y1, (x2-x1), (y2-y1));
fprintf(f, " <g class=\"symbol\" transform=\"translate(%d,%d)\"><rect x=\"0\" y=\"0\" width=\"%d\" height=\"%d\"/><text x=\"0\" y=\"%d\" style=\"font-size: %dpx\"><![CDATA[%s]]></text></g>\n", x1, y1, (x2-x1), (y2-y1), (y2-y1), psize/2, word);
// fprintf(f, " <g class=\"symbol\" transform=\"translate(%d,%d)\"><rect x=\"0\" y=\"0\" width=\"%d\" height=\"%d\"/><text x=\"0\" y=\"%d\" style=\"font-size: %dpx\"><![CDATA[%s]]></text></g>\n", x1, y1, (x2-x1), (y2-y1), (y2-y1), psize/2, word);
fprintf(f, " <g class=\"symbol\" transform=\"translate(%d,%d)\"><rect x=\"0\" y=\"0\" width=\"%d\" height=\"%d\"/><text x=\"0\" y=\"%d\" style=\"font-size: %dpx\"><![CDATA[%s]]></text></g>\n", x1, y1, (x2-x1), (y2-y1), (y2-y1), 12, word);
// printf(" <g class=\"word\" transform=\"translate(%d,%d)\"><rect x=\"0\" y=\"0\" width=\"%d\" height=\"%d\"/><text x=\"0\" y=\"0\"><![CDATA[%s]]></text></g>\n", x1, y1, (x2-x1), (y2-y1), word);
confidence_level+=conf;
j++;
......@@ -110,7 +111,8 @@ int main(int argc, char** argv)
fprintf(stderr, "WORD:'%s', conf: %.2f; BoundingBox: %d,%d,%d,%d; psize: %d, font id: %d\n", word, conf, x1, y1, x2, y2, psize, fid);
if(conf>50) {
// printf(" <rect x=\"%d\" y=\"%d\" width=\"%d\" height=\"%d\"/>\n", x1, y1, (x2-x1), (y2-y1));
fprintf(f, " <g class=\"word\" transform=\"translate(%d,%d)\"><rect x=\"0\" y=\"0\" width=\"%d\" height=\"%d\"/><text x=\"0\" y=\"%d\" style=\"font-size: %dpx\"><![CDATA[%s]]></text></g>\n", x1, y1, (x2-x1), (y2-y1), (y2-y1), psize/2, word);
//fprintf(f, " <g class=\"word\" transform=\"translate(%d,%d)\"><rect x=\"0\" y=\"0\" width=\"%d\" height=\"%d\"/><text x=\"0\" y=\"%d\" style=\"font-size: %dpx\"><![CDATA[%s]]></text></g>\n", x1, y1, (x2-x1), (y2-y1), (y2-y1), psize/2, word);
fprintf(f, " <g class=\"word\" transform=\"translate(%d,%d)\"><rect x=\"0\" y=\"0\" width=\"%d\" height=\"%d\"/><text x=\"0\" y=\"%d\" style=\"font-size: %dpx\"><![CDATA[%s]]></text></g>\n", x1, y1, (x2-x1), (y2-y1), (y2-y1), 12, word);
// printf(" <g class=\"word\" transform=\"translate(%d,%d)\"><rect x=\"0\" y=\"0\" width=\"%d\" height=\"%d\"/><text x=\"0\" y=\"0\"><![CDATA[%s]]></text></g>\n", x1, y1, (x2-x1), (y2-y1), word);
confidence_level+=conf;
j++;
......
#!/usr/bin/env python
from sh import git
import json
import cgi
fs = cgi.FieldStorage()
cmd = fs.getvalue("cmd", "log")
_id = fs.getvalue("id")
if cmd == "log":
log = git("-c", "core.pager=cat", "log", "--oneline", "--no-color")
print "Content-type: application/json"
print
print json.dumps([x.split(" ", 1) for x in reversed(log.splitlines())])
elif cmd == "ls":
log = git("ls-tree", _id)
print "Content-type: application/json"
print
print json.dumps([x.split(None, 3) for x in log.splitlines()])
elif cmd == "cat":
data = git("cat-file", "blob", _id)
print "Content-type: text/plain"
print
sys.stdout.write(data)
# interrogation commands
# git-cat-file
# git-diff-files,index,tre
# git-ls-files,tree (ls-files: index, working tree)
# git-status --porcelain
\ No newline at end of file
from __future__ import print_function
from argparse import ArgumentParser
import sys, os
from srt import srtparse, timecode
from time import sleep
p = ArgumentParser("Draw SVG as a series of stills determined by an SRT format timeline")
p.add_argument("--srt", default="matthew.srt", help="srt file")
p.add_argument("--originals", default="originals/%s", help="pattern to compose the full path of input images referenced in the srt")
p.add_argument("--svg", default="frames_svg/frame%04d.png", help="")
p.add_argument("--warp", default="frames_warp/frame%04d.png", help="")
p.add_argument("--output", default="frames_composed/frame%04d.png", help="")
p.add_argument("--fps", type=float, default=25.0, help="fps, default: 25.0")
p.add_argument("--gap", type=float, default=0.0, help="time to put black between frames, default: 0.0")
p.add_argument("--duration", type=float, default=None, help="total duration, only used with --gap for last title")
p.add_argument("--force", default=False, action="store_true", help="force even if already generated")
p.add_argument("--verbose", default=False, action="store_true")
args = p.parse_args()
with open(args.srt) as f:
titles = srtparse(f.read())
t = 0
curtitle = -1
framenum = 1
fixedpath = None
while True:
t = ((framenum-1) / args.fps)
if t >= args.duration:
break
while ((curtitle+1) < len(titles)) and (titles[curtitle+1]['start'] <= t):
curtitle += 1
fixedpath = args.originals % titles[curtitle]['content'].strip()
outframepath = args.output % framenum
svgpath = args.svg % framenum
warppath = args.warp % framenum
start = titles[curtitle]['start']
if 'end' in titles[curtitle]:
end = titles[curtitle]['end']
else:
end = args.duration
black = (args.gap > 0.0) and (t + args.gap >= end)
# STATUS
use_fixed = fixedpath
if black:
use_fixed = "(Black frame)"
print ("{0}: {1} + {2} + {3} => {4}".format(timecode(t), use_fixed, svgpath, warppath, outframepath), file=sys.stderr)
framenum += 1
# CREATE THE FRAME
if not args.force and os.path.exists(outframepath) and os.path.getsize(outframepath) > 0:
print ("Skipping already generated frame {0}".format(outframepath), file=sys.stderr)
continue
if os.path.exists(fixedpath) and os.path.exists(svgpath) and os.path.exists(warppath):
use_fixed = fixedpath
if black:
use_fixed = "-size 1080x1080 canvas:black"
cmd = """convert {0} {1} -compose Hard_Light -composite {2} -compose Over -composite {3}"""
cmd = cmd.format(use_fixed, svgpath, warppath, outframepath)
os.system(cmd)
else:
print ("Missing input file, skipping", file=sys.stderr)
sleep(0.01)
\ No newline at end of file
mkdir -p output_composed
for i in ls output/*.png
do
# echo $i
b=`basename $i`
b=${b%.*}
convert $i output_contours/$b.png -compose Hard_Light -composite -compose Over -background black -gravity center -extent 1920x1080 output_composed/$b.png
echo $b
done
# convert $i output_contours/$b.png -compose Hard_Light -composite -compose Over -background black -gravity center -extent 1920x1080 output_composed/$b.png
convert originals/S-20-24-1_1080.jpg frames_svg/frame1000.png -compose Hard_Light -composite frames_warp/frame1000.png -compose Over -composite frames_composed/frame1000.png
convert -size 1080x1080 canvas:black frames_svg/frame1000.png -compose Hard_Light -composite frames_warp/frame1000.png -compose Over -composite frames_composed/frame1000b.png
from __future__ import print_function
from xml.etree import cElementTree as ET
import sys, math
def iterparent(tree):
for parent in tree.getiterator():
for child in parent:
yield parent, child
class ContourSVG (object):
"""
Simple class to render simple SVG's composed of only polyline elements with points.
Wnen saving the file, a 'percentage' can be specified to render only a given percentage (from 0.0-1.0)
of the total number of points
"""
def __init__ (self, path):
self.path = path
self.numpoints = None
def open (self):
# http://stackoverflow.com/questions/8983041/saving-xml-files-using-elementtree
ET.register_namespace("","http://www.w3.org/2000/svg")
self.numpoints = 0
with open(self.path) as f:
tree = ET.parse(f)
root = tree.getroot()
for p in root.findall("{http://www.w3.org/2000/svg}polyline"):
pts = p.attrib.get("points").strip().split()
# print len(pts), pts[0], "...", pts[-1]
self.numpoints += len(pts)
# print self.numpoints, "points"
return tree
def output (self, p=1.0, np=None):
tree = self.open()
root = tree.getroot()
if np == None:
np = int(self.numpoints * p)
cp = 0
rem = []
for parent, child in iterparent(root):
if child.tag == "{http://www.w3.org/2000/svg}polyline":
# remove element if we're done
if cp >= np:
# print ("removing element {0} from {1}".format(child, parent), file=sys.stderr)
# parent.remove(child)
rem.append((parent, child))
else:
pts = child.attrib["points"].strip().split()
if cp + len(pts) > np:
# truncate points
pts = pts[:np-cp]
child.attrib['points'] = " ".join(pts)
cp = np
else:
# let pass and count the points
cp += len(pts)
for p, c in rem:
p.remove(c)
return tree
def save (self, name, p=1.0):
t = self.output(p)
with open(name, "wt") as f:
t.write(f)
if __name__ == "__main__":
svg = ContourSVG(sys.argv[1])
# t = svg.output(0, np=2)
# t.write(sys.stdout)
svg.save("test.svg", 1.0)
# print (ET.tostring(t.getroot()))
# print svg.tree
import cv2, sys, os
import numpy as np
width, height = 1080, 1080
black = np.zeros((height, width, 4), 'uint8')
cv2.imwrite("blank.png", black)
sys.exit()
def ensure_alpha (img):
if img.shape[2] < 4:
alpha = np.ones((img.shape[0], img.shape[1], 4), 'uint8') * 255
alpha[:,:,:-1] = img
return alpha
return img
FLANN_INDEX_KDTREE = 0
index_params = dict(algorithm = FLANN_INDEX_KDTREE, trees=5)
search_params = dict(checks = 50)
sift = cv2.xfeatures2d.SIFT_create()
flann = cv2.FlannBasedMatcher(index_params, search_params)
place_image = cv2.imread("frames/frame0001.jpg",0)
place_image_color = ensure_alpha(cv2.imread("frames/frame0001.jpg",-1))
# place_image_color = cv2.imread("frames/frame0001.png",-1)
print place_image_color.shape
fixed_image = cv2.imread("originals/S-20-28-38_1080.jpg",0) # training
kp1, des1 = sift.detectAndCompute(place_image,None)
kp2, des2 = sift.detectAndCompute(fixed_image, None)
matches = flann.knnMatch(des1,des2,k=2) # order: query, training
good = []
for m, n in matches:
if m.distance < 0.7*n.distance:
good.append(m)
if len(good) < 10:
print "not enough matches"
sys.exit()
src_pts = np.float32([ kp1[m.queryIdx].pt for m in good ]).reshape(-1, 1, 2)
dst_pts = np.float32([ kp2[m.trainIdx].pt for m in good ]).reshape(-1, 1, 2)
M, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)
if M == None:
print "findHomography failed"
sys.exit()
if M != None:
# http://docs.opencv.org/modules/imgproc/doc/geometric_transformations.html#warpperspective
warped_image = cv2.warpPerspective(place_image_color, M, fixed_image.shape)
cv2.imwrite("output.png", warped_image)
from PCV.classifiers import *
from PCV.clustering import *
from PCV.geometry import *
from PCV.imagesearch import *
from PCV.localdescriptors import *
from PCV.tools import *
\ No newline at end of file
from numpy import *
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. """
# 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)
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 L2dist(p1,p2):
return sqrt( sum( (p1-p2)**2) )
def L1dist(v1,v2):
return sum(abs(v1-v2))
\ No newline at end of file
from numpy import *
from itertools import combinations
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)
class ClusterLeafNode(object):
def __init__(self,vec,id):
self.vec = vec
self.id = id
def extract_clusters(self,dist):
return [self]
def get_cluster_elements(self):
return [self.id]
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 L2dist(v1,v2):
return sqrt(sum((v1-v2)**2))
def L1dist(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]
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()
from numpy import *
from scipy import linalg
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. """
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
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
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 numpy import *
from scipy import ndimage
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) )
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. """
from PCV.tools 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]
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]
def normalize(points):
""" Normalize a collection of points in
homogeneous coordinates so that last row = 1. """
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]))))
\ No newline at end of file
from pylab import *
from numpy import *
from scipy import linalg
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))
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
U,S,V = linalg.svd(M)
X = V[-1,:4]
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.")
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]
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]
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*')
def skew(a):
""" 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]])
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
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
def compute_fundamental_normalized(x1,x2):
""" 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.")
# 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)
# reverse normalization
F = dot(T1.T,dot(F,T2))
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
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).
input: x1,x2 (3*n arrays) points in hom. coordinates. """
from PCV.tools 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']
import matplotlib.delaunay as md
from scipy import ndimage
from pylab import *
from numpy import *
from PCV.geometry import homography
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
def combine_images(im1,im2,alpha):
""" 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
def triangulate_points(x,y):
""" 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')
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
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
from numpy import *
import pickle
# from pysqlite2 import dbapi2 as sqlite
import sqlite3 as sqlite
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()
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 tf_idf_dist(voc,v1,v2):
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)
# import PIL and pylab for plotting
from PIL import Image
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()
\ No newline at end of file
from numpy import *
from scipy.cluster.vq import *
from PCV.localdescriptors import sift
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 = []