(*
From: levy@math.Princeton.EDU
Date: Thu, 21 Dec 89 13:54:47 -0500
Subject: spline package
Cc: levy@math.Princeton.EDU


I've written a package to draw splines.  It implements John Hobby's
algorithm for finding control points, which I've found to produce very
good results.  It is the algorithm implemented in METAFONT.

The package uses another package I wrote, to solve tridiagonal
(and tridiagonal-plus-corners) matrices.

It also uses ParametricPlot to draw the actual cubic splines.  This takes
the greatest part of the processing time if we have a large number of
curves with few patches.  If we have curves with many patches, on the
other hand, the repeated application of rules in AssignDefaults and
AssignControls makes the time grow quadratically with the number of
patches, which is certainly a nuisance.  Perhaps I'll change the
programming to depend less on pattern matching, and perhaps I won't;
right now I'm drawing mostly short curves.

The Spline.m package doesn't have a "private" environment because I'm
not sure which of the internal constructs should be visible outside.

	SilvioLevy
	
*)
BeginPackage["Spline`", "Tridiagonal`"]

(* 
This package by Silvio Levy is an implementation of certain parts
of Knuth's METAFONT program, concerning the drawing of curves by means
of cubic splines.  The main reason for following METAFONT is its
algorithm for assigning control points when the user does not do so,
and the versatility with which the user can provide partial
information toward the assignment of control points.  The algorithm
itself is due to John Hobby (Discrete and Computational Mathematics 1
(1986), 127-144).  The reader is also referred to Knuth's ``METAFONT: The
Program'' (part D of ``Computers and Typesetting'').
*)

