import bezulate;

int nslice=12;
real camerafactor=1.2;

private real Fuzz=10.0*realEpsilon;
private real nineth=1/9;

struct patch {
  triple[][] P=new triple[4][4];
  triple[] normals; // Optionally specify 4 normal vectors at the corners.
  pen[] colors;     // Optionally specify 4 corner colors.
  bool straight;    // Patch is based on a piecewise straight external path.
  bool3 planar;     // Patch is planar.

  path3 external() {
    return
      P[0][0]..controls P[0][1] and P[0][2]..
      P[0][3]..controls P[1][3] and P[2][3]..
      P[3][3]..controls P[3][2] and P[3][1]..
      P[3][0]..controls P[2][0] and P[1][0]..cycle;
  }

  triple[] internal() {
    return new triple[] {P[1][1],P[1][2],P[2][2],P[2][1]};
  }

  triple cornermean() {
    return 0.25*(P[0][0]+P[0][3]+P[3][3]+P[3][0]);
  }

  triple[] corners() {return new triple[] {P[0][0],P[0][3],P[3][3],P[3][0]};}

  real[] map(real f(triple)) {
    return new real[] {f(P[0][0]),f(P[0][3]),f(P[3][3]),f(P[3][0])};
  }

  triple Bu(int j, real u) {return bezier(P[0][j],P[1][j],P[2][j],P[3][j],u);}
  triple BuP(int j, real u) {return bezierP(P[0][j],P[1][j],P[2][j],P[3][j],u);}
  triple BuPP(int j, real u) {
    return bezierPP(P[0][j],P[1][j],P[2][j],P[3][j],u);
  }
  triple BuPPP(int j) {return bezierPPP(P[0][j],P[1][j],P[2][j],P[3][j]);}

  path3 uequals(real u) {
    return straight ? Bu(0,u)--Bu(3,u) :
      Bu(0,u)..controls Bu(1,u) and Bu(2,u)..Bu(3,u);
  }

  triple Bv(int i, real v) {return bezier(P[i][0],P[i][1],P[i][2],P[i][3],v);}
  triple BvP(int i, real v) {return bezierP(P[i][0],P[i][1],P[i][2],P[i][3],v);}
  triple BvPP(int i, real v) {
    return bezierPP(P[i][0],P[i][1],P[i][2],P[i][3],v);
  }
  triple BvPPP(int i) {return bezierPPP(P[i][0],P[i][1],P[i][2],P[i][3]);}

  path3 vequals(real v) {
    return straight ? Bv(0,v)--Bv(3,v) :
      Bv(0,v)..controls Bv(1,v) and Bv(2,v)..Bv(3,v);
  }

  triple point(real u, real v) {        
    return bezier(Bu(0,u),Bu(1,u),Bu(2,u),Bu(3,u),v);
  }

