代码之家 › 专栏 › 技术社区 › zneak

检查两条三次B_zier曲线是否相交

bezier intersection c++

zneak · 技术社区 · 15 年前

对于一个个人项目,我需要找出两条三次Bzier曲线是否相交。我不需要知道在哪里:我只需要知道他们是否知道。不过,我需要尽快完成。

我一直在清理这个地方,发现了一些资源。大多数情况下,这里有

因此,在我计算出什么是a Sylvester matrix 之后,什么是a constitute ,什么是a resultant”(结果)和 Why it's utility”(为什么它有用),I thought I figured how the solution works.然而,现实开始有所不同,但效果并不那么好。

乱来

我已经使用我的图形计算器绘制了两个相交的bzier样条(我们称之为b _{0
_{和b
_{1
_{)。它们的坐标如下(P
_{0
_{,P
_{1
_{,P
_{2
_{,P
_{3
_):}}}}}}}}}}}

(1,1)(2,4)(3,4)(4,3)
(3,5)(3,6)(0,1)(3,1)

结果如下:b_{0_{being the“horizontal”curve and b_{1_{the other one:。}}}}



按照上述问题的最高投票答案的指示,我减去了b_{0_{to b_{1_{。它给我留下了两个方程(x和y轴),根据我的计算器,它们是:}}}}

x=9t^3-9t^2-3t+2
Y=9t^3-9t^2-6t+4



西尔维斯特矩阵

从这个基础上,我建立了以下西尔维斯特矩阵:

9-9-3 2 0
0 9-9-3 2 0
0 0 9-9-3 2
9-9-6 4 0
0 9-9-6 4 0
0 0 9-9-6 4

之后,利用.HReF=“HTTP://E.WiKiTo.Org/Wik/LaPasePosiod”Re=“NoeFror”> Laplace展开< /A>:,计算了矩阵行列式的C++函数。

template<int size>
浮点行列式(浮点*矩阵)
{
总浮动=0;
浮点数=1;
浮动临时矩阵[大小-1)*(大小-1)]
对于(int i=0;i<size;i++)
{
如果(矩阵[I]!= 0)
{
对于(int j=1;j<size;j++)
{
float*targetOffset=临时矩阵+(j-1)*(大小-1);
float*sourceOffset=matrix+j*大小;
int firstcopysize=i*sizeof*矩阵;
int secondcopysize=(大小-i-1)*sizeof*矩阵;
memcpy(目标偏移量、源偏移量、FirstCopySize);
memcpy(targetoffset+i,sourceoffset+i+1,secondcopysize);
}
浮点子行列式=行列式大小-1>(临时矩阵);
合计+=矩阵[i]*子行列式*符号;
}
符号*=-1;
}
收益总额;
}

模板<gt;
浮点行列式<1>(浮点*矩阵)
{
返回矩阵[0];
}

它在相对较小的矩阵(2x2、3x3和4x4)上似乎工作得很好,所以我希望它也能在6x6矩阵上工作。但是,我没有进行过广泛的测试,而且有可能它已损坏。




问题

如果我正确理解另一个问题的答案,行列式应该是0,因为曲线相交。但是,给我的程序输入我上面制作的西尔维斯特矩阵,它是-2916。


是我的错还是他们的错?找出两条三次Bzier曲线是否相交的正确方法是什么?

我一直在清理这个地方,发现了一些资源。大多数情况下,this question here答案很有希望。

所以当我想到什么是Sylvester matrix,什么是determinant,什么是resultant和why it's useful我想我知道解决方案是如何工作的。然而,现实却不尽相同,而且效果也不太好。



乱来

我用我的绘图计算器画了两个B_)zier样条(我们称之为B_零和B_一)那是相交的。它们的坐标如下(p_零,P_一,P_二,P_三):

(1, 1) (2, 4) (3, 4) (4, 3)
(3, 5) (3, 6) (0, 1) (3, 1)


结果如下,b_零是“水平”曲线和b_一另一个:



按照上述问题的最高投票答案的指示,我减去了b。_零对乙_一. 它给我留下了两个方程(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++函数计算矩阵的行列式。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];
}


它在相对较小的矩阵(2x2、3x3和4x4)上似乎工作得很好,所以我希望它也能在6x6矩阵上工作。不过,我没有进行广泛的测试,而且有可能它坏了。



问题

如果我正确理解另一个问题的答案,行列式应该是0,因为曲线相交。但是,给我的程序输入我上面制作的西尔维斯特矩阵,它是-2916。

