% nurbstobezier.mp
% Troy Henderson
% 2007

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(0);
% 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;

end