
(* :Title: Spline *)

(* :Author: John M. Novak *)

(* :Summary:  Introduces a Spline graphic object and
	a function for generating these objects, as well as
	utilities for rendering them. *)
	
(* :Context: Graphics`Spline` *)

(* :Package Version: 1.1 *)

(* :History: V1.0 by John M. Novak, Dec. 1990 *)

(* :Keywords: splines, curve fitting, graphics *)

(* :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, Dec. 1990.
*)

(* :Warning: adds definitions to the function Display. *)

(* :Mathematica Version: 2.0 *)

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

BeginPackage["Graphics`Spline`"]

Spline::usage = "Spline[points,type] represents a spline object of
kind type through (or controlled by) points.  It produces a spline
graphics primitive of the form Spline[points,type,control]. control
is information relevant to the particular kind of spline, describing
it completely.";
	
Cubic::usage = "
Cubis is a type of spline, which is used in Spline[].  Default
description sets second derivative at endpoints = 0 and unit
parameter spacing between knots.";

CompositeBezier::usage = "
CompositeBezier is a type of spline, which is used in Spline[].
Default description uses unit parameter spacing between sets of
vertices.  This is set to have C1 continuity between joins; thus,
in general, every other point is an interpolating control vertex,
from the first point.  Each segment is defined by four vertices,
the first and last being interpolated, the second being given, and
the third being determined by the continuity condition and derived from
the following noninterpolated control vertex.  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, which is used in Spline[]. In generating,
given n points, it creates a spline of degree n-1.  The spline has
unit parameter spacing."

SplinePoints::usage = "
SplinePoints is an option for Spline, which determines number of
points of interpolation between each control point of a spline.
Default is $SplinePoints.";

SplineDots::usage = "
SplineDots is an option for Spline that specifies a style for a point.
SplineDots=None specifies that nothing should appear at each control
point.  SplineDots=style causes a Point graphic primitive to be placed
at each control point. Default is $SplineDots.";

$SplinePoints::usage = "$SplinePoints gives the default value for 
SplinePoints Option.  The initial value is 15."

$SplineDots::usage = "Default value for SplineDots Option.  The
initial setting is None."
	
Begin["`Private`"]

Spline::cbezlen =
	"You need points to generate a spline!"

Format[Spline[p_,t_,b__]] :=
	SequenceForm["Spline[",p,",",t,",<>]"]

Options[Spline] := {SplinePoints -> $SplinePoints,
	SplineDots-> $SplineDots}

Spline[pts_List,type_Symbol,opts:(_Rule | _RuleDelayed)...] :=
	Spline[pts,type,pts,
		splineinternal[pts,type],
		opts]
		
Spline[pts_List,type_Symbol,cpts_List,
					opts:(_Rule | _RuleDelayed)...] :=
	Spline[pts,type,pts,
		splineinternal[pts,type],
		opts]/; pts != cpts

$SplinePoints = 15;

$SplineDots = None;

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

splineinternal[pts_List,Cubic] :=
	Module[{ln = Length[pts],mat,cpts = Transpose[pts]},
		mat = multmat[ln];
		Transpose[Map[splinecoord[#,mat]&,cpts]]
]

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[Spline::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},
		Table[Binomial[Length[pts] - 1,n],
			{n,0,Length[pts] - 1}]
	]

(* some functions to assist the Cubic splineinternal routine *)
ShiftRight[list_] := Drop[Prepend[list,0],-1]

ShiftLeft[list_] := Drop[Append[list,0],1]

multmat[n_] := Module[{id = IdentityMatrix[n],l,r},
	l = Map[ShiftLeft,id];
	r = Map[ShiftRight,id];
	id = l + r + 4 id;
	id[[1,1]] = 2;
	id[[n,n]] = 2;
	id]
		
splinecoord[vals_,mat_] := 
	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 = LinearSolve[mat,lst];
		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}]]

(* the splinepoints routines.  these routines covert a spline
	into a list of points to actually be connected in generating
	a graphic. *)

splinepoints[pts_,Cubic,internal_,res_] :=
	Module[{rn,drn},
		rn = Range[0,1,1/(res-1)];
		drn = Map[{1,#,#^2,#^3}&,rn];
		Flatten[Map[Transpose,Map[drn.#&,internal,{2}]],1]
	]

splinepoints[pts_,CompositeBezier,internal_,res_] :=
	splinepoints[pts,Cubic,internal,2 res]

splinepoints[pts_,Bezier,internal_,res_] :=
	Module[{rn,eq,n,deg = Length[pts] - 1},
		rn = Range[0,1,1/((res - 1) * deg)];
		eq = Table[#^n (1 - #)^(deg - n),{n,0,deg}];
		Map[Evaluate[Plus @@ (pts internal eq)]&,rn]
	]

(* The following routines handle rendering of splines *)

splinetoline[Spline[pts_List,type_Symbol,pts_List,internal_,
				opts:(_Rule | _RuleDelayed)...]] :=
	Module[{res,dots,spts},
		{res,dots} = {SplinePoints,SplineDots}/.{opts}/.
			Options[Spline];
		spts = splinepoints[pts,type,internal,res];
		If[dots === None,
			Line[spts],
			{Flatten[{dots,Map[Point,pts]}],Line[spts]}]
	]

splinesubsume[shape_] :=
	shape/.{Spline[p_,t_,p_,c_,___,SplinePoints->r_,___] :>
				Sequence @@ splinepoints[p,t,c,r],
			Spline[p_,t_,p_,c_,___] :>
				Sequence @@ splinepoints[p,t,c,SplinePoints/.
								Options[Spline]]}

Unprotect[Display];

Display[f_,gr_?(!FreeQ[#,Spline[___]]&),opts___] :=
	(Display[f,gr/.{p:(Polygon[{___,_Spline,___}] | 
							Line[{___,_Spline,___}]) :>
					splinesubsume[p],
				v:Spline[___]:>splinetoline[v]},
			opts];
			gr)

Protect[Display];

End[]

EndPackage[]

