From e255c1c4ff14246392f049829247d807ec888c6c Mon Sep 17 00:00:00 2001 From: Spagnolo Gasper Date: Sat, 26 Nov 2022 18:02:39 +0100 Subject: [PATCH] hehe --- assignment4/solution.py | 48 ++++++++- assignment4/uz_framework/image.py | 169 +++++++++++++++++++++++------- 2 files changed, 177 insertions(+), 40 deletions(-) diff --git a/assignment4/solution.py b/assignment4/solution.py index 2259d11..90dc64a 100644 --- a/assignment4/solution.py +++ b/assignment4/solution.py @@ -28,7 +28,7 @@ def one_a() -> None: axs[0, i].set_title(f"Sigma: {sigma}") # Plot grayscale image axs[1, i].imshow(img, cmap="gray") - # Plot scatter hessian points + # Plot scatter hessian points (x, y) axs[1, i].scatter(hessian_points[:, 1], hessian_points[:, 0], s=20, c="r", marker="x") plt.show() @@ -51,13 +51,57 @@ def one_b() -> None: axs[1, i].scatter(harris_points[:, 1], harris_points[:, 0], s=20, c="r", marker="x") plt.show() + +def ex2(): + two_a() + +def two_a() -> None: + """ + Hello + """ + graph_a_small = uz_image.imread_gray("data/graf/graf_a_small.jpg", uz_image.ImageType.float64) + graph_b_small = uz_image.imread_gray("data/graf/graf_b_small.jpg", uz_image.ImageType.float64) + + # Get the keypoints + _, graph_a_keypoints = uz_image.harris_detector(graph_a_small, 3, treshold=1e-6) + _, graph_b_keypoints = uz_image.harris_detector(graph_b_small, 3, treshold=1e-6) + + # Get the descriptors + graph_a_descriptors = uz_image.simple_descriptors(graph_a_small, graph_a_keypoints[:,0], graph_a_keypoints[:,1]) + graph_b_descriptors = uz_image.simple_descriptors(graph_b_small, graph_b_keypoints[:,0], graph_b_keypoints[:,1]) + + # Find the correspondences + matches_a = uz_image.find_correspondences(graph_a_descriptors, graph_b_descriptors) + matches_b = uz_image.find_correspondences(graph_b_descriptors, graph_a_descriptors) + + matches_a_coordinates = [] + matche_b_coordinates = [] + for i, match in enumerate(matches_a): + if i % 2 == 0: # plot every second one + if np.flip(match) in matches_b: # Check if the match is reciprocal + matches_a_coordinates.append(np.flip(graph_a_keypoints[match[0]])) + matche_b_coordinates.append(np.flip(graph_b_keypoints[match[1]])) + else: + print("Not reciprocal") + + # Plot the matches + uz_image.display_matches(graph_a_small, matches_a_coordinates, graph_b_small, matche_b_coordinates) + +def two_b() -> None: + """ + jjjjj + """ + + + # ######## # # SOLUTION # # ######## # def main(): - ex1() + #ex1() + ex2() if __name__ == '__main__': main() diff --git a/assignment4/uz_framework/image.py b/assignment4/uz_framework/image.py index 7de1e7e..0442942 100644 --- a/assignment4/uz_framework/image.py +++ b/assignment4/uz_framework/image.py @@ -132,7 +132,7 @@ def calculate_best_treshold_using_otsu_method(image: Union[npt.NDArray[np.float6 im = image.copy() else: raise Exception("Unrecognized image format!") - + treshold_range = np.arange(np.max(im) + 1) criterias = [] @@ -214,7 +214,7 @@ def get_image_bins(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]] for i in range(1, number_of_bins + 1): if i not in bins: empty_bins.append(i) - + # Add empty bins with zeros if empty_bins != []: for i in empty_bins: @@ -271,13 +271,13 @@ def get_image_bins_gradient_magnitude_and_angles(image: Union[npt.NDArray[np.flo gm, da = gradient_magnitude(image, 1 ) gm_blocks = split_into_blocks(gm) da_blocks = split_into_blocks(da) - + histogram = [] for ix in range(64): # Gradient magnitude and angles gm_l = gm_blocks[ix].reshape(-1) da_l = da_blocks[ix].reshape(-1) - + accumulator = np.zeros(8) binned = np.digitize(da_l, angle_linspace) for i in range(binned.size): @@ -309,7 +309,7 @@ def compare_two_histograms(h1: npt.NDArray[np.float64], h2: npt.NDArray[np.float raise Exception('Unsuported method chosen!') return d.astype(float) - + def apply_mask_on_image(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], mask: npt.NDArray[np.uint8]) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]: @@ -351,18 +351,18 @@ def simple_convolution_improved(signal: npt.NDArray[np.float64], kernel: npt.NDA signal_len = signal.size kernel_len = kernel.size signal = signal.copy() - + # Calculate which values to fill in and range EDGE_FRONT = signal[0] EDGE_BACK = signal[-1] EXTEND_RANGE = int(np.floor(kernel.size / 2)) - + # Append end insert edges front = [ EDGE_FRONT for _ in range(EXTEND_RANGE)] back = [ EDGE_BACK for _ in range(EXTEND_RANGE)] signal = np.insert(signal, 1, front, axis=0) signal = np.append(signal, back, axis=0) - + convolved_signal = np.zeros(signal_len) rev_kernel = kernel[::-1].copy() n_conv = signal_len - kernel_len + 1 @@ -400,8 +400,8 @@ def sharpen_image(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], sharpened_image = image.copy() KERNEL = np.array([[-1, -1, -1], - [-1, 17, -1], - [-1, -1,-1]]) * 1./9. * sharpen_factor + [-1, 17, -1], + [-1, -1,-1]]) * 1./9. * sharpen_factor if image.dtype.type == np.float64: sharpened_image = cv2.filter2D(sharpened_image, cv2.CV_64F, KERNEL) @@ -421,7 +421,7 @@ def gaussfilter2D(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], kernel = np.array([get_gaussian_kernel(sigma)]) filtered_image = cv2.filter2D(filtered_image, cv2.CV_64F, kernel) filtered_image = cv2.filter2D(filtered_image, cv2.CV_64F, kernel.T) - + return filtered_image def gaussdx(sigma: float) -> npt.NDArray[np.float64]: @@ -499,7 +499,7 @@ def filter_laplace(image:Union[npt.NDArray[np.float64], # Now apply laplacian filter applied_by_x = cv2.filter2D(image, -1, laplacian_filter) applied_by_y = cv2.filter2D(applied_by_x, -1, laplacian_filter.T) - + return applied_by_y def gauss_noise(I, magnitude=.1) -> npt.NDArray[np.float64]: @@ -570,7 +570,7 @@ def derive_image_by_x(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8 gaussd = np.array([gaussdx(sigma)]) gauss = np.array([get_gaussian_kernel(sigma)]) gaussd = np.flip(gaussd, axis=1) - + applied_by_y = cv2.filter2D(image, cv2.CV_64F, gauss.T) applied_by_x = cv2.filter2D(applied_by_y, cv2.CV_64F, gaussd) @@ -594,7 +594,7 @@ def derive_image_by_y(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8 def derive_image_first_order(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], sigma: float) -> tuple[Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], - Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]]: + Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]]: """ Accepts: image returns: image derived by x, image derived by y @@ -603,8 +603,8 @@ def derive_image_first_order(image: Union[npt.NDArray[np.float64], npt.NDArray[n def derive_image_second_order(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], sigma: float) -> tuple[tuple[Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], - Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]], tuple[Union[npt.NDArray[np.float64], - npt.NDArray[np.uint8]], Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]]]: + Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]], tuple[Union[npt.NDArray[np.float64], + npt.NDArray[np.uint8]], Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]]]: """ Accepts: image Returns: Ixx, Ixy, Iyx, Iyy @@ -616,7 +616,7 @@ def derive_image_second_order(image: Union[npt.NDArray[np.float64], npt.NDArray[ def gradient_magnitude(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], sigma: float) -> tuple[Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], - Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]]: + Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]]: """ Accepts: image Returns: gradient magnitude of image and derivative angles @@ -681,7 +681,7 @@ def find_edges_nms(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], if reduced_magnitude[y, x] < reduced_magnitude[y+gp1[0], x+gp1[1]] or reduced_magnitude[y, x] < reduced_magnitude[y+gp2[0], x+gp2[1]]: reduced_magnitude[y, x] = 0 - + return reduced_magnitude @@ -698,14 +698,14 @@ def find_edges_canny(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8] image_nms = image_nms.astype(np.uint8) new_image = np.zeros_like(image_nms) - + num_labels, labels, _, _ = cv2.connectedComponentsWithStats(image_nms, connectivity=4, ltype=cv2.CV_32S) - + for i in range (1, num_labels): idxs = np.where(image_nms[labels == i] > t_high) if idxs != []: - new_image[labels == i] = 1 + new_image[labels == i] = 1 return new_image @@ -716,7 +716,7 @@ def hough_transform_a_point(x: int, y: int, n_bins: int) -> npt.NDArray[np.float Returns: x, y transformed into hough space """ theta_values = np.linspace(-np.pi/2, np.pi, n_bins) - + cos_theta = np.cos(theta_values) sin_theta = np.sin(theta_values) accumlator = np.zeros((n_bins, n_bins)) @@ -742,7 +742,7 @@ def hough_transform_a_circle(image: Union[npt.NDArray[np.float64] , npt.NDArray[ indices = np.argwhere(edged_image) gm, ga = gradient_magnitude(image, 1) - + # Loop through all nonzero pixels above treshold for i in tqdm(range(len(indices)), desc='Hough transform'): for r in range(0, r_end - r_start): @@ -752,7 +752,7 @@ def hough_transform_a_circle(image: Union[npt.NDArray[np.float64] , npt.NDArray[ b = y - r * np.sin(ga[y, x]) accumlator[int(a), int(b), r ] += 1 return accumlator - + def hough_find_lines(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], n_bins_theta: int, n_bins_rho: int, treshold: float) -> npt.NDArray[np.uint64]: @@ -771,20 +771,20 @@ def hough_find_lines(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8] cos_precalculated = np.cos(theta_values) sin_precalculated = np.sin(theta_values) - + y_s, x_s = np.nonzero(image) - + # Loop through all nonzero pixels above treshold for i in tqdm(range(len(y_s)), desc='Hough transform'): x = x_s[i] y = y_s[i] - + # Precalculate rhos rhos = np.round(x* cos_precalculated + y* sin_precalculated).astype(np.int64) - + # Bin the rhos binned = np.digitize(rhos, rho_values) -1 - + # Add to accumulator for theta in range(n_bins_theta): accumulator[binned[theta], theta] += 1 @@ -792,8 +792,8 @@ def hough_find_lines(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8] return accumulator def hough_find_lines_i(image_with_lines: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]], gradient_angles: npt.NDArray[np.float64], - gradient_magnitude: npt.NDArray[np.float64], - n_bins_theta: int, n_bins_rho: int, treshold: float) -> Union[npt.NDArray[np.uint64], npt.NDArray[np.float64]]: + gradient_magnitude: npt.NDArray[np.float64], + n_bins_theta: int, n_bins_rho: int, treshold: float) -> Union[npt.NDArray[np.uint64], npt.NDArray[np.float64]]: """" Accepts: bw image with lines, n_bins_theta, n_bins_rho, treshold Returns: image points above treshold transformed into hough space @@ -808,7 +808,7 @@ def hough_find_lines_i(image_with_lines: Union[npt.NDArray[np.float64], npt.NDAr cos_precalculated = np.cos(theta_values) sin_precalculated = np.sin(theta_values) - + indices = np.argwhere(image) # Loop through all nonzero pixels above treshold @@ -837,7 +837,7 @@ def nonmaxima_suppression_box(image: Union[npt.NDArray[np.float64], npt.NDArray[ neighbours = image[y - 1 : y + 2, x - 1 : x + 2] if image[y, x] >= np.max(neighbours): nms_image[y, x] = image[y, x] - + return nms_image def retrieve_hough_pairs(original_image: npt.NDArray[np.float64], hough_image: npt.NDArray[np.uint64], @@ -877,14 +877,14 @@ def select_best_pairs(image_line_params, n =10): partition = np.argpartition(image_line_params, kth=len(image_line_params) - n - 1, axis=0)[-n:] image_line_params = image_line_params[partition.T[0]] return image_line_params - + def get_line_to_plot(rho, theta, h, w): """ Accepts: rho, theta, image height h, image width w Returns: line to plot usage example: plt.plot(uz_image.get_line_to_plot(rho, theta, h, w), 'r', linewidth=.7) """ - + c = np.cos(theta) s = np.sin(theta) @@ -925,7 +925,7 @@ def hessian_points(image: Union[npt.NDArray[np.float64], # Calculate determinant determinant = ixx * iyy - ixy * iyx - + # Apply nonmaxima suppression points = nonmaxima_suppression_box(determinant) @@ -965,6 +965,99 @@ def harris_detector(image: Union[npt.NDArray[np.float64], nms_features = nonmaxima_suppression_box(features) points = np.argwhere(nms_features > treshold) - + return features.astype(np.float64), points.astype(np.float64) - + +def simple_descriptors(I, Y, X, n_bins = 16, radius = 40, sigma = 2): + """ + Computes descriptors for locations given in X and Y. + + I: Image in grayscale. + Y: list of Y coordinates of locations. (Y: index of row from top to bottom) + X: list of X coordinates of locations. (X: index of column from left to right) + + Returns: tensor of shape (len(X), n_bins^2), so for each point a feature of length n_bins^2. + """ + + assert np.max(I) <= 1, "Image needs to be in range [0, 1]" + assert I.dtype == np.float64, "Image needs to be in np.float64" + # Additional assertions for dumb programmers as me + assert len(X) == len(Y), "X and Y need to have same length" + assert len(X) > 0, "X and Y need to have at least one element" + + g = get_gaussian_kernel(sigma) + d = gaussdx(sigma) + + Ix = convolve(I, g.T, d) + Iy = convolve(I, g, d.T) + Ixx = convolve(Ix, g.T, d) + Iyy = convolve(Iy, g, d.T) + + mag = np.sqrt(Ix ** 2 + Iy ** 2) + mag = np.floor(mag * ((n_bins - 1) / np.max(mag))) + + feat = Ixx + Iyy + feat += abs(np.min(feat)) + feat = np.floor(feat * ((n_bins - 1) / np.max(feat))) + + desc = [] + + for y, x in zip(Y, X): + miny = max(y - radius, 0) + maxy = min(y + radius, I.shape[0]) + minx = max(x - radius, 0) + maxx = min(x + radius, I.shape[1]) + miny, maxy, minx, maxx = int(miny), int(maxy), int(minx), int(maxx) # Convert to int for indexing + r1 = mag[miny:maxy, minx:maxx].reshape(-1) + r2 = feat[miny:maxy, minx:maxx].reshape(-1) + + a = np.zeros((n_bins, n_bins)) + for m, l in zip(r1, r2): + a[int(m), int(l)] += 1 + + a = a.reshape(-1) + a /= np.sum(a) + + desc.append(a) + + return np.array(desc) + +def find_correspondences(img_a_descriptors: npt.NDArray[np.float64], + img_b_descriptors: npt.NDArray[np.float64]): + correspondances = [] + + # Find img_a correspondences + for idx, descriptor_a in enumerate(img_a_descriptors): + dists = np.sqrt(0.5 * np.sum(np.square(np.sqrt(descriptor_a) - np.sqrt(img_b_descriptors)), axis=1)) + min_dist_idx = np.argmin(dists) + correspondances.append((idx, min_dist_idx)) + + return np.array(correspondances) + + + +def display_matches(I1, pts1, I2, pts2): + """ + Displays matches between images. + + I1, I2: Image in grayscale. + pts1, pts2: Nx2 arrays of coordinates of feature points for each image (first column is x, second is y coordinates) + """ + + assert I1.shape[0] == I2.shape[0] and I1.shape[1] == I2.shape[1], "Images need to be of the same size." + + I = np.hstack((I1, I2)) + w = I1.shape[1] + plt.imshow(I, cmap='gray') + + for p1, p2 in zip(pts1, pts2): + x1 = p1[0] + y1 = p1[1] + x2 = p2[0] + y2 = p2[1] + plt.plot(x1, y1, 'bo', markersize=3) + plt.plot(x2 + w, y2, 'bo', markersize=3) + plt.plot([x1, x2 + w], [y1, y2], 'r', linewidth=.8) + + plt.show() +