Close

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.

agpcooperagp.cooper • 02/12/2024 at 02:11•0 Comments

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;
Previous log

A Non-Suleiman IP Version

01/31/2024 at 13:32 • 0 comments
Next log

Code for Cubic Bezier Version

02/18/2024 at 11:00 • 0 comments

Discussions

Yes, delete it Cancel

Going up?

About Us Contact Hackaday.io Give Feedback Terms of Use Privacy Policy Hackaday API

© 2024 Hackaday