Bezier Curves

Bezier Subdivision

1D Bezier Subdivision

subdivide1D[pts_]:=With[{n=Length[pts]-1},
With[{Q=Table[((((1/2)))^k)Sum[Binomial[k,j]pts[[j+1]],{j,0,k}],{k,0,n}],
R=Table[((((1/2)))^(n-k))Sum[Binomial[n-k,j]pts[[n-j+1]],{j,0,n-k}],{k,0,n}]},
{Q,R}
]
];
subdivide1D[pts_,n_]:=Module[{qs={},rs={},tmp1,tmp2},
If[n==0,{{pts},{pts}},
tmp1=subdivide1D[pts];
qs={tmp1[[1]]};
rs={tmp1[[2]]};
If[n>1,
Do[
tmp1=Table[subdivide1D[qs[[i]]],{i,Length[qs]}];
tmp2=Table[subdivide1D[rs[[i]]],{i,Length[rs]}];
qs=Join[tmp1[[Range[1,Length[tmp1]],1]],tmp1[[Range[1,Length[tmp1]],2]]];
rs=Join[tmp2[[Range[1,Length[tmp2]],1]],tmp2[[Range[1,Length[tmp2]],2]]];
,
{i,2,n}]
];
{qs,rs}
]
]

subdivide1D[{{0,0},{1,1},{2,0}}]

showCurve[p_,q_,r_]:=Graphics[{Black,Thickness[0.005],Line[p],Blue,Thickness[0.001],Line[q],Red,Line[r],PointSize[Large],Blue,Point[q],Red,Point[r]},ImageSize->600];

showCurves[p_,qs_,rs_]:=Graphics[{Black,Thickness[0.005],Line[p],Blue,Thickness[0.001],Table[Line[qs[[i]]],{i,Length[qs]}],Red,Table[Line[rs[[i]]],{i,Length[rs]}],PointSize[Large],Blue,Table[Point[qs[[i]]],{i,Length[qs]}],Red,Table[Point[rs[[i]]],{i,Length[rs]}]},ImageSize->600];

With[{p={{0,0},{1,1},{2,0}}},
With[{qr=subdivide1D[p]},
showCurve[p,qr[[1]],qr[[2]]]
]
]

With[{p={{0,0},{1,1},{2,1},{3,0}}},
With[{qr=subdivide1D[p,0]},
showCurves[p,qr[[1]],qr[[2]]]
]
]

Tensor Product Bezier Subdivision

subdivide2D[pts_]:=Module[{ns=Length[pts]-1,nt=Length[pts[[1]]]-1,ct,tmp,cs},
ct=Table[
With[{bar=subdivide1D[pts[[All,i]]]},
Join[bar[[1]],Rest[bar[[2]]]]]
,{i,nt+1}];
cs=Table[
With[{bar=subdivide1D[ct[[All,i]]]},
Join[bar[[1]],Rest[bar[[2]]]]
]
,{i,Length[ct[[1]]]}];
ns++;nt++;
{{cs[[Range[1,ns],Range[1,nt]]],cs[[Range[1,ns],Range[nt,2nt-1]]]},
{cs[[Range[ns,2ns-1],Range[1,nt]]],cs[[Range[ns,2ns-1],Range[nt,2nt-1]]]}}
];

