# New and Improved Version

A project log for Fitting a Quadratic Bezier Curve to Point Data

Trying to reverse engineer a round bottom sail boat.

Code Rebuild

I rebuilt the code for a different staing point. Instead of:

B(t)=(1-t)^2*B0+2*t*(1-t)*B1+t^2*B2

I refacted the equation to:
B(t)=(B0-2*B1+B2)*t^2+2*(B1-B0)*t+B0

The error squared for each "t" is:

e^2=((Bx0-2*Bx1+Bx2)*t^2+2*(Bx1-Bx0)*t+Bx0-Xi)^2

+((By0-2*By1+By2)*t^2+2*(By1-By0)*t+By0-Yi)^2

This equation (e2) calculates the square of distance between the interpolated point (Xi,Yi) from the closest point on the "trial" Bezier curve.

The updated procedure is:

```void FitBezier(double Xi[],double Yi[],double Ti[],int Ni,double *x1,double *y1)
{
double Bx0=Xi[0];
double By0=Yi[0];
double Bx1=*x1;
double By1=*y1;
double Bx2=Xi[Ni-1];
double By2=Yi[Ni-1];
double Ex,E1x1,E2x1;
double Ey,E1y1,E2y1;
double t;
do {
E1x1=0.0;
E2x1=0.0;
E1y1=0.0;
E2y1=0.0;
// Solve all t
for (int i=1;i<Ni-1;i++) {
// Solve t
double b4=(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)+(By0-2*By1+By2)*(By0-2*By1+By2);
double b3=(Bx0-2*Bx1+Bx2)*(Bx1-Bx0)+(By0-2*By1+By2)*(By1-By0);
double b2=2*(Bx1-Bx0)*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(Bx0-Xi[i])+2*(By1-By0)*(By1-By0)+(By0-2*By1+By2)*(By0-Yi[i]);
double b1=(Bx1-Bx0)*(Bx0-Xi[i])+(By1-By0)*(By0-Yi[i]);
// double b0=(Bx0-Xi[i])*(Bx0-Xi[i])+(By0-Yi[i])*(By0-Yi[i]);
double Et=0.0;
t=Ti[i]; // Initial t
do {
double Et1=4*t*t*t*b4+12*t*t*b3+4*t*b2+4*b1;
double Et2=12*t*t*b4+24*t*b3+4*b2;
Et=Et1/Et2;
t=t-Et;
} while (Et*Et>1e-6); // Error ~0.001
Ti[i]=t;
E1x1=E1x1+(-4*(Bx0-2*Bx1+Bx2))*t*t*t*t+4*(-2*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2))*t*t*t+4*(2*(Bx1-Bx0)-(Bx0-Xi[i]))*t*t+4*(Bx0-Xi[i])*t;
E2x1=E2x1+8*t*t*t*t-16*t*t*t+8*t*t;
E1y1=E1y1+(-4*(By0-2*By1+By2))*t*t*t*t+4*(-2*(By1-By0)+(By0-2*By1+By2))*t*t*t+4*(2*(By1-By0)-(By0-Yi[i]))*t*t+4*(By0-Yi[i])*t;
E2y1=E2y1+8*t*t*t*t-16*t*t*t+8*t*t;
}
Ex=E1x1/E2x1;
Ey=E1y1/E2y1;
Bx1=Bx1-Ex;
By1=By1-Ey;
} while (Ex*Ex+Ey*Ey>1e-6); // Error ~0.001
*x1=Bx1;
*y1=By1;

return;
}```

Derivation

Only interesting if you want to do it yourself!

Refactorise B(t):

B(t)=(1-t)^2*B0+2*t*(1-t)*B1+t^2*B2

=(B0-2*B1+B2)*t^2+2*(B1-B0)*t+B0

Solve Derivatives for e2 with respect to "t":

