Skip to content

Instantly share code, notes, and snippets.

@AsyaMatveeva
Last active June 20, 2023 14:41
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save AsyaMatveeva/cd7833d41e4d71bdef3756e6c0c09646 to your computer and use it in GitHub Desktop.
Save AsyaMatveeva/cd7833d41e4d71bdef3756e6c0c09646 to your computer and use it in GitHub Desktop.
import math
import random
def cos(edge, point):
pe1 = [edge[0][0] - point[0], edge[0][1] - point[1]]
pe2 = [edge[1][0] - point[0], edge[1][1] - point[1]]
return round((pe1[0] * pe2[0] + pe1[1] * pe2[1])/math.dist(edge[0], point)/math.dist(edge[1], point), 14)
def sin(edge, point):
pe1 = [edge[0][0] - point[0], edge[0][1] - point[1]]
pe2 = [edge[1][0] - point[0], edge[1][1] - point[1]]
return round(abs(pe1[1] * pe2[0] - pe1[0] * pe2[1])/math.dist(edge[0], point)/math.dist(edge[1], point), 14)
def delaunay_correct(edge, p1, p2):
sin1 = sin(edge, p1)
cos1 = cos(edge, p1)
sin2 = sin(edge, p2)
cos2 = cos(edge, p2)
return sin1 * cos2 + sin2 * cos1 >= 0
def line(p1, p2):
a = p1[1] - p2[1]
b = p2[0] - p1[0]
c = p1[0] * p2[1] - p2[0] * p1[1]
return [a, b, c]
def naive(points):
points = list(set(points))
points.sort()
p0 = points[0]
edges = []
triangulation = {}
max_abs_slope = [0]
for point in points[1:]:
if point[0] == p0[0]:
slope = float('inf')
else:
slope = (point[1] - p0[1])/(point[0]-p0[0])
if abs(slope) > abs(max_abs_slope[0]):
max_abs_slope = [slope, point]
elif slope == max_abs_slope[0]:
max_abs_slope.append(point)
max_abs_slope[0] = p0
for i in range(1, len(max_abs_slope)):
new_edge = tuple([max_abs_slope[i-1], max_abs_slope[i]])
edges.append(new_edge)
triangulation[new_edge] = []
while len(edges) > 0:
edge = edges.pop(0)
if len(triangulation[edge]) == 2:
continue
coeff = line(edge[0], edge[1])
a, b, c = coeff[0], coeff[1], coeff[2]
if len(triangulation[edge]) == 0:
fixed_point = [edge[0][0] - 1, edge[0][1]]
else:
fixed_point = triangulation[edge][0]
fixed_dist = a * fixed_point[0] + b * fixed_point[1] + c
vertex = [2, []]
for point in points:
if (a * point[0] + b * point[1] + c) * fixed_dist < 0:
cosine = cos(edge, point)
angle = cos((edge[1], point), edge[0])
if cosine < vertex[0]:
vertex[0] = cosine
vertex[1] = [[angle, point]]
elif cosine == vertex[0]:
vertex[1].append([angle, point])
if vertex[0] != 2:
on_circle = vertex[1][::]
on_circle.sort()
on_circle.append([1, edge[1]])
triangulation[edge].append(on_circle[0][1])
new_edge = tuple(sorted([edge[0], on_circle[0][1]]))
if new_edge not in triangulation.keys():
triangulation[new_edge] = [on_circle[1][1]]
edges.append(new_edge)
elif on_circle[1][1] not in triangulation[new_edge]:
triangulation[new_edge].append(on_circle[1][1])
for i in range(len(on_circle)-1):
new_edge = tuple(sorted([on_circle[i][1], on_circle[i+1][1]]))
if new_edge not in triangulation.keys():
triangulation[new_edge] = [edge[0]]
edges.append(new_edge)
elif edge[0] not in triangulation[new_edge]:
triangulation[new_edge].append(edge[0])
for i in range(1, len(on_circle)-1):
new_edge = tuple(sorted([edge[0], on_circle[i][1]]))
triangulation[new_edge] = [on_circle[i-1][1], on_circle[i+1][1]]
return triangulation
flips = 0
all_points = []
flipped = []
k = 5
for x in range(k):
for y in range(k):
all_points.append((x, y))
random.shuffle(all_points)
for i1 in range(k**2-5):
for i2 in range(i1+1, k**2-4):
for i3 in range(i2+1, k**2-3):
for i4 in range(i3+1, k**2-2):
for i5 in range(i4+1, k**2-1):
for i6 in range(i5+1, k**2):
points = [all_points[i1], all_points[i2], all_points[i3], all_points[i4], all_points[i5], all_points[i6]]
triangulation = naive(points)
keys = list(triangulation.keys())
for edge in keys:
val = triangulation[edge]
if len(val) == 2 and not delaunay_correct(edge, val[0], val[1]):
flips += 1
flipped.append(points)
print(flips)
for set_of_points in flipped:
print(set_of_points)
import math
import random
from turtle import *
#returns the cosine of an angle at which an edge is visible from a given point
def cos(edge, point):
pe1 = [edge[0][0] - point[0], edge[0][1] - point[1]]
pe2 = [edge[1][0] - point[0], edge[1][1] - point[1]]
return round((pe1[0] * pe2[0] + pe1[1] * pe2[1])/math.dist(edge[0], point)/math.dist(edge[1], point), 14)
#returns the sine of an angle at which an edge is visible from a given point
def sin(edge, point):
pe1 = [edge[0][0] - point[0], edge[0][1] - point[1]]
pe2 = [edge[1][0] - point[0], edge[1][1] - point[1]]
return round(abs(pe1[1] * pe2[0] - pe1[0] * pe2[1])/math.dist(edge[0], point)/math.dist(edge[1], point), 14)
#returns the coefficients of a line (a * x + b * y + c = 0) through two given points
def line(p1, p2):
a = p1[1] - p2[1]
b = p2[0] - p1[0]
c = p1[0] * p2[1] - p2[0] * p1[1]
return [a, b, c]
points = []
for i in range(7):
x = random.randint(-5, 5)
y = random.randint(-5, 5)
points.append((x, y))
points = list(set(points))
points.sort()
p0 = points[0]
edges = [] #queue with edges
triangulation = {} #dictionary with edges as keys, and points connected to them as values
#First step is to find a line throw p0 and another point from points with the greatest absolute value of slope possible. Every two adjacent points from points[] on the line give a new edge in triangulation. For n points there will be (n-1) edges.
max_abs_slope = [0]
for point in points[1:]:
if point[0] == p0[0]: #Since the slope of a vertical line is undefined as the denominator is equal to zero, it is represented as an infinite value.
slope = float('inf')
else:
slope = (point[1] - p0[1])/(point[0]-p0[0])
if abs(slope) > abs(max_abs_slope[0]):
max_abs_slope = [slope, point]
elif slope == max_abs_slope[0]:
max_abs_slope.append(point)
max_abs_slope[0] = p0
for i in range(1, len(max_abs_slope)):
new_edge = tuple([max_abs_slope[i-1], max_abs_slope[i]])
edges.append(new_edge)
triangulation[new_edge] = []
#Second step is to find the point from which an edge is visible at the greatest angle for every edge in edges[]. Since every edge in edges[] is already connected to a point (except for the edges added in the beginning of the algorithm), the new point is searched in a different half-plane.
while len(edges) > 0:
edge = edges.pop(0)
if len(triangulation[edge]) == 2:
continue
coeff = line(edge[0], edge[1])
a, b, c = coeff[0], coeff[1], coeff[2]
if len(triangulation[edge]) == 0:
fixed_point = [edge[0][0] - 1, edge[0][1]]
else:
fixed_point = triangulation[edge][0]
fixed_dist = a * fixed_point[0] + b * fixed_point[1] + c
vertex = [2, []]
for point in points:
if (a * point[0] + b * point[1] + c) * fixed_dist < 0: #fixed_point and new point are in different half-planes in relation to the edge
cosine = cos(edge, point)
angle = cos((edge[1], point), edge[0])
if cosine < vertex[0]:
vertex[0] = cosine
vertex[1] = [[angle, point]]
elif cosine == vertex[0]:
vertex[1].append([angle, point])
if vertex[0] != 2: #vertex[0] == 2 in case of the absence of points in the half-plane
on_circle = vertex[1][::] #this list contains the vertices of the inscribed polygon except edge[0] and edge[1].
on_circle.sort()
on_circle.append([1, edge[1]])
triangulation[edge].append(on_circle[0][1])
new_edge = tuple(sorted([edge[0], on_circle[0][1]]))
if new_edge not in triangulation.keys(): #new_edge is not connected to any points yet.
triangulation[new_edge] = [on_circle[1][1]]
edges.append(new_edge) #new_edge is added to the queue because one more point could still be connected to this edge.
elif on_circle[1][1] not in triangulation[new_edge]: #connect new_edge to on_circle[1][1], if this has not been done already.
triangulation[new_edge].append(on_circle[1][1]) #new_edge is not added to the queue because it is already connected to two points.
for i in range(len(on_circle)-1): #building outer edges of the inscribed polygon.
new_edge = tuple(sorted([on_circle[i][1], on_circle[i+1][1]]))
if new_edge not in triangulation.keys(): #the same check as for the previous new_edge
triangulation[new_edge] = [edge[0]]
edges.append(new_edge)
elif edge[0] not in triangulation[new_edge]:
triangulation[new_edge].append(edge[0])
for i in range(1, len(on_circle)-1): #building inner edges of the inscribed polygon.
new_edge = tuple(sorted([edge[0], on_circle[i][1]]))
triangulation[new_edge] = [on_circle[i-1][1], on_circle[i+1][1]]
for key, val in triangulation.items():
print(key, ":", val)
def rendering(points, triangulation, k):
speed(1000)
up()
color("red")
for point in points:
goto(point[0] * k, point[1] * k)
dot(5)
up()
color("black")
for key in triangulation.keys():
goto(key[0][0] * k, key[0][1] * k)
down()
goto(key[1][0] * k, key[1][1] * k)
up()
exitonclick()
rendering(points, triangulation, 30)
import math
import random
import time
import matplotlib.pyplot as plt
def cos(edge, point):
pe1 = [edge[0][0] - point[0], edge[0][1] - point[1]]
pe2 = [edge[1][0] - point[0], edge[1][1] - point[1]]
return round((pe1[0] * pe2[0] + pe1[1] * pe2[1])/math.dist(edge[0], point)/math.dist(edge[1], point), 14)
def sin(edge, point):
pe1 = [edge[0][0] - point[0], edge[0][1] - point[1]]
pe2 = [edge[1][0] - point[0], edge[1][1] - point[1]]
return round(abs(pe1[1] * pe2[0] - pe1[0] * pe2[1])/math.dist(edge[0], point)/math.dist(edge[1], point), 14)
def delaunay_correct(edge, p1, p2):
sin1 = sin(edge, p1)
cos1 = cos(edge, p1)
sin2 = sin(edge, p2)
cos2 = cos(edge, p2)
return sin1 * cos2 + sin2 * cos1 >= 0
def line(p1, p2):
a = p1[1] - p2[1]
b = p2[0] - p1[0]
c = p1[0] * p2[1] - p2[0] * p1[1]
return [a, b, c]
def naive(points):
points = list(set(points))
points.sort()
p0 = points[0]
edges = []
triangulation = {}
max_abs_slope = [0]
for point in points[1:]:
if point[0] == p0[0]:
slope = float('inf')
else:
slope = (point[1] - p0[1])/(point[0]-p0[0])
if abs(slope) > abs(max_abs_slope[0]):
max_abs_slope = [slope, point]
elif slope == max_abs_slope[0]:
max_abs_slope.append(point)
max_abs_slope[0] = p0
for i in range(1, len(max_abs_slope)):
new_edge = tuple([max_abs_slope[i-1], max_abs_slope[i]])
edges.append(new_edge)
triangulation[new_edge] = []
while len(edges) > 0:
edge = edges.pop(0)
if len(triangulation[edge]) == 2:
continue
coeff = line(edge[0], edge[1])
a, b, c = coeff[0], coeff[1], coeff[2]
if len(triangulation[edge]) == 0:
fixed_point = [edge[0][0] - 1, edge[0][1]]
else:
fixed_point = triangulation[edge][0]
fixed_dist = a * fixed_point[0] + b * fixed_point[1] + c
vertex = [2, []]
for point in points:
if (a * point[0] + b * point[1] + c) * fixed_dist < 0:
cosine = cos(edge, point)
angle = cos((edge[1], point), edge[0])
if cosine < vertex[0]:
vertex[0] = cosine
vertex[1] = [[angle, point]]
elif cosine == vertex[0]:
vertex[1].append([angle, point])
if vertex[0] != 2:
on_circle = vertex[1][::]
on_circle.sort()
on_circle.append([1, edge[1]])
triangulation[edge].append(on_circle[0][1])
new_edge = tuple(sorted([edge[0], on_circle[0][1]]))
if new_edge not in triangulation.keys():
triangulation[new_edge] = [on_circle[1][1]]
edges.append(new_edge)
elif on_circle[1][1] not in triangulation[new_edge]:
triangulation[new_edge].append(on_circle[1][1])
for i in range(len(on_circle)-1):
new_edge = tuple(sorted([on_circle[i][1], on_circle[i+1][1]]))
if new_edge not in triangulation.keys():
triangulation[new_edge] = [edge[0]]
edges.append(new_edge)
elif edge[0] not in triangulation[new_edge]:
triangulation[new_edge].append(edge[0])
for i in range(1, len(on_circle)-1):
new_edge = tuple(sorted([edge[0], on_circle[i][1]]))
triangulation[new_edge] = [on_circle[i-1][1], on_circle[i+1][1]]
return triangulation
Time = []
for n in range(10, 501, 10):
Time.append(0)
for k in range(10):
points = [(random.randint(-1000, 1000), random.randint(-1000, 1000)) for i in range(n)]
start = time.time()
triangulation = naive(points)
end = time.time()
Time[-1] += end - start
fig, ax = plt.subplots()
result = ax.plot([n for n in range(10, 501, 10)], Time, label='naive')
ax.legend()
plt.xlabel('the number of points')
plt.ylabel('time (sec)')
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment