diff --git a/assignment3/solution.py b/assignment3/solution.py
index 85417bb..33fb9db 100644
--- a/assignment3/solution.py
+++ b/assignment3/solution.py
@@ -192,7 +192,7 @@ def two_b():
     interpolating to more accuracy is not required.
     """
     SIGMA = 1
-    THETA = 0.01
+    THETA = 0.10
 
     museum = uz_image.imread_gray('./images/museum.jpg', uz_image.ImageType.float64)
     museum_edges = uz_image.find_edges_nms(museum, SIGMA, THETA)
@@ -211,7 +211,7 @@ def two_c():
     extract them. Try to avoid explicit for loops as much as possible
     """
     SIGMA = 1
-    THETA = 0.02
+    THETA = 0.10
     T_LOW = 0.04
     T_HIGH = 0.16
 
@@ -230,7 +230,8 @@ def ex3():
     #three_b()
     #three_c()
     #three_d()
-    three_e()
+    #three_e()
+    three_f()
 
 def three_a(): 
     """
@@ -448,7 +449,7 @@ def three_e():
     pier_image_coloured = uz_image.imread('./images/pier.jpg', uz_image.ImageType.float64)
     print('[+] Images loaded')
 
-    bricks_image_edges = uz_image.find_edges_canny(bricks_image_gray, SIGMA, THETA, T_LOW, T_HIGH)
+    bricks_image_edges = uz_image.find_edges_canny(bricks_image_gray, SIGMA, THETA+0.04, T_LOW, T_HIGH)
     pier_image_edges = uz_image.find_edges_canny(pier_image_gray, SIGMA, THETA, T_LOW, T_HIGH)
     print('[+] Edges detected')
 
@@ -508,6 +509,59 @@ def three_e():
         axs[4, 1].plot(xs, ys, 'r', linewidth=0.7)
 
     plt.show()
+
+
+def three_f():
+    """
+    A problem of the Hough transform is that we need a new dimension for
+    each additional parameter in the model, which makes the execution slow for more
+    complex models. We can avoid such parameters if we can reduce the parameter
+    space, e.g. by introducing domain knowledge. Recall from the previous exercise that
+    we can get the local gradient angle besides its magnitude. This angle is perpendicular
+    to the edge and can be used to limit the scope of the parameter ϑ for a specific edge
+    point. We therefore do not have to increase the values of the cells of the entire range
+    of ϑ (calculate multiple values of ρ), but can use the local angle and only work with
+    a single (ρ, ϑ) pair for each edge point.
+    Copy your implementation of the line detector to a new function and modify the
+    algorithm so that it also accepts the matrix of edge angles. Note that the angle values
+    were probably calculated using the np.arctan2(dy, dx) function that returns the
+    values between wider range. You have to adjust the angles so that they are within
+    the [−π/2, π/2] interval. Test the modified function on several images and compare
+    the results with the original implementation.
+    """
+    rectangle_image = uz_image.imread_gray('./images/rectangle.png', uz_image.ImageType.float64)
+
+    image_with_edges_n, derivative_magnitude_n, gradient_angles_n, hough_image_n, hough_image_nms_n, pairs_n, best_pairs_n = uz_image.find_lines_in_image_naive(
+        rectangle_image
+    )
+
+    image_with_edges_i, derivative_magnitude_i, gradient_angles_i, hough_image_i, hough_image_nms_i, pairs_i, best_pairs_i = uz_image.find_lines_in_image_improved(
+        rectangle_image
+    )
+
+    fig, axs = plt.subplots(2, 2)
+
+    axs[0, 0].imshow(hough_image_n)
+    axs[0, 0].set(title='normal')
+    axs[0, 1].imshow(image_with_edges_i)
+    axs[0, 1].set(title='normal')
+
+    axs[1, 0].imshow(rectangle_image, cmap='gray')
+    for param in best_pairs_n:
+        xs, ys = uz_image.get_line_to_plot(param[0], param[1], rectangle_image.shape[0], rectangle_image.shape[1])
+        axs[1, 0].plot(xs, ys, 'r', linewidth=0.7)
+
+    axs[1,1].imshow(rectangle_image, cmap='gray')
+    for param in best_pairs_i:
+        xs, ys = uz_image.get_line_to_plot(param[0], param[1], rectangle_image.shape[0], rectangle_image.shape[1])
+        axs[1, 1].plot(xs, ys, 'r', linewidth=0.7)
+
+    plt.show()
+
+
+
+
+
 # ######## #
 # SOLUTION #
 # ######## #
@@ -515,7 +569,7 @@ def three_e():
 def main():
     #ex1()
     #ex2()
-    ex3()
+    #ex3()
 
 if __name__ == '__main__':
     main()
diff --git a/assignment3/uz_framework/image.py b/assignment3/uz_framework/image.py
index 43e1a49..5ec2f91 100644
--- a/assignment3/uz_framework/image.py
+++ b/assignment3/uz_framework/image.py
@@ -90,7 +90,8 @@ def transform_coloured_image_to_grayscale(image: npt.NDArray[np.float64]) -> npt
 
     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]]:
+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.
     """
