% 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.