diff --git a/assignment5/a5_utils.py b/assignment5/a5_utils.py deleted file mode 100644 index 73224e6..0000000 --- a/assignment5/a5_utils.py +++ /dev/null @@ -1,29 +0,0 @@ -import numpy as np -import cv2 -from matplotlib import pyplot as plt - -def normalize_points(P): - # P must be a Nx2 vector of points - # first coordinate is x, second is y - - # returns: normalized points in homogeneous coordinates and 3x3 transformation matrix - - mu = np.mean(P, axis=0) # mean - scale = np.sqrt(2) / np.mean(np.sqrt(np.sum((P-mu)**2,axis=1))) # scale - T = np.array([[scale, 0, -mu[0]*scale],[0, scale, -mu[1]*scale],[0,0,1]]) # transformation matrix - P = np.hstack((P,np.ones((P.shape[0],1)))) # homogeneous coordinates - res = np.dot(T,P.T).T - return res, T - -def draw_epiline(l,h,w): - # l: line equation (vector of size 3) - # h: image height - # w: image width - - x0, y0 = map(int, [0, -l[2]/l[1]]) - x1, y1 = map(int, [w-1, -(l[2]+l[0]*w)/l[1]]) - - plt.plot([x0,x1],[y0,y1],'r') - - plt.ylim([0,h]) - plt.gca().invert_yaxis() \ No newline at end of file diff --git a/assignment5/solution.py b/assignment5/solution.py index a661209..ea9679a 100644 --- a/assignment5/solution.py +++ b/assignment5/solution.py @@ -4,14 +4,6 @@ from matplotlib import pyplot as plt import cv2 import uz_framework.image as uz_image import uz_framework.text as uz_text -import os - -############################################## -# EXCERCISE 1: Exercise 1: Image derivatives # -############################################## - -def ex1(): - one_b() def one_a() -> None: """ @@ -50,16 +42,114 @@ def one_c() -> None: in x axis in the left camera and at the pixel 300 in the right camera. How far is the object (in meters) in this case? How far is the object if the object is detected at pixel 540 in the right camera? Solve this task analytically and bring your solution - to the presentation of the exercise. + to the presentation of the exercise + + z1 = 0,16216m + z2 = 0.402m """ +def one_d() -> None: + """ + Write a script that calculates the disparity for an image pair. Use + the images in the directory disparity. Since the images were pre-processed we can + limit the search for the most similar pixel to the same row in the other image. Since + just the image intensity carries too little information, we will instead compare small + """ + + left_image = uz_image.imread_gray('/home/gasperspagnolo/Documents/faks_git/uz_assignments/assignment5/data/disparity/office_left.png', uz_image.ImageType.float64) + right_image = uz_image.imread_gray('/home/gasperspagnolo/Documents/faks_git/uz_assignments/assignment5/data/disparity/office_right.png', uz_image.ImageType.float64) + + d = uz_image.get_disparity(left_image, right_image) + + plt.imshow(d, cmap='gray') + plt.show() + + +def two_b() -> None: + """ + mplement a function fundamental_matrix that is given a set of (at least) eight + pairs of points from two images and computes the fundamental matrix using the + eight-point algorithm. + As the eight-point algorithm can be numerically unstable, it is usually not executed + directly on given pairs of points. Instead, the input is first normalized by centering + them to their centroid and scaling their positions so that the average distance to the + centroid is √2. To achieve this, you can use the function normalize_points from + the supplementary material. + Extend the function fundamental_matrix so that it first normalizes the input point- + set of the left camera (we get transformed points and the transformation matrix T1) + and then transform the input point set of the right camera (we get the transformed + points and the transformation matrix T2). Using the transformed points the algo- + rithm computes fundamental matrix ˆF, then transforms it into the original space + using both transformation matrices F = TT + 2 ˆFT1. + Test your function for fundamental matrix estimation using ten correspondence + pairs that you load from the file house_points.txt. The columns are formatted as + follows: x1, y1, x2, y2, i.e. the first column contains the x-coordinates of the points for + the first image etc. Compute the fundamental matrix F and for each point in each + image calculate the corresponding epipolar line in the other image. You can draw the + epipolar lines using draw_epiline from the supplementary material. According to + epipolar geometry the corresponding epipolar line should pass through the point. As + a testing reference the correct fundamental matrix is included in the supplementary + material in file house_fundamental.txt + """ + + points = np.loadtxt('./data/epipolar/house_points.txt') + left_image = uz_image.imread_gray('./data/epipolar/house1.jpg', uz_image.ImageType.float64) + right_image = uz_image.imread_gray('./data/epipolar/house2.jpg', uz_image.ImageType.float64) + + p1 = points[:, :2] + p2 = points[:, 2:] + + lines_p1, lines_p2 = uz_image.get_epipolar_lines(p1, p2) + + plt.imshow(left_image, cmap='gray') + for line, point in zip(lines_p1, p1): + 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, p2): + 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() + +def two_c() -> None: + """ + We use the reprojection error as a quantitative measure of the quality of the estimated fundamental matrix. + Write a function reprojection_error that calculates the reprojection error of a + fundamental matrix F given two matching points. For each point, the function + should calculate the corresponding epipolar line from the point’s match in the other + image, then calculate the perpendicular distance between the point and the line + where a, b and c are the parameters of the epipolar line. Finally, the function should + return the average of the two distances. + Write a script that performs two tests: (1) compute the reprojection error for points + p1 = [85, 233]T + in the left image-plane and p2 = [67, 219]T + in right image-plane using + the fundamental matrix (the error should be approximately 0.15 pixels). (2) Load + the points from the file house_points.txt and compute the average of symmetric + reprojection errors for all pairs of points. If your calculation is correct, the average + error should be approximately 0.33 pixels. + """ + points = np.loadtxt('./data/epipolar/house_points.txt') + p1 = points[:, :2] + p2 = points[:, 2:] + F = uz_image.fundamential_matrix(p1, p2) + 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)) + # ######## # # SOLUTION # # ######## # def main(): - ex1() + #one_d() + two_c() #ex2() #ex3() diff --git a/assignment5/uz_framework/image.py b/assignment5/uz_framework/image.py index 72d60e2..0ff7292 100644 --- a/assignment5/uz_framework/image.py +++ b/assignment5/uz_framework/image.py @@ -1460,3 +1460,160 @@ def sift(grayscale_image: npt.NDArray[np.float64], keypoints = keypoints[:, :-2] return np.array(keypoints), np.array(descriptors) + + +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]: + + if(patch_size % 2 == 0 or window_size % 2 == 0): + raise ValueError("Patch/window size must be odd") + + # Reduce the size of the images + left_image = cv2.resize(left_image, None, fx=reduce_size, fy=reduce_size) + right_image = cv2.resize(right_image, None, fx=reduce_size, fy=reduce_size) + + disparity = np.zeros(left_image.shape) + + def get_patches_per_row(y, image, patch_size): + # Get the patches for the given y value + dd = patch_size//2 + patches = [] + for x in range(dd, image.shape[1] - dd, 1): + patches.append(image[y-dd:y+dd, x-dd:x+dd+1]) + return np.array(patches) + + def compute_ncc(a, bs): + nccs = [] + a = a - np.average(a) + for b in bs: + b = b - np.average(b) + nccs.append(np.sum(a*b) / (np.sqrt(np.sum(a**2)) * np.sqrt(np.sum(b**2)))) + return np.array(nccs) + + # Get the patches + dd = patch_size//2 + dw = window_size//2 + + for y in tqdm(range(dd, left_image.shape[0] - dd, 1), desc="Computing disparity"): + # Get the patches for the right image + right_patches = get_patches_per_row(y, right_image, patch_size) + left_patches = get_patches_per_row(y, left_image, patch_size) + dd = patch_size//2 + dw = window_size//2 + xlen = len(left_patches) + + for ix in range(xlen): + # Compute the NCC + x = ix + dd + patch = left_patches[ix] + + window = right_patches[max(0, ix-dw):min(xlen, ix+dw+1)] + ncc = compute_ncc(patch, window) + + # Get the index of the maximum value + index = np.argmax(ncc) + 1 + + # Compute the disparity + disparity[y, x] = abs(dw - index) + + return disparity + +def normalize_points(P): + # P must be a Nx2 vector of points + # first coordinate is x, second is y + + # returns: normalized points in homogeneous coordinates and 3x3 transformation matrix + + mu = np.mean(P, axis=0) # mean + scale = np.sqrt(2) / np.mean(np.sqrt(np.sum((P-mu)**2,axis=1))) # scale + T = np.array([[scale, 0, -mu[0]*scale],[0, scale, -mu[1]*scale],[0,0,1]]) # transformation matrix + P = np.hstack((P,np.ones((P.shape[0],1)))) # homogeneous coordinates + res = np.dot(T,P.T).T + return res, T + +def draw_epiline(l,h,w): + # l: line equation (vector of size 3) + # h: image height + # w: image width + + x0, y0 = map(int, [0, -l[2]/l[1]]) + x1, y1 = map(int, [w-1, -(l[2]+l[0]*w)/l[1]]) + + plt.plot([x0,x1],[y0,y1],'r') + + plt.ylim([0,h]) + plt.gca().invert_yaxis() + +def fundamential_matrix(P1, P2): + # P1 and P2 are Nx2 vectors of points + + # returns: fundamental matrix + N = P1.shape[0] + # normalize points + P1n, T1 = normalize_points(P1) + P2n, T2 = normalize_points(P2) + + # compute fundamental matrix + A = np.zeros((N,9)) + for i in range(N): + u = P1n[i,0] + v = P1n[i,1] + u_ = P2n[i,0] + v_ = P2n[i,1] + A[i,:] = [u*u_, u_*v, u_, u*v_, v*v_, v_, u, v, 1] + + U,S,V = np.linalg.svd(A) + # F is the last column of V + F = V[-1,:].reshape((3,3)) + + # enforce rank 2 + U,S,V = np.linalg.svd(F) + S[-1] = 0 + F = np.dot(U,np.dot(np.diag(S),V)) + + # denormalize + F = np.dot(T2.T,np.dot(F,T1)) + + return F + + +def get_epipolar_lines(p1, p2): + F = fundamential_matrix(p1, p2) + + p1 = np.hstack((p1, np.ones((p1.shape[0], 1)))) + p2 = np.hstack((p2, np.ones((p2.shape[0], 1)))) + + lines_p2 = np.dot(F, p1.T).T + lines_p1 = np.dot(F.T, p2.T).T + + return lines_p1, lines_p2 + + +def reprojection_error(F, p1, p2): + + # 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)))) + + lines_right = [] + lines_left = [] + for p_left, p_right in zip(p1, p2): + # Left to right + lines_right.append(np.dot(F, p_left).T) + # Right to left + lines_left.append(np.dot(F.T, p_right).T) + + lines_right = np.array(lines_right) + lines_left = np.array(lines_left) + + # Normalize + lines_right = lines_right / lines_right[:, 2][:, np.newaxis] + lines_left = lines_left / lines_left[:, 2][:, np.newaxis] + + # Compute the distances + d_left = np.abs(np.sum(lines_right * p2, axis=1)) / np.sqrt(lines_right[:, 0] ** 2 + lines_right[:, 1] ** 2) + d_right = np.abs(np.sum(lines_left * p1, axis=1)) / np.sqrt(lines_left[:, 0] ** 2 + lines_left[:, 1] ** 2) + + repr_err = np.sum((d_left + d_right) ) / (2 * p1.shape[0]) + + return repr_err