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
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) / 2
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
I3 = numpy.linalg.det([[a, b, d / 2], [b, 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 + 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]]):
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], [b, a]])
w, ww = np.linalg.eig(matrix)
ww = np.squeeze(np.asarray(ww))
y0 = (2 * d * b - 2 * a * e) / (4 * a * c - 4 * b * b) # y-position of the center
x0 = -(e + 2 * c * y0) / (2 * 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 + 2 * b * x0 * y0 + c * y0 * y0 - f
delta_x = r_2 / (a + 2 * b * k + c * k * k)
delta_y = k * k * delta_x
ax = sqrt(delta_x + delta_y)
delta_x = r_2 / (a + 2 * 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], '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:
plt.plot([ellipses[i][2], ellipses[j][2]], [ellipses[i][1], ellipses[j][1]], 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])
return result
def empty_new_ellipse(p1, p2, p3, p4, b, c, p5):
const1 = p1[0] ** 2 + 2 * b * p1[0] * p1[1] + c * p1[1] ** 2
const2 = p2[0] ** 2 + 2 * b * p2[0] * p2[1] + c * p2[1] ** 2
const3 = p3[0] ** 2 + 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, p1[1], 1], [-const2, p2[1], 1], [-const3, p3[1], 1]])
det2 = numpy.linalg.det([[p1[0], -const1, 1], [p2[0], -const2, 1], [p3[0], -const3, 1]])
det3 = numpy.linalg.det([[p1[0], p1[1], -const1], [p2[0], p2[1], -const2], [p3[0], p3[1], -const3]])
if det != 0:
d = det1 / det
e = det2 / det
f = det3 / det
# plt.show()
# 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()
# print(p1, p2, p3, p4)
# x = sympy.symbols('x')
# y = sympy.symbols('y')
# sympy.plot_implicit(sympy.Eq(x*x + 2 * b * x * y + c * y*y + d * x + e * y + f))
#
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
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]
if 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(i[0], i[1], 'bo')
plt.plot(empty_ellipses[0][2], empty_ellipses[0][1], 'go')
plt.show()
# plt.savefig('/home/kate/Рабочий стол/Figure_1.png')