Optimize quadratic curve tracing using numeric methods

Question

I am trying to trace quadratic bezier curves, placing "markers" at a given step length distance. Tried to do it a naive way:

    const p = toPoint(map, points[section + 1]);
    const p2 = toPoint(map, points[section]);
    const {x: cx, y: cy} = toPoint(map, cp);
    const ll1 = toLatLng(map, p),
    ll2 = toLatLng(map, p2),
    llc = toLatLng(map, { x: cx, y: cy });
    const lineLength = quadraticBezierLength(
      ll1.lat,
      ll1.lng,
      llc.lat,
      llc.lng,
      ll2.lat,
      ll2.lng
    );
    for (let index = 0; index < Math.floor(lineLength / distance); index++) {
      const t = distance / lineLength;
      const markerPoint = getQuadraticPoint(
        t * index,
        p.x,
        p.y,
        cx,
        cy,
        p2.x,
        p2.y
      );
      const markerLatLng = toLatLng(map, markerPoint);

      markers.push(markerLatLng);
    }

This approach does not work since the correlation of a quadratic curve between t and L is not linear. I could not find a formula, that would give me a good approximation, so looking at solving this problem using numeric methods [Newton]. One simple option that I am considering is to split the curve into x [for instance 10] times more pieces than needed. After that, using the same quadraticBezierLength() function calculate the distance to each of those points. After this, chose the point so that the length is closest to the distance * index.

This however would be a huge overkill in terms of algorithm complexity. I could probably start comparing points for index + 1 from the subset after/without the point I selected already, thus skipping the beginning of the set. This would lower the complexity some, yet still very inefficient.

Any ideas and/or suggestions?

Ideally, I want a function that would take d - distance along the curve, p0, cp, p1 - three points defining a quadratic bezier curve and return an array of coordinates, implemented with the least complexity possible.

Hypermetropia

Answer 1

OK I found analytic formula for 2D quadratic bezier curve in here:

• Calculate the length of a segment of a quadratic bezier

So the idea is simply binary search the parameter t until analytically obtained arclength matches wanted length...

C++ code:

//---------------------------------------------------------------------------
float x0,x1,x2,y0,y1,y2;    // control points
float ax[3],ay[3];          // coefficients
//---------------------------------------------------------------------------
void get_xy(float &x,float &y,float t)  // get point on curve from parameter t=<0,1>
    {
    float tt=t*t;
    x=ax[0]+(ax[1]*t)+(ax[2]*tt);
    y=ay[0]+(ay[1]*t)+(ay[2]*tt);
    }
//---------------------------------------------------------------------------
float get_l_naive(float t)      // get arclength from parameter t=<0,1>
    {
    // naive iteration
    float x0,x1,y0,y1,dx,dy,l=0.0,dt=0.001;
    get_xy(x1,y1,t);
    for (int e=1;e;)
        {
        t-=dt; if (t<0.0){ e=0; t=0.0; }
        x0=x1; y0=y1; get_xy(x1,y1,t);
        dx=x1-x0; dy=y1-y0;
        l+=sqrt((dx*dx)+(dy*dy));
        }
    return l;
    }
//---------------------------------------------------------------------------
float get_l(float t)        // get arclength from parameter t=<0,1>
    {
    // analytic fomula from: https://mcmap.net/q/541527/-calculate-the-length-of-a-segment-of-a-quadratic-bezier
    float ax,ay,bx,by,A,B,C,b,c,u,k,cu,cb;
    ax=x0-x1-x1+x2;
    ay=y0-y1-y1+y2;
    bx=x1+x1-x0-x0;
    by=y1+y1-y0-y0;
    A=4.0*((ax*ax)+(ay*ay));
    B=4.0*((ax*bx)+(ay*by));
    C=     (bx*bx)+(by*by);
    b=B/(2.0*A);
    c=C/A;
    u=t+b;
    k=c-(b*b);
    cu=sqrt((u*u)+k);
    cb=sqrt((b*b)+k);
    return 0.5*sqrt(A)*((u*cu)-(b*cb)+(k*log(fabs((u+cu))/(b+cb))));
    }
//---------------------------------------------------------------------------
float get_t(float l0)       // get parameter t=<0,1> from arclength
    {
    float t0,t,dt,l;
    for (t=0.0,dt=0.5;dt>1e-10;dt*=0.5)
        {
        t0=t; t+=dt;
        l=get_l(t);
        if (l>l0) t=t0;
        }
    return t;
    }
