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.
OK I found analytic formula for 2D quadratic bezier curve in here:
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:
x0,y0,...t=get_t(wanted_arclength) freelyIn 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;
}
//---------------------------------------------------------------------------