Why two intersected paths can't buildcycle in MetaPost?

Question

Suppose that I have two path:

beginfig(0)
u:=10pt;
vardef randomcirc(expr O, r, w) = 
  save p,i; pair p; numeric i;
    p=O+right*r;
    p for i=1 upto 9: .. (p rotatedaround(O, 36i))+(uniformdeviate w, uniformdeviate w) endfor .. cycle
enddef;
path pat[];
pair p[];
p0:=(0,5u); p3:=(-8u,0);
pat0:=randomcirc(p0,8u,.5u);
pat3:=randomcirc(p3,5u,0.4u);
fill buildcycle(pat0, pat3) withcolor red;
draw pat0; draw pat3;
endfig;

why their intersection are not filled with red color?

I think it's due to the precision of MetaPost too low, how to fix it?

Answer 1

The plain.mp implementation of buildcycle fails to find the correct overlap between two cyclic paths if the beginning (point 0) of one of the paths lies inside the other path. For example:

beginfig(3);
  path A, B; picture p[];
  A = fullcircle scaled 2.5cm; 
  B = fullcircle scaled 2cm shifted (1cm,0);

  p1 = image(fill buildcycle(A,B) withcolor .8[blue,white]; drawarrow A; drawarrow B;);

  A := A rotated 180;
  p2 = image(fill buildcycle(A,B) withcolor .8[blue,white]; drawarrow A; drawarrow B;);

  B := B rotatedabout(center B,180);
  p3 = image(fill buildcycle(A,B) withcolor .8[blue,white]; drawarrow A; drawarrow B;);

  A := A rotated 180;
  p4 = image(fill buildcycle(A,B) withcolor .8[blue,white]; drawarrow A; drawarrow B;);

  for i=1 upto 4: draw p[i] shifted (120i,0); label(decimal i, (7mm+120i,0)); endfor
endfig;

produces

In each subfigure, the arrow head is pointing at the beginning of the path. You can see that buildcycle only works in situation 2 where neither point 0 lies within the other cyclic path. Note that in situation 4 we get the union of the two paths.

Here's a way to correct this fault. First you need a function to determine whether a given point lies within a cyclic path. Following Sedgewick's Algorirthms in C, you can write:

% is point "p" inside cyclic path "ring" ?
vardef inside(expr p, ring) = 
  save t, count, test_line;
  count := 0;
  path test_line;
  test_line = p -- (infinity, ypart p);
  for i = 1 upto length ring:
     t := xpart(subpath(i-1,i) of ring intersectiontimes test_line);
     if ((0<=t) and (t<1)): count := count + 1; fi
  endfor
  odd(count)
  enddef;

And here is an overlap function that replaces buildcycle in this special case of two overlapping cyclic paths.

vardef front_half primary p = subpath(0, 1/2 length p) of p enddef;
vardef back_half  primary p = subpath(1/2 length p, length p) of p enddef;

% a and b should be cyclic paths...
vardef overlap(expr a, b) = 
  save p, q;
  boolean p, q;
  p = inside(point 0 of a, b);
  q = inside(point 0 of b, a);
  if ((not p) and (not q)): 
    buildcycle(a,b)
  elseif (not p): 
    buildcycle(front_half b, a, back_half b)
  elseif (not q): 
    buildcycle(front_half a, b, back_half a)
  else: 
    buildcycle(front_half a, back_half b, front_half b, back_half a)
  fi
  enddef;

The basic idea is that if neither point 0 is inside the other path, then you just call buildcycle, otherwise you split the two cycles up into an appropriate sequence of half cycles. Replacing buildcycle in the above example with overlap gives this:

And replacing buildcycle in the OP example with overlap gives this:

which is probably closer to what was wanted in the first place.

Furthermore, if the two cyclic paths are running in opposite directions the behaviour of buildcycle is different again. Here's an extended version of the first example; in the second row, the larger path is running backwards:

To deal with this possibility you can make overlap more robust by using the useful counterclockwise function (from plain.mp) which returns a copy of a cyclic path running counterclockwise. Here is the improved version of the overlap function.

vardef overlap(expr a, b) = 
  save p, q, A, B;
  boolean p, q;
  p = not inside(point 0 of a, b);
  q = not inside(point 0 of b, a);
  path A, B;
  A = counterclockwise a;
  B = counterclockwise b;
  if (p and q): 
    buildcycle(A,B)
  elseif p: 
    buildcycle(front_half B, A, back_half B)
  elseif q: 
    buildcycle(front_half A, B, back_half A)
  else: 
    buildcycle(front_half A, back_half B, front_half B, back_half A)
  fi
  enddef;

which produces this in my extended example.

Answer 2

