https://mathworld.wolfram.com/SteinerCircumellipse.html

steiner.py

# Author: Pierce Brooks
# https://gist.github.com/PierceLBrooks/d180e8445e3c59fe757e99c932257e9e

import os
import sys
import copy
import math
import cmath
import random
import logging
import functools
import traceback
import numpy as np
import pygame as pg
import pscreen as ps
from datetime import datetime
from scipy.spatial.transform import Rotation as R

def segment(p, q, r):
    if ((q[0] <= max(p[0], r[0])) and (q[0] >= min(p[0], r[0])) and (q[1] <= max(p[1], r[1])) and (q[1] >= min(p[1], r[1]))):
        return True         
    return False
    
def orientation(p, q, r):
    value = (((q[1]-p[1])*(r[0]-q[0]))-((q[0]-p[0])*(r[1]-q[1])))
    if (value == 0):
        return 0
    if (value > 0):
        return 1
    return 2

def intersect(p1, q1, p2, q2):
    o1 = orientation(p1, q1, p2)
    o2 = orientation(p1, q1, q2)
    o3 = orientation(p2, q2, p1)
    o4 = orientation(p2, q2, q1)
    if ((o1 != o2) and (o3 != o4)):
        return True
    if ((o1 == 0) and (segment(p1, p2, q1))):
        return True
    if ((o2 == 0) and (segment(p1, q2, q1))):
        return True
    if ((o3 == 0) and (segment(p2, p1, q2))):
        return True
    if ((o4 == 0) and (segment(p2, q1, q2))):
        return True
    return False
    
def contain(points, p):
    n = len(points)
    if (n < 3):
        return False
    extreme = (10000, p[1])
    count = 0
    i = 0
    while (True):
        next = (i+1)%n
        if (intersect(points[i], points[next], p, extreme)):
            if (orientation(points[i], p, points[next]) == 0):
                return segment(points[i], p, points[next])
            count += 1
        i = next
        if (i == 0):
            break
    if (count % 2 == 1):
        return True
    return False
    
def dot(left, right):
    return (left[0]*right[0])+(left[1]*right[1])

def distance(left, right):
    return (((left[0]-right[0])**2.0)+((left[1]-right[1])**2.0))**0.5
    
def direction(left, right):
    return math.atan2(right[1]-left[1], right[0]-left[0])
    
def magnitude(vector):
    result = (vector[0]**2.0)+(vector[1]**2.0)
    return result**0.5
    
def normalize(vector):
    temp = magnitude(vector)
    return tuple([vector[0]/temp, vector[1]/temp])
    
def determinant(a, b):
    return (a[0]*b[1])-(a[1]*b[0])

def intersection(a, b):
    if not (intersect(a[0], a[1], b[0], b[1])):
        return None
    xdiff = (a[0][0]-a[1][0], b[0][0]-b[1][0])
    ydiff = (a[0][1]-a[1][1], b[0][1]-b[1][1])
    div = determinant(xdiff, ydiff)
    if (div == 0):
       return None
    d = (determinant(a[0], a[1]), determinant(b[0], b[1]))
    x = determinant(d, xdiff)/div
    y = determinant(d, ydiff)/div
    return (x, y)

def trisides(a, b, c):
    sides = []
    sides.append(distance(b, c))
    sides.append(distance(c, a))
    sides.append(distance(a, b))
    return sides
    