  // compute normal vectors for degenerate cases
  private triple normal0(real u, real v, real epsilon) {
    triple n=0.5*(cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u),
                        bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v))+
                  cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u),   
                        bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v)));
    return abs(n) > epsilon ? n :
      0.25*cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u),   
                 bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v))+
      1/6*(cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u),
                 bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v))+
           cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u),   
                 bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v)))+
      1/12*(cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u),
                  bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v))+
            cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u),   
                  bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v)))+
      1/36*cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u),   
                 bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v));
  }

  static real fuzz=1000*realEpsilon;

  triple normal(real u, real v) {
    triple n=cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u),
                   bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v));
    real epsilon=fuzz*change2(P);
    return (abs(n) > epsilon) ? n : normal0(u,v,epsilon);
  }
  
  triple normal00() {
    triple n=9*cross(P[0][1]-P[0][0],P[1][0]-P[0][0]);
    real epsilon=fuzz*change2(P);
    return abs(n) > epsilon ? n : normal0(0,0,epsilon);
  }

  triple normal01() {
    triple n=9*cross(P[0][3]-P[0][2],P[1][3]-P[0][3]);
    real epsilon=fuzz*change2(P);
    return abs(n) > epsilon ? n : normal0(0,1,epsilon);
  }

  triple normal11() {
    triple n=9*cross(P[3][3]-P[3][2],P[3][3]-P[2][3]);
    real epsilon=fuzz*change2(P);
    return abs(n) > epsilon ? n : normal0(1,1,epsilon);
  }

  triple normal10() {
    triple n=9*cross(P[3][1]-P[3][0],P[3][0]-P[2][0]);
    real epsilon=fuzz*change2(P);
    return abs(n) > epsilon ? n : normal0(1,0,epsilon);
  }

  pen[] colors(material m, light light=currentlight) {
    bool nocolors=colors.length == 0;
    if(normals.length > 0)
      return new pen[] {light.color(normals[0],nocolors ? m : colors[0]),
          light.color(normals[1],nocolors ? m : colors[1]),
          light.color(normals[2],nocolors ? m : colors[2]),
          light.color(normals[3],nocolors ? m : colors[3])};
    if(planar) {
      triple normal=normal(0.5,0.5);
      return new pen[] {light.color(normal,nocolors ? m : colors[0]),
          light.color(normal,nocolors ? m : colors[1]),
          light.color(normal,nocolors ? m : colors[2]),
          light.color(normal,nocolors ? m : colors[3])};
    }
    return new pen[] {light.color(normal00(),nocolors ? m : colors[0]),
        light.color(normal01(),nocolors ? m : colors[1]),
        light.color(normal11(),nocolors ? m : colors[2]),
        light.color(normal10(),nocolors ? m : colors[3])};
  }
  
  triple min3,max3;
  bool havemin3,havemax3;

  void init() {
    havemin3=false;
    havemax3=false;
  }

  triple min(triple bound=P[0][0]) {
    if(havemin3) return minbound(min3,bound);
    havemin3=true;
    return min3=minbezier(P,bound);
  }

  triple max(triple bound=P[0][0]) {
    if(havemax3) return maxbound(max3,bound);
    havemax3=true;
    return max3=maxbezier(P,bound);
  }

  triple center() {
    return 0.5*(this.min()+this.max());
  }

  pair min(projection P, pair bound=project(this.P[0][0],P.t)) {
    return minbezier(this.P,P.t,bound);
  }

  pair max(projection P, pair bound=project(this.P[0][0],P.t)) {
    return maxbezier(this.P,P.t,bound);
  }

  void operator init(triple[][] P, triple[] normals=new triple[],
                     pen[] colors=new pen[], bool straight=false,
                     bool3 planar=default) {
    init();
    this.P=copy(P);
    if(normals.length != 0)
      this.normals=copy(normals);
    if(colors.length != 0)
      this.colors=copy(colors);
    this.planar=planar;
    this.straight=straight;
  }

  void operator init(pair[][] P, triple plane(pair)=XYplane,
                     bool straight=false) {
    triple[][] Q=new triple[4][];
    for(int i=0; i < 4; ++i) {
      pair[] Pi=P[i];
      Q[i]=sequence(new triple(int j) {return plane(Pi[j]);},4);
    }
    operator init(Q,straight);
    planar=true;
  }

  void operator init(patch s) {
    operator init(s.P,s.normals,s.colors,s.straight);
  }
  
  // A constructor for a convex cyclic path3 of length <= 4 with optional
  // arrays of 4 internal points, corner normals, and pens.
  void operator init(path3 external, triple[] internal=new triple[],
                     triple[] normals=new triple[], pen[] colors=new pen[],
                     bool3 planar=default) {
    init();

    if(internal.length == 0 && planar == default)
      this.planar=normal(external) != O;
    else this.planar=planar;

    int L=length(external);
    if(L > 4 || !cyclic(external))
      abort("cyclic path3 of length <= 4 expected");
    if(L == 1) {
      external=external--cycle--cycle--cycle;
      if(colors.length > 0) colors.append(array(3,colors[0]));
      if(normals.length > 0) normals.append(array(3,normals[0]));
    } else if(L == 2) {
      external=external--cycle--cycle;
      if(colors.length > 0) colors.append(array(2,colors[0]));
      if(normals.length > 0) normals.append(array(2,normals[0]));
    } else if(L == 3) {
      external=external--cycle;
      if(colors.length > 0) colors.push(colors[0]);
      if(normals.length > 0) normals.push(normals[0]);
    }
    if(normals.length != 0)
      this.normals=copy(normals);
    if(colors.length != 0)
      this.colors=copy(colors);

    if(internal.length == 0) {
      straight=piecewisestraight(external);
      internal=new triple[4];
      for(int j=0; j < 4; ++j)
        internal[j]=nineth*(-4*point(external,j)
                            +6*(precontrol(external,j)+postcontrol(external,j))
                            -2*(point(external,j-1)+point(external,j+1))
                            +3*(precontrol(external,j-1)+
                                postcontrol(external,j+1))
                            -point(external,j+2));
    } else straight=false;

    P=new triple[][] {
      {point(external,0),postcontrol(external,0),precontrol(external,1),
       point(external,1)},
      {precontrol(external,0),internal[0],internal[1],postcontrol(external,1)},
      {postcontrol(external,3),internal[3],internal[2],precontrol(external,2)},
      {point(external,3),precontrol(external,3),postcontrol(external,2),
       point(external,2)}
    };
  }

  // A constructor for a convex quadrilateral.
  void operator init(triple[] external, triple[] internal=new triple[],
                     triple[] normals=new triple[], pen[] colors=new pen[],
                     bool3 planar=default) {
    init();

    if(internal.length == 0 && planar == default)
      this.planar=normal(external) != O;
    else this.planar=planar;

    if(normals.length != 0)
      this.normals=copy(normals);
    if(colors.length != 0)
      this.colors=copy(colors);

    if(internal.length == 0) {
      internal=new triple[4];
      for(int j=0; j < 4; ++j)
        internal[j]=nineth*(4*external[j]+2*external[(j+1)%4]+
                            external[(j+2)%4]+2*external[(j+3)%4]);
    }

    straight=true;

    triple delta[]=new triple[4];
    for(int j=0; j < 4; ++j)
      delta[j]=(external[(j+1)% 4]-external[j])/3;

    P=new triple[][] {
      {external[0],external[0]+delta[0],external[1]-delta[0],external[1]},
      {external[0]-delta[3],internal[0],internal[1],external[1]+delta[1]},
      {external[3]+delta[3],internal[3],internal[2],external[2]-delta[1]},
      {external[3],external[3]-delta[2],external[2]+delta[2],external[2]}
    };
  }
}

patch operator * (transform3 t, patch s)
{ 
  patch S;
  for(int i=0; i < 4; ++i) { 
    triple[] si=s.P[i];
    triple[] Si=S.P[i];
    for(int j=0; j < 4; ++j)
      Si[j]=t*si[j]; 
  }
  
  transform3 t0=shiftless(t);
  for(int i=0; i < s.normals.length; ++i)
    S.normals[i]=t0*s.normals[i];

  S.colors=copy(s.colors);
  S.planar=s.planar;
  S.straight=s.straight;
  return S;
}
 
patch reverse(patch s) 
{
  patch S;
  S.P=transpose(s.P);
  if(s.normals.length > 0) 
    S.normals=
      new triple[] {s.normals[0],s.normals[3],s.normals[2],s.normals[1]};
  if(s.colors.length > 0) 
    S.colors=new pen[] {s.colors[0],s.colors[3],s.colors[2],s.colors[1]};
  S.planar=s.planar;
  S.straight=s.straight;
  return S;
}

struct surface {
  patch[] s;
  
  bool empty() {
    return s.length == 0;
  }

  void operator init(int n) {
    s=new patch[n];
  }

  void operator init(... patch[] s) {
    this.s=s;
  }

  void operator init(surface s) {
    this.s=new patch[s.s.length];
    for(int i=0; i < s.s.length; ++i)
      this.s[i]=patch(s.s[i]);
  }

  void operator init(triple[][][] P, triple[][] normals=new triple[][],
                     pen[][] colors=new pen[][], bool3 planar=default) {
    s=sequence(new patch(int i) {
        return patch(P[i],normals.length == 0 ? new triple[] : normals[i],
                     colors.length == 0 ? new pen[] : colors[i],planar);
      },P.length);
  }

