检查两条三次 Bézier 曲线是否相交
对于个人项目，我需要找出两条三次贝塞尔曲线是否相交。我不需要知道在哪里：我只需要知道他们是否知道。不过，我需要尽快完成。
我一直在这个地方搜寻并找到了一些资源。大多数情况下，这个问题有一个有希望的答案。
因此，在我弄清楚什么是西尔维斯特矩阵之后，什么是行列式，什么是 结果和为什么它有用，我想我知道该解决方案是如何工作的。然而，现实却与此不同，而且效果并不那么好。
闲逛
我使用图形计算器绘制了两条相交的贝塞尔样条线（我们将其称为 B0 和 B1）。它们的坐标如下（P0, P1, P2, P3）：
(1, 1) (2, 4) (3, 4) (4, 3)
(3, 5) (3, 6) (0, 1) (3, 1)
结果是以下，B0 是“水平”曲线，B1 是另一条：
按照上述问题的得票最高的答案的指示，我将 B0 减去 B1。它给我留下了两个方程（X 轴和 Y 轴），根据我的计算器，它们是：
x = 9t^3 - 9t^2 - 3t + 2
y = 9t^3 - 9t^2 - 6t + 4
西尔维斯特矩阵
并据此我构建了以下西尔维斯特矩阵：
9 -9 -3 2 0 0
0 9 -9 -3 2 0
0 0 9 -9 -3 2
9 -9 -6 4 0 0
0 9 -9 -6 4 0
0 0 9 -9 -6 4
之后，我创建了一个 C++ 函数来计算使用拉普拉斯展开的矩阵行列式：
template<int size>
float determinant(float* matrix)
{
float total = 0;
float sign = 1;
float temporaryMatrix[(size - 1) * (size - 1)];
for (int i = 0; i < size; i++)
{
if (matrix[i] != 0)
{
for (int j = 1; j < size; j++)
{
float* targetOffset = temporaryMatrix + (j - 1) * (size - 1);
float* sourceOffset = matrix + j * size;
int firstCopySize = i * sizeof *matrix;
int secondCopySize = (size - i - 1) * sizeof *matrix;
memcpy(targetOffset, sourceOffset, firstCopySize);
memcpy(targetOffset + i, sourceOffset + i + 1, secondCopySize);
}
float subdeterminant = determinant<size - 1>(temporaryMatrix);
total += matrix[i] * subdeterminant * sign;
}
sign *= -1;
}
return total;
}
template<>
float determinant<1>(float* matrix)
{
return matrix[0];
}
它似乎在相对较小的矩阵（2x2、3x3）上工作得很好和 4x4），所以我希望它也适用于 6x6 矩阵。不过，我没有进行广泛的测试，它有可能已损坏。
问题
如果我正确理解了另一个问题的答案，行列式应该为 0，因为曲线相交。然而，将我上面制作的西尔维斯特矩阵输入到我的程序中，它是-2916。
这是我的错误还是他们的错误？确定两条三次贝塞尔曲线是否相交的正确方法是什么？
For a personal project, I'd need to find out if two cubic Bézier curves intersect. I don't need to know where: I just need to know if they do. However, I'd need to do it fast.
I've been scavenging the place and I found several resources. Mostly, there's this question here that had a promising answer.
So after I figured what is a Sylvester matrix, what is a determinant, what is a resultant and why it's useful, I thought I figured how the solution works. However, reality begs to differ, and it doesn't work so well.
Messing Around
I've used my graphing calculator to draw two Bézier splines (that we'll call B0 and B1) that intersect. Their coordinates are as follow (P0, P1, P2, P3):
(1, 1) (2, 4) (3, 4) (4, 3)
(3, 5) (3, 6) (0, 1) (3, 1)
The result is the following, B0 being the "horizontal" curve and B1 the other one:
Following directions from the aforementioned question's top-voted answer, I've subtracted B0 to B1. It left me with two equations (the X and the Y axes) that, according to my calculator, are:
x = 9t^3 - 9t^2 - 3t + 2
y = 9t^3 - 9t^2 - 6t + 4
The Sylvester Matrix
And from that I've built the following Sylvester matrix:
9 -9 -3 2 0 0
0 9 -9 -3 2 0
0 0 9 -9 -3 2
9 -9 -6 4 0 0
0 9 -9 -6 4 0
0 0 9 -9 -6 4
After that, I've made a C++ function to calculate determinants of matrices using Laplace expansion:
template<int size>
float determinant(float* matrix)
{
float total = 0;
float sign = 1;
float temporaryMatrix[(size - 1) * (size - 1)];
for (int i = 0; i < size; i++)
{
if (matrix[i] != 0)
{
for (int j = 1; j < size; j++)
{
float* targetOffset = temporaryMatrix + (j - 1) * (size - 1);
float* sourceOffset = matrix + j * size;
int firstCopySize = i * sizeof *matrix;
int secondCopySize = (size - i - 1) * sizeof *matrix;
memcpy(targetOffset, sourceOffset, firstCopySize);
memcpy(targetOffset + i, sourceOffset + i + 1, secondCopySize);
}
float subdeterminant = determinant<size - 1>(temporaryMatrix);
total += matrix[i] * subdeterminant * sign;
}
sign *= -1;
}
return total;
}
template<>
float determinant<1>(float* matrix)
{
return matrix[0];
}
It seems to work pretty well on relatively small matrices (2x2, 3x3 and 4x4), so I'd expect it to work on 6x6 matrices too. I didn't conduct extensive tests however, and there's a possibility that it's broken.
The Problem
If I understood correctly the answer from the other question, the determinant should be 0 since the curves intersect. However, feeding my program the Sylvester matrix I made above, it's -2916.
Is it a mistake on my end or on their end? What's the correct way to find out if two cubic Bézier curves intersect?
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（8）
贝塞尔曲线的交点是通过（非常酷的）Asymptote矢量图形语言完成的：查找
intersect( )此处。尽管他们没有解释他们实际使用的算法，只是说它来自 p.在“The Metafont Book”的第 137 页中，它的关键似乎是贝塞尔曲线的两个重要属性（尽管我现在找不到该页面，但该网站的其他地方对此进行了解释）：
借助这两个属性和相交多边形的算法，您可以递归到任意精度:
bezInt(B1, B2):
bezInt(B1a, B2b) ||
bezInt(B1b, B2a) ||
bezInt(B1b, B2b).
如果曲线不相交，这会很快——这是通常的情况吗？
[编辑]看起来将贝塞尔曲线一分为二的算法称为de Casteljau 算法。
Intersection of Bezier curves is done by the (very cool) Asymptote vector graphics language: look for
intersect()here.Although they don't explain the algorithm they actually use there, except to say that it's from p. 137 of "The Metafont Book", it appears that the key to it is two important properties of Bezier curves (which are explained elsewhere on that site though I can't find the page right now):
With these two properties and an algorithm for intersecting polygons, you can recurse to arbitrary precision:
bezInt(B1, B2):
bezInt(B1a, B2b) ||
bezInt(B1b, B2a) ||
bezInt(B1b, B2b).
This will be fast if the curves don't intersect -- is that the usual case?
[EDIT] It looks like the algorithm for splitting a Bezier curve in two is called de Casteljau's algorithm.
如果您要为生产代码执行此操作，我建议使用贝塞尔裁剪算法。 此免费在线 CAGD 文本的第 7.7 节（pdf）对此进行了很好的解释，适用于任何贝塞尔曲线的次数，并且快速且鲁棒。
虽然从数学角度来看，使用标准寻根器或矩阵可能更直接，但贝塞尔剪裁相对容易实现和调试，并且实际上具有更少的浮点误差。这是因为每当它创建新数字时，它都会进行加权平均值（凸组合），因此没有机会根据噪声输入进行推断。
If you're doing this for production code, I'd suggest the Bezier clipping algorithm. It's explained well in section 7.7 of this free online CAGD text (pdf), works for any degree of Bezier curve, and is fast and robust.
While using standard rootfinders or matrices might be more straightforward from a mathematical point of view, Bezier clipping is relatively easy to implement and debug, and will actually have less floating point error. This is because whenever it's creating new numbers, it's doing weighted averages (convex combinations) so there's no chance of extrapolating based on noisy inputs.
您对行列式的解释是否基于此答案所附的第四条评论？如果是这样，我相信这就是错误所在。在此转载评论：
我不认为这部分有任何问题，但作者接着说：
我认为这是不正确的。两条曲线完全有可能在 t 值不同的位置相交，在这种情况下，即使矩阵具有非零行列式，也会存在相交。我相信这就是您的情况。
Are you basing your interpretation of the determinant on the 4th comment attached to this answer? If so, I believe that's where the mistake lies. Reproducing the comment here:
I don't see any problems with this part, but the author goes on to say:
I don't think that's correct. It's perfectly possible for the two curves to intersect in a location where the t values differ, and in that case, there will be an intersection even though the matrix has a non-zero determinant. I believe this is what's happening in your case.
这是一个难题。我会将 2 条贝塞尔曲线分成 5-10 条离散线段，然后只进行线与线的相交。
This is a hard problem. I would split each of the 2 Bezier curves into say 5-10 discrete line segments, and just do line-line intersections.
我不知道它有多快，但如果你有两条曲线 C1(t) 和 C2(k)，如果 C1(t) == C2(k)，它们就会相交。因此，对于两个变量 (t, k)，您有两个方程（每个 x 和每个 y）。您可以使用数值方法求解该系统，并且精度足够高。当你找到 t,k 参数时，你应该检查 [0, 1] 上是否有参数。如果是，它们在 [0, 1] 上相交。
I don't know how fast it will be, but if you have two curves C1(t) and C2(k) they intersect if C1(t) == C2(k). So you have two equations (per x and per y) for two variables (t, k). You can solve the system using numerical methods with enough for you accuracy. When you've found t,k parameters you should check if there is a parameter on [0, 1]. If it is they intersects on [0, 1].
我绝不是这类事情的专家，但我遵循一个很好的
http://cagd.cs.byu.edu/~557/text/ch7.pdf(存档副本）I'm by no way an expert on this kind of thing, but I follow a nice blog that talks a lot about curves. He has link to two nice articles talking about your problem (the second link has an interactive demonstration and some source code). Other people may have much better insight into the problem but I hope this helps!
http://cagd.cs.byu.edu/~557/text/ch7.pdf(archived copy)我想说，最简单也可能最快的答案是将它们细分为非常小的线，并找到曲线相交的点（如果它们确实相交）。
显然，蛮力答案很糟糕，请参阅 bo{4} 的答案，并且有很多线性几何和碰撞检测实际上会有很大帮助。
选择曲线所需的切片数。 100 之类的应该就很好了。
我们获取所有段并根据它们具有的 y 最大值对它们进行排序。然后，我们在列表中添加第二个标志，以获取该段的 y 最小值。
我们保留一个活动边的列表。
我们迭代 y 排序的段列表，当遇到前导段时，我们将其添加到活动列表中。当我们达到小 y 标记值时，我们会从活动列表中删除该段。
然后，我们可以简单地迭代整个段集，相当于一条扫描线，当我们简单地迭代列表时，单调增加 y 。我们迭代排序列表中的值，这通常会删除一个段并添加一个新段（或者对于拆分和合并节点，添加两个段或删除两个段）。从而保留相关细分市场的活跃列表。
当我们将新的活动线段添加到活动线段列表中时，我们运行快速失败相交检查，仅针对该线段和当前活动线段。
因此，当我们迭代曲线的采样线段时，我们始终准确地知道哪些线段是相关的。我们知道这些片段在 y 坐标上有重叠。我们可以快速失败任何与 x 坐标不重叠的新段。在极少数情况下，它们在 x 坐标上重叠，然后我们检查这些线段是否相交。
这可能会减少线路交叉点检查的数量，基本上减少交叉点的数量。
checkAgainstActive() 将简单地检查该段和活动列表中的任何段是否与其 x 坐标重叠，如果重叠，则对它们运行线相交检查，并采取适当的操作。
另请注意，这适用于任何数量的曲线、任何类型的曲线、任何顺序的曲线、任何混合形式。如果我们迭代整个段列表，它将找到每个交叉点。它将找到贝塞尔曲线与圆相交的每个点或十几个二次贝塞尔曲线彼此（或自身）的每个交点，并且所有这些都在同一瞬间完成。
对于细分算法，经常引用第7章文档注释：
“一旦一对曲线被足够细分，它们都可以通过线段近似到公差范围内”，
我们实际上可以跳过中间人。我们可以足够快地完成此操作，以便简单地将线段与曲线的误差量进行比较。最后，这就是典型答案的作用。
其次，请注意，这里碰撞检测的大部分速度增加，即基于段的最高 y 排序的有序列表添加到活动列表，以及从活动列表中删除最低 y，同样可以完成直接为贝塞尔曲线的船体多边形。我们的线段只是一个 2 阶多边形，但我们可以简单地处理任意数量的点，并以我们正在处理的曲线的任何顺序检查所有控制点的边界框。因此，我们将拥有一个活动外壳多边形点列表，而不是活动线段列表。我们可以简单地使用 De Casteljau 的算法来分割曲线，从而将这些曲线采样为细分曲线而不是线段。因此，我们不会得到 2 点，而是得到 4 点或 7 点或诸如此类的东西，并且运行相同的程序很容易快速失败。
在最大 y 处添加相关的点组，在最小 y 处将其删除，然后仅比较活动列表。这样我们就可以快速、更好的实现贝塞尔细分算法。通过简单地找到边界框重叠，然后细分那些重叠的曲线，并删除那些不重叠的曲线。同样，我们可以制作任意数量的曲线，甚至是从上一次迭代中的曲线细分出来的曲线。就像我们的线段近似一样，可以非常快速地求解数百条不同曲线（甚至不同阶）之间的每个交点位置。只需检查所有曲线以查看边界框是否重叠，如果重叠，我们将它们细分，直到我们的曲线足够小或者我们用完它们。
I would say that the easiest and likely fastest answer is to subdivide them into very small lines and find the points where the curves intersect, if they actually do.
Obviously the brute force answer is bad See bo{4}'s answer, and there's a lot of linear geometry and collision detection that will actually help quite a lot.
Pick the number of slices you want for the curves. Something like 100 should be great.
We take all the segments and sort them based on the largest value of y they have. We then, add a second flag in the list for the smallest value of y for that segment.
We keep a list of active edges.
We iterate through the y-sorted list of segment, when we encounter a leading segment we add it to the active list. When we hit the small-y flagged value, we remove that segment from the active list.
We then can simply iterate through the entire set of segments with what amounts to a scan line, increasing our y monotonically as we simply iterate the list. We iterate through the values in our sorted list, which will typically remove one segment and add a new segment (or for split and merge nodes, add two segments or remove two segments). Thereby keeping an active list of relevant segments.
We run the fast fail intersect check as we add a new active segment to the list of active segments, against only that segment and and the currently active segments.
So at all times we know exactly which line segments are relevant, as we iterate through the sampled segments of our curves. We know these segments share overlap in the y-coords. We can fast fail any new segment that does not overlap with regard to its x-coords. In the rare case that they overlap in the x-coords, we then check if these segments intersect.
This will likely reduce the number of line intersection checks to basically the number of intersections.
checkAgainstActive() would simply check if this segment and any segment in the active list overlap their x-coords, if they do, then run the line intersect check on them, and take the appropriate action.
Also note, this will work for any number of curves, any type of curves, any order of curves, in any mixture. And if we iterate through the entire list of segments it will find every intersection. It will find every point a Bezier intersects a circle or every intersection that a dozen quadratic Bezier curves have with each other (or themselves), and all in the same split second.
The oft cited Chapter 7 document notes, for the subdivision algorithm:
"Once a pair of curves has been subdivided enough that they can each be approximated by a line segment to within a tolerance"
We can literally just skip the middleman. We can do this fast enough so to simply compare line segments with a tolerable amount of error from the curve. In the end, that's what the typical answer does.
Secondly, note, that the bulk of the speed increase for the collision detection here, namely the ordered list of segments sorted based on their highest y to add to the active list, and lowest y to remove from the active list, can likewise be done for the hull polygons of the Bezier curve directly. Our line segment is simply a polygon of order 2, but we could do any number of points just as trivially, and checking the bounding box of all the control points in whatever order of curve we're dealing with. So rather than a list of active segments, we would have a list of active hull polygons points. We could simply use De Casteljau's algorithm to split the curves, thereby sampling as these as subdivided curves rather than line segments. So rather than 2 points we'd have 4 or 7 or whatnot, and run the same routine being quite apt towards fast failing.
Adding the relevant group of points at max y, removing it at min y, and comparing only the active list. We can thus quickly and better implement the Bezier subdivision algorithm. By simply finding bounding box overlap, and then subdividing those curves which overlapped, and removing those that did not. As, again, we can do any number of curves, even ones subdivided from curves in the previous iteration. And like our segment approximation quickly solve for every intersection position between hundreds of different curves (even of different orders) very quickly. Just check all curves to see if the bounding boxes overlap, if they do, we subdivide those, until our curves are small enough or we ran out of them.
是的，我知道，这个帖子被接受并关闭了很长一段时间，但是......
首先，我要感谢你，zneak，给我灵感。您的努力可以找到正确的方法。
其次，我对所接受的解决方案不太满意。这种类型用在我最喜欢的免费软件 IPE 中，它的 bugzilla 有很多关于交叉点问题的准确性和可靠性低的抱怨，我的也是其中之一。
缺少的技巧允许以zneak提出的方式解决问题：将一条曲线缩短一个因子k< em>>0，则西尔维斯特矩阵的行列式将等于0。显然，如果缩短的曲线会相交，那么原始曲线也会相交。现在，任务变为搜索 k 因子的合适值。这导致了求解 9 次单变量多项式的问题。该多项式的实根和正根负责潜在交点。 （这不足为奇，两条三次贝塞尔曲线最多可以相交 9 次。）执行最终选择以仅查找那些提供 0< 的 k 因子;=t<=1 对于两条曲线。
现在是 Maxima 代码，其中起点是由 zneak 提供的 8 个点：
此时，Maxima 回答：
让 Maxima 求解这个方程：
答案是：
现在是选择正确值的代码k：
马克西玛的回答：
不过，不仅仅是蜂蜜。 所提出的方法不可用
有人可能会问，为什么缩短只执行一次。这已经足够了，因为逆反定律是被发现的en passant，但这是另一个故事了。
Yes, I know, this thread is accepted and closed for a long time, but...
First, I'd like to thank you, zneak, for an inspiration. Your effort allows to find the right way.
Second, I was not quite happy with the accepted solution. This kind is used in my favorite freeware IPE, and its bugzilla is plenty of complains for a low accuracy and reliability about an intersection issue, my among them.
The missing trick which allows to solve the problem in a manner proposed by zneak : It is enough to shorten one of curves by a factor k>0, then the determinant of Sylvester matrix will be equal to zero. It is obvious, if a shortened curve will intersect, then original will too. Now, the task is changed into the search for a suitable value for k factor. This has led to the problem of solving an univariate polynomial of 9 degree. A real and positive roots of this polynomial are responsible for potential points of intersection. (This shouldn't be a surprise, two cubic Bezier curves can intersect up to 9 times.) The final selection is performed to find only those k factors, which provide 0<=t<=1 for both curves.
Now the Maxima code, where the starting point is set of 8 points provided by zneak :
At this point, Maxima answered:
Let Maxima solve this equation:
The answer is:
Now the code to select a right value of k:
Maxima's answer:
Theres not only a honey, though. The presented method is unusable, if:
One can ask, why a shortening is performed only once. It's enough because of reverse-inverse law, which was discovered en passant, but this is an another story.