
(* :Title: Spline *)

(* :Author: John M. Novak *)

(* :Summary:
This package introduces a Spline graphics primitive
and provides utilities for rendering splines.
*)
	
(* :Context: Graphics`Spline` *)

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

(* :Package Version: 2.0.2 *)

(* :History:
	V1.0 by John M. Novak, December 1990.
	V2.0 by John M. Novak, July 1992.
	V2.0.1 by John M. Novak, bug fix to handle DisplayString, February 1997.
	V2.0.2 by John M. Novak, remove default use of RenderSpline, February 1998.
*)

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

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

(* :Mathematica Version: 2.2 *)

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

BeginPackage["Graphics`Spline`",
	"NumericalMath`SplineFit`"]

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.";

SplinePoints::usage =
"SplinePoints is an option for Spline which determines the number of
initial samples for rendering a spline.  The default is 25.";

SplineDots::usage =
"SplineDots is an option for Spline that specifies a style for a point.
A Graphics style primitive, or list of such primitives will render the
point with that style.  None specifies that nothing should appear at each
control point, while Automatic causes a red dot to be placed at each point.";

SplineDivision::usage =
"SplineDivision is an option for Spline which specifies the maximum amount of
subdivision to be used in attempting to generate a smooth spline.";

RenderSpline::usage =
"RenderSpline[pts, type, splinefunction] is an alternate
rendering method for splines.  Rules may be 
attached to this symbol to override the default method that
produces a Line primitive.";

Begin["`Private`"]

Spline::sptt =
"Warning: Value of option SplinePoints -> `1` is not an
integer >= 3, setting to 25.";

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

Options[Spline] =
	{SplinePoints -> 25,
	SplineDots-> None,
	MaxBend -> 10.,
	SplineDivision -> 20.};

Spline[pts_List,type_Symbol,opts:(_Rule | _RuleDelayed)...] :=
	Spline[pts,type,
		SplineFit[pts,type],
		opts]
		
Spline[pts_List,type_Symbol,
		SplineFunction[t_,r_,cpts_,rest___],
		opts:(_Rule | _RuleDelayed)...] :=
	Spline[pts,type,SplineFit[pts,type],
		opts]/; pts != cpts || type != t

(* unfortunately, this definition must be disabled, because
   the PostScript primitive interacts poorly with bounds of
   plots, and line attributes such as thickness. The
   advanced user might wish to reenable this. *)
(*
RenderSpline[{{x1_, y1_}, {x2_, y2_},
			{x3_, y3_}, {x4_, y4_}},
		 Bezier, _] :=
	PostScript[
		"MBeginOrig",
		ToString[N[x1]]<>" "<>ToString[N[y1]]<>
			" moveto",
		ToString[N[x2]]<>" "<>ToString[N[y2]]<>
			" "<>ToString[N[x3]]<>" "<>ToString[N[y3]]<>
			" "<>ToString[N[x4]]<>" "<>ToString[N[y4]]<>
			" curveto",
		"stroke",
		"MEndOrig"
	]
*)

(* The following routines handle rendering of splines *)

Spline::args =
"The spline object `` is not a valid spline of the form
Spline[pts, type, (SplineFunction[...]), (opts)], and cannot be rendered.";

splinetoline[Spline[pts_List,type_Symbol,
		fn_SplineFunction,
		opts:(_Rule | _RuleDelayed)...]] :=
	Module[{res,dots,spts},
		{dots} = {SplineDots}/.{opts}/.
			Options[Spline];
		If[MemberQ[{Automatic, True}, dots],
			dots = {PointSize[.03], Hue[0]}];
		line = RenderSpline[pts, type, fn];
		If[Head[line] === RenderSpline,
			line = Line[splinepoints[pts,
				type, fn, opts]]
		];
		If[dots === None || dots === False,
			line,
			{Flatten[{dots,Map[Point,pts]}],
				line}
		]
	]

splinetoline[sp_] := (Message[Spline::args, sp]; {})

splinepoints[genpts_, type_,
		fn:SplineFunction[_,{min_, max_},___],
		 opts___] :=
	Module[{pts, ipts, rng, spts, sdiv, maxb},
		{spts, sdiv, maxb} = {SplinePoints,
			SplineDivision, MaxBend}/.{opts}/.
			Options[Spline];
		If[!IntegerQ[spts] || spts < 3,
			Message[Spline::sptt, spts];
			spts = 25];
		rng = Range[min, max,(max - min)/(spts-1)];
		ipts = Transpose[{rng, pts = Map[fn,rng]}];
		If[MemberQ[{0,0., Infinity, None, False}, sdiv],
			pts,
			Last[Transpose[sampleall[ipts,fn,maxb,
				(max - min)/(spts sdiv)]]]
		]
	]

splinepoints[args___] := (Message[Spline::args, Spline[args]]; {})

(* following adaptively samples the spline *)

bend[{u1_,pt1_}, {u2_,pt2_}, {u3_, pt3_}, min_] :=
	Module[{v1 = pt1 - pt2, v2 = pt3 - pt2,n},
		If[N[(n = norm[v1] norm[v2]) == 0 ||
				u2 - u1 <= min ||
				u3 - u2 <= min],
			0,
			Re[Pi - ArcCos[(v1 . v2)/n]]
		]
	]

norm[{x_, y_}] := Sqrt[x^2 + y^2]

sampleall[pts_, fun_, maxb_, mini_] :=
	Fold[sample[#1,#2,fun,maxb,mini]&,
		 Take[pts, 2], Drop[pts,2]]

sample[pts_, next_, function_, maxbend_, min_] :=
	With[{first = Drop[pts, -2], last = Sequence @@ Take[pts, -2]},
		If[N[bend[last, next, min] > maxbend Degree],
			Join[first,
				sampleall[interp[last, next, function],
							function, maxbend, min]
			],
		(* else *)
			Append[pts, next]
		]
	]

interp[pt1:{u1_,_},pt2:{u2_,_},pt3:{u3_,_},fun_] :=
	Module[{i1 = (u2 + u1)/2, i2 = (u3 +u2)/2,out},
		{pt1, {i1, fun[i1]}, pt2, {i2, fun[i2]}, pt3}
	]

(* following combines a spline into a polygon or line *)
splinesubsume[shape_] :=
	shape/.sp_Spline :> (Sequence @@ (splinepoints @@ sp))

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];

Unprotect[DisplayString];

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

Protect[DisplayString];

End[]

EndPackage[]
