summaryrefslogtreecommitdiff
path: root/graphics/asymptote/base/three_tube.asy
diff options
context:
space:
mode:
Diffstat (limited to 'graphics/asymptote/base/three_tube.asy')
-rw-r--r--graphics/asymptote/base/three_tube.asy490
1 files changed, 146 insertions, 344 deletions
diff --git a/graphics/asymptote/base/three_tube.asy b/graphics/asymptote/base/three_tube.asy
index 25bb9385ec..86451b97fa 100644
--- a/graphics/asymptote/base/three_tube.asy
+++ b/graphics/asymptote/base/three_tube.asy
@@ -1,41 +1,15 @@
-void render(path3 s, void f(path3, real), render render=defaultrender)
-{
- real granularity=render.tubegranularity;
- void Split(triple z0, triple c0, triple c1, triple z1, real t0=0, real t1=1,
- real depth=mantissaBits) {
- if(depth > 0) {
- real S=straightness(z0,c0,c1,z1);
- if(S > 0) {
- --depth;
- if(S > max(granularity*max(abs(z0),abs(c0),abs(c1),abs(z1)))) {
- 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);
- real tm=0.5*(t0+t1);
- Split(z0,m0,m3,m5,t0,tm,depth);
- Split(m5,m4,m2,z1,tm,t1,depth);
- return;
- }
- }
- }
- f(z0..controls c0 and c1..z1,t0);
- }
- Split(point(s,0),postcontrol(s,0),precontrol(s,1),point(s,1));
-}
-
-struct rmf
-{
+struct rmf {
triple p,r,t,s;
- void operator init(triple p, triple r, triple t)
- {
+ 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
@@ -66,344 +40,172 @@ rmf[] rmf(path3 g, real[] t)
return R;
}
-private real[][][] bispline0(real[][] z, real[][] p, real[][] q, real[][] r,
- real[] x, real[] y, bool[][] cond={})
-{ // z[i][j] is the value at (x[i],y[j])
- // p and q are the first derivatives with respect to x and y, respectively
- // r is the second derivative ddu/dxdy
- int n=x.length-1;
- int m=y.length-1;
-
- bool all=cond.length == 0;
-
- int count;
- if(all)
- count=n*m;
- else {
- count=0;
- for(int i=0; i < n; ++i) {
- bool[] condi=cond[i];
- bool[] condp=cond[i+1];
- for(int j=0; j < m; ++j)
- if(all || (condi[j] && condi[j+1] && condp[j] && condp[j+1]))
- ++count;
- }
- }
-
- real[][][] s=new real[count][][];
- int k=0;
- for(int i=0; i < n; ++i) {
- int ip=i+1;
- real xi=x[i];
- real xp=x[ip];
- real hx=(xp-xi)/3;
- real[] zi=z[i];
- real[] zp=z[ip];
- real[] ri=r[i];
- real[] rp=r[ip];
- real[] pi=p[i];
- real[] pp=p[ip];
- real[] qi=q[i];
- real[] qp=q[ip];
- bool[] condi=all ? null : cond[i];
- bool[] condp=all ? null : cond[i+1];
- for(int j=0; j < m; ++j) {
- if(all || (condi[j] && condi[j+1] && condp[j] && condp[j+1])) {
- real yj=y[j];
- int jp=j+1;
- real yp=y[jp];
- real hy=(yp-yj)/3;
- real hxy=hx*hy;
- real zij=zi[j];
- real zip=zi[jp];
- real zpj=zp[j];
- real zpp=zp[jp];
- real pij=hx*pi[j];
- real ppj=hx*pp[j];
- real qip=hy*qi[jp];
- real qpp=hy*qp[jp];
- real zippip=zip+hx*pi[jp];
- real zppmppp=zpp-hx*pp[jp];
- real zijqij=zij+hy*qi[j];
- real zpjqpj=zpj+hy*qp[j];
-
- s[k]=new real[][] {{zij,zijqij,zip-qip,zip},
- {zij+pij,zijqij+pij+hxy*ri[j],
- zippip-qip-hxy*ri[jp],zippip},
- {zpj-ppj,zpjqpj-ppj-hxy*rp[j],
- zppmppp-qpp+hxy*rp[jp],zppmppp},
- {zpj,zpjqpj,zpp-qpp,zpp}};
- ++k;
- }
+rmf[] rmf(triple z0, triple c0, triple c1, triple z1, real[] t)
+{
+ 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);
}
-
- return s;
-}
-// return the surface values described by a real matrix f, interpolated with
-// xsplinetype and ysplinetype.
-real[][][] bispline(real[][] f, real[] x, real[] y,
- splinetype xsplinetype=null,
- splinetype ysplinetype=xsplinetype, bool[][] cond={})
-{
- real epsilon=sqrtEpsilon*norm(y);
- if(xsplinetype == null)
- xsplinetype=(abs(x[0]-x[x.length-1]) <= epsilon) ? periodic : notaknot;
- if(ysplinetype == null)
- ysplinetype=(abs(y[0]-y[y.length-1]) <= epsilon) ? periodic : notaknot;
- int n=x.length; int m=y.length;
- real[][] ft=transpose(f);
- real[][] tp=new real[m][];
- for(int j=0; j < m; ++j)
- tp[j]=xsplinetype(x,ft[j]);
- real[][] q=new real[n][];
- for(int i=0; i < n; ++i)
- q[i]=ysplinetype(y,f[i]);
- real[][] qt=transpose(q);
- real[] d1=xsplinetype(x,qt[0]);
- real[] d2=xsplinetype(x,qt[m-1]);
- real[][] r=new real[n][];
- real[][] p=transpose(tp);
- for(int i=0; i < n; ++i)
- r[i]=clamped(d1[i],d2[i])(y,p[i]);
- return bispline0(f,p,q,r,x,y,cond);
-}
-
-bool uperiodic(real[][] a) {
- int n=a.length;
- if(n == 0) return false;
- int m=a[0].length;
- real[] a0=a[0];
- real[] a1=a[n-1];
- for(int j=0; j < m; ++j) {
- real norm=0;
- for(int i=0; i < n; ++i)
- norm=max(norm,abs(a[i][j]));
- real epsilon=sqrtEpsilon*norm;
- if(abs(a0[j]-a1[j]) > epsilon) return false;
+ rmf[] R=new rmf[t.length];
+ triple T=c0-z0;
+ if(abs(T) < norm) {
+ T=z0-2*c0+c1;
+ if(abs(T) < norm)
+ T=z1-z0+3.0*(c0-c1);
}
- return true;
-}
-bool vperiodic(real[][] a) {
- int n=a.length;
- if(n == 0) return false;
- int m=a[0].length-1;
- for(int i=0; i < n; ++i) {
- real[] ai=a[i];
- real epsilon=sqrtEpsilon*norm(ai);
- if(abs(ai[0]-ai[m]) > epsilon) return false;
+ T=unit(T);
+ triple Tp=perp(T);
+ 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];
}
- return true;
+ return R;
}
-// return the surface described by a parametric function f evaluated at u and v
-// and interpolated with usplinetype and vsplinetype.
-surface surface(triple f(pair z), real[] u, real[] v,
- splinetype[] usplinetype, splinetype[] vsplinetype=Spline,
- bool cond(pair z)=null)
+surface tube(triple z0, triple c0, triple c1, triple z1, real w)
{
- int nu=u.length-1;
- int nv=v.length-1;
- real[] ipt=sequence(u.length);
- real[] jpt=sequence(v.length);
- real[][] fx=new real[u.length][v.length];
- real[][] fy=new real[u.length][v.length];
- real[][] fz=new real[u.length][v.length];
-
- bool[][] active;
- bool all=cond == null;
- if(!all) active=new bool[u.length][v.length];
-
- for(int i=0; i <= nu; ++i) {
- real ui=u[i];
- real[] fxi=fx[i];
- real[] fyi=fy[i];
- real[] fzi=fz[i];
- bool[] activei=all ? null : active[i];
- for(int j=0; j <= nv; ++j) {
- pair z=(ui,v[j]);
- if(!all) activei[j]=cond(z);
- triple f=f(z);
- fxi[j]=f.x;
- fyi[j]=f.y;
- fzi[j]=f.z;
- }
- }
-
- if(usplinetype.length == 0) {
- usplinetype=new splinetype[] {uperiodic(fx) ? periodic : notaknot,
- uperiodic(fy) ? periodic : notaknot,
- uperiodic(fz) ? periodic : notaknot};
- } else if(usplinetype.length != 3) abort("usplinetype must have length 3");
-
- if(vsplinetype.length == 0) {
- vsplinetype=new splinetype[] {vperiodic(fx) ? periodic : notaknot,
- vperiodic(fy) ? periodic : notaknot,
- vperiodic(fz) ? periodic : notaknot};
- } else if(vsplinetype.length != 3) abort("vsplinetype must have length 3");
-
- real[][][] sx=bispline(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active);
- real[][][] sy=bispline(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active);
- real[][][] sz=bispline(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active);
-
- surface s=surface(sx.length);
- s.index=new int[nu][nv];
- int k=-1;
- for(int i=0; i < nu; ++i) {
- int[] indexi=s.index[i];
- for(int j=0; j < nv; ++j)
- indexi[j]=++k;
- }
-
- for(int k=0; k < sx.length; ++k) {
- triple[][] Q=new triple[4][];
- real[][] Px=sx[k];
- real[][] Py=sy[k];
- real[][] Pz=sz[k];
- for(int i=0; i < 4 ; ++i) {
- real[] Pxi=Px[i];
- real[] Pyi=Py[i];
- real[] Pzi=Pz[i];
- Q[i]=new triple[] {(Pxi[0],Pyi[0],Pzi[0]),
- (Pxi[1],Pyi[1],Pzi[1]),
- (Pxi[2],Pyi[2],Pzi[2]),
- (Pxi[3],Pyi[3],Pzi[3])};
+ 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.s[k]=patch(Q);
+ s.push(patch(P,copy=false));
}
- if(usplinetype[0] == periodic && usplinetype[1] == periodic &&
- usplinetype[1] == periodic) s.ucyclic(true);
-
- if(vsplinetype[0] == periodic && vsplinetype[1] == periodic &&
- vsplinetype[1] == periodic) s.vcyclic(true);
-
+ 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;
}
-path3 interp(path3 a, path3 b, real t)
-{
- int n=size(a);
- return path3(sequence(new triple(int i) {
- return interp(precontrol(a,i),precontrol(b,i),t);},n),
- sequence(new triple(int i) {return interp(point(a,i),point(b,i),t);},n),
- sequence(new triple(int i) {return interp(postcontrol(a,i),
- postcontrol(b,i),t);},n),
- sequence(new bool(int i) {return straight(a,i) && straight(b,i);},n),
- cyclic(a) && cyclic(b));
+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;
+ surface[] s;
path3 center; // tube axis
void Null(transform3) {}
void Null(transform3, bool) {}
- void operator init(path3 p, real width, render render=defaultrender,
- void cylinder(transform3)=Null,
- void sphere(transform3, bool half)=Null,
- void pipe(path3, path3)=null) {
+ 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));
+ surface[] s;
+ void Split(triple z0, triple c0, triple c1, triple z1,
+ real 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);
- if(piecewisestraight(p)) {
- for(int i=0; i < n; ++i) {
+ 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));
- s.append(t*unitcylinder);
- cylinder(t);
- }
- center=center&p;
- } else {
- real[] T;
- path3 G;
- for(int i=0; i < n; ++i)
- render(subpath(p,i,i+1),
- new void(path3 g, real s) {
- G=G&g;
- T.push(i+s);
- },render);
- T.push(n);
- T.cyclic=cyclic(p);
- rmf[] rmf=rmf(p,T);
- triple f(pair t) {
- rmf R=rmf[round(t.x)];
- int n=round(t.y);
- static real[] x={1,0,-1,0};
- static real[] y={0,1,0,-1};
- return point(G,t.x)+r*(R.r*x[n]-R.s*y[n]);
- }
-
- static real[] v={0,1,2,3,0};
- static real[] circular(real[] x, real[] y) {
- static real a=8/3*(sqrt(2)-1);
- return a*periodic(x,y);
- }
-
- static splinetype[] Monotonic={monotonic,monotonic,monotonic};
- static splinetype[] Circular={circular,circular,circular};
- if(T.length > 0) {
- surface S=surface(f,sequence(T.length),v,Monotonic,Circular);
- s.append(S);
-
- // Compute center of tube:
- int n=S.index.length;
- if(T.cyclic) --n;
- triple[] pre=new triple[n+1];
- triple[] point=new triple[n+1];
- triple[] post=new triple[n+1];
-
- int[] index=S.index[0];
- triple Point;
- for(int m=0; m < 4; ++m)
- Point += S.s[index[m]].P[0][0];
- pre[0]=point[0]=0.25*Point;
-
- for(int i=0; i < n; ++i) {
- index=S.index[i];
- triple Pre,Point,Post;
- for(int m=0; m < 4; ++m) {
- triple [][] P=S.s[index[m]].P;
- Post += P[1][0];
- Pre += P[2][0];
- Point += P[3][0];
- }
- post[i]=0.25*Post;
- pre[i+1]=0.25*Pre;
- point[i+1]=0.25*Point;
-
- }
-
- index=S.index[n-1];
- triple Post;
- for(int m=0; m < 4; ++m)
- Post += S.s[index[m]].P[3][0];
- post[n]=0.25*Post;
-
- bool[] b=array(n+1,false);
- path3 Center=path3(pre,point,post,b,T.cyclic);
- center=center&Center;
-
- if(pipe != null) { // Compute path along tube
- triple[] pre=new triple[n+1];
- triple[] point=new triple[n+1];
- triple[] post=new triple[n+1];
- pre[0]=point[0]=S.s[S.index[0][0]].P[0][0];
- for(int i=0; i < n; ++i) {
- triple [][] P=S.s[S.index[i][0]].P;
- post[i]=P[1][0];
- pre[i+1]=P[2][0];
- point[i+1]=P[3][0];
- }
- post[n]=S.s[S.index[n-1][0]].P[3][0];
- pipe(Center,path3(pre,point,post,b,T.cyclic));
- }
- }
+ // 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));
}
}
@@ -416,11 +218,11 @@ struct tube
generate(subpath(p,begin,i));
triple dir=dir(p,i,-1);
transform3 T=t*align(dir);
- s.append(shift(point(p,i))*T*(dir != O ? unithemisphere : unitsphere));
- sphere(shift(point(center,length(center)))*T,
- half=straight(p,i-1) && straight(p,i));
+ s.push(shift(point(p,i))*T*(straight(p,i-1) && straight(p,i) ?
+ unithemisphere : unitsphere));
begin=i;
}
- generate(subpath(p,begin,n));
+ path3 g=subpath(p,begin,n);
+ generate(g);
}
}