diff --git a/.gitignore b/.gitignore index 76aa0f9..4de4a00 100644 --- a/.gitignore +++ b/.gitignore @@ -4,5 +4,5 @@ venv/* env/* -__pycache__/* +*/__pycache__/* *.pyc diff --git a/A1/code.py b/A1/code.py new file mode 100644 index 0000000..a660dfc --- /dev/null +++ b/A1/code.py @@ -0,0 +1,85 @@ +import cv2 +import numpy as np +import matplotlib.pyplot as plt +import pdb +from scipy import linalg + + +def find_corners(img): + """Harris corner detection to find the points in the image.""" + gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) + # find Harris corners + gray = np.float32(gray) + dist = cv2.cornerHarris(gray, 2, 3, 0.04) + dist = cv2.dilate(dist, None) + ret, dist = cv2.threshold(dist, 0.01*dist.max(), 255, 0) + dist = np.uint8(dist) + + # find centroids + ret, labels, stats, centroids = cv2.connectedComponentsWithStats(dist) + + # define the criteria to stop and refine the corners + criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.001) + corners = cv2.cornerSubPix(gray, np.float32(centroids), (5, 5), (-1, -1), criteria) + + # Now draw them + res = np.hstack((centroids, corners)) + res = np.int0(res) + # img[res[:,1],res[:,0]]=[0,0,255] + for i in range(res.shape[0]): + x, y = res[i, 3], res[i, 2] + img[x-2:x+2, y-2:y+2] = [0, 255, 0] + + plt.imshow(img) + plt.show() + pdb.set_trace() + + +def ransac_camera_calibrate(image_points, world_points): + N = image_points.shape[0] + n = 100 + min_error = 10**8 + opt_P = np.zeros((3, 4)) + opt_i = -1 + for i in range(n): + m = np.random.choice(N, 6, replace=False) + i_points = image_points[m] + w_points = world_points[m] + P = camera_calibrate(i_points, w_points) + m_ = np.array(list(set(range(N)) - set(m))) + test_i_points = image_points[m_] + test_w_points = world_points[m_] + test_w_points = np.concatenate([test_w_points, np.ones((len(m_), 1))], axis=-1) + pdb.set_trace() + predicted_i_points = np.dot(P, test_w_points.T).T + predicted_i_points = predicted_i_points/predicted_i_points[:, -1].reshape(len(m_), 1) + error = ((predicted_i_points[:, :2] - test_i_points)**2).sum()/len(m_) + if error < min_error: + min_error = error + opt_P = P + opt_i = i + + print(min_error, opt_P, opt_i) + return opt_P + + +if __name__ == "__main__": + img = cv2.imread('Assignment1_Data/measurements.jpg') + world_points = [[36, 0, 0], [0, 0, 36], [0, 36, 0], + [36, 36, 0], [36, 0, 36], [0, 0, 72]] + image_points = [[396.505, 209.674], [473.951, 246.394], [486.636, 138.904], + [402.514, 132.227], [385.822, 237.047], [465.94, 277.77]] + P = camera_calibrate(image_points, world_points) + K, R, C = decompose_projection(P) + + image_points = [] + world_points = [] + with open('points.txt') as f: + for line in f: + line = line.strip().split(',') + world_points.append([float(i.strip()) for i in line[0:3]]) + image_points.append([float(i.strip()) for i in line[3:]]) + image_points = np.array(image_points) + world_points = np.array(world_points) + ransac_camera_calibrate(image_points, world_points) + # find_corners(img) diff --git a/A1/dlt.py b/A1/dlt.py new file mode 100644 index 0000000..32a5f36 --- /dev/null +++ b/A1/dlt.py @@ -0,0 +1,75 @@ +import numpy as np +import pdb +from scipy import linalg + + +def decompose_projection(P): + """Decompose the projection matrix into K, R, C""" + M = P[:, 0:3] + MC = -1 * P[:, 3] + C = np.dot(np.linalg.inv(M), MC) + K, R = linalg.rq(M) + return K, R, C + + +def reconstruct_projection(K, R, C): + """Create the projection matrix to be multiplied in camera equation.""" + M = np.dot(K, R) + MC = -1 * np.dot(M, C) + P = np.hstack((M, MC)) + return P + + +def DLT_method(image_points, world_points): + assert image_points.shape[0] == world_points.shape[0] + + A = np.zeros((image_points.shape[0]*2, 11)) + U = np.zeros((image_points.shape[0]*2, 1)) + i = 0 + for i_points, w_points in zip(image_points, world_points): + u, v = i_points + x, y, z = w_points + A[i] = [x, y, z, 1, 0, 0, 0, 0, -u*x, -u*y, -u*z] + U[i] = -u + A[i+1] = [0, 0, 0, 0, x, y, z, 1, -v*x, -v*y, -v*z] + U[i+1] = -v + + i += 2 + + # Solving using SVD + H = np.hstack((A, U)) + u_, D, v_T = np.linalg.svd(H) + P = (v_T[-1]).reshape((-1, 1)) + divider = P[-1, 0] + P = P/divider + + P = P.reshape(3, 4) + + return P + + +def find_projection_error(P, image_points, world_points): + """Find the error in predicted points for a projection matrix.""" + total_error = 0 + for img_point, (x, y, z) in zip(image_points, world_points): + A = np.array([[x], [y], [z], [1]]) + i_pred = np.dot(P, A) + i_pred = i_pred.reshape(1, -1)[0] + i_pred = i_pred[:-1] / i_pred[-1] + error = np.linalg.norm(i_pred - img_point, ord=2) + + total_error += error + + return total_error + + +if __name__ == "__main__": + world_points = np.array([[36, 0, 0], [0, 0, 36], [0, 36, 0], + [36, 36, 0], [36, 0, 36], [0, 0, 72]]) + image_points = np.array([[396.505, 209.674], [473.951, 246.394], [486.636, 138.904], + [402.514, 132.227], [385.822, 237.047], [465.94, 277.77]]) + projection_matrix = DLT_method(image_points, world_points) + K, R, C = decompose_projection(projection_matrix) + + error = find_projection_error(projection_matrix, image_points, world_points) + print(error) diff --git a/A1/points.py b/A1/points.py new file mode 100644 index 0000000..992dac7 --- /dev/null +++ b/A1/points.py @@ -0,0 +1,79 @@ +import numpy as np +scale = 36 +world_points = scale * np.array([ + [0, 1, 0], + [0, 2, 0], + [1, 1, 0], + [1, 2, 0], + [2, 1, 0], + [2, 2, 0], + [3, 1, 0], + [3, 2, 0], + [0, 0, 1], + [0, 0, 2], + [1, 0, 1], + [1, 0, 2], + [2, 0, 1], + [2, 0, 2], + [3, 0, 1], + [3, 0, 2], + [4, 1, 0], + [4, 2, 0], + [5, 1, 0], + [5, 2, 0], + [6, 1, 0], + [6, 2, 0], + [4, 0, 1], + [4, 0, 2], + [5, 0, 1], + [5, 0, 2], + [0, 0, 3], + [0, 0, 4], + [1, 0, 3], + [1, 0, 4], + [2, 0, 3], + [2, 0, 4], + [3, 0, 3], + [3, 0, 4], + [4, 0, 3], + [4, 0, 4] +]) + +image_points = np.array([ + [4868.3, 1389.7], + [4916.2, 560.2], + [4030.0, 1333.6], + [4058.9, 514.9], + [3221.3, 1275.5], + [3237.1, 467.6], + [2437.6, 1218.4], + [2434.0, 423.5], + [4726.3, 2456.3], + [4638.2, 2787.5], + [3864.2, 2378.0], + [3723.8, 2711.5], + [3030.9, 2306.0], + [2841.3, 2632.3], + [2222.2, 2233.3], + [1985.3, 2546.7], + [1672.7, 1162], + [1654.0, 374.7], + [927, 1105.5], + [892, 325.5], + [194.2, 1048.8], + [135.3, 282.0], + [1432.5, 2159.0], + [1148.3, 2461.7], + [659.5, 2100.8], + [329.7, 2397.7], + [4539.0, 3161.0], + [4425.7, 3600], + [3564.0, 3076], + [3384.5, 3492.5], + [2625.3, 2992.3], + [2378.7, 3398.7], + [1713.3, 2908.7], + [1409.0, 3308.7], + [831.0, 2813.5], + [466.0, 3217.3] +]) diff --git a/A1/radial.py b/A1/radial.py new file mode 100644 index 0000000..e69de29 diff --git a/A1/ransac.py b/A1/ransac.py new file mode 100644 index 0000000..08a7c8b --- /dev/null +++ b/A1/ransac.py @@ -0,0 +1,33 @@ +import numpy as np +# import pdb + +from dlt import DLT_method, find_projection_error + + +def ransac_estimation(image_points, world_points, iter=500): + """Use the RANSAC algorithm with DLT to find best projection matrix""" + n_points = image_points.shape[0] + best_P = None + min_error = np.inf + for i in range(iter): + print("Iteration number {}".format(i)) + indices = np.random.permutation(n_points) + X, Y = image_points[indices[0:6]], world_points[indices[0:6]] + X_test, Y_test = image_points[indices[6::]], world_points[indices[6::]] + + projection_matrix = DLT_method(X, Y) + error = find_projection_error(projection_matrix, X_test, Y_test) + + if error < min_error: + min_error = error + best_P = projection_matrix + + return best_P + + +if __name__ == "__main__": + from points import image_points, world_points + projection_matrix = ransac_estimation(image_points, world_points) + + error = find_projection_error(projection_matrix, image_points, world_points) + print(error) diff --git a/A1/wireframe.py b/A1/wireframe.py new file mode 100644 index 0000000..d0a5b2a --- /dev/null +++ b/A1/wireframe.py @@ -0,0 +1,55 @@ +import cv2 +import numpy as np +import matplotlib.pyplot as plt +import pdb + +from dlt import DLT_method +from ransac import ransac_estimation + + +def predict_points(projection_matrix, world_points): + """Use the projection matrix to predict image_points""" + points = [] + for (x, y, z) in world_points: + A = np.array([[x], [y], [z], [1]]) + i_pred = np.dot(projection_matrix, A) + i_pred = i_pred.reshape(1, -1)[0] + i_pred = i_pred[:-1] / i_pred[-1] + points.append(i_pred) + + points = np.array(points, dtype=np.int) + return points + + +def plot_wireframe(image, points): + img = image.copy() + for (x, y) in points: + img = cv2.circle(img, (x, y), 50, (255, 0, 0), -1) + + for p in points: + distance = np.sum((points - p)**2, axis=1) + min_index = np.argsort(distance)[1:5] + for nearest_point in points[min_index]: + img = cv2.line(img, tuple(p), tuple(nearest_point), (255, 0, 0), 20) + return img + + +if __name__ == "__main__": + from points import image_points, world_points + + image = cv2.imread('Assignment1_Data/IMG_5455.JPG') + projection_matrix = DLT_method(image_points, world_points) + + pred_points = predict_points(projection_matrix[0:6], world_points[0:6]) + img = plot_wireframe(image, pred_points) + + plt.imshow(img) + plt.show() + + projection_matrix = ransac_estimation(image_points, world_points) + + pred_points = predict_points(projection_matrix, world_points) + img = plot_wireframe(image, pred_points) + + plt.imshow(img) + plt.show()