[Back to GRAPHICS SWAG index] [Back to Main SWAG index] [Original]
{ At last! A couple of spline routines that CONNECT the data points.}
{they are set up in a simple, but probably useless, manner for you to
configure to suit. The reconfigurations are simple, it's the core
routines that seem to be hard to get.}
uses crt,graph;
const maxnodes=4;
sigma:real=1.0; {tension factor}
total=100;
type data=array[0..maxnodes]of real;
var xmin,xmax,sum,xx,yy:real;
i,n,eger:integer;
xc,yc,xp,yp,temp,path:data;
GraphDriver,GraphMode,ErrorCode:Integer;
ch:char;
alternate:boolean;
Procedure Spline1(color:integer);
{ Fits a smooth curve through a given set of points. Small, simple, fast.
No tension factor.
From Algorithms, Robert Sedgewick, sort of...
simple, but limited
Beat beyond recognition by Ron Nossaman June 1996 }
Var i,j,oldy,oldx,x,y:integer;
part,t,xx,yy:real;
d,wy,wx:data;
Function f(g:real):real;
begin
f:=g*g*g-g;
end;
Begin
setcolor(color);
oldx:=999;
x:=round(xc[1]);
y:=round(yc[1]);
{calculate matrix for x & y}
for i:=1 to maxnodes-1 do d[i]:=4; {fake ascending sequence}
for i:=1 to maxnodes-1 do
begin
wy[i]:=6*((yc[i+1]-yc[i])-(yc[i]-yc[i-1]));
wx[i]:=6*((xc[i+1]-xc[i])-(xc[i]-xc[i-1]));
end;
yp[0]:=0.0; xp[0]:=0.0;
yp[maxnodes]:=0.0; xp[maxnodes]:=0.0;
for i:=1 to maxnodes-2 do
begin
wy[i+1]:=wy[i+1]-wy[i]*0.25;
wx[i+1]:=wx[i+1]-wx[i]*0.25;
d[i+1]:=d[i+1]-0.25;
end;
for i:=maxnodes-1 downto 1 do
begin
yp[i]:=(wy[i]-yp[i+1])/d[i];
xp[i]:=(wx[i]-xp[i+1])/d[i];
end;
{ Draw spline }
for i:=0 to maxnodes-1 do
begin
for j:=0 to 30 do {arbitrary number of steps between points}
begin
t:=j/30;
part:=t*yc[i+1]+(1-t)*yc[i]+(f(t)*yp[i+1]+f(1-t)*yp[i])/6.0;
y:=round(part);
part:=t*xc[i+1]+(1-t)*xc[i]+(f(t)*xp[i+1]+f(1-t)*xp[i])/6.0;
x:=round(part);
if oldx<>999 then line(oldx,oldy,x,y);
oldx:=x;
oldy:=y;
end;
end;
end;
Procedure Spline2(color:integer);
{ Fits a smooth curve through a given set of points. Small, simple, fast.
No tension factor. Nicer curve than above Spline1 routine.
From Algorithms, Robert Sedgewick, sort of...
Beat beyond recognition by Ron Nossaman June 1996 }
Var i,oldy,oldx,x,y,j:integer;
part,t,xx,yy:real;
zc,dx,dy,u,wx,wy,px,py:data;
Function f(g:real):real;
begin
f:=g*g*g-g;
end;
Begin
setcolor(color);
oldx:=999;
x:=round(xc[0]);
y:=round(yc[0]);
zc[0]:=0.0;
for i:=1 to maxnodes do
begin
xx:=xc[i-1]-xc[i]; yy:=yc[i-1]-yc[i];
t:=sqrt(xx*xx+yy*yy);
zc[i]:=zc[i-1]+t; {establish a proportional linear progression}
end;
{calculate x & y matrix stuff}
for i:=1 to maxnodes-1 do
begin
dx[i]:=2*(zc[i+1]-zc[i-1]);
dy[i]:=2*(zc[i+1]-zc[i-1]);
end;
for i:=0 to maxnodes-1 do
begin
u[i]:=zc[i+1]-zc[i];
end;
for i:=1 to maxnodes-1 do
begin
wy[i]:=6*((yc[i+1]-yc[i])/u[i]-(yc[i]-yc[i-1])/u[i-1]);
wx[i]:=6*((xc[i+1]-xc[i])/u[i]-(xc[i]-xc[i-1])/u[i-1]);
end;
py[0]:=0.0; px[0]:=0.0; px[1]:=0; py[1]:=0;
py[maxnodes]:=0.0; px[maxnodes]:=0.0;
for i:=1 to maxnodes-2 do
begin
wy[i+1]:=wy[i+1]-wy[i]*u[i]/dy[i];
dy[i+1]:=dy[i+1]-u[i]*u[i]/dy[i];
wx[i+1]:=wx[i+1]-wx[i]*u[i]/dx[i];
dx[i+1]:=dx[i+1]-u[i]*u[i]/dx[i];
end;
for i:=maxnodes-1 downto 1 do
begin
py[i]:=(wy[i]-u[i]*py[i+1])/dy[i];
px[i]:=(wx[i]-u[i]*px[i+1])/dx[i];
end;
{ Draw spline }
for i:=0 to maxnodes-1 do
begin
for j:=0 to 30 do
begin
part:=zc[i]-(((zc[i]-zc[i+1])/30)*j);
t:=(part-zc[i])/u[i];
part:=t*yc[i+1]+(1-t)*yc[i]+u[i]*u[i]*(f(t)*py[i+1]+f(1-t)*py[i])/6.0;
y:=round(part);
part:=zc[i]-(((zc[i]-zc[i+1])/30)*j);
t:=(part-zc[i])/u[i];
part:=t*xc[i+1]+(1-t)*xc[i]+u[i]*u[i]*(f(t)*px[i+1]+f(1-t)*px[i])/6.0;
x:=round(part);
if oldx<>999 then line(oldx,oldy,x,y);
oldx:=x;
oldy:=y;
end;
end;
end;
(* -----------------------------------------------------------------*)
{ Spline under tension}
{ Original Author -- Copyright(c)1985 James R .Van Zandt
Converted to Turbo Pascal, simplified, altered
- Ron Nossaman June 1996 nossaman@southwind.net
Based on algorithms by A. K. Cline }
procedure curv1(var x:data; var y:data; var p:data; n:integer);
var i:Integer;
c1,c2,c3,deln,delnm1,delnn,dels,delx1,delx2,delx12,
diag1,diag2,diagin,dx1,dx2,exps,
sigmap,sinhs,sinhin,slpp1,slppn,spdiag:real;
begin
delx1:=x[1]-x[0];
dx1:=(y[1]-y[0])/delx1;
if sigma>=0.then
begin
slpp1:=0;
slppn:=0;
end else
begin
if n<>1 then
begin
delx2:=x[2]-x[1];
delx12:=x[2]-x[0];
c1:=-(delx12+delx1)/delx12/delx1;
c2:=delx12/delx1/delx2;
c3:=-delx1/delx12/delx2;
slpp1:=c1* y[0]+c2* y[1]+c3* y[2];
deln:=x[n-1]-x[n-2];
delnm1:=x[n-2]-x[n-3];
delnn:=x[n-1]-x[n-3];
c1:=(delnn+deln)/delnn/deln;
c2:=-delnn/deln/delnm1;
c3:=deln/delnn/delnm1;
slppn:=c3* y[n-3]+c2* y[n-2]+c1* y[n-1];
end else
begin
p[0]:=0.;
p[1]:=0.;
end;
end;
(* denormalize tension factor *)
sigmap:=abs(sigma)*(n-1)/(x[n-1]-x[0]);
(* set up right hand side and tridiagonal system for
yp and perform forward elimination *)
dels:=sigmap* delx1;
exps:=exp(dels);
sinhs:=0.5*(exps-1./exps);
sinhin:=1./(delx1* sinhs);
diag1:=sinhin*(dels* 0.5*(exps+1./exps)-sinhs);
diagin:=1./diag1;
p[0]:=diagin*(dx1-slpp1);
spdiag:=sinhin*(sinhs-dels);
temp[0]:=diagin* spdiag;
if(n <> 1) then
begin
for i:=1 to n-2 do
begin
delx2:=x[i+1]-x[i];
dx2:=(y[i+1]-y[i])/delx2;
dels:=sigmap*delx2;
exps:=exp(dels);
sinhs:=0.5*(exps-1./exps);
sinhin:=1./(delx2*sinhs);
diag2:=sinhin*(dels*(0.5*(exps+1./exps))-sinhs);
diagin:=1./(diag1+diag2-spdiag*temp[i-1]);
p[i]:=diagin*(dx2-dx1-spdiag*p[i-1]);
spdiag:=sinhin*(sinhs-dels);
temp[i]:=diagin*spdiag;
dx1:=dx2;
diag1:=diag2;
end;
end;
diagin:=1./(diag1-spdiag*temp[n-2]);
p[n-1]:=diagin*(slppn-dx2-spdiag*p[n-2]);
(* perform back substitution *)
for i:=n-2 downto 0 do p[i]:=p[i]-(temp[i]*p[i+1]);
end;
function curv2(var x:data; var y:data; var p:data; t:real; n:integer):real;
var i,j,i1:integer;
del1,del2,dels,exps,exps1,s,sigmap,sinhd1,sinhd2,sinhs:real;
begin
i1:=1;
s:=x[n-1]-x[0];
sigmap:=abs(sigma)*(n-1)/s;
i:=i1;
while(i<n) and(t>= x[i])do inc(i);
while(i>1) and(x[i-1]>t)do dec(i);
i1:=i;
del1:=t-x[i-1];
del2:=x[i]-t;
dels:=x[i]-x[i-1];
exps1:=exp(sigmap*del1); sinhd1:=0.5*(exps1-1./exps1);
exps:=exp(sigmap*del2); sinhd2:=0.5*(exps-1./exps);
exps:=exps1*exps; sinhs:=0.5*(exps-1./exps);
curv2:=((p[i]*sinhd1+p[i-1]*sinhd2)/sinhs+
((y[i]-p[i])*del1+(y[i-1]-p[i-1])*del2)/dels);
end;
procedure curv0(n,requested:integer);
var i,j,each,stop,seg,regressing:integer;
s,ds,xx,yy,oldx,oldy:real;
begin
oldx:=999;
curv1(path,xc,xp,n);
curv1(path,yc,yp,n);
s:=path[0];
seg:=0;
for j:=0 to n-2 do
begin
stop:=100;
ds:=(path[j+1]-path[j])/stop;
for i:=0 to stop-1 do
begin
xx:=round(curv2(path,xc,xp,s,n));
yy:=round(curv2(path,yc,yp,s,n));
if oldx<>999 then
line(round(oldx),round(oldy),round(xx),round(yy));
oldx:=xx; oldy:=yy;
s:=s+ds;
end;
end;
xx:=xc[n-1];
yy:=yc[n-1];
line(round(oldx),round(oldy),round(xx),round(yy));
end;
procedure tspline(color:integer);
{ Fits a smooth curve through a given set of points, using the splines
under tension introduced by J. Schweikert. They are similar to cubic
splines,but are less likely to introduce spurious inflection points.
Original Author -- Copyright(c)1985James R .Van Zandt
Converted to Turbo Pascal, butchered, simplified, altered --
Ron Nossaman June 1996 nossaman@southwind.net
Based on algorithms by A. K. Cline }
var nn,origin:integer;
begin
sum:=0;
path[0]:=0;
for i:=1 to maxnodes do
begin
xx:=xc[i]-xc[i-1];
yy:=yc[i]-yc[i-1];
sum:=sum+sqrt((xx*xx)+(yy*yy));
path[i]:=sum;
if alternate then path[i]:=i; {my addition, another alternative }
end;
setcolor(color); {display color for spline curve}
curv0(maxnodes+1,100);
end;
procedure loadarrays;
var i:integer;
begin
setcolor(red);
for i:=0 to maxnodes do
begin
xc[i]:=random(540)+50; {try to keep it on the screen for demo}
yc[i]:=random(380)+50;
{show the path, if you want}
if i>0 then line(round(xc[i-1]),round(yc[i-1]),round(xc[i]),round(yc[i]));
end;
setcolor(white);
for i:=0 to maxnodes do circle(round(xc[i]),round(yc[i]),3);
end;
begin
clrscr;
randomize;
ch:=#0;
GraphDriver:=VGA;
Graphmode:=2;
InitGraph(GraphDriver,GraphMode,'');
ErrorCode:=GraphResult;
If ErrorCode<>grOK then Halt(1);
repeat
setfillstyle(solidfill,black);
bar(0,0,639,479);
loadarrays;
{try whatever combinations strike you}
{there are a lot of different ways to get there!}
alternate:=false; sigma:=-2; tspline(green);
alternate:=false; sigma:=0.01; tspline(cyan);
alternate:=true; sigma:=1; tspline(white);
spline1(lightred);
spline2(yellow);
repeat until keypressed;
while keypressed do ch:=readkey;
until ch=#27;
closegraph;
end.
[Back to GRAPHICS SWAG index] [Back to Main SWAG index] [Original]