廉价的计算立方贝塞尔曲线长度的方法

Question

似乎不存在立方贝塞尔曲线长度的解析解，但这并不意味着没有编写廉价解决方案。 廉价指的是大约50-100 ns（或更少）的东西。

有人知道类似的解决方案吗？也许可以分为两类：

1）误差小于1％但代码运行速度较慢。 2）误差高达20％，但运行速度更快？

我在Google上浏览了一下，但没有找到任何看起来像是好的解决方案。 只有将曲线分割成N个线段并对N个平方根求和的方法 - 对于更高的精度来说太慢，并且对于2或3个线段可能太不准确了。

有更好的方法吗？

- user2214913

---

请查看http://stackoverflow.com/a/28764614/107090。 - lhf

大约长度而非精确长度..正在进行中..可用（除了便宜的配方外，它可能适用于许多其他情况）我有时会用它来计算轨迹等，知道给定轨迹的长度很重要，我需要尽可能快地得到它。 - user2214913

作为http://pomax.github.com/bezierinfo的作者，这里最重要的问题是你认为你需要近似而不是正确的弧长是为了什么，因为这决定了可接受的误差范围，以及哪些近似快捷方式甚至可以采用，而不会丢失所需的信息。 - Mike 'Pomax' Kamermans

或许可以请 @MarkRansom 做现实检验，大多数矢量插图软件都内置了弧长计算功能。例如，需要呈现虚线或点划线曲线的任何内容都会实现某种形式的弧长计算。这比你的问题所暗示的要普遍得多。 - Mike 'Pomax' Kamermans

你不会得到一个简单的解决方案。好的估算涉及数值方法，需要一定的数学素养才能理解。你的选择是接受挑战并提高技能，或者雇佣有这些技能的人。 - Codie CodeMonkey

Answer 1

另一种选项是将弧长估计为弦和控制网格的平均值。在实践中：

Bezier bezier = Bezier (p0, p1, p2, p3);

chord = (p3-p0).Length;
cont_net = (p0 - p1).Length + (p2 - p1).Length + (p3 - p2).Length;

app_arc_length = (cont_net + chord) / 2;

您可以将样条曲线段递归地划分为两个部分，并计算直到收敛为止的弧长。我亲自测试过，它实际上收敛得非常快速。我从这个论坛中得到了这个想法。

- Nic

---

我很想看到这种方法和压平方法之间的性能比较...一旦计算出曲线的LUT，压平就不需要重新计算任何东西，而这种递归方法要么必须事先存储每个LUT点处的插值点，占用大量内存，要么在每次迭代时计算这些点以获得新的外壳，使其变得更慢。 - Mike 'Pomax' Kamermans

Answer 2

最简单的算法：将曲线压平并计算欧几里得距离。只要你想要一个近似的弧长，这个解决方案就快速而廉价。给定曲线的坐标LUT（我假设你在谈论速度，所以我认为你使用这些，并且不会不断重新计算坐标），它是一个简单的for循环和累加。在通用代码中，使用一个dist函数来计算两点之间的欧几里得距离：

var arclength = 0,
    last=LUT.length-1,
    i;
for (i=0; i<last; i++) {
  arclength += dist(LUT[i], LUT[i+1]);
}

完成。现在的arclength是基于您的LUT可以形成的曲线最大分段数的近似弧长。想要更快的速度但潜在误差更大？控制分段数。

var arclength = 0,
    segCount = ...,
    last=LUT.length-2,
    step = last/segCount,
    s, i;
for (s=0; s<=segCount; s++) {
  i = (s*step/last)|0;
  arclength += dist(LUT[i], LUT[i+1]);
}

这是可能生成接近真实弧长值的最简单算法。要得到更好的结果，您需要使用更昂贵的数值方法（如Legendre-Gauss积分技术）。

如果想了解原因，请查看“A Primer on Bézier Curves”的弧长部分。

- Mike 'Pomax' Kamermans

Answer 3

在我的情况下，一个快速有效的方法是这样的。（为Unity3d改写成c#）"