是我的错还是他们的错?找出两条三次Bzier曲线是否相交的正确方法是什么?

8 回复 | 直到 8 年前

j_random_hacker 15 年前

贝塞尔曲线的交集是由(非常酷) Asymptote 矢量图形语言:查找 intersect() here .

虽然他们没有解释他们在那里实际使用的算法,除了说它来自“元字体书”的第137页外,似乎它的关键是贝塞尔曲线的两个重要属性(在该网站的其他地方解释,尽管我现在找不到页面):

贝塞尔曲线总是包含在由它的4个控制点定义的边界框中。
贝塞尔曲线总是可以细分为任意的 T 值转换为2个子贝塞尔曲线

有了这两个属性和相交多边形的算法,您可以递归到任意精度:

贝辛特(B _一 ,B _二 ):

B箱(B _一 )相交框(B _二 )?
- 否:返回false。
- 是:继续。
IS区域(B框 _一 ))+区域(B框 _二 ))<阈值?
- 是:返回真。
- 否:继续。
拆分B _一进入B _1A 和B _1B 在 T = 0.5
拆分B _二进入B _2A 和B _2B 在 T = 0.5
返回框(B _1A ,B _2A ) 贝辛特(B _1A ,B _2B ) 贝辛特(B _1B ,B _2A ) 贝辛特(B _1B ,B _2B )

如果曲线不相交,这会很快——这是常见的情况吗?

[编辑] 看起来把贝塞尔曲线一分为二的算法叫做 de Casteljau's algorithm .

tfinniga 12 年前

如果您是为生产代码做这个的,我建议使用贝塞尔裁剪算法。解释得很清楚 section 7.7 of this free online CAGD text (pdf),适用于任何程度的贝塞尔曲线,并且速度快且健壮。

虽然从数学的角度来看,使用标准的rootfinder或矩阵可能更简单,但是贝塞尔剪辑相对容易实现和调试,并且实际上具有更少的浮点误差。这是因为每当它创建新的数字时,它都会进行加权平均(凸组合),因此没有机会根据噪声输入进行外推。

Community CDub 8 年前

是我的错还是他们的错?

你是根据附在第4条评论中的决定因素来解释的吗? this answer ?如果是这样,我相信这就是错误所在。在此处复制评论:

如果行列式为零,则 x和y中的根在*处完全相同 t的值,因此两条曲线的交点。(t) 可能不在间隔0..1内但是)。

我看不出这部分有什么问题,但作者接着说:

如果行列式为零,则可以确保曲线不会在任何地方相交。

我认为那不正确。两条曲线完全有可能在t值不同的位置相交,在这种情况下,即使矩阵具有非零行列式,也会有一个交点。我相信这就是你的案子。

Andrew 15 年前

我不知道有多快,但是如果你有两条曲线c1(t)和c2(k),它们相交,如果c1(t)==c2(k)。对于两个变量(t,k),有两个方程(每x和每y)。你可以用足够精确的数值方法来解这个系统。找到t,k参数后,应检查[0,1]上是否有参数。如果是,它们在[0,1]上相交。

GWW 12 年前

我决不是这类事情的专家,但我喜欢 blog 那是关于曲线的。他有两篇关于你的问题的好文章的链接(第二个链接有一个交互式演示和一些源代码)。其他人可能对这个问题有更好的了解,但我希望这有帮助!

http://cagd.cs.byu.edu/~557/text/ch7.pdf

bobobobo 12 年前

这是个难题。我将把两条贝塞尔曲线中的每一条分割成5-10个离散的线段,然后做直线交叉。

enter image description here

foreach SampledLineSegment line1 in Bezier1
    foreach SampledLineSegment line2 in Bezier2
        if( line1 intersects line2 )
            then Bezier1 intersects Bezier2

Tatarize 9 年前

我会说,最简单也可能最快的答案是将它们细分为非常小的线,并找到曲线相交的点(如果它们确实如此)。

public static void towardsCubic(double[] xy, double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double t) {
    double x, y;
    x = (1 - t) * (1 - t) * (1 - t) * x0 + 3 * (1 - t) * (1 - t) * t * x1 + 3 * (1 - t) * t * t * x2 + t * t * t * x3;
    y = (1 - t) * (1 - t) * (1 - t) * y0 + 3 * (1 - t) * (1 - t) * t * y1 + 3 * (1 - t) * t * t * y2 + t * t * t * y3;
    xy[0] = x;
    xy[1] = y;
}

