Table of Contents
1. Cubic parametric spline
1.1. Definition
1.2 Expressions
1.3 Algorithm
1.4 Code implementation
1.4.1 Read value point
1.4.2 Draw type value points
1.4.3 Drawing of cubic parametric splines
1.4.4 Main function call
2. Cardinal Curve
2.1 Hermite basis matrix
2.2 Cardinal Curve
2.3 Cardinal algorithm
2.4 Code implementation
2.4.1 Read value point
2.4.2 Draw type value points
2.4.3 Calculation and drawing
1. Cubic parametric spline
The only disadvantage of cubic splines is their lack of geometric invariance. That is, when the type value point undergoes geometric transformation, the parameter increment cannot be guaranteed. Therefore, a cubic parametric spline with chord length as parameter was proposed.
1.1. Definition
It is known that n type-valued points Pi(xi, yi), i = 1, 2,…, n and adjacent type-valued points do not overlap; if p(t) satisfies the following conditions:
- The type-valued point Pi is on the function p(t).
- p(t) is second-order continuously differentiable over the entire interval [P1, Pn].
- On each subinterval [Pi, Pi + 1], i = 1, 2,…, n-1, the piecewise function p(t) is a cubic polynomial of parameter t. Then the function is said to be a cubic parametric spline function passing through the value point, and the curve composed of the cubic parametric spline function is called a cubic parametric spline curve.
in,
The chord length Li of the subinterval is:
1.2 Expression
The pi(t) of paragraph i is expressed as: 
In general, a cubic parametric spline is a curve expression that parameterizes each xyz component based on a cubic spline, as follows:
- The three coordinate components x, y, and z are cubic spline functions of parameter t respectively.
- Parameterize the type value points.
- The three coordinate components are processed separately.
1.3 algorithm
(1) Read in the coordinates of n type value points.
(2) Determine the boundary conditions of the cubic parametric spline according to the actual situation.
(3) Calculate the coefficient of the curve and express it as a function of the second derivative of the type value point.
(4) Use the chasing method to solve the three bending moment equations along the x direction and y direction respectively.
(5) Substitute the solved coefficients into the expressions of the x component and y component of the cubic parametric spline function to construct a cubic parametric spline curve.
(6) Loop through each node. In each sub-interval, the parameter t is incremented according to the chord length, and a cubic parametric spline is drawn based on each set of (x, y) coordinates.
1.4 Code Implementation
1.4.1 Read value point
First we need to expand the following two-dimensional point classes p2.h and p2.cpp (mainly as follows):
class CP2
{
public:
CP2();
virtual ~CP2();
CP2(double x, double y);
friend CP2 operator + (const CP2 & amp;p0, const CP2 & amp;p1);//Operator overloading
friend CP2 operator -(const CP2 & amp;p0, const CP2 & amp;p1);
friend CP2 operator -(double scalar, const CP2 & amp;p);
friend CP2 operator -(const CP2 & amp;p, double scalar);
friend CP2 operator *(const CP2 & amp;p, double scalar);
friend CP2 operator *(double scalar, const CP2 & amp;p);
friend CP2 operator /(const CP2 & amp;p0, CP2 & amp;p1);
friend CP2 operator /(const CP2 & amp;p, double scalar);
public:
double x;//X coordinate of the endpoint of the straight line segment
double y;//y coordinate of the endpoint of the straight line segment
};
CP2::CP2()
{
}
CP2::~CP2()
{
}
CP2::CP2(double x, double y)
{
this->x = x;
this->y = y;
}
CP2 operator + (const CP2 & amp;p0, const CP2 & amp;p1)//Sum of objects
{
CP2 result;
result.x = p0.x + p1.x;
result.y = p0.y + p1.y;
return result;
}
CP2 operator -(const CP2 & amp;p0, const CP2 & amp;p1)//Difference of objects
{
CP2 result;
result.x = p0.x - p1.x;
result.y = p0.y - p1.y;
return result;
}
CP2 operator -(double scalar, const CP2 & amp;p)//Difference between constant and object
{
CP2 result;
result.x = scalar - p.x;
result.y = scalar - p.y;
return result;
}
CP2 operator -(const CP2 & amp;p, double scalar)//Difference between object and constant
{
CP2 result;
result.x = p.x - scalar;
result.y = p.y - scalar;
return result;
}
CP2 operator *(const CP2 & amp;p, double scalar)//Product of object and constant
{
return CP2(p.x * scalar, p.y * scalar);
}
CP2 operator *(double scalar, const CP2 & amp;p)//Product of constant and object
{
return CP2(p.x * scalar, p.y * scalar);
}
CP2 operator /(const CP2 & amp;p0, CP2 & amp;p1)//Quotient of the object
{
if(fabs(p1.x)<1e-6)
p1.x = 1.0;
if(fabs(p1.y)<1e-6)
p1.y = 1.0;
CP2 result;
result.x = p0.x / p1.x;
result.y = p0.y / p1.y;
return result;
}
CP2 operator /(const CP2 & amp;p, double scalar)//Division by object number
{
if(fabs(scalar)<1e-6)
scalar = 1.0;
CP2 result;
result.x = p.x / scalar;
result.y = p.y / scalar;
return result;
}
Then define the reading function in the main function:
// CGeometricfiguretestViewmessage handlers
void CGeometricfiguretestView::ReadPoint()
{
P[1].x = -400, P[1].y = -150;
P[2].x = -100, P[2].y = 0;
P[3].x = -300, P[3].y = 100;
P[4].x = 0, P[4].y = 200;
P[5].x = 300, P[5].y = 100;
P[6].x = 100, P[6].y = 0;
P[7].x = 400, P[7].y = -150;
}
1.4.2 Drawing type value points
//Draw type value point
void CGeometricfiguretestView::DrawDataPoint(CDC* pDC)
{
CBrush NewBrush, *OldBrush;
NewBrush.CreateSolidBrush(RGB(0, 0, 0));
OldBrush = pDC->SelectObject( & amp;NewBrush);
for(int i = 1; i < 8; i + + )
pDC->Ellipse(ROUND(P[i].x - 5), ROUND(P[i].y - 5), ROUND(P[i].x + 5), ROUND(P[i].y + 5));
pDC->SelectObject(OldBrush);
}
1.4.3 Cubic Parametric Spline Curve Drawing
Cubic parameter spline calculation and drawing, see the notes for the specific process (the calculation and derivation are the same as those for cubic spline curves):
//Cubic parametric spline curve
void CGeometricfiguretestView::DrawCubicSpline(CDC* pDC)
{
int n = 7;
const int dim = 8;//Dimension of two-dimensional array
double b1 = 0, bn = 0;//The first derivative of the starting point and the end point
double L[dim];//chord length of parameter spline curve
double lambda[dim], mu[dim];
double l[dim], m[dim], u[dim];
CP2 c1, cn;//direction cosine of starting point and end point
CP2 D[dim];
CP2 M[dim], K[dim];//Chasing method transition matrix
CP2 B1[dim], B2[dim], B3[dim], B4[dim];//Coefficients of the function
CP2 delt[dim];
for(int i=1; i<n;i + + )//Calculate the chord length
{
delt[i]=P[i + 1]-P[i];
L[i]=sqrt(delt[i].x*delt[i].x + delt[i].y*delt[i].y);
}
//Projection of boundary conditions
c1.x = b1*cos(delt[1].x/L[1]);//Starting point
c1.y = b1*cos(delt[1].y/L[1]);
cn.x = bn*cos(delt[n-1].x/L[n-1]);//end point
cn.y = bn*cos(delt[n-1].y/L[n-1]);
for(int i = 2; i < n; i + + )
{
lambda[i] = L[i-1]/(L[i-1] + L[i]);//Calculate λ
mu[i]=L[i]/(L[i-1] + L[i]);//Calculate μ
D[i]=6/(L[i-1] + L[i])*((P[i + 1]-P[i])/L[i]-(P[i]-P[i -1])/L[i-1]);//Calculate D
}
D[1]=6*((P[2]-P[1])/L[1]-b1)/L[1];//D[1] at the clamping end
D[n]=6*(bn-(P[n]-P[n-1])/L[n-1])/L[n-1];//D[n] of the clamping end
mu[1]=1;//μ[1] of the clamping end,
lambda[n]=1;//λ[n] of the clamping end
//Catching method to solve three bending moment equations
l[1]=2;
u[1]=mu[1]/l[1];
for(int i=2; i <= n; i + + )
{
m[i]=lambda[i];
l[i]=2-m[i]*u[i-1];
u[i]=mu[i]/l[i];
}
K[1] = D[1]/l[1];//Solution LK=D
for(int i = 2; i <= n;i + + )
{
K[i]=(D[i]-m[i]*K[i-1])/l[i];
}
M[n] = K[n];//Solution UM=K
for(int i = n-1; i >= 1;i--)
{
M[i]=K[i]-u[i]*M[i + 1];
}
//Calculate the coefficients of the cubic spline function
for(int i = 1; i < n; i + + )
{
B1[i]=P[i];
B2[i]=(P[i + 1]-P[i])/L[i]-L[i]*(M[i]/3 + M[i + 1]/6);
B3[i]=M[i]/2;
B4[i]=(M[i + 1]-M[i])/(6*L[i]);
}
CPen pen(PS_SOLID,3,RGB(0,0,0));
pDC->SelectObject( & amp;pen);
pDC->MoveTo(ROUND(P[1].x), ROUND(P[1].y));
double tStep = 0.5;//step size
CP2p;
for(int i = 1; i < n; i + + )//Loop through each node
{
for(double t=0; t <= L[i]; t + = tStep)
{
p =B1[i] + B2[i]*t + B3[i]*t*t + B4[i]*t*t*t;
pDC->LineTo(ROUND(p.x), ROUND(p.y));//Draw parametric spline curve
}
}
1.4.4 Main function call
void CGeometricfiguretestView::OnDraw(CDC* pDC)
{
CTestDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
if (!pDoc)
return;
// TODO: add draw code for native data here
CRect rect;//Define the client area rectangle
GetClientRect( & amp;rect);//Get information about the rectangle in the client area
pDC->SetMapMode(MM_ANISOTROPIC);//Customized two-dimensional coordinate system
pDC->SetWindowExt(rect.Width(), rect.Height());//Set the window range
pDC->SetViewportExt(rect.Width(), -rect.Height());//Set the viewport range, and the x-axis is positive horizontally to the right, and the y-axis is positive vertically upward
pDC->SetViewportOrg(rect.Width() / 2, rect.Height() / 2);//Set the center of the client area as the origin of the two-dimensional coordinate system
rect.OffsetRect(-rect.Width() / 2, -rect.Height() / 2);//The rect rectangle coincides with the client area
ReadPoint();
DrawDataPoint(pDC);
DrawCubicParameterSpline(pDC);
}
Compile and run, as shown below:
We can also change the boundary constraints:
double b1 = 1, bn = 1;//The first derivative of the starting point and the end point
2. Cardinal Curve
2.1 Hermite basis matrix
Each Hermite curve is defined by the coordinates and derivatives of the two endpoints.
Assume that the equation of the Hermite spline can be written in matrix form, as follows:
Substituting the boundary conditions p(0)=yi and p(1)=yi + 1 into the above formula, we get:
And the first-order derivative of the above formula is:
Putting the boundary conditions p’(0)=Ri and p’(1)=Ri + 1 into the above formula, we get:
The matrix representation of Hermite boundary conditions is:
This equation is solved for the polynomial coefficients:
In the formula, Mh is called the Hermite basis matrix, which is the inverse of the boundary constraint matrix matrix. It should be noted that the uniform parameter interpolation cubic spline equation used by Ferguson in 1963 for aircraft design is Mh.
Therefore, we can get:
Finally you can get:
2.2 Cardinal Curve
The Cardinal spline is a deformation of the Hermite spline, and the derivative can be calculated only from the coordinates of adjacent type value points. Let the four adjacent type value points be recorded as Pi-1, Pi, Pi + 1, and Pi + 2 respectively. The Cardinal spline interpolation method stipulates that the boundary conditions of the interpolation polynomial between two types of value points Pi, Pi + 1 are:
In the formula, u is an adjustable parameter, called the tension parameter, which can control the tightness between the value points of the Cardinal spline curve.
Remember s=(1-u)/2, use a method similar to the Hermite curve spline, and substitute the Cardinal boundary condition into the parametric equation below,
You can get the matrix expression:
In the formula, Mc is called the Cardinal matrix.
A Cardinal spline curve is completely given by four continuous type value points. The two middle type value points are the end points of the curve segment, and the other two type value points are used to assist in calculating the end point derivatives. As long as the coordinate values of a set of type-value points are given, a Cardinal spline can be drawn piecewise.
2.3 Cardinal Algorithm
(1) Read in n type value point coordinates Pi, i=1,2,…n.
(2) Add a virtual point P0 and Pn + 1 at both ends.
(3) Set the tension parameter u value and calculate the Cardinal basis matrix Mc.
(4) Algorithm for calculating Mc matrix and boundary condition matrix. Every 4 vertices constitute a condition matrix.
(5) Calculate the value of the interpolation parameter t in each curve according to the number of type value points, and use straight line segments to connect the corresponding points p(x(t), y(t)) to draw the curve.
2.4 Code Implementation
2.4.1 Read-type value point
void CGeometricfiguretestView::ReadPoint()
{
P[0].x =-500, P[0].y =0;
P[1].x =-450, P[1].y =-100;
P[2].x =-350, P[2].y =-250;
P[3].x =-250, P[3].y = 200;
P[4].x =-150, P[4].y = 100;
P[5].x =-50, P[5].y = 150;
P[6].x = 40, P[6].y = 10;
P[7].x = 100, P[7].y = -60;
P[8].x = 150, P[8].y = 80;
P[9].x = 200, P[9].y = 70;
P[10].x =230,P[10].y =-10;
P[11].x =300,P[11].y =-200;
P[12].x =400,P[12].y = 150;
P[13].x =520,P[13].y = 250;
n = 13;
}
2.4.2 Drawing type value points
void CGeometricfiguretestView::DrawDataPoint(CDC* pDC)//Drawing value point
{
CBrush NewBrush, *pOldBrush;
NewBrush.CreateSolidBrush(RGB(0, 0, 0));
pOldBrush = pDC->SelectObject( & amp;NewBrush);
for( int i = 1; i < 13; i + + )
pDC->Ellipse(ROUND(P[i].x - 5), ROUND(P[i].y - 5), ROUND(P[i].x + 5), ROUND(P[i].y + 5));
pDC->SelectObject(pOldBrush);
pDC->Ellipse(ROUND(P[0].x - 5), ROUND(P[0].y - 5), ROUND(P[0].x + 5), ROUND(P[0].y + 5));
pDC->Ellipse(ROUND(P[13].x - 5), ROUND(P[13].y - 5), ROUND(P[13].x + 5), ROUND(P[13].y + 5));
}
2.4.3 Calculation and Drawing
void CGeometricfiguretestView::DrawCardinalCurve(CDC* pDC)//Cardinal Curve
{
double M[4][4];
double u=0.2;
double s=(1-u)/2;
CP2p;
M[0][0]=-s,M[0][1]=2-s,M[0][2]= s-2,M[0][3]= s;//Mc matrix
M[1][0]=2*s,M[1][1]=s-3,M[1][2]=3-2*s,M[1][3]=-s;
M[2][0]=-s ,M[2][1]= 0 ,M[2][2]= s ,M[2][3]= 0;
M[3][0]= 0 ,M[3][1]= 1 ,M[3][2]= 0 ,M[3][3]= 0;
double t3,t2,t1,t0;
CPen pen(PS_SOLID,2,RGB(0,0,0));
CPen* pOldPen=pDC->SelectObject( & amp;pen);
pDC->MoveTo(ROUND(P[1].x),ROUND(P[1].y));
for(int i=0; i< n-2; i + + )
{
Point[0]=P[i],Point[1]=P[i + 1],Point[2]=P[i + 2],Point[3]=P[i + 3];
MultiplyMatrix(M, P, i);
double tStep = 0.001;
for(double t=0.0; t<1; t + = tStep)
{
t3 = t*t*t; t2 = t*t; t1 = t; t0 = 1;
p = t3*Point[0] + t2*Point[1] + t1*Point[2] + t0*Point[3];
pDC->LineTo(ROUND(p.x),ROUND(p.y));
}
}
pDC->SelectObject(pOldPen);
}
void CGeometricfiguretestView::MultiplyMatrix(double M0[][4], CP2 P0[4], int n)
{
for(int i=0;i<4;i + + )
Point[i] = M0[i][0]*P0[n] + M0[i][1]*P0[n + 1] + M0[i][2]*P0[n + 2] + M0[ i][3]*P0[n + 3];
}
The MultiplyMatrix function is defined to realize the multiplication calculation of the 4×4 Mc matrix and the 1×4 boundary constraint matrix. Compile and run as follows:
We can also try changing the tension parameters:
double u=1.0;
The result is as follows:
The above is the code implemented by MFC. For Qt implementation, see: Qt Implementation of Cubic Spline Cardinal Curve
If you want to implement a 3D version of the Cardinal curve, please refer to: osg implements cubic spline Cardinal curve
This article is reproduced from: Computational Geometry 03_Cubic Parametric Splines and Cardinal Curves