summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/base/three_tube.asy
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/base/three_tube.asy')
-rw-r--r--Build/source/utils/asymptote/base/three_tube.asy234
1 files changed, 234 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/base/three_tube.asy b/Build/source/utils/asymptote/base/three_tube.asy
new file mode 100644
index 00000000000..60085a7d512
--- /dev/null
+++ b/Build/source/utils/asymptote/base/three_tube.asy
@@ -0,0 +1,234 @@
+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;
+}
+
+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=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(s),t*max(s),m.opacity == 1);
+ };
+ 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);
+ }
+}