[metapost] envelope approximation

[youping huang]

Hi all,

I have been working on numerical solutions of the 
envelope of a path. My metapost program with a brief
explanation, and a few examples are included. The 
solutions look quite good.

It can be considered as a response to the challenge 
posted by Laurent Siebenmann to tex-fonts on May 6, 
2003.

Youping Huang

**********

Consider drawing a path with an elliptic pen in 
METAPOST, we are interested in approximating the 
envelope with one or two bezier cubics. We assume the 
path is well-behaved.

Let p be the path, Env the envelope, and env the 
approximation of Env.

Because an elliptic pen can be defined as a circular 
pen transformed by a scaling in y direction and a 
rotation, we need to consider only the case of 
circular pen.

In order to make env smooth when p consists of 
smoothly cancatenated bezier cubics, the tangent of 
env and the tangent of p should have the same 
direction at both ends.

Therefore, if it is a bezier cubic, env will have the 
following form.

(z_0+t_0)..controls (a*t_0 shifted z_0+n_0) and 
 (-b*t_1 shifted z_1+n_1)..(z_1+n_1)

where, for i=0 and 1,
  z_i = point i of p
  t_i = direction i of p
  n_i = radius * unitvector (t_i rotated -90)
  a and b = unknown parameters

The error of our approximation is defined as the 
maximum difference between Env and env along normal 
lines of p. In order to calculate the error, pick nsec

points which are equally spaced in time along p. At 
each of these nsec-2 points (excluding endpoints) find

the normal line of p, its intersections with Env and 
env, and the distance between thess two intersection 
points. The approximate error between Env and env is 
calculated as the maximum of these nsec-2 distances.

Downhill simplex method is used to find the optimal a 
and b.

If the error of using a bezier cubic as env is too 
big, then cut p at one of these nsec-2 points, and use

two bezier cubics as env. Search over these nsec-2 
points to find the best env.

**********
prologues:=1;

numeric 
 rpen, % pen radius, major axis
 ysc,  % (minor axis/major axis)
 alpha,% pen angle
 mr,   % ratio of half length of intersecting normal
       % lines to pen radius
 nsec, % number of intersecting normal lines
 tol,  % tolerance for simplex downhill algorithm
 simf[];
pair 
 simx[]; % 0=sum, 1,2,3=simplex, 7=candidate, 8=save,
         % 9=amotry, % 31-33, 41-43, 51-53, 
         % 61-63=starting points, 90-92=best fit
boolean showpen, showfitone, showfittwo;

def showinfo (expr start, ind) =
 show "["& decimal start &"] "& decimal iter &" ("&
  decimal xpart simx[ind] &","& 
  decimal ypart simx[ind] &") -> "&
  decimal simf[ind];
enddef; % showinfo

% 90=best
def fitone (expr isec) =
 simf90 := 1000;
 for start:=30,40,50,60:
  temp := amoeba(afunc) (isec, start, tol);
  if showfitone: showinfo (start, temp); fi
  if simf[temp]<simf90: 
   simx90:=simx[temp]; simf90:=simf[temp]; 
  fi
 endfor
enddef; % fitone

% returns split index, 91,92=best
vardef fittwo =
 save best, diff, cut;

 best = 1000;
 for k:=1 upto nsec-1:
  fitone (k); diff:=simf90; simx93:=simx90;
  fitone (-k);
  if showfittwo: 
   show decimal k &"> "& decimal diff &" "& 
    decimal simf90; 
  fi
  if diff<simf90: diff:=simf90; fi
  if best>diff: 
   best:=diff; cut:=k; 
   simx91:=simx93; simf91:=diff; 
   simx92:=simx90; simf92:=simf90;
  fi
 endfor
 cut
enddef; % fittwo

