MetaPost/Fun question

MetaPost/Fun question

[Giuseppe Bilotta]

Hello,

I would like to know if it is possible in MetaPost (in particular
MetaFun) to get all the possible intersection points of the paths,
with some kind of simple command, with a syntax like, say:

A[] := path1 intersections path2

(where the resulting list could be either a list of points or a
list of time pairs, like for intersectiontimes).

Another possible syntax could be

path3 := path1 intersections path2

which would give a path for which times 0, 1, 2, 3, ... give the
first, second, third, fourth, ... intersection point.

Anything like this?

-- 
Giuseppe "Oblomov" Bilotta

Re: MetaPost/Fun question

[Hans Hagen]

At 09:59 PM 3/22/2003 +0100, you wrote:
>Hello,
>
>I would like to know if it is possible in MetaPost (in particular
>MetaFun) to get all the possible intersection points of the paths,
>with some kind of simple command, with a syntax like, say:
>
>A[] := path1 intersections path2
>
>(where the resulting list could be either a list of points or a
>list of time pairs, like for intersectiontimes).
>
>Another possible syntax could be
>
>path3 := path1 intersections path2
>
>which would give a path for which times 0, 1, 2, 3, ... give the
>first, second, third, fourth, ... intersection point.
>
>Anything like this?

not that i know, but assuming that those path have some distance, i can 
imagine dividing a path in parts and determining the intersectionpoints iof 
the pieces;

delta := .25 ;

for i=0 step delta upto length(p)-delta :
   if intersection_found(1,subpath(i,i+.1) of p) :
     store point
   fi ;
endfor ;

.....


Hans
-------------------------------------------------------------------------
                                   Hans Hagen | PRAGMA ADE | pragma@wxs.nl
                       Ridderstraat 27 | 8061 GH Hasselt | The Netherlands
  tel: +31 (0)38 477 53 69 | fax: +31 (0)38 477 53 74 | www.pragma-ade.com
-------------------------------------------------------------------------
                        information: http://www.pragma-ade.com/roadmap.pdf
                     documentation: http://www.pragma-ade.com/showcase.pdf
-------------------------------------------------------------------------

Re[2]: MetaPost/Fun question

[Giuseppe Bilotta]

Saturday, March 22, 2003 Hans Hagen wrote:

>>I would like to know if it is possible in MetaPost (in particular
>>MetaFun) to get all the possible intersection points of the paths,

[snip]

>>Anything like this?

HH> not that i know, but assuming that those path have some distance, i can
HH> imagine dividing a path in parts and determining the intersectionpoints iof
HH> the pieces;

HH> delta := .25 ;

HH> for i=0 step delta upto length(p)-delta :
HH>    if intersection_found(1,subpath(i,i+.1) of p) :
HH>      store point
HH>    fi ;
HH> endfor ;

I'll see if I can make something of this ... (wow, looks like I'm
gonna need to learn MetaPost too ...)

-- 
Giuseppe "Oblomov" Bilotta

Re: Re[2]: MetaPost/Fun question

[Emil Hedevang Lohse]

Giuseppe Bilotta <gip.bilotta@iol.it> writes:

> Saturday, March 22, 2003 Hans Hagen wrote:
> 
> >>I would like to know if it is possible in MetaPost (in particular
> >>MetaFun) to get all the possible intersection points of the paths,
> 
> [snip]
> 
> >>Anything like this?
> 
> HH> not that i know, but assuming that those path have some distance, i can
> HH> imagine dividing a path in parts and determining the intersectionpoints iof
> HH> the pieces;
> 
> HH> delta := .25 ;
> 
> HH> for i=0 step delta upto length(p)-delta :
> HH>    if intersection_found(1,subpath(i,i+.1) of p) :
> HH>      store point
> HH>    fi ;
> HH> endfor ;
> 
> I'll see if I can make something of this ... (wow, looks like I'm
> gonna need to learn MetaPost too ...)

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

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

When done this way, the intersection times will be sorted according to
the first path.
 

-- 
Emil Hedevang Lohse <http://home.imf.au.dk/emil/> 

Alle spørgsmål er lige dumme. 
Og spørgsmålet "Kan ænder flyve?" er ikke dumt.

Re[4]: MetaPost/Fun question

[Giuseppe Bilotta]

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

Re: Re[4]: MetaPost/Fun question

[Guy Worthington]

Giuseppe Bilotta wrote
> Emil Hedevang Lohse wrote:

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

> [a metapost implementation]

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

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

[Giuseppe Bilotta]

[-- 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 ;