(* ::Package:: *) (* :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-2007, 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. *) Message[General::newpkg, "Graphics`Spline`","Splines Package"] Quiet[ BeginPackage["Graphics`Spline`", "NumericalMath`SplineFit`"] , {General::obspkg, General::newpkg}] 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`"] issueObsoleteFunMessage[fun_, context_] := Message[General::obspkgfn, fun, context]; 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)...] := (issueObsoleteFunMessage[Spline,"Graphics`Spline`"]; Spline[pts,type, SplineFit[pts,type], opts]) Spline[pts_List,type_Symbol, SplineFunction[t_,r_,cpts_,rest___], opts:(_Rule | _RuleDelayed)...] := (issueObsoleteFunMessage[Spline,"Graphics`Spline`"]; 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 2D spline of the form \ Spline[pts, type, (SplineFunction[...]), (opts)], and cannot be rendered."; splinetoline[Spline[pts:{{_?NumericQ,_?NumericQ}..},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)) Typeset`MakeBoxes[sp:Spline[___], fmt_, Graphics] := Typeset`MakeBoxes[#,fmt, Graphics]& @@ {splinetoline[sp]} (* these rules need to be placed before the regular rules for Line/Polygon, or they'll never get hit; thus, we need a by-hand ordering... *) DownValues[Typeset`MakeBoxes] = {HoldPattern[Typeset`MakeBoxes[ Graphics`Spline`Private`p:Polygon[{___, _Spline, ___}], Graphics`Spline`Private`fmt_, Graphics]] :> (Typeset`MakeBoxes[#1, Graphics`Spline`Private`fmt, Graphics] & ) @@ {Graphics`Spline`Private`splinesubsume[ Graphics`Spline`Private`p]}, HoldPattern[Typeset`MakeBoxes[Graphics`Spline`Private`p: Line[{___, _Spline, ___}], Graphics`Spline`Private`fmt_, Graphics]] :> (Typeset`MakeBoxes[#1, Graphics`Spline`Private`fmt, Graphics] & ) @@ {Graphics`Spline`Private`splinesubsume[Graphics`Spline`Private`p]}} ~Join~ DownValues[Typeset`MakeBoxes]; 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[]