//---------------------------------------------------------------------------
void set_coef()             // compute coefficients from control points
    {
    ax[0]=                  (    x0);
    ax[1]=        +(2.0*x1)-(2.0*x0);
    ax[2]=(    x2)-(2.0*x1)+(    x0);
    ay[0]=                  (    y0);
    ay[1]=        +(2.0*y1)-(2.0*y0);
    ay[2]=(    y2)-(2.0*y1)+(    y0);
    }
//---------------------------------------------------------------------------

Usage:

1. set control points x0,y0,...
2. then you can use t=get_t(wanted_arclength) freely

In case you want to use get_t_naive and or get_xy you have to call set_coef first

In case you want to tweak speed/accuracy you can play with the target accuracy of binsearch currently set to1e-10

Here optimized (merged get_l,get_t functions) version:

//---------------------------------------------------------------------------
float get_t(float l0)       // get parameter t=<0,1> from arclength
    {
    float t0,t,dt,l;
    float ax,ay,bx,by,A,B,C,b,c,u,k,cu,cb,cA;
    // precompute get_l(t) constants
    ax=x0-x1-x1+x2;
    ay=y0-y1-y1+y2;
    bx=x1+x1-x0-x0;
    by=y1+y1-y0-y0;
    A=4.0*((ax*ax)+(ay*ay));
    B=4.0*((ax*bx)+(ay*by));
    C=     (bx*bx)+(by*by);
    b=B/(2.0*A);
    c=C/A;
    k=c-(b*b);
    cb=sqrt((b*b)+k);
    cA=0.5*sqrt(A);
    // bin search t so get_l == l0
    for (t=0.0,dt=0.5;dt>1e-10;dt*=0.5)
        {
        t0=t; t+=dt;
        // l=get_l(t);
        u=t+b; cu=sqrt((u*u)+k);
        l=cA*((u*cu)-(b*cb)+(k*log(fabs((u+cu))/(b+cb))));
        if (l>l0) t=t0;
        }
    return t;
    }
//---------------------------------------------------------------------------

Sleekit

I don't really care for cubic, nor distance along straight lines. I am interested in quadratic curves. – Hypermetropia

appreciate the attempt. yet, while you tried to answer a different question, I believe that the solution you offered for the problem you tried to solve is far from optimal. Either way, it is irrelevant to my task, unfortunately. – Hypermetropia

@IgorShmukler I found correct formula with t ... I completely reedited my answer again removed all obsolete stuff and left just binary search based solution ... its working accurately on curves I tried... – Sleekit

Answer 2

For now, I came up with the below:

    for (let index = 0; index < Math.floor(numFloat * times); index++) {
      const t = distance / lineLength / times;
      const l1 = toLatLng(map, p), lcp = toLatLng(map, new L.Point(cx, cy));
      const lutPoint = getQuadraticPoint(
        t * index,
        p.x,
        p.y,
        cx,
        cy,
        p2.x,
        p2.y
      );
      const lutLatLng = toLatLng(map, lutPoint);
      const length = quadraticBezierLength(l1.lat, l1.lng, lcp.lat, lcp.lng, lutLatLng.lat, lutLatLng.lng);
      lut.push({t: t * index, length});
    }
    const lut1 = lut.filter(({length}) => !isNaN(length));
    console.log('lookup table:', lut1);
    for (let index = 0; index < Math.floor(numFloat); index++) {
      const t = distance / lineLength;
      // find t closest to distance * index
      const markerT = lut1.reduce((a, b) => {
        return a.t && Math.abs(b.length - distance * index) < Math.abs(a.length - distance * index) ? b.t : a.t || 0;
      });
      const markerPoint = getQuadraticPoint(
        markerT,
        p.x,
        p.y,
        cx,
        cy,
        p2.x,
        p2.y
      );
      const markerLatLng = toLatLng(map, markerPoint);
    }

I think only that my Bezier curve length is not working as I expected.

function quadraticBezierLength(x1, y1, x2, y2, x3, y3) {
  let a, b, c, d, e, u, a1, e1, c1, d1, u1, v1x, v1y;

  v1x = x2 * 2;
  v1y = y2 * 2;
  d = x1 - v1x + x3;
  d1 = y1 - v1y + y3;
  e = v1x - 2 * x1;
  e1 = v1y - 2 * y1;
  c1 = a = 4 * (d * d + d1 * d1);
  c1 += b = 4 * (d * e + d1 * e1);
  c1 += c = e * e + e1 * e1;
  c1 = 2 * Math.sqrt(c1);
  a1 = 2 * a * (u = Math.sqrt(a));
  u1 = b / u;
  a = 4 * c * a - b * b;
  c = 2 * Math.sqrt(c);
  return (
    (a1 * c1 + u * b * (c1 - c) + a * Math.log((2 * u + u1 + c1) / (u1 + c))) /
    (4 * a1)
  );
}

