diff options
Diffstat (limited to 'Master/texmf-dist/asymptote/three_tube.asy')
-rw-r--r-- | Master/texmf-dist/asymptote/three_tube.asy | 490 |
1 files changed, 146 insertions, 344 deletions
diff --git a/Master/texmf-dist/asymptote/three_tube.asy b/Master/texmf-dist/asymptote/three_tube.asy index 25bb9385ec8..86451b97fa1 100644 --- a/Master/texmf-dist/asymptote/three_tube.asy +++ b/Master/texmf-dist/asymptote/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); } } |