计算 Python 中 SVG 弧的准确长度?
Calculating the exact length of an SVG Arc in Python?
我希望能够计算出 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
,则将其作为选择加入代码,它将使用精确的超级花式解决方案,否则它将退回到较弱的近似代码(逐渐缩小线段,直到误差很小)。问题是这里的数学水平在我之上。而且我不知道是否有一些方法可以为此指定起点和终点。
大多数近似解都是针对椭圆的,但我只想要圆弧。可能还有一个我不知道的解决方案,用于计算椭圆上弧的长度,但因为起点和终点位置可以在任何地方。说后掠角是可能的总角度的 15% 似乎不是立即可行的,因此它是椭圆周长的 15%。
更有效的花式弧近似也可能不错。有越来越好的椭圆近似值,但我不能从椭圆周长到弧长,所以这些目前没有帮助。
假设圆弧参数化是椭圆上的起点和终点。因为这就是 SVG 的参数化方式。但是,任何不像 arc_length 参数化这样的重言式都是正确答案。
如果你想用你的双手和标准库来计算这个,你可以根据 following formula 进行计算。由于 acos
,这仅对椭圆上半部分的两个点有效,但我们将直接将其与角度一起使用。
计算包括以下步骤:
- 从 SVG 数据开始:起点、a、b 旋转、长弧、扫描、终点
- 旋转坐标系以匹配椭圆的水平轴。
- 求解具有 4 个未知数的 4 个方程组以获得中心点和对应于起点和终点的角度
- 通过对小段的离散求和来近似积分。这是您可以使用
scipy.special.ellipeinc
的地方,如评论中所建议的那样。
第二步很简单,用一个旋转矩阵就可以了(注意角度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 段,获取每个段末尾的函数值,按段长度计时,然后将所有这些相加。您可以尝试通过在每个段的 mid-point 处获取函数的值来改进它,但我不确定这样做有什么好处。段的数量可能对精度有更大的影响。
这是上面的一个非常粗略的 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。让我知道这是否值得追求,当我坐在电脑前时,我可以进一步查看。或者也许有人可以同时想出一个干净的方法来做到这一点?
我希望能够计算出 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
,则将其作为选择加入代码,它将使用精确的超级花式解决方案,否则它将退回到较弱的近似代码(逐渐缩小线段,直到误差很小)。问题是这里的数学水平在我之上。而且我不知道是否有一些方法可以为此指定起点和终点。
大多数近似解都是针对椭圆的,但我只想要圆弧。可能还有一个我不知道的解决方案,用于计算椭圆上弧的长度,但因为起点和终点位置可以在任何地方。说后掠角是可能的总角度的 15% 似乎不是立即可行的,因此它是椭圆周长的 15%。
更有效的花式弧近似也可能不错。有越来越好的椭圆近似值,但我不能从椭圆周长到弧长,所以这些目前没有帮助。
假设圆弧参数化是椭圆上的起点和终点。因为这就是 SVG 的参数化方式。但是,任何不像 arc_length 参数化这样的重言式都是正确答案。
如果你想用你的双手和标准库来计算这个,你可以根据 following formula 进行计算。由于 acos
,这仅对椭圆上半部分的两个点有效,但我们将直接将其与角度一起使用。
计算包括以下步骤:
- 从 SVG 数据开始:起点、a、b 旋转、长弧、扫描、终点
- 旋转坐标系以匹配椭圆的水平轴。
- 求解具有 4 个未知数的 4 个方程组以获得中心点和对应于起点和终点的角度
- 通过对小段的离散求和来近似积分。这是您可以使用
scipy.special.ellipeinc
的地方,如评论中所建议的那样。
第二步很简单,用一个旋转矩阵就可以了(注意角度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 段,获取每个段末尾的函数值,按段长度计时,然后将所有这些相加。您可以尝试通过在每个段的 mid-point 处获取函数的值来改进它,但我不确定这样做有什么好处。段的数量可能对精度有更大的影响。
这是上面的一个非常粗略的 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。让我知道这是否值得追求,当我坐在电脑前时,我可以进一步查看。或者也许有人可以同时想出一个干净的方法来做到这一点?