_Algorithm Alley_
by Robert F. Kauffmann

Listing One
{--------------------------TYPES------------------------------------}
type POINT = record                          { define a 2-D point}
        X, Y : integer;
     end;

     PT4 = array[1..5] of POINT; {Array of points controlling 1 seg }

     A4_4 = array[1..4,1..4] of real;     {Basis matrices for cubics}
     A1000_4 = array[1..1000,1..4] of real;  {Lookup table for basis
                                              curve points; up to 1000 points
                                              are pre-calculated in advance
                                              and stored for each of the
                                              four basis equations }
     Curve_Type = ( Cat_Rom,    {Cubic spline which passes through every
                                control point continuity in first derivative}
                    B_Spline,   {Cubic spline which does not pass through
                                control points; continuity in 2nd derivative}
                    Ext_Tspline,  {Triginometric spline which does not pass
                                  through control points; continuity in 
                                  all derivatives}
                    Int_Tspline ); {Cubic spline which passes through all 
                               control points; continuity in all derivatives}

{--------------------------VARIABLES--------------------------------}
var
   Max_Points    : real;    {Max number of points to interpolate per segment}
   Curve         : A1000_4; {Current interpolating curve}
   Curr_Basis    : Curve_Type; {Type of basis function being used}

{-------------------CUBIC CURVE BASIS MATRICES----------------------}
const
    Bspline : A4_4 =
          ( ( -0.166667,  0.5,       -0.5,        0.166667   ),
            (  0.5,      -1.0,        0.0,        0.666667   ),
            ( -0.5,       0.5,        0.5,        0.166667   ),
            (  0.166667,  0.0,        0.0,        0.0        ) );
    Catrom  : A4_4 =
          ( ( -0.5,  1.0,  -0.5,   0.0   ),
            (  1.5, -2.5,   0.0,   1.0   ),
            ( -1.5,  2.0,   0.5,   0.0   ),
            (  0.5, -0.5,   0.0,   0.0   ) );
{--------------------------------------------------------------------}
   function CubicSpline( t : real;
                         n : integer;
                         M : A4_4 ) : real;
  { Calculates value of nth cubic basis function specified by M for parameter
    t; calculates basis curve one point at a time; called by SetBasis }
   var T1, T2, T3 : real;
   begin
      T1 := t; T2 := t*t;T3 := t*t*t;
      CubicSpline := T3*M[n,1] + T2*M[n,2] + T1*M[n,3] + M[n,4];
   end; {CubicSpline}

 {-------------------------------------------------------------------}
   function ExTSpline( t : real;
                      n : integer ) : real;
   { Calculates the value of the n-th basis function for the
     exterpolating triginometric spline for parameter t; used to
     calculate the basis curve one point at a time; Called by
     SetBasis }
   var Sn, Cs : real;
   begin
       Sn := sin( 0.5*Pi*t ); Cs := cos( 0.5*Pi*t );
     case n of
       1 : ExTSpline := 0.25*( 1 - Cs);
       2 : ExTspline := 0.25*( 1 + Sn);
       3 : ExTspline := 0.25*( 1 + Cs);
       4 : ExTspline := 0.25*( 1 - Sn);
     end;
   end; {ExTspline}
 {-------------------------------------------------------------------}
   function InTspline( t : real;
                      n : integer ) : real;
   { Calculates the value of the n-th basis function for the
     interpolating triginometric spline for parameter t; used to
     calculate the basis curve one point at a time; Called by
     SetBasis }
   var Sn, Cs : real;
   begin
       Sn := sin( 0.5*Pi*t ); Cs := cos( 0.5*Pi*t );
     case n of
       1 : InTspline := 0.5*Cs*(Cs - 1 );
       2 : InTspline := 0.5*Sn*(Sn + 1 );
       3 : InTspline := 0.5*Cs*(Cs + 1 );
       4 : InTspline := 0.5*Sn*(Sn - 1 );
     end;
  end; {InTspline}
 {-------------------------------------------------------------------}
   procedure SetBasis( Basis : Curve_Type;
                       Max   : integer);
   {Pre-calculate basis function at specified granularity}
   var t : real;
       i, j : integer;
   begin
      Curr_Basis := Basis;
      if Max > 1000 then Max_Points := 1000  {set granularity <= 1000}
                    else Max_Points := Max;
      for i:= 1 to Max do begin
         t:= i/Max;
         case Basis of       {calculate basis curves on point a time}
            Ext_Tspline :   {exterpolating trig spline}
               for j:= 1 to 4 do Curve[i,j] := ExTspline( t, j );
            Int_Tspline :   {interpolating trig spline}
               for j:= 1 to 4 do Curve[i,j] := InTspline( t, j )
            Cat_Rom :       {Catmull-Rom cubic spline}
               for j:= 1 to 4 do Curve[i,j] := CubicSpline( t, j, CatRom );
            B_Spline :      {BSpline cubic spline}
               for j:= 1 to 4 do Curve[i,j] := CubicSpline( t, j, BSpline );
         end;
      end;
   end; {SetBasis}
 {-------------------------------------------------------------------}
  procedure PlotCurve ( P1,P2,P3,P4 : POINT;
                         C           : integer );
   {Plots a segment of the defined curve in the specified color;
    called by Render_Curve }
   var Xo, X1, X2, X3, X4, t, T1, T2, T3,
       Yo, Y1, Y2, Y3, Y4, X, Y   : real;
       i, j, Max           : integer;
   begin
      Max := round( Max_Points );
      X1 := P1.X; X2 := P2.X; X3 := P3.X; X4 := P4.X;
      Y1 := P1.Y; Y2 := P2.Y; Y3 := P3.Y; Y4 := P4.Y;
      for i:= 1 to Max do begin { vary thru [0,1] & plot basis functn}
          t:= i/Max;
          X := X1*Curve[i,1] + X2*Curve[i,2] +
               X3*Curve[i,3] + X4*Curve[i,4] ;
          Y := Y1*Curve[i,1] + Y2*Curve[i,2] +
               Y3*Curve[i,3] + Y4*Curve[i,4] ;
          put_pixel( X, Y, C );
      end;
   end; {PlotCurve}


3