home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / magazine / pctech / 1986_08 / smooth.pas < prev   
Pascal/Delphi Source File  |  1986-05-19  |  9KB  |  322 lines
  1. program smooth_curve;
  2.  
  3. { Turbo Pascal demonstration of parametric cubic splines }
  4. { Michael A. Covington                              1986 }
  5.  
  6. { This program draws a smooth curve connecting points    }
  7. { marked on the screen by the user.  The points do not   }
  8. { have to define a function; internally, both X and Y    }
  9. { are functions of T, the approximate distance traveled  }
  10. { as the curve moves around the screen.                  }
  11.  
  12. { NOTE: We are using 640x200 graphics but treating it as }
  13. { 320x200.  Thus all X coordinates are multiplied by 2   }
  14. { before being plotted on the screen.  The graphical     }
  15. { cursor occupies the odd-numbered columns in which data }
  16. { is never plotted.                                      }
  17.  
  18. const arraysize  = 100;  { maximum number of known points }
  19.       splinesize = 400;  { number of points that are
  20.                             calculated to draw the curve }
  21.  
  22. type  point_array  = array[1..arraysize] of real;
  23.       spline_array = array[1..splinesize] of real;
  24.  
  25. var x,                   { x-coordinates of known points }
  26.     y,                   { y-coordinates of known points }
  27.     t:     point_array;  { t-coordinates of known points }
  28.  
  29.     count: integer;      { number of known points }
  30.  
  31.     calcx,               { x-coordinates of calculated points }
  32.     calcy,               { y-coordinates of calculated points }
  33.     calct: spline_array; { t-coordinates of calculated points }
  34.       { CALCT is not actually used in this program, but is  }
  35.       { included because other similar programs may use it. }
  36.  
  37.     reply: char;         { user's reply to 'Again?' }
  38.  
  39. procedure get_points;
  40.  
  41.    { Allows the user to move a cursor around the screen and }
  42.    { mark points. X, Y, and T are recorded for each point.  }
  43.  
  44.    const    esc     = #27;
  45.             enter   = ^M;
  46.             up      = 'H';
  47.             down    = 'P';
  48.             left    = 'K';
  49.             right   = 'M';
  50.             beep    = ^G;
  51.  
  52.    var
  53.       key: char;
  54.       ix,              { current screen position }
  55.       iy:  integer;
  56.  
  57.    function inkey: char;
  58.      { Waits for a key to be pressed and returns its code.  }
  59.      { Returns the second byte if a 2-byte key is pressed.  }
  60.    var
  61.      c: char;
  62.    begin
  63.      read(KBD,c);
  64.      if keypressed then read(KBD,c);
  65.      inkey := c
  66.    end;
  67.  
  68.    procedure graphics_cursor(ix,iy,color: integer);
  69.      { Puts the graphics cursor at ix, iy,      }
  70.      { which are x and y expressed as integers. }
  71.    var
  72.      i,dy,dx: integer;
  73.    begin
  74.      for i:=1 to 5 do
  75.        begin
  76.          dy:=i;
  77.          dx:=i+i+1;
  78.          plot(ix+ix+dx,iy+dy,color);
  79.          plot(ix+ix+dx,iy-dy,color);
  80.          plot(ix+ix-dx,iy+dy,color);
  81.          plot(ix+ix-dx,iy-dy,color)
  82.        end
  83.    end;
  84.  
  85.    function dist(x1,y1,x2,y2:real):real;
  86.       { Distance between two points }
  87.    begin
  88.       dist := sqrt(sqr(x1-x2)+sqr(y1-y2))
  89.    end;
  90.  
  91.    function same_point_again:boolean;
  92.       { True if the current point is the same as }
  93.       { the most recently entered point.         }
  94.    begin
  95.       if count < 1 then
  96.          same_point_again := false
  97.       else
  98.          same_point_again :=
  99.             (ix = x[count]) and (iy = y[count])
  100.    end;
  101.  
  102. begin { get_points }
  103.  
  104.    textmode;
  105.    clrscr;
  106.    writeln('Move cursor with arrow keys.');
  107.    writeln('Mark points with Enter key.');
  108.    writeln('Press Esc when finished.');
  109.    writeln;
  110.    writeln('You must mark at least 4 points;');
  111.    writeln('Esc will not work until you do.');
  112.    writeln;
  113.    writeln('You cannot mark more than ',arraysize,' points.');
  114.    writeln('You cannot mark the same point twice unless');
  115.    writeln('you have marked a different point in between.');
  116.    writeln;
  117.    writeln('Erroneous input will result in a beep.');
  118.    writeln;
  119.    writeln('Press Enter to begin...');
  120.    readln;
  121.  
  122.  
  123.    hires;
  124.    hirescolor(yellow);
  125.    { graphbackground(black); Uncomment to run on EGA }
  126.  
  127.    ix:=160;
  128.    iy:=100;
  129.    count:=0;
  130.    key := ' ';
  131.  
  132.    while (count<4) or (key<>esc) do
  133.      begin
  134.        gotoxy(1,25);
  135.        write('Point ',count+1,'  x =',ix,' y =',iy,
  136.                          ' (Enter to mark, Esc to quit)   ');
  137.  
  138.        graphics_cursor(ix,iy,1);
  139.        key := inkey;
  140.        graphics_cursor(ix,iy,0); { erase previous cursor }
  141.  
  142.        case key of
  143.          up:    iy:=iy-1;
  144.          down:  iy:=iy+1;
  145.          left:  ix:=ix-1;
  146.          right: ix:=ix+1;
  147.          enter: if same_point_again or (count=arraysize) then
  148.                   write(beep)
  149.                 else
  150.                   begin
  151.                     count := count+1;
  152.                     x[count] := ix;
  153.                     y[count] := iy;
  154.                     if count=1 then
  155.                         t[count] := 0
  156.                     else
  157.                         t[count] := t[count-1] +
  158.                              dist(x[count],y[count],
  159.                                   x[count-1],y[count-1]);
  160.                     plot(2*ix,iy,1)
  161.                   end
  162.          else  { do nothing };
  163.        end { case }
  164.    end
  165. end;
  166.  
  167.  
  168. procedure spline(f,t:point_array;
  169.                  count:integer;
  170.                  var calcf,calct:spline_array);
  171.  
  172. { Fits a cubic spline to F[1..COUNT] as a function of  }
  173. { T[1..COUNT], then calculates equally spaced values   }
  174. { of T, places them in CALCT[1..SPLINESIZE], computes  }
  175. { the corresponding interpolated values of F, and      }
  176. { returns them in CALCF[1..SPLINESIZE].                }
  177.  
  178. { At least 4 points must be given. The points are }
  179. { assumed to be sorted in order of increasing T.  }
  180.  
  181. { Adapted from the FORTRAN routine CUBSPL (Carl de Boor,   }
  182. { A PRACTICAL GUIDE TO SPLINES, New York: Springer, 1978). }
  183.  
  184. var
  185.    a,b,c: point_array;  { spline coefficients }
  186.    d,e,g,h,             { used in finding spline coefs. }
  187.    wt,dt:    real;      { used in calculating points to plot }
  188.    i,j:      integer;   { loop counters }
  189.  
  190. begin
  191.  
  192.    { Compute first differences of T and F }
  193.  
  194.    for i:=2 to count do
  195.      begin
  196.        b[i] := t[i]-t[i-1];
  197.        c[i] := (f[i]-f[i-1])/b[i]
  198.      end;
  199.  
  200.    { Take care of beginning of curve }
  201.  
  202.    c[1] := b[3];
  203.    b[1] := b[2] + b[3];
  204.    a[1] := ((b[2] + 2*b[1])*c[2]*b[3] + b[2]*b[2]*c[3]) / b[1];
  205.  
  206.    { Forward pass of Gaussian elimination }
  207.  
  208.    for i:=2 to count-1 do
  209.      begin
  210.        g    := -b[i+1] / c[i-1];
  211.        a[i] := g*a[i-1] + 3*(b[i]*c[i+1] + b[i+1]*c[i]);
  212.        c[i] := g*b[i-1] + 2*(b[i]+b[i+1])
  213.      end;
  214.  
  215.    { Take care of end of curve }
  216.  
  217.    g := b[count-1] + b[count];
  218.    a[count] := ((b[count]+g+g) * c[count] * b[count-1]
  219.                 + b[count] * b[count]
  220.                 * (f[count-1]-f[count-2]) / b[count-1]) / g;
  221.    g := -g / c[count-1];
  222.    c[count] := b[count-1];
  223.  
  224.    { Complete the forward pass }
  225.  
  226.    c[count] := g*b[count-1] + c[count];
  227.    a[count] := (g*a[count-1] + a[count])/c[count];
  228.  
  229.    { Back substitution }
  230.  
  231.    for i := count-1 downto 1 do
  232.       a[i] := (a[i]-b[i]*a[i+1]) / c[i];
  233.  
  234.    { Generate cubic coefficients }
  235.  
  236.    for i:=2 to count do
  237.      begin
  238.        d := (f[i]-f[i-1])/b[i];
  239.        e := a[i-1] + a[i] - 2*d;
  240.        b[i-1] := 2*(d-a[i-1]-e)/b[i];
  241.        c[i-1] := (e/b[i])*(6/b[i])
  242.      end;
  243.  
  244.    { Compute the points to be plotted as function of WT }
  245.  
  246.    wt := 0;
  247.    dt := t[count] / (splinesize-1);
  248.    j  := 1;
  249.  
  250.    for i:=1 to splinesize do
  251.      begin
  252.        { Ensure that wt is between t[j] and t[j+1] }
  253.           while t[j+1] < wt do j:=j+1;
  254.        { Calculate a point }
  255.           calct[i] := wt;
  256.           h := wt - t[j];
  257.           calcf[i] := f[j]+h*(a[j]+h*(b[j]+h*c[j]/3)/2);
  258.        { Move to next point }
  259.           wt := wt+dt
  260.      end
  261. end;
  262.  
  263.  
  264. procedure fit_curve;
  265. begin
  266.   gotoxy(60,25);
  267.   write('Calculating...');
  268.   spline(x,t,count,calcx,calct);
  269.   spline(y,t,count,calcy,calct)
  270. end;
  271.  
  272.  
  273. procedure display;
  274.  
  275. { Displays the fitted curve and the original points. }
  276.  
  277. var i,k,ix,iy: integer;
  278.  
  279. begin
  280.  
  281.   hires; hirescolor(yellow);
  282.   { graphbackground(black); Uncomment to run on an EGA }
  283.  
  284.   { Draw crosses at the original points }
  285.      for i:=1 to count do
  286.        begin
  287.          ix := 2*round(x[i]);
  288.          iy := round(y[i]);
  289.          draw(ix-4,iy-2,ix+4,iy+2,1);
  290.          draw(ix+4,iy-2,ix-4,iy+2,1)
  291.        end;
  292.  
  293.   { Plot the spline curve }
  294.      for i:=2 to splinesize do
  295.         draw(round(2*calcx[i-1]),
  296.              round(calcy[i-1]),
  297.              round(2*calcx[i]),
  298.              round(calcy[i]),
  299.              1);
  300.  
  301.   gotoxy(1,25);
  302.   write('  (Press Enter when finished viewing)   ');
  303.   readln
  304. end;
  305.  
  306.  
  307. { main program }
  308.  
  309. begin
  310.     repeat
  311.  
  312.       get_points;
  313.       fit_curve;
  314.       display;
  315.  
  316.       textmode;
  317.       write('Again? (Y/N)  ');
  318.       readln(reply)
  319.  
  320.     until reply in ['N','n']
  321. end.
  322.