I believe that the full curve length is correct, but the partial length that is being calculated for the lookup table is wrong.

Hypermetropia

The length of the entire curve calculation is working fine. The problem, it seems that I am using it to calculate the wrong thing. I find a point on the curve, then assume that the new curve made from start, control and this newly found point will form a shorter part of the same curve. That is wrong, I understand. – Hypermetropia

Did you actually test your code? I built a similar function, and I am not getting the correct results. – Hypermetropia

I get values that are not even close. The length of the entire curve for my case is around 218, while the result for t = 1 is over 50,000. – Hypermetropia

In JavaScript, everything is a float, so there is no need to convert integers to floats. It is a problem for most calculations, but works fine here. BTW, why do you need fabs in your code? – Hypermetropia

this is on what I test i.sstatic.net/ObbQi.png the points are scaled with window size so I can quickly test while resizing window. As you can see the lengths for t=0.75 matches ... also get_l(get_t(200))=200.000005 which is quite good considering I use only floats – Sleekit

the fabs is there because its in the formula I linked ... you know those || near k*log|...| means absolute value in math syntax – Sleekit

I did two calculation, full arc length with: x1: 286, y1: 389, x2: 324.98205596852745, y2: 246.72650759509187, x3: 393.5, y3: 265.5 and the result - line length: 181.71253724569695; then quadraticBezierArcLength(1, 286, 389, 324.98205596852745, 246.72650759509187, 393.5, 265.5) and the result is 41919.46098139183. Wondering how it is possible that you are getting the correct results. – Hypermetropia

I hardcoded your points, used t=1.0 and l=100 ... i.sstatic.net/eHO0k.png ... heh the get_l(200) print should be 100 its hardcoded string :) ... try to copy paste my get_t and get_l functions ... just add the x0,y0,x1,y1,x2,y2 to the calling parameters ... also your quadraticBezierArcLength have 6 parameters how can you call it with 7? that is syntax error ... or you have different code than posted here – Sleekit

@Sleekit The log was fine, but there was an error somewhere. I tried your version, which is very similar to mine. Yet, your version is working correctly, while mine was not. Thank you so much. Will try to figure out what is wrong. – Hypermetropia

btw you can optimize further by merging get_l into get_t as get_l have many constants (not depending on t) that could be computed just once ... will add code to my answer in a minute – Sleekit

yes, a good point. – Hypermetropia

ok I mad a minor change ... there still are few terms that are constant in there but it would be a constants mess ... I think I leave the code as is from now on ... – Sleekit

Answer 3

If I am right, you want points at equally spaced points in terms of curvilinear abscissa (rather than in terms of constant Euclidean distance, which would be a very different problem).

Computing the curvilinear abscissa s as a function of the curve parameter t is indeed an option, but that leads you to the resolution of the equation s(t) = Sk/n for integer k, where S is the total length (or s(t) = kd if a step is imposed). This is not convenient because s(t) is not available as a simple function and is transcendental.

A better method is to solve the differential equation

dt/ds = 1/(ds/dt) = 1/√(dx/dt)²+(dy/dt)²

using your preferred ODE solver (RK4). This lets you impose your fixed step on s and is computationally efficient.

Thibaud

Did not fully understand your post. Hence, I am going to guess what you might have had in mind. s(t) is easy to calculate. I was looking for the closed-form solution for C(s) [in your terms]. While I believe that this problem is theoretically solvable, believe that is no solution except using numeric methods, yet. – Hypermetropia

t(s) depends on an inverse elliptic function, which you must anyway compute numerically. This is easy to achieve as an ODE. – Thibaud

I don't know what is ODE. Don't believe that you MUST solve this numerically. I did use a numeric solution, in fact. Yet, I believe that a closed-form solution is possible. We would get to 4th degree equations, and those are solvable. I personally did not figure out the solution. I believe that someone one day will come with a clever substitution and publish the solution. I will be checking out SIGGRAPH topics from time to time. – Hypermetropia

@IgorShmukler: sorry, I was wrong. For a quadratic Bezier, the arc length has an analytical expression involving an hyperbolic function. This expression cannot be inverted analytically and will never. Just like quintic polynomials will never be solved by radicals. – Thibaud