public static float BezierSingleLength(Vector3[] points){
    var p0 = points[0] - points[1];
    var p1 = points[2] - points[1];
    var p2 = new Vector3();
    var p3 = points[3]-points[2];

    var l0 = p0.magnitude;
    var l1 = p1.magnitude;
    var l3 = p3.magnitude;
    if(l0 > 0) p0 /= l0;
    if(l1 > 0) p1 /= l1;
    if(l3 > 0) p3 /= l3;

    p2 = -p1;
    var a = Mathf.Abs(Vector3.Dot(p0,p1)) + Mathf.Abs(Vector3.Dot(p2,p3));
    if(a > 1.98f || l0 + l1 + l3 < (4 - a)*8) return l0+l1+l3;

    var bl = new Vector3[4];
    var br = new Vector3[4];

    bl[0] = points[0];
    bl[1] = (points[0]+points[1]) * 0.5f;

    var mid = (points[1]+points[2]) * 0.5f;

    bl[2] = (bl[1]+mid) * 0.5f;
    br[3] = points[3];
    br[2] = (points[2]+points[3]) * 0.5f;
    br[1] = (br[2]+mid) * 0.5f;
    br[0] = (br[1]+bl[2]) * 0.5f;
    bl[3] = br[0];

    return BezierSingleLength(bl) + BezierSingleLength(br);
}

- dan_wipf

Answer 4

我已经推导出了三点贝塞尔曲线长度的闭式表达式（如下）。我尚未尝试为4个或更多点计算闭式形式。这很可能难以表示和处理。然而，数值逼近技术，例如龙格-库塔积分算法（请参见此处的我的问答），通过使用弧长公式进行积分可以很好地工作。

这里是一些Java代码，用于计算三点贝塞尔曲线的弧长，其中点为a、b和c。

    v.x = 2*(b.x - a.x);
    v.y = 2*(b.y - a.y);
    w.x = c.x - 2*b.x + a.x;
    w.y = c.y - 2*b.y + a.y;

    uu = 4*(w.x*w.x + w.y*w.y);

    if(uu < 0.00001)
    {
        return (float) Math.sqrt((c.x - a.x)*(c.x - a.x) + (c.y - a.y)*(c.y - a.y));
    }

    vv = 4*(v.x*w.x + v.y*w.y);
    ww = v.x*v.x + v.y*v.y;

    t1 = (float) (2*Math.sqrt(uu*(uu + vv + ww)));
    t2 = 2*uu+vv;
    t3 = vv*vv - 4*uu*ww;
    t4 = (float) (2*Math.sqrt(uu*ww));

    return (float) ((t1*t2 - t3*Math.log(t2+t1) -(vv*t4 - t3*Math.log(vv+t4))) / (8*Math.pow(uu, 1.5)));

- user1300214

---

似乎对于坐标点 (-54.2, 45.5, 84.9)，(-39.2, 45.5, 84.9) 和 (-25.5, 45.5, 92.7) 返回了 NaN。 - Kari

@Kari，可能需要检查一些边缘情况，比如除以零或取零的对数等。 - user1300214

在我的情况下，最后一行是一个大负数的日志（-161）。 - Kari

@Kari 猜测可能需要取绝对值。我明天再确认一下。 - user1300214

Answer 5

public float FastArcLength()
{
    float arcLength = 0.0f;
    ArcLengthUtil(cp0.position, cp1.position, cp2.position, cp3.position, 5, ref arcLength);
    return arcLength;
}

private void ArcLengthUtil(Vector3 A, Vector3 B, Vector3 C, Vector3 D, uint subdiv, ref float L)
{
    if (subdiv > 0)
    {
        Vector3 a = A + (B - A) * 0.5f;
        Vector3 b = B + (C - B) * 0.5f;
        Vector3 c = C + (D - C) * 0.5f;
        Vector3 d = a + (b - a) * 0.5f;
        Vector3 e = b + (c - b) * 0.5f;
        Vector3 f = d + (e - d) * 0.5f;

        // left branch
        ArcLengthUtil(A, a, d, f, subdiv - 1, ref L);
        // right branch
        ArcLengthUtil(f, e, c, D, subdiv - 1, ref L);
    }
    else
    {
        float controlNetLength = (B-A).magnitude + (C - B).magnitude + (D - C).magnitude;
        float chordLength = (D - A).magnitude;
        L += (chordLength + controlNetLength) / 2.0f;
    }
}

- Tyler Heers

Answer 6

首先，你应该了解贝塞尔曲线中使用的算法， 当我在编写一个充满图形材料的C#程序时，我使用了贝塞尔曲线，并且很多时候我需要在贝塞尔曲线中找到一个点的坐标，这似乎在第一眼看起来是不可能的。所以我做的事情是在我的项目中编写了一个自定义数学类中的三次贝塞尔函数。现在我将与您分享代码。

    //--------------- My Costum Power Method ------------------\\

public static float FloatPowerX(float number, int power)
        {
            float temp = number;
            for (int i = 0; i < power - 1; i++)
            {
                temp *= number;
            }
            return temp;
        }

    //--------------- Bezier Drawer Code Bellow ------------------\\

