%
%	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.
%		(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.
%
%	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.
%
%	
%	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:
%
%
%		Spline.SetEvaluation(Spline.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.t" 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.t" 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.
%
%	Read this code from bottom to top to get a idea of the flow.  You
%	will see instructions below on what to add.  Very little modification
%	is required.  Procedures that do not require modification are flagged
%	"OK".
%

unit					% remove this for Turing at home
module Spline
import var Event in "%oot/lib/Event"	% ouch! Turing will have a problem here
export	SetUpSpline,
	SetEvaluation,
	AddControlPoint,
	GetMouseButtonPush,
	point,
	evaluationType,
	CubicCatmullRom, CubicBSpline, CubicBezier, CubicLagrange,
	DirectTM, DirectMP, Horner, ForwardDiff, TableLookup

Event.Set("<Btn1Up>:<Btn1Down>:<Btn1Motion>")

const	NumberOfSplines	:= 4
const	Evaluations	:= 5

type splineType: 1..NumberOfSplines

const	CubicCatmullRom	:splineType := 1
const	CubicBSpline	:splineType := 2
const	CubicBezier	:splineType := 3
const	CubicLagrange	:splineType := 4

const	MultiColour	:= 0
% const	Black		:= 1

type evaluationType: 1..Evaluations
const	DirectTM	:evaluationType := 1	% multiply TMP every time
const	DirectMP	:evaluationType := 2	% premultiply MP to get A
const	Horner		:evaluationType := 3	% use Horner's rule
const	ForwardDiff	:evaluationType := 4	% use forward differencing
const	TableLookup	:evaluationType := 5	% use table lookup

type point:
	record
		x,y: real;
	end record

type Points:
	array 0..3 of point

type matrix:
	array 0..3, 0..3 of real

type spline:
	record
		degree: int			% degree of spline
		overlap: int			% control points shared with
						% previous segment
		basis: matrix			% basis matrix
	end record


var splines: array 1..NumberOfSplines of spline	% records all spline info

%
%	Here are the spline matrix bases
%

const CatRom: matrix := init ( -0.5,  1.5, -1.5,  0.5,
			        1.0, -2.5,  2.0, -0.5,
			       -0.5,  0.0,  0.5,  0.0,
			        0.0,  1.0,  0.0,  0.0
			)


const BSpline: matrix := init (
	-0.16666667,	 0.5,		-0.5,		0.16666667,
	 0.5,		-1.0,	 	0.5,		0.0,
	-0.5,	 	 0.0,	 	0.5,		0.0,
	 0.16666667,	 0.66666667,	0.16666667,	0.0
	)

%
%	Define the bases for the Cubic Bezier and Cubic Lagrange forms here.
%
const Bezier: matrix := init(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
const Lagrange: matrix := init(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)



%
%	Some supported piecewise polynomial curves.  Only Catmull-Rom
%	BSpline and Linear 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.
%
splines(CubicCatmullRom).degree  := 3;
splines(CubicCatmullRom).overlap := 3;
splines(CubicCatmullRom).basis   := CatRom;

splines(CubicBSpline).degree  := 3;
splines(CubicBSpline).overlap := 3;
splines(CubicBSpline).basis   := BSpline;

%
%	Initialise entries for splines(CubicBezier) and
%	splines(CubicLagrange) here.  I've started it off for you.
%
splines(CubicLagrange).basis   := Lagrange;
splines(CubicBezier).basis     := Bezier;


%
%	Basic internal variables.  You don't need to touch these.
%
var last: point 	:= init(-1,-1)
var col: int		:= 1
var pt: Points		:= init( init(0,0), init(0,0), init(0,0), init(0,0))
var ptc: int		:= 0
var whichSpline: int	:= CubicCatmullRom
var basis: matrix
var degree: int
var overlap: int
var dt: real		:= 1/32.0
var colourStyle: int	:= Black
var eval: evaluationType:= DirectTM;


proc SetUpSpline(s: splineType, subdiv: int, style: int)		% OK
	degree      := splines(s).degree
	overlap     := splines(s).overlap
	basis       := splines(s).basis
	whichSpline := s
	colourStyle := style;
	dt	    := 1.0/subdiv	% number of line segs per curve
					% means our step size is 1 / that
end SetUpSpline



proc SetEvaluation(t: evaluationType)					% OK
	if t = TableLookup then
	    put "Table Lookup not implemented.  Using direct evaluation."
	    eval := ForwardDiff
	else
	    eval := t;
	end if
end SetEvaluation


proc drawcross(p: point, c: int)					% OK
	var i,j: int
	i := round(p.x)
	j := round(p.y)
	drawline(i-2,j, i+2,j, c)
	drawline(i,j-2, i,j+2, c)
end drawcross


function GetMouseButtonPush: point					% OK
	var fp, p: point
	var e: Event.Descriptor

	loop
		Event.Get(e)
		if(e.kind = Event.Kind.ButtonPress) then
			drawbox(e.x-2, e.y-2, e.x+2,e.y+2, col)
			p.x := e.x
			p.y := e.y
			exit
		end if
	end loop
	result p
end GetMouseButtonPush



%
%	Draw a line segment between two points on the curve
%
proc PlotPoint(px,py: real, c: int);					% OK

	if (last.x > 0) then
		drawline(round(last.x),round(last.y), round(px), round(py),c);
	end if

	last.x := px;
	last.y := py;

end PlotPoint


%
%	Evaluate the value of the spline at a given t in [0,1],
%	for a given set of control points, and basis.
%
proc EvaluateSpline(TM: array 0..3 of real, pt: Points, c: int)		% OK

	var px, py: real

	px := 0.0
	py := 0.0

	for i: 0..3
		px += TM(i)*pt(i).x;
		py += TM(i)*pt(i).y;
	end for
	PlotPoint(px,py, c)
end EvaluateSpline


%
%	Compute [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].
%
proc ComputeTM(t:real, basis: matrix, var TM: array 0..3 of real)	% OK
	var sum: real;
	var T: array 0..3 of real

	T(3) := 1.0;
	T(2) := t;
	T(1) := t*t;
	T(0) := T(1)*t;

	for i: 0..3
		sum := 0.0
		for j: 0..3
			sum += T(j)*basis(j,i)
		end for
		TM(i) := sum
	end for
end ComputeTM


%
%	Compute 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.
%
proc ComputeA(basis: matrix, P: Points, var A: Points)
	%
	%	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.
	%
end ComputeA


%
%	Plot a segment of the spline.  Assume basis is already set up.
%	You will need to add some code in this procedure.
%
proc PlotSpline(pt: Points)
	var t: real
	var TM, T: array 0..3 of real
	var A: Points
	var c: int
	var px,py: real

	if eval not= DirectTM then
		ComputeA(basis, pt, A)
	end if

	if eval = ForwardDiff then		% set up forward diffs
		%
		%	Set up initial values for forward differences here.
		%	Optional in this assignment.
		%
	end if

	t := 0

	loop
		exit when t > 1
		case eval of
		    label DirectTM:
			ComputeTM(t, basis, TM)
			EvaluateSpline(TM, pt, col)

		    label DirectMP:
			%
			%	Do direct evaluation using T*A.
			%	Use PlotPoint to write a line segment.
			%

		    label Horner:
			%
			%	Rearrange polynomial in Horner's form.
			%	Use PlotPoint to write a line segment.
			%

		    label ForwardDiff:
			%
			%	Plot point and compute forward differences
			%	here.  Optional in this assignment.
			%
		end case

		t += dt

	end loop

	if colourStyle = MultiColour then
		col := col mod 15 + 1
	end if
end PlotSpline


%
%	A local debugging utility that you may or may not find useful.
%
proc PrintControlPoints							% OK
	for i:0..3
		put i,":  (", pt(i).x, ",", pt(i).y, ")"
	end for
end PrintControlPoints


%
%	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.
%
proc AddControlPoint( p: point )
	assert (overlap > 0)
	pt(ptc) := p
	ptc += 1
	if ptc > degree then
		% PrintControlPoints		% uncomment for debugging
		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 into the first overlap slots of pt.
		% Once we set pt up, we make ptc point to next slot in pt.
		%
		for i: 0..overlap-1
			pt(i) := pt(degree-overlap+1 + i)
		end for
		ptc := overlap
	end if
end AddControlPoint


end Spline