python - 在Python中计算SVG弧的精确长度？

Question

我希望能够计算 SVG 弧的精确长度。我可以相当轻松地完成所有操作。但是，我不确定是否有解决方案或解决方案的具体实现。

这是椭圆周长的精确解。使用流行的库就可以了。我完全明白没有简单的解决方案，因为它们都需要超几何函数才能精确。

from scipy import pi, sqrt
from scipy.special import hyp2f1

def exact(a, b):
    t = ((a - b) / (a + b)) ** 2
    return pi * (a + b) * hyp2f1(-0.5, -0.5, 1, t)

a = 2.667950e9
b = 6.782819e8
print(exact(a, b))

我的想法是，如果您碰巧安装了 scipy，则将此作为选择加入代码，它将使用精确的 super 解决方案，否则它将回退到较弱的近似代码(逐渐减小线段，直到误差很小)。问题是这里的数学水平高于我。我不知道是否有一些方法可以指定起点和终点。

大多数近似解都是针对椭圆的，但我只想要圆弧。可能还有一个我不知道的解决方案，用于计算椭圆上的弧长，但由于开始和结束位置可以在任何地方。说后掠角是总可能角度的 15% 因此它是椭圆周长的 15% 似乎并不立即可行。

更有效、更少花哨的弧近似也可能会更好。有逐渐更好的椭圆近似值，但我无法从椭圆周长到弧长，因此这些目前没有帮助。

<小时/>

假设圆弧参数化是椭圆上的起点和终点。因为这就是 SVG 参数化的方式。但是，任何与 arc_length 参数化不同义的内容都是正确的答案。

Answer 1

如果您想徒手使用标准库来计算此值，您可以根据 following formula 进行计算。 。由于 acos，这仅对椭圆上半部分的两个点有效，但我们将直接将其与角度一起使用。

计算包括以下步骤:

1. 从 SVG 数据开始:起点、a、b 旋转、长弧、扫描、终点
2. 旋转坐标系以匹配椭圆的水平轴。
3. 求解包含 4 个未知数的 4 个方程组，得到中心点以及起点和终点对应的角度
4. 通过小段上的离散总和来近似积分。您可以在此处使用 scipy.special.ellipeinc ，正如评论中所建议的。

第 2 步很简单，只需使用旋转矩阵即可(注意角度 rot 顺时针方向为正):

 m = [
        [math.cos(rot), math.sin(rot)], 
        [-math.sin(rot), math.cos(rot)]
   ]

第 3 步在 this answer 中有很好的解释。 。请注意，a1 获得的值是 pi 模，因为它是通过 atan 获得的。这意味着您需要计算两个角度 t1 和 t2 的中心点并检查它们是否匹配。如果没有，请将 pi 添加到 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。让我知道这是否值得追求，当我坐在电脑前时我可以进一步看看。或者也许有人可以同时想出一种干净的方法来做到这一点？