CSE 477
Composite Bezier curves
G. Farin

bern[a_,b_,i_,t_]:=Binomial[n,i]*(b-t)^(n-i)*(t-a)^i/((b-a)^n);

comp[t_]:= If[t< b,bez[bbleft,a,b,n,t],bez[bbright,b,c,n,t]];

bez[pol_,u_,v_,n_,t_]:= Sum[pol[[i+1]]*bern[u,v,i,t],{i,0,n}];

first[t_]:= (comp[t+0.00001]-comp[t])/0.00001;
second[t_]:= (first[t+0.00001]-first[t])/0.00001;
 curvature[t_]:= Norm[Cross[first[t],second[t]]]/Norm[first[t]]^3; 

curvature[1.2]

n=3;
a=0; b=0.5; c=1;

(*G2 and C2: 
bbleft={{-3,2,2},{-2,1,0},{-1,0,0},{0,0,0}};
bbright={{0,0,0},{1.0,0,0},{2,1,0},{3,2,-1}};*)

(*G2 but not C2: *)
bbleft={{-3,2,1},{-2,1,0},{-1,0,0},{0,0,0}};
bbright={{0,0,0},{1.0,0,0},{3,1,0},{3,2,-1}};


curve=ParametricPlot3D[comp[t],{t,a,c},PlotStyle->{Thick},Boxed->False]; 
polyleft=Graphics3D[{Thick,Green,Line[bbleft]}];
polyright=Graphics3D[{Thick,Red,Line[bbright]}];
pointsleft=Graphics3D[{PointSize[Large],Green,Point[bbleft]}];
pointsright=Graphics3D[{PointSize[Large],Red,Point[bbright]}];

Show[{polyleft,polyright,curve,pointsleft,pointsright}]

Plot[curvature[t],{t,a,c},Filling->Axis,AxesOrigin->{0,0},PlotRange->All, PlotLabel->Style["curvature",FontSize->18]]
Plot[Norm[second[t]],{t,a,c},Filling->Axis,AxesOrigin->{0,0},PlotRange->All,PlotLabel->Style["second derivative magnitude",FontSize->18]]