(* :Title: Spline Fit *)

(* :Author: John M. Novak *)

(* :Summary:
This package introduces the SplineFunction object
and a fitting function to generate it, for
curve fitting with splines.
*)
	
(* :Context: NumericalMath`SplineFit` *)

(* :Package Version: 1.0 *)

(* :Copyright: Copyright 1992-1999, Wolfram Research, Inc. *)

(* :History: V1.0 by John M. Novak, July 1992 *)

(* :Keywords: splines, curve fitting *)

(* :Sources:
	Bartels, Beatty, and Barsky: An Introduction to
		Splines for Use in Computer Graphics and
		Geometric Modelling, Morgan Kaufmann, 1987.
	de Boor, Carl: A Practical Guide to Splines,
		Springer-Verlag, 1978.
	Levy, Silvio: Bezier.m: A Mathematica Package for
		Bezier Splines, December 1990.
*)

(* :Mathematica Version: 2.0 *)

(* :Limitation: does not currently handle 3D splines, although
	some spline primitives may produce a curve in space. *)

BeginPackage["NumericalMath`SplineFit`",
				"LinearAlgebra`Tridiagonal`"]

SplineFit::usage = 
"SplineFit[points, type] generates a SplineFunction object
of the given type from the points.  Supported types are Cubic,
CompositeBezier, and Bezier.";

SplineFunction::usage =
"SplineFunction[type, range, internal] represents a spline
created by the SplineFit function. It can accept an
argument in the given range for evaluation at some point
along the curve.  The argument is a parameter along the
curve, with 0 at the first point, 1 at the second, and so
on.";

Cubic::usage = "
Cubic is a type of spline used in SplineFit.  The default
description sets the second derivative at the endpoints to 0 and has unit
parameter spacing between knots.";

CompositeBezier::usage = "
CompositeBezier is a type of spline used in SplineFit.
The default is to use unit parameter spacing between the given knots.
Each set of three knots defines a segment, with the first and third
knots being interpolated (the third knot is also the first knot of the
next segment.)  Each segment is a third-order Bezier spline, by virtue of an
additional knot, which is determined by continuity conditions (C1 continuity
where the segments are joined) and the non-interpolated knot from the
following segment. It is placed between the second and third knots.
Note that this fails at the end;  here, if the number of points given is even,
then the last two points are reversed so that the final point is
interpolated and the next to last is a control point for the final
segment; if odd, then the final vertex is doubled.";

Bezier::usage = " 
Bezier is a type of spline used in SplineFit. Given n points,
a spline of degree n-1 is created.  The spline has unit parameter spacing."

Begin["`Private`"]

SplineFit::cbezlen =
"Points are needed to generate a spline.";

SplineFit[pts_List?(MatrixQ[#, NumberQ[N[#]]&]&),
	type_Symbol?(MemberQ[{Cubic, Bezier, CompositeBezier}, #]&)] :=
	SplineFunction[type, {0., N[Length[pts] - 1]},
		pts,
		splineinternal[pts,type]]

Format[SplineFunction[t_,r_, b__]] :=
	SequenceForm["SplineFunction[",t,", ", r,", <>]"]

SplineFunction::dmval =
"Input value `1` lies outside the domain of the spline
function.";

SplineFunction[type_, {min_, max_}, pts_,
		internal_][in_?(NumberQ[N[#]]&)] :=
	Module[{out},
		If[in < min || in > max,
			Message[SplineFunction::dmval, in];
			out = $Failed,
		(* else *)
			out = evalspline[
				 Which[in == max, Min[max, in],
					   in == min, Max[min, in],
					   True, in], type, pts, internal]
		];
		out/;out =!= $Failed
	]

(* the spline internal routines.  This is where the internal
	forms for various spline types are defined. *)

splineinternal[pts_List,Cubic] :=
		Transpose[Map[splinecoord,Transpose[pts]]]

splineinternal[pts_List,CompositeBezier] :=
	Module[{eqns, gpts = pts,ln = Length[pts],end},
		If[ln < 3 || OddQ[ln],
			Which[ln == 1, gpts = Flatten[Table[gpts,{4}]],
				ln == 2, gpts = {gpts[[1]],gpts[[1]],gpts[[2]],gpts[[2]]},
				OddQ[ln], AppendTo[gpts,Last[gpts]],
				True, Message[SplineFit::cbezlen];
					Return[InString[$Line]]]];
		end = Take[gpts,-4];
		gpts = Partition[Drop[gpts,-2],4,2];
		gpts = Apply[{#1,#2,#3 - (#4 - #3),#3}&,gpts,{1}];
		AppendTo[gpts,end];	
		Apply[Transpose[{#1,3(#2 - #1),
				3(#3 - 2 #2 + #1),#4 - 3 #3 + 3 #2 - #1}]&,
			gpts,{1}]
	]

splineinternal[pts_List,Bezier] :=
	Module[{n, eq, deg = Length[pts] - 1},
		eq = Table[#^n (1 - #)^(deg - n),{n,0,deg}];
		Function[Evaluate[Plus @@ (pts Table[Binomial[deg,n],
			{n,0,deg}] eq)]]
	]

(* some functions to assist the Cubic splineinternal routine *)
trisolve[lst_, ln_] :=
	Module[{a, b, c},
		a = c = Table[1, {ln - 1}];
		b = Join[{2}, Table[4, {ln - 2}], {2}];
		TridiagonalSolve[a, b, c, lst]
	]

splinecoord[vals_] := 
	Module[{lst,ln = Length[vals],d,n},
		lst = Join[{3 (vals[[2]] - vals[[1]])},
			Table[3 (vals[[n + 2]] - vals[[n]]),
					{n,ln - 2}],
			{3 (vals[[ln]] - vals[[ln - 1]])}];
		d = trisolve[lst, ln];
		Table[{vals[[n]],d[[n]],
			3(vals[[n+1]]-vals[[n]])-2 d[[n]]-d[[n+1]],
			2(vals[[n]]-vals[[n+1]])+d[[n]]+d[[n+1]]},
				{n,1,ln - 1}]]

(* routines to evaluate the spline function at particular
	values of the parameter *)

evalspline[pt_?(# == 0 &), Cubic, pts_, internal_] :=
	{internal[[1,1,1]], internal[[1,2,1]]}

evalspline[pt_, Cubic, pts_, internal_] :=
	Module[{tmp},
		({1, #, #^2, #^3}& @
			If[(tmp = Mod[pt,1]) == 0, 1, tmp]) .
		Transpose[internal[[Ceiling[pt]]]]
	]

evalspline[pt_, CompositeBezier, pts_, internal_] :=
	Module[{ln = Length[pts] - 1},
		evalspline[If[ln <= 3,
			pt/ln,
			pt (1/2 - If[OddQ[ln], 1/(2 ln) ,0])],
		 Cubic, pts, internal]
	]

evalspline[pt_, Bezier, pts_, internal_] :=
	internal @ (pt/(Length[pts] - 1))

End[]

EndPackage[]