def triangles(sides):
    angles = [0.0, 0.0, 0.0]
    if (len(sides) < 3):
        return angles
    ab = sides[0]
    bc = sides[1]
    ca = sides[2]
    check = None
    index = 0
    temp = 0.0
    total = 180.0
    king = 0
    major = 1
    minor = 2
    if (ab > bc):
        if (ab > ca):
            index = 0
            king = 0
            if (bc > ca):
                check = True
                major = 1
                minor = 2
            else:
                check = False
                major = 2
                minor = 1
        else:
            check = True
            index = 2
            king = 2
            major = 0
            minor = 1
    else:
        if (bc > ca):
            index = 1
            king = 0
            if (ab > ca):
                check = True
                major = 0
                minor = 2
            else:
                char = False
                major = 2
                minor = 0
        else:
            check = False
            index = 2
            king = 2
            major = 1
            minor = 0
    temp = ((sides[major])**2.0)+((sides[minor])**2.0)
    temp -= (sides[king])**2.0
    temp /= sides[major]*sides[minor]*2.0
    temp = math.acos(temp)
    if (index == 0):
        angles[0] = temp
    elif (index == 1):
        angles[1] = temp
    elif (index == 2):
        angles[2] = temp
    else:
        return angles
    total -= temp
    temp = sides[major]*math.sin(temp)
    temp /= sides[king]
    temp = math.asin(temp)
    while (True):
        if (index%3 == 0):
            if (check):
                angles[1] = temp
            else:
                angles[2] = temp
        elif (index%3 == 1):
            if (check):
                angles[0] = temp
            else:
                angles[2] = temp
        elif (index%3 == 2):
            if (check):
                angles[0] = temp
            else:
                angles[1] = temp
        else:
            return angles
        if (index < 3):
            index += 3
            total -= temp
            temp = total
            if (check):
                check = False
            else:
                check = True
        else:
            break
    return angles
    
def centroid(points):
    if (points == None):
        return (0.0, 0.0)
    if (len(points) == 0):
        return (0.0, 0.0)
    center = complex(0.0, 0.0)
    for i in range(len(points)):
        point = points[i]
        center += complex(point[0], point[1])
    center /= float(len(points))
    return tuple([float(center.real), float(center.imag)])
    
def compare(left, right):
    if (left[0] > right[0]):
        return 1
    if (left[0] < right[0]):
        return -1
    return 0

def sort(points):
    if (points == None):
        return []
    angles = []
    center = centroid(points)
    for i in range(len(points)):
        point = points[i]
        angle = direction(center, point)
        angles.append([angle, i])
    sorting = sorted(angles, key=functools.cmp_to_key(compare))
    result = []
    for i in range(len(sorting)):
        result.append(points[sorting[i][1]])
    return result
    
def interpolate(left, right, key):
    return ((left[0]*(1.0-key))+(right[0]*key), (left[1]*(1.0-key))+(right[1]*key))
    
def split(points, lines):
    if ((len(points) < 2) or (lines < 2)):
        return []
    splits = []
    shares = []
    distances = []
    progresses = []
    progress = 0
    count = 0
    last = 0
    all = 0
    total = 0.0
    for i in range(len(points)):
        if (i == 0):
            continue
        distances.append(distance(points[i-1], points[i]))
        total += distances[len(distances)-1]
    for i in range(len(distances)):
        if not (total > 0.0):
            break
        share = distances[i]/total
        if not (i == 0):
            share += shares[len(shares)-1]
        shares.append(share)
    splits.append(points[0])
    if (lines > 2):
        for i in range(lines):
            line = float(i+1)/float(lines)
            line *= total
            count += 1
            all += 1
            if (line > shares[progress]):
                last = progress
                while (line > shares[progress]):
                    progress += 1
                    if (progress == len(shares)):
                        break
                for i in range(progress-last):
                    progresses.append(int(count/(progress-last)))
                count = 0
                if (progress == len(shares)):
                    break
        if (len(progresses) == 0):
            progresses.append(lines)
            all = lines
        if (all < lines):
            all = int((lines-all)/len(progresses))
            for i in range(len(progresses)):
                progresses[i] += all
        for i in range(len(progresses)):
            for j in range(progresses[i]):
                splits.append(interpolate(points[i], points[i+1], float(j)/float(progresses[i])))
            if (len(splits) == lines):
                break
    splits.append(points[len(points)-1])
    return splits