  void colors(pen[][] palette) {
    for(int i=0; i < s.length; ++i) {
      pen[] palettei=palette[i];
      s[i].colors=new pen[] {palettei[0],palettei[1],palettei[2],palettei[3]};
    }
  }

  triple[][] corners() {
    triple[][] a=new triple[s.length][];
    for(int i=0; i < s.length; ++i)
      a[i]=s[i].corners();
    return a;
  }

  real[][] map(real f(triple)) {
    real[][] a=new real[s.length][];
    for(int i=0; i < s.length; ++i)
      a[i]=s[i].map(f);
    return a;
  }

  triple[] cornermean() {
    return sequence(new triple(int i) {return s[i].cornermean();},s.length);
  }

  // A constructor for a possibly nonconvex cyclic path in a given plane.
  void operator init (path p, triple plane(pair)=XYplane,
                      bool checkboundary=true) {
    if(!cyclic(p))
      abort("cyclic path expected");

    int L=length(p);

    if(L > 4) {
      for(path g : bezulate(p))
        s.append(surface(g,plane,checkboundary).s);
      return;
    }
        
    pair[][] P(path p) {
      if(L == 1)
        p=p--cycle--cycle--cycle;
      else if(L == 2)
        p=p--cycle--cycle;
      else if(L == 3)
        p=p--cycle;
 
      pair[] internal=new pair[4];
      for(int j=0; j < 4; ++j) {
        internal[j]=nineth*(-4*point(p,j)
                            +6*(precontrol(p,j)+postcontrol(p,j))
                            -2*(point(p,j-1)+point(p,j+1))
                            +3*(precontrol(p,j-1)+postcontrol(p,j+1))
                            -point(p,j+2));
      }
    
      return new pair[][] {
        {point(p,0),postcontrol(p,0),precontrol(p,1),point(p,1)},
          {precontrol(p,0),internal[0],internal[1],postcontrol(p,1)},
            {postcontrol(p,3),internal[3],internal[2],precontrol(p,2)},
              {point(p,3),precontrol(p,3),postcontrol(p,2),point(p,2)}
      };
    }

    bool straight=piecewisestraight(p);
    if(L <= 3 && straight) {
      s=new patch[] {patch(P(p),plane,straight)};
      return;
    }
    
    // Split p along the angle bisector at t.
    bool split(path p, real t) {
      pair dir=dir(p,t);
      if(dir != 0) {
        path g=subpath(p,t,t+length(p));
        int L=length(g);
        pair z=point(g,0);
        real[] T=intersections(g,z,z+I*dir);
        for(int i=0; i < T.length; ++i) {
          real cut=T[i];
          if(cut > sqrtEpsilon && cut < L-sqrtEpsilon) {
            pair w=point(g,cut);
            if(!inside(p,0.5*(z+w),zerowinding)) continue;
            pair delta=sqrtEpsilon*(w-z);
            if(intersections(g,z-delta--w+delta).length != 2) continue;
            s=surface(subpath(g,0,cut)--cycle,plane,checkboundary).s;
            s.append(surface(subpath(g,cut,L)--cycle,plane,checkboundary).s);
            return true;
          }
        }
      }
      return false;
    }

    // Ensure that all interior angles are less than 180 degrees.
    real fuzz=1e-4;
    int sign=sgn(windingnumber(p,inside(p,zerowinding)));
    for(int i=0; i < L; ++i) {
      if(sign*(conj(dir(p,i,-1))*dir(p,i,1)).y < -fuzz) {
        if(split(p,i)) return;
      }
    }

    pair[][] P=P(p);

    if(straight) {
      s=new patch[] {patch(P,plane,straight)};
      return;
    }
    
    // Check for degeneracy.
    pair[][] U=new pair[3][4];
    pair[][] V=new pair[4][3];

    for(int i=0; i < 3; ++i) {
      for(int j=0; j < 4; ++j)
        U[i][j]=P[i+1][j]-P[i][j];
    }
      
    for(int i=0; i < 4; ++i) {
      for(int j=0; j < 3; ++j)
        V[i][j]=P[i][j+1]-P[i][j];
    }
      
    int[] choose2={1,2,1};
    int[] choose3={1,3,3,1};

    real T[][]=new real[6][6];
    for(int p=0; p < 6; ++p) {
      int kstart=max(p-2,0);
      int kstop=min(p,3);
      real[] Tp=T[p];
      for(int q=0; q < 6; ++q) {
        real Tpq;
        int jstop=min(q,3);
        int jstart=max(q-2,0);
        for(int k=kstart; k <= kstop; ++k) {
          int choose3k=choose3[k];
          for(int j=jstart; j <= jstop; ++j) {
            int i=p-k;
            int l=q-j;
            Tpq += (conj(U[i][j])*V[k][l]).y*
              choose2[i]*choose3k*choose3[j]*choose2[l];
          }
        }
        Tp[q]=Tpq;
      }
    }

    bool3 aligned=default;
    bool degenerate=false;

    for(int p=0; p < 6; ++p) {
      for(int q=0; q < 6; ++q) {
        if(aligned == default) {
          if(T[p][q] < -sqrtEpsilon) aligned=true;
          if(T[p][q] > sqrtEpsilon) aligned=false;
        } else {
          if((T[p][q] < -sqrtEpsilon && aligned == false) ||
             (T[p][q] > sqrtEpsilon && aligned == true)) degenerate=true;
        }
      }
    }

    if(!degenerate) {
      if(aligned == (sign >= 0))
        s=new patch[] {patch(P,plane)};
      return;
    }

    if(checkboundary) {
      // Polynomial coefficients of (B_i'' B_j + B_i' B_j')/3.
      static real[][][] fpv0={
        {{5, -20, 30, -20, 5},
         {-3, 24, -54, 48, -15},
         {0, -6, 27, -36, 15},
         {0, 0, -3, 8, -5}},
        {{-7, 36, -66, 52, -15},
         {3, -36, 108, -120, 45},
         {0, 6, -45, 84, -45},
         {0, 0, 3, -16, 15}},
        {{2, -18, 45, -44, 15},
         {0, 12, -63, 96, -45},
         {0, 0, 18, -60, 45},
         {0, 0, 0, 8, -15}},
        {{0, 2, -9, 12, -5},
         {0, 0, 9, -24, 15},
         {0, 0, 0, 12, -15},
         {0, 0, 0, 0, 5}}
      };

      // Compute one-ninth of the derivative of the Jacobian along the boundary.
      real[][] c=array(4,array(5,0.0));
      for(int i=0; i < 4; ++i) {
        real[][] fpv0i=fpv0[i];
        for(int j=0; j < 4; ++j) {
          real[] w=fpv0i[j];
          c[0] += w*(conj(P[0][j]-P[1][j])*P[0][i]).y;   // u=0
          c[1] += w*(conj(P[i][3])*(P[j][3]-P[j][2])).y; // v=1
          c[2] += w*(conj(P[3][j]-P[2][j])*P[3][i]).y;   // u=1
          c[3] += w*(conj(P[i][0])*(P[j][1]-P[j][0])).y; // v=0
        }
      }
    
      pair BuP(int j, real u) {
        return bezierP(P[0][j],P[1][j],P[2][j],P[3][j],u);
      }
      pair BvP(int i, real v) {
        return bezierP(P[i][0],P[i][1],P[i][2],P[i][3],v);
      }
      real normal(real u, real v) {
        return (conj(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u))*
                bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v)).y;
      }

