// File: SPLINE2.BCPL // P.Baudelaire // December 5, 1977 4:45 PM // The algorithms implemented by these procedures are described in // "Spline Curve Techniques" // by P.Baudelaire, R.Flegal, & R.Sproull // Xerox Internal Report (May 1977) // Uses MICROCODE floating point routines // outgoing procedures: external [ ParametricSpline CubicSpline PSerror ] // outgoing statics: external [ PSzone ] static [ PSzone=0 // storage zone ] // incoming procedures: external [ FLD; FAD; FML; FST // microFLOAT (Alto floating point package) FLDI; FSB; FDV FCM; FNEG; FSN FPSetup Allocate // Alto SYSTEM Free Zero ] // incoming statics: external [ FPwork // microFLOAT (floating point registers) ] // local definitions: manifest [ naturalSpline=0 periodicSpline=1 // floating point registers: 1 to 4 ac1=1; ac2=2; ac3=3; ac4=4 // constants: zero=5; one=6; two=7; six=8 numFPacs=9 ] structure PSVEC [ FPworkSave word FPworkNew word fpx word fpy word h word a word b word c word r word s word ] manifest lPSVEC= size PSVEC/16 // local statics: static [ PSvec=0 ] let ParametricSpline(n, x, y, d1x, d2x, d3x, d1y, d2y, d3y, splineType, nil; numargs nargs) = valof [ // default arguments, get storage, check various things let tempVec=vec lPSVEC if PSinit(tempVec) eq 0 resultis 0 if n ls 0 then [ n=-n if ConvertToFP(n, lv x, lv y) eq 0 resultis 0 ] switchon nargs into [ case 9: splineType=naturalSpline case 10: case 11: if splineType ne naturalSpline & splineType ne periodicSpline resultis PSquit(PSerror(3)) if n ls 3 then splineType=naturalSpline endcase default: resultis PSquit(PSerror(4)) ] // compute parametrization (polygonal line approximation) let h=PSallocate(lv(PSvec>>PSVEC.h), 2*n) if h eq 0 resultis 0 for i=0 to n-2 do [ FSB(FLD(ac2, x+2*(i+1)), x+2*i) if FSN(ac2) eq -1 then FNEG(ac2) FSB(FLD(ac3, y+2*(i+1)), y+2*i) if FSN(ac3) eq -1 then FNEG(ac3) FDV((FCM(ac2, ac3) eq 1 ? ac3, ac2), two) FST(FAD(ac2, ac3), h+2*i) ] // now, compute the cubic splines x(h) & y(h) let done=ComputeCubicSpline(n, h, x, d1x, d2x, d3x, splineType) if done then done=ComputeCubicSpline(n, h, y, d1y, d2y, d3y, splineType) // reparametrize to the interval [0, 1] if done then for i=0 to n-2 do [ FLD(ac3, FLD(ac1, h+2*i)) FST(FML(FLD(ac2, d1x+2*i), ac1), d1x+2*i) FST(FML(FLD(ac2, d1y+2*i), ac1), d1y+2*i) FML(ac1, ac3) FST(FML(FLD(ac2, d2x+2*i), ac1), d2x+2*i) FST(FML(FLD(ac2, d2y+2*i), ac1), d2y+2*i) FML(ac1, ac3) FST(FML(FLD(ac2, d3x+2*i), ac1), d3x+2*i) FST(FML(FLD(ac2, d3y+2*i), ac1), d3y+2*i) ] resultis PSquit(done) ] and CubicSpline(n,x,y,d1y,d2y,d3y,splineType; numargs nargs) = valof [ // default arguments, get storage, check various things let tempVec=vec lPSVEC if PSinit(tempVec) eq 0 resultis 0 if n ls 0 then [ n=-n if ConvertToFP(n, lv x, lv y) eq 0 resultis 0 ] switchon nargs into [ case 6: splineType=naturalSpline case 7: if splineType ne naturalSpline & splineType ne periodicSpline resultis PSquit(PSerror(3)) if n ls 3 then splineType=naturalSpline endcase default: resultis PSquit(PSerror(4)) ] let h=PSallocate(lv(PSvec>>PSVEC.h), 2*n) if h eq 0 resultis 0 for i=0 to n-2 do FST(FSB(FLD(ac1, x+2*(i+1)), x+2*i), h+2*i) // now, compute let done=ComputeCubicSpline(n, h, y, d1y, d2y, d3y, splineType) if done eq 0 resultis 0 resultis PSquit() ] and ComputeCubicSpline(n, h, y, d1, d2, d3, splineType) = valof [ // IMPORTANT: assume all arguments are right !!!! // ONLY called by CubicSpline & ParametricSpline let c,r,s=0,0,0 if splineType eq periodicSpline then [ if FCM(FLD(ac2,y), y+2*(n-1)) ne 0 resultis PSquit(PSerror(2)) c=PSallocate(lv(PSvec>>PSVEC.c), 2*n) r=PSallocate(lv(PSvec>>PSVEC.r), 2*n) s=PSallocate(lv(PSvec>>PSVEC.s), 2*n) if (c eq 0) % (r eq 0) % (s eq 0) resultis 0 ] let a= PSallocate(lv(PSvec>>PSVEC.a), 2*n) let b= PSallocate(lv(PSvec>>PSVEC.b), 2*n) if (a eq 0) % (b eq 0) resultis 0 // check overlapping knots! for i=0 to n-1 do if FCM(zero, h+2*i) eq 0 then FST(one, h+2*i) test splineType eq naturalSpline ifso [ // a(0)=2*(h(0)+h(1)) FST(FML(FAD(FLD(ac1, h), h+2), two), a) // a(i)=2*(h(i)+h(i+1))-h(i)*h(i)/a(i-1) {i=1,2,...,n-3} for i=1 to n-3 do [ FML(FAD(FLD(ac4, h+2*i), h+2*(i+1)), two) FDV(FML(FLD(ac2,h+2*i), ac2), ac1) FLD(ac1, FST(FSB(ac4, ac2), a+2*i)) ] ] ifnot [ // a(0)=2*(h(0)+h(n-2)) FST(FML(FAD(FLD(ac1, h), h+2*(n-2)), two), a) // a(i)=2*(h(i)+h(i-1))-h(i-1)*h(i-1)/a(i-1) {i=1,2,...,n-3} for i=1 to n-3 do [ FML(FAD(FLD(ac4, h+2*i), h+2*(i-1)), two) FDV(FML(FLD(ac2,h+2*(i-1)), ac2), ac1) FLD(ac1, FST(FSB(ac4, ac2), a+2*i)) ] // c(0)=h(n-2) FST(FLD(ac1, h+2*(n-2)), c) // c(i)=-h(i-1)*c(i-1)/a(i-1) {i=1,2,...,n-3} // Notice: c(i-1) is ac1 for i=1 to n-3 do FST(FNEG(FDV(FML(ac1, h+2*(i-1)), a+2*(i-1))), c+2*i) ] computebc: if n ge 3 then test splineType eq naturalSpline ifso [ // first: D(i)=(y(i+1)-y(i))/h(i) {i=0,1,...,n-2} // and: b(i)=6*(D(i+1)-D(i)) {i=0,1,...,n-3} // Notice: D(i) is ac1, D(i+1) is ac2 FDV(FSB(FLD(ac1, y+2), y), h) for i=0 to n-3 do [ FDV(FSB(FLD(ac2, y+2*(i+2)), y+2*(i+1)), h+2*(i+1)) FST(FML(FSB(FLD(ac3, ac2), ac1), six), b+2*i) FLD(ac1, ac2) ] // then: b(i)=b(i)-h(i)*b(i-1)/a(i-1) {i=1,2,...,n-3} // Notice: b(i-1) is ac1 FLD(ac1, b) for i=1 to n-3 do [ FDV(FML(FLD(ac2, ac1), h+2*i), a+2*(i-1)) FST(FSB(FLD(ac1, b+2*i), ac2), b+i*2) ] ] ifnot [ // first: D(i)=(y(i+1)-y(i))/h(i) {i=0,1,...,n-2} // and: b(i)=6*(D(i)-D(i-1)) {i=1,2,...,n-3} // and: b(0)= 6*(D(0)-D(n-2)) // Notice: D(i-1) is ac1, D(i) is ac2 FDV(FSB(FLD(ac1, y), y+2*(n-2)), h+2*(n-2)) for i=0 to n-3 do [ FDV(FSB(FLD(ac2, y+2*(i+1)), y+2*i), h+2*i) FST(FML(FSB(FLD(ac3, ac2), ac1), six), b+2*i) FLD(ac1, ac2) ] // then: b(i)=b(i)-h(i-1)*b(i-1)/a(i-1) {i=1,2,...,n-3} // Notice: b(i-1) is ac1 FLD(ac1, b) for i=1 to n-3 do [ FDV(FML(FLD(ac2, ac1), h+2*(i-1)), a+2*(i-1)) FST(FSB(FLD(ac1, b+2*i), ac2), b+i*2) ] // r(n-2)=1 and s(n-2)=0 FST(one,r+2*(n-2)) FST(zero,s+2*(n-2)) // r(i)=-(h(i)*r(i+1)+c(i))/a(i) {i=n-3,...,1,0} // s(i)=(b(i)-h(i)*s(i+1))/a(i) {i=n-3,...,1,0} for i=n-3 to 0 by -1 do [ FAD(FML(FLD(ac1, r+2*(i+1)), h+2*i), c+2*i) FST(FDV(FNEG(ac1), a+2*i), r+2*i) FSB(FLD(ac1, b+2*i), FML(FLD(ac2, s+2*(i+1)), h+2*i)) FST(FDV(ac1, a+2*i), s+2*i) ] ] computed2: // COMPUTE SECOND DERIVATIVES test splineType eq naturalSpline ifso [ // d2(0)=d2(n-1)=0 FST(zero,d2); FST(zero,d2+2*(n-1)) // d2(i)=(b(i-1)-h(i)*d2(i+1))/a(i-1) {i=n-2,...,2,1} for i=n-2 to 1 by -1 do [ FML(FLD(ac2, d2+2*(i+1)), h+2*i) FST(FDV(FSB(FLD(ac1, b+2*(i-1)), ac2), a+2*(i-1)), d2+2*i) ] ] ifnot [ // ac1=(6*(D(n-2)-D(n-3))-h(n-2)*s(0)-h(n-3)*s(n-3)) FDV(FSB(FLD(ac1, y+2*(n-1)), y+2*(n-2)), h+2*(n-2)) FSB(ac1, FDV(FSB(FLD(ac2, y+2*(n-2)), y+2*(n-3)), h+2*(n-3))) FSB(FML(ac1, six), FML(FLD(ac3, s), h+2*(n-2))) FSB(ac1, FML(FLD(ac3, s+2*(n-3)), h+2*(n-3))) // ac2=(h(n-2)*r(0)+h(n-3)*r(n-3)+2*(h(n-3)+h(n-2))) FAD(FML(FLD(ac2, r), h+2*(n-2)), FML(FLD(ac3, r+2*(n-3)), h+2*(n-3))) FAD(ac2, FML(FAD(FLD(ac3, h+2*(n-3)), h+2*(n-2)), two)) // ac1=d2(n-2)=ac1/ac2 FST(FDV(ac1, ac2), d2+2*(n-2)) // d2(i)=r(i)*d2(n-2) + s(i) {i=0,1,2,...,n-3} for i=0 to n-3 do FST(FAD(FML(FLD(ac2, ac1), r+2*i), s+2*i), d2+2*i) // d2(n-1)=d2(0) FST(FLD(ac1, d2), d2+2*(n-1)) ] computed1d3: // COMPUTE FIRST & THIRD DERIVATIVES // d1(i)=(y(i+1)-y(i))/h(i) - h(i)*(2*d2(i)+d2(i+1))/6 // d3(i)=(d2(i+1)-d2(i))/h(i) {i=0,1,2,...,n-2} for i=0 to n-2 do [ FDV(FSB(FLD(ac1, y+2*(i+1)), y+2*i), h+2*i) FAD(FML(FLD(ac2, d2+2*i), two), d2+2*(i+1)) FST(FSB(ac1, FDV(FML(ac2, h+2*i), six)), d1+2*i) FST(FDV(FSB(FLD(ac1, d2+2*(i+1)), d2+2*i), h+2*i), d3+2*i) ] resultis true ] and ConvertToFP(n, intXloc, intYloc) = valof [ // convert coordinates from integer to floating point let fpx=PSallocate(lv(PSvec>>PSVEC.fpx), 2*n) let fpy=PSallocate(lv(PSvec>>PSVEC.fpy), 2*n) if (fpx eq 0) % (fpy eq 0) resultis PSquit(PSerror(1)) for i=0 to n-1 do [ FST(FLDI(ac1, (@intXloc)!i), fpx+2*i) FST(FLDI(ac1, (@intYloc)!i), fpy+2*i) ] @intXloc=fpx @intYloc=fpy resultis true ] and PSerror(errorCode,a1,a2,a3,a4) = valof [ (table[#77403; #1401]) ("PS.ERRORS", lv errorCode) resultis 0 ] 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 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 ]