vardef envelope (expr curve) = 
 save temp;
 path curve_tr, nn[];
 numeric time, temp, cutsec;
 pair t_[], n_[], pt_[], pn_[];
 transform tr_;

 tr_ := identity yscaled ysc rotated alpha;
 pickup pencircle scaled 2rpen transformed tr_;
 draw curve withcolor 0.8(1,1,1);
 pickup pencircle scaled 1;
 draw curve;

 curve_tr = curve transformed inverse tr_;
 for i:= 0 upto nsec:
  time := i/nsec;
  pt_[i] := point time of curve_tr;
  t_[i] := direction time of curve_tr;
  n_[i] := rpen*unitvector (t_[i] rotated -90);
  pn_[i] := pt_[i]+n_[i];
  nn[i] := (-mr*n_[i]--mr*n_[i]) shifted
   (pt_[i]+n_[i]);
  if showpen: 
   draw ((fullcircle scaled 2rpen shifted pt_[i])
    transformed tr_) withcolor (0,1,1); 
  fi
 endfor

 showfitone:=true;
 fitone (0);
 show "fitone: "& decimal simf90;
 draw candidate(0, 90) transformed tr_ withcolor red;

 if simf90>1:
  showfitone:=false;
  cutsec:=fittwo;
  show "fittwo: "& decimal cutsec &"> "& 
   decimal simf91 &"~"& decimal simf92;
  draw candidate(cutsec, 91) transformed tr_ 
   withcolor blue;
  draw candidate(-cutsec, 92) transformed tr_ 
   withcolor green;
 fi
enddef; % envelope

% returns a transformed curve
% 0=whole, positive=part1, negative=part2
vardef candidate (expr k, ind) = 
 save temp;

 if k=0: temp=nsec; else: temp=abs(k); fi
 if k>=0:
  pn_0..controls ((xpart simx[ind]*t_0) shifted pn_0) 
   and ((-ypart simx[ind]*t_[temp]) shifted pn_[temp])
   ..pn_[temp]
 else:
  pn_[temp]..controls ((xpart simx[ind]*t_[temp])
   shifted pn_[temp]) and 
   ((-ypart simx[ind]*t_[nsec]) shifted pn_[nsec])
   ..pn_[nsec]
  fi
enddef; % candidate

% returns approximate (transformed) error between 
% envelope and candidate
vardef afunc (expr k, ind) = 
 save diff;
 path q;
 numeric err, diff, ibeg, iend;

 q = candidate (k, ind);
 if k=0: ibeg=1; iend=nsec-1;
 elseif k>0: ibeg=1; iend=k-1;
 else: ibeg=-k+1; iend=nsec-1; fi

 err := 0;
 for i:=ibeg upto iend:
  diff := abs(0.5-xpart(nn[i] intersectiontimes q)); 
  if diff>err: err:=diff; fi
 endfor
 err := 2rpen*mr*err;
 iter:=iter+1;
 err
enddef; % afunc

vardef amotry (suffix $) (expr k, fac) =
 numeric faca, facb, ytry;

 faca = (1-fac)/2;
 facb = faca-fac;
 simx7 := faca*simx0-facb*simx[ihi];
 ytry = $(k, 7);
 if ytry<simf[ihi]:
  simf[ihi] := ytry;
  simx0 := simx1+simx2+simx3;
  simx[ihi] := simx7;
 fi
 ytry
enddef; % amotry

% adapted from numerical recipes (10.4)
% returns index ilo
vardef amoeba (suffix $) (expr k, start, atol) =
 numeric ihi, inhi, ilo, iter, sc;
 boolean calcsum;

 if k=0:     sc=1;
 elseif k>0: sc=k/nsec;
 else:       sc=(nsec+k)/nsec; fi
 for i:=1 upto 3:
  simx[i] := sc*simx[start+i];
  simf[i] := $(k, i);
 endfor

 iter = 0;
 calcsum = true;
 forever:
  if calcsum: simx0 := simx1+simx2+simx3; fi

  if simf1>simf2: ihi:=1; inhi:=ilo:=2;
  else: ihi:=2; inhi:=ilo:=1; fi
  if simf3<=simf[ilo]: ilo:=3; fi
  if simf3>simf[ihi]: inhi:=ihi; ihi:=3;
  elseif simf3>simf[inhi]: inhi:=3; fi
  exitif simf[ihi]-simf[ilo]<atol;

  simf9 := amotry($) (k, -1); %%% reflect
  calcsum := false;
  if simf9<simf[ilo]: simf9 := amotry($) (k, 2);
   %%% expand
  elseif simf9>=simf[inhi]: 
   simf8:=simf[ihi];
   simf9 := amotry($) (k, 0.5); %%% contract
   if simf9>=simf8: %%% multiple contraction
    for i:=ihi,inhi: 
     simx[i] := 0.5[simx[i],simx[ilo]];
     simf[i] := $(k, i);
    endfor
    calcsum := true;
   fi
  fi
  exitif iter>1000;
 endfor
 ilo
