Ransac for computing best fundamental matrix implemented

assignment3/solution.py

@@ -744,7 +744,6 @@ def three_g():
 
     axs[1].imshow(edge_detected_image, cmap='gray')
     axs[1].set(title='Edges')
-    
     plt.show()
 
     for i in range(hugh_transformed_circle.shape[2]):

assignment4/solution.py

@@ -11,9 +11,9 @@ import os
 ##############################################
 
 def ex1():
-	one_a()
+	#one_a()
-	one_b()
+	#one_b()
-	#one_d()
+	one_d()
 
 def one_a() -> None:
 	img = uz_image.imread_gray("data/graf/graf_a.jpg", uz_image.ImageType.float64)
@@ -182,7 +182,7 @@ def two_b() -> None:
 
 
 def ex3():
-	three_a()
+	#three_a()
 	three_b()
 
 def three_a() -> None:
@@ -266,8 +266,8 @@ def three_b() -> None:
 # ######## #
 
 def main():
-	ex1() 
+	#ex1() 
-	ex2()
+	#ex2()
 	ex3()
 
 if __name__ == '__main__':

assignment5/solution.py

@@ -142,6 +142,69 @@ def two_c() -> None:
 	print('Reprojection error for house points', uz_image.reprojection_error(F, np.array([[85, 233], [85, 233]]), np.array([[67, 219], [67, 219]])))
 	print('Reprojection error for house points', uz_image.reprojection_error(F, p1, p2))
 
+def three_a() -> None:
+	"""
+	Implement the function triangulate that accepts a set of correspondence pointsand a 
+	pair of calibration matrices as an input and returns the triangulated 3Dpoints. 
+	Test the triangulation on the ten points from the filehouse_points.txt.Visualize the result usingplt.plotorplt.scatter.
+	Also plot the index of thepoint in 3D space (useplt.text) so the results will be easier to interpret. 
+	Plot thepoints interactively, so that you can rotate the visualization.	
+	"""
+	points = np.loadtxt('./data/epipolar/house_points.txt')
+	p1 = points[:, :2]
+	p2 = points[:, 2:]
+
+	projection_matrix1 = np.loadtxt('./data/epipolar/house1_camera.txt')
+	projection_matrix2 = np.loadtxt('./data/epipolar/house2_camera.txt')
+
+	triangulated_points = uz_image.triangulate(p1, p2, projection_matrix1, projection_matrix2)
+
+	transformation_matrix = np.array([[-1, 0, 0], [0, 0, -1], [0, 1, 0]])
+	
+	triangulated_points = np.dot(transformation_matrix, triangulated_points.T).T
+	
+	fig, axs = plt.subplots(1, 3)
+
+	axs[0].imshow(uz_image.imread_gray('./data/epipolar/house1.jpg', uz_image.ImageType.float64), cmap='gray')
+	axs[0].scatter(p1[:, 0], p1[:, 1], c='r', marker='o', s=10)
+	axs[0].set_title('Left image')
+
+	axs[1].imshow(uz_image.imread_gray('./data/epipolar/house2.jpg', uz_image.ImageType.float64), cmap='gray')
+	axs[1].scatter(p2[:, 0], p2[:, 1], c='r', marker='o', s=10)
+	axs[1].set_title('Right image')
+
+	axs[2] = fig.add_subplot(111, projection='3d')
+	axs[2].scatter(triangulated_points[:, 0], triangulated_points[:, 1], triangulated_points[:, 2])
+	for i, point in enumerate(triangulated_points):
+		axs[2].text(point[0], point[1], point[2], i)
+
+	plt.show()
+
+
+def three_d():
+	left_image = uz_image.imread_gray('./data/desk/DSC02638.JPG', uz_image.ImageType.float64)
+	right_image = uz_image.imread_gray('./data/desk/DSC02639.JPG', uz_image.ImageType.float64)
+
+	keypoints_left, keypoints_right = uz_image.find_matches(left_image, right_image)
+	F, keypoints_left, keypoints_right = uz_image.ransac_fundamental(keypoints_left, keypoints_right, 1000, 1)
+
+	lines_p1, lines_p2 = uz_image.get_epipolar_lines(keypoints_left, keypoints_right)
+
+	plt.imshow(left_image, cmap='gray')
+	for line, point in zip(lines_p1, keypoints_left):
+		uz_image.draw_epiline(line, left_image.shape[0], left_image.shape[1])
+		plt.plot(point[0], point[1], 'r', marker='o', markersize=10)
+	
+	plt.show()
+	
+
+	plt.imshow(right_image, cmap='gray')
+	for line, point in zip(lines_p2, keypoints_right):
+		uz_image.draw_epiline(line, right_image.shape[0], right_image.shape[1])
+		plt.plot(point[0], point[1], 'r', marker='o', markersize=10)
+
+	plt.show()
+
 
 # ######## #
 # SOLUTION #
@@ -149,7 +212,9 @@ def two_c() -> None:
 
 def main():
 	#one_d()
-	two_c()
+	three_d()
+	#three_a()
+	#two_c()
 	#ex2()
 	#ex3()
 

assignment5/uz_framework/image.py

