0%
0%

# Fitting a Quadratic Bezier Curve to Point Data

Trying to reverse engineer a round bottom sail boat.

Similar projects worth following
137 views
A Bezier curve is a family of curves commonly used in CAD programs.
Here we are dealing with simplest (3 point quadratic version).
The (parametric) equations can be described as:
X(t)=Ax*(1-t)^2+Bx*(1-t)*t+Cx*t^2
Y(t)=Ay*(1-t)^2+By*(1-t)*t+Cy*t^2
The points "A", "B" and "C" are the first, control and last points of the Bezier curve.
The problem is to match the "t" parameter to a set of digitised points, then to calculate the control point "B" that minimises the fit error.
This is a non-linear problem so direct calculation is not possible (in this case).
I have solved this problem using the Newton-Raphson method.
As the problem is non-linear there is no guarantee that all problems will return an optimal solution.

Fitting a Quadratice Bezier Curve to Digitised Data

Here is the problem image:

The blue and black digitised lines are the Gunwale and Sheer (beam) of the sailboat design I want to reverse engineer. The image is a plan of the Lynaes 14 sailboat that I found on the Internet.

I have previously tried to model the beam in Excel using polynomial curve fitting but the results were not satisfactory.

However, fitting quadratic Bezier curves was excellent:

The image has four Bezier curves joined at the max beam.

Note how well the Bezier curves meet at the stern. This is a problem area for polynomial curves.

Here is the polynomial 6th order fit. While it looks okay at this scale, the beam wanders and the curve at the stern is not ideal (i.e. perpendicular):

Here is an image of a Lynaes 14 stern which is near perpendular to the keel:

### BezierFitV2.c

x-csrc - 5.67 kB - 01/31/2024 at 13:36

### Results.png

Portable Network Graphics (PNG) - 8.38 kB - 01/31/2024 at 13:36

Suleiman's paper

Adobe Portable Document Format - 154.95 kB - 01/31/2024 at 07:51

### BezierFit.run

A linux program executable

run - 55.88 kB - 01/31/2024 at 07:51

### ezx_graphics.pdf

Manual for EZX

Adobe Portable Document Format - 92.72 kB - 01/31/2024 at 07:51

• ### Code for Cubic Bezier Version

Cubic Code Results

While the results are good it is not exact like the piece-wise quadratic bezier code.

I quess a fourth order bezier would bit better.

Hers is the cubic match:

Note the bezier cuves are too wide just aft of the maxium beam and too narrow near the stern.

Here are the control points:

Here is the full code:

```// Include libraries
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>

// Include in local directory:
#include "ezxdisp.h"
// Library: libezx.a
// Linker: -L. -lezx -lX11 -lm

// EZX window size
static int width=1200;
static int height=700;
static double scale=0.2;
// EZX display
static ezx_t *disp;

void PlotBezier(double Xi[],double Yi[],double Ti[],int Ni,double *Bx1,double *By1,double *Bx2,double *By2)
{
// Plot Fitted Bezier Curve
double X0,Y0,X1,Y1;
X1=Xi[0];
Y1=Yi[0];
for (int i=1;i<Ni;i++) {
double t=1.0*i/(Ni-1);
double s=(1.0-t);
X0=X1;
Y0=Y1;
X1=s*s*s*Xi[0]+3*s*s*t**Bx1+3*s*t*t**Bx2+t*t*t*Xi[Ni-1];
Y1=s*s*s*Yi[0]+3*s*s*t**By1+3*s*t*t**By2+t*t*t*Yi[Ni-1];
ezx_line_2d(disp,scale*X0,height/2-scale*Y0,scale*X1,height/2-scale*Y1,&ezx_red,1);
}

// Plot Fitted Control Point (B1)
ezx_line_2d(disp,scale**Bx1-5,height/2-scale**By1,scale**Bx1+5,height/2-scale**By1,&ezx_black,1);
ezx_line_2d(disp,scale**Bx1,height/2-scale**By1-5,scale**Bx1,height/2-scale**By1+5,&ezx_black,1);

// Plot Fitted Control Point (B2)
ezx_line_2d(disp,scale**Bx2-5,height/2-scale**By2,scale**Bx2+5,height/2-scale**By2,&ezx_black,1);
ezx_line_2d(disp,scale**Bx2,height/2-scale**By2-5,scale**Bx2,height/2-scale**By2+5,&ezx_black,1);

// Plot Input points
for (int i=0;i<Ni;i++) {
ezx_line_2d(disp,scale*Xi[i]-2,height/2-scale*Yi[i]-2,scale*Xi[i]+2,height/2-scale*Yi[i]+2,&ezx_black,1);
ezx_line_2d(disp,scale*Xi[i]+2,height/2-scale*Yi[i]-2,scale*Xi[i]-2,height/2-scale*Yi[i]+2,&ezx_black,1);
}

return;
}

void FitBezier(double Xi[],double Yi[],double Ti[],int Ni,double *x1,double *y1,double *x2,double *y2)
{
double Bx0=Xi[0];
double By0=Yi[0];
double Bx1=*x1;
double By1=*y1;
double Bx2=*x2;
double By2=*y2;
double Bx3=Xi[Ni-1];
double By3=Yi[Ni-1];
double Ex1,E1x1,E2x1;
double Ey1,E1y1,E2y1;
double Ex2,E1x2,E2x2;
double Ey2,E1y2,E2y2;
double t,e2=0.0;
do {
E1x1=0.0;
E2x1=0.0;
E1y1=0.0;
E2y1=0.0;
E1x2=0.0;
E2x2=0.0;
E1y2=0.0;
E2y2=0.0;
e2=0.0;
// Solve all t
for (int i=1;i<Ni-1;i++) {

double Et=0.0;
t=Ti[i]; // Initial t
// printf("i=%3d t=%12.6lf",i,t);
do {
double f1=(t*t*t*(-Bx0+3*Bx1-3*Bx2+Bx3)+t*t*(3*Bx0-6*Bx1+3*Bx2)+t*(-3*Bx0+3*Bx1)+(Bx0-Xi[i]))*(6*t*t*(-Bx0+3*Bx1-3*Bx2+Bx3)+4*t*(3*Bx0-6*Bx1+3*Bx2)+2*(-3*Bx0+3*Bx1))
+(t*t*t*(-By0+3*By1-3*By2+By3)+t*t*(3*By0-6*By1+3*By2)+t*(-3*By0+3*By1)+(By0-Yi[i]))*(6*t*t*(-By0+3*By1-3*By2+By3)+4*t*(3*By0-6*By1+3*By2)+2*(-3*By0+3*By1));
double f2=(3*t*t*(-Bx0+3*Bx1-3*Bx2+Bx3)+2*t*(3*Bx0-6*Bx1+3*Bx2)+(-3*Bx0+3*Bx1))*(6*t*t*(-Bx0+3*Bx1-3*Bx2+Bx3)+4*t*(3*Bx0-6*Bx1+3*Bx2)+2*(-3*Bx0+3*Bx1))
+(t*t*t*(-Bx0+3*Bx1-3*Bx2+Bx3)+t*t*(3*Bx0-6*Bx1+3*Bx2)+t*(-3*Bx0+3*Bx1)+(Bx0-Xi[i]))*(12*t*(-Bx0+3*Bx1-3*Bx2+Bx3)+4*(3*Bx0-6*Bx1+3*Bx2))
+(3*t*t*(-By0+3*By1-3*By2+By3)+2*t*(3*By0-6*By1+3*By2)+(-3*By0+3*By1))*(6*t*t*(-By0+3*By1-3*By2+By3)+4*t*(3*By0-6*By1+3*By2)+2*(-3*By0+3*By1))
+(t*t*t*(-By0+3*By1-3*By2+By3)+t*t*(3*By0-6*By1+3*By2)+t*(-3*By0+3*By1)+(By0-Yi[i]))*(12*t*(-By0+3*By1-3*By2+By3)+4*(3*By0-6*By1+3*By2));
Et=f1/f2;
t=t-Et;
} while (Et*Et>1e-8); // Error ~0.001
Ti[i]=t;
e2=e2+(Bx0*(1-t)*(1-t)*(1-t)+3*Bx1*(1-t)*(1-t)*t+3*Bx2*(1-t)*t*t+Bx3*t*t*t-Xi[i])*(Bx0*(1-t)*(1-t)*(1-t)+3*Bx1*(1-t)*(1-t)*t+3*Bx2*(1-t)*t*t+Bx3*t*t*t-Xi[i])
+(By0*(1-t)*(1-t)*(1-t)+3*By1*(1-t)*(1-t)*t+3*By2*(1-t)*t*t+By3*t*t*t-Yi[i])*(By0*(1-t)*(1-t)*(1-t)+3*By1*(1-t)*(1-t)*t+3*By2*(1-t)*t*t+By3*t*t*t-Yi[i]);
// printf(" -> t=%12.6lf Et*Et=%12.6lf",t,Et*Et);getchar();

// WRT: Bx1
E1x1=E1x1+6*t*(t*t*t*(-Bx0+3*Bx1-3*Bx2+Bx3)+t*t*(3*Bx0-6*Bx1+3*Bx2)+t*(-3*Bx0+3*Bx1)+(Bx0-Xi[i]))*(t*t-2*t+1);
E2x1=E2x1+18*t*t*(t*t-2*t+1)*(t*t-2...```

