计算三次贝塞尔曲线的最快方法？

Posted 2023-02-21

【中文标题】计算三次贝塞尔曲线的最快方法？

【英文标题】：Fastest way to calculate cubic bezier curves?

【发布时间】：2010-07-07 05:17:34

【问题描述】：

现在我是这样计算的：

    double dx1 = a.RightHandle.x - a.UserPoint.x;
    double dy1 = a.RightHandle.y - a.UserPoint.y;
    double dx2 = b.LeftHandle.x - a.RightHandle.x;
    double dy2 = b.LeftHandle.y - a.RightHandle.y;
    double dx3 = b.UserPoint.x - b.LeftHandle.x;
    double dy3 = b.UserPoint.y - b.LeftHandle.y;

    float len = sqrt(dx1 * dx1 + dy1 * dy1) + 
        sqrt(dx2 * dx2 + dy2 * dy2) + 
        sqrt(dx3 * dx3 + dy3 * dy3);




    int NUM_STEPS =  int(len * 0.05);

    if(NUM_STEPS > 55)
    
        NUM_STEPS = 55;
    
    double subdiv_step  = 1.0 / (NUM_STEPS + 1);
    double subdiv_step2 = subdiv_step*subdiv_step;
    double subdiv_step3 = subdiv_step*subdiv_step*subdiv_step;

    double pre1 = 3.0 * subdiv_step;
    double pre2 = 3.0 * subdiv_step2;
    double pre4 = 6.0 * subdiv_step2;
    double pre5 = 6.0 * subdiv_step3;



    double tmp1x = a.UserPoint.x - a.RightHandle.x * 2.0 + b.LeftHandle.x;
    double tmp1y = a.UserPoint.y - a.RightHandle.y  * 2.0 + b.LeftHandle.y;

    double tmp2x = (a.RightHandle.x - b.LeftHandle.x)*3.0 - a.UserPoint.x + b.UserPoint.x;
    double tmp2y = (a.RightHandle.y - b.LeftHandle.y)*3.0 - a.UserPoint.y + b.UserPoint.y;

    double fx = a.UserPoint.x;
    double fy = a.UserPoint.y;

    //a user
    //a right
    //b left
    //b user

    double dfx = (a.RightHandle.x - a.UserPoint.x)*pre1 + tmp1x*pre2 + tmp2x*subdiv_step3;
    double dfy = (a.RightHandle.y - a.UserPoint.y)*pre1 + tmp1y*pre2 + tmp2y*subdiv_step3;

    double ddfx = tmp1x*pre4 + tmp2x*pre5;
    double ddfy = tmp1y*pre4 + tmp2y*pre5;

    double dddfx = tmp2x*pre5;
    double dddfy = tmp2y*pre5;

    int step = NUM_STEPS;



    while(step--)
    


        fx   += dfx;
        fy   += dfy;
        dfx  += ddfx;
        dfy  += ddfy;
        ddfx += dddfx;
        ddfy += dddfy;
        temp[0] = fx;
        temp[1] = fy;
        Contour[currentcontour].DrawingPoints.push_back(temp);
    


    temp[0] = (GLdouble)b.UserPoint.x;
    temp[1] = (GLdouble)b.UserPoint.y;
    Contour[currentcontour].DrawingPoints.push_back(temp);

我想知道是否有更快的方法来插入三次贝塞尔曲线？

谢谢

【参考方案1】：

查看forward differencing 以获得更快的方法。必须小心处理舍入错误。

adaptive subdivision 方法，通过一些检查，可以快速而准确。

【讨论】：

欢迎提供解决方案的链接，但请确保您的答案在没有它的情况下有用：add context around the link 这样您的其他用户就会知道它是什么以及为什么会出现，然后引用最相关的内容您链接到的页面的一部分，以防目标页面不可用。 Answers that are little more than a link may be deleted.

在某些情况下，前向差分在实践中不起作用（特别是在曲线不能用锐角或微小循环微分的点周围）。递归拆分（生成折线）始终有效且速度非常快（并且递归循环自动受平面度参数的约束，该参数也可以使用基本算术计算）。然后，您通过在每个直线段上使用经典的 Bresenham 算法来完成绘图（预置了起点的校正因子，在递归拆分期间计算）。将任何三次贝塞尔曲线压平为 1%：最多得到 628 分。

