diff --git a/RANSAC/RANSAC.py b/RANSAC/RANSAC.py new file mode 100644 index 0000000..31128dc --- /dev/null +++ b/RANSAC/RANSAC.py @@ -0,0 +1,68 @@ +import numpy as np +import matplotlib.pyplot as plt +import skimage.transform as skt +from homography import * + +def apply_homography(X, H): + Xprime = (H.dot(X.T)).T + Xprime/=Xprime[:,2][:,np.newaxis] + return Xprime + +def RANSAC(number_of_iterations,matches,n,r,d): + H_best = np.array([[1,0,0],[0,1,0],[0,0,1]]) + list_of_inliers = np.zeros((2, d, 2)) + # for 3D projection + matches = np.insert(matches, 2, 1, axis=2) + for i in range(number_of_iterations): + # swap the coords so they're xy instead of yx, what a great + matches[:,:,[0, 1]] = matches[:,:,[1, 0]] + # 1. Select a random sample of length n from the matches + indices = np.random.choice(np.arange(0, matches.shape[1]), size=n, replace=False) + sample = matches[:,indices,:] + # 2. Compute a homography based on these points using the methods given above + H = generate_homography(sample[0], sample[1], n=4) + # 3. Apply this homography to the remaining points that were not randomly selected + pred_location = apply_homography(matches[0], H) + # 4. Compute the residual between observed and predicted feature locations + residuals = np.sqrt((matches[1,:,0] - pred_location[:,0])**2 + (matches[1,:,1] - pred_location[:,1])**2) + # 5. Flag predictions that lie within a predefined distance r from observations as inliers + inliers = matches[:,residuals<=r] + # 6. If number of inliers is greater than the previous best + # and greater than a minimum number of inliers d, + # 7. update H_best + # 8. update list_of_inliers + if (inliers.shape[1] > list_of_inliers.shape[1]): + H_best = H + list_of_inliers = inliers + return H_best, list_of_inliers[:,:,:2].astype(int) + +I_1 = plt.imread('photo_1.jpg') +I_2 = plt.imread('photo_2.jpg') + +I_1 = I_1.mean(axis=2) +I_2 = I_2.mean(axis=2) + +corners_1 = harris_corner_detection(I_1) +corners_2 = harris_corner_detection(I_2) + +descriptors_1 = extract_descriptors(I_1, corners_1, 21) +descriptors_2 = extract_descriptors(I_2, corners_2, 21) + +matching_indices = get_matches(descriptors_1, descriptors_2 , .5).astype(int) + +matching_corners_1 = corners_1[matching_indices[:,0]] +matching_corners_2 = corners_2[matching_indices[:,1]] + +matches = np.array([matching_corners_1, matching_corners_2]) + +H_best, inliers = RANSAC(10000, matches, 10, 30, 4) + +# Create a projective transform based on the homography matrix $H$ +proj_trans = skt.ProjectiveTransform(H_best) + +# Warp the image into image 1's coordinate system +I_2_transformed = skt.warp(I_2,proj_trans) + +plt.imshow(I_1) +plt.imshow(I_2_transformed, alpha=.5) +plt.show() diff --git a/RANSAC/RANSAC.py~ b/RANSAC/RANSAC.py~ new file mode 100644 index 0000000..5a644a5 --- /dev/null +++ b/RANSAC/RANSAC.py~ @@ -0,0 +1,68 @@ +import numpy as np +import matplotlib.pyplot as plt +import skimage.transform as skt +from homography import * + +def apply_homography(X, H): + Xprime = (H.dot(X.T)).T + Xprime/=Xprime[:,2][:,np.newaxis] + return Xprime + +def RANSAC(number_of_iterations,matches,n,r,d): + H_best = np.array([[1,0,0],[0,1,0],[0,0,1]]) + list_of_inliers = np.zeros((2, d, 2)) + matches = np.insert(matches, 2, 1, axis=2) + for i in range(number_of_iterations): + # 1. Select a random sample of length n from the matches + indices = np.random.choice(np.arange(0, matches.shape[1]), size=n, replace=False) + sample = matches[:,indices,:] + # for 3D translation + # 2. Compute a homography based on these points using the methods given above + H = generate_homography(sample[0], sample[1], n=4) + # 3. Apply this homography to the remaining points that were not randomly selected + pred_location = apply_homography(matches[0], H) + # 4. Compute the residual between observed and predicted feature locations + residuals = np.sqrt((matches[1,:,0] - pred_location[:,0])**2 + (matches[1,:,1] - pred_location[:,1])**2) + # 5. Flag predictions that lie within a predefined distance r from observations as inliers + inliers = matches[:,residuals<=r] + # 6. If number of inliers is greater than the previous best + # and greater than a minimum number of inliers d, + # 7. update H_best + # 8. update list_of_inliers + if (inliers.shape[1] > list_of_inliers.shape[1]): + H_best = H + list_of_inliers = inliers + return H_best, list_of_inliers + +I_1 = plt.imread('photo_1.jpg') +I_2 = plt.imread('photo_2.jpg') + +I_1 = I_1.mean(axis=2) +I_2 = I_2.mean(axis=2) + +corners_1 = harris_corner_detection(I_1) +corners_2 = harris_corner_detection(I_2) + +descriptors_1 = extract_descriptors(I_1, corners_1, 21) +descriptors_2 = extract_descriptors(I_2, corners_2, 21) + +matching_indices = get_min_descriptor_errors(descriptors_1, descriptors_2 , .7).astype(int) + +matching_corners_1 = corners_1[matching_indices[:,0]] +matching_corners_2 = corners_2[matching_indices[:,1]] + +matches = np.array([matching_corners_1, matching_corners_2]) + +H_best, inliers = RANSAC(1000, matches, 10, 30, 1) + +H_best = np.random.rand(3,3) +# Create a projective transform based on the homography matrix $H$ +proj_trans = skt.ProjectiveTransform(H_best) + +# Warp the image into image 1's coordinate system +I_2_transformed = skt.warp(I_2,proj_trans) +print(I_2_transformed) +print(I_2) +plt.imshow(I_2) +plt.imshow(I_2_transformed, alpha=1) +plt.show() diff --git a/RANSAC/feature_matching.py b/RANSAC/feature_matching.py new file mode 100644 index 0000000..ce0d3ca --- /dev/null +++ b/RANSAC/feature_matching.py @@ -0,0 +1,68 @@ +import numpy as np +from keypoint_detection import * + +def pad_with(vector, pad_width, iaxis, kwargs): + pad_value = kwargs.get('padder', 0) + vector[:pad_width[0]] = pad_value + vector[-pad_width[1]:] = pad_value + return vector + +def get_matches(descriptors_1, descriptors_2, r): + matches = np.zeros((descriptors_1.shape[0], 3)) + for i in range(0, descriptors_1.shape[0]): + error = float('inf') + matches[i][0] = i + for j in range(0, descriptors_2.shape[0]): + test_error = np.sum((descriptors_1[i] - descriptors_2[j])**2) + if (test_error < error): + matches[i][1] = j + # is robust + matches[i][2] = test_error < error * r + error = test_error + return matches[matches[:,2] == 1][:,:2].astype(int) + +def extract_descriptors(I, corners, l): + descriptors = np.zeros((corners.shape[0], l/2*2, l/2*2)) + I = np.pad(I, l/2 + 1, pad_with) + for i in range(0, corners.shape[0]): + # uuh might have indices swapped + x = corners[i][0] + l/2 + y = corners[i][1] + l/2 + descriptor = I[int(x-l/2):int(x+l/2),int(y-l/2):int(y+l/2)] + descriptors[i] = descriptor + return descriptors + +''' +I_1 = plt.imread('photo_1.jpg') +I_2 = plt.imread('photo_2.jpg') + +I_1 = I_1.mean(axis=2) +I_2 = I_2.mean(axis=2) + +corners_1 = harris_corner_detection(I_1) +corners_2 = harris_corner_detection(I_2) + +descriptors_1 = extract_descriptors(I_1, corners_1, 21) +descriptors_2 = extract_descriptors(I_2, corners_2, 21) + +img = np.zeros((I_1.shape[0], I_1.shape[1] * 2)) +img[:,I_1.shape[1]:] = I_2 +img[:,:I_1.shape[1]] = I_1 + +matches = get_matches(descriptors_1, descriptors_2 ,.09) + +for match in matches: + x1 = corners_1[match[0]][1] + y1 = corners_1[match[0]][0] + + x2 = I_2.shape[1] + corners_2[match[1]][1] + y2 = corners_2[match[1]][0] + + plt.plot([x1, x2], [y1, y2]) + +plt.scatter(corners_1[:,1], corners_1[:,0], c='blue') +plt.scatter(I_1.shape[1] + corners_2[:,1], corners_2[:,0], c='white') +plt.imshow(img) +plt.show() +''' + diff --git a/RANSAC/feature_matching.pyc b/RANSAC/feature_matching.pyc new file mode 100644 index 0000000..86d3de4 Binary files /dev/null and b/RANSAC/feature_matching.pyc differ diff --git a/RANSAC/homography.py b/RANSAC/homography.py new file mode 100644 index 0000000..b1db2f4 --- /dev/null +++ b/RANSAC/homography.py @@ -0,0 +1,60 @@ +import numpy as np +import matplotlib.pyplot as plt +from feature_matching import * + +def make_A_row(corner_1, corner_2): + u1 = corner_1[0] + uprime = corner_2[0] + v1 = corner_1[1] + vprime = corner_2[1] + A_row = np.array([ + [0,0,0,-u1,-v1,-1,vprime*u1,vprime*v1,vprime], + [u1,v1,1,0,0,0,-uprime*u1,-uprime*v1,-uprime] + ]) + return A_row + +def generate_A(corners_1, corners_2, n): + A = np.zeros((2*n, 9)) + for i in range(0, 2*n, 2): + A[i:i+2] = make_A_row(corners_1[i/2], corners_2[i/2]) + return A + +def generate_homography(corners_1, corners_2, n): + A = generate_A(corners_1, corners_2, n) + U,Sigma,Vt = np.linalg.svd(A) + return np.reshape(Vt[-1], (3, 3)) + +X = np.array([[0,0,1], + [1,0,1], + [1,1,1], + [0,1,1], + [0,0,1]]) + +H = np.random.rand(3,3) +#H/= H[2,2] + +Xprime = (H.dot(X.T)).T +Xprime/=Xprime[:,2][:,np.newaxis] + +H_gen = generate_homography(X, Xprime, 4) +''' +I_1 = plt.imread('photo_1.jpg') +I_2 = plt.imread('photo_2.jpg') + +I_1 = I_1.mean(axis=2) +I_2 = I_2.mean(axis=2) + +corners_1 = harris_corner_detection(I_1) +corners_2 = harris_corner_detection(I_2) + +descriptors_1 = extract_descriptors(I_1, corners_1, 21) +descriptors_2 = extract_descriptors(I_2, corners_2, 21) + +matches = get_min_descriptor_errors(descriptors_1, descriptors_2, .09).astype(int) + +matching_corners_1 = corners_1[matches[:,0],:] +matching_corners_2 = corners_2[matches[:,1],:] + + +H = generate_H(matching_corners_1, matching_corners_2) +''' diff --git a/RANSAC/homography.pyc b/RANSAC/homography.pyc new file mode 100644 index 0000000..a0b1469 Binary files /dev/null and b/RANSAC/homography.pyc differ diff --git a/RANSAC/homography.py~ b/RANSAC/homography.py~ new file mode 100644 index 0000000..802f07c --- /dev/null +++ b/RANSAC/homography.py~ @@ -0,0 +1,48 @@ +import numpy as np +import matplotlib.pyplot as plt +from feature_matching import * + +def make_A_row(corner_1, corner_2): + + u1 = corner_1[0] + u2 = corner_2[0] + v1 = corner_1[1] + v2 = corner_1[1] + A_row = np.array([ + [0,0,0,-u1,-v1,-1,v2*u1,v2*v1,v1], + [u1,v1,1,0,0,0,-u2*u1,u1*v1,-u1] + ]) + return A_row + +def generate_A(corners_1, corners_2, n): + A = np.zeros((2*n, 9)) + for i in range(0, 2*n, 2): + A[i:i+2] = make_A_row(corners_1[i/2], corners_2[i/2]) + return A + +def generate_homography(corners_1, corners_2, n): + A = generate_A(corners_1, corners_2, n) + U,Sigma,Vt = np.linalg.svd(A) + return np.reshape(Vt[-1], (3, 3)) + +''' +I_1 = plt.imread('photo_1.jpg') +I_2 = plt.imread('photo_2.jpg') + +I_1 = I_1.mean(axis=2) +I_2 = I_2.mean(axis=2) + +corners_1 = harris_corner_detection(I_1) +corners_2 = harris_corner_detection(I_2) + +descriptors_1 = extract_descriptors(I_1, corners_1, 21) +descriptors_2 = extract_descriptors(I_2, corners_2, 21) + +matches = get_min_descriptor_errors(descriptors_1, descriptors_2, .09).astype(int) + +matching_corners_1 = corners_1[matches[:,0],:] +matching_corners_2 = corners_2[matches[:,1],:] + + +H = generate_H(matching_corners_1, matching_corners_2) +''' diff --git a/RANSAC/keypoint_detection.py b/RANSAC/keypoint_detection.py new file mode 100644 index 0000000..370524b --- /dev/null +++ b/RANSAC/keypoint_detection.py @@ -0,0 +1,57 @@ +import numpy as np +from scipy import signal +from matplotlib import pyplot as plt + +def pull_local_maxima(I, evaluation_fraction=.05): + local_maxima = [] + for i in range(1, I.shape[0] - 1): + for j in range(1, I.shape[1] - 1): + if (I[i + 1][j] <= I[i][j] + and I[i - 1][j] <= I[i][j] + and I[i][j + 1] <= I[i][j] + and I[i + 1][j - 1] <= I[i][j]): + local_maxima.append((i, j, I[i, j])) + + local_maxima.sort(key=lambda local_maxima: local_maxima[2], reverse=True) + local_maxima = np.array(local_maxima) + evaluation_slice = int(local_maxima.shape[0] * evaluation_fraction) + local_maxima = local_maxima[:evaluation_slice] + radii = np.zeros((local_maxima.shape[0], 1)) + float("inf") + local_maxima = np.hstack((local_maxima, radii)) + + for i in range(0, local_maxima.shape[0]): + for j in range(0, local_maxima.shape[0]): + distance = np.sqrt( + (local_maxima[i][0] - local_maxima[j][0]) ** 2 + (local_maxima[i][1] - local_maxima[j][1]) ** 2) + if (local_maxima[j][2] > local_maxima[i][2] and distance < local_maxima[i][3]): + local_maxima[i][3] = distance + + indices = np.argsort(-local_maxima[:, 3]) + local_maxima = local_maxima[indices] + return local_maxima[:100] + + +def harris_corner_detection(I, evaluation_fraction=.05): + Su = np.matrix( + [[-1, 0, 1], + [-2, 0, 2], + [-1, 0, 1]]) + w = np.matrix([ + [0.023528, 0.033969, 0.038393, 0.033969, 0.023528], + [0.033969, 0.049045, 0.055432, 0.049045, 0.033969], + [0.038393, 0.055432, 0.062651, 0.055432, 0.038393], + [0.033969, 0.049045, 0.055432, 0.049045, 0.033969], + [0.023528, 0.033969, 0.038393, 0.033969, 0.023528] + ]) + + Iu = signal.convolve2d(I, Su) + Iv = signal.convolve2d(I, Su.T) + + Iuu = signal.convolve2d(np.multiply(Iu, Iu), w) + Ivv = signal.convolve2d(np.multiply(Iv, Iv), w) + Iuv = signal.convolve2d(np.multiply(Iu, Iv), w) + + H = np.divide(np.multiply(Iuu, Ivv) - np.multiply(Iuv, Iuv), Iuu + Ivv + 1e-10) + + local_maxima = pull_local_maxima(H, evaluation_fraction) + return local_maxima[:,:2] diff --git a/RANSAC/keypoint_detection.pyc b/RANSAC/keypoint_detection.pyc new file mode 100644 index 0000000..b23e178 Binary files /dev/null and b/RANSAC/keypoint_detection.pyc differ diff --git a/RANSAC/photo_1.jpg b/RANSAC/photo_1.jpg new file mode 100644 index 0000000..4778723 Binary files /dev/null and b/RANSAC/photo_1.jpg differ diff --git a/RANSAC/photo_2.jpg b/RANSAC/photo_2.jpg new file mode 100644 index 0000000..58a5cbe Binary files /dev/null and b/RANSAC/photo_2.jpg differ diff --git a/RANSAC/test.py b/RANSAC/test.py new file mode 100644 index 0000000..a69bb45 --- /dev/null +++ b/RANSAC/test.py @@ -0,0 +1,15 @@ +import numpy as np +import matplotlib.pyplot as plt + + +X = np.array([0,0,1]) + +H = np.random.rand(3,3) +#H/= H[2,2] + +Xprime = (H.dot(X.T)).T +Xprime/=Xprime[:,2][:,np.newaxis] + +plt.plot(X[:,0],X[:,1],'g-') +plt.plot(Xprime[:,0],Xprime[:,1],'b-') +plt.show() diff --git a/README.md b/README.md index 3703ad8..8223632 100644 --- a/README.md +++ b/README.md @@ -1 +1,2 @@ -# project_2 \ No newline at end of file +# project_2 +This is a panoramic stitcher I made for Computer Vision diff --git a/convolution/convoluution.py b/convolution/convoluution.py new file mode 100644 index 0000000..10d6ade --- /dev/null +++ b/convolution/convoluution.py @@ -0,0 +1,90 @@ +import numpy as np +import matplotlib.pyplot as plt + +def convolve_no_edge(g,h): + gx = g.shape[0] + gy = g.shape[1] + + hx = h.shape[0] + hy = h.shape[1] + padx = hx/2 + pady = hy/2 + + res = np.zeros((gx - (padx - 1), gy - (pady - 1))) + for x in range (padx, gx - padx): + for y in range(pady, gy - pady): + g_local = g[x-padx:x+padx+1,y-pady:y+pady+1] + res_pixel = np.sum(np.multiply(g_local, h)) + res[x-padx,y-pady] = res_pixel + return res + +def convolve_fill(g,h, fill_val): + gx = g.shape[0] + gy = g.shape[1] + + hx = h.shape[0] + hy = h.shape[1] + padx = hx/2 + pady = hy/2 + + res = np.zeros((gx, gy)) + for x in range (0, gx): + for y in range(0, gy): + g_local = np.array([]) + for j in range(x - padx, x + padx + 1): + for k in range(y - pady, y + pady + 1): + if (j < 0 or j >= gx or k < 0 or k >= gy): + g_local = np.append(g_local, fill_val) + else: + g_local = np.append(g_local, g[j,k]) + g_local = g_local.reshape((3,3)) + res_pixel = np.sum(np.multiply(g_local, h)) + res[x-padx,y-pady] = res_pixel + return res + +def convolve_nearest_neighbors(g, h): + gx = g.shape[0] + gy = g.shape[1] + + hx = h.shape[0] + hy = h.shape[1] + padx = hx/2 + pady = hy/2 + + res = np.zeros((gx, gy)) + for x in range (0, gx): + for y in range(0, gy): + g_local = np.array([]) + for j in range(x - padx, x + padx + 1): + for k in range(y - pady, y + pady + 1): + if (j < 0 or j >= gx or k < 0 or k >= gy): + j_copy = j + k_copy = k + while ((j_copy < 0 or j_copy >= gx) or (k_copy < 0 or k_copy >= gy)): + if (j_copy < 0): + j_copy+=1 + elif (j_copy >= gx): + j_copy-=1 + if (k_copy < 0): + k_copy+=1 + elif (k_copy >= gy): + k_copy-=1 + g_local = np.append(g_local, g[j_copy, k_copy]) + else: + g_local = np.append(g_local, g[j,k]) + g_local = g_local.reshape((3,3)) + res_pixel = np.sum(np.multiply(g_local, h)) + print(res_pixel) + res[x-padx,y-pady] = res_pixel + return res + +img_color = plt.imread('stuff.jpeg') +I_gray = img_color.mean(axis=2) +Su = np.matrix( +[[-1, 0, 1], +[-2, 0, 2], +[-1,0,1]]) + +filter = convolve_nearest_neighbors(I_gray, Su.T) +plt.imshow(filter, cmap='gray') +plt.show() diff --git a/convolution/stuff.jpeg b/convolution/stuff.jpeg new file mode 100644 index 0000000..446397b Binary files /dev/null and b/convolution/stuff.jpeg differ diff --git a/feature_match/__init__.py b/feature_match/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/feature_match/feature_matching.py b/feature_match/feature_matching.py new file mode 100644 index 0000000..f11fd60 --- /dev/null +++ b/feature_match/feature_matching.py @@ -0,0 +1,12 @@ +import numpy as np +from matplotlib import pyplot as plt +from feature_match.keypoint_detection import harris_corner_detection, pull_local_maxima + +I = plt.imread('statue_of_liberty.jpg') +I = I.mean(axis=2) + + + + + + diff --git a/feature_match/feature_matching_dgk.py b/feature_match/feature_matching_dgk.py new file mode 100644 index 0000000..a2dc0b5 --- /dev/null +++ b/feature_match/feature_matching_dgk.py @@ -0,0 +1,28 @@ +import numpy as np +from matplotlib import pyplot as plt +from skimage.transform import resize +from feature_match.keypoint_detection import harris_corner_detection + +I_1 = plt.imread('outside_1.jpg') +I_2 = plt.imread('outside_2.jpg') + +I_1 = I_1.mean(axis=2) +I_2 = I_2.mean(axis=2) + +I_1 = resize(I_1, (I_1.shape[0] / 12, I_1.shape[1] / 12), anti_aliasing=True, order=3) +I_2 = resize(I_2, (I_2.shape[0] / 12, I_2.shape[1] / 12), anti_aliasing=True, order=3) + +h_1 = harris_corner_detection(I_1) +h_2 = harris_corner_detection(I_2) + +plt.imshow(I_1) +plt.scatter(h_1[:, 1], h_1[:, 0], c=h_1[:, 2]) +plt.show() + +plt.imshow(I_2) +plt.scatter(h_2[:, 1], h_2[:, 0], c=h_2[:, 2]) +plt.show() + + + + diff --git a/feature_match/keypoint_detection.py b/feature_match/keypoint_detection.py new file mode 100644 index 0000000..d8ac737 --- /dev/null +++ b/feature_match/keypoint_detection.py @@ -0,0 +1,57 @@ +import numpy as np +from scipy import signal +from matplotlib import pyplot as plt + +def pull_local_maxima(I, evaluation_fraction=1.): + local_maxima = [] + for i in range(1, I.shape[0] - 1): + for j in range(1, I.shape[1] - 1): + if (I[i + 1][j] <= I[i][j] + and I[i - 1][j] <= I[i][j] + and I[i][j + 1] <= I[i][j] + and I[i + 1][j - 1] <= I[i][j]): + local_maxima.append((i, j, I[i, j])) + + local_maxima.sort(key=lambda local_maxima: local_maxima[2], reverse=True) + local_maxima = np.array(local_maxima) + evaluation_slice = int(local_maxima.shape[0] * evaluation_fraction) + local_maxima = local_maxima[:evaluation_slice] + radii = np.zeros((local_maxima.shape[0], 1)) + float("inf") + local_maxima = np.hstack((local_maxima, radii)) + + for i in range(0, local_maxima.shape[0]): + for j in range(0, local_maxima.shape[0]): + distance = np.sqrt( + (local_maxima[i][0] - local_maxima[j][0]) ** 2 + (local_maxima[i][1] - local_maxima[j][1]) ** 2) + if (local_maxima[j][2] > local_maxima[i][2] and distance < local_maxima[i][3]): + local_maxima[i][3] = distance + + indices = np.argsort(-local_maxima[:, 3]) + local_maxima = local_maxima[indices] + return local_maxima[:100] + + +def harris_corner_detection(I, evaluation_fraction=1.): + Su = np.matrix( + [[-1, 0, 1], + [-2, 0, 2], + [-1, 0, 1]]) + w = np.matrix([ + [0.023528, 0.033969, 0.038393, 0.033969, 0.023528], + [0.033969, 0.049045, 0.055432, 0.049045, 0.033969], + [0.038393, 0.055432, 0.062651, 0.055432, 0.038393], + [0.033969, 0.049045, 0.055432, 0.049045, 0.033969], + [0.023528, 0.033969, 0.038393, 0.033969, 0.023528] + ]) + + Iu = signal.convolve2d(I, Su) + Iv = signal.convolve2d(I, Su.T) + + Iuu = signal.convolve2d(np.multiply(Iu, Iu), w) + Ivv = signal.convolve2d(np.multiply(Iv, Iv), w) + Iuv = signal.convolve2d(np.multiply(Iu, Iv), w) + + H = np.divide(np.multiply(Iuu, Ivv) - np.multiply(Iuv, Iuv), Iuu + Ivv + 1e-10) + + local_maxima = pull_local_maxima(H, evaluation_fraction) + return local_maxima diff --git a/feature_match/outside_1.jpg b/feature_match/outside_1.jpg new file mode 100644 index 0000000..d159f8e Binary files /dev/null and b/feature_match/outside_1.jpg differ diff --git a/feature_match/outside_2.jpg b/feature_match/outside_2.jpg new file mode 100644 index 0000000..4b30112 Binary files /dev/null and b/feature_match/outside_2.jpg differ diff --git a/feature_match/photo_1.jpg b/feature_match/photo_1.jpg new file mode 100644 index 0000000..4778723 Binary files /dev/null and b/feature_match/photo_1.jpg differ diff --git a/feature_match/photo_2.jpg b/feature_match/photo_2.jpg new file mode 100644 index 0000000..58a5cbe Binary files /dev/null and b/feature_match/photo_2.jpg differ diff --git a/feature_matching/feature_matching.py b/feature_matching/feature_matching.py new file mode 100644 index 0000000..00abec5 --- /dev/null +++ b/feature_matching/feature_matching.py @@ -0,0 +1,66 @@ +import numpy as np +from keypoint_detection import * + +def pad_with(vector, pad_width, iaxis, kwargs): + pad_value = kwargs.get('padder', 0) + vector[:pad_width[0]] = pad_value + vector[-pad_width[1]:] = pad_value + return vector + +def get_min_descriptor_errors(descirptors_1, descriptors_2, r): + descriptor_with_least_error = np.zeros((descriptors_1.shape[0], 2)) + for i in range(0, descriptors_1.shape[0]): + error = float('inf') + for j in range(0, descriptors_2.shape[0]): + if (i != j): + test_error = np.sum((descriptors_1[i] - descriptors_2[j])**2) + if (test_error < error): + descriptor_with_least_error[i][0] = j + # is robust + descriptor_with_least_error[i][1] = test_error < error * r + error = test_error + return descriptor_with_least_error + +def extract_descriptors(I, corners, l): + descriptors = np.zeros((corners.shape[0], l/2*2, l/2*2)) + I = np.pad(I, l/2 + 1, pad_with) + for i in range(0, corners.shape[0]): + # uuh might have indices swapped + x = corners[i][0] + l/2 + y = corners[i][1] + l/2 + descriptor = I[int(x-l/2):int(x+l/2),int(y-l/2):int(y+l/2)] + descriptors[i] = descriptor + return descriptors + +I_1 = plt.imread('photo_1.jpg') +I_2 = plt.imread('photo_2.jpg') + +I_1 = I_1.mean(axis=2) +I_2 = I_2.mean(axis=2) + +corners_1 = harris_corner_detection(I_1) +corners_2 = harris_corner_detection(I_2) + +descriptors_1 = extract_descriptors(I_1, corners_1, 21) +descriptors_2 = extract_descriptors(I_2, corners_2, 21) + +img = np.zeros((I_1.shape[0], I_1.shape[1] * 2)) +img[:,I_1.shape[1]:] = I_2 +img[:,:I_1.shape[1]] = I_1 + +descriptor_with_least_error = get_min_descriptor_errors(descriptors_1, descriptors_2 ,.09) + +for i in range(0, descriptor_with_least_error.shape[0]): + x1 = corners_1[i][1] + y1 = corners_1[i][0] + + x2 = I_2.shape[1] + corners_2[int(descriptor_with_least_error[i][0])][1] + y2 = corners_2[int(descriptor_with_least_error[i][0])][0] + if (descriptor_with_least_error[i][1]): + plt.plot([x1, x2], [y1, y2]) + +plt.scatter(corners_1[:,1], corners_1[:,0], c='blue') +plt.scatter(I_1.shape[1] + corners_2[:,1], corners_2[:,0], c='white') +plt.imshow(img) +plt.show() + diff --git a/feature_matching/keypoint_detection.py b/feature_matching/keypoint_detection.py new file mode 100644 index 0000000..cfeccc2 --- /dev/null +++ b/feature_matching/keypoint_detection.py @@ -0,0 +1,57 @@ +import numpy as np +from scipy import signal +from matplotlib import pyplot as plt + +def pull_local_maxima(I, evaluation_fraction=.05): + local_maxima = [] + for i in range(1, I.shape[0] - 1): + for j in range(1, I.shape[1] - 1): + if (I[i + 1][j] <= I[i][j] + and I[i - 1][j] <= I[i][j] + and I[i][j + 1] <= I[i][j] + and I[i + 1][j - 1] <= I[i][j]): + local_maxima.append((i, j, I[i, j])) + + local_maxima.sort(key=lambda local_maxima: local_maxima[2], reverse=True) + local_maxima = np.array(local_maxima) + evaluation_slice = int(local_maxima.shape[0] * evaluation_fraction) + local_maxima = local_maxima[:evaluation_slice] + radii = np.zeros((local_maxima.shape[0], 1)) + float("inf") + local_maxima = np.hstack((local_maxima, radii)) + + for i in range(0, local_maxima.shape[0]): + for j in range(0, local_maxima.shape[0]): + distance = np.sqrt( + (local_maxima[i][0] - local_maxima[j][0]) ** 2 + (local_maxima[i][1] - local_maxima[j][1]) ** 2) + if (local_maxima[j][2] > local_maxima[i][2] and distance < local_maxima[i][3]): + local_maxima[i][3] = distance + + indices = np.argsort(-local_maxima[:, 3]) + local_maxima = local_maxima[indices] + return local_maxima[:100] + + +def harris_corner_detection(I, evaluation_fraction=.05): + Su = np.matrix( + [[-1, 0, 1], + [-2, 0, 2], + [-1, 0, 1]]) + w = np.matrix([ + [0.023528, 0.033969, 0.038393, 0.033969, 0.023528], + [0.033969, 0.049045, 0.055432, 0.049045, 0.033969], + [0.038393, 0.055432, 0.062651, 0.055432, 0.038393], + [0.033969, 0.049045, 0.055432, 0.049045, 0.033969], + [0.023528, 0.033969, 0.038393, 0.033969, 0.023528] + ]) + + Iu = signal.convolve2d(I, Su) + Iv = signal.convolve2d(I, Su.T) + + Iuu = signal.convolve2d(np.multiply(Iu, Iu), w) + Ivv = signal.convolve2d(np.multiply(Iv, Iv), w) + Iuv = signal.convolve2d(np.multiply(Iu, Iv), w) + + H = np.divide(np.multiply(Iuu, Ivv) - np.multiply(Iuv, Iuv), Iuu + Ivv + 1e-10) + + local_maxima = pull_local_maxima(H, evaluation_fraction) + return local_maxima diff --git a/keypoint_detection/chessboard.png b/keypoint_detection/chessboard.png new file mode 100644 index 0000000..29e51fe Binary files /dev/null and b/keypoint_detection/chessboard.png differ diff --git a/keypoint_detection/convoluution.py b/keypoint_detection/convoluution.py new file mode 100644 index 0000000..d791de9 --- /dev/null +++ b/keypoint_detection/convoluution.py @@ -0,0 +1,90 @@ +import numpy as np +import matplotlib.pyplot as plt + +def convolve_no_edge(g,h): + gx = g.shape[0] + gy = g.shape[1] + + hx = h.shape[0] + hy = h.shape[1] + padx = hx/2 + pady = hy/2 + res = np.zeros((gx - (padx - 1), gy - (pady - 1))) + for x in range (padx, gx - padx): + for y in range(pady, gy - pady): + g_local = g[x-padx:x+padx+1,y-pady:y+pady+1] + res_pixel = np.sum(np.multiply(g_local, h)) + res[x-padx,y-pady] = res_pixel + return res + +def convolve_fill(g,h, fill_val): + gx = g.shape[0] + gy = g.shape[1] + + hx = h.shape[0] + hy = h.shape[1] + padx = hx/2 + pady = hy/2 + + res = np.zeros((gx, gy)) + for x in range (0, gx): + for y in range(0, gy): + g_local = np.array([]) + for j in range(x - padx, x + padx + 1): + for k in range(y - pady, y + pady + 1): + if (j < 0 or j >= gx or k < 0 or k >= gy): + g_local = np.append(g_local, fill_val) + else: + g_local = np.append(g_local, g[j,k]) + g_local = g_local.reshape((3,3)) + res_pixel = np.sum(np.multiply(g_local, h)) + res[x-padx,y-pady] = res_pixel + return res + +def convolve_nearest_neighbors(g, h): + gx = g.shape[0] + gy = g.shape[1] + + hx = h.shape[0] + hy = h.shape[1] + padx = hx/2 + pady = hy/2 + + res = np.zeros((gx, gy)) + for x in range (0, gx): + for y in range(0, gy): + g_local = np.array([]) + for j in range(x - padx, x + padx + 1): + for k in range(y - pady, y + pady + 1): + if (j < 0 or j >= gx or k < 0 or k >= gy): + j_copy = j + k_copy = k + while ((j_copy < 0 or j_copy >= gx) or (k_copy < 0 or k_copy >= gy)): + if (j_copy < 0): + j_copy+=1 + elif (j_copy >= gx): + j_copy-=1 + if (k_copy < 0): + k_copy+=1 + elif (k_copy >= gy): + k_copy-=1 + g_local = np.append(g_local, g[j_copy, k_copy]) + else: + g_local = np.append(g_local, g[j,k]) + g_local = g_local.reshape((3,3)) + res_pixel = np.sum(np.multiply(g_local, h)) + print(res_pixel) + res[x-padx,y-pady] = res_pixel + return res +''' +img_color = plt.imread('stuff.jpeg') +I_gray = img_color.mean(axis=2) +Su = np.matrix( +[[-1, 0, 1], +[-2, 0, 2], +[-1,0,1]]) + +filter = convolve_nearest_neighbors(I_gray, Su.T) +plt.imshow(filter, cmap='gray') +plt.show() +''' diff --git a/keypoint_detection/keypoint_detection.py b/keypoint_detection/keypoint_detection.py new file mode 100644 index 0000000..9f8e2cc --- /dev/null +++ b/keypoint_detection/keypoint_detection.py @@ -0,0 +1,59 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy import signal + +def pull_local_maxima(I): + local_maxima = [] + for i in range(1, I.shape[0]-1): + for j in range(1, I.shape[1]-1): + if (I[i+1][j] <= I[i][j] + and I[i-1][j] <= I[i][j] + and I[i][j+1] <= I[i][j] + and I[i+1][j-1] <= I[i][j]): + local_maxima.append((i, j, I[i,j])) + + local_maxima.sort(key=lambda local_maxima: local_maxima[2], reverse = True) + local_maxima = np.array(local_maxima) + local_maxima = local_maxima[:local_maxima.shape[0]/20] + radii = np.zeros((local_maxima.shape[0],1)) + float("inf") + local_maxima = np.hstack((local_maxima, radii)) + + for i in range(0, local_maxima.shape[0]): + for j in range(0, local_maxima.shape[0]): + distance = np.sqrt((local_maxima[i][0] - local_maxima[j][0])**2 + (local_maxima[i][1] - local_maxima[j][1])**2) + if (local_maxima[j][2] > local_maxima[i][2] and distance < local_maxima[i][3]): + local_maxima[i][3] = r + + indices = np.argsort(-local_maxima[:,3]) + local_maxima = local_maxima[indices] + return local_maxima[:100] + +I = plt.imread('statue_of_liberty.jpg') +I = I.mean(axis=2) +Su = np.matrix( +[[-1, 0, 1], +[-2, 0, 2], +[-1,0,1]]) + +w = np.matrix([ +[0.023528, 0.033969, 0.038393, 0.033969, 0.023528], +[0.033969, 0.049045, 0.055432, 0.049045, 0.033969], +[0.038393, 0.055432, 0.062651, 0.055432, 0.038393], +[0.033969, 0.049045, 0.055432, 0.049045, 0.033969], +[0.023528, 0.033969, 0.038393, 0.033969, 0.023528] +]) + +Iu = signal.convolve2d(I, Su) +Iv = signal.convolve2d(I, Su.T) + +Iuu = signal.convolve2d(np.multiply(Iu, Iu), w) +Ivv = signal.convolve2d(np.multiply(Iv, Iv), w) +Iuv = signal.convolve2d(np.multiply(Iu, Iv), w) + +H = np.divide(np.multiply(Iuu, Ivv) - np.multiply(Iuv, Iuv), Iuu + Ivv + 1e-10) + +local_maxima = pull_local_maxima(H) + +plt.imshow(I) +plt.scatter(local_maxima[:,1], local_maxima[:,0], c=local_maxima[:,2]) +plt.show() diff --git a/keypoint_detection/statue_of_liberty.jpg b/keypoint_detection/statue_of_liberty.jpg new file mode 100644 index 0000000..6dd9f86 Binary files /dev/null and b/keypoint_detection/statue_of_liberty.jpg differ