hvala kurcu za neumna navodila
Spagnolo Gasper
parent
a138b0ffe5
commit
72a58075a3
|
|
@ -0,0 +1,108 @@
|
|||
import numpy as np
|
||||
import cv2
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
|
||||
def gauss(sigma):
|
||||
x = np.arange(np.floor(-3 * sigma), np.ceil(3 * sigma + 1))
|
||||
k = np.exp(-(x ** 2) / (2 * sigma ** 2))
|
||||
k = k / np.sum(k)
|
||||
return np.expand_dims(k, 0)
|
||||
|
||||
|
||||
def gaussdx(sigma):
|
||||
x = np.arange(np.floor(-3 * sigma), np.ceil(3 * sigma + 1))
|
||||
k = -x * np.exp(-(x ** 2) / (2 * sigma ** 2))
|
||||
k /= np.sum(np.abs(k))
|
||||
return np.expand_dims(k, 0)
|
||||
|
||||
|
||||
def convolve(I: np.ndarray, *ks):
|
||||
"""
|
||||
Convolves input image I with all given kernels.
|
||||
|
||||
:param I: Image, should be of type float64 and scaled from 0 to 1.
|
||||
:param ks: 2D Kernels
|
||||
:return: Image convolved with all kernels.
|
||||
"""
|
||||
for k in ks:
|
||||
k = np.flip(k) # filter2D performs correlation, so flipping is necessary
|
||||
I = cv2.filter2D(I, cv2.CV_64F, k)
|
||||
return I
|
||||
|
||||
|
||||
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"
|
||||
|
||||
g = gauss(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])
|
||||
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 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 columnt 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()
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
251 275 342 270
|
||||
539 309 573 314
|
||||
140 517 239 527
|
||||
469 454 524 446
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 256 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 41 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 262 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 55 KiB |
|
|
@ -0,0 +1,3 @@
|
|||
0.76807 -0.63756 108.59988
|
||||
0.64246 0.74179 -33.09045
|
||||
0.00003 -0.00008 1.00000
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
21 96 64 52
|
||||
246 94 238 195
|
||||
25 185 10 122
|
||||
186 207 121 243
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 20 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 20 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 70 KiB |
Binary file not shown.
|
|
@ -0,0 +1,45 @@
|
|||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
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_a()
|
||||
|
||||
def one_a() -> None:
|
||||
img = uz_image.imread_gray("data/graf/graf_a.jpg", uz_image.ImageType.float64)
|
||||
# Get the hessian points
|
||||
sigmas = [3, 6, 9, 12]
|
||||
determinant, hessian_points = uz_image.hessian_points(img, 9, 0.004)
|
||||
# Plot the points
|
||||
fig, axs = plt.subplots(2, len(sigmas))
|
||||
|
||||
for i, sigma in enumerate(sigmas):
|
||||
determinant, hessian_points = uz_image.hessian_points(img, sigma, 0.004)
|
||||
# Plot determinant
|
||||
axs[0, i].imshow(determinant)
|
||||
axs[0, i].set_title(f"Sigma: {sigma}")
|
||||
# Plot grayscale image
|
||||
axs[1, i].imshow(img, cmap="gray")
|
||||
# Plot scatter hessian points
|
||||
axs[1, i].scatter(hessian_points[:, 1], hessian_points[:, 0], s=1, c="r", marker="x")
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
# ######## #
|
||||
# SOLUTION #
|
||||
# ######## #
|
||||
|
||||
def main():
|
||||
ex1()
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
|
|
@ -0,0 +1,944 @@
|
|||
|
||||
import numpy as np
|
||||
import cv2 as cv2
|
||||
from matplotlib import pyplot as plt
|
||||
from PIL import Image
|
||||
from typing import Union
|
||||
import numpy.typing as npt
|
||||
import enum
|
||||
from tqdm import tqdm
|
||||
|
||||
class ImageType(enum.Enum):
|
||||
uint8 = 0
|
||||
float64 = 1
|
||||
|
||||
class DistanceMeasure(enum.Enum):
|
||||
euclidian_distance = 0
|
||||
chi_square_distance = 1
|
||||
intersection_distance = 2
|
||||
hellinger_distance = 3
|
||||
|
||||
def imread(path: str, type: ImageType) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Reads an image in RGB order. Image type is transformed from uint8 to float, and
|
||||
range of values is reduced from [0, 255] to [0, 1].
|
||||
"""
|
||||
I = Image.open(path).convert('RGB') # PIL image.
|
||||
I = np.asarray(I) # Converting to Numpy array.
|
||||
if type == ImageType.float64:
|
||||
I = I.astype(np.float64) / 255
|
||||
return I
|
||||
elif type == ImageType.uint8:
|
||||
return I
|
||||
|
||||
raise Exception(f"Unrecognized image format! {type}")
|
||||
|
||||
|
||||
def imread_gray(path: str, type: ImageType) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Reads an image in gray. Image type is transformed from uint8 to float, and
|
||||
range of values is reduced from [0, 255] to [0, 1].
|
||||
"""
|
||||
|
||||
I = Image.open(path).convert('L') # PIL image opening and converting to gray.
|
||||
I = np.asarray(I) # Converting to Numpy array.
|
||||
|
||||
if type == ImageType.float64:
|
||||
I = I.astype(np.float64) / 255
|
||||
return I
|
||||
elif type == ImageType.uint8:
|
||||
return I
|
||||
|
||||
raise Exception("Unrecognized image format!")
|
||||
|
||||
def signal_show(*signals):
|
||||
"""
|
||||
Plots all given 1D signals in the same plot.
|
||||
Signals can be Python lists or 1D numpy array.
|
||||
"""
|
||||
for s in signals:
|
||||
if type(s) == np.ndarray:
|
||||
s = s.squeeze()
|
||||
plt.plot(s)
|
||||
plt.show()
|
||||
|
||||
|
||||
def convolve(I: np.ndarray, *ks):
|
||||
"""
|
||||
Convolves input image I with all given kernels.
|
||||
|
||||
:param I: Image, should be of type float64 and scaled from 0 to 1.
|
||||
:param ks: 2D Kernels
|
||||
:return: Image convolved with all kernels.
|
||||
"""
|
||||
for k in ks:
|
||||
k = np.flip(k) # filter2D performs correlation, so flipping is necessary
|
||||
I = cv2.filter2D(I, cv2.CV_64F, k)
|
||||
return I
|
||||
|
||||
def transform_coloured_image_to_grayscale(image: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts float64 picture with three colour channels and returns float64 grayscale image
|
||||
with one channel.
|
||||
"""
|
||||
grayscale_image = np.zeros(image.shape[:2])
|
||||
|
||||
for i in range(image.shape[0]):
|
||||
for j in range(image.shape[1]):
|
||||
grayscale_image[i, j] = (image[i, j, 0] + image[i,j, 1] + image[i, j, 2]) / 3
|
||||
|
||||
|
||||
return grayscale_image
|
||||
|
||||
def invert_coloured_image_part(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
startx: int, endx: int, starty: int, endy: int) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts image, starting position end end position for axes x & y. Returns whole image with inverted part.
|
||||
"""
|
||||
inverted_image = image.copy()
|
||||
|
||||
if image.dtype.type == np.float64:
|
||||
for i in range(startx, endx):
|
||||
for j in range(starty, endy):
|
||||
inverted_image[i, j, 0] = 1 - image[i, j, 0]
|
||||
inverted_image[i, j, 1] = 1 - image[i, j, 1]
|
||||
inverted_image[i, j, 2] = 1 - image[i, j, 2]
|
||||
return inverted_image
|
||||
elif image.dtype.type == np.uint8:
|
||||
for i in range(startx, endx):
|
||||
for j in range(starty, endy):
|
||||
inverted_image[i, j, 0] = 255 - image[i, j, 0]
|
||||
inverted_image[i, j, 1] = 255 - image[i, j, 1]
|
||||
inverted_image[i, j, 2] = 255 - image[i, j, 2]
|
||||
return inverted_image
|
||||
|
||||
raise Exception("Unrecognized image format!")
|
||||
|
||||
def invert_coloured_image(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts image and inverts it
|
||||
"""
|
||||
return invert_coloured_image_part(image, 0, image.shape[0], 0, image.shape[1])
|
||||
|
||||
def calculate_best_treshold_using_otsu_method(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]) -> int:
|
||||
"""
|
||||
Accepts image and returns best treshold using otsu method
|
||||
"""
|
||||
|
||||
if image.dtype.type == np.float64:
|
||||
im = image.copy()
|
||||
im = im * (255.0/im.max())
|
||||
elif image.dtype.type == np.uint8:
|
||||
im = image.copy()
|
||||
else:
|
||||
raise Exception("Unrecognized image format!")
|
||||
|
||||
treshold_range = np.arange(np.max(im) + 1)
|
||||
criterias = []
|
||||
|
||||
for treshold in treshold_range:
|
||||
# create the thresholded image
|
||||
thresholded_im = np.zeros(im.shape)
|
||||
|
||||
thresholded_im[im >= treshold] = 1
|
||||
|
||||
# compute weights
|
||||
nb_pixels = im.size
|
||||
nb_pixels1 = np.count_nonzero(thresholded_im)
|
||||
weight1 = nb_pixels1 / nb_pixels
|
||||
weight0 = 1 - weight1
|
||||
|
||||
# if one the classes is empty, eg all pixels are below or above the threshold, that threshold will not be considered
|
||||
# in the search for the best threshold
|
||||
if weight1 == 0 or weight0 == 0:
|
||||
continue
|
||||
|
||||
# find all pixels belonging to each class
|
||||
val_pixels1 = im[thresholded_im == 1]
|
||||
val_pixels0 = im[thresholded_im == 0]
|
||||
|
||||
# compute variance of these classes
|
||||
var0 = np.var(val_pixels0) if len(val_pixels0) > 0 else 0
|
||||
var1 = np.var(val_pixels1) if len(val_pixels1) > 0 else 0
|
||||
|
||||
criterias.append( weight0 * var0 + weight1 * var1)
|
||||
|
||||
best_threshold = treshold_range[np.argmin(criterias)]
|
||||
return best_threshold
|
||||
|
||||
|
||||
def get_image_bins_for_loop(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
number_of_bins: int) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts image in the float64 format or uint8, returns normailzed
|
||||
image bins, histogram
|
||||
"""
|
||||
if image.dtype.type == np.float64 or image.dtype.type == np.uint8:
|
||||
bin_restrictions = np.linspace(np.min(image), np.max(image), num=number_of_bins)
|
||||
else:
|
||||
raise Exception("Unrecognized image format!")
|
||||
|
||||
bins = np.zeros(number_of_bins).astype(np.float64)
|
||||
|
||||
for pixel in image.reshape(-1):
|
||||
# https://stackoverflow.com/a/16244044
|
||||
bins[np.argmax(bin_restrictions > pixel)] += 1
|
||||
|
||||
return bins / np.sum(bins)
|
||||
|
||||
|
||||
# Much faster implementation than for loop
|
||||
def get_image_bins(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
number_of_bins: int) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts image in the float64 format or uint8, returns normailzed
|
||||
image bins, histogram
|
||||
"""
|
||||
if image.dtype.type == np.float64 or image.dtype.type == np.uint8:
|
||||
bins = np.linspace(np.min(image), np.max(image), num=number_of_bins)
|
||||
else:
|
||||
raise Exception("Unrecognized image format!")
|
||||
|
||||
# Put pixels into classes
|
||||
# ex. binsize = 10 then 0.4 would map into 4
|
||||
binarray = np.digitize(image.reshape(-1), bins).astype(np.uint8)
|
||||
|
||||
# Now count those values
|
||||
binarray = np.unique(binarray, return_counts=True)
|
||||
|
||||
counts = binarray[1].astype(np.float64) # Get the counts out of tuple
|
||||
|
||||
# Check if there is any empty bin
|
||||
empty_bins = []
|
||||
bins = binarray[0]
|
||||
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:
|
||||
counts = np.insert(counts, i - 1, 0)
|
||||
|
||||
return counts
|
||||
|
||||
def get_image_bins_ND(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
number_of_bins: int) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts image in the float64 format or uint8 and number of bins
|
||||
Returns normailzed image histogram bins
|
||||
"""
|
||||
bs = []
|
||||
hist = np.zeros((number_of_bins, number_of_bins, number_of_bins))
|
||||
bins = np.linspace(np.min(image), np.max(image), num=number_of_bins)
|
||||
|
||||
for i in range(image.shape[2]):
|
||||
v = image[:, :, i].reshape(-1)
|
||||
bs.append(np.digitize(v, bins).astype(np.uint32))
|
||||
|
||||
for i in range(len(bs[0])):
|
||||
hist[bs[2][i] -1, bs[1][i] -1, bs[0][i] - 1] += 1
|
||||
|
||||
return hist / np.sum(hist)
|
||||
|
||||
def get_image_bins_gradient_magnitude_and_angles(image: Union[npt.NDArray[np.float64],
|
||||
npt.NDArray[np.uint8]]) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: image,
|
||||
Returns: 1D histogram of image using gradient magnitude and gradient angles
|
||||
Works OK on many dimensions
|
||||
"""
|
||||
WIDTH = image.shape[0]
|
||||
HEIGHT = image.shape[1]
|
||||
WIDTH_8 = WIDTH // 8
|
||||
HEIGHT_8 = HEIGHT // 8
|
||||
|
||||
def split_into_blocks(image):
|
||||
blocks = []
|
||||
for y in range(0, 7):
|
||||
for x in range(0, 8):
|
||||
if x == 7:
|
||||
blocks.append(image[y*WIDTH_8 : (y +1)* WIDTH_8, 7*HEIGHT_8 : HEIGHT].copy()) # handle the last column
|
||||
else:
|
||||
blocks.append(image[y*WIDTH_8 : (y +1)* WIDTH_8, x*HEIGHT_8 : (x +1)* HEIGHT_8].copy())
|
||||
|
||||
for x in range(0, 7):
|
||||
blocks.append(image[7*WIDTH_8 : WIDTH, x*HEIGHT_8 : (x +1)* HEIGHT_8].copy()) # handle the last row
|
||||
blocks.append(image[7*WIDTH_8 : WIDTH, 7*HEIGHT_8 : HEIGHT].copy())
|
||||
return blocks
|
||||
|
||||
angle_linspace = np.linspace(-np.pi, np.pi, num=8)
|
||||
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):
|
||||
accumulator[binned[i] - 1] += gm_l[i]
|
||||
|
||||
# Normalize accumulator
|
||||
histogram.extend(accumulator / np.sum(accumulator))
|
||||
|
||||
histogram = np.array(histogram)
|
||||
return histogram / np.sum(histogram)
|
||||
|
||||
|
||||
def compare_two_histograms(h1: npt.NDArray[np.float64], h2: npt.NDArray[np.float64],
|
||||
method: DistanceMeasure) -> float:
|
||||
"""
|
||||
Accepts two histograms and method of comparison
|
||||
Returns distance between them
|
||||
"""
|
||||
|
||||
if method == DistanceMeasure.euclidian_distance:
|
||||
d = np.sqrt(np.sum(np.square(h1 - h2)))
|
||||
elif method == DistanceMeasure.chi_square_distance:
|
||||
d = 0.5 * np.sum(np.square(h1 - h2) / (h1 + h2 + np.finfo(float).eps))
|
||||
elif method == DistanceMeasure.intersection_distance:
|
||||
d = 1 - np.sum(np.minimum(h1, h2))
|
||||
elif method == DistanceMeasure.hellinger_distance:
|
||||
d = np.sqrt(0.5 * np.sum(np.square(np.sqrt(h1) - np.sqrt(h2))))
|
||||
else:
|
||||
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]]:
|
||||
"""
|
||||
Accepts image and applys mask to image
|
||||
"""
|
||||
image = image.copy()
|
||||
|
||||
mask = np.expand_dims(mask, axis=2)
|
||||
image = mask * image
|
||||
|
||||
return image
|
||||
|
||||
def simple_convolution(signal: npt.NDArray[np.float64], kernel: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: signal & kernel
|
||||
|
||||
Returns: convolved signal with a kernel
|
||||
"""
|
||||
N = int(np.ceil(len(kernel) / 2 - 1))
|
||||
n_conv = signal.size - kernel.size + 1
|
||||
print(N, signal.size -N, n_conv)
|
||||
|
||||
convolved_signal = np.zeros(len(signal))
|
||||
rev_kernel = kernel[::-1].copy()
|
||||
|
||||
for i in range(n_conv):
|
||||
convolved_signal[i] = np.dot(signal[i: i+kernel.size], rev_kernel) # Well if you would add i+N then you wuold shift this
|
||||
|
||||
return convolved_signal
|
||||
|
||||
def simple_convolution_improved(signal: npt.NDArray[np.float64], kernel: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: signal & kernel
|
||||
|
||||
Returns: convolved signal with a kernel
|
||||
Improved method replicates edges of an signal
|
||||
"""
|
||||
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
|
||||
|
||||
for i in range(n_conv):
|
||||
convolved_signal[i] = np.dot(signal[i: i+kernel.size], rev_kernel)
|
||||
|
||||
return convolved_signal
|
||||
|
||||
def get_gaussian_kernel(sigma: float) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts sigma
|
||||
Returns gaussian kernel
|
||||
"""
|
||||
# https://github.com/mikepound/convolve/blob/master/run.gaussian.py
|
||||
kernel_size = int(2 * np.ceil(3*sigma) + 1)
|
||||
k_min_max = np.ceil(3*sigma)
|
||||
k_interval = np.arange(-k_min_max, k_min_max + 1.)
|
||||
|
||||
result = (1. / (np.sqrt(2. * np.pi )* sigma)) * np.exp(- (np.square(k_interval)) / (2.*np.square(sigma)))
|
||||
|
||||
assert(kernel_size == len(result))
|
||||
|
||||
return result / np.sum(result)
|
||||
|
||||
def sharpen_image(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
sharpen_factor=1.0) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts: image & sharpen factor
|
||||
|
||||
Returns: sharpened image
|
||||
https://blog.demofox.org/2022/02/26/image-sharpening-convolution-kernels/ <-- sharpening kernel, but also on slides
|
||||
good explanation
|
||||
"""
|
||||
sharpened_image = image.copy()
|
||||
|
||||
KERNEL = np.array([[-1, -1, -1],
|
||||
[-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)
|
||||
elif image.dtype.type == np.uint8:
|
||||
sharpened_image = cv2.filter2D(sharpened_image, cv2.CV_8U, KERNEL)
|
||||
|
||||
return sharpened_image
|
||||
|
||||
def gaussfilter2D(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
sigma: float) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts: image, sigma
|
||||
Applies gaussian noise on image
|
||||
returns: filtered_image
|
||||
"""
|
||||
filtered_image = image.copy()
|
||||
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]:
|
||||
"""
|
||||
Accepts sigma
|
||||
Returns gaussian kernel
|
||||
"""
|
||||
kernel_size = int(2 * np.ceil(3*sigma) + 1)
|
||||
k_min_max = np.ceil(3*sigma)
|
||||
k_interval = np.arange(-k_min_max, k_min_max + 1.)
|
||||
|
||||
result = - (1. / (np.sqrt(2. * np.pi )* np.power(sigma, 3))) * k_interval* np.exp(- (np.square(k_interval)) / (2.*np.square(sigma)))
|
||||
|
||||
assert(kernel_size == len(result))
|
||||
|
||||
return result / np.sum(np.abs(result))
|
||||
|
||||
|
||||
def simple_median(signal: npt.NDArray[np.float64], width: int) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: signal & width
|
||||
returns signal improved using median filter
|
||||
"""
|
||||
if width % 2 == 0:
|
||||
raise Exception('No u won\'t do that')
|
||||
|
||||
signal = signal.copy()
|
||||
for i in range(len(signal) - int(np.ceil(width/2))):
|
||||
middle_element = int(i + np.floor(width/2))
|
||||
signal[middle_element] = np.median(signal[i:i+width])
|
||||
return signal
|
||||
|
||||
def apply_median_method_2D(image:Union[npt.NDArray[np.float64],
|
||||
npt.NDArray[np.uint8]], width: int) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts: image & filter width
|
||||
returns: image with median filter applied
|
||||
"""
|
||||
if width % 2 == 0:
|
||||
raise Exception('No u won\'t do that')
|
||||
|
||||
image = image.copy()
|
||||
W_HALF = int(np.floor(width/2))
|
||||
padded_image = np.pad(image, W_HALF, mode='edge')
|
||||
IMAGE_HEIGHT = image.shape[0] # y
|
||||
IMAGE_WIDTH = image.shape[1] # x
|
||||
|
||||
for x in range(W_HALF, IMAGE_WIDTH): # I think we can start from 0, cuz we padded an image
|
||||
for y in range(W_HALF, IMAGE_HEIGHT):
|
||||
median_filter = np.zeros(0)
|
||||
STARTX = x - W_HALF
|
||||
STARTY = y - W_HALF
|
||||
for m in range(width):
|
||||
median_filter = np.append(median_filter, padded_image[STARTY + m][STARTX: STARTX + width], axis=0)
|
||||
if image.dtype.type == np.uint8:
|
||||
image[y][x] = int(np.mean(median_filter))
|
||||
else:
|
||||
image[y][x] = np.mean(median_filter)
|
||||
return image
|
||||
|
||||
def filter_laplace(image:Union[npt.NDArray[np.float64],
|
||||
npt.NDArray[np.uint8]], sigma: float) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts: image & sigma
|
||||
returns: image with laplace filter applied
|
||||
"""
|
||||
# Prepare unit impulse and gauss kernel
|
||||
unit_impulse = np.zeros((1, 2 * int(np.ceil(3*sigma)) + 1))
|
||||
unit_impulse[0][int(np.ceil(unit_impulse.size /2)) - 1]= 1
|
||||
gauss_kernel = np.array([get_gaussian_kernel(sigma)])
|
||||
assert(len(gauss_kernel[0]) == len(unit_impulse[0]))
|
||||
|
||||
laplacian_filter = unit_impulse - gauss_kernel[0]
|
||||
|
||||
# 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]:
|
||||
"""
|
||||
Accepts: image & magnitude
|
||||
Returns: image with gaussian noise applied
|
||||
"""
|
||||
# input: image, magnitude of noise
|
||||
# output: modified image
|
||||
I = I.copy()
|
||||
|
||||
return I + np.random.normal(size=I.shape) * magnitude
|
||||
|
||||
|
||||
def sp_noise(I, percent=.1) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: image & percent
|
||||
Returns: image with salt and pepper noise applied
|
||||
"""
|
||||
# input: image, percent of corrupted pixels
|
||||
# output: modified image
|
||||
|
||||
res = I.copy()
|
||||
res[np.random.rand(I.shape[0], I.shape[1]) < percent / 2] = 1
|
||||
res[np.random.rand(I.shape[0], I.shape[1]) < percent / 2] = 0
|
||||
|
||||
return res
|
||||
|
||||
def sp_noise1D(signal, percent=.1) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: signal & percent
|
||||
Returns: signal with salt and pepper noise applied
|
||||
"""
|
||||
signal = signal.copy()
|
||||
signal[np.random.rand(signal.shape[0]) < percent / 2] = 2
|
||||
signal[np.random.rand(signal.shape[0]) < percent / 2] = 1
|
||||
signal[np.random.rand(signal.shape[0]) < percent / 2] = 4
|
||||
signal[np.random.rand(signal.shape[0]) < percent / 2] = 0.4
|
||||
return signal
|
||||
|
||||
def sum_two_grayscale_images(image_a: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
image_b :Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts: image_a, image_b
|
||||
Returns: image_a + image_b
|
||||
"""
|
||||
# Merge image_a and image_b
|
||||
return (image_a + image_b)/ 2
|
||||
|
||||
def generate_dirac_impulse(size: int) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: size
|
||||
Returns: dirac impulse of size
|
||||
"""
|
||||
|
||||
dirac_impulse = np.zeros((size, size))
|
||||
dirac_impulse[int(size/2), int(size/2)] = 1
|
||||
|
||||
return dirac_impulse
|
||||
|
||||
def derive_image_by_x(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
sigma: float) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts: image
|
||||
Returns: image derived by x
|
||||
"""
|
||||
image = image.copy()
|
||||
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)
|
||||
|
||||
return applied_by_x
|
||||
|
||||
def derive_image_by_y(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
sigma: float) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts: image
|
||||
Returns: image derived by y
|
||||
"""
|
||||
image = image.copy()
|
||||
gaussd = np.array([gaussdx(sigma)])
|
||||
gauss = np.array([get_gaussian_kernel(sigma)])
|
||||
gaussd = np.flip(gaussd, axis=1)
|
||||
|
||||
applied_by_x = cv2.filter2D(image, cv2.CV_64F, gauss)
|
||||
applied_by_y = cv2.filter2D(applied_by_x, cv2.CV_64F, gaussd.T)
|
||||
|
||||
return applied_by_y
|
||||
|
||||
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]]]:
|
||||
"""
|
||||
Accepts: image
|
||||
returns: image derived by x, image derived by y
|
||||
"""
|
||||
return derive_image_by_x(image, sigma), derive_image_by_y(image, sigma)
|
||||
|
||||
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]]]]:
|
||||
"""
|
||||
Accepts: image
|
||||
Returns: Ixx, Ixy, Iyx, Iyy
|
||||
"""
|
||||
derived_by_x = derive_image_by_x(image, sigma)
|
||||
derived_by_y = derive_image_by_y(image, sigma)
|
||||
|
||||
return derive_image_first_order(derived_by_x, sigma), derive_image_first_order(derived_by_y, sigma)
|
||||
|
||||
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]]]:
|
||||
"""
|
||||
Accepts: image
|
||||
Returns: gradient magnitude of image and derivative angles
|
||||
"""
|
||||
Ix, Iy = derive_image_first_order(image, sigma)
|
||||
return np.sqrt(Ix**2 + Iy**2), np.arctan2(Iy, Ix)
|
||||
|
||||
|
||||
def find_edges_primitive(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
sigma: float, theta: float) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Aceppts: image, sigma & theta
|
||||
Returns: image with edges
|
||||
"""
|
||||
derivative_magnitude, _ = gradient_magnitude(image, sigma)
|
||||
|
||||
binary_mask = np.zeros_like(derivative_magnitude)
|
||||
binary_mask[(derivative_magnitude >= theta)] = 1
|
||||
|
||||
return binary_mask
|
||||
|
||||
|
||||
def find_edges_nms(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
sigma: float, theta: float) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Aceppts: image, sigma & theta
|
||||
Returns: image with edges
|
||||
"""
|
||||
step_size = np.pi/8
|
||||
|
||||
|
||||
def get_gradient_orientation(angle: float) -> tuple[tuple[int, int], tuple[int, int]]:
|
||||
"""
|
||||
Accepts: angle
|
||||
Returns: indexes of gradient orientation (x, y), (x, y)
|
||||
Basically walks around the unit circle and returns the indexes of the closest angle
|
||||
"""
|
||||
angle = angle % np.pi
|
||||
|
||||
for i in range(0, 8):
|
||||
if angle >= i * step_size and angle <= (i+1) * step_size:
|
||||
if i == 0 or i == 7:
|
||||
return (0, 1), (0, -1)
|
||||
elif i == 1 or i == 2:
|
||||
return (-1, -1), (1, 1)
|
||||
elif i == 3 or i == 4:
|
||||
return (-1, 0), (1, 0)
|
||||
elif i == 5 or i == 6:
|
||||
return (-1, 1), (1, -1)
|
||||
raise ValueError(f"Angle {angle} is not in range")
|
||||
|
||||
|
||||
derivative_magnitude, derivative_angles = gradient_magnitude(image, sigma)
|
||||
reduced_magnitude = derivative_magnitude.copy()
|
||||
# Set low elements to 0
|
||||
reduced_magnitude[(derivative_magnitude <= theta)] = 0
|
||||
|
||||
|
||||
for y in range(1, reduced_magnitude.shape[0]-1):
|
||||
for x in range(1, reduced_magnitude.shape[1]-1):
|
||||
gp1, gp2 = get_gradient_orientation(derivative_angles[y, x])
|
||||
|
||||
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
|
||||
|
||||
|
||||
def find_edges_canny(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
|
||||
sigma: float, theta: float, t_low: float, t_high: float) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Aceppts: image, sigma, theta, t_low, t_high
|
||||
Returns: image with edges
|
||||
"""
|
||||
|
||||
image_nms = find_edges_nms(image, sigma, theta)
|
||||
image_nms[image_nms < t_low] = 0
|
||||
image_nms[image_nms >= t_low] = 1
|
||||
|
||||
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
|
||||
|
||||
return new_image
|
||||
|
||||
|
||||
def hough_transform_a_point(x: int, y: int, n_bins: int) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: x, y, n_bins
|
||||
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))
|
||||
|
||||
for i in range(n_bins):
|
||||
r = np.round(x * cos_theta[i] + y * sin_theta[i]) + n_bins/2
|
||||
accumlator[int(r), i] += 1
|
||||
|
||||
return accumlator
|
||||
|
||||
def hough_transform_a_circle(image: Union[npt.NDArray[np.float64] , npt.NDArray[np.uint8]],
|
||||
edged_image: Union[npt.NDArray[np.float64] , npt.NDArray[np.uint8]],
|
||||
r_start: int, r_end: int, treshold: float) -> npt.NDArray[np.float64]:
|
||||
"""
|
||||
Accepts: image, r_start, r_end
|
||||
Returns: hough space
|
||||
"""
|
||||
|
||||
edged_image = edged_image.copy()
|
||||
edged_image[edged_image < treshold] = 0
|
||||
|
||||
accumlator = np.zeros((edged_image.shape[0], edged_image.shape[1], r_end - r_start))
|
||||
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):
|
||||
print(r)
|
||||
y, x = indices[i]
|
||||
a = x - r * np.cos(ga[y, x])
|
||||
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]:
|
||||
""""
|
||||
Accepts: bw image with lines, n_bins_theta, n_bins_rho, treshold
|
||||
Returns: image points above treshold transformed into hough space
|
||||
"""
|
||||
|
||||
image = image.copy()
|
||||
image[image < treshold] = 0
|
||||
theta_values = np.linspace(-np.pi/2, np.pi/2, n_bins_theta)
|
||||
|
||||
D = np.sqrt(image.shape[0]**2 + image.shape[1]**2)
|
||||
rho_values = np.linspace(-D, D, n_bins_rho)
|
||||
accumulator = np.zeros((n_bins_rho, n_bins_theta), dtype=np.uint64)
|
||||
|
||||
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
|
||||
|
||||
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]]:
|
||||
""""
|
||||
Accepts: bw image with lines, n_bins_theta, n_bins_rho, treshold
|
||||
Returns: image points above treshold transformed into hough space
|
||||
"""
|
||||
|
||||
image = image_with_lines.copy()
|
||||
image[image < treshold] = 0
|
||||
theta_values = np.linspace(-np.pi/2, np.pi/2, n_bins_theta)
|
||||
D = np.sqrt(image.shape[0]**2 + image.shape[1]**2)
|
||||
rho_values = np.linspace(-D, D, n_bins_rho)
|
||||
accumulator = np.zeros((n_bins_rho, n_bins_theta), dtype=np.float64)
|
||||
|
||||
cos_precalculated = np.cos(theta_values)
|
||||
sin_precalculated = np.sin(theta_values)
|
||||
|
||||
indices = np.argwhere(image)
|
||||
|
||||
# Loop through all nonzero pixels above treshold
|
||||
for i in tqdm(range(len(indices)), desc='Hough transform'):
|
||||
y, x = indices[i]
|
||||
angle = (np.mod(gradient_angles[y, x] + np.pi/2 , np.pi)) - np.pi/2
|
||||
|
||||
theta = np.digitize(angle, theta_values) -1
|
||||
rho = np.round(x* cos_precalculated[theta] + y* sin_precalculated[theta]).astype(np.float64)
|
||||
binned_rho = np.digitize(rho, rho_values) - 1 # cuz digitize is returning bin number + 1
|
||||
|
||||
# Add to accumulator
|
||||
accumulator[binned_rho, theta] += gradient_magnitude[y, x]
|
||||
|
||||
return accumulator
|
||||
|
||||
def nonmaxima_suppression_box(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]) -> Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]]:
|
||||
"""
|
||||
Accepts: image with sinusoids in hough space
|
||||
Returns: image with sinusoids
|
||||
"""
|
||||
image = image.copy()
|
||||
|
||||
def get_neighbours() -> list[tuple[int, int]]:
|
||||
neighbours = []
|
||||
for i in range(-1, 2):
|
||||
for j in range(-1, 2):
|
||||
if i != 0 or j != 0:
|
||||
neighbours.append((i, j))
|
||||
return neighbours
|
||||
|
||||
neighbours = get_neighbours()
|
||||
|
||||
for y in range(1, image.shape[0]-1):
|
||||
for x in range(1, image.shape[1]-1):
|
||||
for neighbour in neighbours:
|
||||
if image[y, x] < image[y+neighbour[0], x+neighbour[1]]:
|
||||
image[y, x] = 0
|
||||
break
|
||||
return image
|
||||
|
||||
def retrieve_hough_pairs(original_image: npt.NDArray[np.float64], hough_image: npt.NDArray[np.uint64],
|
||||
treshold: int, n_bins_theta: int, n_bins_rho: int) -> list[tuple[int, int]]:
|
||||
"""
|
||||
Accepts: original image, image with sinusoids in hough space, treshold, n_bins theta in h space, n_bins rho in h space
|
||||
Returns: list of pairs of sinusoids
|
||||
"""
|
||||
theta_values = np.linspace(-np.pi/2, np.pi/2, n_bins_theta)
|
||||
D = np.sqrt(original_image.shape[0]**2 + original_image.shape[1]**2)
|
||||
rho_values = np.linspace(-D, D, n_bins_rho)
|
||||
|
||||
hough_image = hough_image.copy()
|
||||
hough_image[hough_image < treshold] = 0
|
||||
|
||||
y_s, x_s = np.nonzero(hough_image)
|
||||
|
||||
pairs = []
|
||||
for i in range(len(x_s)):
|
||||
x = x_s[i] # theta values
|
||||
y = y_s[i] # rho values
|
||||
|
||||
theta = theta_values[x]
|
||||
rho = rho_values[y]
|
||||
|
||||
pairs.append((rho, theta))
|
||||
return pairs
|
||||
|
||||
def select_best_pairs(image_line_params, n =10):
|
||||
"""
|
||||
Accepts: list of pairs of sinusoids
|
||||
Returns: reduced and prioritized list of pairs of sinusoids
|
||||
"""
|
||||
image_line_params = np.array(image_line_params)
|
||||
# Sorts just kth element so every eleement before kth element is lower than kth element
|
||||
# and every element after kth element is higher than kth element
|
||||
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)
|
||||
|
||||
xs = []
|
||||
ys = []
|
||||
if s != 0:
|
||||
y = int(rho / s)
|
||||
if 0 <= y < h:
|
||||
xs.append(0)
|
||||
ys.append(y)
|
||||
|
||||
y = int((rho - w * c) / s)
|
||||
if 0 <= y < h:
|
||||
xs.append(w - 1)
|
||||
ys.append(y)
|
||||
if c != 0:
|
||||
x = int(rho / c)
|
||||
if 0 <= x < w:
|
||||
xs.append(x)
|
||||
ys.append(0)
|
||||
|
||||
x = int((rho - h * s) / c)
|
||||
if 0 <= x < w:
|
||||
xs.append(x)
|
||||
ys.append(h - 1)
|
||||
|
||||
return xs[:2], ys[:2]
|
||||
|
||||
def hessian_points(image: Union[npt.NDArray[np.float64],
|
||||
npt.NDArray[np.uint8]], sigma: float,
|
||||
treshold: float) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
|
||||
"""
|
||||
Accepts: image, sigma, treshold
|
||||
Returns: image with hessian points
|
||||
"""
|
||||
image = image.copy()
|
||||
(ixx, ixy), (iyx, iyy) = derive_image_second_order(image, sigma)
|
||||
|
||||
# Calculate determinant
|
||||
determinant = ixx * iyy - ixy * iyx
|
||||
|
||||
# Apply nonmaxima suppression
|
||||
points = nonmaxima_suppression_box(determinant)
|
||||
|
||||
points = np.argwhere(points > treshold)
|
||||
|
||||
return determinant.astype(np.float64), points.astype(np.float64)
|
||||
|
|
@ -0,0 +1,9 @@
|
|||
import numpy as np
|
||||
|
||||
def read_data(filename: str):
|
||||
# reads a numpy array from a text file
|
||||
with open(filename) as f:
|
||||
s = f.read()
|
||||
|
||||
return np.fromstring(s, sep=' ')
|
||||
|
||||
Loading…
Reference in New Issue