请注意，对于 1% 的任何分段三次贝塞尔曲线，您最多获得的 628 点仅适用于无限分辨率渲染的情况。这些点在曲率最大的区域周围分布得更密集（在曲线的中间）：当它们的距离小于半个像素时，它们中的大部分会被消除。在实践中，您将获得最多包含几十个点的折线，并且将折线展平为 1% 非常好，因为您仍然会在每个像素处获得跳跃（在对角线和正交方向之间）（以及最终的“错误" 保持在半个像素以下）。

上面显示的朴素算法（使用 55 步），仅在曲线段中间（曲率最大）提供约 12% 的展平：如果渲染一个与直线在 9 像素以上，您会看到大于 1 像素的误差！要在任何分辨率下（偏差小于 1/2 像素）获得 1% 的误差，您需要多 12 倍的点（然后通过平均所有点来消除太短的段）。值得注意的是，曲线中间的点是这个平均值最常被“消除”的点

【参考方案2】：

还有一点也很重要，那就是您正在使用大量固定长度的直线段来逼近您的曲线。这在曲线几乎笔直的区域效率低下，并且可能导致曲线非常弯曲的令人讨厌的多角线。对于高曲率和低曲率，没有简单的折衷方案。

要解决这个问题，您可以动态细分曲线（例如，在中点将其分成两部分，然后查看两条线段是否在曲线的合理距离内。如果一个线段是好的适合曲线，停在那里；如果不是，则以相同的方式细分并重复）。以这种方式对曲线进行采样时，您必须小心地对其进行足够的细分，以免错过任何局部（小）特征。

这并不总是“更快”地绘制你的曲线，但它会保证它总是看起来不错，同时使用达到该质量所需的最少线段数。

一旦你“很好”地绘制了曲线，你就可以看看如何“更快”地进行必要的计算。

【参考方案3】：

实际上，您应该继续拆分，直到连接曲线上的点（末端节点）及其最远控制点的两条线“足够平坦”： - 完全对齐或 - 它们的交点位于与两端节点的“平方距离”低于一半“平方像素”的位置） - 请注意，您不需要计算实际距离，因为它需要计算平方根，即慢）

当你遇到这种情况时，忽略控制点并用直线连接两个端点。

这样更快，因为使用经典的 Bresenham 算法，您将很快得到可以直接绘制的直线段，就好像它们是直线一样。

注意：您应该考虑端点的小数位，以正确设置误差变量的初始值，累积差异并由增量 Bresenham 算法使用，以获得更好的结果（尤其是在绘制最终段时）与水平或垂直或两条对角线非常接近）；否则你会得到可见的人工制品。

用于在整数网格上对齐的点之间画线的经典 Bresenham 算法将这个错误变量初始化为零，以表示第一个端节点的位置。但是 Bresenham 算法的一个小修改将两个距离变量和误差值简单地放大了 2 的恒定幂，然后对未缩放的 x 或 y 变量使用 0/+1 增量。

错误变量的高位还允许您计算一个 alpha 值，该值可用于绘制两个具有正确 alpha 阴影的堆叠像素。在大多数情况下，您的图像最多使用 8 位颜色平面，因此您不需要超过 8 位的额外精度来计算误差值，并且放大倍数可以限制为 256 倍：您可以使用它画出“平滑”的线条。

但您可以将自己的缩放因子限制为 16（4 位）：您必须绘制的典型位图图像不是非常宽，它们的分辨率远低于 +/- 20 亿（带符号的 32 位的限制）整数）：当您将​​坐标放大 16 倍时，它将保持 28 位可用，但您应该已经将几何图形“剪切”到要渲染的位图的视图区域，并且错误变量Bresenham 算法在所有情况下都将保持在 56 位以下，并且仍然适合 64 位整数。

如果您的错误变量是 32 位，则在最坏的情况下，您必须将缩放坐标限制在 2^15（不超过 15 位）以下（否则 Bresenham 使用的错误变量的符号测试将不起作用由于最坏情况下的整数溢出），并且使用 16（4 位）的放大系数，您将被限制绘制宽度或高度不大于 11 位的图像，即 2048x2048 图像。

