Mercurial > hg > nsaunier > traffic-intelligence

comparison scripts/compute-homography.py @ 614:5e09583275a4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged Nicolas/trafficintelligence into default
author Mohamed Gomaa <eng.m.gom3a@gmail.com>
date Fri, 05 Dec 2014 12:13:53 -0500
parents 9c429c7efe89
children aded6c1c2ebd
comparison
equal deleted inserted replaced
598:11f96bd08552 614:5e09583275a4
1 #! /usr/bin/env python
2
3 import sys, argparse
4
5 import matplotlib.pyplot as plt
6 import numpy as np
7 import cv2
8
9 import cvutils
10 import utils
11
12 parser = argparse.ArgumentParser(description='The program computes the homography matrix from at least 4 non-colinear point correspondences inputed in the same order in a video frame and a aerial photo/ground map, or from the list of corresponding points in the two planes.', epilog = '''The point correspondence file contains at least 4 non-colinear point coordinates
13 with the following format:
14 - the first two lines are the x and y coordinates in the projected space (usually world space)
15 - the last two lines are the x and y coordinates in the origin space (usually image space)
16
17 If providing video and world images, with a number of points to input
18 and a ration to convert pixels to world distance unit (eg meters per pixel),
19 the images will be shown in turn and the user should click
20 in the same order the corresponding points in world and image spaces.''', formatter_class=argparse.RawDescriptionHelpFormatter,)
21
22 parser.add_argument('-p', dest = 'pointCorrespondencesFilename', help = 'name of the text file containing the point correspondences')
23 parser.add_argument('-i', dest = 'videoFrameFilename', help = 'filename of the video frame')
24 parser.add_argument('-w', dest = 'worldFilename', help = 'filename of the aerial photo/ground map')
25 parser.add_argument('-n', dest = 'nPoints', help = 'number of corresponding points to input', default = 4, type = int)
26 parser.add_argument('-u', dest = 'unitsPerPixel', help = 'number of units per pixel', default = 1., type = float)
27 parser.add_argument('--display', dest = 'displayPoints', help = 'display original and projected points on both images', action = 'store_true')
28 parser.add_argument('--intrinsic', dest = 'intrinsicCameraMatrixFilename', help = 'name of the intrinsic camera file')
29 parser.add_argument('--distortion-coefficients', dest = 'distortionCoefficients', help = 'distortion coefficients', nargs = '*', type = float)
30 parser.add_argument('--undistorted-multiplication', dest = 'undistortedImageMultiplication', help = 'undistorted image multiplication', type = float)
31 parser.add_argument('--undistort', dest = 'undistort', help = 'undistort the video (because features have been extracted that way', action = 'store_true')
32 parser.add_argument('--save', dest = 'saveImages', help = 'save the undistorted video frame (display option must be chosen)', action = 'store_true')
33
34 args = parser.parse_args()
35
36 # TODO process camera intrinsic and extrinsic parameters to obtain image to world homography, taking example from Work/src/python/generate-homography.py script
37 # cameraMat = load(videoFilenamePrefix+'-camera.txt');
38 # T1 = cameraMat[3:6,:].copy();
39 # A = cameraMat[0:3,0:3].copy();
40
41 # # pay attention, rotation may be the transpose
42 # # R = T1[:,0:3].T;
43 # R = T1[:,0:3];
44 # rT = dot(R, T1[:,3]/1000);
45 # T = zeros((3,4),'f');
46 # T[:,0:3] = R[:];
47 # T[:,3] = rT;
48
49 # AT = dot(A,T);
50
51 # nPoints = 4;
52 # worldPoints = cvCreateMat(nPoints, 3, CV_64FC1);
53 # imagePoints = cvCreateMat(nPoints, 3, CV_64FC1);
54
55 # # extract homography from the camera calibration
56 # worldPoints = cvCreateMat(4, 3, CV_64FC1);
57 # imagePoints = cvCreateMat(4, 3, CV_64FC1);
58
59 # worldPoints[0,:] = [[1, 1, 0]];
60 # worldPoints[1,:] = [[1, 2, 0]];
61 # worldPoints[2,:] = [[2, 1, 0]];
62 # worldPoints[3,:] = [[2, 2, 0]];
63
64 # wPoints = [[1,1,2,2],
65 # [1,2,1,2],
66 # [0,0,0,0]];
67 # iPoints = utils.worldToImage(AT, wPoints);
68
69 # for i in range(nPoints):
70 # imagePoints[i,:] = [iPoints[:,i].tolist()];
71
72 # H = cvCreateMat(3, 3, CV_64FC1);
73
74 # cvFindHomography(imagePoints, worldPoints, H);
75
76
77 homography = np.array([])
78 if args.pointCorrespondencesFilename != None:
79 worldPts, videoPts = cvutils.loadPointCorrespondences(args.pointCorrespondencesFilename)
80 homography, mask = cv2.findHomography(videoPts, worldPts) # method=0, ransacReprojThreshold=3
81 elif args.videoFrameFilename != None and args.worldFilename != None:
82 worldImg = plt.imread(args.worldFilename)
83 videoImg = plt.imread(args.videoFrameFilename)
84 if args.undistort:
85 [map1, map2] = cvutils.computeUndistortMaps(videoImg.shape[1], videoImg.shape[0], args.undistortedImageMultiplication, np.loadtxt(args.intrinsicCameraMatrixFilename), args.distortionCoefficients)
86 videoImg = cv2.remap(videoImg, map1, map2, interpolation=cv2.INTER_LINEAR)
87 print('Click on {0} points in the video frame'.format(args.nPoints))
88 plt.figure()
89 plt.imshow(videoImg)
90 videoPts = np.array(plt.ginput(args.nPoints, timeout=3000))
91 print('Click on {0} points in the world image'.format(args.nPoints))
92 plt.figure()
93 plt.imshow(worldImg)
94 worldPts = args.unitsPerPixel*np.array(plt.ginput(args.nPoints, timeout=3000))
95 plt.close('all')
96 homography, mask = cv2.findHomography(videoPts, worldPts)
97 # save the points in file
98 f = open('point-correspondences.txt', 'a')
99 np.savetxt(f, worldPts.T)
100 np.savetxt(f, videoPts.T)
101 f.close()
102
103 if homography.size>0:
104 np.savetxt('homography.txt',homography)
105
106 if args.displayPoints and args.videoFrameFilename != None and args.worldFilename != None and homography.size>0:
107 worldImg = cv2.imread(args.worldFilename)
108 videoImg = cv2.imread(args.videoFrameFilename)
109 if args.undistort:
110 [map1, map2] = cvutils.computeUndistortMaps(videoImg.shape[1], videoImg.shape[0], args.undistortedImageMultiplication, np.loadtxt(args.intrinsicCameraMatrixFilename), args.distortionCoefficients)
111 videoImg = cv2.remap(videoImg, map1, map2, interpolation=cv2.INTER_LINEAR)
112 if args.saveImages:
113 cv2.imwrite(utils.removeExtension(args.videoFrameFilename)+'-undistorted.png', videoImg)
114 invHomography = np.linalg.inv(homography)
115 projectedWorldPts = cvutils.projectArray(invHomography, worldPts.T).T
116 projectedVideoPts = cvutils.projectArray(homography, videoPts.T).T
117 for i in range(worldPts.shape[0]):
118 # world image
119 cv2.circle(worldImg,tuple(np.int32(np.round(worldPts[i]/args.unitsPerPixel))),2,cvutils.cvBlue)
120 cv2.circle(worldImg,tuple(np.int32(np.round(projectedVideoPts[i]/args.unitsPerPixel))),2,cvutils.cvRed)
121 cv2.putText(worldImg, str(i+1), tuple(np.int32(np.round(worldPts[i]/args.unitsPerPixel))+5), cv2.FONT_HERSHEY_PLAIN, 2., cvutils.cvBlue, 2)
122 # video image
123 cv2.circle(videoImg,tuple(np.int32(np.round(videoPts[i]))),2,cvutils.cvBlue)
124 cv2.circle(videoImg,tuple(np.int32(np.round(projectedWorldPts[i]))),2,cvutils.cvRed)
125 cv2.putText(videoImg, str(i+1), tuple(np.int32(np.round(videoPts[i])+5)), cv2.FONT_HERSHEY_PLAIN, 2., cvutils.cvBlue, 2)
126 cv2.imshow('video frame',videoImg)
127 cv2.imshow('world image',worldImg)
128 cv2.waitKey()
129 cv2.destroyAllWindows()