      // Use Rolle's theorem to check for degeneracy on the boundary.
      real M=0;
      real cut;
      for(int i=0; i < 4; ++i) {
        if(!straight(p,i)) {
          real[] ci=c[i];
          pair[] R=quarticroots(ci[4],ci[3],ci[2],ci[1],ci[0]);
          for(pair r : R) {
            if(fabs(r.y) < sqrtEpsilon) {
              real t=r.x;
              if(0 <= t && t <= 1) {
                real[] U={0,t,1,t};
                real[] V={t,1,t,0};
                real[] T={t,t,1-t,1-t};
                real N=sign*normal(U[i],V[i]);
                if(N < M) {
                  M=N; cut=i+T[i];
                }
              }
            }
          }
        }
      }

      // Split at the worst boundary degeneracy.
      if(M < 0 && split(p,cut)) return;
    }
    
    // Split arbitrarily to resolve any remaining (internal) degeneracy.
    checkboundary=false;
    for(int i=0; i < L; ++i)
      if(!straight(p,i) && split(p,i+0.5)) return;

    while(true)
      for(int i=0; i < L; ++i)
        if(!straight(p,i) && split(p,i+unitrand())) return;
  }

  void operator init(explicit path[] g, triple plane(pair)=XYplane) {
    for(path p : bezulate(g))
      s.append(surface(p,plane).s);
  }

  // A general surface constructor for both planar and nonplanar 3D paths.
  void construct(path3 external, triple[] internal=new triple[],
                 triple[] normals=new triple[], pen[] colors=new pen[],
                 bool3 planar=default) {
    int L=length(external);
    if(!cyclic(external)) abort("cyclic path expected");

    if(L <= 3 && piecewisestraight(external)) {
      s.push(patch(external,internal,normals,colors,planar=true));
      return;
    }

    // Construct a surface from a possibly nonconvex planar cyclic path3.
    if(planar != false && internal.length == 0 && normals.length == 0 &&
       colors.length == 0) {
      triple n=normal(external);
      if(n != O) {
        transform3 T=align(n);
        external=transpose(T)*external;
        T *= shift(0,0,point(external,0).z);
        for(patch p : surface(path(external)).s)
          s.push(T*p);
        return;
      }
    }
    
    if(L <= 4 || internal.length > 0) {
      s.push(patch(external,internal,normals,colors,planar));
      return;
    }
      
    // Path is not planar; split into patches.
    real factor=1/L;
    pen[] p;
    triple[] n;
    bool nocolors=colors.length == 0;
    bool nonormals=normals.length == 0;
    triple center;
    for(int i=0; i < L; ++i)
      center += point(external,i);
    center *= factor;
    if(!nocolors)
      p=new pen[] {mean(colors)};
    if(!nonormals)
      n=new triple[] {factor*sum(normals)};
    // Use triangles for nonplanar surfaces.
    int step=normal(external) == O ? 1 : 2;
    int i=0;
    int end;
    while((end=i+step) < L) {
      s.push(patch(subpath(external,i,end)--center--cycle,
                   nonormals ? n : concat(normals[i:end+1],n),
                   nocolors ? p : concat(colors[i:end+1],p),planar));
      i=end;
    }
    s.push(patch(subpath(external,i,L)--center--cycle,
                 nonormals ? n : concat(normals[i:],normals[0:1],n),
                 nocolors ? p : concat(colors[i:],colors[0:1],p),planar));
  }

  void operator init(path3 external, triple[] internal=new triple[],
                     triple[] normals=new triple[], pen[] colors=new pen[],
                     bool3 planar=default) {
    s=new patch[];
    construct(external,internal,normals,colors,planar);
  }

  void operator init(explicit path3[] external,
                     triple[][] internal=new triple[][],
                     triple[][] normals=new triple[][],
                     pen[][] colors=new pen[][], bool3 planar=default) {
    s=new patch[];
    if(planar == true) {// Assume all path3 elements share a common normal.
      if(external.length != 0) {
        triple n=normal(external[0]);
        if(n != O) {
          transform3 T=align(n);
          external=transpose(T)*external;
          T *= shift(0,0,point(external[0],0).z);
          path[] g=sequence(new path(int i) {return path(external[i]);},
                            external.length);
          for(patch p : surface(g).s)
            s.push(T*p);
          return;
        }
      }
    }

    for(int i=0; i < external.length; ++i)
      construct(external[i],
                internal.length == 0 ? new triple[] : internal[i],
                normals.length == 0 ? new triple[] : normals[i],
                colors.length == 0 ? new pen[] : colors[i],planar);
  }

  void push(path3 external, triple[] internal=new triple[],
            triple[] normals=new triple[] ,pen[] colors=new pen[],
            bool3 planar=default) {
    s.push(patch(external,internal,normals,colors,planar));
  }

  // Construct the surface of rotation generated by rotating g
  // from angle1 to angle2 sampled n times about the line c--c+axis.
  // An optional surface pen color(int i, real j) may be specified
  // to override the color at vertex(i,j).
  void operator init(triple c, path3 g, triple axis, int n=nslice,
                     real angle1=0, real angle2= 360,
                     pen color(int i, real j)=null) {
    axis=unit(axis);
    real w=(angle2-angle1)/n;
    int L=length(g);
    s=new patch[L*n];
    int m=-1;
    transform3[] T=new transform3[n+1];
    transform3 t=rotate(w,c,c+axis);
    T[0]=rotate(angle1,c,c+axis);
    for(int k=1; k <= n; ++k)
      T[k]=T[k-1]*t;

    for(int i=0; i < L; ++i) {
      path3 h=subpath(g,i,i+1);
      path3 r=reverse(h);
      path3 H=shift(-c)*h;
      real M=0;
      triple perp;
      void test(real[] t) {
        for(int i=0; i < 3; ++i) {
          triple v=point(H,t[i]);
          triple V=v-dot(v,axis)*axis;
          real a=abs(V);
          if(a > M) {M=a; perp=V;}
        }
      }
      test(maxtimes(H));
      test(mintimes(H));
      
      perp=unit(perp);
      triple normal=unit(cross(axis,perp));
      triple dir(real j) {return Cos(j)*normal-Sin(j)*perp;}
      real j=angle1;
      transform3 Tk=T[0];
      triple dirj=dir(j);
      for(int k=0; k < n; ++k, j += w) {
        transform3 Tp=T[k+1];
        triple dirp=dir(j+w);
        path3 G=Tk*h{dirj}..{dirp}Tp*r{-dirp}..{-dirj}cycle;
        Tk=Tp;
        dirj=dirp;
        s[++m]=color == null ? patch(G) :
          patch(G,new pen[] {color(i,j),color(i+1,j),color(i+1,j+w),
                             color(i,j+w)});
      }
    }
  }

  void push(patch s) {
    this.s.push(s);
  }

  void append(surface s) {
    this.s.append(s.s);
  }

  void operator init(... surface[] s) {
    for(surface S : s)
      this.s.append(S.s);
  }
}