subdivide2D[pts_,n_]:=Module[{s00s={},s01s={},s10s={},s11s={},tmp1,tmp2,tmp3,tmp4},
If[n==0,
{{{pts},{pts}},{{pts},{pts}}}
,
tmp1=subdivide2D[pts];
s00s={tmp1[[1,1]]};s01s={tmp1[[1,2]]};
s10s={tmp1[[2,1]]};s11s={tmp1[[2,2]]};

If[n>1,
Do[
tmp1=Table[subdivide2D[s00s[[i]]],{i,Length[s00s]}];
tmp2=Table[subdivide2D[s01s[[i]]],{i,Length[s01s]}];
tmp3=Table[subdivide2D[s10s[[i]]],{i,Length[s10s]}];
tmp4=Table[subdivide2D[s11s[[i]]],{i,Length[s11s]}];

s00s=Join[tmp1[[All,1,1]],tmp1[[All,1,2]],tmp1[[All,2,1]],tmp1[[All,2,2]]];
s01s=Join[tmp2[[All,1,1]],tmp2[[All,1,2]],tmp2[[All,2,1]],tmp2[[All,2,2]]];
s10s=Join[tmp3[[All,1,1]],tmp3[[All,1,2]],tmp3[[All,2,1]],tmp3[[All,2,2]]];
s11s=Join[tmp4[[All,1,1]],tmp4[[All,1,2]],tmp4[[All,2,1]],tmp4[[All,2,2]]];

,{i,2,n}]
];
{{s00s,s01s},{s10s,s11s}}
]
];

showSurface[p_,{{s00_,s01_},{s10_,s11_}}]:=With[{s=Length[p]-1,t=Length[p[[1]]]-1},Graphics3D[{
Opacity[1.0],Magenta,EdgeForm[],
Table[Polygon[{s00[[i,j]],s00[[i+1,j]],s00[[i+1,j+1]],s00[[i,j+1]]}],{i,s},{j,t}],
Table[Polygon[{s01[[i,j]],s01[[i+1,j]],s01[[i+1,j+1]],s01[[i,j+1]]}],{i,s},{j,t}],
Table[Polygon[{s10[[i,j]],s10[[i+1,j]],s10[[i+1,j+1]],s10[[i,j+1]]}],{i,s},{j,t}],
Table[Polygon[{s11[[i,j]],s11[[i+1,j]],s11[[i+1,j+1]],s11[[i,j+1]]}],{i,s},{j,t}],Cyan,Opacity[0.25],EdgeForm[Thick],Table[Polygon[{p[[i,j]],p[[i+1,j]],p[[i+1,j+1]],p[[i,j+1]]}],{i,s},{j,t}]},ImageSize->300]
];
showSurfaces[p_,{{s00s_,s01s_},{s10s_,s11s_}}]:=With[{s=Length[p]-1,t=Length[p[[1]]]-1},Graphics3D[{
Opacity[1.0],Magenta,EdgeForm[],
Table[Table[Polygon[{s00s[[k,i,j]],s00s[[k,i+1,j]],s00s[[k,i+1,j+1]],s00s[[k,i,j+1]]}],{i,s},{j,t}],{k,Length[s00s]}],
Table[Table[Polygon[{s01s[[k,i,j]],s01s[[k,i+1,j]],s01s[[k,i+1,j+1]],s01s[[k,i,j+1]]}],{i,s},{j,t}],{k,Length[s01s]}],
Table[Table[Polygon[{s10s[[k,i,j]],s10s[[k,i+1,j]],s10s[[k,i+1,j+1]],s10s[[k,i,j+1]]}],{i,s},{j,t}],{k,Length[s10s]}],
Table[Table[Polygon[{s11s[[k,i,j]],s11s[[k,i+1,j]],s11s[[k,i+1,j+1]],s11s[[k,i,j+1]]}],{i,s},{j,t}],{k,Length[s11s]}],EdgeForm[Thick],Cyan,Opacity[0.25],Table[Polygon[{p[[i,j]],p[[i+1,j]],p[[i+1,j+1]],p[[i,j+1]]}],{i,s},{j,t}]},ImageSize->300]
]

