diff --git a/assignment3/solution.py b/assignment3/solution.py index bb971dc..5417319 100644 --- a/assignment3/solution.py +++ b/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]): diff --git a/assignment4/solution.py b/assignment4/solution.py index 3946fbd..f304478 100644 --- a/assignment4/solution.py +++ b/assignment4/solution.py @@ -11,9 +11,9 @@ import os ############################################## def ex1(): - one_a() - one_b() - #one_d() + #one_a() + #one_b() + 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() - ex2() + #ex1() + #ex2() ex3() if __name__ == '__main__': diff --git a/assignment5/solution.py b/assignment5/solution.py index ea9679a..4facc98 100644 --- a/assignment5/solution.py +++ b/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() diff --git a/assignment5/uz_framework/image.py b/assignment5/uz_framework/image.py index 0ff7292..bf15dbd 100644 --- a/assignment5/uz_framework/image.py +++ b/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]