surface operator * (transform3 t, surface s)
{ 
  surface S;
  S.s=new patch[s.s.length];
  for(int i=0; i < s.s.length; ++i)
    S.s[i]=t*s.s[i];
  return S;
}

private string nullsurface="null surface";

triple min(surface s)
{
  if(s.s.length == 0)
    abort(nullsurface);
  triple bound=s.s[0].min();
  for(int i=1; i < s.s.length; ++i)
    bound=s.s[i].min(bound);
  return bound;
}
  
triple max(surface s)
{
  if(s.s.length == 0)
    abort(nullsurface);
  triple bound=s.s[0].max();
  for(int i=1; i < s.s.length; ++i)
    bound=s.s[i].max(bound);
  return bound;
}

pair min(surface s, projection P)
{
  if(s.s.length == 0)
    abort(nullsurface);
  pair bound=s.s[0].min(P);
  for(int i=1; i < s.s.length; ++i)
    bound=s.s[i].min(P,bound);
  return bound;
}
  
pair max(surface s, projection P)
{
  if(s.s.length == 0)
    abort(nullsurface);
  pair bound=s.s[0].max(P);
  for(int i=1; i < s.s.length; ++i)
    bound=s.s[i].max(P,bound);
  return bound;
}

private triple[] split(triple z0, triple c0, triple c1, triple z1, real t=0.5)
{
  triple m0=interp(z0,c0,t);
  triple m1=interp(c0,c1,t);
  triple m2=interp(c1,z1,t);
  triple m3=interp(m0,m1,t);
  triple m4=interp(m1,m2,t);
  triple m5=interp(m3,m4,t);

  return new triple[] {m0,m3,m5,m4,m2};
}

// Return the control points for a subpatch of P on [u,1] x [v,1].
triple[][] subpatchbegin(triple[][] P, real u, real v)
{
  triple[] P0=P[0];
  triple[] P1=P[1];
  triple[] P2=P[2];
  triple[] P3=P[3];

  triple[] c0=split(P0[0],P0[1],P0[2],P0[3],v);
  triple[] c1=split(P1[0],P1[1],P1[2],P1[3],v);
  triple[] c2=split(P2[0],P2[1],P2[2],P2[3],v);
  triple[] c3=split(P3[0],P3[1],P3[2],P3[3],v);

  u=1.0-u;
  
  triple[] c7=split(c3[2],c2[2],c1[2],c0[2],u);
  triple[] c8=split(c3[3],c2[3],c1[3],c0[3],u);
  triple[] c9=split(c3[4],c2[4],c1[4],c0[4],u);
  triple[] c10=split(P3[3],P2[3],P1[3],P0[3],u);

  return new triple[][] {{c7[2],c8[2],c9[2],c10[2]},
      {c7[1],c8[1],c9[1],c10[1]},
        {c7[0],c8[0],c9[0],c10[0]},
          {c3[2],c3[3],c3[4],P3[3]}};
}

// Return the control points for a subpatch of P on [0,u] x [0,v].
triple[][] subpatchend(triple[][] P, real u, real v)
{
  triple[] P0=P[0];
  triple[] P1=P[1];
  triple[] P2=P[2];
  triple[] P3=P[3];

  triple[] c0=split(P0[0],P0[1],P0[2],P0[3],v);
  triple[] c1=split(P1[0],P1[1],P1[2],P1[3],v);
  triple[] c2=split(P2[0],P2[1],P2[2],P2[3],v);
  triple[] c3=split(P3[0],P3[1],P3[2],P3[3],v);

  u=1.0-u;
  
  triple[] c4=split(P3[0],P2[0],P1[0],P0[0],u);
  triple[] c5=split(c3[0],c2[0],c1[0],c0[0],u);
  triple[] c6=split(c3[1],c2[1],c1[1],c0[1],u);
  triple[] c7=split(c3[2],c2[2],c1[2],c0[2],u);

  return new triple[][] {
    {P0[0],c0[0],c0[1],c0[2]},
      {c4[4],c5[4],c6[4],c7[4]},
        {c4[3],c5[3],c6[3],c7[3]},
          {c4[2],c5[2],c6[2],c7[2]}};
}

patch subpatch(patch s, real ua, real va, real ub, real vb)
{
  assert(ua >= 0 && va >= 0 && ub <= 1 && vb <= 1 && ua < ub && va < vb);
  return patch(subpatchbegin(subpatchend(s.P,ub,vb),ua/ub,va/vb),
               s.straight,s.planar);
}

// return an array containing all intersection times of path p and patch s.
real[][] intersections(path3 p, patch s, real fuzz=-1)
{
  return sort(intersections(p,s.P,fuzz));
}

// return an array containing all intersection times of path p and surface s.
real[][] intersections(path3 p, surface s, real fuzz=-1)
{
  real[][] T;
  if(length(p) < 0) return T;
  for(int i=0; i < s.s.length; ++i)
    for(real[] s: intersections(p,s.s[i].P,fuzz))
      T.push(s);

  static real fuzzFactor=10.0;
  static real Fuzz=1000.0*realEpsilon;
  real fuzz=max(fuzzFactor*fuzz,Fuzz)*abs(max(s)-min(s));
  
  // Remove intrapatch duplicate points.
  for(int i=0; i < T.length; ++i) {
    triple v=point(p,T[i][0]);
    for(int j=i+1; j < T.length;) {
      if(abs(v-point(p,T[j][0])) < fuzz)
        T.delete(j);
      else ++j;
    }
  }
  return sort(T);
}

// return an array containing all intersection points of path p and surface s.
triple[] intersectionpoints(path3 p, patch s, real fuzz=-1)
{
  real[][] t=intersections(p,s,fuzz);
  return sequence(new triple(int i) {return point(p,t[i][0]);},t.length);
}

// return an array containing all intersection points of path p and surface s.
triple[] intersectionpoints(path3 p, surface s, real fuzz=-1)
{
  real[][] t=intersections(p,s,fuzz);
  return sequence(new triple(int i) {return point(p,t[i][0]);},t.length);
}

triple point(patch s, real u, real v)
{
  return s.point(u,v);
}

void draw3D(frame f, patch s, material m, light light=currentlight)
{
  if(s.colors.length > 0)
    m=mean(s.colors);
  bool lighton=light.on();
  if(!lighton && !invisible((pen) m))
    m=emissive(m);
  real granularity=m.granularity >= 0 ? m.granularity : defaultgranularity;
  draw(f,s.P,s.straight,m.p,m.opacity,m.shininess,granularity,
       s.planar ? s.normal(0.5,0.5) : O,lighton,s.colors);
}

void tensorshade(transform t=identity(), frame f, patch s,
                 material m, light light=currentlight, projection P)
{
  tensorshade(f,box(t*s.min(P),t*s.max(P)),m.diffuse(),
              s.colors(m,light),t*project(s.external(),P,1),
              t*project(s.internal(),P));
}

restricted pen[] nullpens={nullpen};
nullpens.cyclic(true);

void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
          material[] surfacepen, pen[] meshpen=nullpens,
          light light=currentlight, light meshlight=light,
          projection P=currentprojection)
{
  if(is3D()) {
    for(int i=0; i < s.s.length; ++i)
      draw3D(f,s.s[i],surfacepen[i],light);
    pen modifiers=thin()+squarecap;
    for(int k=0; k < s.s.length; ++k) {
      pen meshpen=meshpen[k];
      if(!invisible(meshpen)) {
        meshpen=modifiers+meshpen;
        real step=nu == 0 ? 0 : 1/nu;
        for(int i=0; i <= nu; ++i)
          draw(f,s.s[k].uequals(i*step),meshpen,meshlight);
        step=nv == 0 ? 0 : 1/nv;
        for(int j=0; j <= nv; ++j)
          draw(f,s.s[k].vequals(j*step),meshpen,meshlight);
      }
    }
  } else {
    begingroup(f);
    // Sort patches by mean distance from camera
    triple camera=P.camera;
    if(P.infinity) {
      triple m=min(s);
      triple M=max(s);
      camera=P.target+camerafactor*(abs(M-m)+abs(m-P.target))*unit(P.vector());
    }

    real[][] depth=new real[s.s.length][];
    for(int i=0; i < depth.length; ++i)
      depth[i]=new real[] {abs(camera-s.s[i].cornermean()),i};

    depth=sort(depth);

    light.T=shiftless(P.modelview());

    // Draw from farthest to nearest
    while(depth.length > 0) {
      real[] a=depth.pop();
      int i=round(a[1]);
      tensorshade(t,f,s.s[i],surfacepen[i],light,P);
      pen meshpen=meshpen[i];
      if(!invisible(meshpen))
        draw(f,t*project(s.s[i].external(),P),meshpen);
    }
    endgroup(f);
  }
}

void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
          material surfacepen=currentpen, pen meshpen=nullpen,
          light light=currentlight, light meshlight=light,
          projection P=currentprojection)
{
  material[] surfacepen={surfacepen};
  pen[] meshpen={meshpen};
  surfacepen.cyclic(true);
  meshpen.cyclic(true);
  draw(t,f,s,nu,nv,surfacepen,meshpen,light,meshlight,P);
}

void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
          material[] surfacepen, pen[] meshpen=nullpens,
          light light=currentlight, light meshlight=light)
{
  if(s.empty()) return;

  pic.add(new void(frame f, transform3 t, picture pic, projection P) {
      surface S=t*s;
      if(is3D()) {
        draw(f,S,nu,nv,surfacepen,meshpen,light,meshlight);
      } else if(pic != null)
        pic.add(new void(frame f, transform T) {
            draw(T,f,S,nu,nv,surfacepen,meshpen,light,meshlight,P);
          },true);
      if(pic != null) {
        pic.addPoint(min(S,P));
        pic.addPoint(max(S,P));
      }
    },true);
  pic.addPoint(min(s));
  pic.addPoint(max(s));

  pen modifiers;
  if(is3D()) modifiers=thin()+squarecap;
  for(int k=0; k < s.s.length; ++k) {
    pen meshpen=meshpen[k];
    if(!invisible(meshpen)) {
      meshpen=modifiers+meshpen;
      real step=nu == 0 ? 0 : 1/nu;
      for(int i=0; i <= nu; ++i)
        addPath(pic,s.s[k].uequals(i*step),meshpen);
      step=nv == 0 ? 0 : 1/nv;
      for(int j=0; j <= nv; ++j)
        addPath(pic,s.s[k].vequals(j*step),meshpen);
    }
  }
}

