Если вы хотите рассчитать это своими руками и стандартной библиотекой, вы можете основывать свои вычисления на следующей формуле . Это верно только для двух точек в верхней половине эллипса из-за acos, но мы собираемся использовать его с углами напрямую.
Расчетсостоит из следующих шагов:
- Начало с данными SVG: начальная точка, вращение a, b, длинная дуга, развертка, конечная точка
- Поворот системы координат в соответствии с горизонтальной осьюэллипса.
- Решите систему из 4 уравнений с 4 неизвестными, чтобы получить центральную точку и углы, соответствующие начальной и конечной точке
- Приблизьте интеграл дискретной суммой по маленьким отрезкам,Здесь вы можете использовать
scipy.special.ellipeinc, как указано в комментариях.
Шаг 2 прост, просто используйте матрицу вращения (обратите внимание, что угол rot положительныйпо часовой стрелке):
m = [
[math.cos(rot), math.sin(rot)],
[-math.sin(rot), math.cos(rot)]
]
Шаг 3 очень хорошо объяснен в этом ответе . Обратите внимание, что значение, полученное для a1, равно модулю пи, потому что оно получено с atan. Это означает, что вам нужно вычислить центральные точки для двух углов t1 и t2 и проверить, совпадают ли они. Если они этого не делают, добавьте пи к a1 и проверьте снова.
Шаг 4 довольно прост. Разделите интервал [t1, t2] на n сегментов, получите значение функции в конце каждого сегмента, рассчитайте время по длине сегмента и суммируйте все это. Вы можете попытаться уточнить это, взяв значение функции в средней точке каждого сегмента, но я не уверен, что это принесет большую пользу. Количество сегментов, скорее всего, окажет большее влияние на точность.
Вот очень грубая версия Python выше (имейте в виду уродливый стиль кодирования, я делал это на своем мобильном телефоне, пока путешествовал ?)
import math
PREC = 1E-6
# matrix vector multiplication
def transform(m, p):
return ((sum(x * y for x, y in zip(m_r, p))) for m_r in m)
# the partial integral function
def ellipse_part_integral(t1, t2, a, b, n=100):
# function to integrate
def f(t):
return math.sqrt(1 - (1 - a**2 / b**2) * math.sin(t)**2)
start = min(t1, t2)
seg_len = abs(t1 - t2) / n
return - b * sum(f(start + seg_len * (i + 1)) * seg_len for i in range(n))
def ellipse_arc_length(x1, y1, a, b, rot, large_arc, sweep, x2, y2):
if abs(x1 - x2) < PREC and abs(y1 - y2) < PREC:
return 0
# get rot in radians
rot = math.pi / 180 * rot
# get the coordinates in the rotated coordinate system
m = [
[math.cos(rot), math.sin(rot)],
[- math.sin(rot), math.cos(rot)]
]
x1_loc, y1_loc, x2_loc, y2_loc = *transform(m, (x1,y1)), *transform(m, (x2,y2))
r1 = (x1_loc - x2_loc) / (2 * a)
r2 = (y2_loc - y1_loc) / (2 * b)
# avoid division by 0 if both points have same y coord
if abs(r2) > PREC:
a1 = math.atan(r1 / r2)
else:
a1 = r1 / abs(r1) * math.pi / 2
if abs(math.cos(a1)) > PREC:
a2 = math.asin(r2 / math.cos(a1))
else:
a2 = math.asin(r1 / math.sin(a1))
# calculate the angle of start and end point
t1 = a1 + a2
t2 = a1 - a2
# calculate centre point coords
x0 = x1_loc - a * math.cos(t1)
y0 = y1_loc - b * math.sin(t1)
x0s = x2_loc - a * math.cos(t2)
y0s = y2_loc - b * math.sin(t2)
# a1 value is mod pi so the centres may not match
# if they don't, check a1 + pi
if abs(x0 - x0s) > PREC or abs(y0 - y0s) > PREC:
a1 = a1 + math.pi
t1 = a1 + a2
t2 = a1 - a2
x0 = x1_loc - a * math.cos(t1)
y0 = y1_loc - b * math.sin(t1)
x0s = x2_loc - a * math.cos(t2)
y0s = y2_loc - b * math.sin(t2)
# get the angles in the range [0, 2 * pi]
if t1 < 0:
t1 += 2 * math.pi
if t2 < 0:
t2 += 2 * math.pi
# increase minimum by 2 * pi for a large arc
if large_arc:
if t1 < t2:
t1 += 2 * math.pi
else:
t2 += 2 * math.pi
return ellipse_part_integral(t1, t2, a, b)
print(ellipse_arc_length(0, 0, 40, 40, 0, False, True, 80, 0))
Хорошей новостью является то, что флаг развертки не имеет значения, пока вы просто ищете длину дуги.
Я не уверен на 100%Проблема по модулю pi обрабатывается правильно, и реализация выше может иметь несколько ошибок. Тем не менее, это дало мне хорошее приближение длины в простом случае полукруга, поэтому я осмелюсь назвать это WIP. Дайте мне знать, если это стоит того, чтобы я посмотрел, когда я буду сидеть за компьютером. Или, может быть, кто-то может придумать чистый способ сделать это тем временем?