struct rmf {
  triple p,r,t,s;
  void operator init(triple p, triple r, triple t) {
    this.p=p;
    this.r=r;
    this.t=t;
    s=cross(t,r);
  }

  transform3 transform() {
    return transform3(r,s,t);
  }
}

// Rotation minimizing frame
// http://www.cs.hku.hk/research/techreps/document/TR-2007-07.pdf
rmf[] rmf(path3 g, real[] t, triple perp=O)
{
  triple T=dir(g,0);
  triple Tp=abs(perp) < sqrtEpsilon ? perp(T) : unit(perp);
  rmf[] R=new rmf[t.length];
  R[0]=rmf(point(g,0),Tp,T);
  for(int i=1; i < t.length; ++i) {
    rmf Ri=R[i-1];
    real t=t[i];
    triple p=point(g,t);
    triple v1=p-Ri.p;
    if(v1 != O) {
      triple r=Ri.r;
      triple u1=unit(v1);
      triple ti=Ri.t;
      triple tp=ti-2*dot(u1,ti)*u1;
      ti=dir(g,t);
      triple rp=r-2*dot(u1,r)*u1;
      triple u2=unit(ti-tp);
      rp=rp-2*dot(u2,rp)*u2;
      R[i]=rmf(p,unit(rp),unit(ti));
    } else
      R[i]=R[i-1];
  }
  return R;
}

rmf[] rmf(triple z0, triple c0, triple c1, triple z1, real[] t, triple perp=O)
{
  static triple s0;

  real norm=sqrtEpsilon*max(abs(z0),abs(c0),abs(c1),abs(z1));

  // Special case of dir for t in (0,1].
  triple dir(real t) {
    if(t == 1) {
      triple dir=z1-c1;
      if(abs(dir) > norm) return unit(dir);
      dir=2.0*c1-c0-z1;
      if(abs(dir) > norm) return unit(dir);
      return unit(z1-z0+3.0*(c0-c1));
    }
    triple a=z1-z0+3.0*(c0-c1);
    triple b=2.0*(z0+c1)-4.0*c0;
    triple c=c0-z0;
    triple dir=a*t*t+b*t+c;
    if(abs(dir) > norm) return unit(dir);
    dir=2.0*a*t+b;
    if(abs(dir) > norm) return unit(dir);
    return unit(a);
  }

  triple T=c0-z0;
  if(abs(T) < norm) {
    T=z0-2*c0+c1;
    if(abs(T) < norm)
      T=z1-z0+3.0*(c0-c1);
  }
  T=unit(T);
  triple Tp=perp == O ? cross(s0,T) : perp;
  Tp=abs(Tp) < sqrtEpsilon ? perp(T) : unit(Tp);
  rmf[] R=new rmf[t.length];
  R[0]=rmf(z0,Tp,T);

  for(int i=1; i < t.length; ++i) {
    rmf Ri=R[i-1];
    real t=t[i];
    triple p=bezier(z0,c0,c1,z1,t);
    triple v1=p-Ri.p;
    if(v1 != O) {
      triple r=Ri.r;
      triple u1=unit(v1);
      triple ti=Ri.t;
      triple tp=ti-2*dot(u1,ti)*u1;
      ti=dir(t);
      triple rp=r-2*dot(u1,r)*u1;
      triple u2=unit(ti-tp);
      rp=rp-2*dot(u2,rp)*u2;
      R[i]=rmf(p,unit(rp),unit(ti));
    } else
      R[i]=R[i-1];
  }
  s0=R[t.length-1].s;
  return R;
}


drawfcn drawTube(triple[] g, real w, triple min, triple max) {
  return new void(frame f, transform3 t=identity4, material[] m,
                  light light=currentlight, render render=defaultrender)
    {
      material m=material(m[0],light);
      drawTube(f,t*g,w,m.p,m.opacity,m.shininess,m.metallic,m.fresnel0,
               t*min,t*max,m.opacity == 1);
    };
}

surface tube(triple z0, triple c0, triple c1, triple z1, real w)
{
  surface s;
  static real[] T={0,1/3,2/3,1};
  rmf[] rmf=rmf(z0,c0,c1,z1,T);

  real aw=a*w;
  triple[] arc={(w,0,0),(w,aw,0),(aw,w,0),(0,w,0)};
  triple[] g={z0,c0,c1,z1};

  void f(transform3 R) {
    triple[][] P=new triple[4][];
    for(int i=0; i < 4; ++i) {
      transform3 T=shift(g[i])*rmf[i].transform()*R;
      P[i]=new triple[] {T*arc[0],T*arc[1],T*arc[2],T*arc[3]};
    }
    s.push(patch(P,copy=false));
  }

  f(identity4);
  f(t1);
  f(t2);
  f(t3);

  s.PRCprimitive=false;
  s.draw=drawTube(g,w,min(s),max(s));
  return s;
}

real tubethreshold=20;

// Note: casting an array of surfaces to a single surface will disable
// primitive compression.
surface operator cast(surface[] s) {
  surface S;
  for(surface p : s)
    S.append(p);
  return S;
}

struct tube
{
  surface[] s;
  path3 center; // tube axis

  void Null(transform3) {}
  void Null(transform3, bool) {}

  surface[] render(path3 g, real r) {
    triple z0=point(g,0);
    triple c0=postcontrol(g,0);
    triple c1=precontrol(g,1);
    triple z1=point(g,1);
    real norm=sqrtEpsilon*max(abs(z0),abs(c0),abs(c1),abs(z1),r);
    surface[] s;
    void Split(triple z0, triple c0, triple c1, triple z1,
               int depth=mantissaBits) {
      if(depth > 0) {
        pair threshold(triple z0, triple c0, triple c1) {
          triple u=c1-z0;
          triple v=c0-z0;
          real x=abs(v);
          return (x,abs(u*x^2-dot(u,v)*v));
        }

        pair a0=threshold(z0,c0,c1);
        pair a1=threshold(z1,c1,c0);
        real rL=r*arclength(z0,c0,c1,z1)*tubethreshold;
        if((a0.x >= norm && rL*a0.y^2 > a0.x^8) ||
           (a1.x >= norm && rL*a1.y^2 > a1.x^8)) {
          triple m0=0.5*(z0+c0);
          triple m1=0.5*(c0+c1);
          triple m2=0.5*(c1+z1);
          triple m3=0.5*(m0+m1);
          triple m4=0.5*(m1+m2);
          triple m5=0.5*(m3+m4);
          --depth;
          Split(z0,m0,m3,m5,depth);
          Split(m5,m4,m2,z1,depth);
          return;
        }
      }

      s.push(tube(z0,c0,c1,z1,r));
    }
    Split(z0,c0,c1,z1);
    return s;
  }

  void operator init(path3 p, real width) {
    center=p;
    real r=0.5*width;

    void generate(path3 p) {
      int n=length(p);
      for(int i=0; i < n; ++i) {
        if(straight(p,i)) {
          triple v=point(p,i);
          triple u=point(p,i+1)-v;
          transform3 t=shift(v)*align(unit(u))*scale(r,r,abs(u));
          // Draw opaque surfaces with core for better small-scale rendering.
          surface unittube=t*unitcylinder;
          unittube.draw=unitcylinderDraw(core=true);
          s.push(unittube);
        } else
          s.append(render(subpath(p,i,i+1),r));
      }
    }

    transform3 t=scale3(r);
    bool cyclic=cyclic(p);
    int begin=0;
    int n=length(p);
    for(int i=cyclic ? 0 : 1; i < n; ++i)
      if(abs(dir(p,i,1)-dir(p,i,-1)) > sqrtEpsilon) {
        generate(subpath(p,begin,i));
        triple dir=dir(p,i,-1);
        transform3 T=t*align(dir);
        s.push(shift(point(p,i))*T*(straight(p,i-1) && straight(p,i) ?
                                    unithemisphere : unitsphere));
        begin=i;
      }
    path3 g=subpath(p,begin,n);
    generate(g);
  }
}