def circumellipse(triangle, lines, polygon):
    result = []
    adjust = []
    if (len(triangle) < 3):
        return result
    try:
        a = triangle[0]
        b = triangle[1]
        c = triangle[2]
        sides = trisides(a, b, c)
        angles = triangles(sides)
        complexA = complex(a[0], a[1])
        complexB = complex(b[0], b[1])
        complexC = complex(c[0], c[1])
        center = (complexA+complexB+complexC)/3.0
        #print(str(center))
        pi = math.pi
        z = (sides[0]**4.0)+(sides[1]**4.0)+(sides[2]**4.0)
        z -= ((sides[0])**2.0)*((sides[1])**2.0)
        z -= ((sides[1])**2.0)*((sides[2])**2.0)
        z -= ((sides[2])**2.0)*((sides[0])**2.0)
        z = z**0.5
        temp = z*2.0
        semimajor = (sides[0]**2.0)+(sides[1]**2.0)+(sides[2]**2.0)
        semiminor = semimajor
        semimajor += temp
        semiminor -= temp
        temp = float(1.0/3.0)
        semimajor = temp*(semimajor**0.5)
        semiminor = temp*(semiminor**0.5)
        focusFirst = center**2.0
        focusFirst -= ((complexA*complexB)+(complexA*complexC)+(complexB*complexC))/3.0
        focusFirst = (focusFirst**0.5)*2.0
        focusSecond = -focusFirst
        focusFirst += center
        focusSecond += center
        dir = direction(tuple([float(focusFirst.real), float(focusFirst.imag)]), tuple([float(focusSecond.real), float(focusSecond.imag)]))*(180.0/pi)
        r = R.from_euler('z', dir, degrees=True)
        for i in range(lines):
            index = float(i*2)
            temp = ((pi/float(lines))*index)-(pi*0.5)
            point = []
            point.append(semimajor*math.cos(temp))
            point.append(semiminor*math.sin(temp))
            point.append(0.0)
            point = r.apply(point)
            result.append(tuple([float(point[0])+float(center.real), float(point[1])+float(center.imag)]))
            adjust.append(result[len(result)-1])
        if not (polygon == None):
            center = tuple([float(center.real), float(center.imag)])
            loop = 0
            exceptions = []
            while (True):
                lastPolygon = -1
                last = -1
                count = 0
                inside = False
                change = False
                for i in range(len(result)*2):
                    if (contain(polygon, result[i%len(result)])):
                        if not (inside):
                            if ((last > -1) and (lastPolygon > -1)):
                                for j in range(len(polygon)):
                                    if (j in exceptions):
                                        continue
                                    point = intersection(tuple([polygon[j], polygon[(j+1)%len(polygon)]]), tuple([result[(i-1)%len(result)], result[i%len(result)]]))
                                    if (point == None):
                                        continue
                                    other = intersection(tuple([polygon[lastPolygon], polygon[(lastPolygon+1)%len(polygon)]]), tuple([result[last%len(result)], result[(last+1)%len(result)]]))
                                    if (other == None):
                                        continue
                                    others = []
                                    if (lastPolygon <= j):
                                        for k in range((j+1)-lastPolygon):
                                            if (k == 0):
                                                others.append(other)
                                                continue
                                            others.append(polygon[(k+lastPolygon)%len(polygon)])
                                        others.append(point)
                                    else:
                                        for k in range((j+len(polygon)+1)-lastPolygon):
                                            if (k == 0):
                                                others.append(other)
                                                continue
                                            others.append(polygon[(k+lastPolygon)%len(polygon)])
                                        others.append(point)
                                    #others = list(reversed(others))
                                    splits = split(others, (i+1)-last)
                                    if (len(splits) > 0):
                                        point = result[i%len(result)]
                                        other = result[last%len(result)]
                                        exceptions.append(j)
                                        if (lastPolygon <= j):
                                            for k in range((j+1)-lastPolygon):
                                                if not ((k+lastPolygon)%len(polygon) in exceptions):
                                                    exceptions.append((k+lastPolygon)%len(polygon))
                                        else:
                                            for k in range((j+len(polygon)+1)-lastPolygon):
                                                if not ((k+lastPolygon)%len(polygon) in exceptions):
                                                    exceptions.append((k+lastPolygon)%len(polygon))
                                        for k in range(len(splits)):
                                            adjust[(k+last)%len(result)] = splits[k]
                                        if (last%len(result) < i%len(result)):
                                            adjust = adjust[:(last%len(result))]+[other]+adjust[(last%len(result)):(i%len(result))]+[point]+adjust[(i%len(result)):]
                                        else:
                                            adjust = adjust[:(i%len(result))]+[point]+adjust[(i%len(result)):(last%len(result))]+[other]+adjust[(last%len(result)):]
                                        change = True
                                    break
                            inside = True
                            lastPolygon = -1
                        last = i
                        count += 1
                    else:
                        if (inside):
                            for j in range(len(polygon)):
                                if (j in exceptions):
                                    continue
                                point = intersection(tuple([polygon[j], polygon[(j+1)%len(polygon)]]), tuple([result[(i-1)%len(result)], result[i%len(result)]]))
                                if (point == None):
                                    continue
                                lastPolygon = j
                                break
                            if (lastPolygon > -1):
                                inside = False
                    if (change):
                        break
                if not (change):
                    if ((count == 0) and (loop == 0)):
                        result = polygon
                    break
                adjust = sort(adjust)
                result = copy.deepcopy(adjust)
                loop += 1
                if (len(exceptions) == len(polygon)):
                    break
    except Exception as exception:
        now = datetime.now()
        print(now.strftime("%d-%m-%Y@%H:%M:%S"))
        logging.error(traceback.format_exc())
    return result

