// File: SPLINE3.BCPL
// P.Baudelaire
// December 5, 1977 4:49 PM
// The procedure ParametricSpline implements the algorithms
// for local cubic splines described in
// "Spline Curve Techniques"
// by P.Baudelaire, R.Flegal, & R.Sproull
// Xerox Internal Report (1977)
// Uses MICROCODED floating point routines
// outgoing procedures:
external [
ParametricSpline
PSerror
]
// outgoing statics:
external [
PSzone
]
static [
PSzone=0 // storage zone
]
// incoming procedures:
external [
FLD; FAD; FML; FST // FLOAT (Alto floating point package)
FLDI; FSB; FDV
FCM; FNEG; FTR
FPSetup
Allocate // Alto SYSTEM
Free
MoveBlock
Zero
]
// incoming statics:
external [
FPwork // FLOAT (floating point accumulator)
]
// local definitions:
manifest [
naturalSpline=0
periodicSpline=1
// floating point registers: 0 to 4
ac0=0; ac1=1; ac2=2; ac3=3; ac4=4
// constants:
zero=5; one=6; two=7; six=8
numFPacs=9
// local spline type
CatmullRom=1
FourPoints=2
Bspline1=3
Bspline2=4
Bspline3=5
Hermit=6
]
structure PSVEC [
FPworkSave word
FPworkNew word
fpx word
fpy word
a word
]
manifest lPSVEC=size PSVEC/16
// local statics:
static [
// working storage
PSvec=0
// various parameters
defaultPStype=Bspline3
PScoefNums=0
PScoefDens=0
PScoefEndNum=1
PScoefEndDen=1
HermitCoef=3
]
let ParametricSpline(n, px, py, p1x, p2x, p3x, p1y, p2y, p3y, endCondition, PStype;
numargs nargs) = valof [
// default arguments, get storage, check various things
let tempVec= vec lPSVEC
if PSinit(tempVec) eq 0 resultis 0
PScoefNums= table [ 1; 1; 1; 1; 1; 1 ]
PScoefDens= table [ 2; 6; 1; 2; 6; 2 ]
let convertXY= n ls 0
if convertXY then n=-n
let x=PSallocate(lv(PSvec>>PSVEC.fpx), 2*n+4)
let y=PSallocate(lv(PSvec>>PSVEC.fpy), 2*n+4)
if (x eq 0) % (y eq 0) resultis 0
test convertXY
ifso [
// convert coordinates from integer to floating point
for i=0 to n-1 do [
FST(FLDI(ac1, x!i), px+2*(i+1))
FST(FLDI(ac1, y!i), py+2*(i+1))
]
]
ifnot [
MoveBlock(x+2, px, 2*n)
MoveBlock(y+2, py, 2*n)
]
switchon nargs into [
case 9:
endCondition=naturalSpline
case 10:
PStype=defaultPStype
case 11:
if endCondition ne naturalSpline &
endCondition ne periodicSpline resultis PSquit(PSerror(3))
if n ls 3 then endCondition=naturalSpline
if PStype eq 0 then PStype=defaultPStype
endcase
default:
resultis PSquit(PSerror(4))
]
test endCondition eq periodicSpline
ifso [
if ((FCM(FLD(ac1, x+2), x+2*n) ne 0) %
(FCM(FLD(ac2, y+2), y+2*n) ne 0)) resultis PSquit(PSerror(2))
MoveBlock(x, x+2*(n-1), 2)
MoveBlock(y, y+2*(n-1), 2)
MoveBlock(x+2*(n+1), x+4, 2)
MoveBlock(y+2*(n+1), y+4, 2)
]
ifnot [
// ac3=b=PScoefEndNum/PScoefEndDen
FDV(FLDI(ac3, PScoefEndNum), FLDI(ac2, PScoefEndDen))
// x(0)=x(1)+b*(x(1)-x(2))
FST(FAD(FML(FSB(FLD(ac1, x+2), x+4), ac3), x+2), x)
FST(FAD(FML(FSB(FLD(ac1, y+2), y+4), ac3), y+2), y)
// x(n+1)=x(n)+b*(x(n)-x(n-1))
FST(FAD(FML(FSB(FLD(ac1, x+2*n), x+2*(n-1)), ac3), x+2*n), x+2*(n+1))
FST(FAD(FML(FSB(FLD(ac1, y+2*n), y+2*(n-1)), ac3), y+2*n), y+2*(n+1))
]
let a=PSallocate(lv(PSvec>>PSVEC.a), 32)
if (a eq 0) resultis 0
Zero(a, 32)
let c=nil
let saveXY=false
FDV(FLDI(ac2, PScoefNums!(PStype-1)), FLDI(ac3, PScoefDens!(PStype-1)))
test PStype eq CatmullRom
ifnot [
switchon PStype into [
case FourPoints:
c=table [
-1; 3; -3; 1;
3; -6; 3; 0;
-2; -3; 6; -1;
0; 6; 0; 0 ]
endcase
case Bspline1:
c=table [
0; 0; 0; 0;
0; 0; 0; 0;
0; -1; 1; 0;
0; 1; 0; 0 ]
endcase
case Bspline2:
c=table [
0; 0; 0; 0;
1; -2; 1; 0;
-2; 2; 0; 0;
1; 1; 0; 0 ]
saveXY=true
endcase
case Bspline3:
c=table [
-1; 3; -3; 1;
3; -6; 3; 0;
-3; 0; 3; 0;
1; 4; 1; 0 ]
saveXY=true
endcase
case Hermit:
c=table [
0; 0; 0; 0;
0; 0; 0; 0;
0; 0; 0; 0;
1; 1; 0; 0 ]
c!0=2-HermitCoef
c!2=HermitCoef-2
c!4=2*HermitCoef-3
c!5=-HermitCoef
c!6=3-HermitCoef
c!8=-HermitCoef
c!9=HermitCoef
saveXY=true
endcase
]
for i=0 to 15 do [
FST(FML(FLDI(ac1, c!i), ac2), a+2*i)
]
]
ifso [
// case CatmullRom:
FST(FST(ac2, a+6), a+20)
FST(FST(FST(FNEG(FLD(ac1, ac2)), a), a+14), a+16)
FST(FNEG(FST(FSB(FLDI(ac1, 2), ac2), a+2)), a+4)
FST(FNEG(FSB(FLDI(ac1, 3), ac2)), a+10)
FST(FSB(FLDI(ac1, 3), FST(FAD(ac2, ac2), a+8)), a+12)
FST(one, a+26)
]
for i=0 to 3 do [
FST(FML(FLD(ac1, a+2*i), six), a+2*i)
FST(FML(FLD(ac1, a+2*(i+4)), two), a+2*(i+4))
]
//do everything twice to get x(t) and y(t).
let p1,p2,p3,p,s=nil,nil,nil,nil,nil
for t=1 to 2 do [
test (t eq 1)
ifso [ s=px; p=x; p1=p1x; p2=p2x; p3=p3x ]
ifnot [ s=py; p=y; p1=p1y; p2=p2y; p3=p3y ]
for i=0 to n-1 do [
FLDI(ac1, 0); FLDI(ac2, 0)
FLDI(ac3, 0); FLDI(ac4, 0)
for j=0 to 3 do [
FAD(ac4, FML(FLD(ac0, a+2*j), p+2*j))
FAD(ac3, FML(FLD(ac0, a+2*(j+4)), p+2*j))
FAD(ac2, FML(FLD(ac0, a+2*(j+8)), p+2*j))
FAD(ac1, FML(FLD(ac0, a+2*(j+12)), p+2*j))
]
FST(ac4, p3+2*i)
FST(ac3, p2+2*i)
FST(ac2, p1+2*i)
if saveXY then test convertXY
ifso s!i=FTR(ac1)
ifnot FST(ac1, s+2*i)
p=p+2
]
]
resultis PSquit(true)
]
and PSinit(psv) = valof [
if PSzone eq 0 resultis PSerror(0)
PSvec=psv
Zero(PSvec, lPSVEC)
// new floating point work area
let FPworkNew= Allocate(PSzone, 4*numFPacs+1)
if FPworkNew eq 0 resultis PSerror(1)
PSvec>>PSVEC.FPworkSave=FPwork
PSvec>>PSVEC.FPworkNew=FPworkNew
FPworkNew!0=numFPacs
FPSetup(FPworkNew)
FLDI(zero,0); FLDI(one, 1); FLDI(two,2); FLDI(six,6)
resultis true
]
and PSallocate(location, m) = valof [
let b=Allocate(PSzone, m)
if b eq 0 resultis PSquit(PSerror(1))
@location=b
resultis b
]
and PSerror(errorCode,a1,a2,a3,a4) = valof [
(table[#77403; #1401]) ("PS.ERRORS", lv errorCode)
resultis 0
]
and PSquit(result) = valof [
if PSvec eq 0 resultis result
FPSetup(PSvec>>PSVEC.FPworkSave)
PSvec>>PSVEC.FPworkSave=0
for i=0 to lPSVEC-1 do if PSvec!i ne 0 then Free(PSzone, PSvec!i)
PSvec=0
resultis result
]