# [metapost] envelope approximation

youping huang yphuang62 at yahoo.com
Fri Mar 25 05:57:57 CET 2005

```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

```