From 61ed0bad099d513326e3435ebc5f9e62b2c4f123 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 20 Dec 2020 12:40:06 +0100 Subject: [PATCH 01/61] - triangle class for easier use of barycentric calculations using 3 points --- curvepy/linear_interpolation.py | 240 ++++++++++++++++++-------------- 1 file changed, 133 insertions(+), 107 deletions(-) diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index 995148c..5afad6e 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -47,7 +47,7 @@ def straight_line_function(a: np.ndarray, b: np.ndarray) -> Callable: def collinear_check(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> bool: """ - Calculates the cross product of d1 and d2 to see if all 3 points are collinear. + Calculates the cross product of (b-a) and (c-a) to see if all 3 points are collinear. Parameters ---------- @@ -97,87 +97,6 @@ def ratio(left: np.ndarray, col_point: np.ndarray, right: np.ndarray) -> float: return 0 -def bary_plane_point(bary_coords: np.ndarray, a: np.ndarray, b: np.ndarray, c: np.ndarray) -> np.ndarray: - """ - Given the barycentric coordinates and three points this method will calculate a new point as a barycentric - combination of a, b, c. - - Parameters - ---------- - bary_coords: np.ndarray - The barycentric coordinates corresponding to a, b, c. Have to sum up to 1. - a: np.ndarray - First point of the triangle. - b: np.ndarray - Second point of the triangle. - c: np.ndarray - Third point of the triangle. - - Returns - ------- - np.ndarray: - Barycentric combination of a, b, c with given coordinates. - """ - if sum(bary_coords) != 1: - raise Exception("The barycentric coordinates don't sum up to 1!") - return np.array((bary_coords[0] * a + bary_coords[1] * b + bary_coords[2] * c)) - - -def triangle_area(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> float: - """ - Calculates the area of a triangle defined by the parameters. - - Parameters - ---------- - a: np.ndarray - First point of the triangle. - b: np.ndarray - Second point of the triangle. - c: np.ndarray - Third point of the triangle. - - Returns - ------- - float: - Area of the triangle. - """ - return 1 / 2 * np.linalg.det(np.array([a, b, c])) - - -def get_bary_coords(p: np.ndarray, a: np.ndarray, b: np.ndarray, c: np.ndarray) -> np.ndarray: - """ - Calculates the barycentric coordinates of p with respect to a, b, c. - - Parameters - ---------- - p: np.ndarray - Point of which we want the barycentric coordinates. - a: np.ndarray - First point of the triangle. - b: np.ndarray - Second point of the triangle. - c: np.ndarray - Third point of the triangle. - - Returns - ------- - np.ndarray: - Barycentric coordinates of p with respect to a, b, c. - """ - # Cramer's Rule for 2d easy - # how do i do Cramer's Rule with 3d points??? - tmp1, tmp2, tmp3, tmp4 = p.copy(), a.copy(), b.copy(), c.copy() - for i in range(len(a)): - if a[i] != b[i] != c[i] and a[i - 1] != b[i - 1] != c[i - 1]: - tmp1[i - 2], tmp2[i - 2], tmp3[i - 2], tmp4[i - 2] = 1, 1, 1, 1 - break - abc_area = triangle_area(tmp2, tmp3, tmp4) - if abc_area == 0: - raise Exception("The area of the triangle defined by a, b, c has to be greater than 0!") - return np.array([triangle_area(tmp1, tmp3, tmp4) / abc_area, triangle_area(tmp2, tmp1, tmp4) / abc_area, - triangle_area(tmp2, tmp3, tmp1) / abc_area]) - - # split in 2d and 3d Polygon similar to 2d and 3d bezier? class Polygon: """ @@ -253,34 +172,138 @@ class Polygon: return self.evaluate(0, ts[0]) return Polygon(np.array([self._piece_funcs[i](ts[0]) for i in range(len(ts))])).blossom(ts[1:]) - # CURRENTLY OUT OF ORDER! - def plot_polygon(self, xs: np.ndarray = np.array([0])) -> None: + # # CURRENTLY OUT OF ORDER! + # def plot_polygon(self, xs: np.ndarray = np.array([0])) -> None: + # """ + # Plots the polygon using matplotlib, either in 2D or 3D. Two Plots are given first one with a given number + # of points which will be highlighted on the function, and second is the function as a whole. + # + # Parameters + # ---------- + # xs: np.ndArray + # the points that may be plotted + # """ + # ep = np.array([self.evaluate(x) for x in xs]) + # np.append(ep, self._points) + # ravel_points = self._points.ravel() # the corner points to plot the function + # if self._dim == 2: + # tmp = ep.ravel() + # xs, ys = tmp[0::2], tmp[1::2] + # func_x, func_y = ravel_points[0::2], ravel_points[1::2] + # plt.plot(func_x, func_y) + # plt.plot(xs, ys, 'o', markersize=3) + # if self._dim == 3: + # tmp = ep.ravel() + # xs, ys, zs = tmp[0::3], tmp[1::3], tmp[2::3] + # func_x, func_y, func_z = ravel_points[0::3], ravel_points[1::3], ravel_points[2::3] + # ax = plt.axes(projection='3d') + # ax.plot(func_x, func_y, func_z) + # ax.plot(xs, ys, zs, 'o', markersize=3) + # plt.show() + + +class Triangle(Polygon): + + def __init__(self, points: np.ndarray) -> None: + points_copy = points.copy() + points_copy = np.append(points_copy, [points_copy[0]], axis=0) + super().__init__(points_copy) + + def bary_plane_point(self, bary_coords: np.ndarray) -> np.ndarray: + """ + Given the barycentric coordinates and three points this method will calculate a new point as a barycentric + combination of the triangle points. + + Parameters + ---------- + bary_coords: np.ndarray + The barycentric coordinates corresponding to a, b, c. Have to sum up to 1. + + Returns + ------- + np.ndarray: + Barycentric combination of a, b, c with given coordinates. + """ + if sum(bary_coords) != 1: + raise Exception("The barycentric coordinates don't sum up to 1!") + return np.array((bary_coords[0] * self._points[0] + bary_coords[1] * self._points[1] + + bary_coords[2] * self._points[2])) + + @staticmethod + def area(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> float: + """ + Calculates the area of a triangle defined by the parameters. + + Parameters + ---------- + a: np.ndarray + First point of the triangle. + b: np.ndarray + Second point of the triangle. + c: np.ndarray + Third point of the triangle. + + Returns + ------- + float: + Area of the triangle. + """ + return 1 / 2 * np.linalg.det(np.array([a, b, c])) + + def map_to_plane(self, p: np.ndarray): """ - Plots the polygon using matplotlib, either in 2D or 3D. Two Plots are given first one with a given number - of points which will be highlighted on the function, and second is the function as a whole. + This method projects p and the points of the triangle on a plane, so that cramer's rule can easily be applied + to them in order to calculate the area of every 3 out of the 4 points. Parameters ---------- - xs: np.ndArray - the points that may be plotted + p: np.ndarray + Additional point that should be on the same plane as the triangle. + + Returns + ------- + np.ndarray: + Copys of p and the triangle points now mapped on to a plane. """ - ep = np.array([self.evaluate(x) for x in xs]) - np.append(ep, self._points) - ravel_points = self._points.ravel() # the corner points to plot the function - if self._dim == 2: - tmp = ep.ravel() - xs, ys = tmp[0::2], tmp[1::2] - func_x, func_y = ravel_points[0::2], ravel_points[1::2] - plt.plot(func_x, func_y) - plt.plot(xs, ys, 'o', markersize=3) - if self._dim == 3: - tmp = ep.ravel() - xs, ys, zs = tmp[0::3], tmp[1::3], tmp[2::3] - func_x, func_y, func_z = ravel_points[0::3], ravel_points[1::3], ravel_points[2::3] - ax = plt.axes(projection='3d') - ax.plot(func_x, func_y, func_z) - ax.plot(xs, ys, zs, 'o', markersize=3) - plt.show() + p_copy, a, b, c = p.copy(), self._points[0].copy(), self._points[1].copy(), self._points[2].copy() + print(self._points[0][-1] != self._points[2][-1]) + # how to map the points on a plane without collapsing points on another? + for i in range(len(self._points)): + if self._points[0][i] != self._points[1][i] != self._points[2][i] \ + and self._points[0][i - 1] != self._points[1][i - 1] != self._points[2][i - 1]: + p_copy[i - 2], a[i - 2], b[i - 2], c[i - 2] = 1, 1, 1, 1 + break + return p_copy, a, b, c + + def get_bary_coords(self, p: np.ndarray) -> np.ndarray: + """ + Calculates the barycentric coordinates of p with respect to the points defining the triangle. + + Parameters + ---------- + p: np.ndarray + Point of which we want the barycentric coordinates. + + Returns + ------- + np.ndarray: + Barycentric coordinates of p with respect to a, b, c. + """ + if np.shape(p)[0] != self._dim: + raise Exception("p has to have the same dimension as the triangle points!") + elif np.shape(p)[0] == 3: + p_copy, a, b, c = self.map_to_plane(p) + elif np.shape(p)[0] == 2: + p_copy, a, b, c = p.copy(), self._points[0].copy(), self._points[1].copy(), self._points[2].copy() + else: + raise Exception("Dimension hast to be 2 or 3!") + + abc_area = self.area(a, b, c) + if abc_area == 0: + raise Exception("The area of the triangle defined by a, b, c has to be greater than 0!") + + return np.array([self.area(p_copy, b, c) / abc_area, self.area(a, p_copy, c) / abc_area, + self.area(a, b, p_copy) / abc_area]) def ratio_test() -> None: @@ -324,10 +347,13 @@ def init() -> None: # Blossom testing # blossom_testing() + # triangle test + t = Triangle(np.array([a, b, c])) + # barycentric coords test - # print(bary_plane_point(coords, a, b, c)) + print(t.bary_plane_point(coords)) # print(np.linalg.det(np.array([a, b, c]))) - print(get_bary_coords(bary_plane_point(coords, a, b, c), a, b, c)) + print(t.get_bary_coords(t.bary_plane_point(coords))) if __name__ == "__main__": -- GitLab From 886cb4b76befaff68553cfca3fd3349a5210cf20 Mon Sep 17 00:00:00 2001 From: Icarus1899 <63288817+Icarus1899@users.noreply.github.com> Date: Fri, 8 Jan 2021 13:34:28 +0100 Subject: [PATCH 02/61] - implemented methods for the calculation of the barycentric coordinates of p work now in 2D and 3D --- curvepy/linear_interpolation.py | 51 ++++++++++++++++++++++----------- 1 file changed, 35 insertions(+), 16 deletions(-) diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index 5afad6e..0433932 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -250,10 +250,12 @@ class Triangle(Polygon): """ return 1 / 2 * np.linalg.det(np.array([a, b, c])) - def map_to_plane(self, p: np.ndarray): + def map_parallel_to_axis_plane(self, p: np.ndarray): """ - This method projects p and the points of the triangle on a plane, so that cramer's rule can easily be applied - to them in order to calculate the area of every 3 out of the 4 points. + This method projects p and the points of the triangle on a plane, for example the y-plane with distance 1 for + all points of the triangle to the plane, so that cramer's rule can easily be applied to them + in order to calculate the area of every 3 out of the 4 points. + But this method does not overwrite the self._points. Parameters ---------- @@ -263,11 +265,10 @@ class Triangle(Polygon): Returns ------- np.ndarray: - Copys of p and the triangle points now mapped on to a plane. + Copy of p and the triangle points now mapped on to a plane. """ p_copy, a, b, c = p.copy(), self._points[0].copy(), self._points[1].copy(), self._points[2].copy() - print(self._points[0][-1] != self._points[2][-1]) - # how to map the points on a plane without collapsing points on another? + # is there another way? for i in range(len(self._points)): if self._points[0][i] != self._points[1][i] != self._points[2][i] \ and self._points[0][i - 1] != self._points[1][i - 1] != self._points[2][i - 1]: @@ -275,35 +276,53 @@ class Triangle(Polygon): break return p_copy, a, b, c - def get_bary_coords(self, p: np.ndarray) -> np.ndarray: + def check_points_for_area_calc(self, p): """ - Calculates the barycentric coordinates of p with respect to the points defining the triangle. + This method checks if the point p and the points of the triangle have the right dimension and will make them so + that cramer's rule can be applied to them. Parameters ---------- p: np.ndarray - Point of which we want the barycentric coordinates. + Additional point that has to be on the same plane as the triangle. Returns ------- - np.ndarray: - Barycentric coordinates of p with respect to a, b, c. + np.ndarrays: + The triangle points and p so that cramer's rule can be used. """ if np.shape(p)[0] != self._dim: raise Exception("p has to have the same dimension as the triangle points!") elif np.shape(p)[0] == 3: - p_copy, a, b, c = self.map_to_plane(p) + return self.map_parallel_to_axis_plane(p) elif np.shape(p)[0] == 2: - p_copy, a, b, c = p.copy(), self._points[0].copy(), self._points[1].copy(), self._points[2].copy() + return np.append(p.copy(), [1]), np.append(self._points[0].copy(), [1]), np.append(self._points[1].copy(), [1]),\ + np.append(self._points[2].copy(), [1]) else: - raise Exception("Dimension hast to be 2 or 3!") + raise Exception("Dimension has to be 2 or 3!") + + def get_bary_coords(self, p: np.ndarray) -> np.ndarray: + """ + Calculates the barycentric coordinates of p with respect to the points defining the triangle. + + Parameters + ---------- + p: np.ndarray + Point of which we want the barycentric coordinates. + + Returns + ------- + np.ndarray: + Barycentric coordinates of p with respect to a, b, c. + """ + p_copy, a, b, c = self.check_points_for_area_calc(p) abc_area = self.area(a, b, c) if abc_area == 0: raise Exception("The area of the triangle defined by a, b, c has to be greater than 0!") return np.array([self.area(p_copy, b, c) / abc_area, self.area(a, p_copy, c) / abc_area, - self.area(a, b, p_copy) / abc_area]) + self.area(a, b, p_copy) / abc_area])[0:np.shape(p)[0]] def ratio_test() -> None: @@ -333,7 +352,7 @@ def blossom_testing() -> None: def init() -> None: coords = np.array([1 / 3, 1 / 3, 1 / 3]) - a, b, c = np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 0, 0]) + a, b, c = np.array([2, 1, 1]), np.array([4, 3, 1]), np.array([5, 1, 1]) # straight_line_point_test() # ratio_test() # test_points = np.array([[0, 0, 0], [1, 1, 1], [3, 4, 4], [5, -2, -2]]) -- GitLab From b1ebc9bb83a1436fb427c302c3ff18424e4e184b Mon Sep 17 00:00:00 2001 From: Icarus1899 <63288817+Icarus1899@users.noreply.github.com> Date: Fri, 8 Jan 2021 13:53:00 +0100 Subject: [PATCH 03/61] - fix in get_bary_coords --- curvepy/linear_interpolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index 0433932..51a1ec3 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -322,7 +322,7 @@ class Triangle(Polygon): raise Exception("The area of the triangle defined by a, b, c has to be greater than 0!") return np.array([self.area(p_copy, b, c) / abc_area, self.area(a, p_copy, c) / abc_area, - self.area(a, b, p_copy) / abc_area])[0:np.shape(p)[0]] + self.area(a, b, p_copy) / abc_area]) def ratio_test() -> None: @@ -352,7 +352,7 @@ def blossom_testing() -> None: def init() -> None: coords = np.array([1 / 3, 1 / 3, 1 / 3]) - a, b, c = np.array([2, 1, 1]), np.array([4, 3, 1]), np.array([5, 1, 1]) + a, b, c = np.array([2, 1]), np.array([4, 3]), np.array([5, 1]) # straight_line_point_test() # ratio_test() # test_points = np.array([[0, 0, 0], [1, 1, 1], [3, 4, 4], [5, -2, -2]]) -- GitLab From 6ec3b146e6b4da3897233e0b36a98b765fa36f02 Mon Sep 17 00:00:00 2001 From: Icarus1899 <63288817+Icarus1899@users.noreply.github.com> Date: Fri, 8 Jan 2021 14:19:19 +0100 Subject: [PATCH 04/61] - addition of a rounding error in bary_plane_point --- curvepy/linear_interpolation.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index 51a1ec3..7361e28 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -224,7 +224,8 @@ class Triangle(Polygon): np.ndarray: Barycentric combination of a, b, c with given coordinates. """ - if sum(bary_coords) != 1: + eps = 0.0000000000000001 # rounding error + if sum(bary_coords) - eps != 1 and sum(bary_coords) + eps != 1: raise Exception("The barycentric coordinates don't sum up to 1!") return np.array((bary_coords[0] * self._points[0] + bary_coords[1] * self._points[1] + bary_coords[2] * self._points[2])) @@ -296,8 +297,9 @@ class Triangle(Polygon): elif np.shape(p)[0] == 3: return self.map_parallel_to_axis_plane(p) elif np.shape(p)[0] == 2: - return np.append(p.copy(), [1]), np.append(self._points[0].copy(), [1]), np.append(self._points[1].copy(), [1]),\ - np.append(self._points[2].copy(), [1]) + return np.append(p.copy(), [1]), np.append(self._points[0].copy(), [1]), np.append(self._points[1].copy(), + [1]), \ + np.append(self._points[2].copy(), [1]) else: raise Exception("Dimension has to be 2 or 3!") @@ -351,8 +353,8 @@ def blossom_testing() -> None: def init() -> None: - coords = np.array([1 / 3, 1 / 3, 1 / 3]) - a, b, c = np.array([2, 1]), np.array([4, 3]), np.array([5, 1]) + coords = np.array([2 / 3, 1 / 6, 1 / 6]) + a, b, c = np.array([2, 1, 1]), np.array([4, 3, 2]), np.array([5, 1, 1]) # straight_line_point_test() # ratio_test() # test_points = np.array([[0, 0, 0], [1, 1, 1], [3, 4, 4], [5, -2, -2]]) -- GitLab From d4bedc1d49062e3ded0b26570aad148f8afd656e Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Fri, 8 Jan 2021 19:05:48 +0100 Subject: [PATCH 05/61] - reformating --- curvepy/linear_interpolation.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index 7361e28..cc5700d 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -97,7 +97,6 @@ def ratio(left: np.ndarray, col_point: np.ndarray, right: np.ndarray) -> float: return 0 -# split in 2d and 3d Polygon similar to 2d and 3d bezier? class Polygon: """ Class for creating a 2D or 3D Polygon. @@ -269,7 +268,6 @@ class Triangle(Polygon): Copy of p and the triangle points now mapped on to a plane. """ p_copy, a, b, c = p.copy(), self._points[0].copy(), self._points[1].copy(), self._points[2].copy() - # is there another way? for i in range(len(self._points)): if self._points[0][i] != self._points[1][i] != self._points[2][i] \ and self._points[0][i - 1] != self._points[1][i - 1] != self._points[2][i - 1]: @@ -292,14 +290,11 @@ class Triangle(Polygon): np.ndarrays: The triangle points and p so that cramer's rule can be used. """ - if np.shape(p)[0] != self._dim: - raise Exception("p has to have the same dimension as the triangle points!") - elif np.shape(p)[0] == 3: + if self._dim == 3: return self.map_parallel_to_axis_plane(p) - elif np.shape(p)[0] == 2: - return np.append(p.copy(), [1]), np.append(self._points[0].copy(), [1]), np.append(self._points[1].copy(), - [1]), \ - np.append(self._points[2].copy(), [1]) + elif self._dim == 2: + return np.append(p.copy(), [1]), np.append(self._points[0].copy(), [1]), \ + np.append(self._points[1].copy(), [1]), np.append(self._points[2].copy(), [1]) else: raise Exception("Dimension has to be 2 or 3!") @@ -354,7 +349,7 @@ def blossom_testing() -> None: def init() -> None: coords = np.array([2 / 3, 1 / 6, 1 / 6]) - a, b, c = np.array([2, 1, 1]), np.array([4, 3, 2]), np.array([5, 1, 1]) + a, b, c = np.array([2, 1, 3]), np.array([4, 3, 3]), np.array([5, 1, 3]) # straight_line_point_test() # ratio_test() # test_points = np.array([[0, 0, 0], [1, 1, 1], [3, 4, 4], [5, -2, -2]]) -- GitLab From 3faea149fc40b8ebc28d1387235602cd3ce8f7b6 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 10 Jan 2021 10:25:40 +0100 Subject: [PATCH 06/61] - small changes in bary_plane_point --- curvepy/linear_interpolation.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index cc5700d..c527fd9 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -3,6 +3,7 @@ from mpl_toolkits.mplot3d import Axes3D from typing import Tuple, Callable import functools import matplotlib.pyplot as plt +import sys def straight_line_point(a: np.ndarray, b: np.ndarray, t: float = 0.5) -> np.ndarray: @@ -223,8 +224,8 @@ class Triangle(Polygon): np.ndarray: Barycentric combination of a, b, c with given coordinates. """ - eps = 0.0000000000000001 # rounding error - if sum(bary_coords) - eps != 1 and sum(bary_coords) + eps != 1: + eps = sys.float_info.epsilon + if np.sum(bary_coords) - eps != 1 and np.sum(bary_coords) + eps != 1: raise Exception("The barycentric coordinates don't sum up to 1!") return np.array((bary_coords[0] * self._points[0] + bary_coords[1] * self._points[1] + bary_coords[2] * self._points[2])) @@ -250,7 +251,7 @@ class Triangle(Polygon): """ return 1 / 2 * np.linalg.det(np.array([a, b, c])) - def map_parallel_to_axis_plane(self, p: np.ndarray): + def squash_parallel_to_axis_plane(self, p: np.ndarray): """ This method projects p and the points of the triangle on a plane, for example the y-plane with distance 1 for all points of the triangle to the plane, so that cramer's rule can easily be applied to them @@ -291,7 +292,7 @@ class Triangle(Polygon): The triangle points and p so that cramer's rule can be used. """ if self._dim == 3: - return self.map_parallel_to_axis_plane(p) + return self.squash_parallel_to_axis_plane(p) elif self._dim == 2: return np.append(p.copy(), [1]), np.append(self._points[0].copy(), [1]), \ np.append(self._points[1].copy(), [1]), np.append(self._points[2].copy(), [1]) -- GitLab From ab7964b990ec8ff08ce5c2f00968317bc3612d44 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Mon, 12 Apr 2021 13:27:17 +0200 Subject: [PATCH 07/61] dirichlet tessellation beginning recursion --- curvepy/linear_interpolation.py | 41 ++++++++++++++++++++++++++------- 1 file changed, 33 insertions(+), 8 deletions(-) diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index c527fd9..0740389 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -1,3 +1,5 @@ +from plistlib import Data + import numpy as np from mpl_toolkits.mplot3d import Axes3D from typing import Tuple, Callable @@ -233,7 +235,8 @@ class Triangle(Polygon): @staticmethod def area(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> float: """ - Calculates the area of a triangle defined by the parameters. + Calculates the "area" of a triangle defined by the parameters. All three points have to be on a plane parallel + to an axis-plane! Parameters ---------- @@ -247,7 +250,7 @@ class Triangle(Polygon): Returns ------- float: - Area of the triangle. + "Area" of the triangle. """ return 1 / 2 * np.linalg.det(np.array([a, b, c])) @@ -255,7 +258,7 @@ class Triangle(Polygon): """ This method projects p and the points of the triangle on a plane, for example the y-plane with distance 1 for all points of the triangle to the plane, so that cramer's rule can easily be applied to them - in order to calculate the area of every 3 out of the 4 points. + in order to calculate the area of the triangle corresponding to every 3 out of the 4 points. But this method does not overwrite the self._points. Parameters @@ -323,6 +326,25 @@ class Triangle(Polygon): self.area(a, b, p_copy) / abc_area]) +class Dirichlet_tessellation: + + def __init__(self, points: np.ndarray) -> None: + self.points_copy = points.copy() + self.tessellation = [] + self.recursive_dirichlet_tessellation(0) + + def recursive_dirichlet_tessellation(self, index: int) -> None: + if index == len(self.points_copy): + return + if not index: + self.tessellation.append(Polygon(np.array([[float('inf'), float('inf')], + [float('inf'), float('-inf')], + [float('-inf'), float('-inf')], + [float('-inf'), float('inf')]]))) + self.recursive_dirichlet_tessellation(1) + + + def ratio_test() -> None: left = np.array([0, 0, 0]) right = np.array([1, 1, 1]) @@ -349,8 +371,9 @@ def blossom_testing() -> None: def init() -> None: - coords = np.array([2 / 3, 1 / 6, 1 / 6]) - a, b, c = np.array([2, 1, 3]), np.array([4, 3, 3]), np.array([5, 1, 3]) + # coords = np.array([2 / 3, 1 / 6, 1 / 6]) + # a, b, c = np.array([2, 1]), np.array([4, 3]), np.array([5, 1]) + # straight_line_point_test() # ratio_test() # test_points = np.array([[0, 0, 0], [1, 1, 1], [3, 4, 4], [5, -2, -2]]) @@ -365,12 +388,14 @@ def init() -> None: # blossom_testing() # triangle test - t = Triangle(np.array([a, b, c])) + # t = Triangle(np.array([a, b, c])) # barycentric coords test - print(t.bary_plane_point(coords)) + # print(t.bary_plane_point(coords)) # print(np.linalg.det(np.array([a, b, c]))) - print(t.get_bary_coords(t.bary_plane_point(coords))) + # print(t.get_bary_coords(t.bary_plane_point(coords))) + + test_tes = Dirichlet_tessellation(np.array([[0, 0]])) if __name__ == "__main__": -- GitLab From b7a4c90fbbd6097666b83d6b3d663d0b9e671204 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 25 Apr 2021 15:29:56 +0200 Subject: [PATCH 08/61] dirichlet tessellation change from recusrion to iterative thought about triangle as way to know when points are neighbours implementation next week --- curvepy/linear_interpolation.py | 58 ++++++++++++++++++++++++++------- 1 file changed, 47 insertions(+), 11 deletions(-) diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index 0740389..8b8d9b9 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -326,23 +326,59 @@ class Triangle(Polygon): self.area(a, b, p_copy) / abc_area]) +class tile_edge: + + def __init__(self) -> None: + """ + start: np.ndarray + start point + end: np.ndarray + end point + """ + self.start = None + self.end = None + + +class tile: + + def __init__(self, center: np.ndarray) -> None: + self.center = center + self.edges = [] + + def add_edge(self, edge: tile_edge): + self.edges.append(edge) + + class Dirichlet_tessellation: def __init__(self, points: np.ndarray) -> None: self.points_copy = points.copy() self.tessellation = [] - self.recursive_dirichlet_tessellation(0) + self.iterative_dirichlet_tessellation() + + def iterative_dirichlet_tessellation(self) -> None: + for p in self.points_copy: + print(f"Punkt zum Einfügen: {p}") + if not self.tessellation: + self.tessellation.append(tile(p)) + continue + + print(f"Liste: {[t.center for t in self.tessellation]}") + min_idx = np.argmin([np.linalg.norm(p - t.center) for t in self.tessellation]) + + print(f"mindinmaler Index: {min_idx}") + self.tessellation.append(tile(p)) + + # das muss eig unter tiles dann können wir uns die schleifen sparen und updaten einfach immer nur die nachbarn + # aber maxi, WaNn siNd TiLeS keInE nAcHbArN MeHR? - def recursive_dirichlet_tessellation(self, index: int) -> None: - if index == len(self.points_copy): - return - if not index: - self.tessellation.append(Polygon(np.array([[float('inf'), float('inf')], - [float('inf'), float('-inf')], - [float('-inf'), float('-inf')], - [float('-inf'), float('inf')]]))) - self.recursive_dirichlet_tessellation(1) + neighbours = [] + for e in self.tessellation[min_idx].edges: + for t in self.tessellation: + if e in t.edges: + neighbours.append(t) + # hier kommt totaler schwachsinn def ratio_test() -> None: @@ -395,7 +431,7 @@ def init() -> None: # print(np.linalg.det(np.array([a, b, c]))) # print(t.get_bary_coords(t.bary_plane_point(coords))) - test_tes = Dirichlet_tessellation(np.array([[0, 0]])) + test_tes = Dirichlet_tessellation(np.array([[0, 0], [1, 2], [-2, -3], [-3, -3]])) if __name__ == "__main__": -- GitLab From db961c08fc9f3f1eb482281a88328c7f5ce9b45e Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 15 May 2021 17:07:18 +0200 Subject: [PATCH 09/61] Changes by lquenti --- curvepy/example_algo.py | 52 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 curvepy/example_algo.py diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py new file mode 100644 index 0000000..8de3cf6 --- /dev/null +++ b/curvepy/example_algo.py @@ -0,0 +1,52 @@ +from dataclasses import dataclass +import numpy as np +from typing import List, Set + +@dataclass +class tile: + center: np.ndarray # 2D float + neighbours: list # of tiles + + +class dirichlet_tessellation: + + def __init__(self): + self.tiles : List[tile] = [] + self.valid_triangulation: List[Set[np.ndarray, np.ndarray]] = [] + + + def append_point(self, p: np.ndarray): + if not self.tiles: + self.tiles.append(tile(p, [])) + return + + min_idx = np.argmin([np.linalg.norm(p - t.center) for t in self.tiles]) + act_tile = self.tiles[min_idx] + + # ERSTELLE TRIANGULATION ZWISCHEN p UND act_tile + + # ZUERST ALLE KOLLISIONSFREIEN + for n in act_tile.neighbours: + if self._intersect(n,p): + continue + # FÜGE TRIANGULATION HINZU + + # DANN ALLE MIT KOLLISIONEN + for n in act_tile.neighbours: + if self._intersect(n,p): + if new_is_more_gleichseitig_sprich_abs_minus_60_minimieren(self, (n,p)): + # replatziere alte Kante mit neuen kanten + # TODO was passiert wenn wir mehr als eine Kollision haben? -> trashen wir die neue Kante dann? + + def _intersect(self, n,p): + return True + +# sum(abs(60-x) for x in [39.80, 71.88, 68.31, 48.63, 68.31, 63.06]) +# sum(abs(60-y) for y in [39.04, 116.94, 24.01, 32.84, 108.11, 39.04]) + +# sum(abs(60-x) for x in [76.95, 85.12, 17.93, 18.19, 140.55, 21.25]) +# sum(abs(60-y) for y in [30.90, 36.12, 112.98, 46.05, 106.37, 27.57]) + +""" + +""" \ No newline at end of file -- GitLab From b6fc9271624e5bb63c6ce2b0b78ee82e4276b01c Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 16 May 2021 16:32:14 +0200 Subject: [PATCH 10/61] =?UTF-8?q?-=20wir=20haben=20dreicke=20und=20nachbar?= =?UTF-8?q?n=20ziemlich=20durch=20-=20n=C3=A4chstes=20mal=20gucken=20wie?= =?UTF-8?q?=20wir=20nachbarn=20vern=C3=BCnftig=20handlen=20->=20datenstruk?= =?UTF-8?q?tur?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- curvepy/example_algo.py | 93 +++++++++++++++++++++++++++++++---------- 1 file changed, 70 insertions(+), 23 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 8de3cf6..3b1690a 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -2,10 +2,11 @@ from dataclasses import dataclass import numpy as np from typing import List, Set + @dataclass class tile: - center: np.ndarray # 2D float - neighbours: list # of tiles + center: np.ndarray # 2D float + neighbours: list # of tiles class dirichlet_tessellation: @@ -14,7 +15,6 @@ class dirichlet_tessellation: self.tiles : List[tile] = [] self.valid_triangulation: List[Set[np.ndarray, np.ndarray]] = [] - def append_point(self, p: np.ndarray): if not self.tiles: self.tiles.append(tile(p, [])) @@ -23,30 +23,77 @@ class dirichlet_tessellation: min_idx = np.argmin([np.linalg.norm(p - t.center) for t in self.tiles]) act_tile = self.tiles[min_idx] - # ERSTELLE TRIANGULATION ZWISCHEN p UND act_tile + # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_tile + self.valid_triangulation.append((p, nearest_tile.center)) - # ZUERST ALLE KOLLISIONSFREIEN - for n in act_tile.neighbours: - if self._intersect(n,p): - continue - # FÜGE TRIANGULATION HINZU + # To make a more educated guess we solve any collisions __after__ adding the trivial connections + collisions_to_check = [] - # DANN ALLE MIT KOLLISIONEN - for n in act_tile.neighbours: - if self._intersect(n,p): - if new_is_more_gleichseitig_sprich_abs_minus_60_minimieren(self, (n,p)): - # replatziere alte Kante mit neuen kanten - # TODO was passiert wenn wir mehr als eine Kollision haben? -> trashen wir die neue Kante dann? + for neighbour in nearest_tile.neighbours: + all_collisions = [x for x in self.valid_triangulation if self._intersect(neighbour.center, p, *x)] + if len(all_collisions) == 0: + self.valid_triangulation.append((p, neighbour.center)) + neighbour.neighbours.append( + p_tile) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind + elif len(all_collisions) == 1: + # 1 collision could be still fine + collisions_to_check.append([p, neighbour.center, *all_collisions[0]]) + # More than 1 collision is always not the best possible triangulation TODO is that true? Yesn't - def _intersect(self, n,p): - return True + for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: + new_triangles = ( + # first triangle + ((p, neighbour), (p, collision_edge_p1)), + ((p, neighbour), (neighbour, collision_edge_p1)), + ((p, neighbour), (neighbour, collision_edge_p1)), + # second triangle + ((p, neighbour), (p, collision_edge_p2)), + ((p, collision_edge_p2), (neighbour, collision_edge_p2)), + ((p, neighbour), (neighbour, collision_edge_p2)) + ) -# sum(abs(60-x) for x in [39.80, 71.88, 68.31, 48.63, 68.31, 63.06]) -# sum(abs(60-y) for y in [39.04, 116.94, 24.01, 32.84, 108.11, 39.04]) + old_triangles = ( + # first triangle + ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, p)), + ((collision_edge_p1, p), (collision_edge_p2, p)), + ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, p)), + # second triangle + ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, neighbour)), + ((collision_edge_p1, neighbour), (collision_edge_p2, neighbour)), + ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, neighbour)) + ) -# sum(abs(60-x) for x in [76.95, 85.12, 17.93, 18.19, 140.55, 21.25]) -# sum(abs(60-y) for y in [30.90, 36.12, 112.98, 46.05, 106.37, 27.57]) + rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) + new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) + if new_is_more_equilateral: + # dann kannte col_edge_p2 - col_edge_p1 entfernen und kannte p - neighbour hinzufügen + collision_edge_p1 + ... -""" + @staticmethod + def _angle(edge1, edge2): + return abs(180 / np.pi * (np.arctan((edge1[1][1] - edge1[0][1]) / (edge1[1][0] - edge1[0][0])) + - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0])))) -""" \ No newline at end of file + @staticmethod + def get_all_edgepairs(A, B, C): + return [(A, B), (B, C), (A, C)] + + @staticmethod + def _rate_triangle(A, B, C): + ... + + @staticmethod + def _intersect(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray) -> bool: + # First we vertical stack the points in an array + vertical_stack = np.vstack([p1, p2, p3, p4]) + # Then we transform them to homogeneous coordinates, to perform a little trick + homogeneous = np.hstack((vertical_stack, np.ones((4, 1)))) + # having our points in this form we can get the lines through the cross product + line_1, line_2 = np.cross(homogeneous[0], homogeneous[1]), np.cross(homogeneous[2], homogeneous[3]) + # when we calculate the cross product of the lines we get intersect point + x, y, z = np.cross(line_1, line_2) + if z == 0: + return False + # we divide with z to turn back to 2D space + return True -- GitLab From ff9cfdaddacddf5bbc176ed697204e51bed79d9e Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 16 May 2021 16:32:51 +0200 Subject: [PATCH 11/61] =?UTF-8?q?-=20n=C3=A4chstes=20mal=20nachbarn=20hand?= =?UTF-8?q?len=20->=20datenstruktur?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- curvepy/example_algo.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 3b1690a..115bc05 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -1,6 +1,6 @@ from dataclasses import dataclass import numpy as np -from typing import List, Set +from typing import List, Set, Union, Dict, Tuple @dataclass @@ -12,8 +12,8 @@ class tile: class dirichlet_tessellation: def __init__(self): - self.tiles : List[tile] = [] - self.valid_triangulation: List[Set[np.ndarray, np.ndarray]] = [] + self.tiles: List[tile] = [] + self.valid_triangulation: List[Tuple[np.ndarray, np.ndarray]] = [] def append_point(self, p: np.ndarray): if not self.tiles: @@ -21,7 +21,11 @@ class dirichlet_tessellation: return min_idx = np.argmin([np.linalg.norm(p - t.center) for t in self.tiles]) - act_tile = self.tiles[min_idx] + nearest_tile = self.tiles[min_idx] + p_tile = tile(center=p, neighbours=[]) + nearest_tile.neighbours.append(p_tile) + p_tile.neighbours.append(nearest_tile) + self.tiles.append(p_tile) # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_tile self.valid_triangulation.append((p, nearest_tile.center)) -- GitLab From 6235a78049817122f9f84b85f30819946c1ec71b Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 14:35:56 +0200 Subject: [PATCH 12/61] =?UTF-8?q?-=20n=C3=A4chstes=20mal=20nachbarn=20hand?= =?UTF-8?q?len=20->=20datenstruktur?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- curvepy/example_algo.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 115bc05..67db006 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -47,13 +47,20 @@ class dirichlet_tessellation: for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: new_triangles = ( # first triangle - ((p, neighbour), (p, collision_edge_p1)), - ((p, neighbour), (neighbour, collision_edge_p1)), - ((p, neighbour), (neighbour, collision_edge_p1)), + (p, neighbour), + (p, collision_edge_p1), + (neighbour, collision_edge_p1), + # ((p, neighbour), (p, collision_edge_p1)), + # ((p, neighbour), (neighbour, collision_edge_p1)), + # ((p, neighbour), (neighbour, collision_edge_p1)), # second triangle - ((p, neighbour), (p, collision_edge_p2)), - ((p, collision_edge_p2), (neighbour, collision_edge_p2)), - ((p, neighbour), (neighbour, collision_edge_p2)) + (p, neighbour), + (p, collision_edge_p2), + (neighbour, collision_edge_p2) + + # ((p, neighbour), (p, collision_edge_p2)), + # ((p, collision_edge_p2), (neighbour, collision_edge_p2)), + # ((p, neighbour), (neighbour, collision_edge_p2)) ) old_triangles = ( @@ -70,7 +77,7 @@ class dirichlet_tessellation: rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) if new_is_more_equilateral: - # dann kannte col_edge_p2 - col_edge_p1 entfernen und kannte p - neighbour hinzufügen + # dann kante col_edge_p2 - col_edge_p1 entfernen und kante p - neighbour hinzufügen collision_edge_p1 ... -- GitLab From 06423ef4c5eeada8045630a74c5dec049b682b51 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 14:36:04 +0200 Subject: [PATCH 13/61] update zwischendurch --- curvepy/testrange.py | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 curvepy/testrange.py diff --git a/curvepy/testrange.py b/curvepy/testrange.py new file mode 100644 index 0000000..80af607 --- /dev/null +++ b/curvepy/testrange.py @@ -0,0 +1,3 @@ + +x = list(itertools.chain(*([(A,B), (B,C), (A,C)] for A,B,C in [((p, neighbour), (p, fst_collision), (neighbour, fst_collision)), + ((p, neighbour), (p, snd_collision), (neighbour, snd_collision))]))) \ No newline at end of file -- GitLab From f96c730c3a1e6bdea35a9fee2d6675dd942adca0 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 14:42:33 +0200 Subject: [PATCH 14/61] update zwischendurch 2 --- curvepy/example_algo.py | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 67db006..031f2ff 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -54,24 +54,14 @@ class dirichlet_tessellation: # ((p, neighbour), (neighbour, collision_edge_p1)), # ((p, neighbour), (neighbour, collision_edge_p1)), # second triangle - (p, neighbour), - (p, collision_edge_p2), - (neighbour, collision_edge_p2) - - # ((p, neighbour), (p, collision_edge_p2)), - # ((p, collision_edge_p2), (neighbour, collision_edge_p2)), - # ((p, neighbour), (neighbour, collision_edge_p2)) + self.get_all_edgepairs(p, neighbour, collision_edge_p2) ) old_triangles = ( # first triangle - ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, p)), - ((collision_edge_p1, p), (collision_edge_p2, p)), - ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, p)), + self.get_all_edgepairs(collision_edge_p1, collision_edge_p2, p), # second triangle - ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, neighbour)), - ((collision_edge_p1, neighbour), (collision_edge_p2, neighbour)), - ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, neighbour)) + self.get_all_edgepairs(collision_edge_p1, collision_edge_p2, neighbour) ) rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) -- GitLab From 1bc33bf51487b32fa3189f647891773309d29733 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 14:42:56 +0200 Subject: [PATCH 15/61] easy dreiecke erstellen --- curvepy/example_algo.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 031f2ff..fb19fc2 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -47,12 +47,7 @@ class dirichlet_tessellation: for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: new_triangles = ( # first triangle - (p, neighbour), - (p, collision_edge_p1), - (neighbour, collision_edge_p1), - # ((p, neighbour), (p, collision_edge_p1)), - # ((p, neighbour), (neighbour, collision_edge_p1)), - # ((p, neighbour), (neighbour, collision_edge_p1)), + self.get_all_edgepairs(p, neighbour, collision_edge_p1), # second triangle self.get_all_edgepairs(p, neighbour, collision_edge_p2) ) -- GitLab From e7a3a47b350a0b8c75152cef1a3133fa8abd64e2 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 14:45:58 +0200 Subject: [PATCH 16/61] einzeilige Dreiecke --- curvepy/example_algo.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index fb19fc2..82a2b4b 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -45,19 +45,8 @@ class dirichlet_tessellation: # More than 1 collision is always not the best possible triangulation TODO is that true? Yesn't for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: - new_triangles = ( - # first triangle - self.get_all_edgepairs(p, neighbour, collision_edge_p1), - # second triangle - self.get_all_edgepairs(p, neighbour, collision_edge_p2) - ) - - old_triangles = ( - # first triangle - self.get_all_edgepairs(collision_edge_p1, collision_edge_p2, p), - # second triangle - self.get_all_edgepairs(collision_edge_p1, collision_edge_p2, neighbour) - ) + new_triangles = (self.get_all_edgepairs(p, neighbour, e) for e in [collision_edge_p1, collision_edge_p2]) + old_triangles = (self.get_all_edgepairs(collision_edge_p1, collision_edge_p2, e) for e in [p, neighbour]) rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) -- GitLab From 524f14cfecb0340689988260d65e9a4d20b92e94 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 15:13:18 +0200 Subject: [PATCH 17/61] replace_valid_triangulation --- curvepy/example_algo.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 82a2b4b..f83664d 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -51,9 +51,7 @@ class dirichlet_tessellation: rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) if new_is_more_equilateral: - # dann kante col_edge_p2 - col_edge_p1 entfernen und kante p - neighbour hinzufügen - collision_edge_p1 - ... + self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) @staticmethod def _angle(edge1, edge2): @@ -61,12 +59,14 @@ class dirichlet_tessellation: - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0])))) @staticmethod - def get_all_edgepairs(A, B, C): - return [(A, B), (B, C), (A, C)] - - @staticmethod - def _rate_triangle(A, B, C): - ... + def get_all_edge_pairs(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray): + return [(p1, p2), (p2, p3), (p1, p3)] + + def replace_valid_triangulation(self, old_pair: Tuple[np.ndarray, np.ndarray], new_pair: Tuple[np.ndarray, np.ndarray]): + # dann kante col_edge_p2 - col_edge_p1 entfernen und kante p - neighbour hinzufügen + # The list comprehension will always just have a single element matching the condition + self.valid_triangulation.remove(*[x for x in self.valid_triangulation if set(*x) == {*old_pair}]) + self.valid_triangulation.append(new_pair) @staticmethod def _intersect(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray) -> bool: -- GitLab From 7b05372caec70b2425d646aa826560ecf2d6963c Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 15:13:23 +0200 Subject: [PATCH 18/61] replace_valid_triangulation --- curvepy/example_algo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index f83664d..d2cd97d 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -45,8 +45,8 @@ class dirichlet_tessellation: # More than 1 collision is always not the best possible triangulation TODO is that true? Yesn't for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: - new_triangles = (self.get_all_edgepairs(p, neighbour, e) for e in [collision_edge_p1, collision_edge_p2]) - old_triangles = (self.get_all_edgepairs(collision_edge_p1, collision_edge_p2, e) for e in [p, neighbour]) + new_triangles = (self.get_all_edge_pairs(p, neighbour, e) for e in [collision_edge_p1, collision_edge_p2]) + old_triangles = (self.get_all_edge_pairs(collision_edge_p1, collision_edge_p2, e) for e in [p, neighbour]) rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) -- GitLab From fd2d1948cf5de4858b153f4182e544914354f350 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 15:24:40 +0200 Subject: [PATCH 19/61] Commit bevor wir die tiles Liste endlich zu einem dict machen --- curvepy/example_algo.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index d2cd97d..f3572a6 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -37,8 +37,9 @@ class dirichlet_tessellation: all_collisions = [x for x in self.valid_triangulation if self._intersect(neighbour.center, p, *x)] if len(all_collisions) == 0: self.valid_triangulation.append((p, neighbour.center)) - neighbour.neighbours.append( - p_tile) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind + # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind + neighbour.neighbours.append(p_tile) + p_tile.neighbours.append(neighbour) elif len(all_collisions) == 1: # 1 collision could be still fine collisions_to_check.append([p, neighbour.center, *all_collisions[0]]) @@ -52,6 +53,8 @@ class dirichlet_tessellation: new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) if new_is_more_equilateral: self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) + self.update_neighbour() + @staticmethod def _angle(edge1, edge2): @@ -62,6 +65,9 @@ class dirichlet_tessellation: def get_all_edge_pairs(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray): return [(p1, p2), (p2, p3), (p1, p3)] + def update_neighbour(self, p_new, p_neighbour, p_out_1, p_out_2): + ... + def replace_valid_triangulation(self, old_pair: Tuple[np.ndarray, np.ndarray], new_pair: Tuple[np.ndarray, np.ndarray]): # dann kante col_edge_p2 - col_edge_p1 entfernen und kante p - neighbour hinzufügen # The list comprehension will always just have a single element matching the condition -- GitLab From cdbc9d4f02ca9260cf74728bca08769b01ca8471 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 15:36:32 +0200 Subject: [PATCH 20/61] Commit nachdem das tile dict entstanden ist --- curvepy/example_algo.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index f3572a6..96eeaf6 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -9,6 +9,7 @@ class tile: neighbours: list # of tiles + class dirichlet_tessellation: def __init__(self): @@ -20,12 +21,9 @@ class dirichlet_tessellation: self.tiles.append(tile(p, [])) return - min_idx = np.argmin([np.linalg.norm(p - t.center) for t in self.tiles]) - nearest_tile = self.tiles[min_idx] - p_tile = tile(center=p, neighbours=[]) - nearest_tile.neighbours.append(p_tile) - p_tile.neighbours.append(nearest_tile) - self.tiles.append(p_tile) + nearest_p = min(self.tiles2.keys(), key=lambda t: np.linalg.norm(t - p)) + self.tiles2[p] = [nearest_p] + self.tiles2[nearest_p].append(p) # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_tile self.valid_triangulation.append((p, nearest_tile.center)) @@ -55,7 +53,6 @@ class dirichlet_tessellation: self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) self.update_neighbour() - @staticmethod def _angle(edge1, edge2): return abs(180 / np.pi * (np.arctan((edge1[1][1] - edge1[0][1]) / (edge1[1][0] - edge1[0][0])) @@ -68,7 +65,8 @@ class dirichlet_tessellation: def update_neighbour(self, p_new, p_neighbour, p_out_1, p_out_2): ... - def replace_valid_triangulation(self, old_pair: Tuple[np.ndarray, np.ndarray], new_pair: Tuple[np.ndarray, np.ndarray]): + def replace_valid_triangulation(self, old_pair: Tuple[np.ndarray, np.ndarray], + new_pair: Tuple[np.ndarray, np.ndarray]): # dann kante col_edge_p2 - col_edge_p1 entfernen und kante p - neighbour hinzufügen # The list comprehension will always just have a single element matching the condition self.valid_triangulation.remove(*[x for x in self.valid_triangulation if set(*x) == {*old_pair}]) -- GitLab From d19063312195440627fac59e75215a5b7bb5d11b Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 15:36:37 +0200 Subject: [PATCH 21/61] Commit nachdem das tile dict entstanden ist --- curvepy/example_algo.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 96eeaf6..0e7ea4f 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -1,24 +1,16 @@ -from dataclasses import dataclass import numpy as np from typing import List, Set, Union, Dict, Tuple -@dataclass -class tile: - center: np.ndarray # 2D float - neighbours: list # of tiles - - - class dirichlet_tessellation: def __init__(self): - self.tiles: List[tile] = [] + self.tiles2: Dict[np.ndarray, List[np.ndarray]] = {} self.valid_triangulation: List[Tuple[np.ndarray, np.ndarray]] = [] def append_point(self, p: np.ndarray): - if not self.tiles: - self.tiles.append(tile(p, [])) + if not self.tiles2: + self.tiles2[p] = [] return nearest_p = min(self.tiles2.keys(), key=lambda t: np.linalg.norm(t - p)) @@ -36,8 +28,8 @@ class dirichlet_tessellation: if len(all_collisions) == 0: self.valid_triangulation.append((p, neighbour.center)) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind - neighbour.neighbours.append(p_tile) - p_tile.neighbours.append(neighbour) + self.tiles2[neighbour].append(p) + self.tiles2[p].append(neighbour) elif len(all_collisions) == 1: # 1 collision could be still fine collisions_to_check.append([p, neighbour.center, *all_collisions[0]]) -- GitLab From 09cce674eb74dbd4a27cd1c7b0e6ea24e4e629d2 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 15:36:47 +0200 Subject: [PATCH 22/61] Commit nachdem das tile dict entstanden ist --- curvepy/example_algo.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 0e7ea4f..18eef55 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -17,23 +17,23 @@ class dirichlet_tessellation: self.tiles2[p] = [nearest_p] self.tiles2[nearest_p].append(p) - # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_tile - self.valid_triangulation.append((p, nearest_tile.center)) + # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_p + self.valid_triangulation.append((p, nearest_p)) # To make a more educated guess we solve any collisions __after__ adding the trivial connections collisions_to_check = [] - for neighbour in nearest_tile.neighbours: - all_collisions = [x for x in self.valid_triangulation if self._intersect(neighbour.center, p, *x)] + for neighbour in self.tiles2[nearest_p]: + all_collisions = [x for x in self.valid_triangulation if self._intersect(neighbour, p, *x)] if len(all_collisions) == 0: - self.valid_triangulation.append((p, neighbour.center)) + self.valid_triangulation.append((p, neighbour)) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind self.tiles2[neighbour].append(p) self.tiles2[p].append(neighbour) elif len(all_collisions) == 1: # 1 collision could be still fine - collisions_to_check.append([p, neighbour.center, *all_collisions[0]]) - # More than 1 collision is always not the best possible triangulation TODO is that true? Yesn't + collisions_to_check.append([p, neighbour, *all_collisions[0]]) + # More than 1 collision is always not the best possible triangulation TODO Is that true? Yesn't for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: new_triangles = (self.get_all_edge_pairs(p, neighbour, e) for e in [collision_edge_p1, collision_edge_p2]) -- GitLab From d8f78a03ddf63f8e9a22275b1311a273411a224c Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 16:24:15 +0200 Subject: [PATCH 23/61] Warum sind np.ndarray unhashable... --- curvepy/example_algo.py | 68 ++++++++++++++++++++++------------------- 1 file changed, 37 insertions(+), 31 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 18eef55..0faa055 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -5,45 +5,51 @@ from typing import List, Set, Union, Dict, Tuple class dirichlet_tessellation: def __init__(self): - self.tiles2: Dict[np.ndarray, List[np.ndarray]] = {} - self.valid_triangulation: List[Tuple[np.ndarray, np.ndarray]] = [] + self.tiles: Dict[np.ndarray, Set[np.ndarray]] = {} + self.valid_triangulation: Set[Tuple[np.ndarray, np.ndarray]] = set() def append_point(self, p: np.ndarray): - if not self.tiles2: - self.tiles2[p] = [] + if not self.tiles: + self.tiles[p] = set() return - nearest_p = min(self.tiles2.keys(), key=lambda t: np.linalg.norm(t - p)) - self.tiles2[p] = [nearest_p] - self.tiles2[nearest_p].append(p) + nearest_p = min(self.tiles.keys(), key=lambda t: np.linalg.norm(t - p)) + self.tiles[p] = {nearest_p} + self.tiles[nearest_p].add(p) # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_p - self.valid_triangulation.append((p, nearest_p)) + self.valid_triangulation.add((p, nearest_p)) # To make a more educated guess we solve any collisions __after__ adding the trivial connections - collisions_to_check = [] - - for neighbour in self.tiles2[nearest_p]: - all_collisions = [x for x in self.valid_triangulation if self._intersect(neighbour, p, *x)] - if len(all_collisions) == 0: - self.valid_triangulation.append((p, neighbour)) - # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind - self.tiles2[neighbour].append(p) - self.tiles2[p].append(neighbour) - elif len(all_collisions) == 1: - # 1 collision could be still fine - collisions_to_check.append([p, neighbour, *all_collisions[0]]) - # More than 1 collision is always not the best possible triangulation TODO Is that true? Yesn't - - for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: - new_triangles = (self.get_all_edge_pairs(p, neighbour, e) for e in [collision_edge_p1, collision_edge_p2]) - old_triangles = (self.get_all_edge_pairs(collision_edge_p1, collision_edge_p2, e) for e in [p, neighbour]) - - rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) - new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) - if new_is_more_equilateral: - self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) - self.update_neighbour() + while True: + collisions_to_check = [] + + for neighbour in self.tiles[nearest_p]: + all_collisions = [x for x in self.valid_triangulation if self._intersect(neighbour, p, *x)] + if not all_collisions: + self.valid_triangulation.add((p, neighbour)) + # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind + self.tiles[neighbour].add(p) + self.tiles[p].add(neighbour) + elif len(all_collisions) == 1: + # 1 collision could be still fine + collisions_to_check.append([p, neighbour, *all_collisions[0]]) + # More than 1 collision is always not the best possible triangulation TODO Is that true? Probably not ;) + + if not collisions_to_check: + return + + for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: + new_triangles = (self.get_all_edge_pairs(p, neighbour, e) for e in + [collision_edge_p1, collision_edge_p2]) + old_triangles = (self.get_all_edge_pairs(collision_edge_p1, collision_edge_p2, e) for e in + [p, neighbour]) + + rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) + new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) + if new_is_more_equilateral: + self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) + self.update_neighbour(p, neighbour, collision_edge_p1, collision_edge_p2) @staticmethod def _angle(edge1, edge2): -- GitLab From e62c2692285f6f07e59439aa58879cb25f68c792 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 16:24:19 +0200 Subject: [PATCH 24/61] Warum sind np.ndarray unhashable... --- curvepy/example_algo.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 0faa055..32e12f7 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -1,3 +1,4 @@ +import pprint import numpy as np from typing import List, Set, Union, Dict, Tuple @@ -60,15 +61,18 @@ class dirichlet_tessellation: def get_all_edge_pairs(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray): return [(p1, p2), (p2, p3), (p1, p3)] - def update_neighbour(self, p_new, p_neighbour, p_out_1, p_out_2): - ... + def update_neighbour(self, p_new, p_neighbour, p_rem_1, p_rem_2): + self.tiles[p_rem_1].remove(p_rem_2) + self.tiles[p_rem_2].remove(p_rem_1) + self.tiles[p_new].add(p_neighbour) + self.tiles[p_neighbour].add(p_new) def replace_valid_triangulation(self, old_pair: Tuple[np.ndarray, np.ndarray], new_pair: Tuple[np.ndarray, np.ndarray]): # dann kante col_edge_p2 - col_edge_p1 entfernen und kante p - neighbour hinzufügen # The list comprehension will always just have a single element matching the condition self.valid_triangulation.remove(*[x for x in self.valid_triangulation if set(*x) == {*old_pair}]) - self.valid_triangulation.append(new_pair) + self.valid_triangulation.add(new_pair) @staticmethod def _intersect(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray) -> bool: @@ -84,3 +88,13 @@ class dirichlet_tessellation: return False # we divide with z to turn back to 2D space return True + + +if __name__ == '__main__': + pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] + d = dirichlet_tessellation() + + for p in pts: + d.append_point(p) + + pprint.pprint(d.valid_triangulation) -- GitLab From c07be16b6616f0fe4829e94891f8cd050fc6f159 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 17:16:50 +0200 Subject: [PATCH 25/61] Wir sind doof (vor allem lars) --- curvepy/example_algo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 32e12f7..be38a5f 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -67,15 +67,15 @@ class dirichlet_tessellation: self.tiles[p_new].add(p_neighbour) self.tiles[p_neighbour].add(p_new) - def replace_valid_triangulation(self, old_pair: Tuple[np.ndarray, np.ndarray], - new_pair: Tuple[np.ndarray, np.ndarray]): + def replace_valid_triangulation(self, old_pair: Tuple[h_tuple, h_tuple], + new_pair: Tuple[h_tuple, h_tuple]): # dann kante col_edge_p2 - col_edge_p1 entfernen und kante p - neighbour hinzufügen # The list comprehension will always just have a single element matching the condition self.valid_triangulation.remove(*[x for x in self.valid_triangulation if set(*x) == {*old_pair}]) self.valid_triangulation.add(new_pair) @staticmethod - def _intersect(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray) -> bool: + def _intersect(p1: h_tuple, p2: h_tuple, p3: h_tuple, p4: h_tuple) -> bool: # First we vertical stack the points in an array vertical_stack = np.vstack([p1, p2, p3, p4]) # Then we transform them to homogeneous coordinates, to perform a little trick -- GitLab From 2034800b0db4ef735edbff005c64907f72cc6390 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 17:16:54 +0200 Subject: [PATCH 26/61] Wir sind doof (vor allem lars) --- curvepy/example_algo.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index be38a5f..1f35a73 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -1,15 +1,26 @@ import pprint import numpy as np -from typing import List, Set, Union, Dict, Tuple +from typing import List, Set, Dict, Tuple + +class h_tuple: + def __init__(self, xs): + self.data: np.ndarray = xs + + def __hash__(self): + return hash((self.data[0], self.data[1])) + + def __iter__(self): + return iter(self.data) class dirichlet_tessellation: def __init__(self): - self.tiles: Dict[np.ndarray, Set[np.ndarray]] = {} - self.valid_triangulation: Set[Tuple[np.ndarray, np.ndarray]] = set() + self.tiles: Dict[h_tuple, Set[h_tuple]] = {} + self.valid_triangulation: Set[Tuple[h_tuple, h_tuple]] = set() def append_point(self, p: np.ndarray): + p = h_tuple(p) if not self.tiles: self.tiles[p] = set() return @@ -61,7 +72,7 @@ class dirichlet_tessellation: def get_all_edge_pairs(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray): return [(p1, p2), (p2, p3), (p1, p3)] - def update_neighbour(self, p_new, p_neighbour, p_rem_1, p_rem_2): + def update_neighbour(self, p_new: h_tuple, p_neighbour: h_tuple, p_rem_1: h_tuple, p_rem_2: h_tuple): self.tiles[p_rem_1].remove(p_rem_2) self.tiles[p_rem_2].remove(p_rem_1) self.tiles[p_new].add(p_neighbour) @@ -77,7 +88,7 @@ class dirichlet_tessellation: @staticmethod def _intersect(p1: h_tuple, p2: h_tuple, p3: h_tuple, p4: h_tuple) -> bool: # First we vertical stack the points in an array - vertical_stack = np.vstack([p1, p2, p3, p4]) + vertical_stack = np.vstack([p1.data, p2.data, p3.data, p4.data]) # Then we transform them to homogeneous coordinates, to perform a little trick homogeneous = np.hstack((vertical_stack, np.ones((4, 1)))) # having our points in this form we can get the lines through the cross product @@ -98,3 +109,6 @@ if __name__ == '__main__': d.append_point(p) pprint.pprint(d.valid_triangulation) + + hash = lambda n: hash((*n)) + np.ndarray -- GitLab From bcf90e5283f21f7166794a81366eb99c0316e70e Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 22 May 2021 17:16:58 +0200 Subject: [PATCH 27/61] Wir sind doof (vor allem lars) --- curvepy/example_algo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 1f35a73..b95a645 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -25,8 +25,8 @@ class dirichlet_tessellation: self.tiles[p] = set() return - nearest_p = min(self.tiles.keys(), key=lambda t: np.linalg.norm(t - p)) - self.tiles[p] = {nearest_p} + nearest_p = min(self.tiles.keys(), key=lambda t: np.linalg.norm(t.data - p.data)) + self.tiles[h_tuple(p)] = {nearest_p} self.tiles[nearest_p].add(p) # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_p @@ -69,7 +69,7 @@ class dirichlet_tessellation: - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0])))) @staticmethod - def get_all_edge_pairs(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray): + def get_all_edge_pairs(p1: h_tuple, p2: h_tuple, p3: h_tuple): return [(p1, p2), (p2, p3), (p1, p3)] def update_neighbour(self, p_new: h_tuple, p_neighbour: h_tuple, p_rem_1: h_tuple, p_rem_2: h_tuple): -- GitLab From f71ad3640f049dad5d57e1e78ac72b5e72ea947d Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 15:22:29 +0200 Subject: [PATCH 28/61] h_tuple funzen --- curvepy/example_algo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index b95a645..c66050f 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -108,7 +108,7 @@ if __name__ == '__main__': for p in pts: d.append_point(p) - pprint.pprint(d.valid_triangulation) + print(d.valid_triangulation) hash = lambda n: hash((*n)) np.ndarray -- GitLab From 18719c126abc616faf8a2780dcc029a6318a7b65 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 15:22:46 +0200 Subject: [PATCH 29/61] h_tuple funzen --- curvepy/example_algo.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index c66050f..1072c6c 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -1,16 +1,15 @@ import pprint import numpy as np -from typing import List, Set, Dict, Tuple +from typing import List, Set, Dict, Tuple, Iterable -class h_tuple: - def __init__(self, xs): - self.data: np.ndarray = xs +class h_tuple(np.ndarray): def __hash__(self): - return hash((self.data[0], self.data[1])) + return hash(tuple(self)) - def __iter__(self): - return iter(self.data) + @staticmethod + def create(xs: Iterable): + return np.array(xs).view(h_tuple) class dirichlet_tessellation: @@ -41,6 +40,8 @@ class dirichlet_tessellation: if not all_collisions: self.valid_triangulation.add((p, neighbour)) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind + print(self.tiles) + print(self.tiles.get(neighbour) is None) self.tiles[neighbour].add(p) self.tiles[p].add(neighbour) elif len(all_collisions) == 1: @@ -88,7 +89,7 @@ class dirichlet_tessellation: @staticmethod def _intersect(p1: h_tuple, p2: h_tuple, p3: h_tuple, p4: h_tuple) -> bool: # First we vertical stack the points in an array - vertical_stack = np.vstack([p1.data, p2.data, p3.data, p4.data]) + vertical_stack = np.vstack([p1, p2, p3, p4]) # Then we transform them to homogeneous coordinates, to perform a little trick homogeneous = np.hstack((vertical_stack, np.ones((4, 1)))) # having our points in this form we can get the lines through the cross product -- GitLab From 441c6bd63de37276e0d1159849392af67e734d1f Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 15:22:50 +0200 Subject: [PATCH 30/61] Wir sind doof (vor allem lars) --- curvepy/example_algo.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 1072c6c..b7b510f 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -19,13 +19,13 @@ class dirichlet_tessellation: self.valid_triangulation: Set[Tuple[h_tuple, h_tuple]] = set() def append_point(self, p: np.ndarray): - p = h_tuple(p) + p = h_tuple.create(p) if not self.tiles: self.tiles[p] = set() return - nearest_p = min(self.tiles.keys(), key=lambda t: np.linalg.norm(t.data - p.data)) - self.tiles[h_tuple(p)] = {nearest_p} + nearest_p = min(self.tiles.keys(), key=lambda t: np.linalg.norm(t - p)) + self.tiles[p] = {nearest_p} self.tiles[nearest_p].add(p) # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_p @@ -36,7 +36,11 @@ class dirichlet_tessellation: collisions_to_check = [] for neighbour in self.tiles[nearest_p]: - all_collisions = [x for x in self.valid_triangulation if self._intersect(neighbour, p, *x)] + print(p) + print(neighbour) + if all(x == y for x, y in zip(p, neighbour)): + continue + all_collisions = [x for x in self.valid_triangulation if self.intersect(neighbour, p, *x)] if not all_collisions: self.valid_triangulation.add((p, neighbour)) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind @@ -87,7 +91,7 @@ class dirichlet_tessellation: self.valid_triangulation.add(new_pair) @staticmethod - def _intersect(p1: h_tuple, p2: h_tuple, p3: h_tuple, p4: h_tuple) -> bool: + def intersect(p1: h_tuple, p2: h_tuple, p3: h_tuple, p4: h_tuple) -> bool: # First we vertical stack the points in an array vertical_stack = np.vstack([p1, p2, p3, p4]) # Then we transform them to homogeneous coordinates, to perform a little trick @@ -103,6 +107,8 @@ class dirichlet_tessellation: if __name__ == '__main__': + #xs = [h_tuple.create(x) for x in ((1,2), (3,4), (1,4), (3,6)) ] + #print(dirichlet_tessellation.intersect(*xs)) pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] d = dirichlet_tessellation() @@ -111,5 +117,5 @@ if __name__ == '__main__': print(d.valid_triangulation) - hash = lambda n: hash((*n)) - np.ndarray + # hash = lambda n: hash(()) + # np.ndarray -- GitLab From 81450632fe9f68bc1f18464095552b6a7f82281c Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 16:11:46 +0200 Subject: [PATCH 31/61] no errors but a loop (no homo) --- curvepy/example_algo.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index b7b510f..585c00c 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -103,15 +103,16 @@ class dirichlet_tessellation: if z == 0: return False # we divide with z to turn back to 2D space - return True if __name__ == '__main__': - #xs = [h_tuple.create(x) for x in ((1,2), (3,4), (1,4), (3,6)) ] - #print(dirichlet_tessellation.intersect(*xs)) + # xs = [h_tuple.create(x) for x in ((1,2), (3,4), (1,4), (3,6)) ] + # print(dirichlet_tessellation.intersect(*xs)) pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] d = dirichlet_tessellation() + # print({2.0 == 2.}) + for p in pts: d.append_point(p) -- GitLab From 0fe7a0421fc802eb1078e51a7f4b095130392bc6 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 16:11:50 +0200 Subject: [PATCH 32/61] no errors but a loop (no homo) --- curvepy/example_algo.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 585c00c..eda9d91 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -9,7 +9,8 @@ class h_tuple(np.ndarray): @staticmethod def create(xs: Iterable): - return np.array(xs).view(h_tuple) + print("xd",xs) + return np.array(xs, dtype=np.float32).view(h_tuple) class dirichlet_tessellation: @@ -40,7 +41,9 @@ class dirichlet_tessellation: print(neighbour) if all(x == y for x, y in zip(p, neighbour)): continue + print("b4") all_collisions = [x for x in self.valid_triangulation if self.intersect(neighbour, p, *x)] + print("after (no homo)") if not all_collisions: self.valid_triangulation.add((p, neighbour)) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind @@ -87,7 +90,7 @@ class dirichlet_tessellation: new_pair: Tuple[h_tuple, h_tuple]): # dann kante col_edge_p2 - col_edge_p1 entfernen und kante p - neighbour hinzufügen # The list comprehension will always just have a single element matching the condition - self.valid_triangulation.remove(*[x for x in self.valid_triangulation if set(*x) == {*old_pair}]) + self.valid_triangulation.remove(*[x for x in self.valid_triangulation if set(x) == set(old_pair)]) self.valid_triangulation.add(new_pair) @staticmethod @@ -100,8 +103,9 @@ class dirichlet_tessellation: line_1, line_2 = np.cross(homogeneous[0], homogeneous[1]), np.cross(homogeneous[2], homogeneous[3]) # when we calculate the cross product of the lines we get intersect point x, y, z = np.cross(line_1, line_2) - if z == 0: - return False + print(f"{h_tuple.create([x/z, y/z])} und die Punkte sind: {p1}, {p2}, {p3}, {p4}") + # print([(h_tuple.create([x/z, y/z]), x) for x in [p1, p2, p3, p4]]) + return not (z == 0 or any([x/z == p[0] and y/z == p[1] for p in [p1, p2, p3, p4]])) # we divide with z to turn back to 2D space -- GitLab From 91b95a3b5b9ea44ecac88cb2fa1aeb2418424297 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 16:11:54 +0200 Subject: [PATCH 33/61] no errors but a loop (no homo) --- curvepy/example_algo.py | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index eda9d91..5e941fd 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -37,8 +37,6 @@ class dirichlet_tessellation: collisions_to_check = [] for neighbour in self.tiles[nearest_p]: - print(p) - print(neighbour) if all(x == y for x, y in zip(p, neighbour)): continue print("b4") @@ -47,8 +45,6 @@ class dirichlet_tessellation: if not all_collisions: self.valid_triangulation.add((p, neighbour)) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind - print(self.tiles) - print(self.tiles.get(neighbour) is None) self.tiles[neighbour].add(p) self.tiles[p].add(neighbour) elif len(all_collisions) == 1: @@ -60,10 +56,27 @@ class dirichlet_tessellation: return for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: - new_triangles = (self.get_all_edge_pairs(p, neighbour, e) for e in - [collision_edge_p1, collision_edge_p2]) - old_triangles = (self.get_all_edge_pairs(collision_edge_p1, collision_edge_p2, e) for e in - [p, neighbour]) + new_triangles = ( + # first triangle + ((p, neighbour), (p, collision_edge_p1)), + ((p, neighbour), (neighbour, collision_edge_p1)), + ((p, neighbour), (neighbour, collision_edge_p1)), + # second triangle + ((p, neighbour), (p, collision_edge_p2)), + ((p, collision_edge_p2), (neighbour, collision_edge_p2)), + ((p, neighbour), (neighbour, collision_edge_p2)) + ) + + old_triangles = ( + # first triangle + ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, p)), + ((collision_edge_p1, p), (collision_edge_p2, p)), + ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, p)), + # second triangle + ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, neighbour)), + ((collision_edge_p1, neighbour), (collision_edge_p2, neighbour)), + ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, neighbour)) + ) rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) -- GitLab From 245d59105474a1708d3f5e0cbd5c70fb0a9eee72 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 17:57:58 +0200 Subject: [PATCH 34/61] eine kannte falsch und error in _angle --- curvepy/example_algo.py | 40 ++++++++++++++++++++++++++++++++-------- 1 file changed, 32 insertions(+), 8 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 5e941fd..d51647f 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -9,7 +9,6 @@ class h_tuple(np.ndarray): @staticmethod def create(xs: Iterable): - print("xd",xs) return np.array(xs, dtype=np.float32).view(h_tuple) @@ -32,6 +31,8 @@ class dirichlet_tessellation: # ERSTELLE TRIANGULATION ZWISCHEN p UND nearest_p self.valid_triangulation.add((p, nearest_p)) + garbage_heap = [] + # To make a more educated guess we solve any collisions __after__ adding the trivial connections while True: collisions_to_check = [] @@ -39,10 +40,11 @@ class dirichlet_tessellation: for neighbour in self.tiles[nearest_p]: if all(x == y for x, y in zip(p, neighbour)): continue - print("b4") - all_collisions = [x for x in self.valid_triangulation if self.intersect(neighbour, p, *x)] - print("after (no homo)") - if not all_collisions: + # print("b4") + all_collisions = [x for x in self.valid_triangulation if self.intersect(neighbour, p, *x) + and x not in garbage_heap and {neighbour, p} not in garbage_heap] + # print("after (no homo)") + if (not all_collisions) and ({p, neighbour} not in garbage_heap): self.valid_triangulation.add((p, neighbour)) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind self.tiles[neighbour].add(p) @@ -83,6 +85,9 @@ class dirichlet_tessellation: if new_is_more_equilateral: self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) self.update_neighbour(p, neighbour, collision_edge_p1, collision_edge_p2) + garbage_heap.append({collision_edge_p1, collision_edge_p2}) + else: + garbage_heap.append({p, neighbour}) @staticmethod def _angle(edge1, edge2): @@ -116,11 +121,24 @@ class dirichlet_tessellation: line_1, line_2 = np.cross(homogeneous[0], homogeneous[1]), np.cross(homogeneous[2], homogeneous[3]) # when we calculate the cross product of the lines we get intersect point x, y, z = np.cross(line_1, line_2) - print(f"{h_tuple.create([x/z, y/z])} und die Punkte sind: {p1}, {p2}, {p3}, {p4}") + # print(f"{h_tuple.create([x/z, y/z])} und die Punkte sind: {p1}, {p2}, {p3}, {p4}") # print([(h_tuple.create([x/z, y/z]), x) for x in [p1, p2, p3, p4]]) - return not (z == 0 or any([x/z == p[0] and y/z == p[1] for p in [p1, p2, p3, p4]])) # we divide with z to turn back to 2D space + # no intersection + if z == 0: + return False + intersection = h_tuple.create([x / z, y / z]) + # if intersection is at one of the points + if any([intersection == p for p in [p1, p2, p3, p4]]): + return False + # if the intersection is in the convex space of both points we good (aka not on the line segment) + return dirichlet_tessellation.intersection_is_between(intersection, p1, p2) + + @staticmethod + def intersection_is_between(intersection: h_tuple, p1: h_tuple, p2: h_tuple): + return abs(intersection - p1) + abs(intersection - p2) == abs(p2 - p1) + if __name__ == '__main__': # xs = [h_tuple.create(x) for x in ((1,2), (3,4), (1,4), (3,6)) ] @@ -128,7 +146,13 @@ if __name__ == '__main__': pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] d = dirichlet_tessellation() - # print({2.0 == 2.}) + # print(h_tuple.create([2,1]) == h_tuple.create([1, 2])) + # cords A, B, D, E + # Winkel altes dreieck1 = 108.44, 26.57, 45 + # Winkel altes dreieck2 = 33.69, 112.62, 33.69 + + # Winkel neues dreieck1 = 49.4, 70.35, 60.26 + # Winkel neues dreieck2 = 59.04, 78.69, 42.27 for p in pts: d.append_point(p) -- GitLab From 1631af5a2717201d7bd283edfd5d51dfe89886be Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 17:58:04 +0200 Subject: [PATCH 35/61] eine kannte falsch und error in _angle --- curvepy/example_algo.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index d51647f..e45712b 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -7,6 +7,9 @@ class h_tuple(np.ndarray): def __hash__(self): return hash(tuple(self)) + def __eq__(self, other): + return len(self) == len(other) and all(x == y for x, y in zip(self, other)) + @staticmethod def create(xs: Iterable): return np.array(xs, dtype=np.float32).view(h_tuple) -- GitLab From 5f776564f43940bf7a3924e9ca9fb69da4e1132b Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 17:58:09 +0200 Subject: [PATCH 36/61] eine kannte falsch und error in _angle --- curvepy/example_algo.py | 1 + 1 file changed, 1 insertion(+) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index e45712b..0ba61de 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -39,6 +39,7 @@ class dirichlet_tessellation: # To make a more educated guess we solve any collisions __after__ adding the trivial connections while True: collisions_to_check = [] + # print(f"collisions_to_check: {collisions_to_check}") for neighbour in self.tiles[nearest_p]: if all(x == y for x, y in zip(p, neighbour)): -- GitLab From ca6e64f2bc6c19af81589c2c14d56a59cb5b26b4 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 30 May 2021 19:46:05 +0200 Subject: [PATCH 37/61] =?UTF-8?q?neue=20Dreiecke=20Winkelberechnung=20sieh?= =?UTF-8?q?t=20richtig=20aus,=20alten=20Dreiecke=20Winkelberechnung=20sieh?= =?UTF-8?q?t=20nicht=20richtig=20aus.=20Warum=20wei=C3=9F=20keiner.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- curvepy/example_algo.py | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 0ba61de..b647bb4 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -61,11 +61,13 @@ class dirichlet_tessellation: if not collisions_to_check: return + # Kantenpaare für rate_tri for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: new_triangles = ( # first triangle ((p, neighbour), (p, collision_edge_p1)), - ((p, neighbour), (neighbour, collision_edge_p1)), + # ((p, neighbour), (neighbour, collision_edge_p1)), <- HEUTE REINKOPIERT + ((p, collision_edge_p1), (neighbour, collision_edge_p1)), # maybe jetzt die richtige kante? ((p, neighbour), (neighbour, collision_edge_p1)), # second triangle ((p, neighbour), (p, collision_edge_p2)), @@ -85,18 +87,37 @@ class dirichlet_tessellation: ) rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) - new_is_more_equilateral = rate_tri(old_triangles) > rate_tri(new_triangles) + a = rate_tri(old_triangles) + b = rate_tri(new_triangles) + new_is_more_equilateral = a > b + print(f"rate_tri(old_triangles) = {a}") + print(f"rate_tri(new_triangles) = {b}") if new_is_more_equilateral: self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) self.update_neighbour(p, neighbour, collision_edge_p1, collision_edge_p2) garbage_heap.append({collision_edge_p1, collision_edge_p2}) else: garbage_heap.append({p, neighbour}) + print(f"garbage_heap = {garbage_heap}") @staticmethod def _angle(edge1, edge2): - return abs(180 / np.pi * (np.arctan((edge1[1][1] - edge1[0][1]) / (edge1[1][0] - edge1[0][0])) - - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0])))) + print(f"edge1 = {edge1}, edge2 = {edge2}") + # print(f"(edge1[1][0] - edge1[0][0]) = {(edge1[1][0] - edge1[0][0])}") + # print(f"(edge2[1][0] - edge2[0][0]) = {(edge2[1][0] - edge2[0][0])}") + # print(f"_angle = {abs(180 / np.pi * (np.arctan((edge1[1][1] - edge1[0][1]) / (edge1[1][0] - edge1[0][0])) - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0]))))}") + slope1 = dirichlet_tessellation.slope(edge1) + slope2 = dirichlet_tessellation.slope(edge2) + print(f"(180 / np.pi) * np.arctan(abs((slope2 - slope1)/(1 + slope1*slope2))) = {(180 / np.pi) * np.arctan(abs((slope2 - slope1)/(1 + slope1*slope2)))}") + return (180 / np.pi) * np.arctan(abs((slope2 - slope1)/(1 + slope1*slope2))) + # return abs(180 / np.pi * (np.arctan((edge1[1][1] - edge1[0][1]) / (edge1[1][0] - edge1[0][0])) + # - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0])))) + + @staticmethod + def slope(edge): + if edge[0][1] == edge[1][1]: + return 0 + return (edge[0][0] - edge[1][0]) / (edge[0][1] - edge[1][1]) @staticmethod def get_all_edge_pairs(p1: h_tuple, p2: h_tuple, p3: h_tuple): @@ -160,6 +181,7 @@ if __name__ == '__main__': for p in pts: d.append_point(p) + print(f"d.valid_triangulation = {d.valid_triangulation}") print(d.valid_triangulation) -- GitLab From 387a721e1d1205564132bdbc9a466e4ecb4db8e8 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 6 Jun 2021 15:12:02 +0200 Subject: [PATCH 38/61] -delauny triangulation mit kreisen -bevor numpy arrays in den Ofen geschossen werden --- curvepy/testrange.py | 194 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 192 insertions(+), 2 deletions(-) diff --git a/curvepy/testrange.py b/curvepy/testrange.py index 80af607..aaa53ad 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -1,3 +1,193 @@ +import numpy as np +from typing import List +import itertools -x = list(itertools.chain(*([(A,B), (B,C), (A,C)] for A,B,C in [((p, neighbour), (p, fst_collision), (neighbour, fst_collision)), - ((p, neighbour), (p, snd_collision), (neighbour, snd_collision))]))) \ No newline at end of file +class Triangle: + def __init__(self, A: np.ndarray, B: np.ndarray, C: np.ndarray): + self.points = [np.array(x) for x in sorted([[*A], [*B], [*C]])] + self.area = self.calc_area(*A, *B, *C) + self.circumcircle = self.calculate_circumcircle() + + def calculate_circumcircle(self): + """ + :return: + + See: https://de.wikipedia.org/wiki/Umkreis#Koordinaten + See: https://de.wikipedia.org/wiki/Umkreis#Radius + """ + A, B, C = self.points + [x1, y1], [x2, y2], [x3, y3] = A, B, C + d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) + xu = ((x1 * x1 + y1 * y1) * (y2 - y3) + (x2 * x2 + y2 * y2) * (y3 - y1) + (x3 * x3 + y3 * y3) * (y1 - y2)) / d + yu = ((x1 * x1 + y1 * y1) * (x3 - x2) + (x2 * x2 + y2 * y2) * (x1 - x3) + (x3 * x3 + y3 * y3) * (x2 - x1)) / d + + lines = [[A, B], [B, C], [A, C]] + c, a, b = [np.linalg.norm(x - y) for x, y in lines] + + R = (a * b * c) / (4 * self.area) + return Circle(center=np.array([xu, yu]), radius=R) + + def get_farthest_point_away_and_nearest_line_to(self, pt: np.ndarray): + points = self.points.copy() + farthest_p = self.points[np.argmax([np.linalg.norm(x - pt) for x in points])] + points.remove(farthest_p) + return farthest_p, points + + def __contains__(self, pt: np.ndarray): + """ + See: https://www.geeksforgeeks.org/check-whether-a-given-point-lies-inside-a-triangle-or-not/ + """ + [x1, y1], [x2, y2], [x3, y3] = self.points + x, y = pt + area1 = self.calc_area(x, y, x2, y2, x3, y3) + area2 = self.calc_area(x1, y1, x, y, x3, y3) + area3 = self.calc_area(x1, y1, x2, y2, x, y) + return self.area == area1 + area2 + area3 + + @staticmethod + def calc_area(x1, y1, x2, y2, x3, y3): + """ + See: https://www.mathopenref.com/coordtrianglearea.html + """ + return abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2.0) + + def check_if_point_is_in_points(self, pt): + for p in self.points: + if [*p] == [*pt]: + return True + return False + + +class Circle: + def __init__(self, center: np.ndarray, radius: float): + self.center = center + self.radius = radius + + def point_in_me(self, pt: np.ndarray) -> bool: + return np.linalg.norm(pt - self.center) <= self.radius + + +class Delaunay_triangulation: + def __init__(self, xmin, xmax, ymin, ymax): + self.xmin, self.xmax, self.ymin, self.ymax = xmin, xmax, ymin, ymax + self.supertriangle = self.create_supertriangle() + self._triangles: List[Triangle] = [self.supertriangle] + self.triangle_queue = [] + + @property + def triangles(self): + # We have to remove everything containing vertices of the supertriangle + rem_if = lambda t: any(pt in t.points for pt in self.supertriangle.points) + return [t for t in self._triangles if rem_if(t)] + + def get_all_points(self): + return np.unique(np.array(sum([t.points for t in self._triangles], [])), axis=0) + + def add_point(self, pt: np.ndarray): + if not self.is_in_range(pt): + raise AssertionError("point not in predefined range") + t = self.find_triangle_containing(pt) + self.add_new_triangles(pt, t) + self._triangles.remove(t) + # new _triangles + self.triangle_queue = self._triangles[-3:] + pts = self.get_all_points() + while self.triangle_queue: + t = self.triangle_queue.pop() + for p in pts: + if t.check_if_point_is_in_points(p): + continue + if t.circumcircle.point_in_me(p): + self.handle_point_in_circumcircle(t, p) + break + + def handle_point_in_circumcircle(self, current_t, p_in_circumcircle): + farthest_pt, nearest_pts = current_t.get_farthest_point_away_and_nearest_line_to(p_in_circumcircle) + t1 = Triangle(p_in_circumcircle, farthest_pt, nearest_pts[0]) + t2 = Triangle(p_in_circumcircle, farthest_pt, nearest_pts[1]) + self.update_triangle_structures(t1, t2, nearest_pts) + + def update_triangle_structures(self, t1_new, t2_new, removed_line): + triangles_to_remove = [t for t in self._triangles if + removed_line[0] in t.points and removed_line[1] in t.points] + for x in triangles_to_remove: + self._triangles.remove(x) + self._triangles += [t1_new, t2_new] + self.triangle_queue += [t1_new, t2_new] + + def create_supertriangle(self) -> Triangle: + # Rectangle CCW + A = np.array([self.xmin, self.ymax]) + B = np.array([self.xmin, self.ymin]) + C = np.array([self.xmax, self.ymin]) + D = np.array([self.xmax, self.ymax]) + + # Those _triangles are chosen "randomly" + # They would look like + # lower_left_triangle = [A, B, A_supert] + # lower_right_triangle = [D, B_supert, C] + A_supert = np.array([self.xmin - (self.xmax - self.xmin) / 2, self.ymin]) + B_supert = np.array([self.xmax + (self.xmax - self.xmin) / 2, self.ymin]) + C_supert = self.intersect_lines(A_supert, A, B_supert, D) + + return Triangle(A_supert, B_supert, C_supert) + + def intersect_lines(self, p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray): + """ + Method checking if line through p1, p2 intersects with line through p3, p4 + + + Parameters + ---------- + p1: np.ndarray: + first point of first line + + p2: np.ndarray: + second point of first line + + p3: np.ndarray: + first point of second line + + p4: np.ndarray: + second point of second line + + Returns + ------- + bool: + True if all point are less than tol away from line otherwise false + """ + # First we vertical stack the points in an array + vertical_stack = np.vstack([p1, p2, p3, p4]) + # Then we transform them to homogeneous coordinates, to perform a little trick + homogeneous = np.hstack((vertical_stack, np.ones((4, 1)))) + # having our points in this form we can get the lines through the cross product + line_1, line_2 = np.cross(homogeneous[0], homogeneous[1]), np.cross(homogeneous[2], homogeneous[3]) + # when we calculate the cross product of the lines we get intersect point + x, y, z = np.cross(line_1, line_2) + if z == 0: + return None + # we divide with z to turn back to 2D space + return np.array([x / z, y / z]) + + def is_in_range(self, pt): + return self.xmin <= pt[0] <= self.xmax and self.ymin <= pt[1] <= self.ymax + + def add_new_triangles(self, pt, t): + for a, b in itertools.product(t.points, t.points): + if [*a] != [*b]: + self._triangles.append(Triangle(a, b, pt)) + + def find_triangle_containing(self, pt) -> Triangle: + for t in self._triangles: + if pt in t: + return t + raise LookupError(f"No triangle containing ({pt})") + + +if __name__ == '__main__': + pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] + d = Delaunay_triangulation(-100, 100, -100, 100) + for p in pts: + d.add_point(p) + + print(d.triangles) -- GitLab From 3ff274e734e8680731713aea6dfcb362b25f5c2f Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 6 Jun 2021 16:16:11 +0200 Subject: [PATCH 39/61] -flipen funktioniert noch nicht --- curvepy/example_algo.py | 144 ++++++++++++++++++++++++++++++---------- 1 file changed, 110 insertions(+), 34 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index b647bb4..ff6365b 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -23,11 +23,13 @@ class dirichlet_tessellation: def append_point(self, p: np.ndarray): p = h_tuple.create(p) + print(f"p = {p}") if not self.tiles: self.tiles[p] = set() return nearest_p = min(self.tiles.keys(), key=lambda t: np.linalg.norm(t - p)) + print(f"nearest_p = {nearest_p}") self.tiles[p] = {nearest_p} self.tiles[nearest_p].add(p) @@ -44,11 +46,13 @@ class dirichlet_tessellation: for neighbour in self.tiles[nearest_p]: if all(x == y for x, y in zip(p, neighbour)): continue - # print("b4") + all_collisions = [x for x in self.valid_triangulation if self.intersect(neighbour, p, *x) and x not in garbage_heap and {neighbour, p} not in garbage_heap] - # print("after (no homo)") + # print(f"all_collisions = {all_collisions}") if (not all_collisions) and ({p, neighbour} not in garbage_heap): + # print(f"p wenn keine collisions = {p}") + # print(f"neighbour wenn keine collisions = {neighbour}") self.valid_triangulation.add((p, neighbour)) # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind self.tiles[neighbour].add(p) @@ -61,44 +65,82 @@ class dirichlet_tessellation: if not collisions_to_check: return - # Kantenpaare für rate_tri for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: - new_triangles = ( - # first triangle - ((p, neighbour), (p, collision_edge_p1)), - # ((p, neighbour), (neighbour, collision_edge_p1)), <- HEUTE REINKOPIERT - ((p, collision_edge_p1), (neighbour, collision_edge_p1)), # maybe jetzt die richtige kante? - ((p, neighbour), (neighbour, collision_edge_p1)), - # second triangle - ((p, neighbour), (p, collision_edge_p2)), - ((p, collision_edge_p2), (neighbour, collision_edge_p2)), - ((p, neighbour), (neighbour, collision_edge_p2)) + new_triangles2 = ( + (p, neighbour, collision_edge_p1), + (p, neighbour, collision_edge_p2) ) - - old_triangles = ( - # first triangle - ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, p)), - ((collision_edge_p1, p), (collision_edge_p2, p)), - ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, p)), - # second triangle - ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, neighbour)), - ((collision_edge_p1, neighbour), (collision_edge_p2, neighbour)), - ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, neighbour)) + old_triangles2 = ( + (p, collision_edge_p1, collision_edge_p2), + (neighbour, collision_edge_p1, collision_edge_p2) ) - rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) - a = rate_tri(old_triangles) - b = rate_tri(new_triangles) - new_is_more_equilateral = a > b - print(f"rate_tri(old_triangles) = {a}") - print(f"rate_tri(new_triangles) = {b}") + print(f"p = {p}") + print(f"neighbour = {neighbour}") + print(f"collision_edge_p1 = {collision_edge_p1}") + print(f"collision_edge_p2 = {collision_edge_p2}") + + angle_list_new_triangle = self._angles_in_triangle(new_triangles2[0]) + # print(f"angle_list_new_triangle = {angle_list_new_triangle}") + angle_rating_new_triangle = sum(abs(60 - angle) for angle in angle_list_new_triangle) + angle_list_new_triangle = self._angles_in_triangle(new_triangles2[1]) + angle_rating_new_triangle += sum(abs(60 - angle) for angle in angle_list_new_triangle) + # print(f"angle_rating_new_triangle = {angle_rating_new_triangle}") + + angle_list_old_triangle = self._angles_in_triangle(old_triangles2[0]) + # print(f"angle_list_old_triangle = {angle_list_old_triangle}") + angle_rating_old_triangle = sum(abs(60 - angle) for angle in angle_list_old_triangle) + angle_list_old_triangle = self._angles_in_triangle(old_triangles2[1]) + angle_rating_old_triangle += sum(abs(60 - angle) for angle in angle_list_old_triangle) + # print(f"angle_rating_new_triangle = {angle_rating_old_triangle}") + + new_is_more_equilateral = angle_rating_old_triangle > angle_rating_new_triangle + if new_is_more_equilateral: self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) self.update_neighbour(p, neighbour, collision_edge_p1, collision_edge_p2) garbage_heap.append({collision_edge_p1, collision_edge_p2}) else: garbage_heap.append({p, neighbour}) - print(f"garbage_heap = {garbage_heap}") + + # Kantenpaare für rate_tri + # for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: + # new_triangles = ( + # # first triangle + # ((p, neighbour), (p, collision_edge_p1)), + # # ((p, neighbour), (neighbour, collision_edge_p1)), <- HEUTE REINKOPIERT + # ((p, collision_edge_p1), (neighbour, collision_edge_p1)), # maybe jetzt die richtige kante? + # ((p, neighbour), (neighbour, collision_edge_p1)), + # # second triangle + # ((p, neighbour), (p, collision_edge_p2)), + # ((p, collision_edge_p2), (neighbour, collision_edge_p2)), + # ((p, neighbour), (neighbour, collision_edge_p2)) + # ) + # + # old_triangles = ( + # # first triangle + # ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, p)), + # ((collision_edge_p1, p), (collision_edge_p2, p)), + # ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, p)), + # # second triangle + # ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, neighbour)), + # ((collision_edge_p1, neighbour), (collision_edge_p2, neighbour)), + # ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, neighbour)) + # ) + # + # rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) + # a = rate_tri(old_triangles) + # b = rate_tri(new_triangles) + # new_is_more_equilateral = a > b + # print(f"rate_tri(old_triangles) = {a}") + # print(f"rate_tri(new_triangles) = {b}") + # if new_is_more_equilateral: + # self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) + # self.update_neighbour(p, neighbour, collision_edge_p1, collision_edge_p2) + # garbage_heap.append({collision_edge_p1, collision_edge_p2}) + # else: + # garbage_heap.append({p, neighbour}) + # print(f"garbage_heap = {garbage_heap}") @staticmethod def _angle(edge1, edge2): @@ -113,6 +155,21 @@ class dirichlet_tessellation: # return abs(180 / np.pi * (np.arctan((edge1[1][1] - edge1[0][1]) / (edge1[1][0] - edge1[0][0])) # - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0])))) + @staticmethod + def _angles_in_triangle(triangle): + # law of cosine + # https://en.wikipedia.org/wiki/Solution_of_triangles + + a = np.linalg.norm(triangle[1]-triangle[2]) + b = np.linalg.norm(triangle[0]-triangle[1]) + c = np.linalg.norm(triangle[2]-triangle[0]) + + alpha = (180 / np.pi) * np.arccos((b**2+c**2-a**2)/(2*b*c)) + beta = (180 / np.pi) * np.arccos((a**2+c**2-b**2)/(2*a*c)) + gamma = 180 - alpha - beta + + return [alpha, beta, gamma] + @staticmethod def slope(edge): if edge[0][1] == edge[1][1]: @@ -158,11 +215,29 @@ class dirichlet_tessellation: if any([intersection == p for p in [p1, p2, p3, p4]]): return False # if the intersection is in the convex space of both points we good (aka not on the line segment) - return dirichlet_tessellation.intersection_is_between(intersection, p1, p2) + return dirichlet_tessellation.intersection_is_between(intersection, p1, p2, p3, p4) @staticmethod - def intersection_is_between(intersection: h_tuple, p1: h_tuple, p2: h_tuple): - return abs(intersection - p1) + abs(intersection - p2) == abs(p2 - p1) + def intersection_is_between(intersection: h_tuple, p1: h_tuple, p2: h_tuple, p3, p4): + # a = np.linalg.norm(intersection - p1) + np.linalg.norm(intersection - p2) + # b = np.linalg.norm(p2 - p1) + # c = a == b + # print(c) + # return c + # return abs(intersection - p1) + abs(intersection - p2) == abs(p2 - p1) + seg1_max_x = max(p1[0], p2[0]) + seg1_min_x = min(p1[0], p2[0]) + seg1_max_y = max(p1[1], p2[1]) + seg1_min_y = min(p1[1], p2[1]) + seg1 = seg1_min_x < intersection[0] < seg1_max_x and seg1_min_y < intersection[1] < seg1_max_y + + seg2_max_x = max(p3[0], p4[0]) + seg2_min_x = min(p3[0], p4[0]) + seg2_max_y = max(p3[1], p4[1]) + seg2_min_y = min(p3[1], p4[1]) + seg2 = seg2_min_x <= intersection[0] <= seg2_max_x and seg2_min_y <= intersection[1] <= seg2_max_y + # print("hello") + return seg1 and seg2 if __name__ == '__main__': @@ -183,7 +258,8 @@ if __name__ == '__main__': d.append_point(p) print(f"d.valid_triangulation = {d.valid_triangulation}") - print(d.valid_triangulation) + print(f"finale d.valid_triangulation = {d.valid_triangulation}") + print(f"len(d.valid_triangulation = {len(d.valid_triangulation)}") # hash = lambda n: hash(()) # np.ndarray -- GitLab From 49061242c529ad0339d37f52c89ac9e02d28031c Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 6 Jun 2021 16:16:40 +0200 Subject: [PATCH 40/61] flippen funktioniert noch nicht --- curvepy/testrange.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/curvepy/testrange.py b/curvepy/testrange.py index aaa53ad..8e95040 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -1,7 +1,8 @@ import numpy as np -from typing import List +from typing import List, Tuple import itertools + class Triangle: def __init__(self, A: np.ndarray, B: np.ndarray, C: np.ndarray): self.points = [np.array(x) for x in sorted([[*A], [*B], [*C]])] @@ -33,7 +34,7 @@ class Triangle: points.remove(farthest_p) return farthest_p, points - def __contains__(self, pt: np.ndarray): + def __contains__(self, pt: Tuple[float, float]): """ See: https://www.geeksforgeeks.org/check-whether-a-given-point-lies-inside-a-triangle-or-not/ """ @@ -51,11 +52,11 @@ class Triangle: """ return abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2.0) - def check_if_point_is_in_points(self, pt): - for p in self.points: - if [*p] == [*pt]: - return True - return False + def __str__(self): + return str(self.points) + + def __repr__(self): + return str(self) class Circle: @@ -81,9 +82,10 @@ class Delaunay_triangulation: return [t for t in self._triangles if rem_if(t)] def get_all_points(self): - return np.unique(np.array(sum([t.points for t in self._triangles], [])), axis=0) + return set(sum([t.points for t in self._triangles], [])) - def add_point(self, pt: np.ndarray): + def add_point(self, pt: Tuple[float, float]): + print(f"pt:{pt}") if not self.is_in_range(pt): raise AssertionError("point not in predefined range") t = self.find_triangle_containing(pt) @@ -173,9 +175,12 @@ class Delaunay_triangulation: return self.xmin <= pt[0] <= self.xmax and self.ymin <= pt[1] <= self.ymax def add_new_triangles(self, pt, t): - for a, b in itertools.product(t.points, t.points): - if [*a] != [*b]: - self._triangles.append(Triangle(a, b, pt)) + for a,b in [[0,1],[0,2],[1,2]]: + self._triangles.append(Triangle(t.points[a], t.points[b], pt)) + # print([*itertools.product(t.points, t.points)]) + # for a, b in itertools.product(t.points, t.points): + # if a != b: + # self._triangles.append(Triangle(a, b, pt)) def find_triangle_containing(self, pt) -> Triangle: for t in self._triangles: -- GitLab From 787a84ead5c4ce7db7328839eaa70a043f0f7295 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 6 Jun 2021 16:16:44 +0200 Subject: [PATCH 41/61] flippen funktioniert noch nicht --- curvepy/testrange.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/curvepy/testrange.py b/curvepy/testrange.py index 8e95040..c58f512 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -4,8 +4,8 @@ import itertools class Triangle: - def __init__(self, A: np.ndarray, B: np.ndarray, C: np.ndarray): - self.points = [np.array(x) for x in sorted([[*A], [*B], [*C]])] + def __init__(self, A: Tuple[float, float], B: Tuple[float, float], C: Tuple[float, float]): + self.points = sorted([A, B, C]) self.area = self.calc_area(*A, *B, *C) self.circumcircle = self.calculate_circumcircle() @@ -23,14 +23,14 @@ class Triangle: yu = ((x1 * x1 + y1 * y1) * (x3 - x2) + (x2 * x2 + y2 * y2) * (x1 - x3) + (x3 * x3 + y3 * y3) * (x2 - x1)) / d lines = [[A, B], [B, C], [A, C]] - c, a, b = [np.linalg.norm(x - y) for x, y in lines] + c, a, b = [np.linalg.norm(np.array(x) - np.array(y)) for x, y in lines] R = (a * b * c) / (4 * self.area) return Circle(center=np.array([xu, yu]), radius=R) - def get_farthest_point_away_and_nearest_line_to(self, pt: np.ndarray): + def get_farthest_point_away_and_nearest_line_to(self, pt): points = self.points.copy() - farthest_p = self.points[np.argmax([np.linalg.norm(x - pt) for x in points])] + farthest_p = self.points[np.argmax([np.linalg.norm(np.array(x) - np.array(pt)) for x in points])] points.remove(farthest_p) return farthest_p, points @@ -65,7 +65,7 @@ class Circle: self.radius = radius def point_in_me(self, pt: np.ndarray) -> bool: - return np.linalg.norm(pt - self.center) <= self.radius + return np.linalg.norm(np.array(pt) - self.center) <= self.radius class Delaunay_triangulation: @@ -95,11 +95,12 @@ class Delaunay_triangulation: self.triangle_queue = self._triangles[-3:] pts = self.get_all_points() while self.triangle_queue: + #print(self.triangle_queue) t = self.triangle_queue.pop() for p in pts: - if t.check_if_point_is_in_points(p): + if p in t.points: continue - if t.circumcircle.point_in_me(p): + if t.circumcircle.point_in_me(np.array(p)): self.handle_point_in_circumcircle(t, p) break @@ -132,7 +133,7 @@ class Delaunay_triangulation: B_supert = np.array([self.xmax + (self.xmax - self.xmin) / 2, self.ymin]) C_supert = self.intersect_lines(A_supert, A, B_supert, D) - return Triangle(A_supert, B_supert, C_supert) + return Triangle(tuple(A_supert), tuple(B_supert), tuple(C_supert)) def intersect_lines(self, p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray): """ @@ -191,8 +192,8 @@ class Delaunay_triangulation: if __name__ == '__main__': pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] - d = Delaunay_triangulation(-100, 100, -100, 100) + d = Delaunay_triangulation(-10, 10, -10, 10) for p in pts: - d.add_point(p) + d.add_point(tuple(p)) print(d.triangles) -- GitLab From 70ccde0f3ef8ca72af4311ade1e2b75b54874c53 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 6 Jun 2021 16:16:51 +0200 Subject: [PATCH 42/61] flippen funktioniert noch nicht --- curvepy/testrange2.py | 56 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 curvepy/testrange2.py diff --git a/curvepy/testrange2.py b/curvepy/testrange2.py new file mode 100644 index 0000000..308c684 --- /dev/null +++ b/curvepy/testrange2.py @@ -0,0 +1,56 @@ +import numpy as np +from collections import namedtuple + +valid_range = namedtuple("valid_range", "xmin xmax ymin ymax") + + +class triangle: + def __init__(self, x: np.ndarray, y: np.ndarray, z: np.ndarray): + self.points = sorted([x, y, z]) + self.circumcircle = self.calculate_circumcircle() + + def calculate_circumcircle(self): + """ + :return: + + See: https://de.wikipedia.org/wiki/Umkreis#Koordinaten + """ + [x1, y1], [x2, y2], [x3, y3] = self.points + d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) + + +class dirchlet_tessellation: + def __init__(self, valid_range: valid_range): + self.valid_range = valid_range + self.triangles = [self.create_supertriangle()] + + def get_all_points(self): + return set(np.ravel([t.points for t in self.triangles])) + + def add_point(self, pt: np.ndarray): + if self.is_not_in_range(): + raise AssertionError("point not in predefined range") + t = self.find_triangle_containing(pt) + self.add_new_triangles(pt, t) + for new_t in self.triangles[-3:]: + self.test_circumcircle_property(new_t) + ... + + def test_circumcircle_property(self, t): + ... + + def create_supertriangle(self) -> triangle: + ... + + def is_not_in_range(self): + ... + + def add_new_triangles(self, pt, t): + for a, b in zip(t.points, t.points): + if a != b: + self.triangles.append(triangle(a, b, pt)) + + def find_triangle_containing(self, pt) -> triangle: + # for each triangle + # if triangle.contains(pt): return triangle + return None -- GitLab From ab3fe04eb797a384a1a78ccc4f13eab59e637776 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 6 Jun 2021 17:11:04 +0200 Subject: [PATCH 43/61] =?UTF-8?q?gro=C3=9Fer=20Umkreis=20macht=20alles=20k?= =?UTF-8?q?aputt?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- curvepy/testrange.py | 32 ++++++++++++++++++++++++++------ 1 file changed, 26 insertions(+), 6 deletions(-) diff --git a/curvepy/testrange.py b/curvepy/testrange.py index c58f512..52d6a42 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -7,6 +7,7 @@ class Triangle: def __init__(self, A: Tuple[float, float], B: Tuple[float, float], C: Tuple[float, float]): self.points = sorted([A, B, C]) self.area = self.calc_area(*A, *B, *C) + print(self) self.circumcircle = self.calculate_circumcircle() def calculate_circumcircle(self): @@ -95,7 +96,7 @@ class Delaunay_triangulation: self.triangle_queue = self._triangles[-3:] pts = self.get_all_points() while self.triangle_queue: - #print(self.triangle_queue) + # print(self.triangle_queue) t = self.triangle_queue.pop() for p in pts: if p in t.points: @@ -104,11 +105,30 @@ class Delaunay_triangulation: self.handle_point_in_circumcircle(t, p) break + def flip(self, current_t, p_in_circumcircle): + # Finde Dreieck, was aus p_in_circumcircle und 2 Punkten aus current_t besteht + triangle_with_problematic_edge = None + for t in self._triangles: + if p_in_circumcircle in t and sum([x in t for x in current_t.points]) == 2: + triangle_with_problematic_edge = t + break + + if not triangle_with_problematic_edge: + raise AssertionError("DEBUG: ich bin doof") + + # Edge herausfinden die blockiert + edge_to_remove = [*set(triangle_with_problematic_edge.points).intersection(set(current_t.points))] + + # Punkt finden zum Verbinden mit p_in_circumcircle + point_to_connect_to = set(current_t.points).difference(edge_to_remove).pop() + + return point_to_connect_to, edge_to_remove + def handle_point_in_circumcircle(self, current_t, p_in_circumcircle): - farthest_pt, nearest_pts = current_t.get_farthest_point_away_and_nearest_line_to(p_in_circumcircle) - t1 = Triangle(p_in_circumcircle, farthest_pt, nearest_pts[0]) - t2 = Triangle(p_in_circumcircle, farthest_pt, nearest_pts[1]) - self.update_triangle_structures(t1, t2, nearest_pts) + point_to_connect_to, edge_to_remove = self.flip(current_t, p_in_circumcircle) + t1 = Triangle(p_in_circumcircle, point_to_connect_to, edge_to_remove[0]) + t2 = Triangle(p_in_circumcircle, point_to_connect_to, edge_to_remove[1]) + self.update_triangle_structures(t1, t2, edge_to_remove) def update_triangle_structures(self, t1_new, t2_new, removed_line): triangles_to_remove = [t for t in self._triangles if @@ -176,7 +196,7 @@ class Delaunay_triangulation: return self.xmin <= pt[0] <= self.xmax and self.ymin <= pt[1] <= self.ymax def add_new_triangles(self, pt, t): - for a,b in [[0,1],[0,2],[1,2]]: + for a, b in [[0, 1], [0, 2], [1, 2]]: self._triangles.append(Triangle(t.points[a], t.points[b], pt)) # print([*itertools.product(t.points, t.points)]) # for a, b in itertools.product(t.points, t.points): -- GitLab From db34eab248565da99faa2405767659745b0381f9 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 19 Jun 2021 17:14:14 +0200 Subject: [PATCH 44/61] bevor das flippen angefasst wird - Idee alle pkte im Umkreis sammeln und dann kleinsten umkreis mit neuem punkt bilden --- curvepy/testrange.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/curvepy/testrange.py b/curvepy/testrange.py index 52d6a42..718afaa 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -98,12 +98,14 @@ class Delaunay_triangulation: while self.triangle_queue: # print(self.triangle_queue) t = self.triangle_queue.pop() - for p in pts: - if p in t.points: - continue - if t.circumcircle.point_in_me(np.array(p)): - self.handle_point_in_circumcircle(t, p) - break + deine_pt_liste = [p for p in pts if (p not in t.points) and (t.circumcircle.point_in_me(np.array(p)))] + # for p in pts: + # if p in t.points: + # continue + # if t.circumcircle.point_in_me(np.array(p)): + # self.handle_point_in_circumcircle(t, p) + # break + def flip(self, current_t, p_in_circumcircle): # Finde Dreieck, was aus p_in_circumcircle und 2 Punkten aus current_t besteht -- GitLab From bb12cb1376249b6a86866c8a11ef9228485c3a05 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 19 Jun 2021 20:21:49 +0200 Subject: [PATCH 45/61] Idee mit den kleinsten Kreisen funktioniert --- curvepy/testrange.py | 66 ++++++++++++++++++++++++++++++++++++-------- 1 file changed, 55 insertions(+), 11 deletions(-) diff --git a/curvepy/testrange.py b/curvepy/testrange.py index 718afaa..0d6e603 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -1,13 +1,14 @@ import numpy as np from typing import List, Tuple import itertools +import matplotlib.pyplot as plt class Triangle: def __init__(self, A: Tuple[float, float], B: Tuple[float, float], C: Tuple[float, float]): self.points = sorted([A, B, C]) self.area = self.calc_area(*A, *B, *C) - print(self) + # print(self) self.circumcircle = self.calculate_circumcircle() def calculate_circumcircle(self): @@ -89,23 +90,58 @@ class Delaunay_triangulation: print(f"pt:{pt}") if not self.is_in_range(pt): raise AssertionError("point not in predefined range") - t = self.find_triangle_containing(pt) + t = self.find_triangle_containing(pt) # FUNKTIONIERT NICHT + print(f"in welchem dreieck ist {pt} = {t}") self.add_new_triangles(pt, t) self._triangles.remove(t) # new _triangles self.triangle_queue = self._triangles[-3:] pts = self.get_all_points() while self.triangle_queue: - # print(self.triangle_queue) + print(f"queue = {self.triangle_queue}") t = self.triangle_queue.pop() + print(f"t = {t}") deine_pt_liste = [p for p in pts if (p not in t.points) and (t.circumcircle.point_in_me(np.array(p)))] - # for p in pts: - # if p in t.points: - # continue - # if t.circumcircle.point_in_me(np.array(p)): - # self.handle_point_in_circumcircle(t, p) - # break - + print(f"deine_pt_liste = {deine_pt_liste}") + if not deine_pt_liste: + continue + print(f"es geht was kaputt") + t_new_1 = self.create_smallest_circumcircle_triangle(pt, deine_pt_liste + t.points) + t_new_2 = self.get_second_new_triangle(t_new_1, t, pt) + print(f"t_new_1 = {t_new_1}") + print(f"t_new_2 = {t_new_2}") + print(f"bevor: {self._triangles}") + self.remove_triangles(t, pt, t_new_1) + print(f"danach: {self._triangles}") + self._triangles.append(t_new_1) + self._triangles.append(t_new_2) + self.triangle_queue.append(t_new_1) + self.triangle_queue.append(t_new_2) + + def remove_triangles(self, t, pt, t_new_1): + self._triangles.remove(t) # WIR ENTFERNEN NUR EIN DREIECK ABER ES GIBT NOCH EIN ZWEITES DREIECK MIT DER KANTE + point_not_in_org_t = list(set(t_new_1.points).difference(set(t.points)))[0] + print(f"point_not_in_org_t = {point_not_in_org_t}") + point_not_in_t_new_1 = list(set(t.points).difference((set(t_new_1.points))))[0] + last_point = (0, 0) + # Wir haben 4 punkte und suchen den letzten den wir brauchen für das andere Dreieck um das zu löschen + for a in t.points+t_new_1.points: + print(f"a = {a}") + if a != pt and a != point_not_in_org_t and a != point_not_in_t_new_1: + last_point = a + to_rem = Triangle(point_not_in_org_t, point_not_in_t_new_1, last_point) + if to_rem in self._triangles: + self._triangles.remove(to_rem) + + def get_second_new_triangle(self, t_new_1, t, pt): + point_not_in_org_t = list(set(t_new_1.points).difference(set(t.points)))[0] + point_not_in_t_new_1 = list(set(t.points).difference((set(t_new_1.points))))[0] + return Triangle(pt, point_not_in_t_new_1, point_not_in_org_t) + + def create_smallest_circumcircle_triangle(self, pt, in_circle_pts): + xs = [Triangle(pt, a, b) for a in in_circle_pts for b in in_circle_pts if a != b and a != pt and b != pt] + # xs = sorted(xs, key=lambda t: t.circumcircle.radius) + return min(xs, key=lambda t: t.circumcircle.radius) def flip(self, current_t, p_in_circumcircle): # Finde Dreieck, was aus p_in_circumcircle und 2 Punkten aus current_t besteht @@ -218,4 +254,12 @@ if __name__ == '__main__': for p in pts: d.add_point(tuple(p)) - print(d.triangles) + for tri in d.triangles: + points = np.ravel(tri.points) + plt.triplot(points[0::2], points[1::2]) + + plt.show() + + print(f"anz Dreiecke = {len(d.triangles)}") + + # print(d.triangles) -- GitLab From 9f85e207a8cf657bd2132419f34d79714d4e08a6 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 19 Jun 2021 22:17:46 +0200 Subject: [PATCH 46/61] =?UTF-8?q?etwas=20aufger=C3=A4umt?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- curvepy/example_algo.py | 70 +++++++---------------------------------- 1 file changed, 12 insertions(+), 58 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index ff6365b..e367019 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -40,61 +40,15 @@ class dirichlet_tessellation: # To make a more educated guess we solve any collisions __after__ adding the trivial connections while True: - collisions_to_check = [] - # print(f"collisions_to_check: {collisions_to_check}") - - for neighbour in self.tiles[nearest_p]: - if all(x == y for x, y in zip(p, neighbour)): - continue - - all_collisions = [x for x in self.valid_triangulation if self.intersect(neighbour, p, *x) - and x not in garbage_heap and {neighbour, p} not in garbage_heap] - # print(f"all_collisions = {all_collisions}") - if (not all_collisions) and ({p, neighbour} not in garbage_heap): - # print(f"p wenn keine collisions = {p}") - # print(f"neighbour wenn keine collisions = {neighbour}") - self.valid_triangulation.add((p, neighbour)) - # TODO das muss auch unten gemacht werden wenn die neuen Dreiecke cooler sind - self.tiles[neighbour].add(p) - self.tiles[p].add(neighbour) - elif len(all_collisions) == 1: - # 1 collision could be still fine - collisions_to_check.append([p, neighbour, *all_collisions[0]]) - # More than 1 collision is always not the best possible triangulation TODO Is that true? Probably not ;) + collisions_to_check, garbage_heap = self.calculate_collisions(p, nearest_p, garbage_heap) if not collisions_to_check: return for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: - new_triangles2 = ( - (p, neighbour, collision_edge_p1), - (p, neighbour, collision_edge_p2) - ) - old_triangles2 = ( - (p, collision_edge_p1, collision_edge_p2), - (neighbour, collision_edge_p1, collision_edge_p2) - ) - - print(f"p = {p}") - print(f"neighbour = {neighbour}") - print(f"collision_edge_p1 = {collision_edge_p1}") - print(f"collision_edge_p2 = {collision_edge_p2}") - - angle_list_new_triangle = self._angles_in_triangle(new_triangles2[0]) - # print(f"angle_list_new_triangle = {angle_list_new_triangle}") - angle_rating_new_triangle = sum(abs(60 - angle) for angle in angle_list_new_triangle) - angle_list_new_triangle = self._angles_in_triangle(new_triangles2[1]) - angle_rating_new_triangle += sum(abs(60 - angle) for angle in angle_list_new_triangle) - # print(f"angle_rating_new_triangle = {angle_rating_new_triangle}") - - angle_list_old_triangle = self._angles_in_triangle(old_triangles2[0]) - # print(f"angle_list_old_triangle = {angle_list_old_triangle}") - angle_rating_old_triangle = sum(abs(60 - angle) for angle in angle_list_old_triangle) - angle_list_old_triangle = self._angles_in_triangle(old_triangles2[1]) - angle_rating_old_triangle += sum(abs(60 - angle) for angle in angle_list_old_triangle) - # print(f"angle_rating_new_triangle = {angle_rating_old_triangle}") - - new_is_more_equilateral = angle_rating_old_triangle > angle_rating_new_triangle + new_triangles, old_triangles = self.create_triangles(p, neighbour, collision_edge_p1, collision_edge_p2) + + new_is_more_equilateral = self.calculate_triangle_rating(old_triangles) > self.calculate_triangle_rating(new_triangles) if new_is_more_equilateral: self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) @@ -243,8 +197,8 @@ class dirichlet_tessellation: if __name__ == '__main__': # xs = [h_tuple.create(x) for x in ((1,2), (3,4), (1,4), (3,6)) ] # print(dirichlet_tessellation.intersect(*xs)) - pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] - d = dirichlet_tessellation() + # pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] + # d = dirichlet_tessellation() # print(h_tuple.create([2,1]) == h_tuple.create([1, 2])) # cords A, B, D, E @@ -254,12 +208,12 @@ if __name__ == '__main__': # Winkel neues dreieck1 = 49.4, 70.35, 60.26 # Winkel neues dreieck2 = 59.04, 78.69, 42.27 - for p in pts: - d.append_point(p) - print(f"d.valid_triangulation = {d.valid_triangulation}") - - print(f"finale d.valid_triangulation = {d.valid_triangulation}") - print(f"len(d.valid_triangulation = {len(d.valid_triangulation)}") + # for p in pts: + # d.append_point(p) + # print(f"d.valid_triangulation = {d.valid_triangulation}") + # + # print(f"finale d.valid_triangulation = {d.valid_triangulation}") + # print(f"len(d.valid_triangulation = {len(d.valid_triangulation)}") # hash = lambda n: hash(()) # np.ndarray -- GitLab From fa4aac2305164fe2e26903f52c76dc756c02890d Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 19 Jun 2021 22:18:12 +0200 Subject: [PATCH 47/61] changes bei lquenti --- curvepy/example_algo.py | 53 +++++++++--- curvepy/testrange.py | 12 ++- curvepy/testrange2.py | 176 +++++++++++++++++++++++++++++++--------- 3 files changed, 190 insertions(+), 51 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index e367019..37e2caf 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -1,4 +1,5 @@ import pprint +import random import numpy as np from typing import List, Set, Dict, Tuple, Iterable @@ -50,12 +51,25 @@ class dirichlet_tessellation: new_is_more_equilateral = self.calculate_triangle_rating(old_triangles) > self.calculate_triangle_rating(new_triangles) - if new_is_more_equilateral: - self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) - self.update_neighbour(p, neighbour, collision_edge_p1, collision_edge_p2) - garbage_heap.append({collision_edge_p1, collision_edge_p2}) - else: - garbage_heap.append({p, neighbour}) + self.update(new_is_more_equilateral, garbage_heap, p, neighbour, collision_edge_p1, collision_edge_p2) + + def update(self, new_is_more_equilateral, garbage_heap, p1, p2, p3, p4): + if new_is_more_equilateral: + self.replace_valid_triangulation((p3, p4), (p1, p2)) + self.update_neighbour(p1, p2, p3, p4) + garbage_heap.append({p3, p4}) + else: + garbage_heap.append({p1, p2}) + + def calculate_triangle_rating(self, triangles): + rating = 0 + for t in triangles: + rating += sum(abs(60 - angle) for angle in self._angles_in_triangle(t)) + + return rating + + def create_triangles(self, p1, p2, p3, p4): + return ((p1, p2, p3), (p1, p2, p4)), ((p1, p2, p4), (p2, p3, p4)) # Kantenpaare für rate_tri # for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: @@ -114,12 +128,12 @@ class dirichlet_tessellation: # law of cosine # https://en.wikipedia.org/wiki/Solution_of_triangles - a = np.linalg.norm(triangle[1]-triangle[2]) - b = np.linalg.norm(triangle[0]-triangle[1]) - c = np.linalg.norm(triangle[2]-triangle[0]) + a = np.linalg.norm(triangle[1] - triangle[2]) + b = np.linalg.norm(triangle[0] - triangle[1]) + c = np.linalg.norm(triangle[2] - triangle[0]) - alpha = (180 / np.pi) * np.arccos((b**2+c**2-a**2)/(2*b*c)) - beta = (180 / np.pi) * np.arccos((a**2+c**2-b**2)/(2*a*c)) + alpha = (180 / np.pi) * np.arccos((b ** 2 + c ** 2 - a ** 2) / (2 * b * c)) + beta = (180 / np.pi) * np.arccos((a ** 2 + c ** 2 - b ** 2) / (2 * a * c)) gamma = 180 - alpha - beta return [alpha, beta, gamma] @@ -195,6 +209,23 @@ class dirichlet_tessellation: if __name__ == '__main__': + _min, _max = -10, 10 + n = 7 + points = np.array([np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))]) + np.save('last_run.npy', points) + tri = Delaunay(points) + plt.triplot(points[:, 0], points[:, 1], tri.simplices) + plt.plot(points[:, 0], points[:, 1], 'o') + plt.show() + + d = dirichlet_tessellation() + + for p in points: + d.append_point(p) + + print(f"finale d.valid_triangulation = {d.valid_triangulation}") + print(f"len(d.valid_triangulation = {len(d.valid_triangulation)}") + # xs = [h_tuple.create(x) for x in ((1,2), (3,4), (1,4), (3,6)) ] # print(dirichlet_tessellation.intersect(*xs)) # pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] diff --git a/curvepy/testrange.py b/curvepy/testrange.py index 0d6e603..a6e5c04 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -60,6 +60,11 @@ class Triangle: def __repr__(self): return str(self) + def __eq__(self, other): + if not isinstance(other, Triangle): + return False + return self.points == other.points + class Circle: def __init__(self, center: np.ndarray, radius: float): @@ -69,6 +74,11 @@ class Circle: def point_in_me(self, pt: np.ndarray) -> bool: return np.linalg.norm(np.array(pt) - self.center) <= self.radius + def __eq__(self, other): + if not isinstance(other, Circle): + return False + return self.center == other.center and self.radius == other.radius + class Delaunay_triangulation: def __init__(self, xmin, xmax, ymin, ymax): @@ -81,7 +91,7 @@ class Delaunay_triangulation: def triangles(self): # We have to remove everything containing vertices of the supertriangle rem_if = lambda t: any(pt in t.points for pt in self.supertriangle.points) - return [t for t in self._triangles if rem_if(t)] + return [t for t in self._triangles if not rem_if(t)] def get_all_points(self): return set(sum([t.points for t in self._triangles], [])) diff --git a/curvepy/testrange2.py b/curvepy/testrange2.py index 308c684..66cb497 100644 --- a/curvepy/testrange2.py +++ b/curvepy/testrange2.py @@ -1,56 +1,154 @@ import numpy as np -from collections import namedtuple +from math import sqrt +from typing import Tuple, List, Dict, Optional -valid_range = namedtuple("valid_range", "xmin xmax ymin ymax") +PointIndex = int +Triangle = Tuple[PointIndex, PointIndex, PointIndex] +Neighbours = List[Optional[Triangle]] +Radius = float +Center = Tuple[float, float] +Circle = Tuple[Center, Radius] -class triangle: - def __init__(self, x: np.ndarray, y: np.ndarray, z: np.ndarray): - self.points = sorted([x, y, z]) - self.circumcircle = self.calculate_circumcircle() - def calculate_circumcircle(self): +class Delaunay: + def __init__(self, center=(0, 0), radius=1000): + center = np.array(center) + + # Rechteck erstellen indem man den Radius mal die Basisvektoren nimmt + self.coords = [center + radius * np.array([i, j]) for i, j in [(-1, -1), (1, -1), (1, 1), (-1, 1)]] + + # Triangles werden wiefolgt gespeichert: + # Der Key ist ein 3-Tuple, welches CCW aus den Indexen der self.coords bestehen. + # Sprich das Tuple + # (2,3,1) + # bedeutet das Dreieck + # (self.coords[2], self.coords[3], self.coords[1]) + self.triangles: Dict[Triangle, Neighbours] = {} + + # Hier speichern wir zu jedem Triangle (Ueber indexe definiert) das 2-Tuple (center, radius) + self.circles: Dict[Triangle, Circle] = {} + + # Hier werden die Indexe der in den coords als Triangles gespeichert + index_lower_left = 0 + index_lower_right = 1 + index_upper_right = 2 + index_upper_left = 3 + + # Wir splitten den durch die self.coords definierten sapce in 2 Triangle auf, indem wir so eine Linie + # durchziehen: + # x────────────────────────────────┐ + # xx │ + # │ xx │ + # │ xxx T2 │ + # │ xxxx │ + # │ xxx │ + # │ xxx │ + # │ xxxx │ + # │ xxxx │ + # │ xxx │ + # │ T1 xxxx │ + # │ xxx │ + # │ xxx │ + # │ xxx │ + # │ xxx│ + # └────────────────────────────────┴ + # Wir definieren sie Counter Clock Wise (CCW), angefangen beim Punkt gegenueber der Hypothenuse + first_super_triangle = (index_lower_left, index_lower_right, index_upper_left) + second_super_triangle = (index_upper_right, index_upper_left, index_lower_right) + + # TODO: Stimmt das, dass wir so die Datenstruktur nutzen? + # Wir speichern die Triangles mit ihren Nachbarn. Offensichtlich sind T1 und T2 oben benachbart und haben + # sonst keine Nachbarn, da sie den validen Bereich begrenzen + self.triangles[first_super_triangle] = [second_super_triangle, None, None] + self.triangles[second_super_triangle] = [first_super_triangle, None, None] + + self.circles[first_super_triangle] = self.circumcenter(first_super_triangle) + self.circles[second_super_triangle] = self.circumcenter(second_super_triangle) + + def circumcenter(self, tri): + """Compute circumcenter and circumradius of a triangle in 2D. + Uses an extension of the method described here: + http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html """ - :return: + pts = np.asarray([self.coords[v] for v in tri]) + pts2 = np.dot(pts, pts.T) + A = np.bmat([[2 * pts2, [[1], + [1], + [1]]], + [[[1, 1, 1, 0]]]]) + + b = np.hstack((np.sum(pts * pts, axis=1), [1])) + x = np.linalg.solve(A, b) + bary_coords = x[:-1] + center = np.dot(bary_coords, pts) - See: https://de.wikipedia.org/wiki/Umkreis#Koordinaten + # radius = np.linalg.norm(pts[0] - center) # euclidean distance + radius = np.sum(np.square(pts[0] - center)) # squared distance + return (center, radius) + + def inCircleFast(self, tri, p): + """Check if point p is inside of precomputed circumcircle of tri. """ - [x1, y1], [x2, y2], [x3, y3] = self.points - d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) + center, radius = self.circles[tri] + return np.sum(np.square(center - p)) <= radius + def add_point(self, p): + """ + TODO: Replace mich + Add point to the current DT, and refine it using Bowyer-Watson + """ + p = np.array(p) + # Wir schieben den Punkt ganz nach hinten (da wir ja alles in indexen speichern und somit darauf angewiesen sind + # dass die vorherigen gleichbleiben). + # (Laenge ist ja 1-indexiert, somit entspricht es dem naechsten Punkt) + idx = len(self.coords) + self.coords.append(p) + # also gilt + # assert self.coords[idx] == self.coords[-1] -class dirchlet_tessellation: - def __init__(self, valid_range: valid_range): - self.valid_range = valid_range - self.triangles = [self.create_supertriangle()] + # Wir nehmen an, dass vor diesem Punkt die Delaunay-Triangulation optimal war + # Somit muessen wir nur diesen Punkt beachten. + # Per Definition ist ein Triangulationsdreieck nur valide wenn in seinem Umkreis kein weiterer Punkt liegt. + bad_trinagles = [tri for tri in self.triangles if self.inCircleFast(tri, p)] - def get_all_points(self): - return set(np.ravel([t.points for t in self.triangles])) + # TODO: Verstehe folgendes + # Find the CCW boundary (star shape) of the bad triangles, + # expressed as a list of edges (point pairs) and the opposite + # triangle to each edge. + # + # Bisherige Vermutung: + # - Star Shape == Triforce mit dem Punkt im Mitteldreieck + # - each edge = Jede Kante des Dreiecks, welches den Punkt im Umkreis beinhaltet + # - opposite triangle = Das Dreieck, was nicht schlecht ist und trotzdem die edge teilt. + # + # UPDATE: Wir sammeln hierbei alle nicht-problematischen Nachbarn (siehe if in der while) + boundaries: List[Tuple[PointIndex, PointIndex, Triangle]] = [] - def add_point(self, pt: np.ndarray): - if self.is_not_in_range(): - raise AssertionError("point not in predefined range") - t = self.find_triangle_containing(pt) - self.add_new_triangles(pt, t) - for new_t in self.triangles[-3:]: - self.test_circumcircle_property(new_t) - ... + # Wir starten bei einem "random" Triangle mit dessen ersten Edge + T = bad_trinagles[0] + edge = 0 - def test_circumcircle_property(self, t): - ... + # Jetzt suchen wir das opposite triangle fuer jede edge, angefangen mit edge = 0 weil irgendwo muss man ja + # + while True: + tri_opposite_to: Triangle = self.triangles[T][edge] + if tri_opposite_to not in bad_trinagles: + # wir speichern die Kante ab an die Tri_op angrenzt + # und zwar die indices der Punkte die die Kante aufspannen + boundaries.append((T[(edge + 1) % 3], T[(edge - 1) % 3], tri_opposite_to)) - def create_supertriangle(self) -> triangle: - ... + edge = (edge + 1) % 3 - def is_not_in_range(self): - ... + # Wenn der naechste Punkt + if boundaries[0][0] == boundaries[-1][1]: + break + else: + # CCW edge nach der edge die wir grade betrachten + # da an dieser noch aktuellen edge das tri_op liegt + # wir wollen aber die ccw nächste da wir die aktuelle abgearbeitet haben + same_edge_other_triangle = self.triangles[tri_opposite_to].index(T) + edge = (same_edge_other_triangle + 1) % 3 + T = tri_opposite_to - def add_new_triangles(self, pt, t): - for a, b in zip(t.points, t.points): - if a != b: - self.triangles.append(triangle(a, b, pt)) - def find_triangle_containing(self, pt) -> triangle: - # for each triangle - # if triangle.contains(pt): return triangle - return None -- GitLab From e7123bfdb89cbcd4f6d728ab1cb14dc742854a6c Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 19 Jun 2021 22:18:21 +0200 Subject: [PATCH 48/61] changes bei timo --- curvepy/example_algo.py | 79 +++++++++++------------------------------ curvepy/testrange.py | 15 ++++++++ 2 files changed, 36 insertions(+), 58 deletions(-) diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py index 37e2caf..9f60ebb 100644 --- a/curvepy/example_algo.py +++ b/curvepy/example_algo.py @@ -2,6 +2,8 @@ import pprint import random import numpy as np from typing import List, Set, Dict, Tuple, Iterable +from scipy.spatial import Delaunay +import matplotlib.pyplot as plt class h_tuple(np.ndarray): @@ -71,57 +73,25 @@ class dirichlet_tessellation: def create_triangles(self, p1, p2, p3, p4): return ((p1, p2, p3), (p1, p2, p4)), ((p1, p2, p4), (p2, p3, p4)) - # Kantenpaare für rate_tri - # for p, neighbour, collision_edge_p1, collision_edge_p2 in collisions_to_check: - # new_triangles = ( - # # first triangle - # ((p, neighbour), (p, collision_edge_p1)), - # # ((p, neighbour), (neighbour, collision_edge_p1)), <- HEUTE REINKOPIERT - # ((p, collision_edge_p1), (neighbour, collision_edge_p1)), # maybe jetzt die richtige kante? - # ((p, neighbour), (neighbour, collision_edge_p1)), - # # second triangle - # ((p, neighbour), (p, collision_edge_p2)), - # ((p, collision_edge_p2), (neighbour, collision_edge_p2)), - # ((p, neighbour), (neighbour, collision_edge_p2)) - # ) - # - # old_triangles = ( - # # first triangle - # ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, p)), - # ((collision_edge_p1, p), (collision_edge_p2, p)), - # ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, p)), - # # second triangle - # ((collision_edge_p1, collision_edge_p2), (collision_edge_p1, neighbour)), - # ((collision_edge_p1, neighbour), (collision_edge_p2, neighbour)), - # ((collision_edge_p1, collision_edge_p2), (collision_edge_p2, neighbour)) - # ) - # - # rate_tri = lambda t: sum(abs(60 - self._angle(*a)) for a in t) - # a = rate_tri(old_triangles) - # b = rate_tri(new_triangles) - # new_is_more_equilateral = a > b - # print(f"rate_tri(old_triangles) = {a}") - # print(f"rate_tri(new_triangles) = {b}") - # if new_is_more_equilateral: - # self.replace_valid_triangulation((collision_edge_p1, collision_edge_p2), (p, neighbour)) - # self.update_neighbour(p, neighbour, collision_edge_p1, collision_edge_p2) - # garbage_heap.append({collision_edge_p1, collision_edge_p2}) - # else: - # garbage_heap.append({p, neighbour}) - # print(f"garbage_heap = {garbage_heap}") + def calculate_collisions(self, p, nearest_p: h_tuple, garbage_heap): + collisions_to_check = [] + # print(f"collisions_to_check: {collisions_to_check}") - @staticmethod - def _angle(edge1, edge2): - print(f"edge1 = {edge1}, edge2 = {edge2}") - # print(f"(edge1[1][0] - edge1[0][0]) = {(edge1[1][0] - edge1[0][0])}") - # print(f"(edge2[1][0] - edge2[0][0]) = {(edge2[1][0] - edge2[0][0])}") - # print(f"_angle = {abs(180 / np.pi * (np.arctan((edge1[1][1] - edge1[0][1]) / (edge1[1][0] - edge1[0][0])) - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0]))))}") - slope1 = dirichlet_tessellation.slope(edge1) - slope2 = dirichlet_tessellation.slope(edge2) - print(f"(180 / np.pi) * np.arctan(abs((slope2 - slope1)/(1 + slope1*slope2))) = {(180 / np.pi) * np.arctan(abs((slope2 - slope1)/(1 + slope1*slope2)))}") - return (180 / np.pi) * np.arctan(abs((slope2 - slope1)/(1 + slope1*slope2))) - # return abs(180 / np.pi * (np.arctan((edge1[1][1] - edge1[0][1]) / (edge1[1][0] - edge1[0][0])) - # - np.arctan((edge2[1][1] - edge2[0][1]) / (edge2[1][0] - edge2[0][0])))) + for neighbour in self.tiles[nearest_p]: + if all(x == y for x, y in zip(p, neighbour)): + continue + + all_collisions = [x for x in self.valid_triangulation if self.intersect(neighbour, p, *x) + and x not in garbage_heap and {neighbour, p} not in garbage_heap] + + if (not all_collisions) and ({p, neighbour} not in garbage_heap): + self.valid_triangulation.add((p, neighbour)) + self.tiles[neighbour].add(p) + self.tiles[p].add(neighbour) + elif len(all_collisions) == 1: + collisions_to_check.append([p, neighbour, *all_collisions[0]]) + + return collisions_to_check, garbage_heap @staticmethod def _angles_in_triangle(triangle): @@ -186,13 +156,7 @@ class dirichlet_tessellation: return dirichlet_tessellation.intersection_is_between(intersection, p1, p2, p3, p4) @staticmethod - def intersection_is_between(intersection: h_tuple, p1: h_tuple, p2: h_tuple, p3, p4): - # a = np.linalg.norm(intersection - p1) + np.linalg.norm(intersection - p2) - # b = np.linalg.norm(p2 - p1) - # c = a == b - # print(c) - # return c - # return abs(intersection - p1) + abs(intersection - p2) == abs(p2 - p1) + def intersection_is_between(intersection: h_tuple, p1: h_tuple, p2: h_tuple, p3: h_tuple, p4: h_tuple): seg1_max_x = max(p1[0], p2[0]) seg1_min_x = min(p1[0], p2[0]) seg1_max_y = max(p1[1], p2[1]) @@ -204,7 +168,6 @@ class dirichlet_tessellation: seg2_max_y = max(p3[1], p4[1]) seg2_min_y = min(p3[1], p4[1]) seg2 = seg2_min_x <= intersection[0] <= seg2_max_x and seg2_min_y <= intersection[1] <= seg2_max_y - # print("hello") return seg1 and seg2 diff --git a/curvepy/testrange.py b/curvepy/testrange.py index a6e5c04..c47549a 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -240,6 +240,21 @@ class Delaunay_triangulation: # we divide with z to turn back to 2D space return np.array([x / z, y / z]) + @staticmethod + def intersection_is_between(intersection, p1, p2, p3, p4): + seg1_max_x = max(p1[0], p2[0]) + seg1_min_x = min(p1[0], p2[0]) + seg1_max_y = max(p1[1], p2[1]) + seg1_min_y = min(p1[1], p2[1]) + seg1 = seg1_min_x < intersection[0] < seg1_max_x and seg1_min_y < intersection[1] < seg1_max_y + + seg2_max_x = max(p3[0], p4[0]) + seg2_min_x = min(p3[0], p4[0]) + seg2_max_y = max(p3[1], p4[1]) + seg2_min_y = min(p3[1], p4[1]) + seg2 = seg2_min_x <= intersection[0] <= seg2_max_x and seg2_min_y <= intersection[1] <= seg2_max_y + return seg1 and seg2 + def is_in_range(self, pt): return self.xmin <= pt[0] <= self.xmax and self.ymin <= pt[1] <= self.ymax -- GitLab From 8af14cc80627b58d10aed859360509f70fe2541f Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 19 Jun 2021 22:18:32 +0200 Subject: [PATCH 49/61] unversioned files --- curvepy/last_run.npy | Bin 0 -> 184 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 curvepy/last_run.npy diff --git a/curvepy/last_run.npy b/curvepy/last_run.npy new file mode 100644 index 0000000000000000000000000000000000000000..dd8a1227ea976ce63992e8ee8e5d61d18515afeb GIT binary patch literal 184 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlWC%^qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I%|ItoUbItsN4WCJcH1_lOZAZ7z%Rw$btNOM5hATdTLA0+no|NsC0fEWOU CL?fdB literal 0 HcmV?d00001 -- GitLab From 153f3ee16e2a99bf939d4a44aaaaba7ce3a385cc Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 19 Jun 2021 23:43:41 +0200 Subject: [PATCH 50/61] halbwegs --- curvepy/testrange.py | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/curvepy/testrange.py b/curvepy/testrange.py index c47549a..df4ba10 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -2,6 +2,8 @@ import numpy as np from typing import List, Tuple import itertools import matplotlib.pyplot as plt +import random as rd +from scipy.spatial import Delaunay class Triangle: @@ -45,7 +47,8 @@ class Triangle: area1 = self.calc_area(x, y, x2, y2, x3, y3) area2 = self.calc_area(x1, y1, x, y, x3, y3) area3 = self.calc_area(x1, y1, x2, y2, x, y) - return self.area == area1 + area2 + area3 + delta = 0.000000000000550 # das ist neu + return area1 + area2 + area3 - delta <= self.area <= area1 + area2 + area3 + delta @staticmethod def calc_area(x1, y1, x2, y2, x3, y3): @@ -100,7 +103,7 @@ class Delaunay_triangulation: print(f"pt:{pt}") if not self.is_in_range(pt): raise AssertionError("point not in predefined range") - t = self.find_triangle_containing(pt) # FUNKTIONIERT NICHT + t = self.find_triangle_containing(pt) print(f"in welchem dreieck ist {pt} = {t}") self.add_new_triangles(pt, t) self._triangles.remove(t) @@ -135,11 +138,12 @@ class Delaunay_triangulation: point_not_in_t_new_1 = list(set(t.points).difference((set(t_new_1.points))))[0] last_point = (0, 0) # Wir haben 4 punkte und suchen den letzten den wir brauchen für das andere Dreieck um das zu löschen - for a in t.points+t_new_1.points: - print(f"a = {a}") + for a in t.points + t_new_1.points: if a != pt and a != point_not_in_org_t and a != point_not_in_t_new_1: last_point = a + continue to_rem = Triangle(point_not_in_org_t, point_not_in_t_new_1, last_point) + print(f"to_rem = {to_rem}") if to_rem in self._triangles: self._triangles.remove(to_rem) @@ -275,7 +279,17 @@ class Delaunay_triangulation: if __name__ == '__main__': pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] + dein_min, dein_max = -5, 5 + n = 7 + pts = np.array([[rd.uniform(dein_min, dein_max), rd.uniform(dein_min, dein_max)] for _ in range(n)]) + + tri = Delaunay(pts) + plt.triplot(pts[:, 0], pts[:, 1], tri.simplices) + plt.plot(pts[:, 0], pts[:, 1], 'o') + plt.show() + d = Delaunay_triangulation(-10, 10, -10, 10) + for p in pts: d.add_point(tuple(p)) -- GitLab From 6716a2b1b6b6e5491109b873938531d1577801c4 Mon Sep 17 00:00:00 2001 From: Lars Quentin Date: Sun, 20 Jun 2021 00:57:25 +0200 Subject: [PATCH 51/61] works, pre refactor (already pretty neat) --- curvepy/my_playground/another_one.py | 175 +++++++++++++++++++++++++++ 1 file changed, 175 insertions(+) create mode 100644 curvepy/my_playground/another_one.py diff --git a/curvepy/my_playground/another_one.py b/curvepy/my_playground/another_one.py new file mode 100644 index 0000000..778f723 --- /dev/null +++ b/curvepy/my_playground/another_one.py @@ -0,0 +1,175 @@ +import numpy as np +from math import sqrt +import random as rd +from functools import cached_property +from typing import List, Tuple, Any +from scipy.spatial import Delaunay +import matplotlib.pyplot as plt + +""" +TODO: Dimension Checks +TODO: Tests +TODO: Voronoi +""" + +Point2D = Tuple[float, float] +Edge2D = Tuple[Point2D, Point2D] + + +class Circle: + def __init__(self, center: Point2D, radius: float): + self.center = np.array(center) + self.radius = radius + + def __contains__(self, pt: Point2D) -> bool: + return np.linalg.norm(np.array(pt) - self.center) <= self.radius + + def __str__(self) -> str: + return f"(CENTER: {self.center}, RADIUS: {self.radius})" + + def __repr__(self) -> str: + return f"" + + def __eq__(self, other: Any) -> bool: + return isinstance(other, Circle) and self.center == other.center and self.radius == other.radius + + +class Triangle: + def __init__(self, a: Point2D, b: Point2D, c: Point2D): + self._points: List[Point2D] = sorted([a, b, c]) + + @property + def points(self) -> List[Point2D]: + return self._points + + @cached_property + def area(self) -> float: + a, b, c = self.points + return self.calc_area(*a, *b, *c) + + @cached_property + def circumcircle(self) -> Circle: + """ + :return: + + See: https://de.wikipedia.org/wiki/Umkreis#Koordinaten + See: https://de.wikipedia.org/wiki/Umkreis#Radius + """ + A, B, C = self.points + [x1, y1], [x2, y2], [x3, y3] = A, B, C + d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) + xu = ((x1 * x1 + y1 * y1) * (y2 - y3) + (x2 * x2 + y2 * y2) * (y3 - y1) + (x3 * x3 + y3 * y3) * (y1 - y2)) / d + yu = ((x1 * x1 + y1 * y1) * (x3 - x2) + (x2 * x2 + y2 * y2) * (x1 - x3) + (x3 * x3 + y3 * y3) * (x2 - x1)) / d + + lines = [[A, B], [B, C], [A, C]] + c, a, b = [np.linalg.norm(np.array(x) - np.array(y)) for x, y in lines] + + R = (a * b * c) / (4 * self.area) + return Circle(center=(xu, yu), radius=R) + + @cached_property + def edges(self) -> List[Edge2D]: + # (<) not only removes same elements but also duplicates in different order + return [(x, y) for x in self.points for y in self.points if x < y] + + @staticmethod + def calc_area(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float) -> float: + """ + See: https://www.mathopenref.com/coordtrianglearea.html + """ + return abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2.0) + + def __str__(self) -> str: + return str(self.points) + + def __repr__(self) -> str: + return f"" + + def __eq__(self, other: Any) -> bool: + return isinstance(other, Triangle) and self.points == other.points + + +class DelaunayTriangulation2D: + def __init__(self, center: Point2D = (0, 0), radius: float = 1000): + self.supertriangles: List[Triangle] = self._create_supertriangles(center, radius) + self._triangles: List[Triangle] = [*self.supertriangles] + + @property + def triangles(self): + # We have to remove everything containing vertices of the supertriangle + all_points_in_supertriangle = sum([x.points for x in self.supertriangles], []) + remove_if = lambda t: any(pt in t.points for pt in all_points_in_supertriangle) + return [t for t in self._triangles if not remove_if(t)] + + def _create_supertriangles(self, center: Point2D, radius): + # Since we have to start with a valid triangulation, we split our allowed range into 2 triangles like that: + # x────────────────────────────────┐ + # xx │ + # │ xx │ + # │ xxx T2 │ + # │ xxxx │ + # │ xxx │ + # │ xxx │ + # │ xxxx │ + # │ xxxx │ + # │ xxx │ + # │ T1 xxxx │ + # │ xxx │ + # │ xxx │ + # │ xxx │ + # │ xxx│ + # └────────────────────────────────┴ + # Those 2 are called the "supertriangles" in most literature + + # for easier calculation + center = np.array(center) + base_rectangle = [center + radius * np.array([i, j]) for i, j in [(-1, -1), (1, -1), (1, 1), (-1, 1)]] + lower_left, lower_right, upper_right, upper_left = [tuple(x) for x in base_rectangle] + lower_triangle = Triangle(lower_left, lower_right, upper_left) + upper_triangle = Triangle(upper_right, upper_left, lower_right) + return [lower_triangle, upper_triangle] + + def add_point(self, p: Point2D): + bad_triangles = [tri for tri in self._triangles if p in tri.circumcircle] + + polygon = [] + for bad_triangle in bad_triangles: + # An edge is part of the boundary iff it doesn't is not part of another bad triangle + other_bad_triangles = list(bad_triangles) + other_bad_triangles.remove(bad_triangle) + for edge in bad_triangle.edges: + if all(edge not in other.edges for other in other_bad_triangles): + polygon.append(edge) + + # set() works since all edges are ordered, we don't care about the order and tuples are hashable + polygon = list(set(polygon)) + + # remove all bad ones + for tri in bad_triangles: + self._triangles.remove(tri) + + for edge in polygon: + self._triangles.append(Triangle(p, *edge)) + + +if __name__ == '__main__': + n = 10 + min, max = -10, 10 + pts = [(rd.uniform(min, max), rd.uniform(min, max)) for _ in range(n)] + d = DelaunayTriangulation2D(radius=n + 5) # Buffer for rounding errors + for p in pts: + d.add_point(p) + + figure, axis = plt.subplots(2) + + axis[0].set_title("meins") + for tri in d.triangles: + points = np.ravel(tri.points) + axis[0].triplot(points[0::2], points[1::2]) + + axis[1].set_title("scipy") + points = np.array([list(x) for x in pts]) + scipy_tri = Delaunay(pts) + axis[1].triplot(points[:,0], points[:,1], scipy_tri.simplices) + axis[1].plot(points[:,0], points[:,1], 'o') + plt.show() -- GitLab From 4f83854b867aae1cbd41eebc49611d8734be231c Mon Sep 17 00:00:00 2001 From: Lars Quentin Date: Sun, 20 Jun 2021 01:10:40 +0200 Subject: [PATCH 52/61] Dafuer haette ich kein Energy trinken muessen --- .../another_one.py => besser_wirds_nicht.py} | 46 +++++++++---------- 1 file changed, 23 insertions(+), 23 deletions(-) rename curvepy/{my_playground/another_one.py => besser_wirds_nicht.py} (87%) diff --git a/curvepy/my_playground/another_one.py b/curvepy/besser_wirds_nicht.py similarity index 87% rename from curvepy/my_playground/another_one.py rename to curvepy/besser_wirds_nicht.py index 778f723..ce7110c 100644 --- a/curvepy/my_playground/another_one.py +++ b/curvepy/besser_wirds_nicht.py @@ -1,17 +1,10 @@ import numpy as np -from math import sqrt import random as rd from functools import cached_property from typing import List, Tuple, Any from scipy.spatial import Delaunay import matplotlib.pyplot as plt -""" -TODO: Dimension Checks -TODO: Tests -TODO: Voronoi -""" - Point2D = Tuple[float, float] Edge2D = Tuple[Point2D, Point2D] @@ -121,7 +114,7 @@ class DelaunayTriangulation2D: # └────────────────────────────────┴ # Those 2 are called the "supertriangles" in most literature - # for easier calculation + # np.array for easier calculation center = np.array(center) base_rectangle = [center + radius * np.array([i, j]) for i, j in [(-1, -1), (1, -1), (1, 1), (-1, 1)]] lower_left, lower_right, upper_right, upper_left = [tuple(x) for x in base_rectangle] @@ -132,25 +125,32 @@ class DelaunayTriangulation2D: def add_point(self, p: Point2D): bad_triangles = [tri for tri in self._triangles if p in tri.circumcircle] - polygon = [] - for bad_triangle in bad_triangles: - # An edge is part of the boundary iff it doesn't is not part of another bad triangle - other_bad_triangles = list(bad_triangles) - other_bad_triangles.remove(bad_triangle) - for edge in bad_triangle.edges: - if all(edge not in other.edges for other in other_bad_triangles): - polygon.append(edge) - - # set() works since all edges are ordered, we don't care about the order and tuples are hashable - polygon = list(set(polygon)) + # An edge is part of the boundary iff it doesn't is not part of another bad triangle + boundaries = self._find_edges_only_used_by_a_single_triangle(bad_triangles) # remove all bad ones for tri in bad_triangles: self._triangles.remove(tri) - - for edge in polygon: + # Replace the hole with the boundaries and our new point + for edge in boundaries: self._triangles.append(Triangle(p, *edge)) + @staticmethod + def _find_edges_only_used_by_a_single_triangle(triangles, no_duplicates=True): + ret = [] + for triangle in triangles: + others = list(triangles) + others.remove(triangle) + for edge in triangle.edges: + if all(edge not in other.edges for other in others): + ret.append(edge) + + # set() works since all edges are ordered, we don't care about the order and tuples are hashable + if no_duplicates: + ret = list(set(ret)) + + return ret + if __name__ == '__main__': n = 10 @@ -170,6 +170,6 @@ if __name__ == '__main__': axis[1].set_title("scipy") points = np.array([list(x) for x in pts]) scipy_tri = Delaunay(pts) - axis[1].triplot(points[:,0], points[:,1], scipy_tri.simplices) - axis[1].plot(points[:,0], points[:,1], 'o') + axis[1].triplot(points[:, 0], points[:, 1], scipy_tri.simplices) + axis[1].plot(points[:, 0], points[:, 1], 'o') plt.show() -- GitLab From 8feb25a0c72e2843595b5fb8a18e4b56c8eb2dfa Mon Sep 17 00:00:00 2001 From: Lars Quentin Date: Sun, 20 Jun 2021 11:56:14 +0200 Subject: [PATCH 53/61] Klappt jetzt gegen Scipy und Reference Implementation --- curvepy/besser_wirds_nicht.py | 18 ++- curvepy/reference_implementation.py | 233 ++++++++++++++++++++++++++++ 2 files changed, 249 insertions(+), 2 deletions(-) create mode 100644 curvepy/reference_implementation.py diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py index ce7110c..e8588a1 100644 --- a/curvepy/besser_wirds_nicht.py +++ b/curvepy/besser_wirds_nicht.py @@ -5,6 +5,8 @@ from typing import List, Tuple, Any from scipy.spatial import Delaunay import matplotlib.pyplot as plt +from curvepy.reference_implementation import Delaunay2D + Point2D = Tuple[float, float] Edge2D = Tuple[Point2D, Point2D] @@ -160,7 +162,8 @@ if __name__ == '__main__': for p in pts: d.add_point(p) - figure, axis = plt.subplots(2) + plt.rcParams["figure.figsize"] = (10,30) + figure, axis = plt.subplots(3) axis[0].set_title("meins") for tri in d.triangles: @@ -172,4 +175,15 @@ if __name__ == '__main__': scipy_tri = Delaunay(pts) axis[1].triplot(points[:, 0], points[:, 1], scipy_tri.simplices) axis[1].plot(points[:, 0], points[:, 1], 'o') - plt.show() + + axis[2].set_title("reference implementation") + d2 = Delaunay2D(radius=n+5) + for p in pts: + d2.addPoint(p) + coord, tris = d2.exportDT() + my_tris = [(coord[a], coord[b], coord[c]) for a,b,c in tris] + for tri in my_tris: + points = np.ravel(tri) + axis[2].triplot(points[0::2], points[1::2]) + + plt.show() \ No newline at end of file diff --git a/curvepy/reference_implementation.py b/curvepy/reference_implementation.py new file mode 100644 index 0000000..7fc3939 --- /dev/null +++ b/curvepy/reference_implementation.py @@ -0,0 +1,233 @@ +# -*- coding: ascii -*- +""" +Simple structured Delaunay triangulation in 2D with Bowyer-Watson algorithm. + +Written by Jose M. Espadero ( http://github.com/jmespadero/pyDelaunay2D ) +Based on code from Ayron Catteau. Published at http://github.com/ayron/delaunay + +Just pretend to be simple and didactic. The only requisite is numpy. +Robust checks disabled by default. May not work in degenerate set of points. +""" + +import numpy as np +from math import sqrt + + +class Delaunay2D: + """ + Class to compute a Delaunay triangulation in 2D + ref: http://en.wikipedia.org/wiki/Bowyer-Watson_algorithm + ref: http://www.geom.uiuc.edu/~samuelp/del_project.html + """ + + def __init__(self, center=(0, 0), radius=9999): + """ Init and create a new frame to contain the triangulation + center -- Optional position for the center of the frame. Default (0,0) + radius -- Optional distance from corners to the center. + """ + center = np.asarray(center) + # Create coordinates for the corners of the frame + self.coords = [center+radius*np.array((-1, -1)), + center+radius*np.array((+1, -1)), + center+radius*np.array((+1, +1)), + center+radius*np.array((-1, +1))] + + # Create two dicts to store triangle neighbours and circumcircles. + self.triangles = {} + self.circles = {} + + # Create two CCW triangles for the frame + T1 = (0, 1, 3) + T2 = (2, 3, 1) + self.triangles[T1] = [T2, None, None] + self.triangles[T2] = [T1, None, None] + + # Compute circumcenters and circumradius for each triangle + for t in self.triangles: + self.circles[t] = self.circumcenter(t) + + def circumcenter(self, tri): + """Compute circumcenter and circumradius of a triangle in 2D. + Uses an extension of the method described here: + http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html + """ + pts = np.asarray([self.coords[v] for v in tri]) + pts2 = np.dot(pts, pts.T) + A = np.bmat([[2 * pts2, [[1], + [1], + [1]]], + [[[1, 1, 1, 0]]]]) + + b = np.hstack((np.sum(pts * pts, axis=1), [1])) + x = np.linalg.solve(A, b) + bary_coords = x[:-1] + center = np.dot(bary_coords, pts) + + # radius = np.linalg.norm(pts[0] - center) # euclidean distance + radius = np.sum(np.square(pts[0] - center)) # squared distance + return (center, radius) + + def inCircleFast(self, tri, p): + """Check if point p is inside of precomputed circumcircle of tri. + """ + center, radius = self.circles[tri] + return np.sum(np.square(center - p)) <= radius + + def inCircleRobust(self, tri, p): + """Check if point p is inside of circumcircle around the triangle tri. + This is a robust predicate, slower than compare distance to centers + ref: http://www.cs.cmu.edu/~quake/robust.html + """ + m1 = np.asarray([self.coords[v] - p for v in tri]) + m2 = np.sum(np.square(m1), axis=1).reshape((3, 1)) + m = np.hstack((m1, m2)) # The 3x3 matrix to check + return np.linalg.det(m) <= 0 + + def addPoint(self, p): + """Add a point to the current DT, and refine it using Bowyer-Watson. + """ + p = np.asarray(p) + idx = len(self.coords) + # print("coords[", idx,"] ->",p) + self.coords.append(p) + + # Search the triangle(s) whose circumcircle contains p + bad_triangles = [] + for T in self.triangles: + # Choose one method: inCircleRobust(T, p) or inCircleFast(T, p) + if self.inCircleFast(T, p): + bad_triangles.append(T) + + # Find the CCW boundary (star shape) of the bad triangles, + # expressed as a list of edges (point pairs) and the opposite + # triangle to each edge. + boundary = [] + # Choose a "random" triangle and edge + T = bad_triangles[0] + edge = 0 + # get the opposite triangle of this edge + while True: + # Check if edge of triangle T is on the boundary... + # if opposite triangle of this edge is external to the list + tri_op = self.triangles[T][edge] + if tri_op not in bad_triangles: + # Insert edge and external triangle into boundary list + boundary.append((T[(edge+1) % 3], T[(edge-1) % 3], tri_op)) + + # Move to next CCW edge in this triangle + edge = (edge + 1) % 3 + + # Check if boundary is a closed loop + if boundary[0][0] == boundary[-1][1]: + break + else: + # Move to next CCW edge in opposite triangle + edge = (self.triangles[tri_op].index(T) + 1) % 3 + T = tri_op + + # Remove triangles too near of point p of our solution + for T in bad_triangles: + del self.triangles[T] + del self.circles[T] + + # Retriangle the hole left by bad_triangles + new_triangles = [] + for (e0, e1, tri_op) in boundary: + # Create a new triangle using point p and edge extremes + T = (idx, e0, e1) + + # Store circumcenter and circumradius of the triangle + self.circles[T] = self.circumcenter(T) + + # Set opposite triangle of the edge as neighbour of T + self.triangles[T] = [tri_op, None, None] + + # Try to set T as neighbour of the opposite triangle + if tri_op: + # search the neighbour of tri_op that use edge (e1, e0) + for i, neigh in enumerate(self.triangles[tri_op]): + if neigh: + if e1 in neigh and e0 in neigh: + # change link to use our new triangle + self.triangles[tri_op][i] = T + + # Add triangle to a temporal list + new_triangles.append(T) + + # Link the new triangles each another + N = len(new_triangles) + for i, T in enumerate(new_triangles): + self.triangles[T][1] = new_triangles[(i+1) % N] # next + self.triangles[T][2] = new_triangles[(i-1) % N] # previous + + def exportTriangles(self): + """Export the current list of Delaunay triangles + """ + # Filter out triangles with any vertex in the extended BBox + return [(a-4, b-4, c-4) + for (a, b, c) in self.triangles if a > 3 and b > 3 and c > 3] + + def exportCircles(self): + """Export the circumcircles as a list of (center, radius) + """ + # Remember to compute circumcircles if not done before + # for t in self.triangles: + # self.circles[t] = self.circumcenter(t) + + # Filter out triangles with any vertex in the extended BBox + # Do sqrt of radius before of return + return [(self.circles[(a, b, c)][0], sqrt(self.circles[(a, b, c)][1])) + for (a, b, c) in self.triangles if a > 3 and b > 3 and c > 3] + + def exportDT(self): + """Export the current set of Delaunay coordinates and triangles. + """ + # Filter out coordinates in the extended BBox + coord = self.coords[4:] + + # Filter out triangles with any vertex in the extended BBox + tris = [(a-4, b-4, c-4) + for (a, b, c) in self.triangles if a > 3 and b > 3 and c > 3] + return coord, tris + + def exportExtendedDT(self): + """Export the Extended Delaunay Triangulation (with the frame vertex). + """ + return self.coords, list(self.triangles) + + def exportVoronoiRegions(self): + """Export coordinates and regions of Voronoi diagram as indexed data. + """ + # Remember to compute circumcircles if not done before + # for t in self.triangles: + # self.circles[t] = self.circumcenter(t) + useVertex = {i: [] for i in range(len(self.coords))} + vor_coors = [] + index = {} + # Build a list of coordinates and one index per triangle/region + for tidx, (a, b, c) in enumerate(sorted(self.triangles)): + vor_coors.append(self.circles[(a, b, c)][0]) + # Insert triangle, rotating it so the key is the "last" vertex + useVertex[a] += [(b, c, a)] + useVertex[b] += [(c, a, b)] + useVertex[c] += [(a, b, c)] + # Set tidx as the index to use with this triangle + index[(a, b, c)] = tidx + index[(c, a, b)] = tidx + index[(b, c, a)] = tidx + + # init regions per coordinate dictionary + regions = {} + # Sort each region in a coherent order, and substitude each triangle + # by its index + for i in range(4, len(self.coords)): + v = useVertex[i][0][0] # Get a vertex of a triangle + r = [] + for _ in range(len(useVertex[i])): + # Search the triangle beginning with vertex v + t = [t for t in useVertex[i] if t[0] == v][0] + r.append(index[t]) # Add the index of this triangle to region + v = t[1] # Choose the next vertex to search + regions[i-4] = r # Store region. + + return vor_coors, regions -- GitLab From a5c30651acd2f8ebd0e5db8bf97f26d77229dec1 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 20 Jun 2021 14:30:58 +0200 Subject: [PATCH 54/61] RIP eigene Kreise --- curvepy/testrange.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/curvepy/testrange.py b/curvepy/testrange.py index df4ba10..2ae817a 100644 --- a/curvepy/testrange.py +++ b/curvepy/testrange.py @@ -87,13 +87,13 @@ class Delaunay_triangulation: def __init__(self, xmin, xmax, ymin, ymax): self.xmin, self.xmax, self.ymin, self.ymax = xmin, xmax, ymin, ymax self.supertriangle = self.create_supertriangle() - self._triangles: List[Triangle] = [self.supertriangle] + self._triangles: List[Triangle] = [*self.supertriangle] self.triangle_queue = [] @property def triangles(self): # We have to remove everything containing vertices of the supertriangle - rem_if = lambda t: any(pt in t.points for pt in self.supertriangle.points) + rem_if = lambda t: any(pt in t.points for pt in self.supertriangle[0].points+self.supertriangle[1].points) return [t for t in self._triangles if not rem_if(t)] def get_all_points(self): @@ -148,7 +148,7 @@ class Delaunay_triangulation: self._triangles.remove(to_rem) def get_second_new_triangle(self, t_new_1, t, pt): - point_not_in_org_t = list(set(t_new_1.points).difference(set(t.points)))[0] + point_not_in_org_t = list(set(t_new_1.points).difference(set(t.points)))[0] # kann immer noch leer sein point_not_in_t_new_1 = list(set(t.points).difference((set(t_new_1.points))))[0] return Triangle(pt, point_not_in_t_new_1, point_not_in_org_t) @@ -190,7 +190,7 @@ class Delaunay_triangulation: self._triangles += [t1_new, t2_new] self.triangle_queue += [t1_new, t2_new] - def create_supertriangle(self) -> Triangle: + def create_supertriangle(self) -> (Triangle, Triangle): # Rectangle CCW A = np.array([self.xmin, self.ymax]) B = np.array([self.xmin, self.ymin]) @@ -205,7 +205,7 @@ class Delaunay_triangulation: B_supert = np.array([self.xmax + (self.xmax - self.xmin) / 2, self.ymin]) C_supert = self.intersect_lines(A_supert, A, B_supert, D) - return Triangle(tuple(A_supert), tuple(B_supert), tuple(C_supert)) + return Triangle(tuple(A), tuple(B), tuple(C)), Triangle(tuple(A), tuple(C), tuple(D)) def intersect_lines(self, p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray): """ @@ -280,7 +280,7 @@ class Delaunay_triangulation: if __name__ == '__main__': pts = [np.array(x) for x in ((2, 3), (6, 5), (3, 7), (8, 3), (5, 1), (8, 8), (-3, -2))] dein_min, dein_max = -5, 5 - n = 7 + n = 20 pts = np.array([[rd.uniform(dein_min, dein_max), rd.uniform(dein_min, dein_max)] for _ in range(n)]) tri = Delaunay(pts) -- GitLab From 243fe0f7fcd4230748d94d4f6918fe841e147c6a Mon Sep 17 00:00:00 2001 From: Lars Quentin Date: Sun, 20 Jun 2021 14:31:27 +0200 Subject: [PATCH 55/61] das isses nun --- curvepy/besser_wirds_nicht.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py index e8588a1..8b296fe 100644 --- a/curvepy/besser_wirds_nicht.py +++ b/curvepy/besser_wirds_nicht.py @@ -155,10 +155,10 @@ class DelaunayTriangulation2D: if __name__ == '__main__': - n = 10 - min, max = -10, 10 + n = 50 + min, max = -100, 100 pts = [(rd.uniform(min, max), rd.uniform(min, max)) for _ in range(n)] - d = DelaunayTriangulation2D(radius=n + 5) # Buffer for rounding errors + d = DelaunayTriangulation2D(radius=max + 5) # Buffer for rounding errors for p in pts: d.add_point(p) @@ -177,7 +177,7 @@ if __name__ == '__main__': axis[1].plot(points[:, 0], points[:, 1], 'o') axis[2].set_title("reference implementation") - d2 = Delaunay2D(radius=n+5) + d2 = Delaunay2D(radius=max+5) for p in pts: d2.addPoint(p) coord, tris = d2.exportDT() -- GitLab From c2943d7ae1d9421768678c1ec511f799a5d351d4 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 20 Jun 2021 15:38:36 +0200 Subject: [PATCH 56/61] vor performanterer Nachbarerstellung --- curvepy/besser_wirds_nicht.py | 90 ++++++++++++++++++++++++++--------- 1 file changed, 68 insertions(+), 22 deletions(-) diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py index 8b296fe..b42da58 100644 --- a/curvepy/besser_wirds_nicht.py +++ b/curvepy/besser_wirds_nicht.py @@ -1,7 +1,7 @@ import numpy as np import random as rd from functools import cached_property -from typing import List, Tuple, Any +from typing import List, Tuple, Any, Set, Dict from scipy.spatial import Delaunay import matplotlib.pyplot as plt @@ -35,6 +35,7 @@ class Triangle: @property def points(self) -> List[Point2D]: + # If it was mutable caching would break return self._points @cached_property @@ -53,6 +54,11 @@ class Triangle: A, B, C = self.points [x1, y1], [x2, y2], [x3, y3] = A, B, C d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) + # xu = (x1 * x1 + y1 * y1) / d * ((y2 - y3) + (x2 * x2 + y2 * y2)) / d * ((y3 - y1) + (x3 * x3 + y3 * y3)) / d * ( + # y1 - y2) / d + # yu = (x1 * x1 + y1 * y1) / d * ((x3 - x2) + (x2 * x2 + y2 * y2)) / d * ((x1 - x3) + (x3 * x3 + y3 * y3)) / d * ( + # x2 - x1) / d + xu = ((x1 * x1 + y1 * y1) * (y2 - y3) + (x2 * x2 + y2 * y2) * (y3 - y1) + (x3 * x3 + y3 * y3) * (y1 - y2)) / d yu = ((x1 * x1 + y1 * y1) * (x3 - x2) + (x2 * x2 + y2 * y2) * (x1 - x3) + (x3 * x3 + y3 * y3) * (x2 - x1)) / d @@ -86,7 +92,8 @@ class Triangle: class DelaunayTriangulation2D: def __init__(self, center: Point2D = (0, 0), radius: float = 1000): - self.supertriangles: List[Triangle] = self._create_supertriangles(center, radius) + self._neighbours = self._create_supertriangles(center, radius) + self.supertriangles: List[Triangle] = [*self.neighbours.keys()] self._triangles: List[Triangle] = [*self.supertriangles] @property @@ -96,6 +103,10 @@ class DelaunayTriangulation2D: remove_if = lambda t: any(pt in t.points for pt in all_points_in_supertriangle) return [t for t in self._triangles if not remove_if(t)] + @property + def neighbours(self): + raise NotImplementedError() + def _create_supertriangles(self, center: Point2D, radius): # Since we have to start with a valid triangulation, we split our allowed range into 2 triangles like that: # x────────────────────────────────┐ @@ -122,7 +133,8 @@ class DelaunayTriangulation2D: lower_left, lower_right, upper_right, upper_left = [tuple(x) for x in base_rectangle] lower_triangle = Triangle(lower_left, lower_right, upper_left) upper_triangle = Triangle(upper_right, upper_left, lower_right) - return [lower_triangle, upper_triangle] + neighbours = {lower_triangle: {upper_triangle}, upper_triangle: {lower_triangle}} + return neighbours def add_point(self, p: Point2D): bad_triangles = [tri for tri in self._triangles if p in tri.circumcircle] @@ -140,12 +152,12 @@ class DelaunayTriangulation2D: @staticmethod def _find_edges_only_used_by_a_single_triangle(triangles, no_duplicates=True): ret = [] - for triangle in triangles: + for t in triangles: others = list(triangles) - others.remove(triangle) - for edge in triangle.edges: - if all(edge not in other.edges for other in others): - ret.append(edge) + others.remove(t) + for e in t.edges: + if all(e not in o.edges for o in others): + ret.append(e) # set() works since all edges are ordered, we don't care about the order and tuples are hashable if no_duplicates: @@ -153,37 +165,71 @@ class DelaunayTriangulation2D: return ret + @property + def voronoi(self): + xs = [(x, y) for x in self.triangles for y in self.triangles if x != y] + print(f"Länge xs: {len(xs)}") + # print("we did it :D") + neighbours = [] + # print("b4 voronoi loop") + for tri1, tri2 in xs: + # Edges are sorted neighbours + if (tri2, tri1) not in neighbours and set(tri1.edges).intersection(tri2.edges): + # print("APPEND!") + neighbours.append((tri1, tri2)) + print("die schleife endet") + neighbours = [(x.circumcircle.center, y.circumcircle.center) for (x, y) in neighbours] + # print("Return") + return neighbours + if __name__ == '__main__': n = 50 min, max = -100, 100 + print("unseres") pts = [(rd.uniform(min, max), rd.uniform(min, max)) for _ in range(n)] d = DelaunayTriangulation2D(radius=max + 5) # Buffer for rounding errors for p in pts: d.add_point(p) - plt.rcParams["figure.figsize"] = (10,30) - figure, axis = plt.subplots(3) + plt.rcParams["figure.figsize"] = (20, 30) + figure, axis = plt.subplots(3, 2) - axis[0].set_title("meins") - for tri in d.triangles: + axis[0, 0].set_title("meins") + print("b4 calc trianglos") + trianglerinos = d.triangles + print("b4 print trianglos") + for tri in trianglerinos: points = np.ravel(tri.points) - axis[0].triplot(points[0::2], points[1::2]) - - axis[1].set_title("scipy") + axis[0, 0].triplot(points[0::2], points[1::2]) + print("b4 calc voronois") + print(f"rawwtf: {len(d._triangles)}") + print(f"wtf: {len(d.triangles)}") + voronois = d.voronoi + print("b4 print voronois") + axis[0, 1].set_title("Voronoi-Schmoronoi sag ich immer!") + for (x, y) in voronois: + axis[0, 1].plot([x[0], y[0]], [x[1], y[1]]) + + print("scipy :)") + axis[2, 0].set_title("scipy") points = np.array([list(x) for x in pts]) scipy_tri = Delaunay(pts) - axis[1].triplot(points[:, 0], points[:, 1], scipy_tri.simplices) - axis[1].plot(points[:, 0], points[:, 1], 'o') + axis[2, 0].triplot(points[:, 0], points[:, 1], scipy_tri.simplices) + axis[2, 0].plot(points[:, 0], points[:, 1], 'o') - axis[2].set_title("reference implementation") - d2 = Delaunay2D(radius=max+5) + print("ugly") + axis[1, 0].set_title("reference implementation") + d2 = Delaunay2D(radius=max + 5) for p in pts: d2.addPoint(p) + print(f"reftris: {d2.exportTriangles()}") coord, tris = d2.exportDT() - my_tris = [(coord[a], coord[b], coord[c]) for a,b,c in tris] + my_tris = [(coord[a], coord[b], coord[c]) for a, b, c in tris] for tri in my_tris: points = np.ravel(tri) - axis[2].triplot(points[0::2], points[1::2]) + axis[1, 0].triplot(points[0::2], points[1::2]) + + axis[1, 1].set_title("ugly but working voronois >:(") - plt.show() \ No newline at end of file + plt.show() -- GitLab From 1febc79cdabe137b1ac4efa5b55748a501a3e444 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 20 Jun 2021 15:56:31 +0200 Subject: [PATCH 57/61] vor performanterer Nachbarerstellung 2 --- curvepy/besser_wirds_nicht.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py index b42da58..0d0f7b6 100644 --- a/curvepy/besser_wirds_nicht.py +++ b/curvepy/besser_wirds_nicht.py @@ -139,6 +139,17 @@ class DelaunayTriangulation2D: def add_point(self, p: Point2D): bad_triangles = [tri for tri in self._triangles if p in tri.circumcircle] + # 1. Walk through full neighbour graph, get subgraphneighburstructure of bad_triangles + # 2. Get boundaries like below, save which triangles the edge belonged (single tri by definition) + # 3. Remove bad_triangles + # 4. + # for edge in boundaries + # for neighbour in neighbours_of_triangle_containing_edge: # (3) + # for other_edge in neighbour: + # if other_edge in bad_triangle_edges: + # diese beiden neuen triangles sind + + # An edge is part of the boundary iff it doesn't is not part of another bad triangle boundaries = self._find_edges_only_used_by_a_single_triangle(bad_triangles) -- GitLab From 74b12c1011f6f1061f3275693e1d9e062333f8d4 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sun, 20 Jun 2021 15:56:36 +0200 Subject: [PATCH 58/61] vor performanterer Nachbarerstellung 2 --- curvepy/besser_wirds_nicht.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py index 0d0f7b6..69475e3 100644 --- a/curvepy/besser_wirds_nicht.py +++ b/curvepy/besser_wirds_nicht.py @@ -133,7 +133,7 @@ class DelaunayTriangulation2D: lower_left, lower_right, upper_right, upper_left = [tuple(x) for x in base_rectangle] lower_triangle = Triangle(lower_left, lower_right, upper_left) upper_triangle = Triangle(upper_right, upper_left, lower_right) - neighbours = {lower_triangle: {upper_triangle}, upper_triangle: {lower_triangle}} + neighbours = {lower_triangle: [upper_triangle], upper_triangle: [lower_triangle]} return neighbours def add_point(self, p: Point2D): -- GitLab From b0dfa7d00b85edaf6afedbfd56048f875daafc07 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 26 Jun 2021 15:10:43 +0200 Subject: [PATCH 59/61] fixes --- curvepy/besser_wirds_nicht.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py index 69475e3..b785380 100644 --- a/curvepy/besser_wirds_nicht.py +++ b/curvepy/besser_wirds_nicht.py @@ -213,14 +213,14 @@ if __name__ == '__main__': for tri in trianglerinos: points = np.ravel(tri.points) axis[0, 0].triplot(points[0::2], points[1::2]) - print("b4 calc voronois") - print(f"rawwtf: {len(d._triangles)}") - print(f"wtf: {len(d.triangles)}") - voronois = d.voronoi - print("b4 print voronois") - axis[0, 1].set_title("Voronoi-Schmoronoi sag ich immer!") - for (x, y) in voronois: - axis[0, 1].plot([x[0], y[0]], [x[1], y[1]]) + # print("b4 calc voronois") + # print(f"rawwtf: {len(d._triangles)}") + # print(f"wtf: {len(d.triangles)}") + # voronois = d.voronoi + # print("b4 print voronois") + # axis[0, 1].set_title("Voronoi-Schmoronoi sag ich immer!") + # for (x, y) in voronois: + # axis[0, 1].plot([x[0], y[0]], [x[1], y[1]]) print("scipy :)") axis[2, 0].set_title("scipy") -- GitLab From 81a28fe82b42003bd4311871fcedc23c54cc89a7 Mon Sep 17 00:00:00 2001 From: Maximilian Winkler Date: Sat, 26 Jun 2021 15:10:49 +0200 Subject: [PATCH 60/61] fixes --- curvepy/besser_wirds_nicht.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py index b785380..9d0391d 100644 --- a/curvepy/besser_wirds_nicht.py +++ b/curvepy/besser_wirds_nicht.py @@ -89,11 +89,14 @@ class Triangle: def __eq__(self, other: Any) -> bool: return isinstance(other, Triangle) and self.points == other.points + def __hash__(self): + return hash(tuple(self.points)) + class DelaunayTriangulation2D: def __init__(self, center: Point2D = (0, 0), radius: float = 1000): self._neighbours = self._create_supertriangles(center, radius) - self.supertriangles: List[Triangle] = [*self.neighbours.keys()] + self.supertriangles: List[Triangle] = [*self._neighbours.keys()] self._triangles: List[Triangle] = [*self.supertriangles] @property @@ -139,15 +142,7 @@ class DelaunayTriangulation2D: def add_point(self, p: Point2D): bad_triangles = [tri for tri in self._triangles if p in tri.circumcircle] - # 1. Walk through full neighbour graph, get subgraphneighburstructure of bad_triangles - # 2. Get boundaries like below, save which triangles the edge belonged (single tri by definition) - # 3. Remove bad_triangles - # 4. - # for edge in boundaries - # for neighbour in neighbours_of_triangle_containing_edge: # (3) - # for other_edge in neighbour: - # if other_edge in bad_triangle_edges: - # diese beiden neuen triangles sind + # An edge is part of the boundary iff it doesn't is not part of another bad triangle -- GitLab From e7bbb681e6d3c48f22d830b8631c4c3296672b5b Mon Sep 17 00:00:00 2001 From: Lars Quentin Date: Sat, 26 Jun 2021 17:49:25 +0200 Subject: [PATCH 61/61] Last update b4 big merge --- curvepy/besser_wirds_nicht.py | 2 +- curvepy/last_run.npy | Bin 184 -> 0 bytes 2 files changed, 1 insertion(+), 1 deletion(-) delete mode 100644 curvepy/last_run.npy diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py index 9d0391d..8d4c971 100644 --- a/curvepy/besser_wirds_nicht.py +++ b/curvepy/besser_wirds_nicht.py @@ -1,7 +1,7 @@ import numpy as np import random as rd from functools import cached_property -from typing import List, Tuple, Any, Set, Dict +from typing import List, Tuple, Any from scipy.spatial import Delaunay import matplotlib.pyplot as plt diff --git a/curvepy/last_run.npy b/curvepy/last_run.npy deleted file mode 100644 index dd8a1227ea976ce63992e8ee8e5d61d18515afeb..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 184 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlWC%^qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I%|ItoUbItsN4WCJcH1_lOZAZ7z%Rw$btNOM5hATdTLA0+no|NsC0fEWOU CL?fdB -- GitLab