Created
March 30, 2014 12:38
-
-
Save eshafeeqe/355c5a783d0b36ca3bae to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/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 = 500 | |
| SURF_MATCH_THRESHOLD = 0.3 | |
| SURF_RATIO_THRESHOLD = 0.9 | |
| INLIER_DISTANCE_THRESHOLD = 20 | |
| 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) | |
| matching_score = np.sort(score,axis=-1)[:,:2] | |
| ratio = np.true_divide(matching_score[:,0],matching_score[:,1]) | |
| matches = np.argmin(sum_square_diff,axis=-1) | |
| invalid_matches = np.logical_or(matching_score[:,0]>SURF_MATCH_THRESHOLD ,ratio>SURF_RATIO_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 | |
| 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 | |
| diff = correspondace_img2 - correspondace_H | |
| 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] | |
| 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 | |
| Homograpgy_best_inliers = HomograpgyCalculation(best_inliers, for_ransac = False) | |
| return best_inliers, outliers, Homograpgy_best_inliers | |
| def homography_transformation(Homograpgy_best_inliers_set): | |
| L = len(Homograpgy_best_inliers_set) | |
| i = int ((L+1)/2) - 2 | |
| while i >= 0: | |
| Homograpgy_best_inliers_set[ i ] = np.dot(Homograpgy_best_inliers_set[ i] ,Homograpgy_best_inliers_set[ i+1 ] ) | |
| i-=1 | |
| i = int ((L+1)/2) + 2 | |
| while i < (L+1): | |
| Homograpgy_best_inliers_set[ i ] = np.dot(Homograpgy_best_inliers_set[ i] , Homograpgy_best_inliers_set[ i-1 ] ) | |
| i+=1 | |
| return 0 | |
| def boundary_calculation(Homograpgy_best_inliers_set, width, height): | |
| homography_transformation(Homograpgy_best_inliers_set) | |
| L = len(Homograpgy_best_inliers_set) | |
| S = np.asarray( Homograpgy_best_inliers_set ) | |
| H = np.array([[0,0,width-1,width-1],[0,height-1,0,height-1],[1,1,1,1]]) | |
| K = np.dot(S,H) | |
| U = (np.hstack((K[:,2,:],K[:,2,:]))).reshape(3,2,4) | |
| boundary_image = (np.divide(K[:,0:2,:],U)) | |
| image_boundaries_min = boundary_image.min(axis=-1) | |
| image_boundaries_max = boundary_image.max(axis=-1) | |
| boundaries = np.hstack((image_boundaries_min,image_boundaries_max)) | |
| image_boundaries = np.insert(boundaries,(L/2)+1,np.array([0,0,width-1,height-1]),0) | |
| cor_min = image_boundaries_min.min(axis=0) | |
| cor_max = image_boundaries_max.max(axis=0) | |
| final_image_boundaries = np.hstack((cor_min,cor_max)) | |
| return final_image_boundaries, image_boundaries | |
| def main ( ) : | |
| image =[] | |
| keypoints_set = [] | |
| descriptors_set = [] | |
| H = 0 | |
| best_inliers_set = [] | |
| Homograpgy_best_inliers_set = [] | |
| #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 | |
| number_of_images = len(image) | |
| for i in range(number_of_images-1): | |
| keypoints = keypoints_set[i:i+2] | |
| descriptors = descriptors_set[i:i+2] | |
| matches, score = 0,0 | |
| best_inliers, outliers, Homograpgy_best_inliers = 0, 0, 0 | |
| if i > ((number_of_images/2)-1): | |
| keypoints[1] , keypoints[0] = keypoints[0] , keypoints[1] | |
| descriptors[1] , descriptors[0] = descriptors[0] , descriptors[1] | |
| matches, score = SURFmatcher(keypoints,descriptors) | |
| best_inliers, outliers, Homograpgy_best_inliers = RANSAC(keypoints, matches, H) | |
| else: | |
| matches, score = SURFmatcher(keypoints,descriptors) | |
| best_inliers, outliers, Homograpgy_best_inliers = RANSAC(keypoints, matches, H) | |
| best_inliers_set.append(best_inliers) | |
| Homograpgy_best_inliers_set.append(Homograpgy_best_inliers) | |
| final_image_boundaries, image_boundaries = boundary_calculation(Homograpgy_best_inliers_set, width, height) | |
| fx_min = final_image_boundaries.item(0) | |
| fx_max = final_image_boundaries.item(2) | |
| fy_min = final_image_boundaries.item(1) | |
| fy_max = final_image_boundaries.item(3) | |
| mosaic_image = np.zeros((int(fy_max -fy_min ), int(fx_max - fx_min) ,3) ,np.uint8) | |
| mosaic_image2 = np.zeros((int(fy_max -fy_min + 20), int(fx_max - fx_min + 20) ,3) ,np.uint8) | |
| H_inverse = [] | |
| count = 0 | |
| for i in range(number_of_images): | |
| if i == number_of_images/2: | |
| H_inv = np.eye(3, dtype=float) | |
| H_inverse.append(H_inv) | |
| else : | |
| H_inv = np.linalg.inv(Homograpgy_best_inliers_set[count]) | |
| H_inverse.append(H_inv) | |
| count += 1 | |
| mosaic = np.indices((int(fy_max - fy_min),int(fx_max - fx_min))).swapaxes(0,2).swapaxes(0,1)[:,:,::-1] | |
| mosaic += [fx_min, fy_min] | |
| ones_to_add = np.ones((int(fy_max - fy_min),int(fx_max - fx_min))) | |
| new_mosaic_set = [] | |
| for i in range(len(image)): | |
| bool_mosaic = np.logical_and ((mosaic < [image_boundaries.item(i,2),image_boundaries.item(i,3)]).all(axis =-1) , \ | |
| (mosaic > [image_boundaries.item(i,0),image_boundaries.item(i,1)]).all(axis=-1)) | |
| bool_mosaic_neg = np.logical_not(bool_mosaic) | |
| new_mosaic = np.dstack((mosaic , ones_to_add)) | |
| desired = new_mosaic[bool_mosaic] | |
| new_mosaic[bool_mosaic_neg] = [0,0,0] | |
| image_points = np.dot(H_inverse[i], desired.T).T | |
| image_points_transformed = (image_points[:,:2].swapaxes(-2,-1)/image_points[:,2]).swapaxes(-1,-2) | |
| condition = np.logical_and((image_points_transformed < [width-1, height-1]).all(axis =-1) , \ | |
| (image_points_transformed > [0,0]).all(axis=-1)) | |
| condition_neg = np.logical_not(condition) | |
| image_points[condition_neg] =[0,0,0] | |
| X_Y = image_points_transformed[condition] | |
| X_Y_floor = np.floor(X_Y) | |
| diff = X_Y - X_Y_floor | |
| diff_X = diff[:,0] | |
| diff_Y = diff[:,1] | |
| indeces_0 = (X_Y_floor[:,::-1]).astype(int) | |
| indeces_1 = np.array([0 , 1]) + indeces_0 | |
| indeces_2 = indeces_0 + np.array([1 , 0]) | |
| indeces_3 = indeces_0 + np.array([1 , 1]) | |
| image_points_0 = image[i][indeces_0[:,0],indeces_0[:,1]] | |
| image_points_1 = image[i][indeces_1[:,0],indeces_1[:,1]] | |
| image_points_2 = image[i][indeces_2[:,0],indeces_2[:,1]] | |
| image_points_3 = image[i][indeces_3[:,0],indeces_3[:,1]] | |
| temp_0 = ((1-diff_X)*(1-diff_Y)*(image_points_1.T)).T | |
| temp_1 = ((diff_X)*(1-diff_Y)*(image_points_1.T)).T | |
| temp_2 = ((1-diff_X)*(diff_Y)*(image_points_2.T)).T | |
| temp_3 = ((diff_X * diff_Y)*(image_points_3.T)).T | |
| temp = temp_0 + temp_1 + temp_2 + temp_3 | |
| image_points[condition] = temp | |
| new_mosaic[bool_mosaic] = image_points | |
| new_mosaic_set.append(new_mosaic) | |
| cv2.imwrite('bool_mosaic'+str(i)+'.png', new_mosaic) | |
| #for i in range(len(new_mosaic_set)): | |
| #mosaic = new_mosaic_set[0] + 0.8*new_mosaic_set[1] + 0.8*new_mosaic_set[2] + 0.8*new_mosaic_set[3] | |
| condition1 = np.logical_and((new_mosaic_set[0]>np.array([0,0,0])).all(axis=-1),(new_mosaic_set[1]>np.array([0,0,0])).all(axis=-1)) | |
| condition2 = np.logical_and((new_mosaic_set[1]>np.array([0,0,0])).all(axis=-1),(new_mosaic_set[2]>np.array([0,0,0])).all(axis=-1)) | |
| condition3 = np.logical_and((new_mosaic_set[2]>np.array([0,0,0])).all(axis=-1),(new_mosaic_set[3]>np.array([0,0,0])).all(axis=-1)) | |
| new_mosaic_set[0][condition1] = 1*new_mosaic_set[0][condition1] | |
| new_mosaic_set[1][condition1] = 0*new_mosaic_set[1][condition1] | |
| new_mosaic_set[1][condition2] = 1*new_mosaic_set[1][condition2] | |
| new_mosaic_set[2][condition2] = 0*new_mosaic_set[2][condition2] | |
| new_mosaic_set[2][condition3] = 1*new_mosaic_set[2][condition3] | |
| new_mosaic_set[3][condition3] = 0*new_mosaic_set[3][condition3] | |
| mosaic = new_mosaic_set[0] + new_mosaic_set[1] + new_mosaic_set[2] + new_mosaic_set[3] | |
| cv2.imwrite('mosaiced_image.png', mosaic) | |
| main() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment