Sweep line Delaunay triangulation algorithm.

correctness test

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 cross_product(vector1, vector2):
    x1, y1 = vector1[1][0] - vector1[0][0], vector1[1][1] - vector1[0][1]
    x2, y2 = vector2[1][0] - vector2[0][0], vector2[1][1] - vector2[0][1]
    return x1 * y2 - x2 * y1


def sweep_line(points):
    points = list(set(points))
    points.sort()
    triangulation = {}
    convex_hull = {}
    p0, p1, p2 = points[0], points[1], points[2]
    convex_hull[p1] = [p0]
    triangulation[(p0, p1)] = []
    prev = p1
    s = 2
    while s < len(points) and (p1[1] - p0[1]) * (points[s][0] - p0[0]) == (points[s][1] - p0[1]) * (p1[0] - p0[0]):
        triangulation[(prev, points[s])] = []
        convex_hull[prev].append(points[s])
        convex_hull[points[s]] = [prev]
        prev = points[s]
        s += 1
    if s == 2:
        triangulation[(p0, p1)] = [p2]
        triangulation[(p0, p2)] = [p1]
        triangulation[(p1, p2)] = [p0]
        if cross_product((p2, p0), (p1, p0)) > 0:
            convex_hull[p0] = [p1, p2]
            convex_hull[p1] = [p2, p0]
            convex_hull[p2] = [p0, p1]
        else:
            convex_hull[p0] = [p2, p1]
            convex_hull[p1] = [p0, p2]
            convex_hull[p2] = [p1, p0]
        s += 1
    else:
        convex_hull[p0] = [points[s-1], p1]
        convex_hull[points[s-1]].append(p0)
        if s < len(points):
            if cross_product((points[s], p0), (points[s], p1)) < 0:
                for key, val in convex_hull.items():
                    convex_hull[key] = val[::-1]
                convex_hull[points[s]] = [p0, points[s-1]]
                convex_hull[p0][1] = points[s]
                convex_hull[points[s-1]][0] = points[s]
            else:
                convex_hull[points[s]] = [points[s-1], p0]
                convex_hull[p0][0] = points[s]
                convex_hull[points[s-1]][1] = points[s]
            edge1, edge2 = (p0, points[s]), (points[s-1], points[s])
            triangulation[edge1] = [p1]
            triangulation[edge2] = [points[s-2]]
            for j in range(1, s - 1):
                triangulation[(points[j], points[s])] = [points[j-1], points[j+1]]
            for k in range(1, s):
                triangulation[(points[k-1], points[k])].append(points[s])
        s += 1    
    
    
    for i in range(s, len(points)):
        point = points[i]
        visible_points = [points[i-1]]
        last = points[i-1]
        right = convex_hull[last][1]
        while right != point and cross_product((last, point), (right, point)) < 0:
            visible_points.append(right)
            last, right = right, convex_hull[right][1]
        last = points[i-1]
        left = convex_hull[last][0]
        while left != point and cross_product((last, point), (left, point)) > 0:
            visible_points = [left] + visible_points[::]
            last, left = left, convex_hull[left][0]
        for elem in visible_points[1:-1]:
            convex_hull.pop(elem)
        convex_hull[visible_points[0]][1] = point
        convex_hull[visible_points[-1]][0] = point
        convex_hull[point] = [visible_points[0], visible_points[-1]]
        edge1, edge2 = tuple(sorted([visible_points[0], point])), tuple(sorted([visible_points[-1], point]))
        triangulation[edge1] = [visible_points[1]]
        triangulation[edge2] = [visible_points[-2]]
        for j in range(1, len(visible_points) - 1):
            triangulation[tuple(sorted([visible_points[j], point]))] = [visible_points[j-1], visible_points[j+1]]
        for k in range(1, len(visible_points)):
            edge = tuple(sorted([visible_points[k-1], visible_points[k]]))
            triangulation[edge].append(point)
            edges = [[edge, point]]
            while len(edges) > 0:
                cur = edges.pop(0)
                cur_edge = cur[0]
                if len(triangulation[cur_edge]) == 2 and not delaunay_correct(cur_edge, triangulation[cur_edge][0], triangulation[cur_edge][1]):
                    e1, e2, p1 = cur_edge[0], cur_edge[1], cur[1]
                    if triangulation[cur_edge][0] == p1:
                        p2 = triangulation[cur_edge][1]
                    else:
                        p2 = triangulation[cur_edge][0]
                    triangulation[tuple(sorted([p1, p2]))] = [e1, e2]
                    triangulation.pop(cur_edge)
                    triangulation[tuple(sorted([e1, p2]))].remove(e2)
                    triangulation[tuple(sorted([e1, p2]))].append(p1)
                    triangulation[tuple(sorted([e1, p1]))].remove(e2)
                    triangulation[tuple(sorted([e1, p1]))].append(p2)
                    triangulation[tuple(sorted([e2, p1]))].remove(e1)
                    triangulation[tuple(sorted([e2, p1]))].append(p2)
                    triangulation[tuple(sorted([e2, p2]))].remove(e1)
                    triangulation[tuple(sorted([e2, p2]))].append(p1)
                    edges.append([tuple(sorted([e1, p2])), p1])
                    edges.append([tuple(sorted([e2, p2])), p1])
    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 = sweep_line(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)

performance test

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 cross_product(vector1, vector2):
    x1, y1 = vector1[1][0] - vector1[0][0], vector1[1][1] - vector1[0][1]
    x2, y2 = vector2[1][0] - vector2[0][0], vector2[1][1] - vector2[0][1]
    return x1 * y2 - x2 * y1


def sweep_line(points):
    points = list(set(points))
    points.sort()
    triangulation = {}
    convex_hull = {}
    p0, p1, p2 = points[0], points[1], points[2]
    convex_hull[p1] = [p0]
    triangulation[(p0, p1)] = []
    prev = p1
    s = 2
    while s < len(points) and (p1[1] - p0[1]) * (points[s][0] - p0[0]) == (points[s][1] - p0[1]) * (p1[0] - p0[0]):
        triangulation[(prev, points[s])] = []
        convex_hull[prev].append(points[s])
        convex_hull[points[s]] = [prev]
        prev = points[s]
        s += 1
    if s == 2:
        triangulation[(p0, p1)] = [p2]
        triangulation[(p0, p2)] = [p1]
        triangulation[(p1, p2)] = [p0]
        if cross_product((p2, p0), (p1, p0)) > 0:
            convex_hull[p0] = [p1, p2]
            convex_hull[p1] = [p2, p0]
            convex_hull[p2] = [p0, p1]
        else:
            convex_hull[p0] = [p2, p1]
            convex_hull[p1] = [p0, p2]
            convex_hull[p2] = [p1, p0]
        s += 1
    else:
        convex_hull[p0] = [points[s-1], p1]
        convex_hull[points[s-1]].append(p0)
        if s < len(points):
            if cross_product((points[s], p0), (points[s], p1)) < 0:
                for key, val in convex_hull.items():
                    convex_hull[key] = val[::-1]
                convex_hull[points[s]] = [p0, points[s-1]]
                convex_hull[p0][1] = points[s]
                convex_hull[points[s-1]][0] = points[s]
            else:
                convex_hull[points[s]] = [points[s-1], p0]
                convex_hull[p0][0] = points[s]
                convex_hull[points[s-1]][1] = points[s]
            edge1, edge2 = (p0, points[s]), (points[s-1], points[s])
            triangulation[edge1] = [p1]
            triangulation[edge2] = [points[s-2]]
            for j in range(1, s - 1):
                triangulation[(points[j], points[s])] = [points[j-1], points[j+1]]
            for k in range(1, s):
                triangulation[(points[k-1], points[k])].append(points[s])
        s += 1    
    
    
    for i in range(s, len(points)):
        point = points[i]
        visible_points = [points[i-1]]
        last = points[i-1]
        right = convex_hull[last][1]
        while right != point and cross_product((last, point), (right, point)) < 0:
            visible_points.append(right)
            last, right = right, convex_hull[right][1]
        last = points[i-1]
        left = convex_hull[last][0]
        while left != point and cross_product((last, point), (left, point)) > 0:
            visible_points = [left] + visible_points[::]
            last, left = left, convex_hull[left][0]
        for elem in visible_points[1:-1]:
            convex_hull.pop(elem)
        convex_hull[visible_points[0]][1] = point
        convex_hull[visible_points[-1]][0] = point
        convex_hull[point] = [visible_points[0], visible_points[-1]]
        edge1, edge2 = tuple(sorted([visible_points[0], point])), tuple(sorted([visible_points[-1], point]))
        triangulation[edge1] = [visible_points[1]]
        triangulation[edge2] = [visible_points[-2]]
        for j in range(1, len(visible_points) - 1):
            triangulation[tuple(sorted([visible_points[j], point]))] = [visible_points[j-1], visible_points[j+1]]
        for k in range(1, len(visible_points)):
            edge = tuple(sorted([visible_points[k-1], visible_points[k]]))
            triangulation[edge].append(point)
            edges = [[edge, point]]
            while len(edges) > 0:
                cur = edges.pop(0)
                cur_edge = cur[0]
                if len(triangulation[cur_edge]) == 2 and not delaunay_correct(cur_edge, triangulation[cur_edge][0], triangulation[cur_edge][1]):
                    e1, e2, p1 = cur_edge[0], cur_edge[1], cur[1]
                    if triangulation[cur_edge][0] == p1:
                        p2 = triangulation[cur_edge][1]
                    else:
                        p2 = triangulation[cur_edge][0]
                    triangulation[tuple(sorted([p1, p2]))] = [e1, e2]
                    triangulation.pop(cur_edge)
                    triangulation[tuple(sorted([e1, p2]))].remove(e2)
                    triangulation[tuple(sorted([e1, p2]))].append(p1)
                    triangulation[tuple(sorted([e1, p1]))].remove(e2)
                    triangulation[tuple(sorted([e1, p1]))].append(p2)
                    triangulation[tuple(sorted([e2, p1]))].remove(e1)
                    triangulation[tuple(sorted([e2, p1]))].append(p2)
                    triangulation[tuple(sorted([e2, p2]))].remove(e1)
                    triangulation[tuple(sorted([e2, p2]))].append(p1)
                    edges.append([tuple(sorted([e1, p2])), p1])
                    edges.append([tuple(sorted([e2, p2])), p1])      
    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 = sweep_line(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='sweep line')
ax.legend()
plt.xlabel('the number of points')
plt.ylabel('time (sec)')
plt.show()

Sweep line Delaunay triangulation algorithm

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)


#checks if the triangulation of four points is a delaunay triangulation
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 cross_product(vector1, vector2):
    x1, y1 = vector1[1][0] - vector1[0][0], vector1[1][1] - vector1[0][1]
    x2, y2 = vector2[1][0] - vector2[0][0], vector2[1][1] - vector2[0][1]
    return x1 * y2 - x2 * y1


points = []
for i in range(6):
    x = random.randint(-7, 7)
    y = random.randint(-7, 7)
    points.append((x, y))
points = list(set(points))
points.sort()
triangulation = {}
convex_hull = {}
p0, p1, p2 = points[0], points[1], points[2]
convex_hull[p1] = [p0]
triangulation[(p0, p1)] = []
prev = p1
s = 2
#building either a triangle, or (in case first three points are on the same line) a line through first three points
while s < len(points) and (p1[1] - p0[1]) * (points[s][0] - p0[0]) == (points[s][1] - p0[1]) * (p1[0] - p0[0]):
    triangulation[(prev, points[s])] = []
    convex_hull[prev].append(points[s])
    convex_hull[points[s]] = [prev]
    prev = points[s]
    s += 1
if s == 2: #if this is true, first three points are not on the same line
    triangulation[(p0, p1)] = [p2]
    triangulation[(p0, p2)] = [p1]
    triangulation[(p1, p2)] = [p0]
    if cross_product((p2, p0), (p1, p0)) > 0:
        convex_hull[p0] = [p1, p2]
        convex_hull[p1] = [p2, p0]
        convex_hull[p2] = [p0, p1]
    else:
        convex_hull[p0] = [p2, p1]
        convex_hull[p1] = [p0, p2]
        convex_hull[p2] = [p1, p0]
    s += 1
else: #the case in which first n points for n > 2 are on the same line is processed separately
    convex_hull[p0] = [points[s-1], p1]
    convex_hull[points[s-1]].append(p0)
    if s < len(points):
        if cross_product((points[s], p0), (points[s], p1)) < 0: #points[s] is above the line
            for key, val in convex_hull.items():
                convex_hull[key] = val[::-1]
            convex_hull[points[s]] = [p0, points[s-1]]
            convex_hull[p0][1] = points[s]
            convex_hull[points[s-1]][0] = points[s]
        else: #points[s] is below the line
            convex_hull[points[s]] = [points[s-1], p0]
            convex_hull[p0][0] = points[s]
            convex_hull[points[s-1]][1] = points[s]
        edge1, edge2 = (p0, points[s]), (points[s-1], points[s])
        triangulation[edge1] = [p1]
        triangulation[edge2] = [points[s-2]]
        for j in range(1, s - 1):
            triangulation[(points[j], points[s])] = [points[j-1], points[j+1]]
        for k in range(1, s):
            triangulation[(points[k-1], points[k])].append(points[s])
    s += 1


for i in range(s, len(points)):
    point = points[i]
    visible_points = [points[i-1]]
    last = points[i-1]
    right = convex_hull[last][1]
    #finding visible points counterclockwise
    while right != point and cross_product((last, point), (right, point)) < 0:
        visible_points.append(right)
        last, right = right, convex_hull[right][1]
    last = points[i-1]
    left = convex_hull[last][0]
    #finding visible points while moving clockwise
    while left != point and cross_product((last, point), (left, point)) > 0:
        visible_points = [left] + visible_points[::]
        last, left = left, convex_hull[left][0]
    #update convex_hull
    for elem in visible_points[1:-1]:
        convex_hull.pop(elem)
    convex_hull[visible_points[0]][1] = point
    convex_hull[visible_points[-1]][0] = point
    convex_hull[point] = [visible_points[0], visible_points[-1]]
    #add new edges
    edge1, edge2 = tuple(sorted([visible_points[0], point])), tuple(sorted([visible_points[-1], point]))
    triangulation[edge1] = [visible_points[1]]
    triangulation[edge2] = [visible_points[-2]]
    for j in range(1, len(visible_points) - 1): #building edges that connect current point and visible points.
        triangulation[tuple(sorted([visible_points[j], point]))] = [visible_points[j-1], visible_points[j+1]]
    for k in range(1, len(visible_points)): #building visible edges.
        edge = tuple(sorted([visible_points[k-1], visible_points[k]]))
        triangulation[edge].append(point)
        edges = [[edge, point]] #the edges of the polygon connected to the point stored in the sublist do not require rebuilding. 
        #sequential rebuilding of the edges
        while len(edges) > 0:
            cur = edges.pop(0)
            cur_edge = cur[0]
            if len(triangulation[cur_edge]) == 2 and not delaunay_correct(cur_edge, triangulation[cur_edge][0], triangulation[cur_edge][1]):
                e1, e2, p1 = cur_edge[0], cur_edge[1], cur[1]
                if triangulation[cur_edge][0] == p1: #edges connected to p1 do not require rebuilding.
                    p2 = triangulation[cur_edge][1]
                else:
                    p2 = triangulation[cur_edge][0]
                #flipping edges
                triangulation[tuple(sorted([p1, p2]))] = [e1, e2] #edge through two points which cur_edge does not contain is added.
                triangulation.pop(cur_edge) #cur_edge is removed
                #since the diagonal in the polygon is changed each outer edge is being connected to another vertex of the polygon.
                triangulation[tuple(sorted([e1, p2]))].remove(e2)
                triangulation[tuple(sorted([e1, p2]))].append(p1)
                triangulation[tuple(sorted([e1, p1]))].remove(e2)
                triangulation[tuple(sorted([e1, p1]))].append(p2)
                triangulation[tuple(sorted([e2, p1]))].remove(e1)
                triangulation[tuple(sorted([e2, p1]))].append(p2)
                triangulation[tuple(sorted([e2, p2]))].remove(e1)
                triangulation[tuple(sorted([e2, p2]))].append(p1)
                #edges not connected to p2 might need to be rebuilt
                edges.append([tuple(sorted([e1, p2])), p1])
                edges.append([tuple(sorted([e2, p2])), p1])
                

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)