In [ ]:

    

import cv2
import numpy as np
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import leastsq
#from multiprocessing import Queue
from queue import PriorityQueue
from scipy.ndimage import filters
import random as rnd
from numpy import *
from pylab import *
from pandas import *
%matplotlib inline


    

In [ ]:

    

def compute_harris_points(img, sigma=2):
    #compute derivates in the image 
    imx = np.zeros(img.size)
    imy = np.zeros(img.size)
    
    imx = filters.gaussian_filter(img, (sigma,sigma), (0,1))
    imy = filters.gaussian_filter(img, (sigma,sigma), (1,0))
    
    # compute the products of derivatives at every pixel
    Sxx = filters.gaussian_filter(imx*imx,sigma)
    Sxy = filters.gaussian_filter(imx*imy,sigma)
    Syy = filters.gaussian_filter(imy*imy,sigma)
    
    # determinant and trace
    Mdet = Sxx*Syy - Sxy**2
    Mtr = Sxx + Syy
    harris = np.divide(Mdet, Mtr)
    harris[np.isposinf(harris)] = 0
    harris[np.isnan(harris)] = 0
    return harris


def doHarrisNonMaxSupression(harrisim,min_dist=10,threshold=0.1):
    #Return corners from a Harris response image
    #min_dist is the minimum number of pixels separating
    #corners and image boundary. 
    global t 
    global dist
    dist=min_dist
    t=threshold
    #print(t)
    # find top corner candidates above a threshold
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # get coordinates of candidates
    coords = array(harrisim_t.nonzero()).T
    # ...and their values
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # sort candidates
    index = argsort(candidate_values)
    # store allowed point locations in array
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # select the best points taking min_distance into account
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
            (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0
    return filtered_coords

def plot_harris_points(image,filtered_coords):
    #""" Plots corners found in image. """
    plt.figure()
    gray()
    plt.imshow(image)
    plt.title('Harris corner detection, dist=%s and threshold=%s'%(dist,t))
    plt.plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*',color = 'r')
    plt.axis('off')
    plt.show()

def get_descriptors(image,filtered_coords,wid=5):

    desc = []
    for coords in filtered_coords:
        patch = image[coords[0]-wid:coords[0]+wid+1, coords[1]-wid:coords[1]+wid+1].flatten()
        desc.append(patch)
    return desc

def match(desc1,desc2,threshold=0.5):
    
    n = len(desc1[0])
    # pair-wise distances
    d = -ones((len(desc1),len(desc2)))
    for i in range(len(desc1)):
        for j in range(len(desc2)):
            d1 = (desc1[i] - mean(desc1[i])) / std(desc1[i])
            d2 = (desc2[j] - mean(desc2[j])) / std(desc2[j])
            ncc_value = sum(d1 * d2) / (n-1)
            if ncc_value > threshold:
                d[i,j] = ncc_value
    ndx = argsort(-d)
    matchscores = ndx[:,0]
    return matchscores

def match_twosided(desc1,desc2,threshold=0.5):
    matches_12 = match(desc1,desc2,threshold)
    matches_21 = match(desc2,desc1,threshold)
    ndx_12 = where(matches_12 >= 0)[0]
    # remove matches that are not symmetric
    for n in ndx_12:
        if matches_21[matches_12[n]] != n:
            matches_12[n] = -1
    return matches_12

def appendimages(im1,im2):

    # select the image with the fewest rows and fill in enough empty rows
    rows1 = im1.shape[0]
    rows2 = im2.shape[0]
    if rows1 < rows2:
        im1 = concatenate((im1,zeros((rows2-rows1,im1.shape[1]))),axis=0)
    elif rows1 > rows2:
        im2 = concatenate((im2,zeros((rows1-rows2,im2.shape[1]))),axis=0)
    # if none of these cases they are equal, no filling needed.
    return concatenate((im1,im2), axis=1)

def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):

    im3 = appendimages(im1,im2)
    if show_below:
        im3 = vstack((im3,im3))
    plt.imshow(im3)
    cols1 = im1.shape[1]
    for i,m in enumerate(matchscores):
        if m>0:
            plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
    axis('off')


    

In [ ]:

    

## harris corner detection
    
    #gray1_1= np.float32(gray1)
    #gray2_2= np.float32(gray2)
    #corners1 = cv2.goodFeaturesToTrack(gray1_1, 100, 0.01, 10)
    #corners1 = np.array(corners1)
    #corners2 = cv2.goodFeaturesToTrack(gray2_2, 100, 0.01, 10)
    #corners2 = np.array(corners2)
    
    
    
    #for corner in corners1:
    #    x,y = corner.ravel()
    #    cv2.circle(gray1_1,(x,y),3,255,-1)
    #plt.imshow(gray1_1)
    #plt.show()
    
    #for corner in corners2:
    #    x,y = corner.ravel()
    #    cv2.circle(gray2_2,(x,y),3,255,-1)
    #plt.imshow(gray2_2)
    #plt.show()
    
    #print(corners1)
    wid =5
    harrisim1 = compute_harris_points(gray1_1)
    filtered_coords1 = doHarrisNonMaxSupression(harrisim1)
    #plot_harris_points(gray1_1, filtered_coords1)
    #print(filtered_coords1)
    d1 = get_descriptors(gray1_1, filtered_coords1, wid)
    
    
    harrisim2 = compute_harris_points(gray2_2)
    filtered_coords2 = doHarrisNonMaxSupression(harrisim2)
    #plot_harris_points(gray2_2, filtered_coords2)
    #print(filtered_coords2)
    d2 = get_descriptors(gray2_2, filtered_coords2, wid)
    
    
    matches = match_twosided(d1,d2)
    
    plt.figure(figsize=(12,12))
    gray()
    plot_matches(gray1_1,gray2_2,filtered_coords1,filtered_coords2,matches)
    plt.show()
    #########