Path::usage = "\n
  Path[...] describes a path `a la METAFONT.  The arguments are knots\n
  (points the path goes through), together with optional expressions\n
  having head Tension, Direction, Curl, Controls and Cycle."

Direction::usage = "\n
  Direction[<angle>] or Direction[<pair>] is an argument to Path\n
  indicating the direction of a path as it goes into a knot (if placed\n
 just before the knot) or as it goes out (if placed just after).\n
  Direction and Curl are mutually exclusive.  If a single Direction\n
  or Curl flanks a knot it is used for both sides."

Curl::usage = "\n
  Curl[<positive number>] is an argument to Path indicating the\n
  ratio between curvatures at the endpoints of a patch (how much the\n
  patch ``curls'').  It applies to the curliness of a path as it goes\n
  into a knot (if placed just before the knot) or as it goes out (if\n
  placed just after).  Direction and Curl are mutually exclusive.\n
  If a single Direction or Curl flanks a knot it is used for both sides."

Tension::usage = "\n
  Tension[<arg>,<arg>] is an argument to Path indicating how flat a\n
  patch of a Bezier spline should be in the middle.  <arg> should\n
  be a number >3/4, or an expression of the form AtLeast[<arg>], where\n
  <arg> is a number >3/4, in which case the tension is increased, if\n
  necessary, until the curve lies within the sector determined by\n
  the direction at the departure point and the line joining the two\n
  endpoints. Tension[x] is equivalent to Tension[x,x]. <arg> defaults to 1."

Controls::usage = "\n
 Controls[<pair>,<pair>] is an argument to Path giving explicit\n
  control points for a patch of a Bezier spline.  It goes between\n
  the two corresponding knots, and excludes of Curl, Tension or Direction."

Cycle::usage = "\n
  Cycle is an argument to Path indicating that the path is a loop.\n
  It should be the last argument."

AssignDefaults::usage = "\n
  AssignDefaults[Path[...]] puts the Path in an internal form suitable\n
  for use with AssignControls.  To look at the internal form, use\n
  FullForm; by default it prints as if it kept the original form."

AssignControls::usage = "\n
  AssignControls[Path[...]] assigns control points to a Path, putting\n
  it into a form suitable for rendering with Render."

Render::usage = "\n
  Render[Path[...]] creates a list, suitable for inclusion in a\n
  Graphics object, containing a representation of the path."

(* SYNTAX *)
[K
(*
After the input conversion a path should have the following syntax:

   <path> = OpenPath[<first_knot>,<knot>...,<last_knot>]
          | CyclePath[<knot>..]

   <knot> = Knot[<dir_info>,<position>,<dir_info>]
   <first_knot> = Knot[EndPoint,<position>,<dir_info>]
   <last_knot> = Knot[<dir_info>,<position>,EndPoint]

   <position> = {_?NumberQ,_?NumberQ}

   <dir_info> = Control{_?NumberQ,_?NumberQ}    (* Knuth's Explicit *)
              | EndPoint
              | Curl[_?CurlyQ,_?TensionQ]
              | Direction[_?NumberQ,_?TensionQ] (* Knuth's Given *)
              | Open[_?TensionQ]

In addition, the following restrictions apply:
  (i) Open can only appear in the same Knot with Open, Control or EndPoint;
 (ii) The right argument of a Knot has head Control if and only if the
     left argument of the next Knot has head Control.
*)

CurlyQ[x_] := NumberQ[x] && x>0
TensionQ[x_] := NumberQ[x] && (x>3/4 || x<-3/4) (* negative = AtLeast *)
PathQ[x_] := SameQ[x,OpenPath] || SameQ[x,CyclePath]
NotControlQ[x_] := !SameQ[Head[x],Control]
OpenQ[Knot[x_,__,y_]] := SameQ[Head[x],Open] && SameQ[Head[y],Open]

GoodPathQ[x_OpenPath] :=
  GoodFirstKnotQ[First[x]] && And @@ GoodKnotQ /@ Take[x,{2,-2}] &&
  GoodLastKnotQ[Last[x]] &&
  And @@ GoodPatchQ /@ Transpose[{List @@ Drop[x,-1],List @@ Drop[x,1]}]
GoodPathQ[x_CyclePath] :=
  And @@ GoodKnotQ /@ x && 
  And @@ GoodPatchQ /@ Transpose[{List @@ x,List @@ RotateLeft[x]}]
GoodPathQ[___] := False

GoodKnotQ[xxx:Knot[xx:(x_[__]),{_?NumberQ,_?NumberQ},yy:(y_[__])]] :=
  GoodSpecsQ[xx] && GoodSpecsQ[yy] &&
  If[(SameQ[x,Open])===SameQ[y,Open] || SameQ[x,Control] || SameQ[y,Control],
    True,Message[Knot::badopen,xxx];False]
GoodKnotQ[x___] := (Message[Knot::garbled,"",x];False)

GoodFirstKnotQ[Knot[EndPoint,{_?NumberQ,_?NumberQ},x_]] := GoodSpecsQ[x]
GoodFirstKnotQ[x___] := (Message[Knot::garbled,"first ",x];False)

GoodLastKnotQ[Knot[x_,{_?NumberQ,_?NumberQ},EndPoint]] := GoodSpecsQ[x]
GoodLastKnotQ[x___] := (Message[Knot::garbled,"last ",x];False)

GoodSpecsQ[Control[x___]] := 
  If[MatchQ[{x},{{_,_}}] && NumberQ[First[x]] && NumberQ[Last[x]],
    True,Message[Knot::badspec,Control[x]];False]
GoodSpecsQ[Curl[x___]] :=
  If[MatchQ[{x},{_,_}] && CurlyQ[First[{x}]] && TensionQ[Last[{x}]],
    True,Message[Knot::badspec,Curl[x]];False]
GoodSpecsQ[Direction[x___]] :=
  If[MatchQ[{x},{_,_}] && NumberQ[First[{x}]] && TensionQ[Last[{x}]],
    True,Message[Knot::badspec,Direction[x]];False]
GoodSpecsQ[Open[x___]] :=
  If[MatchQ[{x},{_}] && TensionQ[x],
    True,Message[Knot::badspec,Open[x]];False]
GoodSpecsQ[x___] := (Message[Knot::badkeyword,x];False)

GoodPatchQ[{xx:Knot[_,_,x_],yy:Knot[y_,_,_]}] := 
  If[SameQ[SameQ[Head[x],Control],SameQ[Head[y],Control]],True,
    Message[Knot::badcontrols,xx,yy];False]

Knot::badcontrols = "Illegal consecutive knots ``and ``"
Knot::badspec = "Illegal specification ``"
Knot::badkeyword = "Illegal keyword in ``"
Knot::badopen = "Illegal Open in ``"
Knot::garbled = "Garbled `` knot ``"

(* INPUT *)

(*
The full specification of a path, knot by knot, is very cumbersome.
In general we want to specify only the positions for the knots,
and perhaps a few directions and curls and tensions -- very rarely
explicit control points.

METAFONT uses a number of rules to specify default values for
specifications not given.  The application of these rules is carried
out in the body of the function AssignDefaults.

In general the input form for a path is very similar to what it would
be in METAFONT, with the following correspondence:

        METAFONT syntax                 Spline.m syntax

        ..                              ,  (argument separator)
        ...                             Dots
        --                              Dash
        ---                             DDash
        {x,y} or {(x,y)}                Direction[x,y] or Direction[Arg[x+Iy]]
        {curl x}                        Curl[x]
        tension x and atleast y         Tension[x,AtLeast[y]]
        controls x and y                Controls[x,y] 
        cycle                           Cycle

If a Curl or Direction is the only directive between two knots, two
commas (that is, a Null argument) must be used to separate it from the
knot that it does not affect. For example, the METAFONT path specification

z0{dir 20}..z1..tension atleast1..{curl 2}z2..z3{-1,-2}..tension 3 and 4..
z4..controls z45 and z54..z5..cycle

can be input as

Path[z0,Dir[20 Degree],,z1,Tension[AtLeast[1]],Curl[2],z2,z3,Direction[-1,-2],
  Tension[3,4],z4,Controls[z45,z54],z5,Cycle].

Its internal representation will be

CyclePath[
  Knot[Open[1],z0,Dir[20 Degree,1]],
  Knot[Open[1],z1,Open[-1]],
  Knot[Curl[2,1],z2,Curl[2,1]],
  Knot[Direction[ArcTan[-1,-2],1],z3,Direction[ArcTan[-1,-2],3]],
  Knot[Open[4.],z4,Control[z45]],
  Knot[Control[z54],z5,Open[1]]
]

TO DO: Combination of paths.
*)

AssignDefaults[path_] := path /.
  Dots -> Tension[-1] /. DDash -> Tension[100000] /. AtLeast[x_] -> -x /.
  Direction[x_,y_] -> Direction[ArcTan[x,y]] /.
 Direction[x_] :> N[Direction[x]] //.
  Path[a___,Dash,b___] -> Path[a,Curl[1],Curl[1],b] //.
  AddDefaultTension //. PropagateCurlOrDirection //. AddOpen //.
  Path[a___,Controls[x_,y_],b___] -> Path[a,Control[x],Control[y],b] /.
  Tension[x_] -> Tension[x,x] //. ContractTension /.
  BoundaryConditions //. FormKnots //. JoinEquals //
  (GoodPathQ[#];#)&  (* give error message if malformed *)

CurlOrDirectionQ[x_] := SameQ[Head[x],Curl] || SameQ[Head[x],Direction]

AddDefaultTension :=
 {Path[a___,x_?CurlOrDirectionQ,,b___] :> Path[a,x,Tension[1],b],
  Path[a___,,x_?CurlOrDirectionQ,b___] :> Path[a,Tension[1],x,b],
  Path[a___,x_List,y_List,b___] :> Path[a,x,Tension[1],y,b],
  Path[a___,x_List,Cycle] :> Path[a,x,Tension[1],Cycle]}

PropagateCurlOrDirection :=
 {Path[a___,y_?CurlOrDirectionQ,x_List,z_Tension,b___] :> Path[a,y,x,y,z,b],
  Path[a___,z_Tension,x_List,y_?CurlOrDirectionQ,b___] :> Path[a,z,y,x,y,b],
  Path[x_List,y_?CurlOrDirectionQ,a___,z_Tension,Cycle] :>
    Path[x,y,a,z,y,Cycle],
  Path[x_List,y_Tension,a___,z_?CurlOrDirectionQ,Cycle] :>
   Path[x,z,y,a,z,Cycle]}

AddOpen :=
 {Path[a___,x_List,y_Tension,b___] :> Path[a,x,Open[],y,b],
  Path[a___,y_Tension,x_List,b___] :> Path[a,y,Open[],x,b],
  Path[a___,y_Tension,Cycle] :> Path[a,y,Open[],Cycle]}

ContractTension :=
  Path[a___,x_[xx___],Tension[y_,z_],w_[ww___],b___] :>
    Path[a,x[xx,y],w[ww,z],b]

BoundaryConditions :=
 {Path[a___,z_,Cycle] :> CyclePath[z,a],
  Path[a___] :> OpenPath[EndPoint,a,EndPoint]}

FormKnots :=
 h_?PathQ[a___,x_,y_List,z_,b___] :> h[a,Knot[x,y,z],b]

JoinEquals =
  h_?PathQ[x___,Knot[l1_,p_,r1_?NotControlQ],Knot[l2_,p_,r2_],y___] :>
    h[x,Knot[l1 /. ToCurl,p,Control[p]],Knot[Control[p],p,r2 /. ToCurl],y]
    
ToCurl = Open[x___] :> Curl[1,x]

(* OUTPUT *)

(* The code to print a path is unreasonably complicated because it
attempts to imitate the way paths are input, or, optionally, the
way METAFONT inputs and ouptuts paths, so results can be compared more easily.

It works as follows: for each knot in the path we look left and right,
except that for the first knot we don't look left, and for the last we
don't look right.
*)

Format[OpenPath[x_,y___,z_]] := SequenceForm["Path[",
  SequenceForm @@ Apply[knot,Transpose[{{Null,x,y},{x,y,z},{y,z,Null}}],{1}],
  "]"]
MFForm[OpenPath[x_,y___,z_]] := ColumnForm[
  Apply[MFForm,Transpose[{{Null,x,y},{x,y,z},{y,z,Null}}],{1}]]

Format[CyclePath[x_,y__]] := SequenceForm["Path[",
  SequenceForm @@ Apply[knot,Transpose[{{Null,x,y},{x,y,x},{y,x,Cycle}}],{1}],
  "]"]
MFForm[CyclePath[x_,y__]] := ColumnForm[
  Apply[MFForm,Transpose[{{Null,x,y},{x,y,x},{y,x,Cycle}}],{1}]]

MFForm[{x_,y_}] := SequenceForm["(",x,",",y,")"]

Format[knot[Null,k2_,k3_]] := halfknot[k2,k3,1]
MFForm[Null,k2_,k3_] := ColumnForm[{" ",MFForm[halfknot[k2,k3,1]]}]

Format[knot[k1_,k2_,Null]] := SequenceForm[halfknot[k2,k1,3],k2[[2]]]
MFForm[k1_,k2_,Null] := 
  SequenceForm[MFForm[halfknot[k2,k1,3]],MFForm[k2[[2]]]]

Format[knot[k1_,k2_,Cycle]] := SequenceForm[halfknot[k2,k1,3],"Cycle"]
MFForm[k1_,k2_,Cycle] := SequenceForm[MFForm[halfknot[k2,k1,3]],"cycle"]

Format[knot[k1_,k2_,k3_]] := SequenceForm[halfknot[k2,k1,3],halfknot[k2,k3,1]]
MFForm[k1_,k2_,k3_] :=
  SequenceForm[MFForm[halfknot[k2,k1,3]],MFForm[halfknot[k2,k3,1]]]

(*
We're on knot k1, looking at k2.
When n=3 we're looking left, when n=1 we're looking right.
For n=1, for example, l1 is the left information of knot k1, and so on.
*)

MFForm[halfknot[k1_,k2_,n_]] := Block[{l1=k1[[n]],r1=k1[[-n]],l2=k2[[n]]},
  SequenceForm[
    If[n==1,MFForm[k1[[2]]]," .."],
    Switch[r1,
      Control[_],If[n==1,SequenceForm["..controls ",
        MFForm[First[r1]]," and ",MFForm[First[l2]]],""],
      Open[_],"",
      Curl[__],SequenceForm["{curl ",First[r1],"}"],
      Direction[__],
        SequenceForm["{",4096 Cos[First[r1]],",",4096 Sin[First[r1]],"}"]],
    If[n==1,MFForm[tension[gettension[r1],gettension[l2]]],""]
  ]
]

Format[halfknot[k1_,k2_,n_]] := Block[{l1=k1[[n]],r1=k1[[-n]],l2=k2[[n]]},
  SequenceForm[
    If[n==1,SequenceForm[k1[[2]],", "],""],
    Switch[r1,
      Control[_],
       If[n==1,SequenceForm["Controls[",First[r1],", ",First[l2],"], "],""],
      Open[_],"",
      Curl[__],SequenceForm[Drop[r1,-1],", "],
      Direction[__],SequenceForm[Drop[r1,-1],", "]],
    If[n==1,tension[gettension[r1],gettension[l2]],""]
  ]
]

Format[tension[t_,u_]] := 
  If[t==1 && u==1,"",
    SequenceForm["Tension[",tension[t],
      If[!SameQ[t,u],SequenceForm[", ",tension[u]],""],"], "]]
Format[tension[t_]] := If[t>0,t,SequenceForm["AtLeast[",-t,"]"]]

MFForm[tension[t_,u_]] := 
  If[t==1 && u==1,"",
    SequenceForm["..tension ",MFForm[tension[t]],
      If[!SameQ[t,u],SequenceForm[" and ",MFForm[tension[u]]],""]]]
MFForm[tension[t_]] := If[t>0,t,SequenceForm["atleast",-t]]

gettension[Open[t_]] := t
gettension[Direction[_,t_]] := t
gettension[Curl[_,t_]] := t
gettension[_] := 1
  
(* CHOOSING CONTROL POINTS *)
(*
We now get to Hobby's algorithm proper.
*)

AssignControls[x_Path] := AssignControls[AssignDefaults[x]]
AssignControls[path:(_?PathQ)[__]] := 
path //. {
  h_?PathQ[a___,x:Knot[_,_,_?NotControlQ]?((!OpenQ[#])&),y:(_?OpenQ...),
      z:Knot[_?NotControlQ,_,_]?((!OpenQ[#])&),b___] :>
    h @@ Join[{a},OpenFillIn[{x,y,z}],{b}],

  CyclePath[y:(_?OpenQ...),z:Knot[_?NotControlQ,_,_]?((!OpenQ[#])&),
      yy:(_?OpenQ...)] :>
    RotateRight[CyclePath @@ (OpenFillIn[{z,yy,y,z}]/.MergeEnds),Length[{y}]],

  CyclePath[y:(_?OpenQ...),z:Knot[_?NotControlQ,_,_]?((!OpenQ[#])&),a___,
      x:Knot[_,_,_?NotControlQ]?((!OpenQ[#])&),yy:(_?OpenQ...)] :>
    RotateRight[CyclePath @@ Join[{a},OpenFillIn[{x,yy,y,z}]],Length[{y}]+1],
[K
  CyclePath[y:(_?OpenQ..)] :>
    CyclePath @@ CyclicFillIn[{y}]
}

MergeEnds := {Knot[_,x_,a_],y___,Knot[b_,x_,_]} :> {Knot[b,x,a],y}

OpenFillIn[{a_,b___,c_}] := Block[
  {x={RemoveOpen[1,a],b,Reverse[RemoveOpen[-1,Reverse[c]]]},
   z,deltaz,psi,theta,phi,alpha,beta},
  z=Map[#[[2,1]]+I #[[2,2]]&,x];
  deltaz=N[Drop[RotateLeft[z]-z,-1]];
  psi=Append[Arg[deltaz/RotateRight[deltaz]],0]; (* ignore first *)
  alpha=(1/gettension[#[[3]]])& /@ Drop[x,-1];
  beta=(1/gettension[#[[1]]])& /@ Drop[x,1]; (* shifted *)
  theta=If[Length[z]==2 && Head[x[[1,-1]]]===Head[x[[-1,1]]]===Curl,
    {0,0},
    Tridiagonal @@ 
      MakeCoefficients[Abs[alpha],Abs[beta],deltaz,psi,x[[1,-1]],x[[-1,1]]]];
  phi= -psi-theta;
  ApplyHobby[z,Append[deltaz,0],
    Append[alpha,0],Prepend[beta,0],theta,phi,x[[1,1]],x[[-1,-1]]]
][K

MakeCoefficients[alpha_,beta_,deltaz_,psi_] :=
    {#1/(#2^2 #3),(3-#1)/(#2^2 #3),(3-#2)/(#1^2 #3),#2/(#1^2 #3)}& @@
      {alpha,RotateLeft[beta],Abs[deltaz]} //
    {Transpose[{RotateRight[#[[1]]],RotateRight[#[[2]]]+#[[3]],#[[4]]}],
      -psi RotateRight[#[[2]]] -RotateLeft[psi] #[[4]]}&

MakeCoefficients[alpha_,beta_,deltaz_,psi_,init_,final_] :=
    {#1/(#2^2 #3),(3-#1)/(#2^2 #3),(3-#2)/(#1^2 #3),#2/(#1^2 #3)}& @@
      {alpha,beta,Abs[deltaz]} //
    {Transpose[{RotateRight[#[[1]]],RotateRight[#[[2]]]+#[[3]],#[[4]]}],
      -Drop[psi,-1] RotateRight[#[[2]]] -Drop[psi,1] #[[4]]}& // Transpose //
    MapAt[FirstRow[init,alpha[[1]],beta[[1]],psi[[2]],deltaz[[1]]]&,#,{1}]& //
    Append[#,LastRow[final,alpha[[-1]],beta[[-1]],deltaz[[-1]]]]& //
    Transpose

FirstRow[Direction[d_,___],_,_,_,deltaz_] := {{0,1,0},Arg[E^(I d)/deltaz]}
FirstRow[Curl[gamma_,___],alpha_,beta_,psi_,_] :=
  Block[{chi=alpha^2gamma/beta^2},
    {{0,chi alpha+3-beta,(3-alpha)chi+beta},-((3-alpha)chi + beta)psi}]

LastRow[Direction[d_,___],_,_,deltaz_] := {{0,1,0},Arg[E^(I d)/deltaz]}
LastRow[Curl[gamma_,___],alpha_,beta_,_] := 
  Block[{chi=beta^2gamma/alpha^2},{{(3-beta)chi+alpha,chi beta+3-alpha,0},0}]

RemoveOpen[sign_,Knot[Control[z_],p_,Open[x___]]] :=
  Knot[Control[z],p,If[p==z,Curl[1,x],Direction[ArcTan@@(sign(p-z)),x]]]
RemoveOpen[_,x_] := x

CyclicFillIn[x_] := Block[
  {z,deltaz,psi,theta,phi,alpha,beta},
  z=Map[#[[2,1]]+I #[[2,2]]&,x];
  deltaz=N[RotateLeft[z]-z];
  psi=Arg[deltaz/RotateRight[deltaz]];
  alpha=(1/gettension[#[[3]]])& /@ x;
  beta=(1/gettension[#[[1]]])& /@ x;
  theta=CyclicTridiagonal @@
    MakeCoefficients[Abs[alpha],Abs[beta],deltaz,psi];
  phi= -psi-theta;
  ApplyHobby[z,deltaz,alpha,beta,theta,phi]
]

ApplyHobby[z_,deltaz_,alpha_,beta_,theta_,phi_] :=
 ({z-HobbyF[phi,RotateRight[theta],beta] E^(-I phi) RotateRight[deltaz],
      z,z+HobbyF[theta,RotateLeft[phi],alpha] E^(I theta) deltaz} ~
    Map[#2,#1,{2}]& ~ ({Re[#],Im[#]}&) //
    Transpose // Apply[Knot,#,{1}]&) /. 
    Knot[a_,b_,c_] -> Knot[Control[a],b,Control[c]]

ApplyHobby[z_,deltaz_,alpha_,beta_,theta_,phi_,init_,final_] :=
  ApplyHobby[z,deltaz,alpha,beta,theta,phi] //
    MapAt[init&,#,{{1,1}}]& // MapAt[final&,#,{{-1,-1}}]&

HobbyF[theta_,phi_,alpha_] := 
  Block[{f,st=Sin[theta],sp=Sin[phi],ct=Cos[theta],cp=Cos[phi],ssum},
    f=Min[4,Abs[alpha]*(2.+Sqrt[2.] (st-sp/16.)(sp-st/16.)(ct-cp))/
      (3.(1.+.5(Sqrt[5.]-1)ct+.5(3.-Sqrt[5.])cp))];
    If[alpha<0 && st sp>=0,
      Block[{ssum=Abs[st]cp+Abs[sp]ct},If[ssum>0,Min[f,Abs[sp]/ssum],f]],f]
  ]
SetAttributes[HobbyF,Listable]

(* DRAWING THE CURVE *)
(*
This is a slow implementation since we draw each Bezier patch at top
level.  If Mathematica had a bicubic spline primitive, or better
interface with PostScript, the drawing would be infinitely faster.
We could also change the call to ParametricPlot into something faster
(e.g., taking into account the resolution).
*)

Render[x_Path] := Render[AssignControls[x]]
Render[x:(_?PathQ)[__]] := Render[AssignControls[x],OK]
Render[OpenPath[x_,y___,z_],OK] := Patch /@ Transpose[{{x,y},{y,z}}]
Render[CyclePath[x_,y__],OK] := Patch /@ Transpose[{{x,y},{y,x}}]

Patch[{Knot[_,x_,Control[y_]],Knot[Control[z_],w_,_]}] := Patch[x,y,z,w]
Patch[a_,b_,c_,d_] := Block[{aa=a, bb=3(b-a), cc=3(c-2b+a), dd=d-3c+3b-a},
  ParametricPlot[aa+t(bb+t(cc+t dd)),{t,0,1},DisplayFunction->Identity] //
  First // First]

EndPackage[]

