1. CUBIC SPLINE:      Approach - 3rd degree polynomial

PseudoCode CubicSpline

Y=F(X)=A+B(X-X0)1+C(X-X0)2+D(X-X0)3
Transformation of coordinates:
  y=(Y-Y0)/(Y1-Y0) and x=(X-X0)/(X1-X0)
=>Y=Y0+(Y1-Y0)y    y=f(x)=a+bx1+cx2+dx3
                   (f'(x)=b+2cx1+3dx2)

=>F(X)=Y0+(Y1-Y0)f((X-X0)/(X1-X0))

M=m(Y1-Y0)/(X1-X0) =>m=M(X1-X0)/(Y1-Y0)

F(X0)=Y0  =>f(0)=0 (1)
F'(X0)=M0 =>f'(0)=m0 (2)
F(X1)=Y1  =>f(1)=1 (3)
F'(X1)=M1 =>f'(1)=m1 (4)

(1) f(0)=a+b(0)+c(0)+d(0)=0 =>a=0
(2) f'(0)=b+2c(0)+3d(0)=m0  =>b=m0
(3) f(1)=m0+c+d=1     =>1c+1d=1-m0 (5)
(4) f'(1)=m0+2c+3d=m1 =>2c+3d=m1-m0 (6)
(5)×3: 3c+3d=3(1-m0)
         (5)×3-(6)×1:   =>c=3-2m0-m1
(6)×1: 2c+3d=m1-m0
         (6)×1-(5)×2:   =>d=m0+m1-2
(5)×2: 2c+2d=2(1-m0)
=>y=f(x)=m0x1+(3-2m0-m1)x2-(2-m0-m1)x3
f'(x)=m0+2(3-2m0-m1)x1-3(2-m0-m1)x2
f''(x)=2(3-2m0-m1)-6(2-m0-m1)x1
Inflection point? 
=>x=(2m0+m1-3)/(3m0+3m1-6)

Input: X0,Y0,X1,Y1,M0,M1
m0=M0*(X1-X0)/(Y1-Y0);
m1=M1*(X1-X0)/(Y1-Y0);
for (X=X0; X<=X1; X++)
{ x=(X-X0)/(X1-X0);
  y=m0*x
   +(3-2*m0-m1)*x*x
   -(2-m0-m1)*x*x*x;
  Y=Y0+(Y1-Y0)*y;
  SetPixel(X,Y);
}
//Caution: 
//if |m0|>3 or |m1|>3
//then y could go
//outside of [0..1]

  2. CUBIC SPLINE ROTATED:   Approach - 3rd degree polynomial (rotated 45°)

PseudoCode CubicSpline

Y=F(X)=A+B(X-X0)1+C(X-X0)2+D(X-X0)3
Transformation of coordinates:
  xt=(X-X0)/(X1-X0) and yt=(Y-Y0)/(Y1-Y0)
  x=(yt+xt)/2 and y=(yt-xt)/2
  xt=x-y and yt=x+y
=>X=X0+(X1-X0)(x-y)   Y=Y0+(Y1-Y0)(x+y)
y=f(x)=a+bx1+cx2+dx3 (f'(x)=b+2cx1+3dx2)

M=mt(Y1-Y0)/(X1-X0) =>mt=M(X1-X0)/(Y1-Y0)
m=(mt-1)/(mt+1)   mt=(1+m)/(1-m)

F(X0)=Y0  =>f(0)=0 (1)
F'(X0)=M0 =>f'(0)=m0 (2)
F(X1)=Y1  =>f(1)=0 (3)
F'(X1)=M1 =>f'(1)=m1 (4)

(1) f(0)=a+b(0)+c(0)+d(0)=0 =>a=0
(2) f'(0)=b+2c(0)+3d(0)=m0  =>b=m0
(3) f(1)=m0+c+d=0     =>1c+1d=0-m0 (5)
(4) f'(1)=m0+2c+3d=m1 =>2c+3d=m1-m0 (6)
(5)×3: 3c+3d=3(0-m0)
         (5)×3-(6)×1:   =>c=0-2m0-m1
(6)×1: 2c+3d=m1-m0
         (6)×1-(5)×2:   =>d=m0+m1-0
(5)×2: 2c+2d=2(0-m0)
=>y=f(x)=m0x1+(0-2m0-m1)x2-(0-m0-m1)x3
f'(x)=m0+2(0-2m0-m1)x1-3(0-m0-m1)x2
f''(x)=2(0-2m0-m1)-6(0-m0-m1)x1
Inflection point? 
=>x=(2m0+m1-0)/(3m0+3m1-0)

Input: X0,Y0,X1,Y1,M0,M1
mt0=M0*(X1-X0)/(Y1-Y0);
mt1=M1*(X1-X0)/(Y1-Y0);
m0=(mt0-1)/(mt0+1);
m1=(mt1-1)/(mt1+1);
dx=Math.abs(X1-X0);
dx+=Math.abs(Y1-Y0);
dx=1/dx;
for (x=0; x<=1; x+=dx)
{ y=m0*x
   +(0-2*m0-m1)*x*x
   -(0-m0-m1)*x*x*x;
  X=X0+(X1-X0)*(x-y); 
  Y=Y0+(Y1-Y0)*(x+y);
  SetPixel(X,Y);
}

  3. BEZIER SPLINE:      Approach - Parametric curve

PseudoCode BezierSpline

P=[Px,Py]=[X,Y]
P(t)=( ((1-t)P0 + tL0))(1-t)
      +((1-t)L0 + tL1))t    )(1-t)
    +( ((1-t)L0 + tL1))(1-t)
      +((1-t)L1 + tP1))t    )t
                             t∈[0..1]