e2=((Bx0-2*Bx1+Bx2)*t^2+2*(Bx1-Bx0)*t+Bx0-Xi)^2+((By0-2*By1+By2)*t^2+2*(By1-By0)*t+By0-Yi)^2
=((Bx0-2*Bx1+Bx2)*t^2+2*(Bx1-Bx0)*t+Bx0-Xi)^2
+((By0-2*By1+By2)*t^2+2*(By1-By0)*t+By0-Yi)^2
=(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)*t^4
+4*(Bx1-Bx0)*(Bx1-Bx0)*t^2
+  (Bx0-Xi)*(Bx0-Xi)
+2*(Bx0-2*Bx1+Bx2)*t^2*2*(Bx1-Bx0)*t
+2*(Bx0-2*Bx1+Bx2)*t^2  *(Bx0-Xi)
+2*2*(Bx1-Bx0)*t*(Bx0-Xi)
+...
=t^4  *(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)
+t^2*4*(Bx1-Bx0)*(Bx1-Bx0)
+      (Bx0-Xi)*(Bx0-Xi)
+t^3*4*(Bx0-2*Bx1+Bx2)*(Bx1-Bx0)
+t^2*2*(Bx0-2*Bx1+Bx2)*(Bx0-Xi)
+t  *4*(Bx1-Bx0)*(Bx0-Xi)
+...
=t^4  *(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)
+t^3*4*(Bx0-2*Bx1+Bx2)*(Bx1-Bx0)
+t^2*2*(2*(Bx1-Bx0)*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(Bx0-Xi))
+t  *4*(Bx1-Bx0)*(Bx0-Xi)
+      (Bx0-Xi)*(Bx0-Xi)
+...
=t^4  *(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)
+t^3*4*(Bx0-2*Bx1+Bx2)*(Bx1-Bx0)
+t^2*2*(2*(Bx1-Bx0)*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(Bx0-Xi))
+t  *4*(Bx1-Bx0)*(Bx0-Xi)
+      (Bx0-Xi)*(Bx0-Xi)
+t^4  *(By0-2*By1+By2)*(By0-2*By1+By2)
+t^3*4*(By0-2*By1+By2)*(By1-By0)
+t^2*2*(2*(By1-By0)*(By1-By0)+(By0-2*By1+By2)*(By0-Yi))
+t  *4*(By1-By0)*(By0-Yi)
+      (By0-Yi)*(By0-Yi)

=t^4  *((Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)+(By0-2*By1+By2)*(By0-2*By1+By2))
+t^3*4*((Bx0-2*Bx1+Bx2)*(Bx1-Bx0)+(By0-2*By1+By2)*(By1-By0))
+t^2*2*(2*(Bx1-Bx0)*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(Bx0-Xi)+2*(By1-By0)*(By1-By0)+(By0-2*By1+By2)*(By0-Yi))
+t  *4*((Bx1-Bx0)*(Bx0-Xi)+(By1-By0)*(By0-Yi))
+      ((Bx0-Xi)*(Bx0-Xi)+(By0-Yi)*(By0-Yi))

b4=(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)+(By0-2*By1+By2)*(By0-2*By1+By2)
b3=(Bx0-2*Bx1+Bx2)*(Bx1-Bx0)+(By0-2*By1+By2)*(By1-By0)
b2=2*(Bx1-Bx0)*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(Bx0-Xi)+2*(By1-By0)*(By1-By0)+(By0-2*By1+By2)*(By0-Yi)
b1=(Bx1-Bx0)*(Bx0-Xi)+(By1-By0)*(By0-Yi)
b0=(Bx0-Xi)*(Bx0-Xi)+(By0-Yi)*(By0-Yi)

Summary for Derivatives for e(t) with respect for "t":

e2(t)=t*t*t*t*b4+4*t*t*t*b3+2*t*t*b2+4*t*b1+b0
e2'(t)=4*t*t*t*b4+12*t*t*b3+4*t*b2+4*b1
e2"(t)=12*t*t*b4+24*t*b3+4*b2
t=t-e2'(t)/e2"(t)

Solve Derivatives with respect to "B1":

b4(Bx1)=(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)+(By0-2*By1+By2)*(By0-2*By1+By2)
b3(Bx1)=(Bx0-2*Bx1+Bx2)*(Bx1-Bx0)+(By0-2*By1+By2)*(By1-By0)
b2(Bx1)=2*(Bx1-Bx0)*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(Bx0-Xi)+2*(By1-By0)*(By1-By0)+(By0-2*By1+By2)*(By0-Yi)
b1(Bx1)=(Bx1-Bx0)*(Bx0-Xi)+(By1-By0)*(By0-Yi)
b0(Bx1)=(Bx0-Xi)*(Bx0-Xi)+(By0-Yi)*(By0-Yi)