subdivisionApp[pts_]:=DynamicModule[{showHull=True,cpts=pts,s=Length[pts]-1,t=Length[pts[[1]]]-1,n=0,s00s,s01s,s10s,s11s,reset,mins,maxs},
mins={Min[pts[[All,All,1]]],Min[pts[[All,All,2]]],Min[pts[[All,All,3]]]};
maxs={Max[pts[[All,All,1]]],Max[pts[[All,All,2]]],Max[pts[[All,All,3]]]};

reset=Function[{tn},{{s00s,s01s},{s10s,s11s}}=subdivide2D[cpts,tn];n];
{{s00s,s01s},{s10s,s11s}}=subdivide2D[cpts,n];
Panel[Column[{Row[{EventHandler[Dynamic@Graphics3D[{
Opacity[1.0],Magenta,EdgeForm[Thin],Specularity[White,20],
Table[Table[Polygon[{s00s[[k,i,j]],s00s[[k,i+1,j]],s00s[[k,i+1,j+1]],s00s[[k,i,j+1]]}],{i,s},{j,t}],{k,Length[s00s]}],
Table[Table[Polygon[{s01s[[k,i,j]],s01s[[k,i+1,j]],s01s[[k,i+1,j+1]],s01s[[k,i,j+1]]}],{i,s},{j,t}],{k,Length[s01s]}],
Table[Table[Polygon[{s10s[[k,i,j]],s10s[[k,i+1,j]],s10s[[k,i+1,j+1]],s10s[[k,i,j+1]]}],{i,s},{j,t}],{k,Length[s10s]}],
Table[Table[Polygon[{s11s[[k,i,j]],s11s[[k,i+1,j]],s11s[[k,i+1,j+1]],s11s[[k,i,j+1]]}],{i,s},{j,t}],{k,Length[s11s]}],If[showHull,{Specularity[0,1],EdgeForm[Thick],Cyan,Opacity[0.0],Table[Polygon[{cpts[[i,j]],cpts[[i+1,j]],cpts[[i+1,j+1]],cpts[[i,j+1]]}],{i,s},{j,t}]},{}]},PlotRange->{{Min[mins],Max[maxs]},{Min[mins],Max[maxs]},{Min[mins],Max[maxs]}},Boxed->False,(*RotationAction->"Clip",*)ImageSize->300],
{(*"MouseClicked":>(col=Red)*)}
]}],
Row[Table[With[{i=i},Button["Level "<>ToString[i],n=i;reset[n]]],{i,0,4}]]
(*{Slider[Dynamic[n],{1,4,1}],Dynamic[reset[n]]}*),
Row[{Labeled[RadioButton[Dynamic[showHull],False],"Hide Control Polygon"],Labeled[RadioButton[Dynamic[showHull],True],"Show Control Polygon"]}],
}],ImageSize->300]]

Lagrange Interpolation and Neville's Algorithm

neville[cpts_,nodes_,level_,t_]:=If[Length[cpts]==1,cpts[[1]],neville[Table[((nodes[[i+level]]-t)cpts[[i]]+(t-nodes[[i]]) cpts[[i+1]])/(nodes[[i+level]]-nodes[[i]]),{i,(Length[cpts]-1)}],nodes,(level+1),t]]

Linear

Manipulate[Module[{nodes={0,1}},Show[{Graphics[{{Thickness[.004],Line[cpts]},{Blue,Thickness[.004],Line[cpts]}}],ParametricPlot[neville[cpts,nodes,1,t],{t,0,1},PlotStyle->Thickness[.009],DisplayFunction->Identity]},Axes->False,ImageSize->250,AspectRatio->Automatic,PlotRange->{{-2,4},{-2,4}}]],{{cpts,{{-1,0},{3,0}}},Locator}]

Quadratic

Manipulate[Module[{nodes={0,1,2}},Show[{Graphics[{{Thickness[.004],Line[cpts]},{Blue,Thickness[.004],Line[cpts]}}],ParametricPlot[neville[cpts,nodes,1,t],{t,0,2},PlotStyle->Thickness[.009],DisplayFunction->Identity]},Axes->False,ImageSize->250,AspectRatio->Automatic,PlotRange->{{-2,4},{-2,4}}]],{{cpts,{{-1,0},{1,2},{3,0}}},Locator}]

Cubic

Manipulate[Module[{nodes={0,1,2,3}},Show[{Graphics[{{Thickness[.004],Line[cpts]},{Blue,Thickness[.004],Line[cpts]}}],ParametricPlot[neville[cpts,nodes,1,t],{t,0,3},PlotStyle->Thickness[.009],DisplayFunction->Identity]},Axes->False,ImageSize->250,AspectRatio->Automatic,PlotRange->{{-2,6},{-2,4}}]],{{cpts,{{-1,0},{1,2},{3,2},{5,0}}},Locator}]

Quartic

Manipulate[Module[{nodes={0,1,2,3,4}},Show[{Graphics[{{Thickness[.004],Line[cpts]},{Blue,Thickness[.004],Line[cpts]}}],ParametricPlot[neville[cpts,nodes,1,t],{t,0,4},PlotStyle->If[display,Thickness[.008],PlotStyle->Thickness[.000001]],DisplayFunction->Identity]},Axes->False,ImageSize->250,AspectRatio->Automatic,PlotRange->{{-1,5},{-2,4}}]],{{cpts,{{0,0},{1,0},{2,2},{3,0},{4,0}}},Locator},{display,{True,False}},TrackedSymbols->{display,cpts}]

Sextic

Manipulate[Module[{nodes={0,1,2,3,4,5,6}},Show[{Graphics[{{Thickness[.004],Line[cpts]},{Blue,Thickness[.004],Line[cpts]}}],ParametricPlot[neville[cpts,nodes,1,t],{t,0,6},PlotStyle->If[display,Thickness[.008],PlotStyle->Thickness[.0000001]],DisplayFunction->Identity]},Axes->False,ImageSize->250,AspectRatio->Automatic,PlotRange->{{-2,12},{-4,8}}]],{{cpts,{{-1,0},{1,0},{3,0},{5,5},{7,0},{9,0},{11,0}}},Locator},{display,{True,False}},TrackedSymbols->{display,cpts}]

Bezier Approximation and the De Casteljau Algorithm

Setup

Clear[average,deCast,leftsubd,rightsubd,subd];

average[pt_,qt_]:=(pt+qt)/2

deCast[cpts_]:=Module[{dpts={cpts}},If[Length[cpts]==1,dpts,Do[AppendTo[dpts,Table[average[dpts[[k,i]], dpts[[k,i+1]]],{i,Length[dpts[[k]]]-1}]],{k,1,Length[cpts]-1}]];dpts];

leftsubd[cpts_]:=Table[First[deCast[cpts][[k]]],{k,1,Length[cpts]}];

rightsubd[cpts_]:=Table[Last[deCast[cpts][[k]]],{k,Length[cpts],1,-1}]

subd[cpts_]:=Flatten[Table[{leftsubd[cpts[[k]]],rightsubd[cpts[[k]]]},{k,1,Length[cpts]}],1]

Quadratic

Manipulate[Module[{dpts={cpts}},dpts=Nest[subd,dpts,level];epts=Nest[subd,dpts,4];Graphics[{{Thickness[.004],Line[cpts]},{Darker[Red],Thickness[.008],Line[epts]},{Blue,Thickness[.004],Line[dpts]}},Axes->False,ImageSize->250,AspectRatio->Automatic,PlotRange->{{-2,4},{-2,4}}]],{{cpts,{{-1,0},{1,2},{3,0}}},Locator},{{level,0},{0,1,2,3,4}}]

Cubic

Manipulate[Module[{dpts={cpts}},dpts=Nest[subd,dpts,level];epts=Nest[subd,dpts,4];Graphics[{{Thickness[.004],Line[cpts]},{Darker[Red],Thickness[.008],Line[epts]},{Blue,Thickness[.004],Line[dpts]}},Axes->False,ImageSize->250,AspectRatio->Automatic,PlotRange->{{-1,4},{-1,3}}]],{{cpts,{{0,0},{.5,1.5},{2.5,1.5},{3,0}}},Locator},{{level,0},{0,1,2,3,4}}]

Sextic

Manipulate[Module[{dpts={cpts}},dpts=Nest[subd,dpts,level];epts=Nest[subd,dpts,4];Graphics[{{Thickness[.004],Line[cpts]},{Darker[Red],Thickness[.008],Line[epts]},{Blue,Thickness[.004],Line[dpts]}},Axes->False,ImageSize->250,AspectRatio->Automatic,PlotRange->{{-2,12},{-2,6}}]],{{cpts,{{-1,0},{1,0},{3,0},{5,5},{7,0},{9,0},{11,0}}},Locator},{{level,0},{0,1,2,3,4}}]