@@ -1463,7 +1463,7 @@ def sift(grayscale_image: npt.NDArray[np.float64],
 
 
 def get_disparity(left_image: npt.NDArray[np.float64], right_image: npt.NDArray[np.float64],\
-             window_size: int = 35, patch_size: int = 11, reduce_size: float = 0.5) -> npt.NDArray[np.float64]:
+             window_size: int = 35, patch_size: int = 11, reduce_size: float = 0.3) -> npt.NDArray[np.float64]:
 
     if(patch_size % 2 == 0 or window_size % 2 == 0):
         raise ValueError("Patch/window size must be odd")
@@ -1484,9 +1484,9 @@ def get_disparity(left_image: npt.NDArray[np.float64], right_image: npt.NDArray[
 
     def compute_ncc(a, bs):
         nccs = []
-        a = a - np.average(a)
+        a = a - np.mean(a)
         for b in bs: 
-            b = b - np.average(b)
+            b = b - np.mean(b)
             nccs.append(np.sum(a*b) / (np.sqrt(np.sum(a**2)) * np.sqrt(np.sum(b**2))))
         return np.array(nccs)
 
@@ -1514,7 +1514,7 @@ def get_disparity(left_image: npt.NDArray[np.float64], right_image: npt.NDArray[
             index = np.argmax(ncc) + 1
 
             # Compute the disparity
-            disparity[y, x] = abs(dw - index)
+            disparity[y, x] =  np.abs(dw - index)
 
     return disparity
 
@@ -1617,3 +1617,84 @@ def reprojection_error(F, p1, p2):
     repr_err  = np.sum((d_left + d_right) ) / (2 * p1.shape[0])
 
     return repr_err
+
+
+def triangulate(p1, p2, P1, P2):
+    # p1, p2: Nx2 vectors of points
+    # P1, P2: 3x4 projection matrices
+
+    # returns: 3D points
+
+    # Add a column of ones to the points, even if there is single point
+    p1 = np.hstack((p1, np.ones((p1.shape[0], 1))))
+    p2 = np.hstack((p2, np.ones((p2.shape[0], 1))))
+
+    D_points = []
+
+    def get_matrix_from_vec(a):
+        return np.array([[0, -a[2], a[1]], [a[2], 0, -a[0]], [-a[1], a[0], 0]])
+    
+    for p_left, p_right in zip(p1, p2):
+        # Build A1 and A2
+        X1 = get_matrix_from_vec(p_left)
+        X2 = get_matrix_from_vec(p_right)
+        A1 = np.dot(X1, P1)
+        A2 = np.dot(X2, P2)
+        # Construct A
+        A = np.vstack((A1, A2))        # Compute the SVD
+        U, S, V = np.linalg.svd(A)
+        # Get the last column of V
+        X = V[-1, :]
+        # Normalize
+        X = X / X[-1]
+    
+        D_points.append(X[:-1])
+
+    return np.array(D_points)
+
+def ransac_fundamental(correspondences_a: npt.NDArray[np.float64], correspondences_b: npt.NDArray[np.float64],
+            iterations: int = 5000,
+                       threshold: float = 3):
+    """
+    RANSAC algorithm for estimating homography.
+    Accepts two images and their corresponding keypoints.
+    Returns the best homography matrix and the inliers.
+    """
+    # Find the best homography
+    best_inliers = [] 
+    best_fund_matrix = None
+
+    for i in tqdm(range(iterations), desc='RANSAC'):
+        # Randomly sample 4 correspondences
+        sample = np.random.choice(correspondences_a.shape[0], 8, replace=False)
+        # correspondance_a[0] (x,y) -> correspondance_b[0] (x,y)
+        sample_a = correspondences_a[sample]
+        sample_b = correspondences_b[sample]
+        # Make it (x,y) -> (x,y) = x,y,x,y matrix
+        # Estimate homography
+        fundamential_matrix_x = fundamential_matrix(sample_a, sample_b)
+        # Compute the inliers
+        inlier_indices = []
+        # Calculate the distance between the transformed points and the actual points
+        # and count the number of inliers
+        for i, correspondence in enumerate(zip(correspondences_a, correspondences_b)):
+            p1, p2 = correspondence
+            if reprojection_error(fundamential_matrix_x, np.array([p1]), np.array([p2])) < threshold:
+                inlier_indices.append(i)
+
+        inliers_proportion = len(inlier_indices)/len(correspondences_a)
+        print(inliers_proportion)
+        # Check if we have a new best homography
+        if inliers_proportion > 0.15:
+            fundamential_matrix_x_2 = fundamential_matrix(correspondences_a[inlier_indices], correspondences_b[inlier_indices])
+            inlier_indices_2 = []
+            for i, correspondence in enumerate(zip(correspondences_a, correspondences_b)):
+                p1, p2 = correspondence
+                if reprojection_error(fundamential_matrix_x_2, np.array([p1]), np.array([p2])) < threshold:
+                    inlier_indices_2.append(i)
+
+            if len(inlier_indices_2) / len(correspondences_a) > len(best_inliers) / len(correspondences_a):
+                best_inliers = inlier_indices_2
+                best_fund_matrix = fundamential_matrix_x_2 
+
+    return best_fund_matrix, correspondences_a[best_inliers], correspondences_b[best_inliers]