A1 parts.

.gitignore

@@ -4,5 +4,5 @@
 venv/*
 env/*
 
-__pycache__/*
+*/__pycache__/*
 *.pyc

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)

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)

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]
+])

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)

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()