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.

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://stackoverflow.com/a/11857788/2521214
    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;
    }
//---------------------------------------------------------------------------

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.

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.