import numpy.linalg
import numpy.random
from matplotlib import pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull
import sympy
import matplotlib.path as path
from math import pi, cos, sin, atan, asin, fabs, sqrt
def calculate_b_c(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 = (-1) * numpy.linalg.det(B)
c = numpy.linalg.det(C)
d = (-1) * numpy.linalg.det(D)
e = numpy.linalg.det(E)
f = (-1) * numpy.linalg.det(F)
return a, b, c, d, e, f
def in_ellipse(p, koef):
if koef[0] * p[0] ** 2 + koef[1] * p[0] * p[1] + koef[2] * p[1] ** 2 + koef[3] * p[0] + koef[4] * p[1] + koef[
5] < 0:
return True
return False
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))]
return H
def empty_ellipse(p, points):
bbPath = path.Path(p)
a, b, c, d, e, f = calculate_b_c(p)
for i in points:
if [i[0], i[1]] not in p:
if bbPath.contains_point(i, radius=0.0):
return False
else:
if in_ellipse(i, [a, b, c, d, e, f]):
return False
return [a, b, c, d, e, f]
def equation(a, b, c):
discriminant = b ** 2 - 4 * a * c
if discriminant < 0:
print('Корней нет')
return None
elif discriminant == 0:
x = -b / (2 * a)
return x
else:
x1 = (-b + discriminant ** 0.5) / (2 * a)
x2 = (-b - discriminant ** 0.5) / (2 * a)
return max(abs(x1), abs(x2))
points = []
while len(points) < 7:
new_p = [numpy.random.randint(-10, 10), numpy.random.randint(-10, 10)]
if new_p not in points:
points.append(new_p)
print(points)
# points = [[-5, -3], [-4, 4], [2, 9], [-1, 8], [1, -4], [3, 4], [7, 4]]
empty_ellipses = []
ee = []
for i1 in range(len(points)):
for i2 in range(i1 + 1, len(points)):
for i3 in range(i2 + 1, len(points)):
for i4 in range(i3 + 1, len(points)):
for i5 in range(i4 + 1, len(points)):
conv = jarvismarch([points[i1], points[i2], points[i3], points[i4], points[i5]])
if (len(conv) == 5):
flag = empty_ellipse(conv, points)
if flag:
# print(i1, i2, i3, i4, i5)
empty_ellipses.append(flag)
ee.append([points[i1], points[i2], points[i3], points[i4], points[i5]])
k = 0
for i in points:
if i in ee[0]:
plt.plot(i[0], i[1], 'ro')
else:
plt.plot(i[0], i[1], 'bo')
plt.annotate(str("%d" % k), [i[0], i[1]])
k = k + 1
points = np.array(points)
plt.grid()
plt.axis([-11, 11, -11, 11], 'equal')
plt.gca().set_aspect("equal")
# x = sympy.symbols('x')
# y = sympy.symbols('y')
# p0 = sympy.plot_implicit(sympy.Eq(empty_ellipses[0][0] * x*x + empty_ellipses[0][1] * x * y + empty_ellipses[0][2] * y*y + empty_ellipses[0][3] * x + empty_ellipses[0][4] * y + empty_ellipses[0][5]), show=False)
# for i in empty_ellipses:
# p0.extend(sympy.plot_implicit(sympy.Eq(i[0] * (x ** 2) + i[1] * x * y + i[2] * (y ** 2) + i[3] * x + i[4] * y + i[5]), show=False))
# p0.show()
a = empty_ellipses[0][0]
b = empty_ellipses[0][1]
c = empty_ellipses[0][2]
d = empty_ellipses[0][3]
e = empty_ellipses[0][4]
f = empty_ellipses[0][5]
matrix = np.mat([[c, b/2], [b/2, a]])
print (matrix)
w, ww = np.linalg.eig(matrix)
print(ww)
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]
xrange = [-10, 0, 10]
yrange = [y0 + k * (_x - x0) for _x in xrange]
plt.plot(xrange, yrange)
yrange = [y0 + k1 * (_x - x0) for _x in xrange]
plt.plot(xrange, yrange)
# plt.plot([ww[0][0]+x0, x0, ww[1][0]+x0], [ww[0][1]+y0, y0, ww[1][1]+y0])
#
# _A = a + b * k + c * k * k
# _B = b * y0 - b * k * x0 - 2 * c * k * x0 + 2 * k * c * y0 + d + e * k
# _C = c * y0 * y0 + c * k * k * x0 * x0 - 2 * k * c * y0 * x0 + e * y0 - e * k * x0 + f
#
# print(_A, _B, _C)
#
# point_a = equation(_A, _B, _C)
# point_b = point_a * k + y0
#
# plt.plot(point_a, point_b, 'yo')
# ax = sqrt((point_a - x0) ** 2 + (point_b - y0) ** 2)
#
# _A = a + b * k1 + c * k1 * k1
# _B = b * y0 - b * k1 * x0 - 2 * c * k1 * x0 + 2 * k1 * c * y0 + d + e * k1
# _C = c * y0 * y0 + c * k1 * k1 * x0 * x0 - 2 * k1 * c * y0 * x0 + e * y0 - e * k1 * x0 + f
#
# print(_A, _B, _C)
#
# point_a = equation(_A, _B, _C)
# point_b = point_a * k1 + y0
#
# plt.plot(point_a, point_b, 'yo')
#
# by = sqrt((point_a - x0) ** 2 + (point_b - y0) ** 2)
#
# # radius on the x-axis
# # radius on the y-axis
#
# 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)]])
# # 2-D rotation matrix
#
# 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
# plt.grid(color='lightgray', linestyle='--')
# # plt.show()
plt.savefig('/home/kate/Рабочий стол/Figure_1.png')
# x = sympy.symbols('x')
# y = sympy.symbols('y')
# sympy.plot_implicit(sympy.Eq(a * x * x + b * x * y + c * y * y + d * x + e * y + f))