// March 31, 1978 2:02 PM *** RESIDENT ***
get "zpDefs.bcpl"
// outgoing procedures:
external [
newSplineID
splineType
makeSpline
remakeSpline
computeSpline
]
// outgoing statics:
external [
@computeRange
@numRate
@denRate
@pseudoCyclic
@splineTension
]
static [
// spline computation parameters
@computeRange=64
@numRate=3
@denRate=2
// flag for pseudo cyclic hack vs periodic splines
@pseudoCyclic=false
@splineTension=false
]
// incoming procedures:
external [
MoveBlock // SYSTEM
SetBlock
Zero
Gets
typeForm // ZPUTIL
giveUp
ParametricSpline // PSPLINE
FLD; FST; FTR; FLDI // FLOAT
FNEG; FAD; FSB; FML
FDV; FCM; FSN
FSTDP; FLDDP
DPAD; DPSB
obtainBlock // ZPBLOCK
putBlock
checkPoint // ZPDRAW
startCode
endCode
checkSplineID
curve
splineColorSymbol
giveMeXY
]
// incoming statics:
external [
keys // SYSTEM
@splineTable // ZPINIT
@maxSplineID
// globals for encoding:
@newX // ZPDRAW
@newY
]
// local statics:
static [
@debugSpline=0
]
// local definitions
//*****************************************************************
// First: make spline block
//*****************************************************************
let makeSpline(n, xTable, yTable, brush, color, cyclic; numargs m) = valof [makeSpline
// xTable & yTable are FLOATING POINT
//returns id of new spline, 0 if error
unless n gr 0 resultis 0
if m ls 6 then cyclic=false
// check for false curves that actually are lines or dots
let minVal=vec 2
let maxVal=vec 2
if n gr 2 & TableConstant(n, yTable) then [
TableMax(n, xTable, maxVal)
TableMin(n, xTable, minVal)
n=2
FST(FLD(0, minVal), xTable)
FST(FLD(0, maxVal), xTable+2)
]
if n gr 2 & TableConstant(n, xTable) then [
TableMax(n, yTable, maxVal)
TableMin(n, yTable, minVal)
n=2
FST(FLD(0, minVal), yTable)
FST(FLD(0, maxVal), yTable+2)
]
if n eq 2 & TableConstant(2, xTable) & TableConstant(2, yTable) then n=1
// now make the spline
let splinePointer=obtainBlock(SPLINEknotBase+4*n)
unless splinePointer resultis giveUp("[makeSpline]")
Zero(splinePointer, SPLINEknotBase+4*n)
//make spline block
let knotTable=splinePointer+SPLINEknotBase
let dn=2*n
MoveBlock(knotTable, xTable, dn)
MoveBlock(knotTable+dn, yTable, dn)
splinePointer>>SPLINE.nKnots=n
splinePointer>>SPLINE.brush=brush
splinePointer>>SPLINE.color=color
splinePointer>>SPLINE.cyclic=cyclic
splineType(splinePointer)
giveMeXY(splinePointer, lv splinePointer>>SPLINE.xColor,
lv splinePointer>>SPLINE.yColor)
let id=remakeSpline(splinePointer)
unless id then putBlock(splinePointer)
resultis id
]makeSpline
and splineType(splinePointer) be [splineType
let n=splinePointer>>SPLINE.nKnots
let xTable=splinePointer+SPLINEknotBase
let yTable=xTable+2*n
let x1=FTR(FLD(0, xTable))
let y1=FTR(FLD(0, yTable))
let x2, y2=nil, nil
test n eq 1
ifso [ x2=x1; y2=y1 ]
ifnot [
x2=FTR(FLD(0, xTable+2))
y2=FTR(FLD(0, yTable+2))
]
splinePointer>>SPLINE.left= (x1 ls x2) ? x1, x2
splinePointer>>SPLINE.right= (x1 gr x2) ? x1, x2
splinePointer>>SPLINE.top= (y1 gr y2) ? y1, y2
splinePointer>>SPLINE.bottom= (y1 ls y2) ? y1, y2
splinePointer>>SPLINE.type = selecton n into [
case 1: dotSpline;
case 2: ((y1 eq y2) ? horSpline, ((x1 eq x2) ? verSpline, regSpline));
default: regSpline ]
]splineType
//*****************************************************************
// Second: get spline id
//*****************************************************************
and remakeSpline(splinePointer) = valof [remakeSpline
let id=newSplineID()
unless id resultis 0
test computeSpline(splinePointer)
ifso [
splineTable!0=splineTable!0+1
splineTable!id=splinePointer
resultis id
]
ifnot resultis 0
]remakeSpline
and newSplineID() = valof [newSplineID
for id=1 to maxSplineID do unless splineTable!id resultis id
typeForm(0, "Sorry, no room for more than ", 10, maxSplineID, 0, "curves*N",
0, "To get a larger work space for curves, start DRAW with switch /S (e.g.: DRAW ", 10, 2*maxSplineID, 0, "/S )*N")
resultis 0
]newSplineID
//*****************************************************************
// Third: compute [ & generate chain encoding ]
//*****************************************************************
and computeSpline(splinePointer) = valof [
let r=((splinePointer>>SPLINE.type eq regSpline) ?
computeRegularSpline(splinePointer), computeFalseSpline(splinePointer))
splineColorSymbol(splinePointer)
resultis r
]
and computeFalseSpline(splinePointer) = valof [computeFalseSpline
// i.e. single point, or horizontal or vertical line
splinePointer>>SPLINE.chain=0
splinePointer>>SPLINE.nBeads=0
splinePointer>>SPLINE.xStart=splinePointer>>SPLINE.left
splinePointer>>SPLINE.yStart=splinePointer>>SPLINE.top
curve(splinePointer, drawMode)
resultis true
]computeFalseSpline
and computeRegularSpline(splinePointer) = valof [computeRegularSpline
// returns 0 if no storage
let procedureName="[computeRegularSpline]"
let n=splinePointer>>SPLINE.nKnots
let dn=2*n
// 4 extra knots for pseudo-cyclic curve
let xn=(splinePointer>>SPLINE.cyclic) & pseudoCyclic ? 4, 0
let nt=n+xn
let dnt=2*nt
// table for derivatives
let dTable=obtainBlock(nt*12)
unless dTable resultis giveUp(procedureName)
let knotTable=splinePointer+SPLINEknotBase
let kTable=obtainBlock(2*dnt)
unless kTable resultis giveUp(procedureName, dTable)
let xTable=kTable+xn
let yTable=xTable+dnt
MoveBlock(xTable, knotTable, dn)
MoveBlock(yTable, knotTable+dn, dn)
if xn ne 0 then [
// extra knots
MoveBlock(kTable, xTable+dn-6, 4)
MoveBlock(xTable+dn, xTable+2, 4)
MoveBlock(yTable-4, yTable+dn-6, 4)
MoveBlock(yTable+dn, yTable+2, 4)
]
manifest [
//local floating Point registers
t0=0; t1=1; t2=2; t3=3
//caution: global floating Point registers
X=8; Y=9
d1X=10; d2X=11; d3X=12
d1Y=13; d2Y=14; d3Y=15
]
// get storage for chain encoding
// (try to estimate block size...)
let s=1
let kx=FTR(FLD(t0,xTable))
let ky=FTR(FLD(t1,yTable))
for i=1 to n-1 do [
let kx1=FTR(FLD(t0,xTable+2*i))
let ky1=FTR(FLD(t1,yTable+2*i))
s=s+1+max(abs(kx-kx1),abs(ky-ky1))/16
kx,ky=kx1,ky1
]
let nBeads=(((s*4)/3)/BEADsize)+1
let maxByteCount=16*nBeads
let chainSize=nBeads*(BEADsize+2)
let chainBlock=obtainBlock(chainSize+maxByteCount/2)
unless chainBlock resultis giveUp(procedureName, dTable, (xn ? kTable, 0))
let wTable=0
if splineTension then [
wTable=obtainBlock(nt)
if wTable ne 0 then [
SetBlock(wTable, 0, nt)
typeForm(0,"Input spline weigths now @", 8, wTable, 1, $*N)
Gets(keys)
]
]
// spline computation:
// caution: ParametricSpline(n,X,Y,X',X'',X''',Y',Y'',Y''')
let psDone=ParametricSpline(nt, kTable, kTable+dnt, dTable, dTable+dnt, dTable+2*dnt,
dTable+3*dnt, dTable+4*dnt, dTable+5*dnt,
((splinePointer>>SPLINE.cyclic & not pseudoCyclic) ? periodicSpline, naturalSpline),
wTable)
putBlock(wTable)
unless psDone resultis giveUp(procedureName, dTable, chainBlock, (xn ? kTable, 0))
newX=FTR(FLD(t0,xTable))
newY=FTR(FLD(t1,yTable))
splinePointer>>SPLINE.chain=chainBlock
chainBlock!(chainSize-1)=maxByteCount
splinePointer>>SPLINE.nBeads=nBeads
splinePointer>>SPLINE.xStart=newX
splinePointer>>SPLINE.yStart=newY
startCode(splinePointer,newX,newY)
let two=vec 2; FST(FLDI(t0,2),two)
let three=vec 2; FST(FLDI(t0,3),three)
let six=vec 2; FST(FLDI(t0,6),six)
for k=0 to n-2 do [
//knot k
let kX,kY=xTable+2*k,yTable+2*k
//stepping parameters
let m=max(abs(FTR(FLD(t0,kX))-FTR(FLD(t1,kX+2))),
abs(FTR(FLD(t2,kY))-FTR(FLD(t3,kY+2))))
m=(m*numRate)/denRate
let ni,r=1,m
if m gr computeRange then [ r=computeRange; ni=m/r; m=ni*r ]
//constants
let fni=vec 2; FST(FLDI(t0,ni), fni)
test m eq 0
ifso FLDI(t1,0)
ifnot FDV(FLDI(t1,1), FLDI(t0,m))
let delta=vec 2; FST(t1,delta)
let delta2=vec 2; FST(FML(t1,delta), delta2)
let delta3=vec 2; FST(FML(t1,delta), delta3)
//derivatives - floating point -
let kX1=dTable+xn+2*k
let kX2=kX1+dnt
let kX3=kX2+dnt
let kY1=kX3+dnt
let kY2=kY1+dnt
let kY3=kY2+dnt
// show values
if debugSpline then [
typeForm(0, "m=", 10, m, 0, ", delta=", -2, delta, 0, ", X, X', X'', X''' =*N", -2, kX)
for i=0 to 2 do typeForm(1, $|, -2, kX1+dnt*i)
typeForm(1, $*N, -2, kY)
for i=0 to 2 do typeForm(1, $|, -2, kY1+dnt*i)
Gets(keys)
typeForm(1, $*N)
]
//starting derivatives of subintervals - floating point -
let X1=vec 2; let X2=vec 2; let X3=vec 2
let Y1=vec 2; let Y2=vec 2; let Y3=vec 2
//start at knot k
FLD(X,kX); FLD(Y,kY)
X1!0=kX1!0; X1!1=kX1!1
Y1!0=kY1!0; Y1!1=kY1!1
X2!0=kX2!0; X2!1=kX2!1
Y2!0=kY2!0; Y2!1=kY2!1
X3!0=kX3!0; X3!1=kX3!1
Y3!0=kY3!0; Y3!1=kY3!1
//3rd derivatives & differences are constant
FML(FLD(d3X,delta3), X3)
FML(FLD(d3Y,delta3), Y3)
//forward difference computation in subinterval
for i=1 to ni do [
//floating point computation of initial values (blahh!)
FML(FLD(d2X,delta2),X2)
FML(FLD(d1X,delta),X1)
FAD(d1X,FDV(FLD(t0,d2X),two))
FAD(d1X,FDV(FLD(t0,d3X),six))
FML(FLD(d2Y,delta2),Y2)
FML(FLD(d1Y,delta),Y1)
FAD(d1Y,FDV(FLD(t0,d2Y),two))
FAD(d1Y,FDV(FLD(t0,d3Y),six))
FAD(d2X,d3X); FAD(d2Y,d3Y)
//double precision fixed point variables
let dpX=vec 2; let dpd1X=vec 2; let dpd2X=vec 2; let dpd3X=vec 2
let dpY=vec 2; let dpd1Y=vec 2; let dpd2Y=vec 2; let dpd3Y=vec 2
//initial values
FSTDP(X,dpX); FSTDP(d1X,dpd1X); FSTDP(d2X,dpd2X); FSTDP(d3X,dpd3X)
FSTDP(Y,dpY); FSTDP(d1Y,dpd1Y); FSTDP(d2Y,dpd2Y); FSTDP(d3Y,dpd3Y)
//compute
for j=1 to r do [
newX=DPAD(dpX,dpd1X)
newY=DPAD(dpY,dpd1Y)
DPAD(dpd1X,dpd2X); DPAD(dpd1Y,dpd2Y)
DPAD(dpd2X,dpd3X); DPAD(dpd2Y,dpd3Y)
checkPoint()
]
//now, next interval or next knot
test i eq ni
ifso [
//next knot
newX=FTR(FLD(t0,kX+2))
newY=FTR(FLD(t0,kY+2))
]
ifnot [
//starting point of next interval
FDV(FLDI(t1,i),fni)
FDV(FLD(t2,t1),two)
FDV(FLD(t3,t1),three)
FST(FAD(FML(FLD(t0,kX3),t1),kX2),X2)
FST(FAD(FML(FLD(t0,kY3),t1),kY2),Y2)
FST(FAD(FML(FAD(FML(FLD(t0,kX3),t2),kX2),t1),kX1),X1)
FST(FAD(FML(FAD(FML(FLD(t0,kY3),t2),kY2),t1),kY1),Y1)
newX=FTR(FAD(FML(FAD(FML(FAD(FML(FLD(X,kX3),t3),kX2),t2),kX1),t1),kX))
newY=FTR(FAD(FML(FAD(FML(FAD(FML(FLD(Y,kY3),t3),kY2),t2),kY1),t1),kY))
]
checkPoint()
]
]
//end of chain encoding
endCode(splinePointer)
//return storage [ chain block is trimmed by endCode() ]
putBlock(dTable)
putBlock(kTable)
resultis true
]computeRegularSpline
//*****************************************************************
// miscellaneous useful stuff
//*****************************************************************
and abs(i) = ((i gr 0) ? i, -i)
and max(i1,i2) = ((i1 ge i2) ? i1, i2)
and TableMax(n, fpTable, fpMax) be [
manifest [ t0=0 ]
FST(FLD(t0, fpTable), fpMax)
for k=1 to n-1 do if FCM(t0, fpTable+2*k) eq 1 then
FST(FLD(t0, fpTable+2*k), fpMax)
]
and TableMin(n, fpTable, fpMin) be [
manifest [ t0=0 ]
FST(FLD(t0, fpTable), fpMin)
for k=1 to n-1 do if FCM(t0, fpTable+2*k) eq -1 then
FST(FLD(t0, fpTable+2*k), fpMin)
]
and TableConstant(n, fpTable) = valof [
manifest [ t0=0 ]
FLD(t0, fpTable)
for k=1 to n-1 do if FCM(t0, fpTable+2*k) ne 0 then resultis false
resultis true
]