b4(Bx1)=(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)+(By0-2*By1+By2)*(By0-2*By1+By2)
b4'(Bx1)=2*(Bx0-2*Bx1+Bx2)*(Bx0-2*Bx1+Bx2)'
b4'(Bx1)=2*(Bx0-2*Bx1+Bx2)*(-2)
b4'(Bx1)=-4*(Bx0-2*Bx1+Bx2)
b4"(Bx1)=-4*(-2)
b4"(Bx1)=8

b3(Bx1)=(Bx0-2*Bx1+Bx2)*(Bx1-Bx0)+(By0-2*By1+By2)*(By1-By0)
b3'(Bx1)=(Bx0-2*Bx1+Bx2)'*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(Bx1-Bx0)'
b3'(Bx1)=(-2)*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(1)
b3'(Bx1)=-2*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)
b3"(Bx1)=-2*(1)+(-2)
b3"(Bx1)=-4

b2(Bx1)=2*(Bx1-Bx0)*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)*(Bx0-Xi)+2*(By1-By0)*(By1-By0)+(By0-2*By1+By2)*(By0-Yi)
b2'(Bx1)=2*(Bx1-Bx0)'*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)'*(Bx0-Xi)+2*(Bx1-Bx0)*(Bx1-Bx0)'+(Bx0-2*Bx1+Bx2)*(Bx0-Xi)'
b2'(Bx1)=2*(1)*(Bx1-Bx0)+(-2)*(Bx0-Xi)+2*(Bx1-Bx0)*(1)+(Bx0-2*Bx1+Bx2)*(0)
b2'(Bx1)=2*(Bx1-Bx0)-2*(Bx0-Xi)+2*(Bx1-Bx0)
b2'(Bx1)=4*(Bx1-Bx0)-2*(Bx0-Xi)
b2"(Bx1)=4*(1)
b2"(Bx1)=4

b1(Bx1)=(Bx1-Bx0)*(Bx0-Xi)+(By1-By0)*(By0-Yi)
b1'(Bx1)=(1)*(Bx0-Xi)
b1'(Bx1)=(Bx0-Xi)
b1"(Bx1)=0

b0(Bx1)=(Bx0-Xi)*(Bx0-Xi)+(By0-Yi)*(By0-Yi)
b0'(Bx1)=0
b0"(Bx1)=0

b4'(Bx1)=-4*(Bx0-2*Bx1+Bx2)
b4"(Bx1)=8
b3'(Bx1)=-2*(Bx1-Bx0)+(Bx0-2*Bx1+Bx2)
b3"(Bx1)=-4
b2'(Bx1)=4*(Bx1-Bx0)-2*(Bx0-Xi)
b2"(Bx1)=4
b1'(Bx1)=(Bx0-Xi)
b1"(Bx1)=0
b0'(Bx1)=0
b0"(Bx1)=0

Summary of derivatives with respect to "B1":

b4'(By1)=-4*(By0-2*By1+By2)
b4"(By1)=8
b3'(By1)=-2*(By1-By0)+(By0-2*By1+By2)
b3"(By1)=-4
b2'(By1)=4*(By1-By0)-2*(By0-Yi)
b2"(By1)=4
b1'(By1)=(By0-Yi)
b1"(By1)=0
b0'(By1)=0
b0"(By1)=0

E1x1=Sum(t*t*t*t*b4'(Bx1)+4*t*t*t*b3'(Bx1)+2*t*t*b2'(Bx1)+4*t*b1'(Bx1));
E2x1=Sum(t*t*t*t*b4"(Bx1)+4*t*t*t*b3"(Bx1)+2*t*t*b2"(Bx1));
Bx1=Bx1-E1x1/E2x1;

E1y1=Sum(t*t*t*t*b4'(By1)+4*t*t*t*b3'(By1)+2*t*t*b2'(By1)+4*t*b1'(By1));
E2y1=Sum(t*t*t*t*b4"(By1)+4*t*t*t*b3"(By1)+2*t*t*b2"(By1));
By1=By1-E1y1/E2y1;