diff --git a/curvepy/besser_wirds_nicht.py b/curvepy/besser_wirds_nicht.py new file mode 100644 index 0000000000000000000000000000000000000000..8d4c971f374773365cd86ce6ff3fe041fe47cb4f --- /dev/null +++ b/curvepy/besser_wirds_nicht.py @@ -0,0 +1,241 @@ +import numpy as np +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 + +from curvepy.reference_implementation import Delaunay2D + +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]: + # If it was mutable caching would break + 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) / 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 + + 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 + + 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._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)] + + @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────────────────────────────────┐ + # 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 + + # 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] + 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]} + return neighbours + + def add_point(self, p: Point2D): + bad_triangles = [tri for tri in self._triangles if p in tri.circumcircle] + + + + + # 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) + # 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 t in triangles: + others = list(triangles) + 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: + ret = list(set(ret)) + + 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"] = (20, 30) + figure, axis = plt.subplots(3, 2) + + 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, 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[2, 0].triplot(points[:, 0], points[:, 1], scipy_tri.simplices) + axis[2, 0].plot(points[:, 0], points[:, 1], 'o') + + 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] + for tri in my_tris: + points = np.ravel(tri) + axis[1, 0].triplot(points[0::2], points[1::2]) + + axis[1, 1].set_title("ugly but working voronois >:(") + + plt.show() diff --git a/curvepy/example_algo.py b/curvepy/example_algo.py new file mode 100644 index 0000000000000000000000000000000000000000..9f60ebb40b9570121d27961f566f4ee7842f8420 --- /dev/null +++ b/curvepy/example_algo.py @@ -0,0 +1,213 @@ +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): + 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) + + +class dirichlet_tessellation: + + def __init__(self): + 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.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) + + # 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, 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_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) + + 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)) + + def calculate_collisions(self, p, nearest_p: h_tuple, garbage_heap): + 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] + + 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): + # 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]: + 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): + 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): + 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[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) == set(old_pair)]) + self.valid_triangulation.add(new_pair) + + @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]) + # 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) + # 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]]) + # 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, p3, p4) + + @staticmethod + 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]) + 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 + + +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))] + # d = dirichlet_tessellation() + + # 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) + # 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 diff --git a/curvepy/linear_interpolation.py b/curvepy/linear_interpolation.py index 995148c7f7e060dd19ff7d26b9cb6dd181e175df..8b8d9b906e243a9b5c607a45a46a6f2522323224 100644 --- a/curvepy/linear_interpolation.py +++ b/curvepy/linear_interpolation.py @@ -1,8 +1,11 @@ +from plistlib import Data + import numpy as np 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: @@ -47,7 +50,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,88 +100,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: """ Class for creating a 2D or 3D Polygon. @@ -253,34 +174,211 @@ 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. + """ + 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])) + + @staticmethod + def area(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> float: """ - 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. + 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 ---------- - xs: np.ndArray - the points that may be plotted + 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 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 + 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 + ---------- + p: np.ndarray + Additional point that should be on the same plane as the triangle. + + Returns + ------- + np.ndarray: + 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() + 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 check_points_for_area_calc(self, p): + """ + 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 + Additional point that has to be on the same plane as the triangle. + + Returns + ------- + np.ndarrays: + The triangle points and p so that cramer's rule can be used. """ - 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() + 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]) + else: + 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]) + + +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.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? + + 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: @@ -309,8 +407,9 @@ 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]) + # 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]]) @@ -324,10 +423,15 @@ 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))) + + test_tes = Dirichlet_tessellation(np.array([[0, 0], [1, 2], [-2, -3], [-3, -3]])) if __name__ == "__main__": diff --git a/curvepy/reference_implementation.py b/curvepy/reference_implementation.py new file mode 100644 index 0000000000000000000000000000000000000000..7fc393985edb177f1408442efce45f0b371b63f0 --- /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 diff --git a/curvepy/testrange.py b/curvepy/testrange.py new file mode 100644 index 0000000000000000000000000000000000000000..2ae817acb3df7ccab34b23f85d9cf7bfa5a86ab4 --- /dev/null +++ b/curvepy/testrange.py @@ -0,0 +1,304 @@ +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: + 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): + """ + :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=np.array([xu, yu]), radius=R) + + def get_farthest_point_away_and_nearest_line_to(self, pt): + points = self.points.copy() + 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 + + def __contains__(self, pt: Tuple[float, float]): + """ + 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) + 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): + """ + See: https://www.mathopenref.com/coordtrianglearea.html + """ + return abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2.0) + + def __str__(self): + return str(self.points) + + 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): + self.center = center + self.radius = radius + + 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): + 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[0].points+self.supertriangle[1].points) + 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], [])) + + 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) + 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(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)))] + 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: + 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) + + 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] # 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) + + 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 + 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): + 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 + 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, 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(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): + """ + 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]) + + @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 + + def add_new_triangles(self, pt, t): + 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: + 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))] + dein_min, dein_max = -5, 5 + n = 20 + 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)) + + 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) diff --git a/curvepy/testrange2.py b/curvepy/testrange2.py new file mode 100644 index 0000000000000000000000000000000000000000..66cb4974788d6dfb2b4b78e57154ffdb54d80d5d --- /dev/null +++ b/curvepy/testrange2.py @@ -0,0 +1,154 @@ +import numpy as np +from math import sqrt +from typing import Tuple, List, Dict, Optional + +PointIndex = int +Triangle = Tuple[PointIndex, PointIndex, PointIndex] +Neighbours = List[Optional[Triangle]] + +Radius = float +Center = Tuple[float, float] +Circle = Tuple[Center, Radius] + + +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 + """ + 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 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] + + # 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)] + + # 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]] = [] + + # Wir starten bei einem "random" Triangle mit dessen ersten Edge + T = bad_trinagles[0] + edge = 0 + + # 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)) + + edge = (edge + 1) % 3 + + # 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 + +