' ' -------------------------INTERPOLATION MINI GUIDE------------------------- ' ' This is not a tutorial on interpolation. It shows a very common algorithm used ' for interpolation which I made general purpose to compare several ' interpolation types. I added comments to explain some subtle aspects. ' ' I included some common types: Bézier, B-spline and Catmul-Rom. Some others I ' ran across are: parabolic, simple cubic, Hermite and beta with tension & bias. ' I completely derived the linear basis matrix myself (big deal). ' ' I also derived the Kochanek-Bartels from their 84 SIGGRAPH paper to provide ' global control of tension, bias, continuity (it's been 25 years since I've ' done matrix math so it took a few passes). NO NURBS, I haven't gotten in that ' deep. ' ' I am not a math major nor an expert on splines and interpolation. For the ' last 2-1/2 years I read comp.graphics.algorithms, did some searches and ' grabbed related things. In Jan. '96, 4 months ago, I decided to finally code ' something and found that it was rather simple, thanks mostly to Leon DeBoer's ' posting. ' ' This routine does interpolation on a 2D list in the array: IT(c,n) where c=1 ' for the X coordinate, c=2 for Y and "n" the vector number. The same type of ' interpolation is done on both X and Y coordinates though this is not required ' in general and some interesting effects may be possible by using different ' interpolation types for the dimensions - room for a little art. There are ' comments for going to more dimensions later. I set up this routine to study ' the interpolation types and debug the code with the goal of moving to key ' frame animation. ' ' There is a problem with terminology which I will now define. Since Laser ' (light show) Graphics is "vector" graphics, the coordinate sets (X,Y here) are ' called vectors. As this code is in my image editor, I use that term. I think ' I use both vector and control point in my comments. They are the same. This ' algorithm considers all input image vectors as control points, even for ' Hermite (once you understand the Hermite, you'll see why this is incorrect for ' it). For key frame animation, the input control points are called key frames ' or "keys". ' ' This is parametric interpolation, which means that the coordinate values are ' interpolated vs. a parameter, in this case J. Y is NOT interpolated as a ' function of X as we might think to do. ' ' Because some types need data points before and after an interval in order to ' calculate a curve, the first and last vectors need extra data. A common ' solution is simply to repeat the end vectors to provide the data. This is ' done here. This forces the curve to always hit the 1st and last vectors. ' This along with the way each type works can confuse what is going on. The ' extra values could be made changeable to gain more control of the curve at the ' end vectors. ' ' This is a general purpose routine, a single interpolation type won't need all ' the terms to be calculated. Any basis matrix entry of zero means that the ' corresponding term in the coef. equation can be removed. An entire zero row ' means that the corresponding power of T calculation can be dropped. This ' routine calculates the polynomial equation in matrix form: ' ' | n-1 | ' 1 | n | Steve Noskowicz ' P = --- * [T^3 T^2 T 1] * [basis] * | n+1 | Q10706@email.mot.com ' D | n+2 | May-13:Nov-25 1996 ' ' Entry conditions (if I remember everything): ' * Arrays subscripts start at 1, not zero. ' * The first vector is in IT(c,2). IT(c,1) is the extra n-1 data. ' * The last vector is in IT(c,P). ' * InterPts is the number of NEW points added "between" key points as T goes ' 0-->1. ' ' NOTE: ' For types which hit the control points (linear, cubic, Catmul-Rom), when T=0 ' (or 1) control vectors are being "re-calculated". For these types, time can ' be saved by eliminating those calculations and passing the original vectors ' through. Also, if T goes to 1, you are double calculating. ' ' IType is used for selecting the basis matrix. Spline1: IType=1 : GOTO Spline ' CatMul Spline2: IType=2 : GOTO Spline ' B-Spline Spline3: IType=3 : GOTO Spline ' Parabolic Spline4: IType=4 : GOTO Spline ' Bézier (cubic) Spline5: IType=5 : GOTO Spline ' Linear Interp Spline6: IType=6 : GOTO Spline ' Hermite Spline7: IType=7 : GOTO Spline ' Beta with Tension & Bias Spline8: IType=8 : GOTO Spline ' Simple Cubic Spline9: IType=9 : GOTO Spline ' Kochanek-Bartels, the Holy Grail. Spline10: IType=10 : GOTO Spline ' Bezier (quadratic) Spline: ' This section adds the extra data at the ends by repeating the end vectors. ' Making these accessible would allow customizing the end behavior. IT(1,1) = IT(1,2) ' Make the n-1 data for first X IT(2,1) = IT(2,2) ' for Y ' Look where your type starts & stops to see if you need to do this. IT(1,P+1) = IT(1,P) ' Make the n+1 for IT(2,P+1) = IT(2,P) ' the last vector IT(1,P+2) = IT(1,P) ' Make the n+2 for last vector IT(2,P+2) = IT(2,P) ' This is only necessary for some conditions. 'You may want different start/stop points (NStart, NEnd) from what I have here.. ' ' In this long section we select the desired basis matrix. IF IType=1 THEN ' CatMul Rom basis Matrix (hits every input vector) ' ' As T goes from 0 to 1 you go from n to n+1 ' n-1 n n+1 n+2 ' T^3 -1 +3 -3 +1 / ' T^2 +2 -5 4 -1 / ' T^1 -1 0 1 0 / 2 ' T^0 0 2 0 0 / ' DD!=2 'divisor This allows the Basis matrix to be all integers. Ft1=-1/DD! : Ft2=+3/DD! : Ft3=-3/DD! : Ft4=+1/DD! Fu1=+2/DD! : Fu2=-5/DD! : Fu3=+4/DD! : Fu4=-1/DD! Fv1=-1/DD! : Fv2=0 : Fv3=+1/DD! : Fv4=0 Fw1=0 : Fw2=+2/DD! : Fw3=0 : Fw4=0 NStart=2 ' Start at the first vector NEnd=P-1 ' End at next to last NStep=1 ' Step one at a time JS=0 ' Interpolate from zero ' ELSEIF IType=2 THEN ' B-Spline basis Matrix ' ' Doesn't hit any input vectors ' As T goes from 0 to 1 you go from near n to near n+1 ' n-1 n n+1 n+2 ' T^3 -1 +3 -3 +1 / ' T^2 +3 -6 3 0 / ' T^1 -3 0 3 0 / 6 ' T^0 1 4 1 0 / ' DD!=6 ' divisor Ft1=-1/DD! : Ft2=+3/DD! : Ft3=-3/DD! : Ft4=+1/DD! Fu1=+3/DD! : Fu2=-6/DD! : Fu3=+3/DD! : Fu4=0 Fv1=-3/DD! : Fv2=0 : Fv3=3/DD! : Fv4=0 Fw1=+1/DD! : Fw2=+4/DD! : Fw3=1/DD! : Fw4=0 NStart=2 ' Start at the first vector NEnd=P-1 ' End at next to last NStep=1 ' Step one at a time JS=0 ' Interpolate from zero ' ELSEIF IType=3 THEN ' Parabolic basis Matrix from Leon de Boer. ' ' Doesn't hit key vectors. As T goes from 0 to 1 you go from ' the midpoint of (n-1) & n to the midpoint n & (n+1). ' n-1 n n+1 n+2 ' T^3 0 0 0 0 / ' T^2 +1 -2 1 0 / ' T^1 -2 +2 0 0 / 2 ' T^0 1 1 0 0 / ' DD!=2 ' The matrix divisor Ft1=0 : Ft2=0 : Ft3=0 : Ft4=0 Fu1=+1/DD! : Fu2=-2/DD! : Fu3=+1/DD! : Fu4=0 Fv1=-2/DD! : Fv2=+2/DD! : Fv3=0 : Fv4=0 Fw1=+1/DD! : Fw2=+1/DD! : Fw3=0 : Fw4=0 NStart=2 ' Start at the first vector. NEnd=P ' End at next to last NStep=1 ' Step one at a time JS=0 ' Interpolate from zero ' ELSEIF IType=4 THEN ' Bézier ' ' Hits every third vector so STEP=3 & start@ second vector (#3) ' As T goes from 0 TO 1 you go from n-1 TO n+2 ' n-1 n n+1 n+2 ' T^3 -1 +3 -3 +1 / ' T^2 +3 -6 +3 +0 / ' T^1 -3 +3 +0 +0 / 1 ' T^0 +1 +0 +0 +0 / ' DD!=1 'divisor Ft1=-1 : Ft2=+3 : Ft3=-3 : Ft4=+1 Fu1=+3 : Fu2=-6 : Fu3=+3 : Fu4=0 Fv1=-3 : Fv2=+3 : Fv3=0 : Fv4=0 Fw1=+1 : Fw2=0 : Fw3=0 : Fw4=0 NStart=3 ' Start at the second vector (the first is its n-1) NEnd=P-2 ' End down two since it hits n+2 NStep=3 ' Step three at a time JS=0 ' Interpolate from zero ' ELSEIF IType=5 THEN ' Linear Interpolation ' ' As T goes from 0 to 1 you go from n to n+1 ' n-1 n n+1 n+2 ' T^3 0 0 0 0 / ' T^2 0 0 0 0 / ' T^1 0 -1 1 0 / 1 ' T^0 0 1 0 0 / ' DD!=1 ' divisor Ft1=0 : Ft2=0 : Ft3=0 : Ft4=0 Fu1=0 : Fu2=0 : Fu3=0 : Fu4=0 Fv1=0 : Fv2=-1 : Fv3=+1 : Fv4=0 Fw1=0 : Fw2=+1 : Fw3=0 : Fw4=0 NStart=2 ' Start at the first vector NEnd=P-1 ' End at next to last NStep=1 ' Step one at a time JS=0 ' Interpolate from zero ' ELSEIF IType=6 THEN ' Hermite (non-spline) From Llew Mason ' ' Bézier is similar and easier to use for graphics because the extra ' control points are relative to the adjacent vectors. ' I arranged the basis matrix columns to be similar to the Bézier. ' The two tangent vector controls are in the middle (R1 & R2) ' However, they must have values near zero, not near the adjacent vectors. ' As T goes from 0 TO 1 you go from n-1 to n+2 ' The R1 & R2 are the velocity from/to end vectors respectively. ' R1 R2 ' n-1 n n+1 n+2 ' T^3 2 +1 +1 -2 / The order of columns changed ' T^2 -3 -2 -1 +3 / from regular Hermite. ' T^1 0 1 0 0 / 1 ' T^0 +1 0 0 0 / ' DD!=1 ' divisor Ft1=+2 : Ft2=+1 : Ft3=+1 : Ft4=-2 Fu1=-3 : Fu2=-2 : Fu3=-1 : Fu4=+3 Fv1=0 : Fv2=+1 : Fv3=0 : Fv4=0 Fw1=+1 : Fw2=0 : Fw3=0 : Fw4=0 NStart=3 ' Start at the second vector (the first is its n-1) NEnd=P-2 ' End down two since it hits n+2 NStep=3 ' Step three at a time JS=0 ' Interpolate from zero ' ELSEIF IType=7 THEN ' Beta Spline with Tension & Bias from Llew Mason. ' ' As tension goes from 0 to big, attraction is toward control point n. ' As tension goes to about -6, curve "repels" from control point n. ' As Bias goes from 1 to infinity attraction moves from point n toward n-1. ' As bias goes from 1 to 0, attraction moves toward n+1. ' As T goes from 0 to 1 you go from near n to near n+1 ' n-1 n n+1 n+2 ' T^3 -2B^3 2(T+B^3+B^2+B) -2(T+B^2+B+1) 2 / ' T^2 +6B^3 -3(T+2B^3+2B^2) 3(T+2B^2) 0 / ' T^1 -6B^3 6(B^3-B) 6B 0 / T+2B^3+4B^2+4B+2 ' T^0 2B^3 T+4(B^2+B) 2 0 / ' For B=1 & T=0 this traces the B-Spline ' FB=1' Bias FT=0' Tension ' DD!=FT+2*FB^3+4*FB^2+4*FB+2 ' The matrix divisor for Beta Ft1=-2*FB^3/DD! : Ft2=+2*(FT+FB^3+FB^2+FB)/DD! : Ft3=-2*(FT+FB^2+FB+1)/DD! : Ft4=2/DD! Fu1=+6*FB^3/DD! : Fu2=-3*(FT+2*FB^3+2*FB^2)/DD! : Fu3=+3*(FT+2*FB^2)/DD! : Fu4=0 Fv1=-6*FB^3/DD! : Fv2=+6*(FB^3-FB)/DD! : Fv3=+6*FB/DD! : Fv4=0 Fw1=+2*FB^3/DD! : Fw2=(FT+4*(FB^2+FB))/DD! : Fw3=+2/DD! : Fw4=0 NStart=2 ' Start at the first vector NEnd=P-1 ' End at next to last NStep=1 ' Step one at a time JS=0 ' Interpolate from zero ' ' ELSEIF IType=8 THEN ' Simple Cubic ' ' Derived from code found on the net written by Toby Orloff & Jim Larson ' U of Minn Geometry Supercomputer Project "omni_interp" source code. ' The derivative=0 at control points which gives straight ' lines between control points,. ' Hits every vector, but slows down at them. ' As T goes from 0 TO 1 you go from n TO n+1 ' n-1 n n+1 n+2 ' T^3 0 +2 -2 0 / ' T^2 0 -3 +3 0 / ' T^1 0 0 0 0 / 1 ' T^0 0 +1 0 0 / ' DD!=1 ' divisor Ft1=0 : Ft2=+2 : Ft3=-2 : Ft4=0 Fu1=0 : Fu2=-3 : Fu3=+3 : Fu4=0 Fv1=0 : Fv2=0 : Fv3=0 : Fv4=0 Fw1=0 : Fw2=+1 : Fw3=0 : Fw4=0 NStart=2 ' Start at the first vector NEnd=P-1 ' End at next to last NStep=1 ' Step one at a time JS=0 ' Interpolate from zero ' ELSEIF IType=9 THEN ' Kochanek-Bartels ' ' Basis matrix for global Tension, Continuity & Bias ' T H I S I S I T I figured it out from the SIGGRAPH paper ! ' As T goes from 0 to 1 you go from n to n+1 ' n-1 n n+1 n+2 ' T^3 -A 4+A-B-C -4+B+C-D D / ' T^2 +2A -6-2A+2B+C 6-2B-C+D -D / ' T^1 -A A-B B 0 / 2 ' T^0 0 2 0 0 / ' FT=0 ' Tension T=+1-->Tight T=-1--> Round FB=0 ' Bias B=+1-->Post Shoot B=-1--> Pre shoot FC=-1 ' Continuity C=+1-->Inverted corners C=-1--> Box corners ' ' When T=B=C=0 this is the Catmul-Rom. ' When T=1 & B=C=0 this is the Simple Cubic. ' When T=B=0 & C=-1 this is the linear interp. ' ' Here, the three parameters are folded into the basis matrix. If you want ' independent control of the segment start and end, make two T, C & Bs. ' One for A & B (beginning) and one for C & D (end). For local control of ' each point, you'll need an array of T, C & Bs for each individual point. ' If you want the local control as shown on the video or in the paper, you use ' the "A" & "B" for the current segment and the "C" & "D" from the previous ' segment. ' FFA=(1-FT)*(1+FC)*(1+FB) FFB=(1-FT)*(1-FC)*(1-FB) FFC=(1-FT)*(1-FC)*(1+FB) FFD=(1-FT)*(1+FC)*(1-FB) ' Here, the basis matrix coefficients are defined DD!=2 ' divisor for K-B Ft1=-FFA/DD! :Ft2=(+4+FFA-FFB-FFC)/DD! :Ft3=(-4+FFB+FFC-FFD)/DD! : Ft4=FFD/DD! Fu1=+2*FFA/DD! :Fu2=(-6-2*FFA+2*FFB+FFC)/DD! :Fu3=(+6-2*FFB-FFC+FFD)/DD! : Fu4=-FFD/DD! Fv1=-FFA/DD! :Fv2=(FFA-FFB)/DD! :Fv3=FFB/DD! : Fv4=0 Fw1=0 :Fw2=+2/DD! :Fw3=0 : Fw4=0 NStart=2 ' Start at the first vector NEnd=P-1 ' End at next to last NStep=1 ' Step one at a time JS=0 ' Interpolate from zero ELSEIF IType=9 THEN ' Bezier (quadratic) ' ' It uses the second degree Berstein basis function. ' There is ONE control point between interpolated points. ' The end velocities are both determined by the end-to-center vector. ' As T goes from 0 to 1 you go from n to n+2 ' n-1 n n+1 n+2 ' T^3 +0 -0 +0 +0 / ' T^2 -0 +1 -2 +1 / ' T^1 +0 -2 +2 -0 / 1 ' T^0 +0 +1 +0 +0 / DD!=1 ' divisor Ft1=0 : FT2=0 : Ft3=0 : Ft4=0 Fu1=0 : Fu2=+1 : Fu3=-2 : Fu4=+1 Fv1=0 : Fv2=-2 : Fv3=+2 : Fv4=0 Fw1=0 : Fw2=+1 : Fw3=0 : Fw4=0 NStart=2 ' Start at the first vector NEnd=P-2 ' End at next to last NStep=2 ' Step one at a time JS=0 ' Interpolate from zero END IF SplineDo: ' Finally, this is where the rubber meets the road. ' Derived from source code posted by Leon de Boer. ' delta2!=delta! * delta! ' Pre-compute delta squared and cubed delta3!=delta2! * delta! ' These two lines for forward differences only. ' ' Step through the vectors FOR n = NStart TO NEnd STEP NStep ' ' These are the coefficients a, b, c, d, in aT^3 + bT2 + cT + d ' NOTE: All terms are here, so all types use same code. ' To conserve calculations, terms with a zero in the basis matrix can be ' removed. FAX = Ft1*IT(1,n-1) + Ft2*IT(1,n) + Ft3*IT(1,n+1) + Ft4*IT(1,n+2) FBX = Fu1*IT(1,n-1) + Fu2*IT(1,n) + Fu3*IT(1,n+1) + Fu4*IT(1,n+2) FCX = Fv1*IT(1,n-1) + Fv2*IT(1,n) + Fv3*IT(1,n+1) + Fv4*IT(1,n+2) FDX = Fw1*IT(1,n-1) + Fw2*IT(1,n) + Fw3*IT(1,n+1) + Fw4*IT(1,n+2) ' FAY = Ft1*IT(2,n-1) + Ft2*IT(2,n) + Ft3*IT(2,n+1) + Ft4*IT(2,n+2) FBY = Fu1*IT(2,n-1) + Fu2*IT(2,n) + Fu3*IT(2,n+1) + Fu4*IT(2,n+2) FCY = Fv1*IT(2,n-1) + Fv2*IT(2,n) + Fv3*IT(2,n+1) + Fv4*IT(2,n+2) FDY = Fw1*IT(2,n-1) + Fw2*IT(2,n) + Fw3*IT(2,n+1) + Fw4*IT(2,n+2) ' For 3D you will have four more equations here for Z ' ' Calculate one segment of the interpolation here. ' NOTE: For types which pass through the keys, J can stop at InterPts ' Then just pass the key through (and the first one). WINDOW OUTPUT 2 ' Drawing window ' ' ================= Start of replacable section =========================== ' FX=FDX ' The first point in the segment FY=FDY ' And another for "Z" ' ' Start drawing at the very first point IF n = NStart THEN PSET(FX,255-FY) ' Accent the beginning of eachsegment (knot) LINE-STEP(2,2):LINE-STEP(-4,0):LINE-STEP(2,-2) ' Note, My display is 0-255 high & wide & the origin is in the lower left. FOR J = JS+1 TO InterPts+1 ' Interpolate one segment T=J/(InterPts+1) ' The interpolation parameter steps from delta to ' 1.0 ' The simple polynomial evaluation first pre-computes T^2 & T^3 ' T2 = T * T : T3 = T2 * T ' Then the polynomials ' FX = FAX*T3 + FBX*T2 + FCX*T + FDX ' FY = FAY*T3 + FBY*T2 + FCY*T + FDY ' And another for "Z" ' ' However, using Horner's Rule saves calculations & time. FX = ((FAX*T + FBX)*T + FCX)*T +FDX 'interpolated x value (using Horner) FY = ((FAY*T + FBY)*T + FCY)*T +FDY 'interpolated y value ' And another for "Z" ' ================ End of replacable section ========================= ' IF ILines=1 THEN LINE-(FX,255-FY) ELSE PSET(FX,255-FY) 'Draw lines/dots NEXT ' J sub-step "between" vectors NEXT ' n Move to the next vector set. ' Accent the very last point LINE-STEP(2,2):LINE-STEP(-4,0):LINE-STEP(2,-2) SplineEnd: RETURN ' Interpolation complete. Back to Kansas ' ' ============== Forward Difference version ===================== ' This goes between the two === lines above. ' Derived from Foley, van Damm, Feiner & Hughes ' Forward Difference "derivatives" calculated here. ' FX=FDX : FY=FDY ' Initial position ' And another for "Z" ' FD1X=FAX*delta3!+FBX*delta2!+FCX*delta! ' Initial velocity FD1Y=FAY*delta3!+FBY*delta2!+FCY*delta! ' And another for "Z" ' FD2X=6*FAX*delta3!+2*FBX*delta2! ' Initial acceleration FD2Y=6*FAY*delta3!+2*FBY*delta2! ' And another for "Z" ' FD3X=6*FAX*delta3! ' Acceleration change FD3Y=6*FAY*delta3! ' And another for "Z" ' ' Start drawing at the very first point IF n = NStart THEN PSET(FX,255-FY) ' ' Accent the beginning of each segment (knot) LINE-STEP(2,2):LINE-STEP(-4,0):LINE-STEP(2,-2) ' FOR J = JS TO InterPts ' Interpolate one segment FX = FX+FD1X 'interpolated x value FY = FY+FD1Y 'interpolated y value ' And another for "Z" ' FD1X=FD1X+FD2X : FD1Y=FD1Y+FD2Y ' Next speed ' And another for "Z" ' FD2X=FD2X+FD3X : FD2Y=FD2Y+FD3Y ' Next acceleration ' And another for "Z" ' '================END Forward Difference version =================