Оптимизируйте трассировку квадратичных кривых с помощью численных методов

ОК, я нашел аналитическую формулу для двумерной квадратичной кривой Безье здесь:

• Вычислите длину отрезка квадратичного безье

Итак, идея заключается в простом бинарном поиске параметра t до тех пор, пока аналитически полученная длина дуги не совпадет с желаемой длиной...

Код С++:

//---------------------------------------------------------------------------
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);
    }
//---------------------------------------------------------------------------

Применение:

1. установить контрольные точки x0,y0,...
2. тогда вы можете свободно использовать t=get_t(wanted_arclength)

Если вы хотите использовать get_t_naive и или get_xy, вы должны сначала вызвать set_coef

Если вы хотите настроить скорость / точность, вы можете играть с целевой точностью binsearch, установленной в настоящее время на1e-10

Здесь оптимизированная (объединенная get_l,get_t функция) версия:

//---------------------------------------------------------------------------
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;
    }
//---------------------------------------------------------------------------

На данный момент я придумал следующее:

    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.info('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);
    }

Я думаю только о том, что моя длина кривой Безье работает не так, как я ожидал.

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)
  );
}

Я считаю, что полная длина кривой верна, но частичная длина, вычисляемая для таблицы поиска, неверна.

Если я прав, вам нужны точки в равноотстоящих точках с точки зрения криволинейной абсциссы (а не с точки зрения постоянного евклидова расстояния, что было бы совсем другой проблемой).

Вычисление криволинейной абсциссы s как функции параметра кривой t действительно вариант, но это приводит к решению уравнения s(t) = Sk/n для целого числа k, где S — общая длина (или s(t) = kd, если наложен шаг). Это неудобно, потому что s(t) недоступна как простая функция и является трансцендентной.

Лучшим методом является решение дифференциального уравнения

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

используя предпочитаемый вами решатель ODE (RK4). Это позволяет вам наложить фиксированный шаг на s и эффективно с точки зрения вычислений.