Skip to content

Instantly share code, notes, and snippets.

@eshafeeqe
Created March 25, 2014 12:37
Show Gist options
  • Save eshafeeqe/8d505a052b41830ff197 to your computer and use it in GitHub Desktop.
Save eshafeeqe/8d505a052b41830ff197 to your computer and use it in GitHub Desktop.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys , os , random
import cv, cv2
import numpy as np
from numpy.lib.stride_tricks import as_strided
SURF_THRESHOLD = 1000
SURF_MATCH_THRESHOLD = 0.3
INLIER_DISTANCE_THRESHOLD = 30
NUMBER_OF_INLIERS_THRESHOLD = 0.8
def SURFdetector ( image ):
global SURF_THRESHOLD
detector = cv2.SURF(SURF_THRESHOLD, 10, 10)
keypoints, descriptors = detector.detectAndCompute(image, None)
return keypoints, descriptors
def SURFmatcher(keypoints_set , descriptors_set):
keypoint_1 = keypoints_set[0]
keypoint_2 = keypoints_set[1]
descriptor_1 = descriptors_set[0]
descriptor_2 = descriptors_set[1]
diff = descriptor_2 - descriptor_1[:,None]
squre_of_diff = diff ** 2
sum_square_diff = squre_of_diff.sum(axis=-1)
score = np.sqrt(sum_square_diff.min(axis=-1))
matches = np.argmin(sum_square_diff,axis=-1)
invalid_matches = score > SURF_MATCH_THRESHOLD
score[invalid_matches] = -1
matches[invalid_matches] = -1
return matches , score
def correspondace_map(features_1, features_2, matches):
def feature_to_point(keypoint):
x,y = keypoint.pt
return x , y
valid_matches_bool = matches != -1
valid_matches = matches[valid_matches_bool]
feature_to_array = np.vectorize(feature_to_point)
features1 = (np.asarray(feature_to_array(features_1))).T
features2 = (np.asarray(feature_to_array(features_2))).T
valid_features1 = features1[valid_matches_bool]
valid_features2 = features2[valid_matches,:]
correspondace = np.hstack((valid_features1, valid_features2))
return correspondace
def HomograpgyCalculation(correspondace, for_ransac=True):
def Making_A(point):
x,y,x1,y1 = point.item(0), point.item(1), point.item(2), point.item(3)
return x,y,1,0,0,0,-(x*x1),-(x1*y),0,0,0,x,y,1,-(x*y1),-(y*y1)
desired_correspondaces = 0
if for_ransac == True:
random_indices = np.array(random.sample(range(correspondace.shape[0]), 8))
desired_correspondaces = correspondace[random_indices,:]
else :
desired_correspondaces = correspondace
#print desired_correspondaces.shape
B = desired_correspondaces[:,2:4].flatten()
A = np.apply_along_axis(Making_A,1,desired_correspondaces).reshape(-1,8)
#Computing over determined solution using psudo inverse of A
A_psudo_inverse = np.dot ( np.linalg.inv( np.dot ( A.T , A ) ) , A.T )
x = np.dot ( A_psudo_inverse , B )
H = (np.append(x,1)).reshape(3,3)
return H
def RANSAC ( features_set, matches, H ):
features1 = features_set[0]
features2 = features_set[1]
correspondace = correspondace_map(features1 , features2 , matches)
number_of_correspondences = correspondace.shape[0]
inlier_set = 0
outliers = 0
best_inliers = 0
trials = 0
N = 1e3
max_inliers = 10
min_variance = 1e10
while trials < N :
temp_H = HomograpgyCalculation(correspondace)
correspondace_img1 = correspondace[:,0:2]
correspondace_img2 = correspondace[:,2:4]
#padding one at the end
correspondace_img = np.ones((correspondace_img1.shape[0],correspondace_img1.shape[1]+1))
correspondace_img[:,:-1] = correspondace_img1
#multiplying with temp
correspondace_mul_H = np.dot(temp_H,correspondace_img.T).T
correspondace_H = (np.divide(correspondace_mul_H[:,0:2].T,correspondace_mul_H[:,2])).T
#print np.dot(H,correspondace_img.T),T
diff = correspondace_img2 - correspondace_H
#print diff
error = (((diff[:,0]**2) + (diff[:,1]**2))**0.5)
inlier_indices = error < INLIER_DISTANCE_THRESHOLD
inlier_set = correspondace[inlier_indices,:]
number_of_inliers = inlier_set.shape[0]
#print number_of_inliers
if number_of_inliers > max_inliers :
error_mean = (error.sum(axis=-1))/number_of_inliers
variance = ((error**2).sum(axis=-1)) - error_mean
if variance < min_variance :
max_inliers = number_of_inliers
min_variance = variance
H = temp_H
best_inliers = inlier_set
outliers_indices = np.logical_not(inlier_indices)
outliers = correspondace[outliers_indices,:]
#Update N and no of trials
trials +=1
if number_of_inliers > 0 :
e = 1.0 - float ( number_of_inliers )/ float ( number_of_correspondences )
e_1 = 1.0 - e
if e_1 == 1:
break
if np . log (1.0 - e_1 * e_1 * e_1 * e_1 * e_1 * e_1 * e_1 * e_1 ) !=0:
N = int (np. log (1.0-0.99) /np . log (1.0- e_1 * e_1 * e_1 * e_1 * e_1 * e_1 * e_1 * e_1 ) )
if float ( number_of_inliers ) / float ( number_of_correspondences ) < NUMBER_OF_INLIERS_THRESHOLD \
and trials > N:
trials = 1
#print H
#print correspondace
#print best_inliers
Homograpgy_best_inliers = HomograpgyCalculation(best_inliers,for_ransac = False)
return best_inliers, outliers, Homograpgy_best_inliers
def main ( ) :
image =[]
keypoints_set = []
descriptors_set = []
H = 0
#loading images
for i in range (len ( sys . argv ) - 1 ):
filename = sys . argv [ i + 1 ]
image.append(cv2.imread (filename))
width = image[0].shape[1]
height = image[0].shape[0]
surf_features = np.zeros((height , len(image)*width ,3) ,np.uint8)
surf_matches = np.zeros((height , len(image)*width ,3) ,np.uint8)
for i in range(len(image)):
keypoints, descriptors = SURFdetector( image[i])
keypoints_set.append(keypoints)
descriptors_set.append(descriptors)
surf_features[0:height , i*width : (i + 1)*width , :] = image[i]
for j in range ( len ( keypoints ) ) :
x , y = keypoints[j].pt
cv2.circle ( surf_features , ( (i*width) + int (x) , int (y) ) , 0 , (255 , 0 , 0) , 4)
cv2.imwrite('surf_features.png', surf_features)
surf_best_inliers = surf_features
outliers_image = surf_features
matches, score = SURFmatcher(keypoints_set,descriptors_set)
best_inliers, outliers, Homograpgy_best_inliers = RANSAC(keypoints_set, matches, H)
for i in range ( len ( best_inliers ) ) :
cv2.line (surf_best_inliers , ((int(best_inliers[i][0])) , (int(best_inliers[i][1]))), \
((int(best_inliers[i][2]))+width , (int(best_inliers[i][3]))), \
(255*( i%4) ,255*(( i+1)%4) , 255*(( i+2)%4) ) , 1 , cv2.CV_AA, 0)
cv2.imwrite('surf_best_inliers.png', surf_best_inliers)
for i in range ( len ( outliers ) ) :
cv2.line (outliers_image , ((int(outliers[i][0])) , (int(best_inliers[i][1]))), \
((int(outliers[i][2]))+width , (int(outliers[i][3]))), \
(255 ,255 , 255 ) , 1 , cv2.CV_AA, 0)
cv2.imwrite('surf_outliers.png', outliers_image)
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment