summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/doc/metapost/featpost/example/splineperspective.mp
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/doc/metapost/featpost/example/splineperspective.mp')
-rw-r--r--Master/texmf-dist/doc/metapost/featpost/example/splineperspective.mp550
1 files changed, 550 insertions, 0 deletions
diff --git a/Master/texmf-dist/doc/metapost/featpost/example/splineperspective.mp b/Master/texmf-dist/doc/metapost/featpost/example/splineperspective.mp
new file mode 100644
index 00000000000..c65cd5b2ce4
--- /dev/null
+++ b/Master/texmf-dist/doc/metapost/featpost/example/splineperspective.mp
@@ -0,0 +1,550 @@
+% splineperspective.mp
+% L. Nobre G. and Troy Henderson
+% 2007
+
+input featpost3Dplus2D;
+
+prologues := 1;
+
+% Evaluate a cubic polynomial of the "standard" Bezier form at t
+vardef evalbezier(expr p,t) =
+ save _a,_b,_c,_d;
+ numeric _a,_b,_c,_d;
+ _a:=(1-t)**3;
+ _b:=3*((1-t)**2)*t;
+ _c:=3*(1-t)*(t**2);
+ _d:=t**3;
+ (point 0 of p)*_a + (postcontrol 0 of p)*_b + (precontrol 1 of p)*_c +
+ (point 1 of p)*_d
+enddef;
+
+% Evaluate the derivative of a cubic polynomial of the "standard"
+% Bezier form at t
+vardef evalbezierderivative(expr p,t) =
+ save _a,_b,_c;
+ pair _a,_b,_c;
+ _a:=3*((point 1 of p) - 3*(precontrol 1 of p) + 3*(postcontrol 0 of p)
+ -(point 0 of p));
+ _b:=6*((precontrol 1 of p) - 2*(postcontrol 0 of p) + (point 0 of p));
+ _c:=3*((postcontrol 0 of p) - (point 0 of p));
+ _a*(t**2) + _b*t + _c
+enddef;
+
+% Evaluate a rational function of the "standard" cubic NURBS form at t
+vardef evalnurbs(expr p,w,t) =
+ save _q,_r;
+ path _q,_r;
+ _q:=((cyanpart w)*(point 0 of p))..
+ controls ((magentapart w)*(postcontrol 0 of p))
+ and ((yellowpart w)*(precontrol 1 of p)) .. ((blackpart w)*(point 1 of p));
+ _r:=(cyanpart w,0) ..
+ controls (magentapart w,0) and (yellowpart w,0) .. (blackpart w,0);
+ evalbezier(_q,t)/(xpart evalbezier(_r,t))
+enddef;
+
+% Evaluate the derivative of a rational function of the "standard"
+% cubic NURBS form at t
+vardef evalnurbsderivative(expr p,w,t) =
+ save _a,_b,_c,_d,_q,_r;
+ pair _a,_b;
+ numeric _c,_d;
+ path _q,_r;
+ _q:=((cyanpart w)*(point 0 of p)) ..
+ controls ((magentapart w)*(postcontrol 0 of p))
+ and ((yellowpart w)*(precontrol 1 of p)) .. ((blackpart w)*(point 1 of p));
+ _r:=(cyanpart w,0) ..
+ controls (magentapart w,0) and (yellowpart w,0) .. (blackpart w,0);
+ _a:=evalbezier(_q,t);
+ _b:=evalbezierderivative(_q,t);
+ _c:=xpart evalbezier(_r,t);
+ _d:=xpart evalbezierderivative(_r,t);
+ (_b*_c-_a*_d)/(_c**2)
+enddef;
+
+% Fit a cubic polynomial of the "standard" Bezier form to a
+% rational function of the
+% "standard" cubic NURBS form with function and derivative agreement
+% at tmin and tmax
+vardef nurbstobezier(expr p,w,tmin,tmax) =
+ save _a,_b,_c,_d,_e;
+ pair _a,_b,_c,_d;
+ numeric _e;
+ _e:=(tmax-tmin)/3;
+ _a:=evalnurbs(p,w,tmin);
+ _b:=_a + _e*evalnurbsderivative(p,w,tmin);
+ _d:=evalnurbs(p,w,tmax);
+ _c:=_d - _e*evalnurbsderivative(p,w,tmax);
+ _a .. controls _b and _c .. _d
+enddef;
+
+% Reparameterize a cubic polynomial of the "standard" Bezier form by mapping
+% the interval [tmin,tmax] to [0,1]
+vardef beziertobezier(expr p,tmin,tmax) =
+ nurbstobezier(p,(1,1,1,1),tmin,tmax)
+enddef;
+
+% Evalute the L^2[0,1] norm of a cubic polynomial of the "standard"
+% Bezier form
+vardef beziernorm(expr p) =
+ save _a,_b,_c,_d,_i,_xabs,_yabs,_A,_B,_C,_D,_I;
+ numeric _a,_b,_c,_d,_i,_xabs,_yabs,_A,_B,_C,_D,_I;
+ _xabs:=max(
+ abs(xpart point 0 of p),
+ abs(xpart postcontrol 0 of p),
+ abs(xpart precontrol 1 of p),
+ abs(xpart point 1 of p));
+ _yabs:=max(
+ abs(ypart point 0 of p),
+ abs(ypart postcontrol 0 of p),
+ abs(ypart precontrol 1 of p),
+ abs(ypart point 1 of p));
+ if (_xabs > 0):
+ _a:=xpart((point 1 of p) - 3*(precontrol 1 of p)
+ + 3*(postcontrol 0 of p) - (point 0 of p))/_xabs;
+ _b:=3*xpart((precontrol 1 of p) - 2*(postcontrol 0 of p)
+ + (point 0 of p))/_xabs;
+ _c:=3*xpart((postcontrol 0 of p) - (point 0 of p))/_xabs;
+ _d:=xpart(point 0 of p)/_xabs;
+ _i:=(_a**2)/7 + ((_b)**2 + 2*_a*_c)/5 + (_a*_b + 2*_b*_d + (_c**2))/3 + (_a*_d + _b*_c)/2 + (_c*_d + (_d**2));
+ else:
+ _i:=0;
+ fi;
+ if (_yabs > 0):
+ _A:=ypart((point 1 of p) - 3*(precontrol 1 of p)
+ + 3*(postcontrol 0 of p) - (point 0 of p))/_yabs;
+ _B:=3*ypart((precontrol 1 of p) - 2*(postcontrol 0 of p)
+ + (point 0 of p))/_yabs;
+ _C:=3*ypart((postcontrol 0 of p) - (point 0 of p))/_yabs;
+ _D:=ypart(point 0 of p)/_yabs;
+ _I:=(_A**2)/7 + ((_B)**2 + 2*_A*_C)/5
+ + (_A*_B + 2*_B*_D + (_C**2))/3 + (_A*_D + _B*_C)/2 + (_C*_D + (_D**2));
+ else:
+ _I:=0;
+ fi;
+ (_xabs*sqrt(_i)) ++ (_yabs*sqrt(_I))
+enddef;
+
+% Fit a cubic Bezier spline to a rational function of the "standard"
+% cubic NURBS form by iteratively refining the Bezier curve.
+% p is a 4 point path containing the 4 cubic NURBS (2D) control points
+% w is a cmykcolor containing the 4 cubic NURBS weights
+% EPS is the tolerance to stop refining each branch of the Bezier spline
+vardef fitnurbswithbezier(expr p,w,EPS) =
+ save _a,_b,_c,_e,_error,_k,_q;
+ numeric _a,_b,_c,_error,_k;
+ path _q,_q[],_e;
+ _a:=0;
+ _b:=1;
+ _k:=1/sqrt(2);
+ _q:=(point 0 of p);
+ _q[4]:=nurbstobezier(p,w,_a,_b);
+ forever:
+ exitunless(_a<1);
+ _q[1]:=_q[4];
+ _c:=_b-_k*((_b-_a)**2);
+ _q[2]:=beziertobezier(_q[1],_a,_c);
+ _q[3]:=nurbstobezier(p,w,_a,_c);
+ _q[4]:=_q[3];
+ _e:=((point 0 of _q[2])-(point 0 of _q[3])) ..
+ controls ((postcontrol 0 of _q[2])-(postcontrol 0 of _q[3]))
+ and ((precontrol 1 of _q[2])-(precontrol 1 of _q[3])) ..
+ ((point 1 of _q[2])-(point 1 of _q[3]));
+ _error:=beziernorm(_e)/beziernorm(_q[3]);
+% show _error;
+ if (_error > EPS):
+ _b:=_c;
+ else:
+ _q[2]:=beziertobezier(_q[1],_c,_b);
+ _q[3]:=nurbstobezier(p,w,_c,_b);
+ _e:=((point 0 of _q[2])-(point 0 of _q[3])) ..
+ controls ((postcontrol 0 of _q[2])-(postcontrol 0 of _q[3]))
+ and ((precontrol 1 of _q[2])-(precontrol 1 of _q[3])) ..
+ ((point 1 of _q[2])-(point 1 of _q[3]));
+ _error:=beziernorm(_e)/beziernorm(_q[3]);
+ if (_error > EPS):
+ _q:=_q .. controls (postcontrol 0 of _q[4])
+ and (precontrol 1 of _q[4]) .. (point 1 of _q[4]);
+ _a:=_c;
+ _q[4]:=_q[3];
+ else:
+ _q:=_q .. controls (postcontrol 0 of _q[1])
+ and (precontrol 1 of _q[1]) .. (point 1 of _q[1]);
+ _a:=_b;
+ _q[4]:=nurbstobezier(p,w,_a,1);
+ fi;
+ _b:=1;
+ fi;
+ endfor;
+ _q
+enddef;
+
+% This macro is used to provide a path to draw the NURBS
+% It returns a path of length N passing through N+1 equally spaced
+% (in time) points along the NURBS connected by line segments
+vardef samplednurbs(expr p,w,N) =
+ save _a,_b,_c,_d,_n,_t,_q;
+ numeric _a,_b,_c,_d,_n,_t;
+ path _q;
+ _q:=(point 0 of p);
+ for _n=1 upto N:
+ _t:=_n/N;
+ _a:=(cyanpart w)*((1-_t)**3);
+ _b:=3*(magentapart w)*((1-_t)**2)*_t;
+ _c:=3*(yellowpart w)*(1-_t)*(_t**2);
+ _d:=(blackpart w)*(_t**3);
+ _q:=_q .. ((_a*(point 0 of p)+_b*(postcontrol 0 of p)
+ +_c*(precontrol 1 of p)+_d*(point 1 of p))/(_a+_b+_c+_d));
+ endfor;
+ ( _q )
+enddef;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Here's where the fun begins %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+beginfig(4);
+% p contains the 4 control points of the rational function of the
+% "standard" cubic NURBS form
+ path p;
+ p:=(297.63725,297.63725) .. controls (132.98871,286.67885) and (180.62535,152.16249) .. (429.54399,226.31157);
+
+% w contains the 4 weights for the rational function of the
+% "standard" cubic NURBS form
+ cmykcolor w;
+ w:=(2.15756,1.6709,0.8598,1.34647);
+
+% EPS represents the minimum "acceptable error" to stop refining any
+% given branch of the Bezier
+ Err:=0.040;
+
+% q represents the Bezier spline fit to the rational function of the
+% "standard" cubic NURBS form
+ path q;
+ q:=fitnurbswithbezier(p,w,Err);
+% q:=fitnurbswithbezier(reverse p,(blackpart w,yellowpart w,magentapart w,cyanpart w),Err);
+
+% draw the NURBS by sampling it at many points and connecting the
+% samples via line segments
+ draw samplednurbs(p,w,20) withcolor red withpen pencircle scaled 2bp;
+
+% draw the Bezier spline and its knots
+ draw q;
+ for n=0 upto length q:
+ draw fullcircle scaled 2 shifted point n of q withcolor blue;
+ endfor;
+endfig;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+def casteljau( expr Za, Zb, Zc, Zd, Pt ) = %%%%%%%%%%%%%%%%%%% 2D or 3D
+ begingroup
+ save A, B, C, D;
+ numeric A, B, C, D;
+ A = (1-Pt)**3;
+ B = 3*((1-Pt)**2)*Pt;
+ C = 3*(1-Pt)*(Pt**2);
+ D = Pt**3;
+ ( (A*Za+B*Zb+C*Zc+D*Zd) )
+ endgroup
+enddef;
+
+def twothr( expr Z ) = ( xpart Z, ypart Z, 0 ) enddef;
+
+def twotwo( expr Z ) = rp( twothr( Z ) ) enddef;
+
+def xoy( expr Z ) = rp( ( X(Z), Y(Z), 0 ) ) enddef;
+
+def yoz( expr W ) = rp( ( 0, Y(W), Z(W) ) ) enddef;
+
+def xoz( expr W ) = rp( ( X(W), 0, Z(W) ) ) enddef;
+
+f := 0.35*(3,5,2);
+Spread := 160;
+
+def nextthirty( expr Za, Zb, Zc, Zd, Pt ) = %%% input 3D and return 2D
+ begingroup
+ save A, B, C, D, Tot, P;
+ numeric A, B, C, D, Tot;
+ color P;
+ P = N( f - viewcentr );
+ A = ((1-Pt)**3)*cdotprod( P, f-Za );
+ B = 3*((1-Pt)**2)*Pt*cdotprod( P, f-Zb );
+ C = 3*(1-Pt)*(Pt**2)*cdotprod( P, f-Zc );
+ D = (Pt**3)*cdotprod( P, f-Zd );
+ Tot = A+B+C+D;
+ ( (A*rp(Za)+B*rp(Zb)+C*rp(Zc)+D*rp(Zd))/Tot )
+ endgroup
+enddef;
+
+vardef nurbstobezierold (expr p,w) =
+ save _a,_b,_c,_d,_j,_n,_r,_s,_t,_A,_B,_Aold,_Bold,_C,_D,_EPS,_J,_N;
+ _EPS:=0.00001;
+ _J:=10;
+
+ _Aold:=0;
+ _Bold:=0;
+ _A:=1;
+ _B:=1;
+ _s:=((_A-_Aold)++(_B-_Bold))/(_A++_B);
+ _j:=1; _r:=0;
+ forever:
+ exitunless((_s>_EPS) and (_j<_J));
+ _j:=_j+1;
+ _N:=2**_j;
+ _Aold:=_A;
+ _Bold:=_B;
+ _D:=_N+1/_N-21/_N/_N/_N-1/_N/_N/_N/_N/_N+20/_N/_N/_N/_N/_N/_N/_N;
+ _C:=120*(2+2/_N/_N-5/_N/_N/_N/_N)/_D;
+ _D:=60*(3+3/_N/_N+10/_N/_N/_N/_N)/_D;
+ _c:=5/_N/_N/_N/_N;
+ _a:=2+2/_N/_N-_c;
+ _b:=2-3/_N/_N+_c;
+ _c:=1+6/_N/_N+_c;
+ _A:=(-2*(cyanpart p)*_a+(blackpart p)*_b)/_c;
+ _B:=((cyanpart p)*_b-2*(blackpart p)*_a)/_c;
+ for _n=0 upto _N:
+ _t:=_n/_N;
+ _a:=(1-_t)**3;
+ _b:=((1-_t)**2)*_t;
+ _c:=(1-_t)*(_t**2);
+ _d:=_t**3;
+ _r:=((cyanpart w)*(cyanpart p)*_a + 3*(magentapart
+ w)*(magentapart p)*_b + 3*(yellowpart w)*(yellowpart p)*_c +
+ (blackpart w)*(blackpart p)*_d)/((cyanpart w)*_a + 3*(magentapart
+ w)*_b + 3*(yellowpart w)*_c + (blackpart w)*_d);
+ _A:=_A+(_C*_b-_D*_c)*_r;
+ _B:=_B+(_C*_c-_D*_b)*_r;
+ endfor;
+ _s:=((_A-_Aold)++(_B-_Bold))/(_A++_B);
+ endfor;
+ (_A,_B)/3
+enddef;
+
+def nurbsapprox( expr Pa, Pb, Pc, Pd ) =
+ begingroup
+ color Pn;
+ numeric wa, wb, wc, wd;
+ path returnpath;
+ pair xpair, ypair, ba, bb, bc, bd;
+ cmykcolor xcontrols, ycontrols;
+ Pn = N( f - viewcentr );
+ wa = cdotprod( Pn, f-Pa );
+ wb = cdotprod( Pn, f-Pb );
+ wc = cdotprod( Pn, f-Pc );
+ wd = cdotprod( Pn, f-Pd );
+ xcontrols = (xpart rp(Pa),xpart rp(Pb),xpart rp(Pc),xpart rp(Pd));
+ ycontrols = (ypart rp(Pa),ypart rp(Pb),ypart rp(Pc),ypart rp(Pd));
+ xpair = nurbstobezierold( xcontrols, (wa,wb,wc,wd) );
+ ypair = nurbstobezierold( ycontrols, (wa,wb,wc,wd) );
+ ba = rp( Pa );
+ bb = (xpart xpair, xpart ypair);
+ bc = (ypart xpair, ypart ypair);
+ bd = rp( Pd );
+ %show ba;
+ %show bb;
+ %show bc;
+ %show bd;
+ %show wa;
+ %show wb;
+ %show wc;
+ %show wd;
+ returnpath = ba..controls bb and bc..bd;
+ (returnpath)
+ endgroup
+enddef;
+
+def fitthreednurbswithtwodbezier( expr pa, pb, pc, pd, EPS ) =
+ begingroup
+ save _a,_b,_c,_e,_error,_k,_q,w,wa,wb,wc,wd,pn,za, zb, zc, zd;
+ numeric _a,_b,_c,_error,_k,wa,wb,wc,wd;
+ path p,_q,_q[],_e;
+ color pn;
+ pair za, zb, zc, zd;
+ cmykcolor w;
+ za = rp(pa); show za;
+ zb = rp(pb); show zb;
+ zc = rp(pc); show zc;
+ zd = rp(pd); show zd;
+ p = za .. controls zb and zc .. zd;
+ pn = N( f - viewcentr );
+ wa = cdotprod( pn, f-pa ); show wa;
+ wb = cdotprod( pn, f-pb ); show wb;
+ wc = cdotprod( pn, f-pc ); show wc;
+ wd = cdotprod( pn, f-pd ); show wd;
+ w = ( wa, wb, wc, wd );
+ _a:=0;
+ _b:=1;
+ _k:=1/sqrt(2);
+ _q:=(point 0 of p);
+ _q[4]:=nurbstobezier(p,w,_a,_b);
+ forever:
+ exitunless(_a<1);
+ _q[1]:=_q[4];
+ _c:=_b-_k*((_b-_a)**2);
+ _q[2]:=beziertobezier(_q[1],_a,_c);
+ _q[3]:=nurbstobezier(p,w,_a,_c);
+ _q[4]:=_q[3];
+ _e:=((point 0 of _q[2])-(point 0 of _q[3])) ..
+ controls ((postcontrol 0 of _q[2])-(postcontrol 0 of _q[3]))
+ and ((precontrol 1 of _q[2])-(precontrol 1 of _q[3])) ..
+ ((point 1 of _q[2])-(point 1 of _q[3]));
+ _error:=beziernorm(_e)/beziernorm(_q[3]);
+ if (_error > EPS):
+ _b:=_c;
+ else:
+ _q[2]:=beziertobezier(_q[1],_c,_b);
+ _q[3]:=nurbstobezier(p,w,_c,_b);
+ _e:=((point 0 of _q[2])-(point 0 of _q[3])) ..
+ controls ((postcontrol 0 of _q[2])-(postcontrol 0 of _q[3]))
+ and ((precontrol 1 of _q[2])-(precontrol 1 of _q[3])) ..
+ ((point 1 of _q[2])-(point 1 of _q[3]));
+ _error:=beziernorm(_e)/beziernorm(_q[3]);
+ if (_error > EPS):
+ _q:=_q .. controls (postcontrol 0 of _q[4])
+ and (precontrol 1 of _q[4]) .. (point 1 of _q[4]);
+ _a:=_c;
+ _q[4]:=_q[3];
+ else:
+ _q:=_q .. controls (postcontrol 0 of _q[1])
+ and (precontrol 1 of _q[1]) .. (point 1 of _q[1]);
+ _a:=_b;
+ _q[4]:=nurbstobezier(p,w,_a,1);
+ fi;
+ _b:=1;
+ fi;
+ endfor;
+ ( _q )
+ endgroup
+enddef;
+
+beginfig(1);
+ numeric tu, num, i, fac;
+ pen pencontrol, penspline, penalytic;
+ color colcontrol, colspline, colorytic, colormark;
+ color w[];
+ pencontrol = pencircle scaled 4pt;
+ penspline = pencircle scaled 2pt;
+ penalytic = pencircle scaled 1pt;
+ colcontrol = black;
+ colspline = red;
+ colorytic = blue+green;
+ colormark = (0.8,0.8,0.1);
+ tu = 6cm;
+ num = 50;
+ fac = 1.2;
+ transform T;
+ T = identity scaled tu;
+ z21 = origin;
+ z22 = (1,0);
+ z23 = (1,1);
+ z24 = (0,1);
+ z1 = z21 transformed T;
+ z2 = z22 transformed T;
+ z3 = z23 transformed T;
+ z4 = z24 transformed T;
+ z6 = (fac,0) transformed T;
+ z8 = (0,fac) transformed T;
+ drawarrow z1--z6;
+ drawarrow z1--z8;
+ label.lrt( "x", z6 );
+ label.ulft( "y", z8 );
+ dotlabels.urt(1,2,3,4);
+ z11 = twotwo( z21 );
+ z12 = twotwo( z22 );
+ z13 = twotwo( z23 );
+ z14 = twotwo( z24 );
+ w1 = twothr( z21 );
+ w2 = twothr( z22 );
+ w3 = twothr( z23 );
+ w4 = twothr( z24 );
+ cartaxes( fac, fac, 0.3*fac );
+ draw z12--z13--z14 dashed evenly;
+ draw z2--z3--z4 dashed evenly;
+ % 1) Next line: MetaPost intrinsic path.
+ draw z1..controls z2 and z3..z4 withpen penspline withcolor colspline;
+ % 2) Next line: my implementation of the MetaPost intrinsic path.
+ draw z1 for i=1 upto num: ..casteljau(z1,z2,z3,z4,i/num) endfor
+ withpen penalytic withcolor colorytic;
+ % 5) Next line: hopefully, how it should be done. Yeah! Way to go!
+ draw z11 for i=1 upto num: ..nextthirty(w1,w2,w3,w4,i/num) endfor
+ withpen pencontrol withcolor colcontrol;
+ % 4) Next line: MetaPost intrinsic path of perspectived control points.
+ draw z11..controls z12 and z13..z14
+ withpen penspline withcolor colspline;
+ % 3) Next line: what should be drawn in perspective.
+ draw z11 for i=1 upto num: ..rp(casteljau(w1,w2,w3,w4,i/num)) endfor
+ withpen penalytic withcolor colorytic;
+ % 6) Next line: Troy's approximation
+ draw nurbsapprox(w1,w2,w3,w4) withcolor colormark;
+
+
+endfig;
+
+f := 1.05*(3,5,2);
+Spread := 160;
+
+beginfig(2);
+ color w[];
+ w1 = (1,0,0);
+ w2 = (0,0,1);
+ w3 = (0,1,0);
+ w4 = (1,1,1);
+ w5 = (1,1,0);
+ w6 = (1,0,1);
+ w7 = (0,1,1);
+ cartaxes( fac, fac, fac );
+ draw rp(w1)--rp(w2)--rp(w3)--rp(w4)--rp(w5)--rp(w1)--rp(w6)--
+ rp(w2)--rp(w7)--rp(w3)--rp(w5) dashed withdots;
+ draw rp(w6)--rp(w4)--rp(w7) dashed withdots;
+ draw xoy(w1) for i=1 upto num:
+ ..xoy(casteljau(w1,w3,black,w5,i/num))
+ endfor withcolor colormark;
+ draw xoy(w1) for i=1 upto num: ..xoy(casteljau(w1,w2,w3,w4,i/num)) endfor;
+ draw xoz(w1) for i=1 upto num: ..xoz(casteljau(w1,w2,w3,w4,i/num)) endfor;
+ draw rp(w1) for i=1 upto num: ..nextthirty(w1,w2,w3,w4,i/num) endfor
+ withpen pencontrol withcolor colcontrol;
+ draw rp(w1) for i=1 upto num: ..rp(casteljau(w1,w2,w3,w4,i/num)) endfor
+ withpen penalytic withcolor colorytic;
+ draw nextthirty(w1,w2,w3,w4,0.5) withpen pencontrol withcolor red;
+endfig;
+
+beginfig(5);
+ color w[];
+ for i=1 upto 4:
+ w[i]=(uniformdeviate(1),uniformdeviate(1),uniformdeviate(1));
+ draw rp(w[i]) withpen pencontrol;
+ draw xoy(w[i]) withpen penspline;
+ draw xoz(w[i]) withpen penspline;
+ draw yoz(w[i]) withpen penspline;
+ endfor;
+ cartaxes( fac, fac, fac );
+ draw rp(w1)--rp(w2)--rp(w3)--rp(w4) withpen penspline dashed evenly;
+ draw xoy(w1)--xoy(w2)--xoy(w3)--xoy(w4) withpen penalytic dashed evenly;
+ draw xoz(w1)--xoz(w2)--xoz(w3)--xoz(w4) withpen penalytic dashed evenly;
+ draw yoz(w1)--yoz(w2)--yoz(w3)--yoz(w4) withpen penalytic dashed evenly;
+ draw xoy(w1) for i=1 upto num: ..xoy(casteljau(w1,w2,w3,w4,i/num)) endfor;
+ draw xoz(w1) for i=1 upto num: ..xoz(casteljau(w1,w2,w3,w4,i/num)) endfor;
+ draw yoz(w1) for i=1 upto num: ..yoz(casteljau(w1,w2,w3,w4,i/num)) endfor;
+ draw rp(w1) for i=1 upto num: ..nextthirty(w1,w2,w3,w4,i/num) endfor
+ withpen pencontrol withcolor colcontrol;
+ draw fitthreednurbswithtwodbezier(w1,w2,w3,w4,0.005)
+ withpen penalytic withcolor red;
+endfig;
+
+beginfig(3);
+ path xyp, xzp;
+ xyp = origin..tension 2 and 0.75..right...(right+up+right+up);
+ xzp = origin..tension 2 and 0.75..(right+up)...(right+up+right);
+ z1 = postcontrol 0 of xyp;
+ z2 = postcontrol 0 of xzp;
+ z3 = precontrol 1 of xyp;
+ z4 = precontrol 1 of xzp;
+% show (xpart z1);
+% show (xpart z2);
+% show (xpart z3);
+% show (xpart z4);
+ draw xyp;
+ draw xzp;
+endfig;
+
+end.