% splineperspective.mp % L. Nobre G. and Troy Henderson % 2007 -- 2012 prologues := 1; defaultfont := "cmss17"; color f, viewcentr; boolean ParallelProj; f := (3,5,4); % This f is the point of view in 3D viewcentr := black; % This is the aim of the view ParallelProj := false; % Kind of perspective % def X(expr A) = if color A: redpart A else: cyanpart A fi enddef; def Y(expr A) = if color A: greenpart A else: magentapart A fi enddef; def Z(expr A) = if color A: bluepart A else: yellowpart A fi enddef; def conorm(expr A) = ( X(A) ++ Y(A) ++ Z(A) ) enddef; def N(expr A) = begingroup save M, exitcolor; numeric M; color exitcolor; M = conorm( A ); if M > 0: exitcolor = ( X(A)/M, Y(A)/M, Z(A)/M ); else: exitcolor := black; fi; ( exitcolor ) endgroup enddef; def cdotprod(expr A, B) = ( X(A)*X(B) + Y(A)*Y(B) + Z(A)*Z(B) ) enddef; def ccrossprod(expr A, B) = ( Y(A)*Z(B) - Z(A)*Y(B), Z(A)*X(B) - X(A)*Z(B), X(A)*Y(B) - Y(A)*X(B) ) enddef; % The dotproduct of two normalized vectors is the cosine of the angle % they form. def ndotprod(expr A, B) = begingroup save a, b; color a, b; a = N(A); b = N(B); ( ( X(a)*X(b) + Y(a)*Y(b) + Z(a)*Z(b) ) ) endgroup enddef; % The normalized crossproduct of two vectors. % Also check getangle below. def ncrossprod(expr A, B) = N( ccrossprod( A, B ) ) enddef; % Haahaa! Trigonometry. def getangle(expr A, B) = begingroup save coss, sine; numeric coss, sine; coss := cdotprod( A, B ); sine := conorm( ccrossprod( A, B ) ); ( angle( coss, sine ) ) endgroup enddef; def rp(expr R) = begingroup save v, u; save verti, horiz, eta, squarf, radio; color v, u; numeric verti, horiz, eta, squarf, radio; v = N( (-Y(f-viewcentr), X(f-viewcentr), 0) ); u = ncrossprod( f-viewcentr, v ); horiz = cdotprod( R-viewcentr, v ); verti = cdotprod( R-viewcentr, u ); if ParallelProj: eta = 1; else: squarf = cdotprod( f-viewcentr, f-viewcentr ); radio = cdotprod( R-viewcentr, f-viewcentr ); eta = 1 - radio/squarf; fi; ( 150*(horiz,verti)/eta ) endgroup enddef; def cartaxes(expr axex, axey, axez) = begingroup save orig, axxc, ayyc, azzc; color orig, axxc, ayyc, azzc; orig = (0,0,0); axxc = (axex,0,0); ayyc = (0,axey,0); azzc = (0,0,axez); drawarrow rp(orig)..rp(axxc); drawarrow rp(orig)..rp(ayyc); drawarrow rp(orig)..rp(azzc); label.bot( "x" ,rp(axxc)); %%%%%%%%%%%%%%%%%%%%%%%%% label.bot( "y" ,rp(ayyc)); %% Some Labels... %% label.lft( "z" ,rp(azzc)); %%%%%%%%%%%%%%%%%%%%%%%%% endgroup enddef; def line( expr Ang ) = begingroup numeric a, b, c; a = (2-(1 ++ cosd(Ang))*cosd(3*Ang))*cosd(Ang); b = (2-(1 ++ cosd(Ang))*cosd(3*Ang))*sind(Ang); c =1.5+(1 ++ cosd(Ang))*sind(3*Ang); ( (a,b,c) ) endgroup enddef; % 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; % The code below is a development by Troy from an original by Przemek Koprowski vardef rationalbezier(expr A,B,C,D) = begingroup save P,Q,E,a,b,c,d,r,EPS; color P[],Q[],E; pair a,b,c,d; EPS:=1/10; a:=(redpart A,greenpart A)/(bluepart A); b:=(redpart B,greenpart B)/(bluepart B); c:=(redpart C,greenpart C)/(bluepart C); d:=(redpart D,greenpart D)/(bluepart D); r:=max(abs(a-b),abs(a-c),abs(a-d),abs(b-c),abs(b-d),abs(c-d)); if (max(abs(b-1/2[a,c]),abs(c-1/2[b,d])) > EPS*r): P[0]:=A; P[1]:=1/2[A,B]; E:=1/2[B,C]; Q[2]:=1/2[C,D]; Q[3]:=D; P[2]:=1/2[P[1],E]; Q[1]:=1/2[E,Q[2]]; P[3]:=1/2[P[2],Q[1]]; Q[0]:=P[3]; rationalbezier(P[0],P[1],P[2],P[3]) & rationalbezier(Q[0],Q[1],Q[2],Q[3]) else: a .. controls b and c .. d fi endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Here's where the fun begins % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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; 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(0); draw for i=1 step 6 until 360: xoy(line(i)).. endfor cycle; draw for i=1 step 6 until 360: rp(line(i)).. endfor cycle; for i=1 step 15 until 360: draw rp(line(i)) withpen pencircle scaled 2mm; endfor; endfig; 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.067; % 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f := 0.35*(3,5,2); 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); 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; %%%%%%%%%%%%% the following line may crash; %% draw fitthreednurbswithtwodbezier(w1,w2,w3,w4,0.2) %% 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.