main
Spagnolo Gasper 2022-11-26 18:02:39 +01:00
parent 1922be45e2
commit e255c1c4ff
2 changed files with 177 additions and 40 deletions

View File

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

View File

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