path: root/Build/source/utils/asymptote/base/bezulate.asy

// Bezier triangulation routines written by Orest Shardt, 2008.

private real fuzz=sqrtEpsilon;
real duplicateFuzz=1e-3; // Work around font errors.
real maxrefinements=7;

private real[][] intersections(pair a, pair b, path p)
{
  pair delta=fuzz*unit(b-a);
  return intersections(a-delta--b+delta,p,fuzz);
}

int countIntersections(path[] p, pair start, pair end)
{
  int intersects=0;
  for(path q : p)
    intersects += intersections(start,end,q).length;
  return intersects;
}

path[][] containmentTree(path[] paths)
{
  path[][] result;
  for(path g : paths) {
    // check if current curve contains or is contained in a group of curves
    int j;
    for(j=0; j < result.length; ++j) {
      path[] resultj=result[j];
      int test=inside(g,resultj[0],zerowinding);
      if(test == 1) {
        // current curve contains group's toplevel curve;
        // replace toplevel curve with current curve
        resultj.insert(0,g);
        // check to see if any other groups are contained within this curve
        for(int k=j+1; k < result.length;) {
          if(inside(g,result[k][0],zerowinding) == 1) {
            resultj.append(result[k]);
            result.delete(k);
          } else ++k;
        }
        break;
      } else if(test == -1) {
        // current curve contained within group's toplevel curve
        resultj.push(g);
        break;
      }
    }
    // create a new group if this curve does not belong to another group
    if(j == result.length)
      result.push(new path[] {g});
  }
  return result;
}

bool isDuplicate(pair a, pair b, real relSize)
{
  return abs(a-b) <= duplicateFuzz*relSize;
}

path removeDuplicates(path p)
{
  real relSize = abs(max(p)-min(p));
  bool cyclic=cyclic(p);
  for(int i=0; i < length(p); ++i) {
    if(isDuplicate(point(p,i),point(p,i+1),relSize)) {
      p=subpath(p,0,i)&subpath(p,i+1,length(p));
      --i;
    }
  }
  return cyclic ? p&cycle : p;
}

path section(path p, real t1, real t2, bool loop=false)
{
  if(t2 < t1 || loop && t1 == t2)
    t2 += length(p);
  return subpath(p,t1,t2);
}

path uncycle(path p, real t)
{
  return subpath(p,t,t+length(p));
}