Tensor Product Bezier Patches

Bernstein Basis Functions

bernsteinbasis[k_,n_,t_]:=Binomial[n,k](t^k)(((1-t))^(n-k))


bernsteinbasissurf[j_,k_,m_,n_,s_,t_]:=bernsteinbasis[j,m,s]*bernsteinbasis[k,n,t]

Bezier Curves and Surfaces

bezier[cpts_,t_]:=(∑_(k=0)^(Length[cpts]-1))bernsteinbasis[k,Length[cpts]-1,t]*cpts[[k+1]];

beziersurf[cpts_,s_,t_]:=(∑_(j=0)^(Length[cpts[[1]]]-1))(∑_(k=0)^(Length[cpts]-1))bernsteinbasissurf[j,k,Length[cpts[[1]]]-1,Length[cpts]-1,s,t]*cpts[[j+1,k+1]]

Example

Clear[cpts,cpoly,bsurf];

cpts={{{0,0,0},{1,0,1},{2,0,1},{3,0,0}},{{0,1,2},{1,1,3},{2,1,3},{3,1,2}},{{0,2,2},{1,2,3},{2,2,3},{3,2,2}},{{0,3,0},{1,3,1},{2,3,1},{3,3,0}}};

cpoly[cpts_]:=Graphics3D[{Thickness[.006],Table[{Line[cpts[[i]]],Line[Transpose[cpts][[i]]]},{i,1,4}]},PlotRange->{{-1,4},{-1,4},{-1,4}},Boxed->False,DisplayFunction->Identity]

bsurf[cpts_]:=ParametricPlot3D[Evaluate[beziersurf[cpts,s,t]],{s,1.`*^-8,0.99999999`},{t,1.`*^-8,0.99999999`},PlotRange->{{-1,4},{-1,4},{-1,4}},Boxed->False,Axes->False]

Module[{z=3},Show[{bsurf[cpts],cpoly[cpts]},AspectRatio->Automatic,ImageSize->300]]

bsurf[cpts]

Clear[scurve,tcurve,bcurv];

tcurve=Table[ParametricPlot3D[Evaluate[bezier[cpts[[i]],t]],{t,.00000001,.99999999},DisplayFunction->Identity],{i,1,4}];

scurve=Table[ParametricPlot3D[Evaluate[bezier[Transpose[cpts][[i]],s]],{s,.00000001,.99999999},DisplayFunction->Identity],{i,1,4}];

bcurv=Show[Graphics3D[Thickness[.01]],scurve,DisplayFunction->$DisplayFunction];

Show[Graphics3D[Thickness[.01]],scurve,DisplayFunction->$DisplayFunction]

Show[bsurf[cpts],{Graphics3D[Thickness[.01]],scurve},DisplayFunction->$DisplayFunction]

Show[bsurf[cpts],{Graphics3D[Thickness[.01]],scurve,tcurve},DisplayFunction->$DisplayFunction]

Module[{z=3},Show[{bsurf[cpts],cpoly[cpts],Graphics3D[Thickness[.01]],scurve,tcurve},AspectRatio->Automatic,ImageSize->300]]

Application

quadraticPatch={{{0,0,0},{1,1,0},{2,0,0}},{{0,0,1},{1,1,1},{2,0,1}},{{0,0,2},{1,1,2},{2,0,2}}};
cubicPatch={{{0,0,0},{1,0,0},{2,0,0},{3,0,0}},{{0,0,1},{1,1,1},{2,1,1},{3,0,1}},{{0,0,2},{1,1,2},{2,1,2},{3,0,2}},{{0,0,3},{1,0,3},{2,0,3},{3,0,3}}};

subdivisionApp[quadraticPatch]

subdivisionApp[cubicPatch]