public static void towardsQuad(double[] xy, double x0, double y0, double x1, double y1, double x2, double y2, double t) {
    double x, y;
    x = (1 - t) * (1 - t) * x0 + 2 * (1 - t) * t * x1 + t * t * x2;
    y = (1 - t) * (1 - t) * y0 + 2 * (1 - t) * t * y1 + t * t * y2;
    xy[0] = x;
    xy[1] = y;
}

很明显,蛮力的答案是不好的,看bo 4的答案,有很多线性几何和碰撞检测实际上会有很大帮助。

选择曲线所需的切片数。像100这样的应该很好。

我们取所有的段,并根据它们所拥有的y的最大值对它们进行排序。然后,在列表中为该段的最小值y添加第二个标志。

我们保留一个活动边的列表。

我们迭代Y排序的段列表,当我们遇到一个前导段时,我们将它添加到活动列表中。当我们点击小的Y标记值时,我们从活动列表中删除该段。

然后我们可以简单地用相当于一条扫描线的数据迭代整个段集,在简单地迭代列表时单调地增加y。我们迭代排序列表中的值,这通常会删除一个段并添加一个新段(或者对于拆分和合并节点,添加两个段或删除两个段)。从而保持相关段的活动列表。

当我们将一个新的活动段添加到活动段列表中时,只针对该段和当前活动段运行快速失败交叉检查。

因此,当我们在曲线的采样段中迭代时,我们始终确切地知道哪些线段是相关的。我们知道这些段在y坐标中共享重叠。我们可以快速失败任何新的部分不重叠的X坐标。在极少数情况下,它们在X坐标中重叠,然后我们检查这些段是否相交。

这可能会将线路交叉口检查的数量减少到基本上的交叉口数量。

foreach(segment in sortedSegmentList) {
    if (segment.isLeading()) {
        checkAgainstActives(segment);
        actives.add(segment);
    }
    else actives.remove(segment)
}

checkAgInstactive()只需检查此段和活动列表中的任何段是否与它们的X坐标重叠,如果重叠,则运行行交叉检查,并采取适当的操作。

还要注意,这对于任何数量的曲线、任何类型的曲线、任何顺序的曲线、任何混合曲线都有效。如果我们遍历整个段列表,它会找到每个交叉点。它会发现一个贝塞尔曲线上的每一个点都与一个圆相交,或是一打二次贝塞尔曲线上的每一个相交点(或它们本身),并且都在同一个时间间隔内。

经常引用的 Chapter 7 document 注:对于细分算法:

“一旦对一对曲线进行了足够的细分,每一条曲线都可以通过一个线段近似到公差范围内。”

我们可以跳过中间人。我们可以以足够快的速度来比较直线段和曲线的误差。最后,这就是典型的答案。

其次,请注意,在这里碰撞检测的速度增加的大部分,即根据要添加到活动列表中的最高Y和要从活动列表中删除的最低Y排序的分段列表,同样可以直接对贝塞尔曲线的外壳多边形进行。我们的直线段只是一个2级的多边形,但是我们可以做任何数量的点,就像琐碎一样,并检查所有控制点的边界框,不管我们处理的是什么样的曲线顺序。因此,我们将有一个活动的船体多边形点列表,而不是一个活动段列表。我们可以简单地使用de Casteljau的算法来分割曲线,从而作为细分曲线而不是直线段进行采样。因此,我们不需要2分,而是4分或7分或其他什么分数,并且运行相同的程序,很容易快速失败。

在最大Y处添加相关的点组,在最小Y处删除它,并仅比较活动列表。从而可以快速、更好地实现贝塞尔细分算法。通过简单地查找边界框重叠,然后细分那些重叠的曲线,并删除那些不重叠的曲线。同样,我们可以做任意数量的曲线,甚至是在前一次迭代中从曲线细分出来的曲线。像我们的分段近似法一样,快速求解数百条不同曲线(甚至不同阶数)之间的每个交叉点位置。只需检查所有曲线,看看边界框是否重叠,如果重叠,我们就细分这些曲线,直到我们的曲线足够小或用完它们为止。

BiÅy 8 年前

是的,我知道,这个线程被接受并关闭了很长一段时间,但是…

首先,我要感谢你, 泽纳克 为了灵感。你的努力可以找到正确的方法。

第二,我对被接受的解决方案不太满意。这种是在我最喜欢的免费软件IPE中使用的,它的Bugzilla对交叉点问题的准确性和可靠性都很低,这其中就包括。

