求闭合二维均匀三次 B 样条的面积

Question

我有一个二维点列表，它们是闭合均匀三次 B 样条的控制顶点 (Dx)。我假设一条简单的曲线（非自相交，所有控制点都是不同的）。

我试图找到曲线所包围的区域：

如果我计算结点（Px），我可以将曲线视为多边形；然后我“只”需要在实际曲线和连接结点的直线之间找到每个段的剩余增量区域。

我知道 Bspline 的形状（以及因此面积）在旋转和平移下是不变的 - 因此对于每个段，我可以找到一个平移来将 t=0 结放在原点，并找到一个旋转来放置 t =+x 轴上的 1 个结：

我可以通过插入点并重新分组来找到曲线的方程：

P(t) = (
    (t**3)*(-Dm1 + 3*D0 - 3*D1 + D2)
    + (t**2)*(3*Dm1 - 6*D0 + 3*D1)
    + t*(-3*Dm1 + 3*D1)
    + (Dm1 + 4*D0 + D1)
) / 6.

但我正在努力整合它 - 我可以做到

 1
/
| Py(t) dt
/
t=0

但这并没有给我空间。我认为我需要的是

 Px(t=1)
/
| Py(t) (dPx(t) / dt) dt
/
x = Px(t=0)

但在我走得更远之前，我真的很想知道：

1. 这是对面积的正确计算吗？理想情况下，分析解决方案会让我很开心！

2. 找到这个区域后，如何判断是否需要从基础多边形中添加或减去它（第一张图中的红色与绿色区域）？

3. 是否有任何 Python 模块可以为我进行此计算？ Numpy 有一些评估三次 Bsplines 的方法，但似乎没有一种方法可以处理面积。

4. 有没有更简单的方法来做到这一点？我正在考虑可能在一堆点上评估 P(t) - 比如t = numpy.arange(0.0, 1.0, 0.05) - 并将整个事物视为一个多边形。知道需要多少细分才能保证给定的准确度（我希望误差

Answer 1

就我个人而言，我会充分利用样条曲线，并使用Green's theorem 将面积积分重写为轮廓积分。

由于您已经知道曲线，因此使用足够阶的高斯求积进行积分将是一件容易的事。没有恶作剧来估计所需的额外区域。我敢打赌它的计算效率也会很高，因为多项式上的高斯求积表现得非常好。三次 B 样条将很好地集成。

我会以这样一种方式编写您的代码，即第一个点和最后一个点必须重合。这是问题的不变量。

这种方法甚至适用于有洞的区域。你沿着外曲线积分，制作一条连接外曲线和内曲线的假想直线，沿着内曲线积分，然后沿着到达那里的直线返回到外曲线。

Answer 2

1. 选择任意点作为枢轴p₀（例如原点 (0,0)）
2. 沿曲线选取一点p₁ = (x,y)
3. 对该点的曲线求微分，得到速度v = vx,vy>
4. 从三个点和calculate the area 组成一个三角形。使用向量 p₀p₁ 和 v 之间的叉积最容易完成，然后除以二。
5. 在 t 上积分此区域，从 0 到 1。

你得到的结果是整条曲线的一段的面积。有些人会是消极的，因为他们面对相反的方向。如果将所有线段的面积相加，则得到整条曲线的面积。

结果是：

面积_i = (31 x_i-1 y_i + 28 x_i-1 y_i+1 + x_{i- 1}y_i+2 - 31 x_iy_i-1 + 183 x_iy_i+1 + 28 x_iy_i+2 - 28 x_i+1y_i-1 - 183 x_i+1y_i + 31 x_i+1y_i+2 - x_{i +2}y_i-1 - 28 x_i+2y_i - 31 x_i+2y_i+1) / 720

如果你把它转换成矩阵形式，你会得到：

面积_i = x_{i-1   xi   xi+1   xi+2> · P · yi-1   yi   y i+1   yi+2>^T}

P在哪里

[    0   31   28    1]
[  -31    0  183   28] / 720
[  -28 -183    0   31]
[   -1  -28  -31    0]

---

如果控制点是[(0,0) (1,0) (1,1) (0,1)]，则结果区域变为：

[(0,0), (1,0), (1,1), (0,1)] -> 242/720
[(1,0), (1,1), (0,1), (0,0)] -> 242/720
[(1,1), (0,1), (0,0), (1,0)] ->   2/720
[(0,1), (0,0), (1,0), (1,1)] ->   2/720

总和是 488/720 = 61/90。

Answer 3

嗯……看来你已经知道该怎么做了。

1. 是的，这确实是计算线段下面积的正确方法

   dS=Y*dX=Y*(dX/dt)*dt
   

2. 您无需关心加法或减法。您需要一直添加，但是某些积分（红色段）将为负数（如果您始终设置将 P_n 放在坐标的开头，并且 P_{n+1 沿 X 轴正方向）。因此，对于每个段，您都需要计算平移和旋转，然后进行积分（全部通过解析完成）。}

3. 我不知道 Python 模块，但是制作这样一个模块似乎最多需要一天时间。
4. 我认为解析解会更好，也不会那么难。

无论如何，图形绝对令人惊叹