def run(size):
    lines = 50
    radius = 3.0
    triangle = []
    ellipse = []
    polygon = []
    buttonsNext = [False, False, False]
    buttonsPrevious = [False, False, False]
    triangleTarget = -1
    polygonTarget = -1
    triangleColor = (255, 0, 0, 255)
    ellipseColor = (0, 255, 0, 255)
    polygonColor = (0, 0, 255, 128)
    backgroundColor = (128, 128, 128)
    move = False
    quit = False
    next = None
    previous = None
    mouse = None
    triangleCenter = None
    ellipseCenter = None
    polygonCenter = None
    ps.LoadScreen(size)
    while (True):
        mouse = ps.MouseGetPosition()
        mouse = tuple([float(mouse[0]), float(mouse[1])])
        buttonsPrevious[0] = buttonsNext[0]
        buttonsPrevious[1] = buttonsNext[1]
        buttonsPrevious[2] = buttonsNext[2]
        buttonsNext[0] = ps.MouseGetButtonL()
        buttonsNext[1] = ps.MouseGetButtonM()
        buttonsNext[2] = ps.MouseGetButtonR()
        if (ps.KeyIsPressed("escape")):
            print("quit")
            quit = True
        if ((triangleTarget < 0) or (polygonTarget < 0)):
            move = False
        ps.ClearScreen(backgroundColor)
        for i in range(len(triangle)):
            next = triangle[i]
            dist = distance(next, mouse)
            #print(str(i)+" "+str(dist))
            if ((triangleTarget == i) or (dist < radius*2.0)):
                if (buttonsNext[1]):
                    if (len(triangle) == 0):
                        triangle = []
                    elif (i == 0):
                        triangle = triangle[1:]
                    elif (i < len(triangle)-1):
                        triangle = triangle[0:i]+triangle[(i+1):]
                    else:
                        triangle = triangle[0:i]
                    ellipse = []
                    ellipseCenter = None
                    triangleCenter = None
                    break
                elif (buttonsNext[0]):
                    if not ((triangleTarget > -1) and not (triangleTarget == i)):
                        if ((triangleTarget == i) or not (buttonsPrevious[0])):
                            move = True
                            triangleTarget = i
                            triangle[i] = mouse
                            next = mouse
                            if (len(triangle) >= 3):
                                old = triangleCenter
                                ellipse = circumellipse(triangle, lines, None)
                                ellipseCenter = centroid(ellipse)
                                triangleCenter = centroid(triangle)
            color = triangleColor
            ps.Circle(int(next[0]-(radius*0.5)), int(next[1]-(radius*0.5)), int(radius*2.0), color=color)
            if (i == 0):
                previous = triangle[len(triangle)-1]
            else:
                previous = triangle[i-1]
            ps.Line(int(previous[0]), int(previous[1]), int(next[0]), int(next[1]), color=color)
        if not (triangleCenter == None):
            ps.Circle(int(triangleCenter[0]-(radius*0.5)), int(triangleCenter[1]-(radius*0.5)), int(radius*2.0), color=triangleColor)
        for i in range(len(ellipse)):
            next = ellipse[i]
            color = ellipseColor
            if (i == 0):
                color = (255, 0, 255, 255)
            elif (i == 1):
                color = (255, 255, 0, 255)
            ps.Circle(int(next[0]-(radius*0.5)), int(next[1]-(radius*0.5)), int(radius*2.0), color=color)
            if (i == 0):
                previous = ellipse[len(ellipse)-1]
            else:
                previous = ellipse[i-1]
            ps.Line(int(previous[0]), int(previous[1]), int(next[0]), int(next[1]), color=color)
        if not (ellipseCenter == None):
            ps.Circle(int(ellipseCenter[0]-(radius*0.5)), int(ellipseCenter[1]-(radius*0.5)), int(radius*2.0), color=ellipseColor)
        for i in range(len(polygon)):
            next = polygon[i]
            dist = distance(next, mouse)
            #print(str(i)+" "+str(dist))
            if ((polygonTarget == i) or (dist < radius*2.0)):
                if (buttonsNext[1]):
                    if (len(polygon) == 0):
                        polygon = []
                    elif (i == 0):
                        polygon = polygon[1:]
                    elif (i < len(polygon)-1):
                        polygon = polygon[0:i]+polygon[(i+1):]
                    else:
                        polygon = polygon[0:i]
                    if (len(polygon) >= 3):
                        old = polygonCenter
                        polygonCenter = centroid(polygon)
                    break
                elif (buttonsNext[2]):
                    if not ((polygonTarget > -1) and not (polygonTarget == i)):
                        if ((polygonTarget == i) or not (buttonsPrevious[2])):
                            move = True
                            polygonTarget = i
                            polygon[i] = mouse
                            next = mouse
                            if (len(polygon) >= 3):
                                old = polygonCenter
                                polygonCenter = centroid(polygon)
            color = polygonColor
            ps.Circle(int(next[0]-(radius*0.5)), int(next[1]-(radius*0.5)), int(radius*2.0), color=color)
            if (i == 0):
                previous = polygon[len(polygon)-1]
            else:
                previous = polygon[i-1]
            ps.Line(int(previous[0]), int(previous[1]), int(next[0]), int(next[1]), color=color)
        if not (polygonCenter == None):
            ps.Circle(int(polygonCenter[0]-(radius*0.5)), int(polygonCenter[1]-(radius*0.5)), int(radius*2.0), color=polygonColor)
        ps.UpdateScreen()
        if not (move):
            if (triangleTarget < 0):
                if (len(triangle) < 3):
                    if ((buttonsNext[0]) and not (buttonsPrevious[0])):
                        triangle.append(mouse)
                        triangleTarget = len(triangle)
                        if (len(triangle) >= 3):
                            polygonReal = None
                            if (len(polygon) >= 3):
                                polygonReal = polygon
                            ellipse = circumellipse(triangle, lines, sort(polygonReal))
                            ellipseCenter = centroid(ellipse)
                            triangleCenter = centroid(triangle)
            else:
                if not (buttonsNext[0]):
                    triangleTarget = -1
                    if (len(triangle) >= 3):
                        polygonReal = None
                        if (len(polygon) >= 3):
                            polygonReal = polygon
                        ellipse = circumellipse(triangle, lines, sort(polygonReal))
                        ellipseCenter = centroid(ellipse)
                        triangleCenter = centroid(triangle)
            if (polygonTarget < 0):
                if ((buttonsNext[2]) and not (buttonsPrevious[2])):
                    polygon.append(mouse)
                    polygonTarget = len(polygon)
                    if (len(polygon) >= 3):
                        old = polygonCenter
                        polygonCenter = centroid(polygon)
            else:
                if not (buttonsNext[2]):
                    polygonTarget = -1
                    if (len(triangle) >= 3):
                        polygonReal = None
                        if (len(polygon) >= 3):
                            polygonReal = polygon
                        ellipse = circumellipse(triangle, lines, sort(polygonReal))
                        ellipseCenter = centroid(ellipse)
                        triangleCenter = centroid(triangle)
        if (quit):
            break
    ps.UnloadScreen()
    return 0

def launch(arguments):
    result = run((1280, 720))
    print(str(result))
    if not (result == 0):
        return False
    return True

if (__name__ == "__main__"):
    print(str(launch(sys.argv)))