缺失的技巧,允许以提出的方式解决问题 泽纳克 :只需将其中一条曲线缩短一个系数即可。 K >0,那么西尔维斯特矩阵的行列式等于零。很明显,如果缩短的曲线会相交,那么原始曲线也会相交。现在,任务变为搜索合适的值 K 因素。这就导致了求解9度单变量多项式的问题。此多项式的实根和正根负责 潜在的 交叉点。(这并不奇怪,两条三次贝塞尔曲线最多可以相交9次。)最后的选择是只找到那些 K 系数,提供0<。= T &两条曲线lt;=1。

现在是maxima代码,其中起始点是由 泽纳克 :

p0x:1; p0y:1;
p1x:2; p1y:4;
p2x:3; p2y:4;
p3x:4; p3y:3;

q0x:3; q0y:5;
q1x:3; q1y:6;
q2x:0; q2y:1;
q3x:3; q3y:1;

c0x:p0x;
c0y:p0y;
c1x:3*(p1x-p0x);
c1y:3*(p1y-p0y);
c2x:3*(p2x+p0x)-6*p1x;
c2y:3*(p2y+p0y)-6*p1y;
c3x:3*(p1x-p2x)+p3x-p0x;
c3y:3*(p1y-p2y)+p3y-p0y;

d0x:q0x;
d0y:q0y;
d1x:3*(q1x-q0x);
d1y:3*(q1y-q0y);
d2x:3*(q2x+q0x)-6*q1x;
d2y:3*(q2y+q0y)-6*q1y;
d3x:3*(q1x-q2x)+q3x-q0x;
d3y:3*(q1y-q2y)+q3y-q0y;

x:c0x-d0x + (c1x-d1x*k)*t+ (c2x-d2x*k^2)*t^2+ (c3x-d3x*k^3)*t^3;
y:c0y-d0y + (c1y-d1y*k)*t+ (c2y-d2y*k^2)*t^2+ (c3y-d3y*k^3)*t^3;

z:resultant(x,y,t);

此时,马克西玛回答说:

(%o35)â2*(1004*k^9â5049*k^8+5940*k^7â1689*k^6+10584*k^5â8235*k^4â2307*k^3+1026*k^2+108*k+76)

让马克西玛解这个方程:

rr: float( realroots(z,1e-20))

答案是:

(%o36) [k=â0.40256438624399,k=0.43261490325108,k=0.84718739982868,k=2.643321910825066,k=2.71772491293651]

现在选择正确值的代码 K :

for item in rr do ( 
  evk:ev(k,item),
  if evk>0 then (  
    /*print("k=",evk),*/ 
    xx:ev(x,item),  rx:float( realroots(xx,1e-20)),/*print("x(t)=",xx,"   roots: ",rx),*/
    yy:ev(y,item),  ry:float( realroots(yy,1e-20)),/*print("y(t)=",yy,"   roots: ",ry),*/
    for it1 in rx do (  t1:ev(t,it1),
    for it2 in ry do (  t2:ev(t,it2),
       dt:abs(t1-t2),
       if dt<1e-10 then (
         /*print("Common root=",t1,"  delta t=",dt),*/
         if (t1>0) and (t1<=1) then ( t2:t1*evk,
         if (t2>0) and (t2<=1) then (                           
                 x1:c0x + c1x*t1+ c2x*t1^2+ c3x*t1^3,
                 y1:c0y + c1y*t1+ c2y*t1^2+ c3y*t1^3,
                 print("Intersection point:     x=",x1, "      y=",y1)
)))))/*,disp ("-----")*/
));

马克西玛的回答:

"Intersection point:     x="1.693201254437358"      y="2.62375005067273
(%o37) done

不过,这里不仅有蜂蜜。如果出现以下情况,则所述方法不可用:

如果p0位于第二条曲线(或其延伸)上,则p0=q0,或者更一般地说。可以尝试交换曲线。
当两条曲线都属于一个K族时(例如,它们的无限延伸是相同的),这种情况极为罕见。
注意 (sqr(c3x)+sqr(c3y))=0 或 (sqr(d3x)+sqr(3y))=0 在这种情况下,二次曲线假装为三次贝塞尔曲线。

人们可以问,为什么缩短只执行一次。这就足够了,因为反逆定律被发现了 在过去 但这是另一个故事。

检查两条三次B_zier曲线是否相交

乱来

西尔维斯特矩阵

问题

乱来

西尔维斯特矩阵

问题

贝辛特(B 一 ,B 二 ):

贝辛特(B _一 ,B _二 ):