• ### New and Improved Version

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))...

• ### A Non-Suleiman IP Version

Non-Suleiman Version

After my last post I realised I did not need Suleiman's code. The Newton-Rapson method could solve the "t" values directly, but it would likely need more iterations.

Here are the results:

So the results are indentical which means my derivation (based on first and second derivates) is but a refactorisation of Suleiman's.

BezierFitV2.c is in the files directory.

• ### Coding a Solution

Mapping Digitised Points to a Trial Quadratic Bezier Curve

While the quadratic Bezier curve equation is quite simple, it evaded me how to convert a digitised point to the Bezier curve parameter "t".

Upon research I came upon a paper by Nouri A. Suleiman (2013), "Least Squares Data Fitting with Quadratic Bezier Curves".

Although I did not follow the Least Squares route, the paper presents a method for mapping a digitised point  (i.e. finding the t parameter) to a trial quadratic bezier curve.

For each digitised point Xi[i] and Yi[i], and the trial Bezier curve:

A*(1-t)^2+2*B*(1-t)*t+C*t*t

Calculate the following (Suleiman) values:

d0=(Ax-Xi[i])*(Ax-Xi[i])+(Ay-Yi[i])*(Ay-Yi[i]);

d1=(*Bx-Ax)*(Ax-Xi[i])+(*By-Ay)*(Ay-Yi[i]);

d2=(Ax-2*(*Bx)+Cx)*(Ax-Xi[i])+2*(*Bx-Ax)*(*Bx-Ax)+(Ay-2*(*By)+Cy)*(Ay-Yi[i])

+2*(*By-Ay)*(*By-Ay);

d3=(Ax-2*(*Bx)+Cx)*(*Bx-Ax)+(Ay-2*(*By)+Cy)*(*By-Ay);

d4=Ax*(Ax-2*(*Bx)+Cx)-2*(*Bx)*(Ax-2*(*Bx)+Cx)+Cx*(Ax-2*(*Bx)+Cx)

+Ay*(Ay-2*(*By)+Cy)-2*(*By)*(Ay-2*(*By)+Cy)+Cy*(Ay-2*(*By)+Cy);

Then calculate the (Suleiman) parameter mapping error:

e2=d4*t*t*t*t+4*d3*t*t*t+2*d2*t*t+4*d1*t+d0;

Then solve the minimisation problem of e2 for t. We can use the Newton-Raphson method:

t=t-e2'/e2"

Where:

e2'=4*(d4*t*t*t+3*d3*t*t+d2*t+d1);
e2"=4*(3*d4*t*t+6*d3*t+d2);

After finding the trial Bezier curve parameter t for each digitized point, we need to solve the minimisation problem of e2 for Bx of the control point (i.e calculate the first and second derivatives for the sum of e2):

Sum[e2]=Sum[(Ax*s*s+2*Bx*s*t+Cx*t*t-Xi[i])^2], for each i
Sum[e2']=Sum[4*s*t*(Ax*s*s+2*Bx*s*t+Cx*t*t-Xi[i])], for each i
=Sum[4*(Ax*s*s*s*s+2*Bx*s*s*t*t+Cx*s*t*t*t-Xi[i]*s*t)], for each i
Sum[e2"]=Sum[8*s*s*t*t], for each i

Then apply Newton-Raphson:
Bx=Bx-e2'/e2''

Repeat for By.

In the files directory is some C code implimenting the psudo code above.

Share

## Similar Projects

Project Owner Contributor

### DXF Laser/Mill Box Maker

agp.cooper

Project Owner Contributor

### TinyPong

Saul

Project Owner Contributor

### Mechatronic Ears

deʃhipu

Project Owner Contributor

bret.wiebe

# Does this project spark your interest?

Become a member to follow this project and never miss any updates