Re[4]: MetaPost/Fun question

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

Sunday, March 23, 2003 Emil Hedevang Lohse wrote:

EHL> I implemented it by finding some intersectiontime and the recursively
EHL> finding the intersectionstimes of the two subpaths (each with some
EHL> neighbourhood of the first intersection point cut off). For the time
EHL> being I cannot find my code, but here it is in pseude code.

EHL>   let p and q be paths;
EHL>   let (s,t) be a pair of intersection times of p and q;
EHL>   recursively find the intersection times of the 
EHL>     subpath (0, s - small number) of p and q;
EHL>   store (s,t);
EHL>   recursively find the intersection times of the 
EHL>     subpath (s + small number, infinity) of p and q;
EHL>   feel happy;

EHL> When done this way, the intersection times will be sorted according to
EHL> the first path.
 
Ok, I came up with something which might be useful; Hans, you may
consider adding it to the MetaFun core:

newinternal intersectiontolerance ; intersectiontolerance := 4*eps ;
def findallintersections (expr p, q) = 
pair intersections[], intersectionpoints[] ; intersectionsfound := 0 ; % reset intersections
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 ;
% watch the trick to ensure that the times are
% relative to the original paths
intersections[j] := tp intersectiontimes q ;
intersections[j] := (xpart (p intersectiontimes ((point xpart intersections[j] of tp)--cycle)),
                     ypart intersections[j]);
endfor ;
intersectionsfound := j ;
endgroup ;
fi ; 
enddef ;


Some notes:

* this is my first MetaPost macro, so I expect it to be slow,
inefficient, and buggy; feel free to comment

* USAGE: findallintersections(path, path) reinitializes the
following global expressions:

+ intersectionsfound, a numeric containing the number of
intersections

+ intersections[0] .. intersections[intersectionsfound-1], pairs
containing the time of the intersections on the two paths

+ intersectionpoints[0] ..
intersectionpoints[intersectionsfound-1], points of intersection
of the two paths

* PARAMETER: there is a global parameter intersectiontolerance:
setting it too low will send the macro in a neverending loop in case of
tangent curves (example:

beginfig(1)
path c[] ;
c1:= fullcircle scaled 10cm ;
c2:= fullsquare scaled 10cm ;

draw c1 withcolor blue; draw c2 withcolor blue;

loggingall ;
findallintersections (c1, c2) ;
show intersectionsfound ;
for i= 0 upto intersectionsfound-1 : drawpoint intersectionpoints[i] ; endfor ;

endfig

will fail for intersectiontolerance:=eps ;

Comments? Ideas? Suggestions?

-- 
Giuseppe "Oblomov" Bilotta