It works here if I change slightly the definition of pat3, making pat3 start at the "left" (so to speak) of the circle and not at the "right" (to the contrary of pat0). For this, I defined another vardef macro (randomcirc_bis)

beginfig(0)
    u:=10pt;
    vardef randomcirc(expr O, r, w) = 
        save p,i; pair p; numeric i;
         p=O+right*r;
        p for i=1 upto 9: .. (p rotatedaround(O, 36i))+(uniformdeviate w, uniformdeviate w)                
        endfor .. cycle
    enddef;
    %
    vardef randomcirc_bis(expr O, r, w) = 
        save p,i; pair p; numeric i;
        p=O+left*r;
        p for i=1 upto 9: .. (p rotatedaround(O, 36i))+(uniformdeviate w, uniformdeviate w)  
        endfor .. cycle
    enddef;
%
    path pat[];
    pair p[];
    p0:=(0,5u); p3:=(-8u,0);
    pat0:=randomcirc(p0,8u,.5u);
    pat3:=randomcirc_bis(p3,5u,0.4u);
    pat4 = buildcycle(pat0, pat3);
    fill pat4 withcolor red;
    draw pat0; draw pat3 dashed evenly;
    %
    dotlabel.ulft("1", pat0 intersectionpoint pat3);
    dotlabel.lrt("2", pat3 intersectionpoint pat0);
    label.rt("aa", point 0 of pat0);
    label.lft("b", point 0 of pat3);
endfig;
end.

In this drawing I tried to reproduce exactly the modus operandi of buildcycle as explained in the MetaPost manual, p. 30 (section 9.1, "Building Cycles"). Citation:

The buildcycle macro detects the two intersections labeled 1 and 2 in Figure 24b. Then it constructs the cyclic path shown in bold in the figure by going forward along path aa from intersection 1 to intersection 2 and then forward around the counter-clockwise path b back to intersection 1.

If I understand the macro intersectionpoint correctly, point 1 is in fact aa intersectionpoint b and it is the first intersection that aa (pat0) encounters with b (pat3).

Similarly point 2 is b intersectiopoint aa, and it is also the first intersection that b (pat3) encounters with aa (pat0).

To reproduce this scheme it has been necessary to start the path pat3 at the left. Otherwise the first intersection that pat3=b would have met with pat0=aa would have been point 1. The two intersection points would have been the same, and in that case buildcycle fails.

I hope I have not been too confusing: English is not my native language…

UPDATE

I have simplified the code. My new vardef macro randomcircbis wasn't necessary. It was enough to replace the line pat3:=randomcirc(p3,5u,0.4u);in the original code by pat3:=randomcirc(p3,-5u,0.4u);. The resulting figure is the same, of course.

beginfig(0)
    u:=10pt;
    vardef randomcirc(expr O, r, w) = 
        save p,i; pair p; numeric i;
         p=O+right*r;
        p for i=1 upto 9: .. (p rotatedaround(O, 36i))+(uniformdeviate w, uniformdeviate w)                
        endfor .. cycle
    enddef;
%
    path pat[];
    pair p[];
    p0:=(0,5u); p3:=(-8u,0);
    pat0:=randomcirc(p0,8u,.5u);
    pat3:=randomcirc(p3,-5u,0.4u);
    pat4 = buildcycle(pat0, pat3);
    fill pat4 withcolor red;
    draw pat0; draw pat3 dashed evenly;
    %
    dotlabel.ulft("1", pat0 intersectionpoint pat3);
    dotlabel.lrt("2", pat3 intersectionpoint pat0);
    label.rt("aa", point 0 of pat0);
    label.lft("b", point 0 of pat3);
endfig;
end.

Answer 3

Here is another implementation of the same idea as in fpast's answer. Instead of starting with a point on the left, you can just rotate the original circle by 180 degrees. (I also simplified the code for randomized circle a bit).

\startMPdefinitions
  vardef randomcircle(expr O, r, w, a) = 
    ((fullcircle rotated  a scaled 2r) randomized w shifted O)
  enddef;
\stopMPdefinitions
\starttext

\startMPpage[offset=2mm]
  u:=10pt;
  path pat[];
  pair p[];

  p0:=(0,5u); p3:=(-8u,0);

  pat0 = randomcircle(p0, 8u, 0.4u, 0);
  pat3 = randomcircle(p3, 5u, 0.4u, 180);

  fill buildcycle(pat0, pat3) withcolor red;

  draw pat0;
  draw pat3;

  dotlabel.ulft("1", pat0 intersectionpoint pat3);
  dotlabel.lrt("2", pat3 intersectionpoint pat0);
  label.rt("aa", point 0 of pat0);
  label.lft("b", point 0 of pat3);
\stopMPpage

\stoptext

which gives