enddef; % amoeba

 mr = 4; 
 nsec=40; 
 tol = 0.0001;
 simx31 = (1.2,1.2);
 simx32 = (1.0,0.8); 
 simx33 = (0.8,1.0);
 simx41 = (1.1,1.1); 
 simx42 = (1.0,0.7); 
 simx43 = (0.7,1.0);
 simx51 = (1.1,1.2); 
 simx52 = (1.0,0.7); 
 simx53 = (0.8,1.0);
 simx61 = (1.2,1.1); 
 simx62 = (1.0,0.8); 
 simx63 = (0.7,1.0);
 showpen = false;
 showfittwo = true; 

path p;

beginfig(1);
 rpen := 50;
 ysc := 0.3;
 alpha := 30;
 z1=(0,0); z2=(20,150); z3=(100,200); z4=(300,200);
 p := z1..controls z2 and z3..z4;
 envelope(p);
endfig;

beginfig(2);
 rpen := 35;
 ysc := 0.3;
 alpha := 30;
 z1=(0,0); z2=(0,250); z3=(0,250); z4=(350,250);
 p := z1..controls z2 and z3..z4;
 envelope(p);
endfig;

beginfig(3);
 rpen:=30;
 ysc := 0.4;
 z1=(0,0); z2=(140,300); z3=(250,-100); z4=(300,200);
 p := z1..controls z2 and z3..z4;
 envelope(p);
endfig;

% Lars Hellstr\"om's example
beginfig(4);
 z1=(50,300); z2=(50,150); z3=(150,50); z4=(300,50);
 p := z1..controls z2 and z3..z4;

 transform T, Tinv;
 (0,0) transformed T = (0,0);
 (0.5,0) transformed T = (20,5);
 (0,0.5) transformed T = (0,10);
 Tinv = inverse T;
 rsc = 1/(abs(xxpart Tinv)+abs(yypart Tinv));
 show rsc;
 Tinv := (inverse T) scaled 2rsc;
 show "choose sc to have moderate numbers in..";
 show Tinv;

 a = (xxpart Tinv);
 b = (xypart Tinv);
 c = (yxpart Tinv);
 d = (yypart Tinv);

 alpha := 0.5*angle(a*a-b*b+c*c-d*d,2*(a*b+c*d));
 sina = sind(alpha);
 cosa = cosd(alpha);
 xs = cosa*cosa*(a*a+c*c)+2*sina*cosa*(a*b+c*d)
  +sina*sina*(b*b+d*d);
 ys = sina*sina*(a*a+c*c)-2*sina*cosa*(a*b+c*d)
  +cosa*cosa*(b*b+d*d);
 if xs<ys: 
  rpen := sqrt(1/xs); 
  ysc := sqrt(xs/ys);
 else:     
  rpen := sqrt(1/ys); 
  ysc := sqrt(ys/xs);     alpha:=alpha+90;
 fi
 show rpen; show ysc; show alpha;
 rpen := rsc*rpen;
 showpen:=true;
 envelope(p);

% double check tr_
 draw (fullcircle transformed T shifted (200,200))
  withcolor red;
 draw fullcircle scaled 2rpen transformed tr_ 
  shifted (200,200) dashed evenly withcolor blue;
endfig;
end


__________________________________________________
Do You Yahoo!?
Tired of spam?  Yahoo! Mail has the best spam protection around 
http://mail.yahoo.com