@@ -166,7 +167,8 @@ def calculate_best_treshold_using_otsu_method(image: Union[npt.NDArray[np.float6
     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]:
+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
@@ -186,7 +188,8 @@ def get_image_bins_for_loop(image: Union[npt.NDArray[np.float64],  npt.NDArray[n
 
 
 # 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]:
+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
@@ -219,7 +222,8 @@ def get_image_bins(image: Union[npt.NDArray[np.float64],  npt.NDArray[np.uint8]]
 
     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]:
+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
@@ -237,7 +241,8 @@ def get_image_bins_ND(image: Union[npt.NDArray[np.float64],  npt.NDArray[np.uint
 
     return hist / np.sum(hist)
 
-def compare_two_histograms(h1: npt.NDArray[np.float64], h2: npt.NDArray[np.float64], method: DistanceMeasure) -> float:
+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
@@ -257,7 +262,8 @@ def compare_two_histograms(h1: npt.NDArray[np.float64], h2: npt.NDArray[np.float
     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]]:
+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
     """
@@ -333,7 +339,8 @@ def get_gaussian_kernel(sigma: float) -> npt.NDArray[np.float64]:
 
     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]]:
+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 
 
@@ -354,7 +361,8 @@ def sharpen_image(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
 
     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]]:
+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
@@ -397,7 +405,8 @@ def simple_median(signal: npt.NDArray[np.float64], width: int) -> npt.NDArray[np
         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]]:
+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
@@ -424,7 +433,8 @@ def apply_median_method_2D(image:Union[npt.NDArray[np.float64], npt.NDArray[np.u
                 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]]:
+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
@@ -481,7 +491,8 @@ def sp_noise1D(signal, percent=.1) -> npt.NDArray[np.float64]:
     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]]:
+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
@@ -500,7 +511,8 @@ def generate_dirac_impulse(size: int) -> npt.NDArray[np.float64]:
 
     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]]:
+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
@@ -515,7 +527,8 @@ def derive_image_by_x(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8
 
     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]]:
+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
@@ -530,14 +543,19 @@ def derive_image_by_y(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8
 
     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]]]:
+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]]]]:
+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 
@@ -547,7 +565,9 @@ def derive_image_second_order(image: Union[npt.NDArray[np.float64], npt.NDArray[
 
     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]]]:
+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
@@ -556,7 +576,8 @@ def gradient_magnitude(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint
     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]]:
+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
@@ -569,7 +590,8 @@ def find_edges_primitive(image: Union[npt.NDArray[np.float64], npt.NDArray[np.ui
     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]]:
+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
@@ -590,11 +612,11 @@ def find_edges_nms(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
                 if i == 0 or i == 7:
                     return (0, 1), (0, -1)
                 elif i == 1 or i == 2:
-                    return (-1, 0), (1, 0)
-                elif i == 3 or i == 4:
-                    return (-1, 1), (1, -1)
-                elif i == 5 or i == 6:
                     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")
 
 
@@ -614,7 +636,8 @@ def find_edges_nms(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]],
     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]]:
+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
@@ -656,7 +679,8 @@ def hough_transform_a_point(x: int, y: int, n_bins: int) -> npt.NDArray[np.float
     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]:
+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
@@ -691,6 +715,39 @@ def hough_find_lines(image: Union[npt.NDArray[np.float64], npt.NDArray[np.uint8]
 
     return accumulator
 
+def hough_find_lines_i(image_with_lines: npt.NDArray[np.float64], gradient_angles: npt.NDArray[np.float64],
+                       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_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.uint64)
+
+    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]
+              
+        theta = np.digitize(gradient_angles[y, x] / 2, theta_values) -1
+        rho = np.round(x* cos_precalculated[theta] + y* sin_precalculated[theta]).astype(np.int64)
+        binned_rho = np.digitize(rho, rho_values) - 1
+
+        # Add to accumulator
+        accumulator[binned_rho, theta] += 1
+
+    return accumulator
+
 def nonmaxima_suppression_box(image: npt.NDArray[np.uint64]) -> npt.NDArray[np.uint64]:
     """
     Accepts: image with sinusoids in hough space
@@ -717,7 +774,8 @@ def nonmaxima_suppression_box(image: npt.NDArray[np.uint64]) -> npt.NDArray[np.u
     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]]:
+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
@@ -742,7 +800,73 @@ def retrieve_hough_pairs(original_image: npt.NDArray[np.float64], hough_image: n
         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 find_lines_in_image_naive(image: npt.NDArray[np.float64], SIGMA=1, THETA=0.02, T_LOW=0.04,
+                              T_HIGH=0.16, N_BINS_THETA=360, N_BINS_RHO=360, TRESHOLD=0.2):
+    """
+    Aplies all methods to transform image into hough space and find lines
+    """
+
+    image = image.copy()
+
+    # First step: apply canny edge detector
+    image_with_edges = find_edges_canny(image, SIGMA, THETA, T_LOW, T_HIGH)
+
+    # Second step: Retrieve gradient angles
+    derivative_magnitude, gradient_angles = gradient_magnitude(image, SIGMA)
+
+    # Third step: Transform image into hough space
+    hough_image = hough_find_lines(image_with_edges, int(N_BINS_THETA), int(N_BINS_RHO), TRESHOLD)
+
+    # Fourth step: Apply nonmaxima suppression
+    hough_image_nms = nonmaxima_suppression_box(hough_image)
+
+    # Fifth step: Retrieve sigma and theta pairs
+    pairs = retrieve_hough_pairs(image, hough_image_nms, np.max(hough_image_nms) *0.6, int(N_BINS_THETA), int(N_BINS_RHO))
+
+    # Sixth step: select best pairs
+    best_pairs = select_best_pairs(pairs, 10)
+
+    return image_with_edges, derivative_magnitude, gradient_angles, hough_image, hough_image_nms, pairs, best_pairs
+
+def find_lines_in_image_improved(image: npt.NDArray[np.float64], SIGMA=1, THETA=0.02, T_LOW=0.04, T_HIGH=0.16, N_BINS_THETA=360, N_BINS_RHO=360, TRESHOLD=0.2):
+    """
+    Aplies all methods to transform image into hough space and find lines
+    """
+
+    image = image.copy()
+
+    # First step: apply canny edge detector
+    image_with_edges = find_edges_canny(image, SIGMA, THETA, T_LOW, T_HIGH)
+
+    # Second step: Retrieve gradient angles
+    derivative_magnitude, gradient_angles = gradient_magnitude(image, SIGMA)
+
+    # Third step: Transform image into hough space
+    hough_image = hough_find_lines_i(image_with_edges, gradient_angles, N_BINS_THETA, N_BINS_RHO, TRESHOLD)
+
+    # Fourth step: Apply nonmaxima suppression
+    hough_image_nms = nonmaxima_suppression_box(hough_image)
+
+    # Fifth step: Retrieve sigma and theta pairs
+    pairs = retrieve_hough_pairs(image, hough_image_nms, np.max(hough_image_nms) *0.001, N_BINS_THETA, N_BINS_RHO)
+    # Sixth step: select best pairs
+    best_pairs = select_best_pairs(pairs, 10)
+    best_pairs = pairs
+
+    return image_with_edges, derivative_magnitude, gradient_angles, hough_image, hough_image_nms, pairs, best_pairs
+    
 def get_line_to_plot(rho, theta, h, w):
     """
     Accepts: rho, theta, image height h, image width w