void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
          material surfacepen=currentpen, pen meshpen=nullpen,
          light light=currentlight, light meshlight=light)
{
  material[] surfacepen={surfacepen};
  pen[] meshpen={meshpen};
  surfacepen.cyclic(true);
  meshpen.cyclic(true);
  draw(pic,s,nu,nv,surfacepen,meshpen,light,meshlight);
}

void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
          material[] surfacepen, pen meshpen,
          light light=currentlight, light meshlight=light)
{
  pen[] meshpen={meshpen};
  meshpen.cyclic(true);
  draw(pic,s,nu,nv,surfacepen,meshpen,light,meshlight);
}

surface extrude(path p, triple axis=Z)
{
  static patch[] allocate;
  path3 G=path3(p);
  path3 G2=shift(axis)*G;
  return surface(...sequence(new patch(int i) {
        return patch(subpath(G,i,i+1)--subpath(G2,i+1,i)--cycle);
      },length(G)));
}

surface extrude(explicit path[] p, triple axis=Z)
{
  surface s;
  for(path g:p)
    s.append(extrude(g,axis));
  return s;
}

triple rectify(triple dir) 
{
  real scale=max(abs(dir.x),abs(dir.y),abs(dir.z));
  if(scale != 0) dir *= 0.5/scale;
  dir += (0.5,0.5,0.5);
  return dir;
}

path3[] align(path3[] g, transform3 t=identity4, triple position,
              triple align, pen p=currentpen)
{
  if(determinant(t) == 0) return g;
  triple m=min(g);
  triple dir=rectify(inverse(t)*-align);
  triple a=m+realmult(dir,max(g)-m);
  return shift(position+align*labelmargin(p))*t*shift(-a)*g;
}

surface align(surface s, transform3 t=identity4, triple position,
              triple align, pen p=currentpen)
{
  if(determinant(t) == 0) return s;
  triple m=min(s);
  triple dir=rectify(inverse(t)*-align);
  triple a=m+realmult(dir,max(s)-m);
  return shift(position+align*labelmargin(p))*t*shift(-a)*s;
}

surface surface(Label L, triple position=O)
{
  surface s=surface(texpath(L));
  return L.align.is3D ? align(s,L.T3,position,L.align.dir3,L.p) :
    shift(position)*L.T3*s;
}

path[] path(Label L, pair z=0, projection P)
{
  path[] g=texpath(L);
  if(L.defaulttransform3) {
    return L.align.is3D ? align(g,z,project(L.align.dir3,P)-project(O,P),L.p) :
      shift(z)*g;
  } else {
    path3[] G=path3(g);
    return L.align.is3D ? shift(z)*project(align(G,L.T3,O,L.align.dir3,L.p),P) :
      shift(z)*project(L.T3*G,P);
  }
}

void label(frame f, Label L, triple position, align align=NoAlign,
           pen p=currentpen, bool targetsize=false, light light=nolight,
           projection P=currentprojection)
{
  Label L=L.copy();
  L.align(align);
  L.p(p);
  if(targetsize && settings.render != 0)
    L.T=L.T*scale(abs(P.camera-position)/abs(P.vector()));
  if(L.defaulttransform3)
    L.T3=transform3(P);
  if(is3D()) {
    for(patch S : surface(L,position).s)
      draw3D(f,S,L.p,light);
  } else {
    if(L.filltype == NoFill)
      fill(f,path(L,project(position,P.t),P),
           light.color(L.T3*Z,L.p,shiftless(P.modelview())));
    else {
      frame d;
      fill(d,path(L,project(position,P.t),P),
           light.color(L.T3*Z,L.p,shiftless(P.modelview())));
      add(f,d,L.filltype);
    }
  }
}

void label(picture pic=currentpicture, Label L, triple position,
           align align=NoAlign, pen p=currentpen, bool targetsize=false,
           light light=nolight)
{
  Label L=L.copy();
  L.align(align);
  L.p(p);
  L.position(0);
  pic.add(new void(frame f, transform3 t, picture pic, projection P) {
      // Handle relative projected 3D alignments.
      Label L=L.copy();
      triple v=t*position;
      if(!align.is3D && L.align.relative && L.align.dir3 != O &&
         determinant(P.t) != 0)
          L.align(L.align.dir*unit(project(v+L.align.dir3,P.t)-project(v,P.t)));
      
      if(targetsize && settings.render != 0)
        L.T=L.T*scale(abs(P.camera-v)/abs(P.vector()));
      if(L.defaulttransform3)
        L.T3=transform3(P);
      if(is3D())
        for(patch S : surface(L,v).s)
          draw3D(f,S,L.p,light);
      if(pic != null) {
        if(L.filltype == NoFill)
          fill(project(v,P.t),pic,path(L,P),
               light.color(L.T3*Z,L.p,shiftless(P.modelview())));
        else {
          picture d;
          fill(project(v,P.t),d,path(L,P),
               light.color(L.T3*Z,L.p,shiftless(P.modelview())));
          add(pic,d,L.filltype);
        }
      }
    },!L.defaulttransform3);

  Label L=L.copy();
  if(targetsize && settings.render != 0)
    L.T=L.T*scale(abs(currentprojection.camera-position)/
                  abs(currentprojection.vector()));
  path[] g=texpath(L);
  if(g.length == 0 || (g.length == 1 && size(g[0]) == 0)) return;
  if(L.defaulttransform3)
    L.T3=transform3(currentprojection);
  path3[] G=path3(g);
  G=L.align.is3D ? align(G,L.T3,O,L.align.dir3,L.p) : L.T3*G;
  pic.addBox(position,position,min(G),max(G));
}

