Код: Выделить всё
import numpy as np
import matplotlib.pyplot as plt
def conventional_form(a, b, c):
# Find longest side and scale to length 1
side1 = np.abs(a - b)
side2 = np.abs(b - c)
side3 = np.abs(c - a)
if side1 >= side2 and side1 >= side3:
longest_side = side1
v1, v2, v3 = a, b, c
elif side2 >= side1 and side2 >= side3:
longest_side = side2
v1, v2, v3 = b, c, a
else:
longest_side = side3
v1, v2, v3 = c, a, b
scaling_factor = 1 / longest_side
scaled_v1 = v1 * scaling_factor
scaled_v2 = v2 * scaling_factor
scaled_v3 = v3 * scaling_factor
# Rotating scaled_vertices so the longest side is horizontal and the other vertex is above the horizontal line
theta = np.angle(v2 - v1)
rotation_factor = np.exp(-1j * theta)
rotated_v1 = scaled_v1 * rotation_factor
rotated_v2 = scaled_v2 * rotation_factor
rotated_v3 = scaled_v3 * rotation_factor
if rotated_v3.imag < rotated_v1.imag:
rotated_v3 = complex(rotated_v3.real, 2 * rotated_v1.imag - rotated_v3.imag)
# Translate rotated_vertices so v1 = 0 and v2 = 0 + 1j
translation = -rotated_v1
a_new = rotated_v1 + translation
b_new = rotated_v2 + translation
c_new = rotated_v3 + translation
# Moving top vertex to left-side if needed
if c_new.real >= 0.50:
c_new = complex(1 - c_new.real, c_new.imag)
return a_new, b_new, c_new
def plot_point(point, color='black', markersize=2):
plt.scatter(point.real, point.imag, color=color, s=markersize)
def barycentric_subdivision(a, b, c, subdivisions):
if subdivisions == 0:
return
# Transform the vertices to the conventional form
a_new, b_new, c_new = conventional_form(a, b, c)
# Plot the c_new vertex
plot_point(c_new)
# Define the 6 new triangles
ab_mid = (a + b) / 2
bc_mid = (b + c) / 2
ca_mid = (c + a) / 2
centroid = (a + b + c) / 3
triangles = [
(a, ab_mid, centroid),
(ab_mid, b, centroid),
(b, bc_mid, centroid),
(bc_mid, c, centroid),
(c, ca_mid, centroid),
(ca_mid, a, centroid)
]
# Recursively apply barycentric subdivision to each new triangle
for tri in triangles:
barycentric_subdivision(*tri, subdivisions - 1)
def calculate_plotted_points(subdivisions):
return 1 + 6 ** (subdivisions - 1)
# Input
a = complex(1, 1.5)
b = complex(1.5, 2.5)
c = complex(0.5, 2.5)
subdivisions = 7
plt.figure()
barycentric_subdivision(a, b, c, subdivisions)
plt.xlim(0, 1)
plt.ylim(0, 1)
num_plotted_points = calculate_plotted_points(subdivisions)
plt.title(f"Subdivisions: {subdivisions}\nPlotted Points: {num_plotted_points}")
plt.show()
Сначала дан треугольник масштабируйте его так, чтобы его самая длинная сторона имела длину один. Затем поверните его так, чтобы его длинная сторона была горизонтальной. При необходимости поверните его еще раз (на 180$^\circ$), чтобы треугольник оказался на вершине самой длинной стороны. При необходимости переверните его вокруг вертикальной оси так, чтобы его самая короткая сторона оказалась слева. Треугольник, который мы сейчас рассматриваем, очень особенный: его нижняя сторона представляет собой горизонтальный сегмент длиной 1 (потому что мы соответствующим образом масштабировали его); его верхняя вершина лежит на расстоянии не более одной от его правой вершины (поскольку длина правой стороны не превышает длины нижней); а верхняя вершина лежит слева от центра нижней (из-за нашего соглашения относительно зеркальных изображений). Если нам задана система координат, мы можем даже разместить нижнюю часть отрезка от (0, 0) до (1,0).
При использовании 6 подразделений эта программа занимает около 20 секунд, но 7 так и не завершились. Спасибо всем, кто может помочь.
Подробнее здесь: https://stackoverflow.com/questions/790 ... ubdivision
