ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng)/ProfitPress-MegaCDROM2.B6I/MAGAZINE/MISC/AUG86.ZIP/SMOOTH.PAS

Pascal/Delphi Source File | 1986-05-19 | 8.7 KB | 322 lines

program smooth_curve;

{ Turbo Pascal demonstration of parametric cubic splines }
{ Michael A. Covington                              1986 }

{ This program draws a smooth curve connecting points    }
{ marked on the screen by the user.  The points do not   }
{ have to define a function; internally, both X and Y    }
{ are functions of T, the approximate distance traveled  }
{ as the curve moves around the screen.                  }

{ NOTE: We are using 640x200 graphics but treating it as }
{ 320x200.  Thus all X coordinates are multiplied by 2   }
{ before being plotted on the screen.  The graphical     }
{ cursor occupies the odd-numbered columns in which data }
{ is never plotted.                                      }

const arraysize  = 100;  { maximum number of known points }
      splinesize = 400;  { number of points that are
                            calculated to draw the curve }

type  point_array  = array[1..arraysize] of real;
      spline_array = array[1..splinesize] of real;

var x,                   { x-coordinates of known points }
    y,                   { y-coordinates of known points }
    t:     point_array;  { t-coordinates of known points }

    count: integer;      { number of known points }

    calcx,               { x-coordinates of calculated points }
    calcy,               { y-coordinates of calculated points }
    calct: spline_array; { t-coordinates of calculated points }
      { CALCT is not actually used in this program, but is  }
      { included because other similar programs may use it. }

    reply: char;         { user's reply to 'Again?' }

procedure get_points;

   { Allows the user to move a cursor around the screen and }
   { mark points. X, Y, and T are recorded for each point.  }

   const    esc     = #27;
            enter   = ^M;
            up      = 'H';
            down    = 'P';
            left    = 'K';
            right   = 'M';
            beep    = ^G;

   var
      key: char;
      ix,              { current screen position }
      iy:  integer;

   function inkey: char;
     { Waits for a key to be pressed and returns its code.  }
     { Returns the second byte if a 2-byte key is pressed.  }
   var
     c: char;
   begin
     read(KBD,c);
     if keypressed then read(KBD,c);
     inkey := c
   end;

   procedure graphics_cursor(ix,iy,color: integer);
     { Puts the graphics cursor at ix, iy,      }
     { which are x and y expressed as integers. }
   var
     i,dy,dx: integer;
   begin
     for i:=1 to 5 do
       begin
         dy:=i;
         dx:=i+i+1;
         plot(ix+ix+dx,iy+dy,color);
         plot(ix+ix+dx,iy-dy,color);
         plot(ix+ix-dx,iy+dy,color);
         plot(ix+ix-dx,iy-dy,color)
       end
   end;

   function dist(x1,y1,x2,y2:real):real;
      { Distance between two points }
   begin
      dist := sqrt(sqr(x1-x2)+sqr(y1-y2))
   end;

   function same_point_again:boolean;
      { True if the current point is the same as }
      { the most recently entered point.         }
   begin
      if count < 1 then
         same_point_again := false
      else
         same_point_again :=
            (ix = x[count]) and (iy = y[count])
   end;

begin { get_points }

   textmode;
   clrscr;
   writeln('Move cursor with arrow keys.');
   writeln('Mark points with Enter key.');
   writeln('Press Esc when finished.');
   writeln;
   writeln('You must mark at least 4 points;');
   writeln('Esc will not work until you do.');
   writeln;
   writeln('You cannot mark more than ',arraysize,' points.');
   writeln('You cannot mark the same point twice unless');
   writeln('you have marked a different point in between.');
   writeln;
   writeln('Erroneous input will result in a beep.');
   writeln;
   writeln('Press Enter to begin...');
   readln;


   hires;
   hirescolor(yellow);
   { graphbackground(black); Uncomment to run on EGA }

   ix:=160;
   iy:=100;
   count:=0;
   key := ' ';

   while (count<4) or (key<>esc) do
     begin
       gotoxy(1,25);
       write('Point ',count+1,'  x =',ix,' y =',iy,
                         ' (Enter to mark, Esc to quit)   ');

       graphics_cursor(ix,iy,1);
       key := inkey;
       graphics_cursor(ix,iy,0); { erase previous cursor }

       case key of
         up:    iy:=iy-1;
         down:  iy:=iy+1;
         left:  ix:=ix-1;
         right: ix:=ix+1;
         enter: if same_point_again or (count=arraysize) then
                  write(beep)
                else
                  begin
                    count := count+1;
                    x[count] := ix;
                    y[count] := iy;
                    if count=1 then
                        t[count] := 0
                    else
                        t[count] := t[count-1] +
                             dist(x[count],y[count],
                                  x[count-1],y[count-1]);
                    plot(2*ix,iy,1)
                  end
         else  { do nothing };
       end { case }
   end
end;


procedure spline(f,t:point_array;
                 count:integer;
                 var calcf,calct:spline_array);

{ Fits a cubic spline to F[1..COUNT] as a function of  }
{ T[1..COUNT], then calculates equally spaced values   }
{ of T, places them in CALCT[1..SPLINESIZE], computes  }
{ the corresponding interpolated values of F, and      }
{ returns them in CALCF[1..SPLINESIZE].                }

{ At least 4 points must be given. The points are }
{ assumed to be sorted in order of increasing T.  }

{ Adapted from the FORTRAN routine CUBSPL (Carl de Boor,   }
{ A PRACTICAL GUIDE TO SPLINES, New York: Springer, 1978). }

var
   a,b,c: point_array;  { spline coefficients }
   d,e,g,h,             { used in finding spline coefs. }
   wt,dt:    real;      { used in calculating points to plot }
   i,j:      integer;   { loop counters }

begin

   { Compute first differences of T and F }

   for i:=2 to count do
     begin
       b[i] := t[i]-t[i-1];
       c[i] := (f[i]-f[i-1])/b[i]
     end;

   { Take care of beginning of curve }

   c[1] := b[3];
   b[1] := b[2] + b[3];
   a[1] := ((b[2] + 2*b[1])*c[2]*b[3] + b[2]*b[2]*c[3]) / b[1];

   { Forward pass of Gaussian elimination }

   for i:=2 to count-1 do
     begin
       g    := -b[i+1] / c[i-1];
       a[i] := g*a[i-1] + 3*(b[i]*c[i+1] + b[i+1]*c[i]);
       c[i] := g*b[i-1] + 2*(b[i]+b[i+1])
     end;

   { Take care of end of curve }

   g := b[count-1] + b[count];
   a[count] := ((b[count]+g+g) * c[count] * b[count-1]
                + b[count] * b[count]
                * (f[count-1]-f[count-2]) / b[count-1]) / g;
   g := -g / c[count-1];
   c[count] := b[count-1];

   { Complete the forward pass }

   c[count] := g*b[count-1] + c[count];
   a[count] := (g*a[count-1] + a[count])/c[count];

   { Back substitution }

   for i := count-1 downto 1 do
      a[i] := (a[i]-b[i]*a[i+1]) / c[i];

   { Generate cubic coefficients }

   for i:=2 to count do
     begin
       d := (f[i]-f[i-1])/b[i];
       e := a[i-1] + a[i] - 2*d;
       b[i-1] := 2*(d-a[i-1]-e)/b[i];
       c[i-1] := (e/b[i])*(6/b[i])
     end;

   { Compute the points to be plotted as function of WT }

   wt := 0;
   dt := t[count] / (splinesize-1);
   j  := 1;

   for i:=1 to splinesize do
     begin
       { Ensure that wt is between t[j] and t[j+1] }
          while t[j+1] < wt do j:=j+1;
       { Calculate a point }
          calct[i] := wt;
          h := wt - t[j];
          calcf[i] := f[j]+h*(a[j]+h*(b[j]+h*c[j]/3)/2);
       { Move to next point }
          wt := wt+dt
     end
end;


procedure fit_curve;
begin
  gotoxy(60,25);
  write('Calculating...');
  spline(x,t,count,calcx,calct);
  spline(y,t,count,calcy,calct)
end;


procedure display;

{ Displays the fitted curve and the original points. }

var i,k,ix,iy: integer;

begin

  hires; hirescolor(yellow);
  { graphbackground(black); Uncomment to run on an EGA }

  { Draw crosses at the original points }
     for i:=1 to count do
       begin
         ix := 2*round(x[i]);
         iy := round(y[i]);
         draw(ix-4,iy-2,ix+4,iy+2,1);
         draw(ix+4,iy-2,ix-4,iy+2,1)
       end;

  { Plot the spline curve }
     for i:=2 to splinesize do
        draw(round(2*calcx[i-1]),
             round(calcy[i-1]),
             round(2*calcx[i]),
             round(calcy[i]),
             1);

  gotoxy(1,25);
  write('  (Press Enter when finished viewing)   ');
  readln
end;


{ main program }

begin
    repeat

      get_points;
      fit_curve;
      display;

      textmode;
      write('Again? (Y/N)  ');
      readln(reply)

    until reply in ['N','n']
end.