summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/metafont/roex
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2006-01-12 23:43:37 +0000
committerKarl Berry <karl@freefriends.org>2006-01-12 23:43:37 +0000
commit9199479b7f615c58da260c7ded57b63aed5cf799 (patch)
tree2a9a4068b7a757edd5db10b94946b6c90404b4d6 /Master/texmf-dist/metafont/roex
parent2828582f70c71954ec767bead50118e55e23ceae (diff)
trunk/Master/texmf-dist/metafont
git-svn-id: svn://tug.org/texlive/trunk@604 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/metafont/roex')
-rw-r--r--Master/texmf-dist/metafont/roex/roex.mf1756
1 files changed, 1756 insertions, 0 deletions
diff --git a/Master/texmf-dist/metafont/roex/roex.mf b/Master/texmf-dist/metafont/roex/roex.mf
new file mode 100644
index 00000000000..dbf6f8e35dc
--- /dev/null
+++ b/Master/texmf-dist/metafont/roex/roex.mf
@@ -0,0 +1,1756 @@
+%%%%
+%%%% This file belongs to the ROEX package.
+%%%%
+%%%% ---------------------------------------------------------------------
+%%%% MFT formatting commands
+%%%% ---------------------------------------------------------------------
+%%% length quicksort
+%%% length cycle zang pos_turn neg_turn
+%%% good enc
+%%% labels makelabel
+%%% length make_cycle make_join make_cyclic_join make_end make_edge
+%%% point predir postdir upredir upostdir udir
+%%% dotprod det
+%%% subpath pos_subpath neg_subpath
+%%% message info_ro info_es
+%%% draw roex_default
+%%% -- &&
+%%%% ---------------------------------------------------------------------
+%%%% \TeX formatting commands
+%%%% ---------------------------------------------------------------------
+%%\vsize245mm
+%%\font\titfnt cmtt10 at 48 pt
+%%{\let\makefootline\empty \let\makeheadline\empty
+%%\vglue0ptplus1fill
+%%\centerline{\titfnt ROEX.MF}
+%%\bigskip
+%%\centerline{ver. 0.56 (Wednesday, October 25th, 1995)}
+%%\vfill\vfill\eject}
+%% % ---
+%% \vsize 245mm
+%% % an innocent formatting trick: the underscore character ending a name
+%% % will be typeset as an superscript asterisk
+%% \let\oriunderscore\_
+%% \newif\ifbgroupopen\bgroupopenfalse
+%% \def\altdblbackslash#1{\bgroup\bgroupopentrue\it#1}
+%% \def\optegroup{\ifbgroupopen\egroup\fi}
+%% \def\underscoreasasterisk#1{%
+%% \ifx#1\relax\optegroup^*\else\oriunderscore#1\fi}
+%% \def\\#1{%
+%% \let\_\underscoreasasterisk
+%% \altdblbackslash{#1\relax}\optegroup
+%% \let\_\oriunderscore}
+%% % ---
+%% \def\dblhyph{--}
+%% \def\8#1{\def\eightparm{#1}\mathrel{\mathcode`\.="8000 \mathcode`\-="8000
+%% \ifx\eightparm\dblhyph\setbox\shorthyf\hbox{\bf -\kern-.05em}\fi%
+%% #1\unkern}} % `..' and `--'
+%% % ---
+%% \def\MP{{\tenlogo META}\-{\tenlogo POST}}
+% ------------------------------------------------------------------------
+% This is ROEX.MF file containing \MF definitions implementing
+% operations known as `remove overlap' and `expand stroke'.
+% ------------------------------------------------------------------------
+% Authors: \bf{}B. Jackowski, P. Pianowski, M. Ry\'cko \& S. Soko\l{}owski
+% ------------------------------------------------------------------------
+% H I S T O R Y
+% ver. 0.1 (1 / 9 VI 1994):
+% * incunabula version
+% ver. 0.5 (15 VIII / 1 IX 1995):
+% * pioneer version, released during the 9th Euro\TeX conference in Arnhem
+% ver. 0.55 (26 IX 1995):
+% * if a single path is an argument to |remove_overlap|, removing
+% of self-overlaps is performed, hence several adjustments, most
+% significant changes were introduced in |is_far_enough| and
+% |intersect_two_segments|; this ismprovement is, in fact, a prelude
+% towards a more general approach
+% * a bug trap added in |clean_path|
+% * positioning of labels not forced in |mark_nodes|
+% * |quicksort| more flexible
+% * more statistics available (optionally) in |find_minimal_secant|
+% * displaying information changed
+% * \TeX formatting comments collected at the beginning of the file
+% * a result of mental laps corrected in |build_node_structure|
+% (minimal secant has nothing to do with minimal distance between nodes)
+% * a silly bug removed in |prepare_input_data| (|W| instead of |W_|)
+% This version was released during the CyrTUG-95 meeting in Moscow
+% ver. 0.56 (27 X 1995):
+% * comments adjusted to a new distribution
+% * the name |miter_limit| changed to |miter_size| in order to avoid
+% misunderstanding, as in this implementation it is a dimen, while
+% in PostScript it is a dimensionless quantity
+% ------------------------------------------------------------------------
+% S Y N O P S I S
+% ------------------------------------------------------------------------
+%
+% Such operations as `remove overlap' and `expand stroke' are perhaps
+% particularly useful in the contex of exporting data from \MF to other
+% languages, e.g., to PostScript or HP-GL. Therefore the file ROEX.MF has
+% been included into the MFTOEPS package (which accomplishes export from \MF
+% to Encapsulated PostScript), although it can be used by ``normal'' \MF
+% users, too. Therefore our favourite macros (e.g., |pos_turn|, |neg_turn|,
+% |make_list|), are defined identically in both ROEX and MFTOEPS.
+%
+% We hope that tiny adjustments, if any, should be sufficient for transforming
+% the macros to the form accepted by both \MF and \MP.
+%
+% Sample \MF programs (i.e., simple examples) illustrating the use of the
+% ``interface'' macros, i.e., |remove_overlap|, |expand_stroke| and
+% |change_weight|, can be found in a subdirectory ROEXSAMP. It is instructive
+% to generate EPS files and then to play around with the results using
+% CorelDRAW! or Adobe Illustrator.
+%
+%
+% REMOVING OVERLAPS
+%
+% The command |remove_overlap_| requires three parameters. The first
+% parameter is a list of paths to be processed; the paths are assumed to
+% have a non-zero |turningnumber| and no self intersecting points (no
+% checking is performed, except that non-cyclic paths are ignored).
+% The second parameter is a list (possibly empty) of weights assigned
+% to paths; more exactly, it is a list of pairs |(i,w.i)|, where |i| is
+% the order number of a path and |w.i| is the respective weight.
+% If the weight is not specified it is assumed to be equal to |1|.
+% The last parameter is a suffix, i.e., the name of the resulting data
+% structure; given a suffix is |R|, |R.num| is the number of the resulting
+% paths, and |R1|, |R2|, ..., |R[R.num]| are the paths. If the suffix
+% contains an index, e.g., |P[x]q|, the user is responsible for providing
+% appropriate declarations prior to calling |remove_overlap|,
+% in this case: |numeric P[\\]q.num; path P[\\]q[\\]|. If a variable
+% |append_results| is assigned a definite value (by default it is undefined),
+% |R.num| is not zeroed at the stage of initialisation, thus the results
+% are accumulated (see example ROES-02.MF in the subdirectory ROEXSAMP).
+%
+% The algorithm assumes that a path |p| fills its interior with the colour
+% |w*turningnumber(p)|, where |w| is the weight assigned to |p|. If an
+% area is coloured by several paths, the colors are summed up. The user
+% decides which areas are the resulting ones. By default, these are areas
+% which have the interior painted with a colour $\ge1$ and the exterior
+% painted with a colour $\le0$. There is a two-parameter function
+% (parameters are numbers) that can be redefined by a user, |good_colors|,
+% which governs the decision. The user is responsible for a proper definition
+% of this function (the formula |good_colors(m,n) and good_colors(n,m)| must
+% be false; cf. the default definition of |good_colors| at the end of this
+% file). There is also a variable |background_color|, by default equal to |0|,
+% which determines the colour of the Euclidian plane. One more function that
+% is meant to be redefined by the user, if needed, is one-parameter function
+% |touch_path|; the function is applied to every input path at the stage of
+% initialisation, and can be used, e.g., for adjusting the direction of paths
+% (cf. example RO-04.MF in the subdirectory ROEXSAMP).
+%
+% The orientation of paths generated by the |remove_overlap| macro is defined
+% in such a way that in order to fill the resulting figure the internal
+% variable |turningcheck| should be set to zero prior to using the |fill|
+% command.
+%
+% Examples:
+% Assume that paths |A|, |B| and |C| are defined as follow (say, |w=h=1cm|):
+% |A=unitsquare xscaled 1/5w yscaled h shifted (2/5w,0);|
+% |B=A rotatedaround((1/2w,1/2h), 60);|
+% |C=B rotatedaround((1/2w,1/2h), 60);|
+% Calling
+% |remove_overlap (A,B,C) () R;|
+% will result in generating a single path |R1| (|R.num=1|) of a six-arm
+% propeller shape. Adding one more path:
+% |D=reverse fullcircle scaled 3/4w shifted (1/2w,1/2h);|
+% and calling
+% |remove_overlap (A,B,C,D) ((4,2)) R;|
+% (|D| has weight |=2|) will result in generating seven objects: six
+% ``tips'' of a propeler and a regular hexagon in the center. Try to guess
+% what would be the result of
+% |remove_overlap (A,B,C,D) () R;|
+% Not all paths need to intersect. For example, the following set of paths
+% |A=fullcircle scaled w shifted (1/2w,1/2h);|
+% |B=reverse unitsquare xscaled 1/5w yscaled 3/5w shifted (2/5w,1/5w);|
+% |C=B rotatedaround((1/2w,1/2h), 90);|
+% after calling
+% |remove_overlap (A,B,C) () R;|
+% will yield a circle surrounding a white cross. Since the orientation
+% of the resulting paths is important here, the |fill| commands should be
+% preceded by
+% |interim turningcheck:=0;|
+% assignment.
+%
+%
+% EXPANDING STROKES
+%
+% Expanding stroke means finding the trace of the outline of an imaginary pen
+% moving along a path. There are two commands accomplishing expanding stroke:
+% |expand_stroke| and |change_weight|. Both make use of the essentially the
+% same algorithm, except that the latter finds only one edge and ignores
+% non-cyclic paths. Both commands require three parameters: first and third
+% are analogous to the parameters of the |expand_stroke| macro (see above),
+% the second denotes the radius (not diameter) of the circular pen.
+% The algorithm works in such a way that the result of the |expand_stroke|
+% does not depend on the direction of a path for cyclic paths, namely,
+% the outer edge is always positively directed and the inner is negatively
+% directed, provided the radius is positive; if the radius is negative,
+% the outer edge is negatively directed and the inner one---positively.
+% For non-cyclic paths positive radius yields positively directed resulting
+% paths, negative radius---negatively oriented paths. Although the macro
+% |change_weight| is subdued to the same rules, the result depends both
+% on the direction of a path and on the sign of a radius. Let |t| and |r|
+% denote the turning number and the radius, respectively; there are
+% four cases:
+% 1) |t>0| and |r>0|: the resulting path is an outer edge positively
+% directed,
+% 2) |t>0| and |r<0|: the resulting path is an inner edge positively
+% directed,
+% 3) |t<0| and |r>0|: the resulting path is an inner edge negatively
+% directed,
+% 4) |t<0| and |r<0|: the resulting path is an outer edge negatively
+% directed.
+% Following PostScript, we introduced three variables which govern the shape
+% of joins and ends: |join_kind| (corresponds to |setlinejoin|), |end_kind|
+% (corresponds to |setlinecap|) and |miter_size| (corresponds, as the name
+% suggest, to |setmiterlimit|; however, here |miter_size| is a dimen,
+% while in PostScript miter limit is a dimensionless quantity).
+% Currently both |join_kind| and |end_kind| can receive value |0| or |1|,
+% while in PostScript value |2| is also admissible. (The latter option
+% will perhaps be included also into the ROEX package some day, but more
+% tempting is the implementation of extrapolated non-linear joins.)
+% Since the interpretation of |miter_size| (dimen) is slightly different than
+% the interpretation of |miter_limit| (a number), |miter_size| must merely
+% be non-negative, while |miter_limit| must be greater than or equal to $1$.
+% Roughly speaking, value |0| for |join_kind| and |end_kind| denotes cusp
+% joins, cut if necessary at miter limit; value |1| denotes rounded rounded
+% joins (for details see, e.g., ``PostScript Language Reference Manual,''
+% second ed., Addison-Wesley Publishing Company, Ltd.).
+%
+% Example:
+% Assume that a path |A| is simply a square (say, |w=h=1cm|):
+% |A=unitsquare scaled w;|
+% After calling
+% |expand_stroke(A)(1mm)R;|
+% |R1| is a positively directed square of side |12mm|, and |R2| is
+% negatively directed square of side |8mm|.
+% ------------------------------------------------------------------------
+% C A V E A T S , H I N T S A N D C O M M E N T S
+% * The employed algorithms expect that the results are well defined;
+% if the data are weird (e.g., self-loooping path are supplied)
+% the results, if any, may be weird as well.
+% * The case of curves partially overlapping is not handled and, frankly
+% speaking, we have no idea how to implement it efficiently and robustly;
+% if there are such pairs of paths in the input data, the algorithm almost
+% certainly will not produce good results.
+% * Only circular pens are implemented so far.
+% * Be aware of rounding errors, they may cause unpredictable results;
+% in some cases increasing accuracy by using a higher resolution may
+% help, but more adequate seems to be preparing better data (cf. the
+% program RO-07.MF in the subdirectory ROEXSAMP).
+% * Comments in the code are meant primarily for the authors; the user
+% is kindly requested not to complain fiercely if they are of a little
+% use to her/him.
+% * Unfortunately, \MF has no error-handling facility, hence a lot of
+% ``bug traps'' can be found in the code; messages issued in the case
+% of falling into such a trap are rather useless if you don't know the
+% details of the algorithm; this part of the program is certainly to be
+% improved; usually the error help says ``Better stop now! Algorithm
+% failed'' and this advice should be followed; in practice this means that
+% \MF is not able to recognize the details of the picture because of
+% very close nodes (intersection points).
+% * Usually, the first stage of removing overlaps (finding all intersection
+% points) is the longest one, the more segments paths have the longer it
+% lasts; a pity that \MF has no built-in function informing about all
+% intersection points/times of two B\'ezier curves.
+% * Improper definition of |good_colors| may result in erroneous behaviour of
+% the algorithm.
+% * One peculiar case is considered by the expanding stroke algorithm,
+% namely cyclic path of length 2; some more cases might have been taken
+% into account...
+% * There remain a lot of unsolved problems with numerical instability
+% connected with detecting tangent and close points.
+% * In The \MF{}book, p. 229, D. E. Knuth writes:
+% ``...tiny little loops won't hurt anything if you are filling cycles
+% in the correct direction.''
+% Cf. also preceding dangerous band paragraph and exercise on pp. 228--229.
+% ROEX does much more complex things with paths than merely filling them,
+% hence tiny loops may cause some mess, the more so as the built-in
+% function |turningnumber| is very sensitive to such loops, e.g., it may
+% happen that |turningnumber(p)=1| while |turningnumber(reverse p)=0|
+% (cf. example RO-6.MF in the subdirectory ROEXSAMP); hence a hopefully
+% more robust (from the point of view of this application) function
+% is used, |check_turn|, which makes use of \MF's |fill| operation.
+% * There remain several suboptimal algorithms employed, partially on
+% purpose: less efficient algorithms are usually (although not necessarily)
+% more comprehensible and flexible (easier to modify), which is important
+% at the stage of developing a program.
+% * Parameters that may have influence on the process of removing
+% overlaps are |epsil.time|, |epsil.ang| (in degrees), |epsil.dist| (in
+% resolution-dependent units), and |max_idx|; the choice of good default
+% values will need some practice.
+% * Incompatible modifications may come, although we shall do our best
+% to avoid them.
+% ------------------------------------------------------------------------
+% We follow the naming convention of The \MF{}book:
+% ``Private tokens always end with the underscore character.''
+% Since the underscore is a rather illegible character, in a ``neat''
+% printing (using MFT utility) it will appear as an superscript asterisk.
+% ------------------------------------------------------------------------
+% D E F I N I T I O N S
+% ------------------------------------------------------------------------
+% UNIVERSAL MACROS:
+% ---
+% Without the following redefinition:
+def -- = {curl 1} .. tension (1+eps) .. {curl 1} enddef;
+% the result of |p intersectiontimes reverse p|, where |p=(a,b)--(a+3c,b+3d)|,
+% |a|, |b|, |c|, |d| are arbitrary (sic!) \MF's numbers, yields the result
+% |(1/2,1/2)|, which contradicts the statement preceding the exercise 14.17
+% on the page 137 of The \MF{}book. Since it is no longer a ``standard''
+% macro, its formatting is slightly modified.
+% ---
+%%% length ]]] ]]]]
+%%% ; ]
+def ]]] = ] ] ] enddef;
+def ]]]] = ] ] ] ] enddef; % right brackets should be loners, indeed
+%%% ) ] ]] ]]] ]]]]
+% ---
+vardef distance(expr za,zb) = length(za-zb) enddef; % in fact, an alias
+% ---
+vardef interval(expr ta,tb,p) =
+ save ta_,tb_;
+ if cycle p:
+ ta_:=ta mod length(p); tb_:=tb mod length(p);
+ min(length(p)-abs(ta_-tb_), abs(ta_-tb_))
+ else:
+ ta_:=min(max(0,ta),length(p)); tb_:=min(max(0,tb),length(p));
+ abs(ta_-tb_)
+ fi
+enddef;
+% ---
+def make_list(expr k,l) suffix s =
+ for i_:=k upto l: if i_>k: , fi \\ s[i_] endfor
+enddef;
+% ---
+vardef dec_pair(expr z) =
+ "(" & decimal(xpart z) & "," & decimal(ypart z) & ")"
+enddef;
+% ---
+primarydef u det v = % dual operation to |dotprod|
+ (xpart u * ypart v - xpart v * ypart u)
+enddef;
+% ---
+vardef zang(expr u,v) = % useful during testing
+% computes the angle form |u| to |v| (useful for testing)
+ angle(u dotprod v,u det v) mod 360 % CAVEAT! rounding errors
+enddef;
+% ---
+vardef turn_ang(expr za,zb) = % more robust version of |zang|
+% The idea of computing the turn angle is based on the following observation:
+% |z reflectedabout (origin,right)=1/z| for a complex number |z| such that
+% |abs(z)=1|; recall also that multiplication of complex numbers
+% (|zscaled| operation) implies addition of their angle arguments.
+ if (abs(za)>=epsil.len) and (abs(zb)>=epsil.len): % |eps| may be not enough
+ angle(unitvector(za) zscaled (unitvector(zb) reflectedabout (origin,right)))
+ else: whatever fi
+enddef;
+% ---
+def predir expr t of p = ((point t of p)-(precontrol t of p)) enddef;
+def postdir expr t of p = ((postcontrol t of p)-(point t of p)) enddef;
+def udir expr t of p = unitvector(direction t of p) enddef;
+def upredir expr t of p = unitvector(predir t of p) enddef;
+def upostdir expr t of p = unitvector(postdir t of p) enddef;
+% ---
+vardef pos_turn primary p =
+ interim autorounding:=0;
+ if check_turn(p)=0: show p;
+ errhelp "I will leave the path intact, continue with crossed fingers.";
+ errmessage "Cannot make positive turn (check_turn=0)";
+ elseif check_turn(p)<0: reverse fi \\ p
+enddef;
+% ---
+vardef neg_turn primary p =
+ interim autorounding:=0;
+ if check_turn(p)=0: show p;
+ errhelp "I will leave the path intact, continue with crossed fingers.";
+ errmessage "Cannot make negative turn (check_turn=0)";
+ elseif check_turn(p)>0: reverse fi \\ p
+enddef;
+% ---
+vardef check_turn primary p = % seems more adequate than |turningnumber|
+% |epsilon|=|totalweight currentpicture| after |fill unitsquare|,
+% |eps/epsilon=32|, i.e., we admit accuracy of 32 pixels (isn't it too many?)
+ save r_,currentpicture; picture currentpicture;
+ interim turningcheck:=0; interim autorounding:=0;
+ currentpicture:=nullpicture; fill p; r_:=totalweight(currentpicture);
+ if r_>eps: 1 elseif r_<-eps: -1 else: turningnumber(p) fi
+enddef;
+% ---
+def check_embedding(expr a,b)(suffix res) =
+begingroup
+% see comment in |check_turn|
+ save napb_,panb_,currentpicture; picture currentpicture;
+ interim turningcheck:=0; interim autorounding:=0;
+ currentpicture:=nullpicture; fill pos_turn a; fill neg_turn b; cullit;
+ panb_:=totalweight currentpicture;
+ currentpicture:=nullpicture; fill neg_turn a; fill pos_turn b; cullit;
+ napb_:=totalweight currentpicture;
+ if (panb_<eps) and (napb_<>0): res:=1; % $a \subset b$
+ elseif (panb_<>0) and (napb_<eps): res:=2; % $b \subset a$
+ else: res:=0; fi % undefined result
+endgroup
+enddef;
+% ---
+vardef pos_subpath expr z of p =
+ if not cycle p: subpath z of p
+ else:
+ if xpart(z)<=ypart(z): subpath z of p
+ else: subpath (xpart(z),ypart(z)+length(p)) of p
+ fi
+ fi
+enddef;
+% ---
+vardef neg_subpath expr z of p =
+ if not cycle p: subpath z of p
+ else: reverse(pos_subpath (ypart z,xpart z) of p) fi
+enddef;
+% ---
+tertiarydef p && q = % |length(p)>0|
+ (subpath(0,length(p)-1) of p) ..
+ controls (postcontrol length(p)-1 of p) and (precontrol length(p) of p)
+ .. q
+enddef;
+% ---
+def make_cycle expr p = % |length(p)>0|
+ (subpath(0,length(p)-1) of p) ..
+ controls (postcontrol length(p)-1 of p) and (precontrol length(p) of p)
+ .. cycle
+enddef;
+% ---
+vardef is_line(expr B) =
+% checks if a B\'ezier segment |B| is an almost straight line;
+% recall that |z reflectedabout (origin,right)=1/z| for a complex
+% number |z| such that |length(z)=1|; recall also that the multiplication
+% of complex numbers (|zscale| operation) implies the addition of
+% their angle arguments
+ save pa_,pb_,pc_,pd_,ba_,da_,dc_; pair pa_,pb_,pc_,pd_,ba_,da_,dc_;
+ pa_:=point 0 of B; pd_:=point 1 of B;
+ if distance(pa_,pd_)<epsil.dist:
+ false % either really not a line or an uncertain situation (rounding errors)
+ else:
+ da_=unitvector(pd_-pa_) reflectedabout (origin,right);
+ pb_:=postcontrol 0 of B; if distance(pa_,pb_)<epsil.dist: pb_:=pa_; fi
+ pc_:=precontrol 1 of B; if distance(pd_,pc_)<epsil.dist: pc_:=pd_; fi
+ if (pa_=pb_) and (pc_=pd_): true
+ elseif (pa_=pb_):
+ dc_=unitvector(pd_-pc_); abs(angle(dc_ zscaled da_))<epsil.ang
+ elseif (pc_=pd_):
+ ba_=unitvector(pb_-pa_); abs(angle(ba_ zscaled da_))<epsil.ang
+ else:
+ ba_=unitvector(pb_-pa_); dc_=unitvector(pd_-pc_);
+ (abs(angle(ba_ zscaled da_))<epsil.ang)
+ and (abs(angle(dc_ zscaled da_))<epsil.ang)
+ fi
+ fi
+enddef;
+% ---
+vardef is_tiny_bez(expr B) =
+% checks if B\'ezier segment |B| is negligibly small
+ (distance((postcontrol 0 of B),(point 0 of B))<epsil.dist)
+ and (distance((precontrol 1 of B),(point 0 of B))<epsil.dist)
+ and (distance((point 1 of B),(point 0 of B))<epsil.dist)
+enddef;
+% ---
+vardef are_parallel(expr B,C) =
+% checks if B\'ezier segments |B| and |C| are linear and parallel
+save a_;
+ if is_line(B) and is_line(C):
+ a_:=turn_ang((point 0 of B)-(point 1 of B),(point 0 of C)-(point 1 of C));
+ (if known a_: abs(a_)<epsil.ang else: false fi)
+ else: false fi
+enddef;
+% ---
+vardef tidy_lines(expr P) =
+% converts almost linear segments of a path |P| into a ``tidy'' lines (|--|)
+ save B_; path B_;
+ for i_:=1 upto length(P): if i_>1: & fi
+ hide(B_:=subpath (i_-1,i_) of P)
+ if is_line(B_): ((point 0 of B_)--(point 1 of B_)) else: B_ fi
+ endfor if cycle P: & cycle fi
+enddef;
+% ---
+def add_bez(expr ta,tb, p) =
+ .. controls (postcontrol ta of p) and (precontrol tb of p) .. (point tb of p)
+enddef;
+% ---
+vardef clean_path(expr P) =
+% this routine joins together colinear neighbouring segments and removes
+% ``tiny'' edges of a cyclic path |P| (performed at the end of removing
+% overlaps and expanding stroke); since some nodes may become ``midline''
+% ones after cleaning, the operation is performed twice
+ if cycle P:
+ save P_,for_del_,not_del_,i_,j_; path P_;
+% mark all deletable nodes and one non-deletable node:
+ for i_:=0 upto length(P)-1:
+ if are_parallel(subpath (i_-1,i_) of P,subpath (i_,i_+1) of P)
+ or is_tiny_bez(subpath(i_-1,i_) of P):
+ for_del_[i_]:=1;
+ else:
+ not_del_:=i_;
+ fi
+ endfor;
+% BUG TRAP:
+ if unknown not_del_:
+ err_helpless;
+ errmessage "ROEX ERROR: all nodes deleted during path cleaning";
+ fi
+% delete nodes:
+ i_:=j_:=not_del_; % we start with |not_del_|: one of not deleted points
+ P_:=(point j_ of P)
+ forever:
+% invariant: |i_| recent not deleted point, |j_| current point
+ hide(j_:=(j_+1) mod length(P))
+ if unknown for_del_[j_]: add_bez(i_,j_,P) \\ hide(i_:=j_) fi
+ exitif j_=not_del_;
+ endfor & cycle;
+ tidy_lines(P_)
+ else: P fi
+enddef;
+% ---
+vardef is_less(expr a,b) = (a<b) enddef;
+vardef quicksort@#(expr ii,jj)(suffix s)(text t) =
+% sorts |@#.s[ii..jj]| along with |@#.$[ii..jj]| for |$| in |t|,
+% using Tony Hoare's ``quick sort'' method; suffix |s| must must not occur
+% in the |t| list (no checking is performed); if both |s| and |t| are empty,
+% |t| is ignored.
+% REMARK 1: the algorithm has no explicit recursion, because of \MF's limits
+% on recursion level.
+% REMARK 2: the algorithm, of course, is not stable, i.e., it does not
+% preserve the order of equal items, but it does not matter here
+ save i_,j_,k_,l_,cell_,stack_,incl_t_; boolean incl_t_;
+ pair stack_[\\]; stack_.lev:=0; stack_[incr stack_.lev]:=(ii,jj);
+ i_:=0; for $:=t: i_:=i_+1; endfor % ``measure'' |t|-list
+ incl_t_:=(str s <> "") or ((str s = "") and (i_<>0));
+ forsuffixes $:= s if incl_t_: , t fi:
+ if numeric @#.$[ii]: numeric cell_.$;
+ elseif string @#.$[ii]: string cell_.$;
+ elseif boolean @#.$[ii]: boolean cell_.$;
+ fi
+ endfor
+ forever:
+ exitif stack_.lev<=0;
+ numeric i_,j_; (i_,j_)=stack_[stack_.lev]; stack_.lev:=stack_.lev-1;
+ if i_<j_:
+ forsuffixes $:= s if incl_t_: , t fi: cell_.$:=@#.$[i_]; endfor
+ l_:=i_;
+ for k_:=i_+1 upto j_:
+ if is_less(@#.s[k_],cell_.s):
+ forsuffixes $:=s if incl_t_: , t fi:
+ @#.$[l_]:=@#.$[k_]; @#.$[k_]:=@#.$[l_+1];
+ endfor
+ l_:=l_+1;
+ fi
+ endfor
+ forsuffixes $:= s if incl_t_: , t fi: @#.$[l_]:=cell_.$; endfor
+ stack_[incr stack_.lev]:=(i_,l_-1); stack_[incr stack_.lev]:=(l_+1,j_);
+ fi
+ endfor
+enddef;
+% ---
+% R-O MACROS:
+% ---
+% visualising macros (useful for testing):
+% ---
+def mark_nodes =
+ if proofing>0:
+ for i_:=1 upto NODE_.num:
+ makelabel(decimal(i_) & ":" & decimal(NODE_.pth[i_]),
+ point TIME_[NODE_.pth[i_]]tim[NODE_.tim[i_]] of PATH_[NODE_.pth[i_]]);
+ endfor
+ fi
+enddef;
+% ---
+def mark_area(expr i) =
+begingroup
+ save j_,v_; j_:=i; mark_edge(j_); v_[j_]:=0;
+ forever: j_:=EDGE_.out[j_]; exitif (j_=i) or (known v_.emerg);
+ if known v_[j_]: v_.emerg:=0; else: mark_edge(j_); v_[j_]:=0; fi
+ endfor
+endgroup
+enddef;
+% ---
+def mark_edge(expr i) =
+begingroup
+ if proofing>0:
+ save currentpen, currentpen_path; pen currentpen; path currentpen_path;
+ makelabel(decimal(i),
+ (point .5length(the_edge(i)) of the_edge(i))+
+ 1pt*(udir .5length(the_edge(i)) of the_edge(i)) rotated 90);
+ pickup pencircle scaled 1;
+ draw (point .5length(the_edge(i)) of the_edge(i))--
+ ((point .5length(the_edge(i)) of the_edge(i))+
+ (1pt*(udir .5length(the_edge(i)) of the_edge(i)) rotated 90));
+ makelabel("", point 0 of the_edge(i));
+ fi
+endgroup
+enddef;
+% ---
+def mark_edges =
+ for i_:=-EDGE_.num upto EDGE_.num: if i_<>0: mark_edge(i_); fi endfor
+enddef;
+% ---
+def show_area(expr i) =
+begingroup
+ save j_,v_; j_:=i;
+ message "EDGE " & decimal(j_) & "/" & decimal(EDGE_.pth[j_]) & ":";
+ message "color " &
+ if known EDGE_.col[j_]: decimal(EDGE_.col[j_]) else: "???" fi;
+ v_[j_]:=0;
+ forever: j_:=EDGE_.out[j_]; exitif (j_=i) or (known v_.emerg);
+ if known v_[j_]: v_.emerg:=0; fi
+ v_[j_]:=0; message " " & decimal(j_) & "/" & decimal(EDGE_.pth[j_]);
+ endfor
+endgroup
+enddef;
+% ---
+def show_areas =
+ for i_:=-EDGE_.num upto EDGE_.num: if i_<>0: show_area(i_); fi endfor
+enddef;
+% ---
+def err_helpless =
+ errhelp "Better stop now! Algorithm failed.";
+enddef;
+% ---
+def err_extra_info(expr i,j) =
+ message
+"========================== BEGIN OF ERROR INFO: ==========================";
+ for k_:=i,j:
+ if known k_:
+ message "Edge " & decimal(k_) &
+ " (a subpath of the path " & decimal(EDGE_.pth[k_]) & "):";
+ message "Color:"; show EDGE_.col[k_]; show the_edge(k_);
+ fi
+ endfor;
+enddef;
+% ---
+% principal macros:
+% ---
+vardef edge_path(expr i) = PATH_[EDGE_.pth[i]] enddef;
+vardef first_time(expr i) =
+ TIME_[NODE_.pth[EDGE_.fnd[i]]]tim[NODE_.tim[EDGE_.fnd[i]]]
+enddef;
+vardef last_time(expr i) =
+ TIME_[NODE_.pth[EDGE_.lnd[i]]]tim[NODE_.tim[EDGE_.lnd[i]]]
+enddef;
+% ---
+vardef the_edge(expr i) =
+ if i>0: pos_subpath else: neg_subpath fi
+ (first_time(i), last_time(i)) of edge_path(i)
+enddef;
+% ---
+vardef make_area(expr i) =
+ save j_,q_,v_; path q_; j_:=i; v_[j_]:=0; q_:=the_edge(j_);
+ forever: j_:=EDGE_.out[j_]; exitif (j_=i) or (known v_.emerg);
+ if known v_[j_]:
+ show_area(i); err_helpless;
+ errmessage "RO ERROR: Edge " & decimal(j_) & " revisited";
+ v_.emerg:=0;
+ fi
+ v_[j_]:=0; q_:=q_ && the_edge(j_);
+ endfor
+ make_cycle(q_)
+enddef;
+% ---
+vardef is_tangent(expr i,j,k,l) =
+ save e_,d_,pi_,pj_,ti_,tj_; path e_,pi_,pj_;
+ if (TIME_[i]num=0) or (TIME_[j]num=0): true
+ else:
+ ti_.loc:=TIME_[i]tim[k];
+ ti_.prv:=TIME_[i]tim[(k-1) mod (TIME_[i]num+1)];
+ ti_.nxt:=TIME_[i]tim[(k+1) mod (TIME_[i]num+1)];
+ tj_.loc:=TIME_[j]tim[l];
+ tj_.prv:=TIME_[j]tim[(l-1) mod (TIME_[j]num+1)];
+ tj_.nxt:=TIME_[j]tim[(l+1) mod (TIME_[j]num+1)];
+ pi_:=PATH_[i] shifted (-point ti_.loc of PATH_[i]);
+ pj_:=PATH_[j] shifted (-point tj_.loc of PATH_[j]);
+ d_:=min(
+ distance(point ti_.loc of pi_, point ti_.prv of pi_),
+ distance(point ti_.loc of pi_, point ti_.nxt of pi_),
+ distance(point tj_.loc of pj_, point tj_.prv of pj_),
+ distance(point tj_.loc of pj_, point tj_.nxt of pj_));
+% BUG TRAP 1:
+ if d_<epsil.dist:
+ err_helpless;
+ errmessage "RO ERROR: Cannot check tangency (too short secants)";
+ fi
+ e_:=enc.pth scaled (1/2[epsil.dist,d_]);
+ save ta_,tb_,tc_,td_;
+ save tt_;
+ (tt_,ta_)=(pos_subpath (ti_.prv,ti_.loc) of pi_) intersectiontimes e_;
+ save tt_;
+ (tt_,tb_)=(pos_subpath (ti_.loc,ti_.nxt) of pi_) intersectiontimes e_;
+ save tt_;
+ (tt_,tc_)=(pos_subpath (tj_.prv,tj_.loc) of pj_) intersectiontimes e_;
+ save tt_;
+ (tt_,td_)=(pos_subpath (tj_.loc,tj_.nxt) of pj_) intersectiontimes e_;
+% BUG TRAP 2:
+ if (ta_<0) or (tb_<0) or (tc_<0) or (td_<0):
+ err_helpless; errmessage "RO ERROR: Cannot check tangency";
+ fi
+ forsuffixes tt_:=tb_,tc_,td_: tt_:=(tt_-ta_) mod enc.len; endfor
+ ((tc_>=tb_) and (td_>=tb_)) or ((tc_<=tb_) and (td_<=tb_))
+ fi
+enddef;
+% ---
+vardef multi_path_case = PATH_.num>1 enddef;
+def prepare_input_data(text P)(text W) =
+% |P|: list of paths to be processed (non-cyclic paths are ignored);
+% |W|: list of weights given as pairs: (index, value)
+ PATH_.num:=0;
+ for P_:=P: if cycle P_: PATH_[incr PATH_.num]:=touch_path(P_); fi endfor
+ for W_:=W: PATH_.wei[xpart W_]:=ypart W_; endfor
+ for i_:=1 upto PATH_.num:
+ if unknown PATH_.wei[i_]: PATH_.wei[i_]:=1; fi
+ endfor
+enddef;
+% ---
+def initialise_removing_overlaps =
+% Given paths are |PATH_1|, |PATH_2|, ..., |PATH_[P.num]|;
+% if |PATH_[i][j]| is known, paths |PATH_[i]| and |PATH_[j]| at least touch
+% each other; |PATH_.wei[i]| is a weight of a path (corresponds to
+% multiplying a turning number by this value or, in other words, to
+% applying |PATH_.wei[i]| times a fill operation to the path |PATH_[i]|).
+numeric PATH_.num, PATH_[\\][\\], PATH_.wei[\\]; path PATH_[\\];
+%
+% Lone paths are stored in variable |LONE_|; |LONE_.col[i]| determines
+% the color (being an integer number) of the plane surrounding the path
+% |LONE_[i]|; |LONE_.wei| is a weight inherited from |PATH_.wei| (see above).
+numeric LONE_.num, LONE_.col[\\], LONE_.wei[\\]; path LONE_[\\];
+%
+% |TIME_[i]num| is the number of intersection points for paths |PATH_[i]|,
+% |TIME_[i]tim[j]| is the time of intersection of the |j|-th point of path
+% |PATH_[i]| (points are sorted with respect to time), |TIME_[i]ntp[j]|
+% marks non-tangent points (if known), |TIME_[i]nod[j]| is the node number
+% of |j|-th point of path |PATH_[i]| (only non-tangent points are considered
+% to be nodes, points on a path are numbered from |0|).
+numeric TIME_[\\]num, TIME_[\\]tim[\\], TIME_[\\]ntp[\\], TIME_[\\]nod[\\];
+%
+% Variables with prefix |EDGE_| describe the edge structure that results from
+% intersecting process; the data structure is similar to Dijkstra's data
+% structure for the algorithm finding the convex hull of for a given
+% set of points (E. W. Dijkstra, ``A Discipline of programming'',
+% Prentice-Hall, Inc., 1976): the edges (edsges?) are numbered
+% |-EDGE_.num|, |-EDGE_.num+1|, ..., |-1|, |1|, ... |EDGE_.num-1|,
+% |EDGE_.num|; edges |i| and |-i| are in fact the same edge but
+% differently oriented, the positive value denotes the edge which
+% direction is consistent with the direction of the original path;
+% |EDGE_.out[i]| is the number of the leftmost edge outcoming from
+% the last node, i.e., |EDGE_.lnd[i]|; the number of the first node
+% of |i|-th edge is |EDGE_.fnd[i]|; color of |i|-th edge, i.e., the color
+% of the area surrounded by the edge and its leftmost successors is stored
+% in |EDGE_.col[i]|; |i|-th edge belongs to the path of |PATH_[EDGE_.pth[i]]|;
+% |EDGE_.aux[\\]| is an auxiliary variable; all intersecting paths
+% can be grouped into |SPOT_.num| of disjoint ``spots''; for |i=1|, |2|, ...,
+% |SPOT_.num|, |EDGE_.bed[i]| is the number of the edge which leftmost
+% successors form the area being a boundary of the intersecting paths
+% for a given spot and |EDGE_.bpa[i]| is the boundary (since there is
+% one-to-one correspondence between boundaries and spots, boundaries are
+% pairwise disjoint, too).
+numeric EDGE_.num, EDGE_.pth[\\], EDGE_.out[\\], EDGE_.aux[\\], EDGE_.col[\\],
+ EDGE_.fnd[\\], EDGE_.lnd[\\], EDGE_.are[\\], EDGE_.bed[\\];
+path EDGE_.bpa[\\];
+%
+% Variables with prefix |NODE_| describe the node structure and are related
+% to the edge structure; the node is not a point on a plane but a point on a
+% path, hence several nodes may correspond to one Euclidian point;
+% |NODE_.num| is the number of nodes, |NODE_.pth[i]| is the number of a path
+% to which the node |i| belongs, |NODE_.tim[i]| is the corresponding time on
+% path |PATH_[NODE_.pth[i]]|, |NODE_.ped[i]| is the ordering number of a
+% positively-numbered edge leaving the node |i|, |NODE_.ned[i]| is the
+% ordering number of negatitively-numbered edge leaving node |i|,
+% |NODE_.nod[i]num| is the number of nodes coinciding with node |i| and
+% these are nodes |NODE_.nod[i]1|, |NODE_.nod[i]2|, ...,
+% |NODE_.nod[i][NODE_.nod[i]num]|.
+numeric NODE_.num, NODE_.pth[\\], NODE_.tim[\\], NODE_.ned[\\],
+ NODE_.ped[\\], NODE_.nod[\\]num, NODE_.nod[\\][\\];
+%
+% It often happens that intersecting paths form disjoint areas or that
+% there are paths that do not intersect; it is crucial for the colouring
+% algorithm to know the embedding ``hierarchy''; the hierarchy is stored in
+% a tree structure: the links (suffix |emb|) point upward (from leaves to
+% the root), moreover, with each leaf (node) is associated a ``level,''
+% i.e., the number of leaves beneath this leaf; the information stored in
+% a leaf is either the number of a lone path or the number of an area being
+% the result of the intersecting process (negative value marks the former
+% case); the tree is built by adding at first lone paths and next boundary
+% paths, i.e., negatively oriented paths surrounding groups of paths (areas,
+% see below) resulting from the intersecting process; if there is a lone path
+% or a boundary path |q| embedded in a boundary path |p|, there must be also
+% an area which belongs to the group of areas surrounded by |p|, which apears
+% in the tree between |q| and |p|; such a structure is convenient at the
+% stage of finding colors of areas and lone paths.
+numeric TREE_.num,TREE_.pth[\\],TREE_.emb[\\],TREE_.lev[\\];
+%
+% Finally there are areas which arise during intersecting process;
+% |AREA_1|, |AREA_2|, ..., |AREA_[AREA_.num]| are the ordering numbers
+% of edges which leftmost successors form areas such that areas arising
+% from |AREA_[i]| and |AREA_[j]| are either disjoint or have at most a common
+% edge, and, moreover, areas arising from |AREA_1|, |AREA_2|, ...,
+% |AREA_[AREA_.num]| exhaust the list of all possible areas in question
+% (plus lone paths, i.e., with no intersecting points); |AREA_.spt[i]| is
+% a spot number (areas of the same spot number are subsets of the same
+% boundary, different boundaries are disjoint); areas are sorted wrt spot
+% numbers, moreover, |AREA_[SPOT_[s-1]+1]| thru |AREA_[SPOT_[s]]| are areas
+% belonging to a spot |s|, |s=1|, |2|, ..., |SPOT_.num|.
+numeric AREA_.num, AREA_[\\], AREA_.spt[\\]; numeric SPOT_.num, SPOT_[\\];
+enddef;
+% ---
+vardef is_far_enough(expr i,k,dk) =
+ if act_idx_>=max_idx: were_more_:=1; false
+ else:
+ save z_; pair z_; z_:=point k+dk of PATH_[i];
+ true
+ for j_:=0 upto TIME_[i]num:
+ and
+ if multi_path_case:
+ (distance(point TIME_[i]tim[j_] of PATH_[i],z_)>=epsil.dist)
+ else:
+ (interval(TIME_[i]tim[j_],(k+dk),PATH_[i])>=epsil.time)
+ fi
+ endfor
+ for j_:=0 upto ignored_.num:
+ and
+ if multi_path_case:
+ (distance(point ignored_[j_] of PATH_[i],z_)>=epsil.dist)
+ else:
+ (interval(ignored_[j_],(k+dk),PATH_[i])>=epsil.time)
+ fi
+ endfor
+ fi
+enddef;
+% ---
+def intersect_two_segments(expr i,j,k,l) =
+begingroup
+ save pi_,pj_,stack_; path pi_,pj_,stack_[\\]; numeric stack_.lev;
+ pi_:=subpath (k,k+1) of PATH_[i]; pj_:=subpath (l,l+1) of PATH_[j];
+ stack_.lev:=1; stack_[stack_.lev]:=pj_;
+ forever:
+ exitif stack_.lev<=0;
+ pj_:=stack_[stack_.lev]; stack_.lev:=stack_.lev-1;
+ save dk_,dl_; (dk_,dl_)=pi_ intersectiontimes pj_;
+ if dk_>=0:
+ if is_far_enough(i,k,dk_):
+ act_idx_:=act_idx_+1; TIME_[i].tim[incr TIME_[i].num]:=k+dk_;
+ else: ignored_[incr ignored_.num]:=k+dk_;
+ fi
+ if (dl_+epsil.time)<length(pj_):
+ stack_[incr stack_.lev]:=subpath (dl_+epsil.time,length pj_) of pj_;
+ fi
+ if (dl_-epsil.time)>0:
+ stack_[incr stack_.lev]:=subpath (0,dl_-epsil.time) of pj_;
+ fi
+ fi
+ endfor
+endgroup
+enddef;
+% ---
+def intersect_two_paths(expr i,j) =
+begingroup
+ save ignored_,were_more_,act_idx_;
+ act_idx_:=0; ignored_.num:=-1;
+ for k_:=0 upto length(PATH_[i])-1:
+ for l_:=0 upto length(PATH_[j])-1:
+ if (i<>j) or (k_<>l_):
+ intersect_two_segments(i,j,k_,l_);
+ fi
+ endfor
+ endfor
+ if known were_more_:
+ errhelp "Dangerous situation: rounding errors may screw up results.";
+ errmessage "RO ERROR: there were more than "
+ & decimal(max_idx) & " intersections (thus some were ignored)";
+ fi
+ quicksort TIME_[i](0,TIME_[i].num)(tim)();
+endgroup
+enddef;
+% ---
+def intersect_all_paths =
+ for i_:=1 upto PATH_.num: TIME_[i_]num:=-1; endfor
+ for i_:=1 upto PATH_.num:
+ for j_:=i_+1 upto PATH_.num:
+ if xpart(PATH_[i_] intersectiontimes PATH_[j_])>-1:
+ PATH_[i_][j_]:=0;
+% the process is repeated twice (for both paths in turn) because we haven't
+% invented an efficient soultion to the following problem:
+% given a subpath |S| of a path |P| and a time |t.S|; find a time |t.P| such
+% that |point t.S of S=point t.P of P|
+ intersect_two_paths(i_,j_); intersect_two_paths(j_,i_);
+ fi
+ endfor
+ endfor
+ if not multi_path_case:
+ for i_:=1 upto PATH_.num:
+ PATH_[i_][i_]:=0; intersect_two_paths(i_,i_);
+ endfor
+ fi
+enddef;
+% ---
+vardef find_minimal_secant =
+ save secants_, intervals_;
+ secants_.num:=0; intervals_.num:=0;
+ minimal_secant:=minimal_interval:=infinity;
+ for i_:=1 upto PATH_.num:
+ for j_:=0 upto TIME_[i_].num:
+ if TIME_[i_].num>0:
+ secants_[if tracingremoving>1: incr fi \\ secants_.num]:=
+ distance(point TIME_[i_]tim[j_] of PATH_[i_],
+ point TIME_[i_]tim[(j_+1) mod (TIME_[i_].num+1)] of PATH_[i_]);
+ if tracingremoving>1:
+ secants_.pth[secants_.num]:=i_; secants_.tim[secants_.num]:=j_;
+ fi
+ minimal_secant:=min(minimal_secant,secants_[secants_.num]);
+ intervals_[if tracingremoving>1: incr fi \\ intervals_.num]:=
+ interval(TIME_[i_]tim[j_],
+ TIME_[i_]tim[(j_+1) mod (TIME_[i_].num+1)], PATH_[i_]);
+ if tracingremoving>1:
+ intervals_.pth[intervals_.num]:=i_; intervals_.tim[intervals_.num]:=j_;
+ fi
+ minimal_interval:=min(minimal_interval,intervals_[intervals_.num]);
+ fi
+ endfor;
+ endfor;
+ if minimal_secant<>infinity:
+ info_ro "Minimal secant = " & decimal(minimal_secant/pt)
+ & "pt, i.e., " & decimal(minimal_secant) & "pxl, " &
+ if minimal_secant<4/3epsil.dist: "CAVEAT!" else: "seems OK" fi
+ & " (bound=" & decimal(epsil.dist) & "pxl)";
+ fi
+ if minimal_interval<>infinity:
+ info_ro "Minimal interval = " & decimal(minimal_interval) & ", " &
+ if minimal_interval<4/3epsil.time: "CAVEAT!" else: "seems OK" fi
+ & " (bound=" & decimal(epsil.time) & ")";
+ fi
+ if tracingremoving>1:
+ quicksort secants_(1,secants_.num)()(tim,pth);
+ quicksort intervals_(1,intervals_.num)()(tim,pth);
+ for i_:=1 upto secants_.num:
+ info_ro "secant=" & decimal(secants_[i_])
+ & " path=" & decimal(secants_.pth[i_])
+ & " time=" & decimal(secants_.tim[i_]);
+ endfor
+ for i_:=1 upto intervals_.num:
+ info_ro "interval=" & decimal(intervals_[i_])
+ & " path=" & decimal(intervals_.pth[i_])
+ & " time=" & decimal(intervals_.tim[i_]);
+ endfor
+ fi
+enddef;
+% ---
+def build_node_structure =
+begingroup
+ save n_,Tik_,Tjl_;
+ NODE_.num:=0;
+ for i_:=1 upto PATH_.num:
+ for j_:=i_ if multi_path_case: +1 fi upto PATH_.num:
+ if known PATH_[i_][j_]:
+ for k_:=0 upto TIME_[i_]num:
+ for l_:=if i_=j_: k_ else: 0 fi upto TIME_[j_]num:
+ if distance(point TIME_[i_]tim[k_] of PATH_[i_],
+ point TIME_[j_]tim[l_] of PATH_[j_])<epsil.dist:
+ if if multi_path_case: not is_tangent(i_,j_,k_,l_) else: true fi:
+ TIME_[i_]ntp[k_]:=1; TIME_[j_]ntp[l_]:=1;
+ if unknown TIME_[i_]nod[k_]:
+ NODE_.num:=NODE_.num+1; TIME_[i_]nod[k_]:=NODE_.num;
+ NODE_.pth[NODE_.num]:=i_; NODE_.tim[NODE_.num]:=k_;
+ NODE_.nod[NODE_.num]num:=0;
+ fi
+ if unknown TIME_[j_]nod[l_]:
+ NODE_.num:=NODE_.num+1; TIME_[j_]nod[l_]:=NODE_.num;
+ NODE_.pth[NODE_.num]:=j_; NODE_.tim[NODE_.num]:=l_;
+ NODE_.nod[NODE_.num]num:=0;
+ fi
+ Tik_:=TIME_[i_]nod[k_]; Tjl_:=TIME_[j_]nod[l_];
+ NODE_.nod[Tik_][incr NODE_.nod[Tik_]num]:=Tjl_;
+ if (i_<>j_) or (k_<>l_):
+ NODE_.nod[Tjl_][incr NODE_.nod[Tjl_]num]:=Tik_;
+ fi
+ fi
+ fi
+ endfor
+ endfor
+ fi
+ endfor
+ endfor
+% BUG TRAP:
+ for i_:=1 upto PATH_.num:
+ n_:=0;
+ for j_:=0 upto TIME_[i_]num: if known TIME_[i_]ntp[j_]: n_:=n_+1; fi endfor
+ if n_=1:
+ err_helpless;
+ errmessage "RO ERROR: Number of non-tangent points must not be 1 (path "
+ & decimal(i_) & ")";
+ fi
+ endfor
+endgroup
+enddef;
+% ---
+def identify_close_nodes =
+begingroup
+% It is assumed that `close points' and `coinciding points' means the same,
+% hence we make a transitive closure of the relation of `being close'
+% (in ``normal'' cases the relation is transitive, although from
+% a mathematical point of view it is obviously not):
+ forever:
+ % Is this loop really needed? I (BJ) could not devise a case where
+ % more than two turns would be necessary and where the algorithm
+ % would still work properly
+ save changed_;
+ for i_:=1 upto NODE_.num:
+ save v_;
+ for j_:=1 upto NODE_.nod[i_]num: v_[NODE_.nod[i_][j_]]:=0; endfor
+ for j_:=1 upto NODE_.nod[i_]num:
+ k_:=NODE_.nod[i_][j_];
+ for l_:=1 upto NODE_.nod[k_]num:
+ if NODE_.nod[k_][l_]<>i_:
+ if unknown v_[NODE_.nod[k_][l_]]:
+ NODE_.nod[i_][incr NODE_.nod[i_]num]:=NODE_.nod[k_][l_];
+ changed_:=v_[NODE_.nod[k_][l_]]:=0;
+ fi
+ fi
+ endfor
+ endfor
+ endfor
+ exitif unknown changed_;
+ endfor
+endgroup
+enddef;
+% ---
+def build_edge_structure =
+begingroup
+ numeric min_sec_[\\],min_sec_.tmp;
+% |min_sec_[i]| is a minimal secant for |i|-th path, |i=1|, |2|, ..., |PATH_.num|
+ save i_,j_,k_;
+ EDGE_.num:=0;
+ if NODE_.num>0:
+ for i_:=1 upto PATH_.num:
+ for j_:=0 upto TIME_[i_]num:
+ if known TIME_[i_]ntp[j_]:
+ EDGE_.num:=EDGE_.num+1; EDGE_.pth[EDGE_.num]=i_;
+ EDGE_.pth[-EDGE_.num]=i_; EDGE_.fnd[EDGE_.num]:=TIME_[i_]nod[j_];
+% each path should contain at least two non-tangent nodes
+ k_:=j_+1;
+ forever: exitif known TIME_[i_]ntp[k_ mod (TIME_[i_]num+1)];
+ k_:=k_+1;
+ endfor;
+ EDGE_.lnd[EDGE_.num]:=TIME_[i_]nod[k_ mod (TIME_[i_]num+1)];
+ min_sec_.tmp:=
+ distance(point first_time(EDGE_.num) of edge_path(EDGE_.num),
+ point last_time(EDGE_.num) of edge_path(EDGE_.num));
+ if if unknown min_sec_[i_]: true else: min_sec_[i_]>min_sec_.tmp fi:
+ min_sec_[i_]:=min_sec_.tmp;
+ fi
+ EDGE_.fnd[-EDGE_.num]:=EDGE_.lnd[EDGE_.num];
+ EDGE_.lnd[-EDGE_.num]:=EDGE_.fnd[EDGE_.num];
+ fi
+ endfor
+ endfor
+ for i_:=-EDGE_.num upto EDGE_.num:
+ if i_>0: NODE_.ped[EDGE_.fnd[i_]]:=i_;
+ elseif i_<0: NODE_.ned[EDGE_.fnd[i_]]:=i_;
+ fi
+ endfor
+ else:
+ info_ro "RO WARNING: no intersections detected.";
+ fi
+endgroup
+enddef;
+% ---
+def find_leftmost_edges =
+% a simple method is used: a tiny circle is drawn in a node and its
+% intersection points with all edges leaving the node are examined
+begingroup
+ save ei_,ej_,i_,j_,k_,leftmost_; path ei_,ej_; numeric min_sec_.loc;
+ for i_:=-EDGE_.num upto EDGE_.num: if i_<>0:
+ if tracingleftmost>0:
+ message "@@@ " & decimal(i_) & "/" & decimal(EDGE_.pth[i_])
+ & " (" & decimal(EDGE_.fnd[i_]) & "," & decimal(EDGE_.lnd[i_]) & "):";
+ fi
+ numeric leftmost_.edg,leftmost_.tim,leftmost_.tmp;
+ min_sec_.loc:=min_sec_[EDGE_.pth[i_]];
+ for k_:=1 upto NODE_.nod[EDGE_.lnd[i_]]num:
+ forsuffixes $:=ped,ned:
+ j_:=NODE_$[NODE_.nod[EDGE_.lnd[i_]][k_]];
+ min_sec_.loc:=min(min_sec_.loc,min_sec_[EDGE_.pth[j_]]);
+ endfor
+ endfor
+% BUG TRAP 1: (should not happen, see |find_minimal_secant|);
+ if min_sec_.loc<epsil.dist:
+ err_extra_info(i_,whatever); showvariable min_sec_; err_helpless;
+ errmessage "RO ERROR: Cannot continue searching for the leftmost edge";
+ fi
+ ei_:=enc.pth scaled (1/2[epsil.dist,min_sec_.loc])
+ shifted point infinity of the_edge(i_);
+ save tei_,tt_; (tei_,tt_)=ei_ intersectiontimes the_edge(i_);
+ for k_:=1 upto NODE_.nod[EDGE_.lnd[i_]]num:
+ forsuffixes $:=ped,ned:
+ j_:=NODE_$[NODE_.nod[EDGE_.lnd[i_]][k_]];
+ ej_:=enc.pth scaled (1/2[epsil.dist,min_sec_.loc])
+ shifted point 0 of the_edge(j_);
+ if tracingleftmost>0:
+ message " " & decimal(j_) & "/" & decimal(EDGE_.pth[j_])
+ & " (" & decimal(EDGE_.fnd[j_]) & "," & decimal(EDGE_.lnd[j_]) & "):";
+ fi
+ save tej_,tt_; (tej_,tt_)=ej_ intersectiontimes the_edge(j_);
+% BUG TRAP 2:
+ if (tei_<0) or (tej_<0):
+ err_extra_info(i_,j_); showvariable min_sec_;
+ message "Times: " & decimal(tei_) & " " & decimal(tej_);
+ err_helpless;
+ errmessage "RO ERROR: Unsuccesful search for the leftmost edge";
+ fi
+% it happens that |i_=j_| if |multi_path_case=false|
+ leftmost_.tmp:=if (i_=-j_): 0 else: (tej_-tei_) mod enc.len fi;
+ if tracingleftmost>0:
+ message " " & decimal(leftmost_.tmp) & " " & dec_pair((tei_,tej_));
+ fi
+ if if unknown leftmost_.tim: true else: leftmost_.tmp>leftmost_.tim fi:
+ leftmost_.edg:=j_; leftmost_.tim:=leftmost_.tmp;
+ fi
+ endfor
+ endfor
+ EDGE_.out[i_]:=leftmost_.edg;
+ if tracingleftmost>0: j_:=EDGE_.out[i_];
+ message ">>> " & decimal(j_) & "/" & decimal(EDGE_.pth[j_])
+ & " (" & decimal(EDGE_.fnd[j_]) & "," & decimal(EDGE_.lnd[j_]) & "):";
+ fi
+ fi endfor
+endgroup
+enddef;
+% ---
+def build_area_structure =
+begingroup
+ save i_,j_,v_;
+ AREA_.num:=0;
+ for i_:=1 upto EDGE_.num: % \MF's linear equation solver employed
+ EDGE_.col[i_]-PATH_.wei[EDGE_.pth[i_]]=EDGE_.col[-i_];
+ endfor
+% split to areas (edges surrounding the same area are assigned the same colour):
+ for i_:=-EDGE_.num upto EDGE_.num:
+ if (i_<>0) and (unknown EDGE_.are[i_]):
+ AREA_[incr AREA_.num]:=i_;
+ save v_; j_:=i_; v_[j_]:=0; EDGE_.are[j_]:=AREA_.num;
+ forever: j_:=EDGE_.out[j_]; exitif (j_=i_) or (known v_.emerg);
+% BUG TRAP 1:
+ if j_=-i_:
+ err_extra_info(i_,whatever); show_area(i_); err_helpless;
+ errmessage "RO ERROR: strange area";
+ fi
+% BUG TRAP 2:
+ if known v_[j_]:
+ err_extra_info(i_,j_); show_area(i_); err_helpless;
+ errmessage "RO ERROR: Edge " & decimal(j_) & " revisited";
+ v_.emerg:=0;
+ fi
+ v_[j_]:=0; EDGE_.are[j_]:=AREA_.num;
+ if known (EDGE_.col[i_]-EDGE_.col[j_]):
+% BUG TRAP 3:
+ if (EDGE_.col[i_]-EDGE_.col[j_])<>0:
+ err_extra_info(i_,j_); show_area(i_); err_helpless;
+ errmessage "RO ERROR: Edges " & decimal(i_) & " and " & decimal(j_)
+ & " have inconsistent colors";
+ fi
+ else: EDGE_.col[i_]=EDGE_.col[j_]; fi
+ endfor
+ fi
+ endfor
+endgroup
+enddef;
+% ---
+def build_spot_structure =
+begingroup
+ save i_,j_;
+% areas having a common edge belong to the same spot (\MF's linear
+% equation solver employed):
+ for i_:=1 upto EDGE_.num: if (i_<>0):
+ if unknown (AREA_.spt[EDGE_.are[i_]]-AREA_.spt[EDGE_.are[-i_]]):
+ AREA_.spt[EDGE_.are[i_]]=AREA_.spt[EDGE_.are[-i_]];
+ fi
+ fi endfor
+% count different spots:
+ SPOT_.num:=0;
+ for i_:=1 upto AREA_.num:
+ if unknown AREA_.spt[i_]: AREA_.spt[i_]=incr SPOT_.num; fi
+ endfor;
+% sort areas wrt spot numbers:
+ quicksort AREA_(1,AREA_.num)(spt)();
+% define |SPOT_[s]|, |s=1|, |2|, ..., |SPOT_.num|, such that
+% |AREA_[SPOT_[s-1]+1]| thru |AREA_[SPOT_[s]]| are areas having the
+% same spot number |s|:
+ SPOT_0=0; for i_:=1 upto AREA_.num: SPOT_[AREA_.spt[i_]]:=i_; endfor
+% identify paths for which |check_turn=-1| (such areas are boundaries and
+% should be unique for each spot):
+ i_:=0;
+ for j_:=1 upto AREA_.num:
+ if check_turn(make_area(AREA_[j_]))<=0:
+ i_:=i_+1;
+ EDGE_.bar[i_]:=j_; EDGE_.bed[i_]:=AREA_[j_];
+ EDGE_.bpa[i_]:=make_area(AREA_[j_]);
+ fi
+ endfor
+% BUG TRAP 4:
+ if i_<>SPOT_.num:
+ message "Number of spots=" & decimal(SPOT_.num) &
+ ", number of boundaries=" & decimal(i_);
+ err_helpless;
+ errmessage "RO ERROR: Inconsistent number of spots and boundaries";
+ fi
+endgroup
+enddef;
+% ---
+def update_tree_levels(expr n) =
+% update levels above the inserted leaf (|n|):
+begingroup
+ save i_,j_;
+ i_:=TREE_.emb[n]; j_:=TREE_.lev[n];
+ forever:
+ exitif i_=0;
+ j_:=TREE_.lev[i_]:=max(j_+1,TREE_.lev[i_]); i_:=TREE_.emb[i_];
+ endfor
+endgroup
+enddef;
+% ---
+vardef embedding_pair(expr p,n) =
+% returns a pair of numbers, |(out_,in_)|, such that |out_| is either a tree
+% address of an area surrounded |p| or zero if not found, and |in_| is a
+% tree address of an area surrounding |p| or zero if not found; one branch
+% is searched, starting from a zero level leaf, |n|.
+ save in_,out_,q_; path q_; in_:=out_:=0; n_:=n;
+ forever:
+ q_:=if TREE_.pth[n_]<0: LONE_[-TREE_.pth[n_]]
+ else: make_area(AREA_[TREE_.pth[n_]]) fi;
+ check_embedding(p_,q_,r_);
+ if r_=1: % |p_| $\subset$ |q_|
+ if in_=0: in_:=n_; fi % ``minimal'' surrounding path is to be found
+ elseif r_=2: % |q_| $\subset$ |p_|
+ out_:=n_; % ``maximal'' surrounded path is to be found
+ fi
+ exitif TREE_.emb[n_]=0;
+ n_:=TREE_.emb[n_];
+ endfor
+ (out_,in_)
+enddef;
+% ---
+def add_to_queue (expr p,q) =
+% updates a queue, i.e., adds to a queue an area belonging to the spot
+% of a boundary area, i.e., |AREA_[q]|, surrounding a path described by
+% a tree address |p|.
+begingroup
+ save p_,r_,s_,found_; path p_; boolean found_;
+% |AREA_[q]| is a boundary, i.e., it is negatively oriented,
+ p_:=if TREE_.pth[p]<0: LONE_[-TREE_.pth[p]]
+ else: make_area(AREA_[TREE_.pth[p]]) fi;
+% there must be a unique positively oriented area containing |p_|, belonging
+% to the ``spot'' of the boundary |AREA_[q]|.
+ found_:=false; s_:=SPOT_[AREA_.spt[q]-1];
+% |s_| is the last area from a previous spot, |s_+1| will be the first
+% area of the current spot
+ forever: s_:=s_+1;
+ if s_<>q:
+ check_embedding(make_area(AREA_[s_]),p_,r_);
+ found_:=(r_=2); % |r_=2| implies |p_| $\subset$ |make_area(AREA_[s_])|
+ fi
+ if found_:
+ if unknown LVQ_.inq[s_]: LVQ_[incr LVQ_.num]:=s_; LVQ_.inq[s_]:=1; fi
+ elseif s_=SPOT_[AREA_.spt[q]]:
+% the list of candidates has been exhausted without a success, hence
+% BUG TRAP:
+ err_helpless;
+ errmessage "RO ERROR: cannot build embedding tree (boundary " &
+ decimal(q) & ")";
+ found_:=true;
+ fi
+ exitif found_;
+ endfor;
+endgroup
+enddef;
+% ---
+def add_to_tree (expr leaf) =
+begingroup
+ save boundary_,found_,N_,p_; path p_;
+ TREE_.num:=TREE_.num+1;
+ N_:=TREE_.num; % abbreviation
+ TREE_.pth[N_]:=leaf; TREE_.emb[N_]:=0; TREE_.lev[N_]:=0;
+ p_:=if leaf<0: LONE_[-leaf] else: make_area(AREA_[leaf]) fi;
+ if (leaf>0) and (check_turn(p_)<0): boundary_:=1; fi
+ if N_>1:
+ for l_:=1 upto LVZ_.num: % climbing up from level zero
+ save out_,in_; (out_,in_)=embedding_pair(p_,LVZ_[l_]);
+ if (out_<>0) or (in_<>0): % a feasible branch found
+ if (out_=0) and (in_<>0): % to be added at the bottom, certain
+ TREE_.emb[N_]:=in_;
+ if in_=LVZ_[l_]: LVZ_[l_]:=N_; % replace bottom leaf
+ else: LVZ_[incr LVZ_.num]:=N_; fi % add new bottom leaf
+ elseif (out_<>0) and (in_=0): % to be added at the top, optional
+ if TREE_.emb[out_]<>N_: % we weren't here, add
+% invariant: |TREE_.emb[out_]=0|
+ TREE_.emb[out_]:=N_;
+ TREE_.lev[N_]:=max(TREE_.lev[N_],TREE_.lev[out_]+1);
+ if known boundary_: add_to_queue(out_,leaf); fi
+ fi
+ else: % to be added in the midst, optional
+ if TREE_.emb[out_]<>N_: % we weren't here, add
+% invariant: |TREE_.emb[out_]=in_|
+ TREE_.emb[out_]:=N_; TREE_.emb[N_]:=in_;
+ TREE_.lev[N_]:=max(TREE_.lev[N_],TREE_.lev[out_]+1);
+ if known boundary_: add_to_queue(out_,leaf); fi
+ fi
+ fi
+ found_:=1; update_tree_levels(N_);
+ fi
+ endfor;
+ fi
+ if unknown found_: LVZ_[incr LVZ_.num]:=N_; fi % a ``separate'' leaf appeared
+endgroup
+enddef;
+% ---
+def build_embedding_tree =
+begingroup
+ save LVZ_,QUE_;
+% |LVZ_1|, |LVZ_2|, ..., |LVZ_[LVZ_.num]| is the list of zero-level leaves,
+% (a temporary data structure, used during building a tree), |LVQ_1|, |LVQ_2|,
+% ..., |LVQ_[LVQ_.num]| is the list of leaves waiting in a queue (also
+% a temporary data structure); if |LVQ_.inq[i]| is known, |i=1|, |2|, ...,
+% |AREA_.num|, area |i| is already in a queue.
+ TREE_.num:=0; LONE_.num:=0; LVZ_.num:=0; LVQ_.num:=0;
+% identify lone paths:
+ for i_:=1 upto PATH_.num:
+ if true for j_:=0 upto TIME_[i_]num: and (unknown TIME_[i_]ntp[j_]) endfor:
+ LONE_[incr LONE_.num]:=PATH_[i_]; LONE_.wei[LONE_.num]:=PATH_.wei[i_];
+ fi
+ endfor
+% build the tree:
+ for i_:=1 upto LONE_.num: add_to_tree(-i_); endfor
+ for i_:=1 upto SPOT_.num: add_to_tree(EDGE_.bar[i_]); endfor
+ for i_:=1 upto LVQ_.num: add_to_tree(LVQ_[i_]); endfor
+endgroup
+enddef;
+% ---
+def color_paths = % \MF's linear equation solver heavily exploited
+begingroup
+ save i_,j_;
+ for i_:=1 upto TREE_.num:
+ if TREE_.emb[i_]=0: % outer path
+ if TREE_.pth[i_]<0: LONE_.col[-TREE_.pth[i_]]=background_color;
+ else: EDGE_.col[AREA_[TREE_.pth[i_]]]=background_color; fi
+ else: % inner path, inherits color from the surrounding path
+ j_:=TREE_.emb[i_];
+ if TREE_.pth[i_]<0:
+ if TREE_.pth[j_]<0:
+ LONE_.col[-TREE_.pth[i_]]=LONE_.col[-TREE_.pth[j_]]
+ +LONE_.wei[-TREE_.pth[j_]]*check_turn(LONE_[-TREE_.pth[j_]]);
+ else:
+ LONE_.col[-TREE_.pth[i_]]=EDGE_.col[AREA_[TREE_.pth[j_]]];
+ fi
+ else:
+ if TREE_.pth[j_]<0:
+ EDGE_.col[AREA_[TREE_.pth[i_]]]=LONE_.col[-TREE_.pth[j_]]
+ +LONE_.wei[-TREE_.pth[j_]]*check_turn(LONE_[-TREE_.pth[j_]]);
+ else:
+ if AREA_.spt[TREE_.pth[i_]]<>AREA_.spt[TREE_.pth[j_]]:
+ EDGE_.col[AREA_[TREE_.pth[i_]]]=EDGE_.col[AREA_[TREE_.pth[j_]]];
+ fi
+ fi
+ fi
+ fi
+ endfor;
+endgroup
+enddef;
+% ---
+def recombine_edges(suffix R) =
+% this routine can be used several times (after completing the process
+% of finding the structure of paths after intersecting) with various
+% definitions of |good_color| function in order to select various
+% sets of areas
+ if not path R0: numeric R.num; path R[\\]; fi
+ if (unknown R.num) or (unknown append_results): R.num:=0; fi
+% |R|: resulting data structure, namely, |R.num| is the number of output
+% paths, |R1|, |R2|, ..., |R[R.num]| are the resulting paths
+begingroup
+ save i_,j_,out_,in_;
+ for i_:=1 upto LONE_.num:
+ out_:=LONE_.col[i_]; in_:=out_+LONE_.wei[i_]*check_turn(LONE_[i_]);
+ if good_colors(in_,out_) or good_colors(out_,in_):
+ R[incr R.num]:=LONE_[i_];
+ R[R.num]:=if good_colors(in_,out_): pos_turn else: neg_turn fi \\ R[R.num];
+ fi
+ endfor
+ for i_:=-EDGE_.num upto EDGE_.num: if i_<>0:
+ EDGE_.aux[i_]:=whatever;
+ fi endfor
+ for i_:=-EDGE_.num upto EDGE_.num: if i_<>0:
+% BUG TRAP 1:
+ if unknown EDGE_.col[i_]:
+ err_extra_info(i_,whatever); err_helpless;
+ errmessage "RO ERROR: Edge " & decimal(j_) & " not colored";
+ fi
+ if good_colors(EDGE_.col[i_],EDGE_.col[-i_]) and (unknown EDGE_.aux[i_]):
+ save v_;
+ R.num:=R.num+1; j_:=i_; v_[j_]:=0;
+ EDGE_.aux[j_]:=0; R[R.num]:=the_edge(j_);
+ forever: j_:=EDGE_.out[j_]; exitif (j_=i_) or (known v_.emerg);
+% BUG TRAP 2:
+ if known v_[j_]:
+ err_extra_info(i_,j_); err_helpless;
+ errmessage "RO ERROR: Edge " & decimal(j_) & " revisited";
+ v_.emerg:=0;
+ fi
+ v_[j_]:=0;
+ if good_colors(EDGE_.col[j_],EDGE_.col[-j_]):
+ EDGE_.aux[j_]:=0; R[R.num]:=R[R.num] && the_edge(j_);
+ else: j_:=-j_;
+ fi
+ endfor
+ R[R.num]:=clean_path(clean_path(make_cycle(R[R.num])));
+ fi
+ fi endfor
+endgroup
+enddef;
+% ---
+def remove_overlap (text P)(text W) suffix R =
+begingroup interim autorounding:=0;
+% |P|: list of paths to be processed (non-cyclic paths are ignored);
+% |W|: list of weights given as pairs: (index, value)
+% |R|: resulting data structure, i.e., |R.num| is the number of output paths,
+% |R1|, |R2|, ..., |R[R.num]| are the resulting paths
+ info_ro "initialise_removing_overlaps"; initialise_removing_overlaps;
+ info_ro "prepare_input_data"; prepare_input_data(P)(W);
+ info_ro "intersect_all_paths"; intersect_all_paths;
+ info_ro "find_minimal_secant"; find_minimal_secant;
+ info_ro "build_node_structure"; build_node_structure;
+ info_ro "identify_close_nodes"; identify_close_nodes;
+ info_ro "build_edge_structure"; build_edge_structure;
+ info_ro "find_leftmost_edges"; find_leftmost_edges;
+ info_ro "build_area_structure"; build_area_structure;
+ info_ro "build_spot_structure"; build_spot_structure;
+ info_ro "build_embedding_tree"; build_embedding_tree;
+ info_ro "color_paths"; color_paths;
+ info_ro "recombine_edges"; recombine_edges(R);
+endgroup
+enddef;
+% ---
+% E-S MACROS:
+% ---
+vardef make_join@#(expr pa,pb)=
+ save kind_; string kind_; kind_:=str @#; if kind_="": kind_:="0" fi;
+ if (kind_<>"0") and (kind_<>"1"):
+ errhelp "Will use default.";
+ errmessage "ES ERROR: don't know how to join";
+ kind_:="0";
+ fi
+ if distance(point length(pa) of pa,point 0 of pb)<epsil.dist:
+ if (point length(pa) of pa)<>(point 0 of pb):
+ info_es "Points " & dec_pair(point length(pa) of pa) &
+ " and " & dec_pair(point 0 of pb) & " joined";
+ if (tracingexpanding>0) and (proofing>0):
+ makelabel.lft.nodot("joined",point length(pa) of pa);
+ fi
+ fi
+ pa && pb
+ elseif kind_="0":
+ if miter_size<=0: % a special case, isn't it?
+ pa--pb
+ else:
+ save ta_,tb_,za_,da_,zb_,db_,zc_,zd_,ze_,zf_;
+ pair za_,da_,zb_,db_,zc_,zd_,ze_,zf_;
+ za_=point length(pa) of pa;
+ da_=direction length(pa) of pa;
+ zb_=point 0 of pb;
+ db_=direction 0 of pb;
+ zc_=whatever[za_,za_+da_]=whatever[ze_,ze_+(zb_-za_)];
+ zd_=whatever[zb_,zb_+db_]=whatever[ze_,ze_+(zb_-za_)];
+ ze_=.5[za_,zb_]+miter_size*(unitvector(da_-db_));
+% we used to check |turningnumber(za_--zc_--zd_--zb_--cycle)|, but it was
+% not sufficiently robust
+ (ta_,tb_)=(za_--zc_) intersectiontimes (zd_--zb_);
+ if ta_<0: % |miter_size| in force:
+ pa
+ if distance(point length(pa) of pa,zc_)>=epsil.dist: --zc_ fi
+ if (distance(zc_,zd_)>=epsil.dist)
+ and (distance(point 0 of pb,zd_)>=epsil.dist): --zd_ fi
+ --pb
+ else:
+ zf_:=point ta_ of (za_--zc_);
+ if abs(zf_-.5[za_,zb_])>abs(ze_-.5[za_,zb_]): % |miter_size| in force:
+ pa
+ if distance(point length(pa) of pa,zc_)>=epsil.dist: --zc_ fi
+ if (distance(zc_,zd_)>=epsil.dist)
+ and (distance(point 0 of pb,zd_)>=epsil.dist): --zd_ fi
+ --pb
+ else:
+ pa
+ if (distance(point length(pa) of pa,zf_)>=epsil.dist)
+ and (distance(point 0 of pb,zf_)>=epsil.dist): --zf_ fi
+ --pb
+ fi
+ fi
+ fi
+ elseif kind_="1":
+ pa{direction length(pa) of pa}..{direction 0 of pb}pb
+ fi
+enddef;
+% ---
+vardef make_cyclic_join@#(expr p)=
+ save kind_; string kind_; kind_:=str @#; if kind_="": kind_:="0" fi;
+ if (kind_<>"0") and (kind_<>"1"):
+ errhelp "Will use default.";
+ errmessage "ES ERROR: don't know how to join";
+ kind_:="0";
+ fi
+ if distance(point length(p) of p,point 0 of p)<epsil.dist:
+ if (point length(p) of p)<>(point 0 of p):
+ info_es "Points " & dec_pair(point length(p) of p) &
+ " and " & dec_pair(point 0 of p) & " joined (cycle)";
+ if (tracingexpanding>0) and (proofing>0):
+ makelabel.lft.nodot("joined (cycle)",point length(p) of p);
+ fi
+ fi
+ make_cycle(p)
+ elseif kind_="0":
+ if miter_size<=0: % a special case, isn't it?
+ p--cycle
+ else:
+ save ta_,tb_,za_,da_,zb_,db_,zc_,zd_,ze_,zf_; pair za_,da_,zb_,db_,zc_,zd_,ze_,zf_;
+ za_=point length(p) of p; da_=direction length(p) of p;
+ zb_=point 0 of p; db_=direction 0 of p;
+ zc_=whatever[za_,za_+da_]=whatever[ze_,ze_+(zb_-za_)];
+ zd_=whatever[zb_,zb_+db_]=whatever[ze_,ze_+(zb_-za_)];
+ ze_=.5[za_,zb_]+miter_size*(unitvector(da_-db_));
+% we used to check |turningnumber(za_--zc_--zd_--zb_--cycle)|, but it was
+% not sufficiently robust
+ (ta_,tb_)=(za_--zc_) intersectiontimes (zd_--zb_);
+ if ta_<0: % |miter_size| in force:
+ p
+ if distance(point length(p) of p,zc_)>=epsil.dist: --zc_ fi
+ if (distance(zc_,zd_)>=epsil.dist)
+ and (length((point 0 of p)-zd_)>=epsil.dist): --zd_ fi
+ --cycle
+ else:
+ zf_:=point ta_ of (za_--zc_);
+ if abs(zf_-.5[za_,zb_])>abs(ze_-.5[za_,zb_]): % |miter_size| in force:
+ p
+ if distance(point length(p) of p,zc_)>=epsil.dist: --zc_ fi
+ if (distance(zc_,zd_)>=epsil.dist)
+ and (distance(point 0 of p,zd_)>=epsil.dist): --zd_ fi
+ --cycle
+ else:
+ p
+ if (distance(point length(p) of p,zf_)>=epsil.dist)
+ and (distance(point 0 of p,zf_)>=epsil.dist): --zf_ fi
+ --cycle
+ fi
+ fi
+ fi
+ elseif kind_="1":
+ p{direction length(p) of p}..{direction 0 of p}cycle
+ fi
+enddef;
+% ---
+vardef make_end@#(expr pr,pl) =
+ save kind_; string kind_; kind_:=str @#; if kind_="": kind_:="0" fi;
+ if (kind_<>"0") and (kind_<>"1"):
+ errhelp "Will use default.";
+ errmessage "ES ERROR: don't know how to end";
+ kind_:="0";
+ fi
+ if kind_="0": pr--pl--cycle
+ elseif kind_="1":
+ save za_,zb_; pair za_,zb_;
+ za_=1/2[point length(pr) of pr,point 0 of pl]
+ +(1/2((point length(pr) of pr)-(point 0 of pl)) rotated 90);
+ zb_=1/2[point length(pl) of pl,point length 0 of pr]
+ +(1/2((point length(pl) of pl)-(point 0 of pr)) rotated 90);
+ pr{direction length(pr) of pr}..za_..{direction 0 of pl}pl
+ {direction length(pl) of pl}..zb_..{direction 0 of pr}cycle
+ fi
+enddef;
+% ---
+vardef opt_tensions(expr p,b) =
+% for a given B\'ezier segment |p| and a distance |b|, an optimal pair of
+% `tensions' $(\alpha,\beta)$ is found using least square method such that
+% |bez_edge|$(p,b,\alpha,\beta)$ (see below) approximates the edge of
+% a circular pen of diameter |b| traversing |p| (more on the employed
+% method be found in the article of B. Jackowski and M. Ry\'cko:
+% ``Labyrinth of \MF paths in outline,'' proceedings of the 8th European
+% \TeX Conference, Sept. 26--30, 1994, Gda\'nsk, Poland)
+%
+ save alpha_,beta_,gx_,gy_,n_,t_,ta_,tb_,tc_,td_,u_,v_,nu_,nv_,x_,y_;
+ numeric alpha_,beta_,n_,ta_,tb_,tc_,td_,
+ gx_[\\],gy_[\\],gx_.alpha[\\],gy_.alpha[\\],gx_.beta[\\],gy_.beta[\\],
+ u_.x,u_.y,v_.x,v_.y,nu_.x,nu_.y,nv_.x,nv_.y,
+ x_[\\],y_[\\];
+ n_:=5; % perhaps for |n_|$=\infty$ algebraic formulas can be derived, but...
+ (u_.x,u_.y)=(postcontrol 0 of p)-(point 0 of p);
+ (v_.x,v_.y)=(precontrol 1 of p)-(point 1 of p);
+ (nu_.x,nu_.y)=unitvector(u_.x,u_.y); (nv_.x,nv_.y)=unitvector(v_.x,v_.y);
+ for t_:=0 upto n_:
+ (x_[t_],y_[t_])=(point t_/n_ of p)+b*((udir t_/n_ of p) rotated -90);
+ endfor
+ for t_:=1 upto n_-1:
+ td_:=t_/n_; ta_:=1-td_; tb_:=3ta_*ta_*td_; tc_:=3ta_*td_*td_;
+ ta_:=ta_*ta_*ta_; td_:=td_*td_*td_;
+ gx_[t_]=ta_*x_0+tb_*(x_0+alpha_*u_.x)+tc_*(x_[n_]+beta_*v_.x)+td_*x_[n_];
+ gx_.alpha[t_]=tb_*nu_.x; gx_.beta[t_]=tc_*nv_.x;
+ gy_[t_]=ta_*y_0+tb_*(y_0+alpha_*u_.y)+tc_*(y_[n_]+beta_*v_.y)+td_*y_[n_];
+ gy_.alpha[t_]=tb_*nu_.y; gy_.beta[t_]=tc_*nv_.y;
+ endfor
+ 0=0 for t_:=1 upto n_-1:
+ +((gx_[t_]-x_[t_])*gx_.alpha[t_]+(gy_[t_]-y_[t_])*gy_.alpha[t_])/n_
+ endfor;
+ 0=0 for t_:=1 upto n_-1:
+ +((gx_[t_]-x_[t_])*gx_.beta[t_]+(gy_[t_]-y_[t_])*gy_.beta[t_])/n_
+ endfor;
+%| (u_.x,u_.y)=(postcontrol 0 of p)-(point 0 of p);|
+%| (v_.x,v_.y)=(precontrol 1 of p)-(point 1 of p);|
+%| ta_:=1/4length((u_.x,u_.y))+1/4length((v_.x,v_.y))|
+%| +1/4length((postcontrol 0 of p)-(precontrol 1 of p))|
+%| +1/4length((point 0 of p)-(point 1 of p));|
+%| message "accuracy=" & decimal|
+%| (0+for t_:=1 upto n_-1:+(((gx_[t_]-x_[t_])++(gy_[t_]-y_[t_]))/ta_)/n_|
+%| endfor);|
+%| message " alpha=" & decimal(alpha_) & " beta=" & decimal(beta_);|
+%| for t_:=0 upto n_: fill fullcircle scaled 3 shifted (x_[t_],y_[t_]); endfor|
+%| for t_:=1 upto n_-1: makelabel("g" & decimal(t_),(gx_[t_],gy_[t_])); endfor|
+ (alpha_,beta_)
+enddef;
+% ---
+vardef bez_edge(expr p,b,uv) =
+ save za_,zb_,u_,v_; pair za_,zb_; u_:=xpart(uv); v_:=ypart(uv);
+ za_=b*((udir 0 of p) rotated -90); zb_=b*((udir 1 of p) rotated -90);
+ ((point 0 of p)+za_) .. controls (u_[point 0 of p,postcontrol 0 of p]+za_)
+ and (v_[point 1 of p,precontrol 1 of p]+zb_) .. ((point 1 of p)+zb_)
+enddef;
+% ---
+def remove_global_loops(suffix E) =
+begingroup
+% warning: we don't trust too much in the results of ex. 14.17 from
+% The \MF{}book, hence a ``par force'' approach; there still exist
+% weird cases (e.g., local loops) which remain unsolved, but in practice
+% the following algorithm should suffice:
+ save opt_,ta_,tb_; pair opt_;
+ opt_:=(0,length(E));
+ for i_:=0 upto length(E)-1:
+ for j_:=i_+2 upto length(E)-1:
+ numeric ta_,tb_;
+ (ta_,tb_)=(subpath (i_,i_+1) of E)
+ intersectiontimes (subpath (j_,j_+1) of E);
+ if (ta_>0) and ((ta_+i_)>xpart(opt_)) and ((tb_+j_)<ypart(opt_)):
+ opt_:=(ta_+i_,tb_+j_);
+ fi
+ endfor
+ endfor
+ if xpart(opt_)>0:
+ E:=make_cycle(subpath(xpart(opt_),ypart(opt_)) of E);
+ fi
+endgroup
+enddef;
+% ---
+vardef make_edge@#(expr p,b)=
+ save E_,e_,ta_,tb_,tc_,td_; path E_,e_[\\];
+ for i_:=0 upto length(p)-1:
+ E_:=subpath (i_,i_+1) of p; e_[i_]=bez_edge(E_,b,opt_tensions(E_,b));
+ endfor
+ E_:=e_0;
+ for i_:=1 upto length(p)-1:
+ if (length(p)=2) and (cycle p): % this is a peculiar case, indeed!
+ numeric ta_,tb_,tc_,td_;
+ (ta_,tb_)=E_ intersectiontimes e_[i_];
+ (1-tc_,1-td_)=reverse(E_) intersectiontimes reverse(e_[i_]);
+ if ta_>=0:
+ E_:=(subpath(min(ta_,tc_),max(ta_,tc_)) of E_)
+ && (subpath(min(tb_,td_),max(tb_,td_)) of e_[i_]);
+ else: E_:=make_join@#(E_,e_[i_]);
+ fi
+ else:
+ numeric ta_,tb_;
+ (ta_,tb_)=(subpath(length(E_)-1,length(E_)) of E_)
+ intersectiontimes e_[i_];
+ if ta_>=0:
+ E_:=(subpath (0,length(E_)-1+ta_) of E_)
+ && (subpath(tb_,1) of e_[i_]);
+ else: E_:=make_join@#(E_,e_[i_]); fi
+ fi
+ endfor
+ if cycle p:
+ remove_global_loops(E_);
+ if not (cycle E_): E_:=make_cyclic_join@#(E_); fi
+ fi
+ E_
+enddef;
+% ---
+def expand_stroke(text P)(expr b) suffix R =
+begingroup interim autorounding:=0;
+ numeric PATH_.num; path PATH_[\\];
+ PATH_.num:=0; for P_:=P: PATH_[incr PATH_.num]:=touch_path(P_); endfor
+ if not path R0: numeric R.num; path R[\\]; fi
+ if (unknown R.num) or (unknown append_results): R.num:=0; fi
+ if unknown join_kind: save join_kind; join_kind=0; fi
+ if unknown end_kind: save end_kind; end_kind=0; fi
+ for i_:=1 upto PATH_.num:
+ if not cycle PATH_[i_]:
+ R[incr R.num]:=make_end[end_kind]
+ (make_edge[join_kind](PATH_[i_],b),
+ reverse make_edge[join_kind](PATH_[i_],-b));
+ else:
+ R[incr R.num]:=make_edge[join_kind](PATH_[i_],b);
+ R[incr R.num]:=reverse make_edge[join_kind](PATH_[i_],-b);
+ fi
+ endfor
+ for i_:=1 upto R.num: R[i_]:=clean_path(clean_path(R[i_])); endfor
+endgroup
+enddef;
+% ---
+def change_weight(text P)(expr b) suffix R =
+begingroup interim autorounding:=0;
+ numeric PATH_.num; path PATH_[\\];
+ PATH_.num:=0; for P_:=P: PATH_[incr PATH_.num]:=touch_path(P_); endfor
+ if not path R0: numeric R.num; path R[\\]; fi
+ if (unknown R.num) or (unknown append_results): R.num:=0; fi
+ if unknown join_kind: save join_kind; join_kind=0; fi
+ for i_:=1 upto PATH_.num:
+% non-cyclic paths are ignored
+ if cycle PATH_[i_]: R[incr R.num]:=make_edge[join_kind](PATH_[i_],b); fi
+ endfor
+endgroup
+enddef;
+% ---
+def info_ro expr s = if tracingremoving>0: message s; message ""; fi enddef;
+def info_es expr s = if tracingexpanding>0: message s; message ""; fi enddef;
+% ---
+% DEFAULTS:
+% ---
+def roex_default text t =
+ forsuffixes S_:=t:
+ if str S_ = "good_colors":
+% the formula |good_colors(p,q) and good_colors(q,p)| must be |false|!
+ vardef good_colors(expr i,o) = ((i>=1) and (o<=0)) enddef;
+ elseif str S_ = "touch_path":
+ vardef touch_path(expr p) = p enddef;
+ elseif str S_ = "background_color": background_color:=0;
+ elseif str S_ = "miter_size":
+ miter_size:=10pixels_per_inch/72; % i.e., 10bp
+% incidentally, |10bp| would convert to |10.00002| during export at |300dpi|
+ elseif str S_ = "epsil.ang": epsil.ang:=1/10; % in degrees
+ elseif str S_ = "epsil.dist": epsil.dist:=1/10pt; % ca |2/5|pxl at |300dpi|
+ elseif str S_ = "epsil.time": epsil.time:=1/100;
+ elseif str S_ = "epsil.len": epsil.len:=1/1000; % used in |turn_ang|
+ elseif str S_ = "max_idx": max_idx:=125;
+ elseif str S_ = "enc":
+% |enc| is a prefix of a data structure used in checking tangent
+% points and searching for the leftmost edge; |enc.pth| is in both
+% cases scaled differently
+ vardef enc.pth = fullcircle enddef; enc.len:=length(enc.pth);
+ fi
+ endfor
+enddef;
+%
+roex_default good_colors, touch_path, background_color, miter_size,
+ epsil.ang, epsil.dist, epsil.time, epsil.len, max_idx, enc;
+% ---
+numeric append_results; % initially unknown
+newinternal tracingleftmost; tracingleftmost:=0;
+newinternal tracingremoving; tracingremoving:=0;
+newinternal tracingexpanding; tracingexpanding:=0;
+% ---
+endinput
+%%\end
+