Commit 067702c6 authored by lars.quentin's avatar lars.quentin 💬
Browse files

Merge branch 'linear_interpolation_branch' into 'master'

branch differences degraded too much, lets work on master until we have a bigger and more structured code base

See merge request !5
parents 3445025b e7bbb681
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"<CIRCLE: {str(self)}>"
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"<TRIANLGE: {str(self)}>"
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()
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
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()