但是，如果您的绘图区域实际上低于 2048x2048 像素，则绘制由绘图颜色的 16 个 alpha 阴影值平滑的线条是没有问题的（要绘制 alpha 阴影像素，您需要读取图像中的原始像素值在混合 alpha 阴影颜色之前，除非计算的阴影对于您不需要绘制的第一个堆叠像素为 0%，而对于您可以直接用普通绘制颜色覆盖的第二个堆叠像素为 100%）

如果您的计算图像还包含一个 alpha 通道，则您的绘制颜色也可以有自己的 alpha 值，您需要对它进行着色并与要绘制的像素的 alpha 值相结合。但是您不需要任何中间缓冲区来绘制线条，因为您可以直接在目标缓冲区中绘制。

使用 Bresenham 算法使用的误差变量，完全没有因舍入误差引起的问题，因为它们已被此变量考虑在内。所以要正确设置它的初始值（或者，在开始递归细分之前，只需将所有坐标放大 16 倍，Bresenham 算法本身的样条线要慢 16 倍）。

【参考方案4】：

请注意如何计算“足够平坦”。 “平面度”是两个连续段之间的最小绝对角度（0 到 180° 之间）的度量，但您不需要计算实际角度，因为此平面度也相当于为 余弦设置最小值 他们的相对角度。

该余弦值也不需要直接计算，因为您实际上只需要两个向量的向量积，并将其与它们长度最大值的平方进行比较。

还要注意，“余弦的平方”也是“1 减去正弦的平方”。最大平方余弦值也是最小平方正弦值......现在您知道使用哪种“向量积”：计算最快和最简单的是标量积，其平方与两者的平方正弦成正比向量和两个向量的平方长度的乘积。

因此，检查曲线是否“足够平坦”很简单：计算两个标量积之间的比率，并查看该比率是否低于“平坦度”常数值（最小平方正弦）。没有除以零，因为您将确定两个向量中的哪一个是最长的，如果这个向量的平方长度低于 1/4，那么您的曲线已经足够平坦，可以用于渲染分辨率；否则检查最长和最短向量之间的比率（由包含端点和控制点的凸包的交叉对角线形成）：

二次贝塞尔曲线，凸包是一个三角形，你选择两对

对于三次贝塞尔曲线，凸包是一个 4 边凸多边形，对角线可以将端点与两个控制点之一连接起来，也可以将两个端点和两个控制点连接在一起你有六种可能

使用组合给出第一个端点和其他三个点之一之间的第一个向量的最大长度，第二个向量连接另外两个点）：

您需要确定从一个端点或控制点开始到序列中下一个控制点或端点的分段的“最小平方长度”。 （在二次贝塞尔曲线中，您只需比较两个线段，在二次贝塞尔曲线中，您检查 3 段）

如果这个“最小平方长度”低于 1/4，您可以停在那里，曲线“足够平坦”。

然后确定从一个端点开始到另一个端点或控制点中的任何一个的线段的“最大平方长度”（使用二次贝塞尔曲线，您可以安全地使用与上面，使用三次贝塞尔曲线，您丢弃了上面使用的连接 2 个控制点的 3 个线段之一，但您添加了连接两个端节点的线段。

然后检查“最小平方长度”是否低于常数“平坦度”（最小平方正弦）乘以“最大平方长度”的乘积（如果是，则曲线“足够平坦”。

在这两种情况下，当您的曲线“足够平坦”时，您只需绘制连接两个端点的线段。否则，您将递归地拆分样条曲线。

您可以对递归进行限制，但在实践中，除非曲线的凸包在非常大的绘制区域中覆盖了非常大的区域，否则永远不会达到此限制；即使有 32 级回避，它也永远不会在对角线短于 2^32 像素的矩形绘制区域中爆炸（只有当您在具有浮点的几乎无限空间中分割“虚拟贝塞尔曲线”时才会达到限制坐标，但您不打算直接绘制它，因为在这样的空间中您不会有 1/2 像素的限制，并且只有当您为“平面度”设置了 extreme 值时您的“最小平方正弦”常数参数为 1/2^32 或更低）。