Re: Re: Re[4]: MetaPost/Fun question

[Giuseppe Bilotta <gip.bilotta@iol.it>]

[-- Attachment #1: Type: text/plain, Size: 456 bytes --]

Monday, March 24, 2003 Guy Worthington wrote:

GW> Giuseppe Bilotta wrote
>> Emil Hedevang Lohse wrote:

>>> [an algorithm for finding the intersection points of two paths]

>> [a metapost implementation]

GW> Thanks Emil and Giuseppe, a nice example that was well worth the read.

You may want to have a look at the stuff in the attached MetaPost
macroset, which I'm preparing for a better support of the
Eukleides program.

-- 
Giuseppe "Oblomov" Bilotta

[-- Attachment #2: eukleides.mp --]
[-- Type: application/octet-stream, Size: 4613 bytes --]

%D This \MP\ source is supposed to help in the conversion of Eukleides
%D datafiles to \MP\ format.

%D The first macro is findallintersections ; given two paths,
%D it sets one global variable (intersectionsfound) containing the number of
%D found intersections, and two sets of global pairs:
%D \startitemize
%D \item intersections[0] \dots\ intersections[intersectionsfound-1] which
%D contain the times of the intersections point (xpart for ther first curve,
%D ypart for the second curve;
%D \item intersectionpointsp[0] \dots\ intersectionpoints[intersectionsfound-1]
%D which contain the actual points;
%D \stopitemize
%D the behaviour of this macro is controlled by the global parameter
%D intersectiontolerance: setting it too high will skip some intersections,
%D setting it too low will find too many, possibly entering an infinite loop;
%D this is especially the case for tangent, osculant or hyperosculant paths.

%D The macro works by finding an intersection point, removing the slice of the
%D first path upto the first intersection point (plus tolerance), and reiterate
%D until either no more intersection points are found or the first paths gets
%D exhausted (this check is necessary because of infinite loops that can
%D happen otherwise.

newinternal intersectiontolerance ; intersectiontolerance := eps ;

def findallintersections (expr p, q) = 
% we start by resetting the results
pair intersections[], intersectionpoints[] ; intersectionsfound := 0 
if (path p and path q):
begingroup ;
save j, tp ; path tp; tp := p ;
j:= 0 ;
intersections[j] = tp intersectiontimes q ;
forever:
exitif (xpart intersections[j] = -1) or (ypart intersections[j] = -1) or 
       (xpart intersections[j] + intersectiontolerance > length p) ;
intersectionpoints[j] = point ypart intersections[j] of q ;
j:=j+1 ;
tp := subpath (xpart intersections[j-1] + intersectiontolerance, length p) of p ;
intersections[j] := tp intersectiontimes q ;
% watch the trick to ensure the xpart will be the time on the real path
intersections[j] := (xpart (p intersectiontimes (point xpart intersections[j] of tp)),
                     ypart intersections[j]);
endfor ;
intersectionsfound := j ;
endgroup ;
fi ; 
enddef ;

%D We then create an auxiliary macro to determine the time of a point
%D on a curve; I'm not aware of any \MP\ primitives that do this, so \dots

primarydef p timedon q =
(xpart (q intersectiontimes p))
enddef ;

%D I hate the counter-intuitiveness of fullcircle not having radius 1, so:

path trigcircle ; trigcircle := fullcircle scaled 2 ;

%D A useful shorthand:

def point_at_angle (expr degrees) =
(cosd(degrees), sind(degrees))
enddef ;

%D The next macro returns a circular arc. It accepts two parameters
%D (starting angle, ending angle) in degrees. Angles are assumed to be
%D oriented.

def arc (expr starting, ending) =
begingroup 
save startingpoint, endingpoint, startingtime, endingtime ;
pair startingpoint, endingpoint ;
startingpoint := point_at_angle(starting) ;
endingpoint   := point_at_angle(ending) ;
startingtime  := startingpoint timedon trigcircle ;
endingtime    := endingpoint   timedon trigcircle ;
% watch the trick to ensure that arcs crossing the origin
% are drawn correctly
if startingtime > endingtime:
arc (starting+10, ending+10) rotated -10
else:
subpath (startingtime, endingtime) of trigcircle 
fi
endgroup
enddef ;

%D These are really not needed by Eukleides, but anyway \dots

def radiusat (expr degrees) =
(origin -- point_at_angle(degrees))
enddef ;

def circularsector (expr starting, ending) =
(origin -- arc(starting, ending) -- cycle)
enddef ;

%D Back to Eukleides; a dotmark with a plus shape

picture dot_mark_plus; dot_mark_plus := nullpicture ;
addto dot_mark_plus doublepath (-3,0)--(3,0) ;
addto dot_mark_plus doublepath (0,3)--(0,-3) ;

%D In Eukleides we can set labels centered at a given distance and angle
%D in proximity of a given point (aka psfig's \tex{uput})

color labelcolor ; labelcolor := black ;

def pslabel( expr t,  p,  d,  a) =
% t: text; p: point; d: distance; a: angle
begingroup
save thislabel, thispoint, drawingcolor ;
picture thislabel ; pair thispoint ; color drawingcolor ;
if picture t:  thislabel = t
else:    thislabel = t infont defaultfont scaled defaultscale
fi; 
xpart thispoint = d*cosd(a) - xpart center thislabel * cosd(a) ;
ypart thispoint = d*sind(a) - ypart center thislabel * sind(a) ;
drawingcolor := labelcolor ;
label(thislabel, p + thispoint) ;
endgroup
enddef ;


endinput ;