Transformation of coordinates:
original CS => orthonormal CS =>
normalize M-vector => back to o. CS
y=(Y-Y0)/(Y1-Y0) and x=(X-X0)/(X1-X0)
M=m(Y1-Y0)/(X1-X0) =>m=M(X1-X0)/(Y1-Y0)
n(m)=[nx,ny]=[1,m]/√(1+m2)
(n is vector of length 1)
N(M)=[Nx,Ny]=[(X1-X0),m(Y1-Y0)]/√(1+m2)
=>N0=[(X1-X0),M0(X1-X0)]/√(1+m02)
=>N1=[(X1-X0),M1(X1-X0)]/√(1+m12)

=>L0=P0+kN0 and L1=P1-kN1
(k ≈ 1/3 is free parameter to choose)
=>
P(t)=(1-t)3P0+3(1-t)2tL0+3(1-t)t2L1+t3P1

M0-M1-Intersection? 0+m0x+m1(1-x)=1 =>
x=(1-m1)/(m0-m1) and y=m0(1-m1)/(m0-m1)

Inflection point? 
P'(t) × P''(t) = 0 (cross product)
=>P'x(t) P''y(t) = P'y(t) P''x(t)
Q0(t)=P0+t(L0-P0)
Q1(t)=L0+t(L1-L0)
Q2(t)=L1+t(P1-L1)
R0(t)=Q0(t)+t(Q1(t)-Q0(t)) 
R1(t)=Q1(t)+t(Q2(t)-Q1(t))
P(t)=(R0(t)+t(R1(t)-R0(t)))
=> P(t)=
P0+3t(L0-P0)+3t2(L1-2L0+P0)+t3(P1-3L1+3L0-P0)
=> P'(t)=
3(L0-P0)+6t(L1-2L0+P0)+3t2(P1-3L1+3L0-P0)
Substitution: 
A=(L0-P0), B=(L1-2L0+P0), C=(P1-3L1+3L0-P0)
=> P'(t)=3A+6Bt+3Ct2 and P''(t)=6B+6Ct

(3Ax+6Bxt+3Cxt2)(6By+6Cyt)=
(3Ay+6Byt+3Cyt2)(6Bx+6Cxt)

=> (BxCy-ByCx)t2+(AxCy-AyCx)t+(AxBy-AyBx)=0
Substitution: a t2 + b t + c = 0
=> t=-b/(2a)±√(b2/(2a)2-c/a) (a≠0)
or t=-c/b (a=0, b≠0)              t∈[0..1]

Input: X0,Y0,X1,Y1,M0,M1
k=1/3;
m0=M0*(X1-X0)/(Y1-Y0);
m1=M1*(X1-X0)/(Y1-Y0);
n0=Math.sqrt(1+m0*m0);
n1=Math.sqrt(1+m1*m1);
L0X=X0+k*(X1-X0)/n0;
L0Y=Y0+k*M0*(X1-X0)/n0;
L1X=X1-k*(X1-X0)/n1;
L1Y=Y1-k*M1*(X1-X0)/n1;
dt=Math.abs(X1-X0);
dt+=Math.abs(Y1-Y0);
dt=1/dt;
CX=3*(L0X-X0);
BX=3*(L1X-2*L0X+X0);
AX=(X1-3*L1X+3*L0X-X0);
CY=3*(L0Y-Y0);
BY=3*(L1Y-2*L0Y+Y0);
AY=(Y1-3*L1Y+3*L0Y-Y0);
for (t=0; t<=1; t+=dt)
{ t2=t*t; t3=t2*t;
  X=X0+t*CX+t2*BX+t3*AX;
  Y=Y0+t*CY+t2*BY+t3*AY;
  SetPixel(X,Y);
}