// returns outer paths
void connect(path[] paths, path[] result, path[] patch)
{
  path[][] tree=containmentTree(paths);
  for(path[] group : tree) {
    path outer = group[0];
    group.delete(0);
    path[][] innerTree = containmentTree(group);
    path[] remainingCurves;
    path[] inners;
    for(path[] innerGroup:innerTree)
      {
        inners.push(innerGroup[0]);
        if(innerGroup.length>1)
          remainingCurves.append(innerGroup[1:]);
      }
    connect(remainingCurves,result,patch);
    real d=2*abs(max(outer)-min(outer));
    while(inners.length > 0) {
      int curveIndex = 0;
      //pair direction=I*dir(inners[curveIndex],0,1); // Use outgoing direction
      //if(direction == 0) // Try a random direction
      //  direction=expi(2pi*unitrand());
      //pair start=point(inners[curveIndex],0);
      
      // find shortest distance between a node on the inner curve and a node
      // on the outer curve
      
      real mindist = d;
      int inner_i = 0;
      int outer_i = 0;
      for(int ni = 0; ni < length(inners[curveIndex]); ++ni)
      {
          for(int no = 0; no < length(outer); ++no)
          {
              real dist = abs(point(inners[curveIndex],ni)-point(outer,no));
              if(dist < mindist)
              {
                  inner_i = ni;
                  outer_i = no;
                  mindist = dist;
              }
          }
      }
      pair start=point(inners[curveIndex],inner_i);
      pair end = point(outer,outer_i);  
      
      // find first intersection of line segment with outer curve
      //real[][] ints=intersections(start,start+d*direction,outer);
      real[][] ints=intersections(start,end,outer);
      assert(ints.length != 0);
      real endtime=ints[0][1]; // endtime is time on outer
      end = point(outer,endtime);
      // find first intersection of end--start with any inner curve
      real starttime=inner_i; // starttime is time on inners[curveIndex]
      real earliestTime=1;
      for(int j=0; j < inners.length; ++j) {
        real[][] ints=intersections(end,start,inners[j]);
        
        if(ints.length > 0 && ints[0][0] < earliestTime) {
          earliestTime=ints[0][0]; // time on end--start
          starttime=ints[0][1]; // time on inner curve
          curveIndex=j;
        }
      }
      start=point(inners[curveIndex],starttime);
      
      
      bool found_forward = false;
      real timeoffset_forward = 2;
      path portion_forward;
      path[] allCurves = {outer};
      allCurves.append(inners);

      while(!found_forward && timeoffset_forward > fuzz) {
        timeoffset_forward /= 2;
        if(countIntersections(allCurves,start,
                              point(outer,endtime+timeoffset_forward)) == 2)
          {
            portion_forward = subpath(outer,endtime,endtime+timeoffset_forward)--start--cycle;
            
            found_forward=true;
            // check if an inner curve is inside the portion
            for(int k = 0; found_forward && k < inners.length; ++k)
              {
                if(k!=curveIndex &&
                   inside(portion_forward,point(inners[k],0),zerowinding))
                  found_forward = false;
              }
          }
      }
      
      bool found_backward = false;
      real timeoffset_backward = -2;
      path portion_backward;
      while(!found_backward && timeoffset_backward < -fuzz) {
        timeoffset_backward /= 2;
        if(countIntersections(allCurves,start,
                              point(outer,endtime+timeoffset_backward))==2)
          {
            portion_backward = subpath(outer,endtime+timeoffset_backward,endtime)--start--cycle;
            found_backward = true;
            // check if an inner curve is inside the portion
            for(int k = 0; found_backward && k < inners.length; ++k)
              {
                if(k!=curveIndex &&
                   inside(portion_backward,point(inners[k],0),zerowinding))
                  found_backward = false;
              }
          }
      }
      assert(found_forward || found_backward);
      real timeoffset;
      path portion;
      if(found_forward && !found_backward)
      {
        timeoffset = timeoffset_forward;
        portion = portion_forward;
      }
      else if(found_backward && !found_forward)
      {
        timeoffset = timeoffset_backward;
        portion = portion_backward;
      }
      else // assert handles case of neither found
      {
        if(timeoffset_forward > -timeoffset_backward)
        {
          timeoffset = timeoffset_forward;
          portion = portion_forward;
        }
        else
        {
          timeoffset = timeoffset_backward;
          portion = portion_backward;
         }
      }
      
      endtime=min(endtime,endtime+timeoffset);
      // or go from timeoffset+timeoffset_backward to timeoffset+timeoffset_forward?
      timeoffset=abs(timeoffset);

      // depends on the curves having opposite orientations
      path remainder=section(outer,endtime+timeoffset,endtime)
        --uncycle(inners[curveIndex],
                  starttime)--cycle;
      inners.delete(curveIndex);
      outer = remainder;
      patch.append(portion);
    }
    result.append(outer);
  }
}

bool checkSegment(path g, pair p, pair q)
{
  pair mid=0.5*(p+q);
  return intersections(p,q,g).length == 2 &&
    inside(g,mid,zerowinding) && intersections(g,mid).length == 0;
}

path subdivide(path p)
{
  path q;
  int l=length(p);
  for(int i=0; i < l; ++i)
    q=q&subpath(p,i,i+0.5)&subpath(p,i+0.5,i+1);
  return cyclic(p) ? q&cycle : q;
}

path[] bezulate(path[] p)
{
  if(p.length == 1 && length(p[0]) <= 4) return p;
  path[] patch;
  path[] result;
  connect(p,result,patch);
  for(int i=0; i < result.length; ++i) {
    path p=result[i];
    int refinements=0;
    if(size(p) <= 1) return p;
    if(!cyclic(p))
      abort("path must be cyclic and nonselfintersecting.");
    p=removeDuplicates(p);
    if(length(p) > 4) {
      static real SIZE_STEPS=10;
      static real factor=1.05/SIZE_STEPS;
      for(int k=1; k <= SIZE_STEPS; ++k) {
        real L=factor*k*abs(max(p)-min(p));
        for(int i=0; length(p) > 4 && i < length(p); ++i) {
          bool found=false;
          pair start=point(p,i);
          //look for quadrilaterals and triangles with one line, 4 | 3 curves
          for(int desiredSides=4; !found && desiredSides >= 3;
              --desiredSides) {
            if(desiredSides == 3 && length(p) <= 3)
              break;
            pair end;
            int endi=i+desiredSides-1;
            end=point(p,endi);
            found=checkSegment(p,start,end) && abs(end-start) < L;
            if(found) {
              path p1=subpath(p,endi,i+length(p))--cycle;
              patch.append(subpath(p,i,endi)--cycle);
              p=removeDuplicates(p1);
              i=-1; // increment will make i be 0
            }
          }
          if(!found && k == SIZE_STEPS && length(p) > 4 && i == length(p)-1) {
            // avoid infinite recursion
            ++refinements;
            if(refinements > maxrefinements) {
              warning("subdivisions","too many subdivisions",position=true);
            } else {
              p=subdivide(p);
              i=-1;
            }
          }
        }
      }
    }
    if(length(p) <= 4)
      patch.append(p);
  }
  return patch;
}