基于python计算曲线的曲率
文章目录一、实现原理1.1、计算点到直线的距离——海伦公式1.2、弓高和弦长计算半径二、python实现曲率计算最近需要对曲线的曲率做一个粗略的估计,在此记录下。其实计算曲率就是为了求这段弧长对应的半径,也就是说,我们把曲线看成圆的弧长就行,那么问题就简单了。一、实现原理1.1、计算点到直线的距离——海伦公式如下图所示,要计算A到CB的长度。设Δ\DeltaΔABC的三条边分别为a,b,c,那么海
最近需要对曲线的曲率做一个粗略的估计,在此记录下。其实计算曲率就是为了求这段弧长对应的半径,也就是说,我们把曲线看成圆的弧长就行,那么问题就简单了。
一、实现原理
1.1、计算点到直线的距离——海伦公式
如下图所示,要计算A到CB的长度。
设
Δ
\Delta
ΔABC的三条边分别为a,b,c,那么海伦公式计算面积S如下:
S
=
p
(
p
−
a
)
(
p
−
b
)
(
p
−
c
)
其
中
:
p
=
1
2
(
a
+
b
+
c
)
S=\sqrt{p(p-a)(p-b)(p-c)} \\ 其中: p=\frac{1}{2}(a+b+c)
S=p(p−a)(p−b)(p−c)其中:p=21(a+b+c)
将p带入海伦公式,得:
S
=
1
4
(
a
+
b
+
c
)
(
a
+
b
−
c
)
(
a
+
c
−
b
)
(
b
+
c
−
a
)
S=\frac{1}{4}\sqrt{(a+b+c)(a+b-c)(a+c-b)(b+c-a)}
S=41(a+b+c)(a+b−c)(a+c−b)(b+c−a)
如下图所示,AD是A到直线BC的距离,即是AD是
Δ
\Delta
ΔABC边BC上的高
S
=
a
∗
A
D
2
S=\frac{a*AD}{2}
S=2a∗AD
即:
1
4
(
a
+
b
+
c
)
(
a
+
b
−
c
)
(
a
+
c
−
b
)
(
b
+
c
−
a
)
=
a
∗
A
D
2
\frac{1}{4}\sqrt{(a+b+c)(a+b-c)(a+c-b)(b+c-a)}=\frac{a*AD}{2}
41(a+b+c)(a+b−c)(a+c−b)(b+c−a)=2a∗AD
合并可以得:
A
D
=
1
2
a
(
a
+
b
+
c
)
(
a
+
b
−
c
)
(
a
+
c
−
b
)
(
b
+
c
−
a
)
AD=\frac{1}{2a}\sqrt{(a+b+c)(a+b-c)(a+c-b)(b+c-a)}
AD=2a1(a+b+c)(a+b−c)(a+c−b)(b+c−a)
1.2、弓高和弦长计算半径
如图所示,我们已知了弦长L和弓高H,需要求R。我们可以根据三角形得勾股定理求得。
对于直角三角形OPQ,OP=R-H,OQ=R,PQ=
L
2
\frac{L}{2}
2L,那么计算如下:
R
2
=
(
L
2
)
2
+
(
R
−
H
)
2
R
2
=
L
4
4
+
R
2
−
2
R
H
+
H
2
2
R
H
=
L
4
4
+
H
2
R
=
L
2
4
+
H
2
2
H
R^{2}=(\frac{L}{2})^{2}+(R-H)^{2} \\ R^{2}= \frac{L^{4}}{4}+R^{2}-2RH+H^{2} \\ 2RH=\frac{L^{4}}{4}+H^{2} \\ R=\frac{\frac{L^{2}}{4}+H^{2}}{2H}
R2=(2L)2+(R−H)2R2=4L4+R2−2RH+H22RH=4L4+H2R=2H4L2+H2
二、python实现曲率计算
import numpy as np
def get_arc_curve(pts):
'''
获取弧度值
:param pts:
:return:
'''
# 计算弦长
start = np.array(pts[0])
end = np.array(pts[len(pts) - 1])
l_arc = np.sqrt(np.sum(np.power(end - start, 2)))
# 计算弧上的点到直线的最大距离
# 计算公式:\frac{1}{2a}\sqrt{(a+b+c)(a+b-c)(a+c-b)(b+c-a)}
a = l_arc
b = np.sqrt(np.sum(np.power(pts - start, 2), axis=1))
c = np.sqrt(np.sum(np.power(pts - end, 2), axis=1))
dist = np.sqrt((a + b + c) * (a + b - c) * (a + c - b) * (b + c - a)) / (2 * a)
h = dist.max()
# 计算曲率
r = ((a * a) / 4 + h * h) / (2 * h)
return r
if __name__ == '__main__':
x = np.linspace(1, 100, 99).astype(np.int64)
y = (x ** 2)
xy = list(zip(x, y)) # list of points in 2D space
print(get_arc_curve(xy))
参考连接:
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)