public static void CubicBezierDrawer(Graphics graphics, Pen pen, float[] startPointPixel, float[] firstControlPointPixel
                , float[] secondControlPointPixel, float[] endPointPixel)
        {
            float[] px = new float[1111], py = new float[1111];
            float[] x = new float[4] { startPointPixel[0], firstControlPointPixel[0], secondControlPointPixel[0], endPointPixel[0] };
            float[] y = new float[4] { startPointPixel[1], firstControlPointPixel[1], secondControlPointPixel[1], endPointPixel[1] };
        int i = 0;

        for (float t = 0; t <= 1F; t += 0.001F)
        {
            px[i] = FloatPowerX((1F - t), 3) * x[0] + 3 * t * FloatPowerX((1F - t), 2) * x[1] + 3 * FloatPowerX(t, 2) * (1F - t) * x[2] + FloatPowerX(t, 3) * x[3];
            py[i] = FloatPowerX((1F - t), 3) * y[0] + 3 * t * FloatPowerX((1F - t), 2) * y[1] + 3 * FloatPowerX(t, 2) * (1F - t) * y[2] + FloatPowerX(t, 3) * y[3];
            graphics.DrawLine(pen, px[i - 1], py[i - 1], px[i], py[i]);
            i++;
        }
    }

如您所见，这是贝塞尔函数的工作方式，并且它绘制与 Microsoft 贝塞尔函数相同的 Bezier（我已经测试过）。您可以通过增加数组大小和计数器大小或绘制椭圆而使其更加精确，而不是绘制线条等...。所有这些都取决于您的需求和所需精度级别等等。

回到主要目标，问题是如何计算长度？

好吧，答案是我们有大量的点，每个点都有 x 坐标和 y 坐标，这使我们想起了一个三角形形状，特别是一个直角三角形形状。因此，如果我们有点 p1 和 p2，我们可以将它们之间的距离计算为直角三角形弦长。正如我们在学校的数学课堂上记得的那样，在 ABC 类型的直角三角形中，弦长是 -> Sqrt(Angle's FrontCostalLenght ^ 2 + Angle's SideCostalLeghth ^ 2);

并且在当前点和当前点之前的最后一个点（例如 p[i-1] 和 p[i]）之间计算长度并将它们全部存储在一个变量中。让我们在下面的代码中显示它。

//--------------- My Costum Power Method ------------------\\

public static float FloatPower2(float number)
        {
            return number * number;
        }

//--------------- My Bezier Lenght Calculator Method ------------------\\

public static float CubicBezierLenghtCalculator(float[] startPointPixel
            , float[] firstControlPointPixel, float[] secondControlPointPixel, float[] endPointPixel)
        {
            float[] tmp = new float[2];
            float lenght = 0;
            float[] px = new float[1111], py = new float[1111];

            float[] x = new float[4] { startPointPixel[0], firstControlPointPixel[0]
                , secondControlPointPixel[0], endPointPixel[0] };

            float[] y = new float[4] { startPointPixel[1], firstControlPointPixel[1]
                , secondControlPointPixel[1], endPointPixel[1] };

            int i = 0;

            for (float t = 0; t <= 1.0; t += 0.001F)
            {
                px[i] = FloatPowerX((1.0F - t), 3) * x[0] + 3 * t * FloatPowerX((1.0F - t), 2) * x[1] + 3F * FloatPowerX(t, 2) * (1.0F - t) * x[2] + FloatPowerX(t, 3) * x[3];
                py[i] = FloatPowerX((1.0F - t), 3) * y[0] + 3 * t * FloatPowerX((1.0F - t), 2) * y[1] + 3F * FloatPowerX(t, 2) * (1.0F - t) * y[2] + FloatPowerX(t, 3) * y[3];
                if (i > 0)
                {
                    tmp[0] = Math.Abs(px[i - 1] - px[i]);// calculating costal lenght
                    tmp[1] = Math.Abs(py[i - 1] - py[i]);// calculating costal lenght
                    lenght += (float)Math.Sqrt(FloatPower2(tmp[0]) + FloatPower2(tmp[1]));// calculating the lenght of current RightTriangle Chord  & add it each time to variable
                }
                i++;
            }
            return lenght;
        }

如果您希望进行更快的计算，只需减少px和py数组长度和循环次数即可。

我们还可以通过将px和py减少到数组长度为1或创建一个简单的双精度变量来减少内存需求，但由于出现了条件情况，增加了我们的大O值，我没有这样做。

希望它对您有所帮助。如果有其他问题，请随时提问。 此致，敬礼！Heydar - 伊朗伊斯兰共和国。

- Heydar