%%%% %%%% 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_0): res:=1; % $a \subset b$ elseif (panb_<>0) and (napb_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_)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 "") 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_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_=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)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_])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_.loc0: 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)(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)(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_)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