void label(picture pic=currentpicture, Label L, path3 g, align align=NoAlign,
           bool targetsize=false, pen p=currentpen)
{
  Label L=L.copy();
  L.align(align);
  L.p(p);
  bool relative=L.position.relative;
  real position=L.position.position.x;
  if(L.defaultposition) {relative=true; position=0.5;}
  if(relative) position=reltime(g,position);
  if(L.align.default) {
    align a;
    a.init(-I*(position <= sqrtEpsilon ? S :
              position >= length(g)-sqrtEpsilon ? N : E),relative=true);
    a.dir3=dir(g,position); // Pass 3D direction via unused field.
    L.align(a);             
  }
  label(pic,L,point(g,position),targetsize);
}

surface extrude(Label L, triple axis=Z)
{
  Label L=L.copy();
  path[] g=texpath(L);
  surface S=extrude(g,axis);
  surface s=surface(g);
  S.append(s);
  S.append(shift(axis)*s);
  return S;
}

restricted surface nullsurface;

private real a=4/3*(sqrt(2)-1);
private transform3 t=rotate(90,O,Z);
private transform3 t2=t*t;
private transform3 t3=t2*t;
private transform3 i=xscale3(-1)*zscale3(-1);

restricted patch octant1=patch(X{Z}..{-X}Z..Z{Y}..{-Z}Y{X}..{-Y}cycle,
                               new triple[] {(1,a,a),(a,a^2,1),(a^2,a,1),
                                             (a,1,a)});

restricted surface unithemisphere=surface(octant1,t*octant1,t2*octant1,
                                          t3*octant1);
restricted surface unitsphere=surface(octant1,t*octant1,t2*octant1,t3*octant1,
                                      i*octant1,i*t*octant1,i*t2*octant1,
                                      i*t3*octant1);

restricted patch unitfrustum(real t1, real t2)
{
  real s1=interp(t1,t2,1/3);
  real s2=interp(t1,t2,2/3);
  return patch(interp(Z,X,t2)--interp(Z,X,t1){Y}..{-X}interp(Z,Y,t1)--
               interp(Z,Y,t2){X}..{-Y}cycle,
               new triple[] {(s2,s2*a,1-s2),(s1,s1*a,1-s1),(s1*a,s1,1-s1),
                                          (s2*a,s2,1-s2)});
}

// Return a unitcone constructed from n frusta (the final one being degenerate)
surface unitcone(int n=6)
{
  surface unitcone;
  unitcone.s=new patch[4*n];
  real r=1/3;
  for(int i=0; i < n; ++i) {
    patch s=unitfrustum(i < n-1 ? r^(i+1) : 0,r^i);
    unitcone.s[i]=s;
    unitcone.s[n+i]=t*s;
    unitcone.s[2n+i]=t2*s;
    unitcone.s[3n+i]=t3*s;
  }
  return unitcone;
}

restricted surface unitcone=unitcone();
restricted surface unitsolidcone=surface(patch(unitcircle3)...unitcone.s);

private patch unitcylinder1=patch(X--X+Z{Y}..{-X}Y+Z--Y{X}..{-Y}cycle);

restricted surface unitcylinder=surface(unitcylinder1,t*unitcylinder1,
                                        t2*unitcylinder1,t3*unitcylinder1);

private patch unitplane=patch(new triple[] {O,X,X+Y,Y});
restricted surface unitcube=surface(reverse(unitplane),
                                    rotate(90,O,X)*unitplane,
                                    rotate(-90,O,Y)*unitplane,
                                    shift(Z)*unitplane,
                                    rotate(90,X,X+Y)*unitplane,
                                    rotate(-90,Y,X+Y)*unitplane);
restricted surface unitplane=surface(unitplane);
restricted surface unitdisk=surface(unitcircle3);

void dot(frame f, triple v, material p=currentpen,
         light light=nolight, projection P=currentprojection)
{
  pen q=(pen) p;
  if(is3D()) {
    material m=material(p,p.granularity >= 0 ? p.granularity : dotgranularity);
    for(patch s : unitsphere.s)
      draw3D(f,shift(v)*scale3(0.5*linewidth(dotsize(q)+q))*s,m,light);
  } else dot(f,project(v,P.t),q);
}

void dot(frame f, path3 g, material p=currentpen,
         projection P=currentprojection)
{
  for(int i=0; i <= length(g); ++i) dot(f,point(g,i),p,P);
}

void dot(frame f, path3[] g, material p=currentpen,
         projection P=currentprojection)
{
  for(int i=0; i < g.length; ++i) dot(f,g[i],p,P);
}

void dot(picture pic=currentpicture, triple v, material p=currentpen,
         light light=nolight)
{
  pen q=(pen) p;
  real size=0.5*linewidth(dotsize(q)+q);
  pic.add(new void(frame f, transform3 t, picture pic, projection P) {
      if(is3D()) {
        material m=material(p,p.granularity >= 0 ? p.granularity :
                            dotgranularity);
        for(patch s : unitsphere.s)
          draw3D(f,shift(t*v)*scale3(size)*s,m,light);
      }
      if(pic != null)
        dot(pic,project(t*v,P.t),q);
    },true);
  triple R=size*(1,1,1);
  pic.addBox(v,v,-R,R);
}

void dot(picture pic=currentpicture, triple[] v, material p=currentpen)
{
  for(int i=0; i < v.length; ++i) dot(pic,v[i],p);
}

void dot(picture pic=currentpicture, explicit path3 g, material p=currentpen)
{
  for(int i=0; i <= length(g); ++i) dot(pic,point(g,i),p);
}

void dot(picture pic=currentpicture, path3[] g, material p=currentpen)
{
  for(int i=0; i < g.length; ++i) dot(pic,g[i],p);
}

void dot(picture pic=currentpicture, Label L, triple v, align align=NoAlign,
         string format=defaultformat, material p=currentpen)
{
  Label L=L.copy();
  if(L.s == "") {
    if(format == "") format=defaultformat;
    L.s="("+format(format,v.x)+","+format(format,v.y)+","+
      format(format,v.z)+")";
  }
  L.align(align,E);
  L.p((pen) p);
  dot(pic,v,p);
  label(pic,L,v);
}