import numpy.linalg
import numpy.random
from matplotlib import pyplot as plt
import numpy as np
import matplotlib.path as path
from math import pi, cos, sin, atan, asin, fabs, sqrt
import sympy
def ellipse_by_5_points(p):
A = [[p[0][0] * p[0][1], p[0][1] * p[0][1], p[0][0], p[0][1], 1],
[p[1][0] * p[1][1], p[1][1] * p[1][1], p[1][0], p[1][1], 1],
[p[2][0] * p[2][1], p[2][1] * p[2][1], p[2][0], p[2][1], 1],
[p[3][0] * p[3][1], p[3][1] * p[3][1], p[3][0], p[3][1], 1],
[p[4][0] * p[4][1], p[4][1] * p[4][1], p[4][0], p[4][1], 1]]
B = [[p[0][0] * p[0][0], p[0][1] * p[0][1], p[0][0], p[0][1], 1],
[p[1][0] * p[1][0], p[1][1] * p[1][1], p[1][0], p[1][1], 1],
[p[2][0] * p[2][0], p[2][1] * p[2][1], p[2][0], p[2][1], 1],
[p[3][0] * p[3][0], p[3][1] * p[3][1], p[3][0], p[3][1], 1],
[p[4][0] * p[4][0], p[4][1] * p[4][1], p[4][0], p[4][1], 1]]
C = [[p[0][0] * p[0][0], p[0][0] * p[0][1], p[0][0], p[0][1], 1],
[p[1][0] * p[1][0], p[1][0] * p[1][1], p[1][0], p[1][1], 1],
[p[2][0] * p[2][0], p[2][0] * p[2][1], p[2][0], p[2][1], 1],
[p[3][0] * p[3][0], p[3][0] * p[3][1], p[3][0], p[3][1], 1],
[p[4][0] * p[4][0], p[4][0] * p[4][1], p[4][0], p[4][1], 1]]
D = [[p[0][0] * p[0][0], p[0][0] * p[0][1], p[0][1] * p[0][1], p[0][1], 1],
[p[1][0] * p[1][0], p[1][0] * p[1][1], p[1][1] * p[1][1], p[1][1], 1],
[p[2][0] * p[2][0], p[2][0] * p[2][1], p[2][1] * p[2][1], p[2][1], 1],
[p[3][0] * p[3][0], p[3][0] * p[3][1], p[3][1] * p[3][1], p[3][1], 1],
[p[4][0] * p[4][0], p[4][0] * p[4][1], p[4][1] * p[4][1], p[4][1], 1]]
E = [[p[0][0] * p[0][0], p[0][0] * p[0][1], p[0][1] * p[0][1], p[0][0], 1],
[p[1][0] * p[1][0], p[1][0] * p[1][1], p[1][1] * p[1][1], p[1][0], 1],
[p[2][0] * p[2][0], p[2][0] * p[2][1], p[2][1] * p[2][1], p[2][0], 1],
[p[3][0] * p[3][0], p[3][0] * p[3][1], p[3][1] * p[3][1], p[3][0], 1],
[p[4][0] * p[4][0], p[4][0] * p[4][1], p[4][1] * p[4][1], p[4][0], 1]]
F = [[p[0][0] * p[0][0], p[0][0] * p[0][1], p[0][1] * p[0][1], p[0][0], p[0][1]],
[p[1][0] * p[1][0], p[1][0] * p[1][1], p[1][1] * p[1][1], p[1][0], p[1][1]],
[p[2][0] * p[2][0], p[2][0] * p[2][1], p[2][1] * p[2][1], p[2][0], p[2][1]],
[p[3][0] * p[3][0], p[3][0] * p[3][1], p[3][1] * p[3][1], p[3][0], p[3][1]],
[p[4][0] * p[4][0], p[4][0] * p[4][1], p[4][1] * p[4][1], p[4][0], p[4][1]]]
a = numpy.linalg.det(A)
b = -numpy.linalg.det(B)
c = numpy.linalg.det(C)
d = -numpy.linalg.det(D)
e = numpy.linalg.det(E)
f = -numpy.linalg.det(F)
b = b / a
c = c / a
d = d / a
e = e / a
f = f / a
a = 1
I1 = a + c
I2 = a * c - (b / 2) ** 2
I3 = numpy.linalg.det([[a, b / 2, d / 2], [b / 2, c, e / 2], [d / 2, e / 2, f]])
if I2 > 0 and I1 * I3 < 0:
return [a, b, c, d, e, f]
else:
return None
def in_ellipse(p, k):
if k[0] * p[0] ** 2 + k[1] * p[0] * p[1] + k[2] * p[1] ** 2 + k[3] * p[0] + k[4] * p[1] + k[5] < 0:
return True
return False
def empty_ellipse(k, points, used_points):
for i in points:
if i not in used_points and in_ellipse(i, k):
return False
return True
def empty_ellipses_into_5_points(points):
empty_ellipses = []
e_numbers = []
for p1 in range(len(p)):
for p2 in range(p1 + 1, len(p)):
for p3 in range(p2 + 1, len(p)):
for p4 in range(p3 + 1, len(p)):
for p5 in range(p4 + 1, len(p)):
if points_for_ellipse([p[p1], p[p2], p[p3], p[p4], p[p5]]):
k = ellipse_by_5_points([points[p1], points[p2], points[p3], points[p4], points[p5]])
if k is not None:
if empty_ellipse(k, points, [p[p1], p[p2], p[p3], p[p4], p[p5]]):
# draw_ellipse(k)
print(p1, p2, p3, p4, p5)
empty_ellipses.append(k)
e_numbers.append([p1, p2, p3, p4, p5])
return empty_ellipses, e_numbers
def draw_ellipse(ellipse):
a = ellipse[0]
b = ellipse[1]
c = ellipse[2]
d = ellipse[3]
e = ellipse[4]
f = ellipse[5]
matrix = np.mat([[c, b / 2], [b / 2, a]])
w, ww = np.linalg.eig(matrix)
ww = np.squeeze(np.asarray(ww))
y0 = (d * b - 2 * a * e) / (4 * a * c - b * b) # y-position of the center
x0 = -(e + 2 * c * y0) / b # x-position of the center
k = ww[0][1] / ww[0][0]
k1 = ww[1][1] / ww[1][0]
r_2 = a * x0 * x0 + b * x0 * y0 + c * y0 * y0 - f
delta_x = r_2 / (a + b * k + c * k * k)
delta_y = k * k * delta_x
ax = sqrt(delta_x + delta_y)
delta_x = r_2 / (a + b * k1 + c * k1 * k1)
delta_y = k1 * k1 * delta_x
by = sqrt(delta_x + delta_y)
t_rot = atan(k) # rotation angle
t = np.linspace(0, 2 * pi, 100)
Ell = np.array([ax * np.cos(t), by * np.sin(t)])
R_rot = np.array([[cos(t_rot), -sin(t_rot)], [sin(t_rot), cos(t_rot)]])
Ell_rot = np.zeros((2, Ell.shape[1]))
for i in range(Ell.shape[1]):
Ell_rot[:, i] = np.dot(R_rot, Ell[:, i])
plt.plot(x0 + Ell_rot[0, :], y0 + Ell_rot[1, :], 'darkorange') # rotated ellipse
def create_points(N):
points = []
while len(points) < N:
new_p = [numpy.random.randint(-10, 10), numpy.random.randint(-10, 10)]
if new_p not in points:
points.append(new_p)
return points
def rotate(A, B, C):
return (B[0] - A[0]) * (C[1] - B[1]) - (B[1] - A[1]) * (C[0] - B[0])
def jarvismarch(A):
n = len(A)
P = [i for i in range(n)]
for i in range(1, n):
if A[P[i]][0] < A[P[0]][0]:
P[i], P[0] = P[0], P[i]
H = [P[0]]
del P[0]
P.append(H[0])
while True:
right = 0
for i in range(1, len(P)):
if rotate(A[H[-1]], A[P[right]], A[P[i]]) < 0:
right = i
if P[right] == H[0]:
break
else:
H.append(P[right])
del P[right]
H = [A[H[i]] for i in range(len(H))]
for i in range(len(H)):
if angle_180(H[i - 2], H[i - 1], H[i]):
return None
return H
def angle_180(a, b, c):
A = [a[0] - b[0], a[1] - b[1]]
B = [c[0] - b[0], c[1] - b[1]]
skal = A[0] * B[0] + A[1] * B[1]
d_A = A[0] * A[0] + A[1] * A[1]
d_B = B[0] * B[0] + B[1] * B[1]
if skal / (d_A * d_B) == -1:
return True
else:
return False
def points_for_ellipse(p):
conv = jarvismarch(p)
if (conv is not None and len(conv) == 5):
for i in range(len(p)):
return True
else:
return False
def print_omega(ellipses, numbers):
y = np.linspace(-5, 5, 1000)
x = y ** 2
plt.grid()
plt.xlabel('c')
plt.ylabel('b')
plt.plot(x, y)
for i in ellipses:
plt.plot(i[2], i[1] / 2, 'ro')
for i in range(len(numbers)):
for j in range(i + 1, len(numbers)):
c = list(set(numbers[i]) & set(numbers[j]))
if len(c) == 4:
print(c, numbers[i], numbers[j])
plt.plot([ellipses[i][2], ellipses[j][2]], [ellipses[i][1] / 2, ellipses[j][1] / 2], color='orange')
def koordinates_by_bases(A, B, C):
det = B[0] * C[1] - C[0] * B[1]
if det != 0:
det_1 = A[0] * C[1] - C[0] * A[1]
det_2 = B[0] * A[1] - A[0] * B[1]
return [det_1 / det, det_2 / det]
def solve_equation(_A, _B, _C):
D = _B ** 2 - 4 * _A * _C
if D == 0:
return [-_B / (2 * _A)]
if D < 0:
return None
if D > 0:
return [(-_B - sqrt(D)) / (2 * _A), (-_B + sqrt(D)) / (2 * _A)]
def line_for_4_poitns(p1, p2, p3, p4):
P21 = [p2[0] - p1[0], p2[1] - p1[1]]
P31 = [p3[0] - p1[0], p3[1] - p1[1]]
P41 = [p4[0] - p1[0], p4[1] - p1[1]]
a = koordinates_by_bases(P41, P21, P31)
k1 = a[0] * P21[0] ** 2 + a[1] * P31[0] ** 2 - P41[0] ** 2
k2 = 2 * (a[0] * P21[0] * P21[1] + a[1] * P31[0] * P31[1] - P41[0] * P41[1])
k3 = a[0] * P21[1] ** 2 + a[1] * P31[1] ** 2 - P41[1] ** 2
result = []
b = solve_equation(k3, k2, k1)
for i in b:
result.append([i ** 2, i/2])
return result
def empty_new_ellipse(p1, p2, p3, p4, b, c, p5):
const1 = p1[0] ** 2 + b * p1[0] * p1[1] + c * p1[1] ** 2
const2 = p2[0] ** 2 + b * p2[0] * p2[1] + c * p2[1] ** 2
const3 = p3[0] ** 2 + b * p3[0] * p3[1] + c * p3[1] ** 2
det = numpy.linalg.det([[p1[0], p1[1], 1], [p2[0], p2[1], 1], [p3[0], p3[1], 1]])
det1 = numpy.linalg.det([[-const1, -const2, -const3], [p2[0], p2[1], 1], [p3[0], p3[1], 1]])
det2 = numpy.linalg.det([[p1[0], p1[1], 1], [-const1, -const2, -const3], [p3[0], p3[1], 1]])
det3 = numpy.linalg.det([[p1[0], p1[1], 1], [p2[0], p2[1], 1], [-const1, -const2, -const3]])
if det != 0:
d = det1 / det
e = det2 / det
f = det3 / det
plt.plot(p1[0], p1[1], 'go')
plt.plot(p2[0], p2[1], 'go')
plt.plot(p3[0], p3[1], 'go')
plt.plot(p4[0], p4[1], 'go')
plt.plot(p5[0], p5[1], 'ro')
draw_ellipse([1, b, c, d, e, f])
plt.show()
return in_ellipse(p5, [1, b, c, d, e, f])
# p = create_points(8)
# print(p)
# p = [[-3, -6], [6, 8], [-8, 8], [-4, 8], [-8, 7], [0, 0], [0, 9], [2, -6]] divide by zero
p = [[6, -2], [-10, 8], [-1, -8], [-4, -4], [4, -2], [-7, -9], [-4, 1], [8, -6]]
empty_ellipses, e_numbers = empty_ellipses_into_5_points(p)
# k = 0
# for i in p:
# plt.plot(i[0], i[1], 'bo')
# plt.annotate(str("%d" % k), [i[0], i[1]])
# k = k + 1
# plt.show()
# plt.grid()
#
# plt.gca().set_aspect("equal")
# plt.show()
# print_omega(empty_ellipses, e_numbers)
points_on_parab = line_for_4_poitns(p[e_numbers[0][0]], p[e_numbers[0][1]], p[e_numbers[0][2]], p[e_numbers[0][3]])
for i in points_on_parab:
new_point = [(i[0] + empty_ellipses[0][2]) / 2, (i[1] + empty_ellipses[0][1])/2]
empty_new_ellipse(p[e_numbers[0][0]], p[e_numbers[0][1]], p[e_numbers[0][2]], p[e_numbers[0][3]], new_point[1], new_point[0], p[e_numbers[0][4]])
# plt.plot(new_point[0], new_point[1], 'bo')
plt.show()
# plt.savefig('/home/kate/Рабочий стол/Figure_1.png')