/* % % Polynomial Spline Package--Stub % % See the drivers for examples of use. The basic structure is as % follows: % % SetUpSpline: first called to indicate choice of basis etc. % % SetEvaluation: called to choose the evaluation type (default = DirectTM) % % AddControlPoint: adds a vertex to be interpolated to the list. % It automatically decides when to plot to the screen. % % This package supports several splines. You are encouraged to try both % Catmull-Rom and B-spline. % % % You will be doing two things in this part of the assignment: % * Implement two other interpolation techniques, namely % the 4-point (Lagrange) interpolant and cubic Bezier. % * Implement two polynomial evaluation schemes. % % % Implementing New Splines % % The system is laid out so that to define a spline curve technique, % you must: % (a) provide a 4*4 basis matrix for the spline family. % Put the basis matrices in Bases.c. % (b) indicate the degree and the number of points shared % by adjacent curve segments (see the splines definitions % below. % % If you do (a) and (b) properly, then implementing the cubic Lagrange % basis should require NO further changes to the code. I.e., two lines % of code in InitSplines and four lines for the matrix definition in % the file Bases.h. % % However, Bezier curves will require a small modification. We want % the Bezier curves to be piecewise continuous in the first derivative % across curve segments. As you will see, what this means is that if % A,B,C,D define Bezier curve segment i, then D,E,F,G will define % curve segment i+1, where E is colinear with the line segment CD % and moreover is the same distance from D as C is from D. You will % need to be able to generate such points for the Bezier case. % A probable place to put this logic is in AddControlPoint below. % % % Evaluation of Polynomials % % One evaluation technique is already implemented for you: DirectTM. % This is the slowest possible technique, in that to get the value on % a curve, it multiplies T*M*P, where T is the row vector of the % value of the parameter, M is a 4*4 basis matrix, and P is a % 4-element list of values in each of x and y. Run the drivers, % first making sure you invoke: % % % SetEvaluation(DirectTM) % % and study the result. % % You can improve things for each curve segment by precomputing % A = M*P. A is simply the coefficients of the resulting polynomial % in monomial form. % % Once you have A, you can compute T*A for each value of t in [0,1] % desired. Furthermore, it is also easy to implement Horner's once % you have A. % % You must add code to evaluate the polynomial using % premultipling T*P to get A (DirectMP) % Horner's rule (Horner) % % % Implement these things one at a time. They're all fairly easy to % do. Suggested order of implementation: % % * first add code to ComputeA (3 or 4 lines). This will give % you the coefficients of the polynomials to render. % * use these coefficients to implement DirectMP (3 lines) % * then use these coefficients to implement Horner (3 lines) % % % The driver "interp.c" will invoke routines from this module to track % your mouse. Try it immediately with the DirectTM evaluation. % Put your mouse in the graphics window and hit the left mouse button. % After a while, you should get a trajectory to follow your mouse % pressings, albeit slowly. % % The driver "interp1.c" runs standalone. It takes in a sequence of % x,y pairs and plots a trajectory through them. You will be using % this routine for benchmarks. % % I repeat: both of these drivers can be run immediately so that you % can get an idea of how things work. % */ #include #include "X.h" #include "Spline.h" #include "Bases.h" static spline splines[NumberOfSplines]; static int init = FALSE; static RealPoint2D last = {-1.0,-1.0}; static int col = 1; static Points pt = {{0,0}, {0,0}, {0,0}, {0,0}}; static int ptc = 0; static matrix basis; static int degree; static int overlap; static double dt = 1/32.0; static int colourStyle = Black; static int eval = DirectTM; static int whichSpline = CubicCatmullRom; static int checking = CheckOff; /* Assignment 3 */ static double checkval = 0.0; /* Assignment 3 */ static void CopyMatrix (s,d) matrix s,d; { register int i,j; for (i=0; i<4; i++) for (j=0; j<4; j++) d[i][j] = s[i][j]; } /* % % Initialise the supported piecewise polynomial curves. Only % Catmull-Rom and BSpline make any sense right now. Notice that you % need to define the degree and the number of points shared by % adjacent curve segments. The assumption here is that if % A,B,C,D define curve segment i, and if, for example, the "overlap" % is 2, then points C,D,E,F would define curve segment i+1. % */ static void InitSplines() { CopyMatrix(CatRom,splines[CubicCatmullRom].basis); CopyMatrix(BSpline, splines[CubicBSpline].basis); splines[CubicCatmullRom].degree = 3; splines[CubicCatmullRom].overlap = 3; splines[CubicBSpline].degree = 3; splines[CubicBSpline].overlap = 3; /* * Set up spline values for CubicLagrange and CubicBezier here. * I've started things off. */ CopyMatrix(Bezier,splines[CubicBezier].basis); CopyMatrix(Lagrange,splines[CubicLagrange].basis); init = TRUE; } void SetUpSpline(s, subdiv, style) /* OK */ int s; int subdiv; int style; { if (!init) InitSplines(); whichSpline = s; degree = splines[s].degree; overlap = splines[s].overlap; dt = 1.0/subdiv; colourStyle = style; CopyMatrix(splines[s].basis,basis); /* for compatiblity */ } /* SetUpSpline */ void SetEvaluation(t) /* OK */ int t; { if (t == TableLookup) { printf("Table Lookup not implemented. Using direct evaluation\n"); eval = DirectTM; } else eval = t; } /* SetEvaluation */ /* This is relevant to Assignment 3 only */ void ChooseChecking(c,v) /* OK */ int c; double v; { checking = c; checkval = v; } /* ChooseChecking */ static void drawcross(p, c) /* OK */ RealPoint2D p; int c; { int i,j; i = round(p.x); j = round(p.y); drawline(i-2,j, i+2,j, c); drawline(i,j-2, i,j+2, c); } /* drawcross */ static void drawbox(lx,ly, tx, ty, col) /* OK */ int lx, ly, tx,ty, col; { drawline(lx,ly,lx,ty, col); drawline(lx,ty,tx,ty, col); drawline(tx,ty,tx,ly, col); drawline(tx,ly,lx,ly, col); } /* drawbox */ RealPoint2D GetMouseButtonPush() /* OK */ { RealPoint2D p; int done = FALSE; int x,y; XButtonEvent *button; do { button = GetEvent(); switch (button->button) { case 1: case 2: x = button->x; y = button->y; drawbox(x-2, y-2, x+2, y+2, col); p.x = x; p.y = y; done = TRUE; break; case 3: p.x = p.y = -1.0; done = TRUE; break; } } while (!done); return(p); } /* GetMouseButtonPush */ /* * Draw a line segment between two points on the curve */ static void PlotPoint(px,py, c) /* OK */ double px,py; int c; { if (last.x > 0) drawline(round(last.x),round(last.y), round(px), round(py),c); last.x = px; last.y = py; } /* PlotPoint */ /* * Compute the value of the spline at a given t in [0,1], * for a given set of control points, and basis. This is slow. */ static void EvaluateSpline(TM, pt, c) /* OK */ double TM[4]; Points pt; int c; { double px, py; int i; px = 0.0; py = 0.0; for (i=0; i<4; i++) { px += TM[i]*pt[i].x; py += TM[i]*pt[i].y; } PlotPoint(px,py, c); } /* EvaluateSpline */ /* * Compute TM = [t^3 t^2 t 1]M, where M is the basis matrix * for the spline and t is a specific value in [0,1]. */ static void ComputeTM(t, basis, TM) /* OK */ double t; matrix basis; double TM[4]; { double sum; int i,j; double T[4]; T[3] = 1.0; T[2] = t; T[1] = t*t; T[0] = T[1]*t; for (i=0; i<4; i++) { sum = 0.0; for (j=0; j<4; j++) sum += T[j]*basis[j][i]; TM[i] = sum; } } /* ComputeTM */ /* * Compute A = MP, where M is the basis matrix and P is a set of * points. This gives us the coefficients of the curves in x,y. */ static void ComputeA(basis, P, A) matrix basis; Points P; Points A; { /* * Add code here to compute the A coefficients for the x and y * polynomials. The loop looks similar to ComputeTM except that * the vector is on the other side of the matrix. */ } /* ComputeA */ /* * Plot a segment of the spline. Assumes the basis has already been set up. */ static void PlotSpline(pt) Points pt; { double t; double TM[4], T[4]; Points A; int c; double px,py; double h, hh, hhh; /* dt, dt^2 and dt^3 */ RealPoint2D D, DD, DDD; /* forward differences */ if (eval != DirectTM) ComputeA(basis, pt, A); if (eval == ForwardDiff) { /* set up forward diffs */ /* * Set up initial values for forward differences here. * Optional in this assignment. */ } t = 0; while (t <= 1.0) { switch (eval) { case DirectTM: ComputeTM(t, basis, TM); EvaluateSpline(TM, pt, col); break; case DirectMP: /* * Do direct evaluation using T*A. * Use PlotPoint to write a line segment. */ break; case Horner: /* * Rearrange polynomial in Horner's form. * Use PlotPoint to write a line segment. */ break; case ForwardDiff: /* * Plot point and compute forward differences * here. Optional in this assignment. */ break; } t += dt; } if (colourStyle == MultiColour) col++; } /* PlotSpline */ /* * A debugging utility to print current control points. */ void PrintControlPoints() { int i; for (i=0; i<4; i++) { printf("%d: (%f,%f)\n", i,pt[i].x,pt[i].y); } } /* * * Add a control point to list. We buffer the number of points until * we exceed the degree of the curve. Then we plot the curve and reset * our array pt of control points to be ready for the next curve segment. * * Note: you may have to change the code here slightly to implement * continuous Bezier curves. */ void AddControlPoint( p ) RealPoint2D p; { int i; pt[ptc] = p; ptc++; if (ptc > degree) { /* Uncomment the following for debugging purposes. */ /* PrintControlPoints(); */ PlotSpline(pt); /* * The following code is hard to understand. Simulate it to * see how it works. The idea is this. We are done with * the current curve segment. So, we shift exactly overlap * points from the end of pt into the first overlap slots of pt. * Once we set pt up, we make ptc point to next slot in pt. */ for (i=0; i