diff options
author | Karl Berry <karl@freefriends.org> | 2013-04-07 18:19:31 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2013-04-07 18:19:31 +0000 |
commit | 752012c605d34cd943795527a9738475a6958fcc (patch) | |
tree | 4ee06acdd8333a662c2d6f6ef716235053468f55 /Master/texmf/asymptote/three.asy | |
parent | 9789d09132f18a838e84f041b4b3aff28d3426ec (diff) |
texmf -> texmf-dist: start with unique dirs from texmf
git-svn-id: svn://tug.org/texlive/trunk@29712 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf/asymptote/three.asy')
-rw-r--r-- | Master/texmf/asymptote/three.asy | 3319 |
1 files changed, 0 insertions, 3319 deletions
diff --git a/Master/texmf/asymptote/three.asy b/Master/texmf/asymptote/three.asy deleted file mode 100644 index 5d7ae1a8ba5..00000000000 --- a/Master/texmf/asymptote/three.asy +++ /dev/null @@ -1,3319 +0,0 @@ -private import math; - -if(inXasyMode) settings.render=0; - -if(prc0()) { - if(settings.tex == "context") settings.prc=false; - else { - access embed; - Embed=embed.embed; - Link=embed.link; - } -} - -// Useful lossy compression values. -restricted real Zero=0; -restricted real Low=0.0001; -restricted real Medium=0.001; -restricted real High=0.01; - -restricted int PRCsphere=0; // Renders slowly but produces smaller PRC files. -restricted int NURBSsphere=1; // Renders fast but produces larger PRC files. - -struct render -{ - // PRC parameters: - real compression; // lossy compression parameter (0=no compression) - real granularity; // PRC rendering granularity - - bool closed; // use one-sided rendering? - bool tessellate; // use tessellated mesh to store straight patches? - - bool3 merge; // merge nodes before rendering, for faster but - // lower quality PRC rendering (the value default means - // merge opaque patches only). - - int sphere; // PRC sphere type (PRCsphere or NURBSsphere). - - // General parameters: - real margin; // shrink amount for rendered openGL viewport, in bp. - real tubegranularity; // granularity for rendering tubes - bool labelfill; // fill subdivision cracks in unlighted labels - - static render defaultrender; - - void operator init(real compression=defaultrender.compression, - real granularity=defaultrender.granularity, - bool closed=defaultrender.closed, - bool tessellate=defaultrender.tessellate, - bool3 merge=defaultrender.merge, - int sphere=defaultrender.sphere, - real margin=defaultrender.margin, - real tubegranularity=defaultrender.tubegranularity, - bool labelfill=defaultrender.labelfill) - { - this.compression=compression; - this.granularity=granularity; - this.closed=closed; - this.tessellate=tessellate; - this.merge=merge; - this.sphere=sphere; - this.margin=margin; - this.tubegranularity=tubegranularity; - this.labelfill=labelfill; - } -} - -render operator init() {return render();} - -render defaultrender=render.defaultrender=new render; -defaultrender.compression=High; -defaultrender.granularity=Medium; -defaultrender.closed=false; -defaultrender.tessellate=false; -defaultrender.merge=false; -defaultrender.margin=0.02; -defaultrender.tubegranularity=0.005; -defaultrender.sphere=NURBSsphere; -defaultrender.labelfill=true; - -real defaultshininess=0.25; - -real angleprecision=1e-5; // Precision for centering perspective projections. -int maxangleiterations=25; - -string defaultembed3Doptions; -string defaultembed3Dscript; -real defaulteyetoview=63mm/1000mm; - -string partname(int i=0) -{ - return string(i+1); -} - -triple O=(0,0,0); -triple X=(1,0,0), Y=(0,1,0), Z=(0,0,1); - -// A translation in 3D space. -transform3 shift(triple v) -{ - transform3 t=identity(4); - t[0][3]=v.x; - t[1][3]=v.y; - t[2][3]=v.z; - return t; -} - -// Avoid two parentheses. -transform3 shift(real x, real y, real z) -{ - return shift((x,y,z)); -} - -transform3 shift(transform3 t) -{ - transform3 T=identity(4); - T[0][3]=t[0][3]; - T[1][3]=t[1][3]; - T[2][3]=t[2][3]; - return T; -} - -// A 3D scaling in the x direction. -transform3 xscale3(real x) -{ - transform3 t=identity(4); - t[0][0]=x; - return t; -} - -// A 3D scaling in the y direction. -transform3 yscale3(real y) -{ - transform3 t=identity(4); - t[1][1]=y; - return t; -} - -// A 3D scaling in the z direction. -transform3 zscale3(real z) -{ - transform3 t=identity(4); - t[2][2]=z; - return t; -} - -// A 3D scaling by s in the v direction. -transform3 scale(triple v, real s) -{ - v=unit(v); - s -= 1; - return new real[][] { - {1+s*v.x^2, s*v.x*v.y, s*v.x*v.z, 0}, - {s*v.x*v.y, 1+s*v.y^2, s*v.y*v.z, 0}, - {s*v.x*v.z, s*v.y*v.z, 1+s*v.z^2, 0}, - {0, 0, 0, 1}}; -} - -// A transformation representing rotation by an angle in degrees about -// an axis v through the origin (in the right-handed direction). -transform3 rotate(real angle, triple v) -{ - if(v == O) abort("cannot rotate about the zero vector"); - v=unit(v); - real x=v.x, y=v.y, z=v.z; - real s=Sin(angle), c=Cos(angle), t=1-c; - - return new real[][] { - {t*x^2+c, t*x*y-s*z, t*x*z+s*y, 0}, - {t*x*y+s*z, t*y^2+c, t*y*z-s*x, 0}, - {t*x*z-s*y, t*y*z+s*x, t*z^2+c, 0}, - {0, 0, 0, 1}}; -} - -// A transformation representing rotation by an angle in degrees about -// the line u--v (in the right-handed direction). -transform3 rotate(real angle, triple u, triple v) -{ - return shift(u)*rotate(angle,v-u)*shift(-u); -} - -// Reflects about the plane through u, v, and w. -transform3 reflect(triple u, triple v, triple w) -{ - triple n=unit(cross(v-u,w-u)); - if(n == O) - abort("points determining reflection plane cannot be colinear"); - - return new real[][] { - {1-2*n.x^2, -2*n.x*n.y, -2*n.x*n.z, u.x}, - {-2*n.x*n.y, 1-2*n.y^2, -2*n.y*n.z, u.y}, - {-2*n.x*n.z, -2*n.y*n.z, 1-2*n.z^2, u.z}, - {0, 0, 0, 1} - }*shift(-u); -} - -// Project u onto v. -triple project(triple u, triple v) -{ - v=unit(v); - return dot(u,v)*v; -} - -// Return a unit vector perpendicular to a given unit vector v. -triple perp(triple v) -{ - triple u=cross(v,Y); - if(abs(u) > sqrtEpsilon) return unit(u); - u=cross(v,Z); - return (abs(u) > sqrtEpsilon) ? unit(u) : X; -} - -// Return the transformation corresponding to moving the camera from the target -// (looking in the negative z direction) to the point 'eye' (looking at target), -// orienting the camera so that direction 'up' points upwards. -// Since, in actuality, we are transforming the points instead of the camera, -// we calculate the inverse matrix. -// Based on the gluLookAt implementation in the OpenGL manual. -transform3 look(triple eye, triple up=Z, triple target=O) -{ - triple f=unit(target-eye); - if(f == O) - f=-Z; // The eye is already at the origin: look down. - - triple s=cross(f,up); - - // If the eye is pointing either directly up or down, there is no - // preferred "up" direction. Pick one arbitrarily. - s=s != O ? unit(s) : perp(f); - - triple u=cross(s,f); - - transform3 M={{ s.x, s.y, s.z, 0}, - { u.x, u.y, u.z, 0}, - {-f.x, -f.y, -f.z, 0}, - { 0, 0, 0, 1}}; - - return M*shift(-eye); -} - -// Return a matrix to do perspective distortion based on a triple v. -transform3 distort(triple v) -{ - transform3 t=identity(4); - real d=length(v); - if(d == 0) return t; - t[3][2]=-1/d; - t[3][3]=0; - return t; -} - -projection operator * (transform3 t, projection P) -{ - projection P=P.copy(); - if(!P.absolute) { - P.camera=t*P.camera; - triple target=P.target; - P.target=t*P.target; - if(P.infinity) - P.normal=t*(target+P.normal)-P.target; - else - P.normal=P.vector(); - P.calculate(); - } - return P; -} - -// With this, save() and restore() in plain also save and restore the -// currentprojection. -addSaveFunction(new restoreThunk() { - projection P=currentprojection.copy(); - return new void() { - currentprojection=P; - }; - }); - -pair project(triple v, projection P=currentprojection) -{ - return project(v,P.t); -} - -pair dir(triple v, triple dir, projection P) -{ - return unit(project(v+0.5dir,P)-project(v-0.5*dir,P)); -} - -// Uses the homogenous coordinate to perform perspective distortion. -// When combined with a projection to the XY plane, this effectively maps -// points in three space to a plane through target and -// perpendicular to the vector camera-target. -projection perspective(triple camera, triple up=Z, triple target=O, - real zoom=1, real angle=0, pair viewportshift=0, - bool showtarget=true, bool autoadjust=true, - bool center=autoadjust) -{ - if(camera == target) - abort("camera cannot be at target"); - return projection(camera,up,target,zoom,angle,viewportshift, - showtarget,autoadjust,center, - new transformation(triple camera, triple up, triple target) - {return transformation(look(camera,up,target), - distort(camera-target));}); -} - -projection perspective(real x, real y, real z, triple up=Z, triple target=O, - real zoom=1, real angle=0, pair viewportshift=0, - bool showtarget=true, bool autoadjust=true, - bool center=autoadjust) -{ - return perspective((x,y,z),up,target,zoom,angle,viewportshift,showtarget, - autoadjust,center); -} - -projection orthographic(triple camera, triple up=Z, triple target=O, - real zoom=1, pair viewportshift=0, - bool showtarget=true, bool center=false) -{ - return projection(camera,up,target,zoom,viewportshift,showtarget, - center=center,new transformation(triple camera, triple up, - triple target) { - return transformation(look(camera,up,target));}); -} - -projection orthographic(real x, real y, real z, triple up=Z, - triple target=O, real zoom=1, pair viewportshift=0, - bool showtarget=true, bool center=false) -{ - return orthographic((x,y,z),up,target,zoom,viewportshift,showtarget, - center=center); -} - -// Compute camera position with x axis below the horizontal at angle alpha, -// y axis below the horizontal at angle beta, and z axis up. -triple camera(real alpha, real beta) -{ - real denom=Tan(alpha+beta); - real Tanalpha=Tan(alpha); - real Tanbeta=Tan(beta); - return (sqrt(Tanalpha/denom),sqrt(Tanbeta/denom),sqrt(Tanalpha*Tanbeta)); -} - -projection oblique(real angle=45) -{ - transform3 t=identity(4); - real c2=Cos(angle)^2; - real s2=1-c2; - t[0][2]=-c2; - t[1][2]=-s2; - t[2][2]=1; - t[2][3]=-1; - return projection((c2,s2,1),up=Y,normal=(0,0,1), - new transformation(triple,triple,triple) { - return transformation(t);}); -} - -projection obliqueZ(real angle=45) {return oblique(angle);} - -projection obliqueX(real angle=45) -{ - transform3 t=identity(4); - real c2=Cos(angle)^2; - real s2=1-c2; - t[0][0]=-c2; - t[1][0]=-s2; - t[1][1]=0; - t[0][1]=1; - t[1][2]=1; - t[2][2]=0; - t[2][0]=1; - t[2][3]=-1; - return projection((1,c2,s2),normal=(1,0,0), - new transformation(triple,triple,triple) { - return transformation(t);}); -} - -projection obliqueY(real angle=45) -{ - transform3 t=identity(4); - real c2=Cos(angle)^2; - real s2=1-c2; - t[0][1]=c2; - t[1][1]=s2; - t[1][2]=1; - t[2][1]=-1; - t[2][2]=0; - t[2][3]=-1; - return projection((c2,-1,s2),normal=(0,-1,0), - new transformation(triple,triple,triple) { - return transformation(t);}); -} - -projection oblique=oblique(); -projection obliqueX=obliqueX(), obliqueY=obliqueY(), obliqueZ=obliqueZ(); - -projection LeftView=orthographic(-X,showtarget=true); -projection RightView=orthographic(X,showtarget=true); -projection FrontView=orthographic(-Y,showtarget=true); -projection BackView=orthographic(Y,showtarget=true); -projection BottomView=orthographic(-Z,up=-Y,showtarget=true); -projection TopView=orthographic(Z,up=Y,showtarget=true); - -currentprojection=perspective(5,4,2); - -projection projection() -{ - projection P; - real[] a=_projection(); - if(a.length == 0 || a[10] == 0.0) return currentprojection; - int k=0; - return a[0] == 1 ? - orthographic((a[++k],a[++k],a[++k]),(a[++k],a[++k],a[++k]), - (a[++k],a[++k],a[++k]),a[++k],(a[k += 2],a[++k])) : - perspective((a[++k],a[++k],a[++k]),(a[++k],a[++k],a[++k]), - (a[++k],a[++k],a[++k]),a[++k],a[++k],(a[++k],a[++k])); -} - -// Map pair z to a triple by inverting the projection P onto the -// plane perpendicular to normal and passing through point. -triple invert(pair z, triple normal, triple point, - projection P=currentprojection) -{ - transform3 t=P.t; - real[][] A={{t[0][0]-z.x*t[3][0],t[0][1]-z.x*t[3][1],t[0][2]-z.x*t[3][2]}, - {t[1][0]-z.y*t[3][0],t[1][1]-z.y*t[3][1],t[1][2]-z.y*t[3][2]}, - {normal.x,normal.y,normal.z}}; - real[] b={z.x*t[3][3]-t[0][3],z.y*t[3][3]-t[1][3],dot(normal,point)}; - real[] x=solve(A,b,warn=false); - return x.length > 0 ? (x[0],x[1],x[2]) : P.camera; -} - -// Map pair to a triple on the projection plane. -triple invert(pair z, projection P=currentprojection) -{ - return invert(z,P.normal,P.target,P); -} - -// Map pair dir to a triple direction at point v on the projection plane. -triple invert(pair dir, triple v, projection P=currentprojection) -{ - return invert(project(v,P)+dir,P.normal,v,P)-v; -} - -pair xypart(triple v) -{ - return (v.x,v.y); -} - -struct control { - triple post,pre; - bool active=false; - bool straight=true; - void operator init(triple post, triple pre, bool straight=false) { - this.post=post; - this.pre=pre; - active=true; - this.straight=straight; - } -} - -control nocontrol; - -control operator * (transform3 t, control c) -{ - control C; - C.post=t*c.post; - C.pre=t*c.pre; - C.active=c.active; - C.straight=c.straight; - return C; -} - -void write(file file, control c) -{ - write(file,".. controls "); - write(file,c.post); - write(file," and "); - write(file,c.pre); -} - -struct Tension { - real out,in; - bool atLeast; - bool active; - void operator init(real out=1, real in=1, bool atLeast=false, - bool active=true) { - real check(real val) { - if(val < 0.75) abort("tension cannot be less than 3/4"); - return val; - } - this.out=check(out); - this.in=check(in); - this.atLeast=atLeast; - this.active=active; - } -} - -Tension operator init() -{ - return Tension(); -} - -Tension noTension; -noTension.active=false; - -void write(file file, Tension t) -{ - write(file,"..tension "); - if(t.atLeast) write(file,"atleast "); - write(file,t.out); - write(file," and "); - write(file,t.in); -} - -struct dir { - triple dir; - real gamma=1; // endpoint curl - bool Curl; // curl specified - bool active() { - return dir != O || Curl; - } - void init(triple v) { - this.dir=v; - } - void init(real gamma) { - if(gamma < 0) abort("curl cannot be less than 0"); - this.gamma=gamma; - this.Curl=true; - } - void init(dir d) { - dir=d.dir; - gamma=d.gamma; - Curl=d.Curl; - } - void default(triple v) { - if(!active()) init(v); - } - void default(dir d) { - if(!active()) init(d); - } - dir copy() { - dir d=new dir; - d.init(this); - return d; - } -} - -void write(file file, dir d) -{ - if(d.dir != O) { - write(file,"{"); write(file,unit(d.dir)); write(file,"}"); - } else if(d.Curl) { - write(file,"{curl "); write(file,d.gamma); write(file,"}"); - } -} - -dir operator * (transform3 t, dir d) -{ - dir D=d.copy(); - D.init(unit(shiftless(t)*d.dir)); - return D; -} - -void checkEmpty(int n) { - if(n == 0) - abort("nullpath3 has no points"); -} - -int adjustedIndex(int i, int n, bool cycles) -{ - checkEmpty(n); - if(cycles) - return i % n; - else if(i < 0) - return 0; - else if(i >= n) - return n-1; - else - return i; -} - -struct flatguide3 { - triple[] nodes; - bool[] cyclic; // true if node is really a cycle - control[] control; // control points for segment starting at node - Tension[] Tension; // Tension parameters for segment starting at node - dir[] in,out; // in and out directions for segment starting at node - - bool cyclic() {int n=cyclic.length; return n > 0 ? cyclic[n-1] : false;} - bool precyclic() {int i=find(cyclic); return i >= 0 && i < cyclic.length-1;} - - int size() { - return cyclic() ? nodes.length-1 : nodes.length; - } - - void node(triple v, bool b=false) { - nodes.push(v); - control.push(nocontrol); - Tension.push(noTension); - in.push(new dir); - out.push(new dir); - cyclic.push(b); - } - - void control(triple post, triple pre) { - if(control.length > 0) { - control c=control(post,pre,false); - control[control.length-1]=c; - } - } - - void Tension(real out, real in, bool atLeast) { - if(Tension.length > 0) - Tension[Tension.length-1]=Tension(out,in,atLeast,true); - } - - void in(triple v) { - if(in.length > 0) { - in[in.length-1].init(v); - } - } - - void out(triple v) { - if(out.length > 0) { - out[out.length-1].init(v); - } - } - - void in(real gamma) { - if(in.length > 0) { - in[in.length-1].init(gamma); - } - } - - void out(real gamma) { - if(out.length > 0) { - out[out.length-1].init(gamma); - } - } - - void cycleToken() { - if(nodes.length > 0) - node(nodes[0],true); - } - - // Return true if outgoing direction at node i is known. - bool solved(int i) { - return out[i].active() || control[i].active; - } -} - -void write(file file, string s="", explicit flatguide3 x, suffix suffix=none) -{ - write(file,s); - if(x.size() == 0) write(file,"<nullpath3>"); - else for(int i=0; i < x.nodes.length; ++i) { - if(i > 0) write(file,endl); - if(x.cyclic[i]) write(file,"cycle"); - else write(file,x.nodes[i]); - if(i < x.nodes.length-1) { - // Explicit control points trump other specifiers - if(x.control[i].active) - write(file,x.control[i]); - else { - write(file,x.out[i]); - if(x.Tension[i].active) write(file,x.Tension[i]); - } - write(file,".."); - if(!x.control[i].active) write(file,x.in[i]); - } - } - write(file,suffix); -} - -void write(string s="", flatguide3 x, suffix suffix=endl) -{ - write(stdout,s,x,suffix); -} - -// A guide3 is most easily represented as something that modifies a flatguide3. -typedef void guide3(flatguide3); - -restricted void nullpath3(flatguide3) {}; - -guide3 operator init() {return nullpath3;} - -guide3 operator cast(triple v) -{ - return new void(flatguide3 f) { - f.node(v); - }; -} - -guide3 operator cast(cycleToken) { - return new void(flatguide3 f) { - f.cycleToken(); - }; -} - -guide3 operator controls(triple post, triple pre) -{ - return new void(flatguide3 f) { - f.control(post,pre); - }; -}; - -guide3 operator controls(triple v) -{ - return operator controls(v,v); -} - -guide3 operator cast(tensionSpecifier t) -{ - return new void(flatguide3 f) { - f.Tension(t.out, t.in, t.atLeast); - }; -} - -guide3 operator cast(curlSpecifier spec) -{ - return new void(flatguide3 f) { - if(spec.side == JOIN_OUT) f.out(spec.value); - else if(spec.side == JOIN_IN) f.in(spec.value); - else - abort("invalid curl specifier"); - }; -} - -guide3 operator spec(triple v, int side) -{ - return new void(flatguide3 f) { - if(side == JOIN_OUT) f.out(v); - else if(side == JOIN_IN) f.in(v); - else - abort("invalid direction specifier"); - }; -} - -guide3 operator -- (... guide3[] g) -{ - return new void(flatguide3 f) { - if(g.length > 0) { - for(int i=0; i < g.length-1; ++i) { - g[i](f); - f.out(1); - f.in(1); - } - g[g.length-1](f); - } - }; -} - -guide3 operator .. (... guide3[] g) -{ - return new void(flatguide3 f) { - for(int i=0; i < g.length; ++i) - g[i](f); - }; -} - -typedef guide3 interpolate3(... guide3[]); - -interpolate3 join3(tensionSpecifier t) -{ - return new guide3(... guide3[] a) { - if(a.length == 0) return nullpath3; - guide3 g=a[0]; - for(int i=1; i < a.length; ++i) - g=g..t..a[i]; - return g; - }; -} - -interpolate3 operator ::=join3(operator tension(1,true)); -interpolate3 operator ---=join3(operator tension(infinity,true)); - -flatguide3 operator cast(guide3 g) -{ - flatguide3 f; - g(f); - return f; -} - -flatguide3[] operator cast(guide3[] g) -{ - flatguide3[] p=new flatguide3[g.length]; - for(int i=0; i < g.length; ++i) { - flatguide3 f; - g[i](f); - p[i]=f; - } - return p; -} - -// A version of asin that tolerates numerical imprecision -real asin1(real x) -{ - return asin(min(max(x,-1),1)); -} - -// A version of acos that tolerates numerical imprecision -real acos1(real x) -{ - return acos(min(max(x,-1),1)); -} - -struct Controls { - triple c0,c1; - - // 3D extension of John Hobby's control point formula - // (cf. The MetaFont Book, page 131), - // as described in John C. Bowman and A. Hammerlindl, - // TUGBOAT: The Communications of the TeX Users Group 29:2 (2008). - - void operator init(triple v0, triple v1, triple d0, triple d1, real tout, - real tin, bool atLeast) { - triple v=v1-v0; - triple u=unit(v); - real L=length(v); - d0=unit(d0); - d1=unit(d1); - real theta=acos1(dot(d0,u)); - real phi=acos1(dot(d1,u)); - if(dot(cross(d0,v),cross(v,d1)) < 0) phi=-phi; - c0=v0+d0*L*relativedistance(theta,phi,tout,atLeast); - c1=v1-d1*L*relativedistance(phi,theta,tin,atLeast); - } -} - -private triple cross(triple d0, triple d1, triple reference) -{ - triple normal=cross(d0,d1); - return normal == O ? reference : normal; -} - -private triple dir(real theta, triple d0, triple d1, triple reference) -{ - triple normal=cross(d0,d1,reference); - if(normal == O) return d1; - return rotate(degrees(theta),dot(normal,reference) >= 0 ? normal : -normal)* - d1; -} - -private real angle(triple d0, triple d1, triple reference) -{ - real theta=acos1(dot(unit(d0),unit(d1))); - return dot(cross(d0,d1,reference),reference) >= 0 ? theta : -theta; -} - -// 3D extension of John Hobby's angle formula (The MetaFont Book, page 131). -// Notational differences: here psi[i] is the turning angle at z[i+1], -// beta[i] is the tension for segment i, and in[i] is the incoming -// direction for segment i (where segment i begins at node i). - -real[] theta(triple[] v, real[] alpha, real[] beta, - triple dir0, triple dirn, real g0, real gn, triple reference) -{ - real[] a,b,c,f,l,psi; - int n=alpha.length; - bool cyclic=v.cyclic; - for(int i=0; i < n; ++i) - l[i]=1/length(v[i+1]-v[i]); - int i0,in; - if(cyclic) {i0=0; in=n;} - else {i0=1; in=n-1;} - for(int i=0; i < in; ++i) - psi[i]=angle(v[i+1]-v[i],v[i+2]-v[i+1],reference); - if(cyclic) { - l.cyclic=true; - psi.cyclic=true; - } else { - psi[n-1]=0; - if(dir0 == O) { - real a0=alpha[0]; - real b0=beta[0]; - real chi=g0*(b0/a0)^2; - a[0]=0; - b[0]=3a0-a0/b0+chi; - real C=chi*(3a0-1)+a0/b0; - c[0]=C; - f[0]=-C*psi[0]; - } else { - a[0]=c[0]=0; - b[0]=1; - f[0]=angle(v[1]-v[0],dir0,reference); - } - if(dirn == O) { - real an=alpha[n-1]; - real bn=beta[n-1]; - real chi=gn*(an/bn)^2; - a[n]=chi*(3bn-1)+bn/an; - b[n]=3bn-bn/an+chi; - c[n]=f[n]=0; - } else { - a[n]=c[n]=0; - b[n]=1; - f[n]=angle(v[n]-v[n-1],dirn,reference); - } - } - - for(int i=i0; i < n; ++i) { - real in=beta[i-1]^2*l[i-1]; - real A=in/alpha[i-1]; - a[i]=A; - real B=3*in-A; - real out=alpha[i]^2*l[i]; - real C=out/beta[i]; - b[i]=B+3*out-C; - c[i]=C; - f[i]=-B*psi[i-1]-C*psi[i]; - } - - return tridiagonal(a,b,c,f); -} - -triple reference(triple[] v, int n, triple d0, triple d1) -{ - triple[] V=sequence(new triple(int i) { - return cross(v[i+1]-v[i],v[i+2]-v[i+1]); - },n-1); - if(n > 0) { - V.push(cross(d0,v[1]-v[0])); - V.push(cross(v[n]-v[n-1],d1)); - } - - triple max=V[0]; - real M=abs(max); - for(int i=1; i < V.length; ++i) { - triple vi=V[i]; - real a=abs(vi); - if(a > M) { - M=a; - max=vi; - } - } - - triple reference; - for(int i=0; i < V.length; ++i) { - triple u=unit(V[i]); - reference += dot(u,max) < 0 ? -u : u; - } - - return reference; -} - -// Fill in missing directions for n cyclic nodes. -void aim(flatguide3 g, int N) -{ - bool cyclic=true; - int start=0, end=0; - - // If the cycle contains one or more direction specifiers, break the loop. - for(int k=0; k < N; ++k) - if(g.solved(k)) {cyclic=false; end=k; break;} - for(int k=N-1; k >= 0; --k) - if(g.solved(k)) {cyclic=false; start=k; break;} - while(start < N && g.control[start].active) ++start; - - int n=N-(start-end); - if(n <= 1 || (cyclic && n <= 2)) return; - - triple[] v=new triple[cyclic ? n : n+1]; - real[] alpha=new real[n]; - real[] beta=new real[n]; - for(int k=0; k < n; ++k) { - int K=(start+k) % N; - v[k]=g.nodes[K]; - alpha[k]=g.Tension[K].out; - beta[k]=g.Tension[K].in; - } - if(cyclic) { - v.cyclic=true; - alpha.cyclic=true; - beta.cyclic=true; - } else v[n]=g.nodes[(start+n) % N]; - int final=(end-1) % N; - - triple d0=g.out[start].dir; - triple d1=g.in[final].dir; - - triple reference=reference(v,n,d0,d1); - - real[] theta=theta(v,alpha,beta,d0,d1,g.out[start].gamma,g.in[final].gamma, - reference); - - v.cyclic=true; - theta.cyclic=true; - - for(int k=1; k < (cyclic ? n+1 : n); ++k) { - triple w=dir(theta[k],v[k]-v[k-1],v[k+1]-v[k],reference); - g.in[(start+k-1) % N].init(w); - g.out[(start+k) % N].init(w); - } - - if(g.out[start].dir == O) - g.out[start].init(dir(theta[0],v[0]-g.nodes[(start-1) % N],v[1]-v[0], - reference)); - if(g.in[final].dir == O) - g.in[final].init(dir(theta[n],v[n-1]-v[n-2],v[n]-v[n-1],reference)); -} - -// Fill in missing directions for the sequence of nodes i...n. -void aim(flatguide3 g, int i, int n) -{ - int j=n-i; - if(j > 1 || g.out[i].dir != O || g.in[i].dir != O) { - triple[] v=new triple[j+1]; - real[] alpha=new real[j]; - real[] beta=new real[j]; - for(int k=0; k < j; ++k) { - v[k]=g.nodes[i+k]; - alpha[k]=g.Tension[i+k].out; - beta[k]=g.Tension[i+k].in; - } - v[j]=g.nodes[n]; - - triple d0=g.out[i].dir; - triple d1=g.in[n-1].dir; - - triple reference=reference(v,j,d0,d1); - - real[] theta=theta(v,alpha,beta,d0,d1,g.out[i].gamma,g.in[n-1].gamma, - reference); - - for(int k=1; k < j; ++k) { - triple w=dir(theta[k],v[k]-v[k-1],v[k+1]-v[k],reference); - g.in[i+k-1].init(w); - g.out[i+k].init(w); - } - if(g.out[i].dir == O) { - triple w=dir(theta[0],g.in[i].dir,v[1]-v[0],reference); - if(i > 0) g.in[i-1].init(w); - g.out[i].init(w); - } - if(g.in[n-1].dir == O) { - triple w=dir(theta[j],g.out[n-1].dir,v[j]-v[j-1],reference); - g.in[n-1].init(w); - g.out[n].init(w); - } - } -} - -private real Fuzz=10*realEpsilon; - -triple XYplane(pair z) {return (z.x,z.y,0);} -triple YZplane(pair z) {return (0,z.x,z.y);} -triple ZXplane(pair z) {return (z.y,0,z.x);} - -bool cyclic(guide3 g) {flatguide3 f; g(f); return f.cyclic();} -int size(guide3 g) {flatguide3 f; g(f); return f.size();} -int length(guide3 g) {flatguide3 f; g(f); return f.nodes.length-1;} - -triple dir(path3 p) -{ - return dir(p,length(p)); -} - -triple dir(path3 p, path3 h) -{ - return 0.5*(dir(p)+dir(h)); -} - -// return the point on path3 p at arclength L -triple arcpoint(path3 p, real L) -{ - return point(p,arctime(p,L)); -} - -// return the direction on path3 p at arclength L -triple arcdir(path3 p, real L) -{ - return dir(p,arctime(p,L)); -} - -// return the time on path3 p at the relative fraction l of its arclength -real reltime(path3 p, real l) -{ - return arctime(p,l*arclength(p)); -} - -// return the point on path3 p at the relative fraction l of its arclength -triple relpoint(path3 p, real l) -{ - return point(p,reltime(p,l)); -} - -// return the direction of path3 p at the relative fraction l of its arclength -triple reldir(path3 p, real l) -{ - return dir(p,reltime(p,l)); -} - -// return the initial point of path3 p -triple beginpoint(path3 p) -{ - return point(p,0); -} - -// return the point on path3 p at half of its arclength -triple midpoint(path3 p) -{ - return relpoint(p,0.5); -} - -// return the final point of path3 p -triple endpoint(path3 p) -{ - return point(p,length(p)); -} - -path3 path3(triple v) -{ - triple[] point={v}; - return path3(point,point,point,new bool[] {false},false); -} - -path3 path3(path p, triple plane(pair)=XYplane) -{ - int n=size(p); - return path3(sequence(new triple(int i) {return plane(precontrol(p,i));},n), - sequence(new triple(int i) {return plane(point(p,i));},n), - sequence(new triple(int i) {return plane(postcontrol(p,i));},n), - sequence(new bool(int i) {return straight(p,i);},n), - cyclic(p)); -} - -path3[] path3(explicit path[] g, triple plane(pair)=XYplane) -{ - return sequence(new path3(int i) {return path3(g[i],plane);},g.length); -} - -path3 invert(path p, triple normal, triple point, - projection P=currentprojection) -{ - return path3(p,new triple(pair z) {return invert(z,normal,point,P);}); -} - -path3 invert(path p, triple point, projection P=currentprojection) -{ - return path3(p,new triple(pair z) {return invert(z,P.normal,point,P);}); -} - -path3 invert(path p, projection P=currentprojection) -{ - return path3(p,new triple(pair z) {return invert(z,P.normal,P.target,P);}); -} - -// Construct a path from a path3 by applying P to each control point. -path path(path3 p, pair P(triple)=xypart) -{ - int n=length(p); - if(n < 0) return nullpath; - guide g=P(point(p,0)); - if(n == 0) return g; - for(int i=1; i < n; ++i) - g=straight(p,i-1) ? g--P(point(p,i)) : - g..controls P(postcontrol(p,i-1)) and P(precontrol(p,i))..P(point(p,i)); - - if(straight(p,n-1)) - return cyclic(p) ? g--cycle : g--P(point(p,n)); - - pair post=P(postcontrol(p,n-1)); - pair pre=P(precontrol(p,n)); - return cyclic(p) ? g..controls post and pre..cycle : - g..controls post and pre..P(point(p,n)); -} - -void write(file file, string s="", explicit path3 x, suffix suffix=none) -{ - write(file,s); - int n=length(x); - if(n < 0) write("<nullpath3>"); - else { - for(int i=0; i < n; ++i) { - write(file,point(x,i)); - if(i < length(x)) { - if(straight(x,i)) write(file,"--"); - else { - write(file,".. controls "); - write(file,postcontrol(x,i)); - write(file," and "); - write(file,precontrol(x,i+1),newl); - write(file," .."); - } - } - } - if(cyclic(x)) - write(file,"cycle",suffix); - else - write(file,point(x,n),suffix); - } -} - -void write(string s="", explicit path3 x, suffix suffix=endl) -{ - write(stdout,s,x,suffix); -} - -void write(file file, string s="", explicit path3[] x, suffix suffix=none) -{ - write(file,s); - if(x.length > 0) write(file,x[0]); - for(int i=1; i < x.length; ++i) { - write(file,endl); - write(file," ^^"); - write(file,x[i]); - } - write(file,suffix); -} - -void write(string s="", explicit path3[] x, suffix suffix=endl) -{ - write(stdout,s,x,suffix); -} - -path3 solve(flatguide3 g) -{ - int n=g.nodes.length-1; - - // If duplicate points occur consecutively, add dummy controls (if absent). - for(int i=0; i < n; ++i) { - if(g.nodes[i] == g.nodes[i+1] && !g.control[i].active) - g.control[i]=control(g.nodes[i],g.nodes[i],straight=true); - } - - // Fill in empty direction specifiers inherited from explicit control points. - for(int i=0; i < n; ++i) { - if(g.control[i].active) { - g.out[i].init(g.control[i].post-g.nodes[i]); - g.in[i].init(g.nodes[i+1]-g.control[i].pre); - } - } - - // Propagate directions across nodes. - for(int i=0; i < n; ++i) { - int next=g.cyclic[i+1] ? 0 : i+1; - if(g.out[next].active()) - g.in[i].default(g.out[next]); - if(g.in[i].active()) { - g.out[next].default(g.in[i]); - g.out[i+1].default(g.in[i]); - } - } - - // Compute missing 3D directions. - // First, resolve cycles - int i=find(g.cyclic); - if(i > 0) { - aim(g,i); - // All other cycles can now be reduced to sequences. - triple v=g.out[0].dir; - for(int j=i; j <= n; ++j) { - if(g.cyclic[j]) { - g.in[j-1].default(v); - g.out[j].default(v); - if(g.nodes[j-1] == g.nodes[j] && !g.control[j-1].active) - g.control[j-1]=control(g.nodes[j-1],g.nodes[j-1]); - } - } - } - - // Next, resolve sequences. - int i=0; - int start=0; - while(i < n) { - // Look for a missing outgoing direction. - while(i <= n && g.solved(i)) {start=i; ++i;} - if(i > n) break; - // Look for the end of the sequence. - while(i < n && !g.solved(i)) ++i; - - while(start < i && g.control[start].active) ++start; - - if(start < i) - aim(g,start,i); - } - - // Compute missing 3D control points. - for(int i=0; i < n; ++i) { - int next=g.cyclic[i+1] ? 0 : i+1; - if(!g.control[i].active) { - control c; - if((g.out[i].Curl && g.in[i].Curl) || - (g.out[i].dir == O && g.in[i].dir == O)) { - // fill in straight control points for path3 functions - triple delta=(g.nodes[i+1]-g.nodes[i])/3; - c=control(g.nodes[i]+delta,g.nodes[i+1]-delta,straight=true); - } else { - Controls C=Controls(g.nodes[i],g.nodes[next],g.out[i].dir,g.in[i].dir, - g.Tension[i].out,g.Tension[i].in, - g.Tension[i].atLeast); - c=control(C.c0,C.c1); - } - g.control[i]=c; - } - } - - // Convert to Knuth's format (control points stored with nodes) - int n=g.nodes.length; - bool cyclic; - if(n > 0) { - cyclic=g.cyclic[n-1]; - if(cyclic) --n; - } - triple[] pre=new triple[n]; - triple[] point=new triple[n]; - triple[] post=new triple[n]; - bool[] straight=new bool[n]; - if(n > 0) { - for(int i=0; i < n-1; ++i) { - point[i]=g.nodes[i]; - post[i]=g.control[i].post; - pre[i+1]=g.control[i].pre; - straight[i]=g.control[i].straight; - } - point[n-1]=g.nodes[n-1]; - if(cyclic) { - pre[0]=g.control[n-1].pre; - post[n-1]=g.control[n-1].post; - straight[n-1]=g.control[n-1].straight; - } else { - pre[0]=point[0]; - post[n-1]=point[n-1]; - straight[n-1]=false; - } - } - - return path3(pre,point,post,straight,cyclic); -} - -path nurb(path3 p, projection P, int ninterpolate=P.ninterpolate) -{ - triple f=P.camera; - triple u=unit(P.normal); - transform3 t=P.t; - - path nurb(triple v0, triple v1, triple v2, triple v3) { - return nurb(project(v0,t),project(v1,t),project(v2,t),project(v3,t), - dot(u,f-v0),dot(u,f-v1),dot(u,f-v2),dot(u,f-v3),ninterpolate); - } - - path g; - - if(straight(p,0)) - g=project(point(p,0),t); - - int last=length(p); - for(int i=0; i < last; ++i) { - if(straight(p,i)) - g=g--project(point(p,i+1),t); - else - g=g&nurb(point(p,i),postcontrol(p,i),precontrol(p,i+1),point(p,i+1)); - } - - int n=length(g); - if(cyclic(p)) g=g&cycle; - - return g; -} - -path project(path3 p, projection P=currentprojection, - int ninterpolate=P.ninterpolate) -{ - guide g; - - int last=length(p); - if(last < 0) return g; - - transform3 t=P.t; - - if(ninterpolate == 1 || piecewisestraight(p)) { - g=project(point(p,0),t); - // Construct the path. - int stop=cyclic(p) ? last-1 : last; - for(int i=0; i < stop; ++i) { - if(straight(p,i)) - g=g--project(point(p,i+1),t); - else { - g=g..controls project(postcontrol(p,i),t) and - project(precontrol(p,i+1),t)..project(point(p,i+1),t); - } - } - } else return nurb(p,P); - - if(cyclic(p)) - g=straight(p,last-1) ? g--cycle : - g..controls project(postcontrol(p,last-1),t) and - project(precontrol(p,last),t)..cycle; - return g; -} - -pair[] project(triple[] v, projection P=currentprojection) -{ - return sequence(new pair(int i) {return project(v[i],P.t);},v.length); -} - -path[] project(explicit path3[] g, projection P=currentprojection) -{ - return sequence(new path(int i) {return project(g[i],P);},g.length); -} - -guide3 operator cast(path3 p) -{ - int last=length(p); - - bool cyclic=cyclic(p); - int stop=cyclic ? last-1 : last; - return new void(flatguide3 f) { - if(last >= 0) { - f.node(point(p,0)); - for(int i=0; i < stop; ++i) { - if(straight(p,i)) { - f.out(1); - f.in(1); - } else - f.control(postcontrol(p,i),precontrol(p,i+1)); - f.node(point(p,i+1)); - } - if(cyclic) { - if(straight(p,stop)) { - f.out(1); - f.in(1); - } else - f.control(postcontrol(p,stop),precontrol(p,last)); - f.cycleToken(); - } - } - }; -} - -// Return a unit normal vector to a planar path p (or O if the path is -// nonplanar). -triple normal(path3 p) -{ - triple normal; - real fuzz=sqrtEpsilon*abs(max(p)-min(p)); - real absnormal; - real theta; - - bool Cross(triple a, triple b) { - if(abs(a) >= fuzz && abs(b) >= fuzz) { - triple n=cross(unit(a),unit(b)); - real absn=abs(n); - if(absn < sqrtEpsilon) return false; - n=unit(n); - if(absnormal > 0 && - abs(normal-n) > sqrtEpsilon && abs(normal+n) > sqrtEpsilon) - return true; - else { - int sign=dot(n,normal) >= 0 ? 1 : -1; - theta += sign*asin1(absn); - if(absn > absnormal) { - absnormal=absn; - normal=n; - theta=sign*theta; - } - } - } - return false; - } - - int L=length(p); - if(L <= 0) return O; - - triple zi=point(p,0); - triple v0=zi-precontrol(p,0); - for(int i=0; i < L; ++i) { - triple c0=postcontrol(p,i); - triple c1=precontrol(p,i+1); - triple zp=point(p,i+1); - triple v1=c0-zi; - triple v2=c1-c0; - triple v3=zp-c1; - if(Cross(v0,v1) || Cross(v1,v2) || Cross(v2,v3)) return O; - v0=v3; - zi=zp; - } - return theta >= 0 ? normal : -normal; -} - -// Return a unit normal vector to a polygon with vertices in p. -triple normal(triple[] p) -{ - triple normal; - real fuzz=sqrtEpsilon*abs(maxbound(p)-minbound(p)); - real absnormal; - real theta; - - bool Cross(triple a, triple b) { - if(abs(a) >= fuzz && abs(b) >= fuzz) { - triple n=cross(unit(a),unit(b)); - real absn=abs(n); - n=unit(n); - if(absnormal > 0 && absn > sqrtEpsilon && - abs(normal-n) > sqrtEpsilon && abs(normal+n) > sqrtEpsilon) - return true; - else { - int sign=dot(n,normal) >= 0 ? 1 : -1; - theta += sign*asin1(absn); - if(absn > absnormal) { - absnormal=absn; - normal=n; - theta=sign*theta; - } - } - } - return false; - } - - if(p.length <= 0) return O; - - triple zi=p[0]; - triple v0=zi-p[p.length-1]; - for(int i=0; i < p.length-1; ++i) { - triple zp=p[i+1]; - triple v1=zp-zi; - if(Cross(v0,v1)) return O; - v0=v1; - zi=zp; - } - return theta >= 0 ? normal : -normal; -} - -// Transforms that map XY plane to YX, YZ, ZY, ZX, and XZ planes. -restricted transform3 XY=identity4; -restricted transform3 YX=rotate(-90,O,Z); -restricted transform3 YZ=rotate(90,O,Z)*rotate(90,O,X); -restricted transform3 ZY=rotate(-90,O,X)*YZ; -restricted transform3 ZX=rotate(-90,O,Z)*rotate(-90,O,Y); -restricted transform3 XZ=rotate(-90,O,Y)*ZX; - -private transform3 flip(transform3 t, triple X, triple Y, triple Z, - projection P) -{ - static transform3 flip(triple v) { - static real s(real x) {return x > 0 ? -1 : 1;} - return scale(s(v.x),s(v.y),s(v.z)); - } - - triple u=unit(P.normal); - triple up=unit(perp(P.up,u)); - bool upright=dot(Z,u) >= 0; - if(dot(Y,up) < 0) { - t=flip(Y)*t; - upright=!upright; - } - return upright ? t : flip(X)*t; -} - -restricted transform3 XY(projection P=currentprojection) -{ - return flip(XY,X,Y,Z,P); -} - -restricted transform3 YX(projection P=currentprojection) -{ - return flip(YX,Y,X,Z,P); -} - -restricted transform3 YZ(projection P=currentprojection) -{ - return flip(YZ,Y,Z,X,P); -} - -restricted transform3 ZY(projection P=currentprojection) -{ - return flip(ZY,Z,Y,X,P); -} - -restricted transform3 ZX(projection P=currentprojection) -{ - return flip(ZX,Z,X,Y,P); -} - -restricted transform3 XZ(projection P=currentprojection) -{ - return flip(XZ,X,Z,Y,P); -} - -// Transform3 that projects in direction dir onto plane with normal n -// through point O. -transform3 planeproject(triple n, triple O=O, triple dir=n) -{ - real a=n.x, b=n.y, c=n.z; - real u=dir.x, v=dir.y, w=dir.z; - real delta=1.0/(a*u+b*v+c*w); - real d=-(a*O.x+b*O.y+c*O.z)*delta; - return new real[][] { - {(b*v+c*w)*delta,-b*u*delta,-c*u*delta,-d*u}, - {-a*v*delta,(a*u+c*w)*delta,-c*v*delta,-d*v}, - {-a*w*delta,-b*w*delta,(a*u+b*v)*delta,-d*w}, - {0,0,0,1} - }; -} - -// Transform3 that projects in direction dir onto plane defined by p. -transform3 planeproject(path3 p, triple dir=O) -{ - triple n=normal(p); - return planeproject(n,point(p,0),dir == O ? n : dir); -} - -// Transform for projecting onto plane through point O with normal cross(u,v). -transform transform(triple u, triple v, triple O=O, - projection P=currentprojection) -{ - transform3 t=P.t; - static real[] O={0,0,0,1}; - real[] tO=t*O; - real tO3=tO[3]; - real factor=1/tO3^2; - real[] x=(tO3*t[0]-tO[0]*t[3])*factor; - real[] y=(tO3*t[1]-tO[1]*t[3])*factor; - triple x=(x[0],x[1],x[2]); - triple y=(y[0],y[1],y[2]); - u=unit(u); - v=unit(v); - return (0,0,dot(u,x),dot(v,x),dot(u,y),dot(v,y)); -} - -// Project Label onto plane through point O with normal cross(u,v). -Label project(Label L, triple u, triple v, triple O=O, - projection P=currentprojection) { - Label L=L.copy(); - L.position=project(O,P.t); - L.transform(transform(u,v,O,P)); - return L; -} - -path3 operator cast(guide3 g) {return solve(g);} -path3 operator cast(triple v) {return path3(v);} - -guide3[] operator cast(triple[] v) -{ - return sequence(new guide3(int i) {return v[i];},v.length); -} - -path3[] operator cast(triple[] v) -{ - return sequence(new path3(int i) {return v[i];},v.length); -} - -path3[] operator cast(guide3[] g) -{ - return sequence(new path3(int i) {return solve(g[i]);},g.length); -} - -guide3[] operator cast(path3[] g) -{ - return sequence(new guide3(int i) {return g[i];},g.length); -} - -void write(file file, string s="", explicit guide3[] x, suffix suffix=none) -{ - write(file,s,(path3[]) x,suffix); -} - -void write(string s="", explicit guide3[] x, suffix suffix=endl) -{ - write(stdout,s,(path3[]) x,suffix); -} - -triple point(explicit guide3 g, int t) { - flatguide3 f; - g(f); - int n=f.size(); - return f.nodes[adjustedIndex(t,n,f.cyclic())]; -} - -triple[] dirSpecifier(guide3 g, int t) -{ - flatguide3 f; - g(f); - int n=f.size(); - checkEmpty(n); - if(f.cyclic()) t=t % n; - else if(t < 0 || t >= n-1) return new triple[]; - return new triple[] {f.out[t].dir,f.in[t].dir}; -} - -triple[] controlSpecifier(guide3 g, int t) { - flatguide3 f; - g(f); - int n=f.size(); - checkEmpty(n); - if(f.cyclic()) t=t % n; - else if(t < 0 || t >= n-1) return new triple[]; - control c=f.control[t]; - if(c.active) return new triple[] {c.post,c.pre}; - else return new triple[]; -} - -tensionSpecifier tensionSpecifier(guide3 g, int t) -{ - flatguide3 f; - g(f); - int n=f.size(); - checkEmpty(n); - if(f.cyclic()) t=t % n; - else if(t < 0 || t >= n-1) return operator tension(1,1,false); - Tension T=f.Tension[t]; - return operator tension(T.out,T.in,T.atLeast); -} - -real[] curlSpecifier(guide3 g, int t) -{ - flatguide3 f; - g(f); - int n=f.size(); - checkEmpty(n); - if(f.cyclic()) t=t % n; - else if(t < 0 || t >= n-1) return new real[]; - return new real[] {f.out[t].gamma,f.in[t].gamma}; -} - -guide3 reverse(guide3 g) -{ - flatguide3 f; - bool cyclic=cyclic(g); - g(f); - - if(f.precyclic()) - return reverse(solve(g)); - - int n=f.size(); - checkEmpty(n); - guide3 G; - if(n >= 0) { - int start=cyclic ? n : n-1; - for(int i=start; i > 0; --i) { - G=G..f.nodes[i]; - control c=f.control[i-1]; - if(c.active) - G=G..operator controls(c.pre,c.post); - else { - dir in=f.in[i-1]; - triple d=in.dir; - if(d != O) G=G..operator spec(-d,JOIN_OUT); - else if(in.Curl) G=G..operator curl(in.gamma,JOIN_OUT); - dir out=f.out[i-1]; - triple d=out.dir; - if(d != O) G=G..operator spec(-d,JOIN_IN); - else if(out.Curl) G=G..operator curl(out.gamma,JOIN_IN); - } - } - if(cyclic) G=G..cycle; - else G=G..f.nodes[0]; - } - return G; -} - -triple intersectionpoint(path3 p, path3 q, real fuzz=-1) -{ - real[] t=intersect(p,q,fuzz); - if(t.length == 0) abort("paths do not intersect"); - return point(p,t[0]); -} - -// return an array containing all intersection points of p and q -triple[] intersectionpoints(path3 p, path3 q, real fuzz=-1) -{ - real[][] t=intersections(p,q,fuzz); - return sequence(new triple(int i) {return point(p,t[i][0]);},t.length); -} - -triple[] intersectionpoints(explicit path3[] p, explicit path3[] q, - real fuzz=-1) -{ - triple[] v; - for(int i=0; i < p.length; ++i) - for(int j=0; j < q.length; ++j) - v.append(intersectionpoints(p[i],q[j],fuzz)); - return v; -} - -path3 operator &(path3 p, cycleToken tok) -{ - int n=length(p); - if(n < 0) return nullpath3; - triple a=point(p,0); - triple b=point(p,n); - return subpath(p,0,n-1)..controls postcontrol(p,n-1) and precontrol(p,n).. - cycle; -} - -// return the point on path3 p at arclength L -triple arcpoint(path3 p, real L) -{ - return point(p,arctime(p,L)); -} - -// return the point on path3 p at arclength L -triple arcpoint(path3 p, real L) -{ - return point(p,arctime(p,L)); -} - -// return the direction on path3 p at arclength L -triple arcdir(path3 p, real L) -{ - return dir(p,arctime(p,L)); -} - -// return the time on path3 p at the relative fraction l of its arclength -real reltime(path3 p, real l) -{ - return arctime(p,l*arclength(p)); -} - -// return the point on path3 p at the relative fraction l of its arclength -triple relpoint(path3 p, real l) -{ - return point(p,reltime(p,l)); -} - -// return the direction of path3 p at the relative fraction l of its arclength -triple reldir(path3 p, real l) -{ - return dir(p,reltime(p,l)); -} - -// return the point on path3 p at half of its arclength -triple midpoint(path3 p) -{ - return relpoint(p,0.5); -} - -real relative(Label L, path3 g) -{ - return L.position.relative ? reltime(g,L.relative()) : L.relative(); -} - -// return the linear transformation that maps X,Y,Z to u,v,w. -transform3 transform3(triple u, triple v, triple w=cross(u,v)) -{ - return new real[][] { - {u.x,v.x,w.x,0}, - {u.y,v.y,w.y,0}, - {u.z,v.z,w.z,0}, - {0,0,0,1} - }; -} - -// return the rotation that maps Z to a unit vector u about cross(u,Z), -transform3 align(triple u) -{ - real a=u.x; - real b=u.y; - real c=u.z; - real d=a^2+b^2; - - if(d != 0) { - d=sqrt(d); - real e=1/d; - return new real[][] { - {-b*e,-a*c*e,a,0}, - {a*e,-b*c*e,b,0}, - {0,d,c,0}, - {0,0,0,1}}; - } - return c >= 0 ? identity(4) : diagonal(1,-1,-1,1); -} - -// return a rotation that maps X,Y to the projection plane. -transform3 transform3(projection P=currentprojection) -{ - triple w=unit(P.normal); - triple v=unit(perp(P.up,w)); - if(v == O) v=cross(perp(w),w); - triple u=cross(v,w); - return u != O ? transform3(u,v,w) : identity(4); -} - -triple[] triples(real[] x, real[] y, real[] z) -{ - if(x.length != y.length || x.length != z.length) - abort("arrays have different lengths"); - return sequence(new triple(int i) {return (x[i],y[i],z[i]);},x.length); -} - -path3[] operator cast(path3 p) -{ - return new path3[] {p}; -} - -path3[] operator cast(guide3 g) -{ - return new path3[] {(path3) g}; -} - -path3[] operator ^^ (path3 p, path3 q) -{ - return new path3[] {p,q}; -} - -path3[] operator ^^ (path3 p, explicit path3[] q) -{ - return concat(new path3[] {p},q); -} - -path3[] operator ^^ (explicit path3[] p, path3 q) -{ - return concat(p,new path3[] {q}); -} - -path3[] operator ^^ (explicit path3[] p, explicit path3[] q) -{ - return concat(p,q); -} - -path3[] operator * (transform3 t, explicit path3[] p) -{ - return sequence(new path3(int i) {return t*p[i];},p.length); -} - -triple[] operator * (transform3 t, triple[] v) -{ - return sequence(new triple(int i) {return t*v[i];},v.length); -} - -triple[][] operator * (transform3 t, triple[][] v) -{ - triple[][] V=new triple[v.length][]; - for(int i=0; i < v.length; ++i) { - triple[] vi=v[i]; - V[i]=sequence(new triple(int j) {return t*vi[j];},vi.length); - } - return V; -} - -triple min(explicit path3[] p) -{ - checkEmpty(p.length); - triple minp=min(p[0]); - for(int i=1; i < p.length; ++i) - minp=minbound(minp,min(p[i])); - return minp; -} - -triple max(explicit path3[] p) -{ - checkEmpty(p.length); - triple maxp=max(p[0]); - for(int i=1; i < p.length; ++i) - maxp=maxbound(maxp,max(p[i])); - return maxp; -} - -path3 randompath3(int n, bool cumulate=true, interpolate3 join=operator ..) -{ - guide3 g; - triple w; - for(int i=0; i <= n; ++i) { - triple z=(unitrand()-0.5,unitrand()-0.5,unitrand()-0.5); - if(cumulate) w += z; - else w=z; - g=join(g,w); - } - return g; -} - -path3[] box(triple v1, triple v2) -{ - return - (v1.x,v1.y,v1.z)-- - (v1.x,v1.y,v2.z)-- - (v1.x,v2.y,v2.z)-- - (v1.x,v2.y,v1.z)-- - (v1.x,v1.y,v1.z)-- - (v2.x,v1.y,v1.z)-- - (v2.x,v1.y,v2.z)-- - (v2.x,v2.y,v2.z)-- - (v2.x,v2.y,v1.z)-- - (v2.x,v1.y,v1.z)^^ - (v2.x,v2.y,v1.z)-- - (v1.x,v2.y,v1.z)^^ - (v1.x,v2.y,v2.z)-- - (v2.x,v2.y,v2.z)^^ - (v2.x,v1.y,v2.z)-- - (v1.x,v1.y,v2.z); -} - -restricted path3[] unitbox=box(O,(1,1,1)); -restricted path3 unitcircle3=X..Y..-X..-Y..cycle; -restricted path3 unitsquare3=O--X--X+Y--Y--cycle; - -path3 circle(triple c, real r, triple normal=Z) -{ - path3 p=normal == Z ? unitcircle3 : align(unit(normal))*unitcircle3; - return shift(c)*scale3(r)*p; -} - -// return an arc centered at c from triple v1 to v2 (assuming |v2-c|=|v1-c|), -// drawing in the given direction. -// The normal must be explicitly specified if c and the endpoints are colinear. -path3 arc(triple c, triple v1, triple v2, triple normal=O, bool direction=CCW) -{ - v1 -= c; - real r=abs(v1); - v1=unit(v1); - v2=unit(v2-c); - - if(normal == O) { - normal=cross(v1,v2); - if(normal == O) abort("explicit normal required for these endpoints"); - } - - transform3 T; - bool align=normal != Z; - if(align) { - T=align(unit(normal)); - transform3 Tinv=transpose(T); - v1=Tinv*v1; - v2=Tinv*v2; - } - - string invalidnormal="invalid normal vector"; - real fuzz=sqrtEpsilon*max(abs(v1),abs(v2)); - if(abs(v1.z) > fuzz || abs(v2.z) > fuzz) - abort(invalidnormal); - - real[] t1=intersect(unitcircle3,O--2*(v1.x,v1.y,0)); - real[] t2=intersect(unitcircle3,O--2*(v2.x,v2.y,0)); - - if(t1.length == 0 || t2.length == 0) - abort(invalidnormal); - - real t1=t1[0]; - real t2=t2[0]; - int n=length(unitcircle3); - if(direction) { - if (t1 >= t2) t1 -= n; - } else if(t2 >= t1) t2 -= n; - - path3 p=subpath(unitcircle3,t1,t2); - if(align) p=T*p; - return shift(c)*scale3(r)*p; -} - -// return an arc centered at c with radius r from c+r*dir(theta1,phi1) to -// c+r*dir(theta2,phi2) in degrees, drawing in the given direction -// relative to the normal vector cross(dir(theta1,phi1),dir(theta2,phi2)). -// The normal must be explicitly specified if c and the endpoints are colinear. -path3 arc(triple c, real r, real theta1, real phi1, real theta2, real phi2, - triple normal=O, bool direction) -{ - return arc(c,c+r*dir(theta1,phi1),c+r*dir(theta2,phi2),normal,direction); -} - -// return an arc centered at c with radius r from c+r*dir(theta1,phi1) to -// c+r*dir(theta2,phi2) in degrees, drawing drawing counterclockwise -// relative to the normal vector cross(dir(theta1,phi1),dir(theta2,phi2)) -// iff theta2 > theta1 or (theta2 == theta1 and phi2 >= phi1). -// The normal must be explicitly specified if c and the endpoints are colinear. -path3 arc(triple c, real r, real theta1, real phi1, real theta2, real phi2, - triple normal=O) -{ - return arc(c,r,theta1,phi1,theta2,phi2,normal, - theta2 > theta1 || (theta2 == theta1 && phi2 >= phi1) ? CCW : CW); -} - -private real epsilon=1000*realEpsilon; - -// Return a representation of the plane through point O with normal cross(u,v). -path3 plane(triple u, triple v, triple O=O) -{ - return O--O+u--O+u+v--O+v--cycle; -} - -// PRC/OpenGL support - -include three_light; - -void draw(frame f, path3 g, material p=currentpen, light light=nolight, - string name="", render render=defaultrender, - projection P=currentprojection); - -void begingroup3(frame f, string name="", render render=defaultrender, - triple center=O, int interaction=0) -{ - _begingroup3(f,name,render.compression,render.granularity,render.closed, - render.tessellate,render.merge == false, - render.merge == true,center,interaction); -} - -void begingroup3(picture pic=currentpicture, string name="", - render render=defaultrender, - triple center=O, int interaction=0) -{ - pic.add(new void(frame f, transform3, picture pic, projection) { - if(is3D()) - begingroup3(f,name,render,center,interaction); - if(pic != null) - begingroup(pic); - },true); -} - -void endgroup3(picture pic=currentpicture) -{ - pic.add(new void(frame f, transform3, picture pic, projection) { - if(is3D()) - endgroup3(f); - if(pic != null) - endgroup(pic); - },true); -} - -void addPath(picture pic, path3 g, pen p) -{ - if(size(g) > 0) - pic.addBox(min(g),max(g),min3(p),max3(p)); -} - -include three_surface; -include three_margins; - -pair min(path3 p, projection P) -{ - path3 q=P.T.modelview*p; - if(P.infinity) - return xypart(min(q)); - return maxratio(q)/P.T.projection[3][2]; -} - -pair max(path3 p, projection P) -{ - path3 q=P.T.modelview*p; - if(P.infinity) - return xypart(max(q)); - return minratio(q)/P.T.projection[3][2]; -} - -pair min(frame f, projection P) -{ - frame g=P.T.modelview*f; - if(P.infinity) - return xypart(min3(g)); - return maxratio(g)/P.T.projection[3][2]; -} - -pair max(frame f, projection P) -{ - frame g=P.T.modelview*f; - if(P.infinity) - return xypart(max3(g)); - return minratio(g)/P.T.projection[3][2]; -} - -void draw(picture pic=currentpicture, Label L="", path3 g, - align align=NoAlign, material p=currentpen, margin3 margin=NoMargin3, - light light=nolight, string name="", render render=defaultrender) -{ - pen q=(pen) p; - pic.add(new void(frame f, transform3 t, picture pic, projection P) { - path3 G=margin(t*g,q).g; - if(is3D()) { - draw(f,G,p,light,name,render,null); - if(pic != null && size(G) > 0) - pic.addBox(min(G,P),max(G,P),min(q),max(q)); - } - if(pic != null) - draw(pic,project(G,P),q); - },true); - Label L=L.copy(); - L.align(align); - if(L.s != "") { - L.p(q); - label(pic,L,g); - } - addPath(pic,g,q); -} - -include three_tube; - -draw=new void(frame f, path3 g, material p=currentpen, - light light=nolight, string name="", - render render=defaultrender, - projection P=currentprojection) { - pen q=(pen) p; - if(is3D()) { - p=material(p); - real width=linewidth(q); - void drawthick(path3 g) { - if(settings.thick) { - if(width > 0) { - bool prc=prc(); - void cylinder(transform3) {}; - void sphere(transform3, bool half) {}; - void disk(transform3) {}; - void pipe(path3, path3); - if(prc) { - cylinder=new void(transform3 t) {drawPRCcylinder(f,t,p,light);}; - sphere=new void(transform3 t, bool half) - {drawPRCsphere(f,t,half,p,light,render);}; - disk=new void(transform3 t) {draw(f,t*unitdisk,p,light,render);}; - pipe=new void(path3 center, path3 g) - {drawPRCtube(f,center,g,p,light);}; - } - real linecap=linecap(q); - real r=0.5*width; - bool open=!cyclic(g); - int L=length(g); - triple g0=point(g,0); - triple gL=point(g,L); - if(open && L > 0) { - if(linecap == 2) { - g0 -= r*dir(g,0); - gL += r*dir(g,L); - g=g0..g..gL; - L += 2; - } - } - tube T=tube(g,width,render,cylinder,sphere,pipe); - path3 c=T.center; - if(L >= 0) { - if(open) { - int Lc=length(c); - triple c0=point(c,0); - triple cL=point(c,Lc); - triple dir0=dir(g,0); - triple dirL=dir(g,L); - triple dirc0=dir(c,0); - triple dircL=dir(c,Lc); - transform3 t0=shift(g0)*align(-dir0); - transform3 tL=shift(gL)*align(dirL); - transform3 tc0=shift(c0)*align(-dirc0); - transform3 tcL=shift(cL)*align(dircL); - if(linecap == 0 || linecap == 2){ - transform3 scale2r=scale(r,r,1); - T.s.append(t0*scale2r*unitdisk); - disk(tc0*scale2r); - if(L > 0) { - T.s.append(tL*scale2r*unitdisk); - disk(tcL*scale2r); - } - } else if(linecap == 1) { - transform3 scale3r=scale3(r); - T.s.append(t0*scale3r* - (dir0 != O ? unithemisphere : unitsphere)); - sphere(tc0*scale3r,half=straight(c,0)); - if(L > 0) { - T.s.append(tL*scale3r* - (dirL != O ? unithemisphere : unitsphere)); - sphere(tcL*scale3r,half=straight(c,Lc-1)); - } - } - } - if(opacity(q) == 1) - _draw(f,c,q); - } - for(patch s : T.s.s) - draw3D(f,s,p,light,prc=false); - } else _draw(f,g,q); - } else _draw(f,g,q); - } - if(q != nullpen) - begingroup3(f,name == "" ? "curve" : name,render); - if(linetype(q).length == 0) drawthick(g); - else { - real[] dash=linetype(adjust(q,arclength(g),cyclic(g))); - if(sum(dash) > 0) { - dash.cyclic=true; - real offset=offset(q); - real L=arclength(g); - int i=0; - real l=offset; - while(l <= L) { - real t1=arctime(g,l); - l += dash[i]; - real t2=arctime(g,min(l,L)); - drawthick(subpath(g,t1,t2)); - ++i; - l += dash[i]; - ++i; - } - } - } - if(q != nullpen) - endgroup3(f); - } else draw(f,project(g,P),q); -}; - -void draw(frame f, explicit path3[] g, material p=currentpen, - light light=nolight, string name="", - render render=defaultrender, projection P=currentprojection) -{ - if(g.length > 1) - begingroup3(f,name == "" ? "curve" : name,render); - for(int i=0; i < g.length; ++i) - draw(f,g[i],p,light,partname(i),render,P); - if(g.length > 1) - endgroup3(f); -} - -void draw(picture pic=currentpicture, explicit path3[] g, - material p=currentpen, margin3 margin=NoMargin3, light light=nolight, - string name="", render render=defaultrender) -{ - if(g.length > 1) - begingroup3(pic,name == "" ? "curves" : name,render); - for(int i=0; i < g.length; ++i) - draw(pic,g[i],p,margin,light,partname(i),render); - if(g.length > 1) - endgroup3(pic); -} - -include three_arrows; - -void draw(picture pic=currentpicture, Label L="", path3 g, - align align=NoAlign, material p=currentpen, arrowbar3 arrow, - arrowbar3 bar=None, margin3 margin=NoMargin3, light light=nolight, - light arrowheadlight=currentlight, string name="", - render render=defaultrender) -{ - bool group=arrow != None || bar != None; - if(group) - begingroup3(pic,name,render); - bool drawpath=arrow(pic,g,p,margin,light,arrowheadlight); - if(bar(pic,g,p,margin,light,arrowheadlight) && drawpath) - draw(pic,L,g,align,p,margin,light,render); - if(group) - endgroup3(pic); - if(L.s != "") - label(pic,L,g,align,(pen) p); -} - -void draw(frame f, path3 g, material p=currentpen, arrowbar3 arrow, - light light=nolight, light arrowheadlight=currentlight, - string name="", render render=defaultrender, - projection P=currentprojection) -{ - picture pic; - bool group=arrow != None; - if(group) - begingroup3(f,name,render); - if(arrow(pic,g,p,NoMargin3,light,arrowheadlight)) - draw(f,g,p,light,render,P); - add(f,pic.fit()); - if(group) - endgroup3(f); -} - -void add(picture pic=currentpicture, void d(picture,transform3), - bool exact=false) -{ - pic.add(d,exact); -} - -// Fit the picture src using the identity transformation (so user -// coordinates and truesize coordinates agree) and add it about the point -// position to picture dest. -void add(picture dest, picture src, triple position, bool group=true, - bool above=true) -{ - dest.add(new void(picture f, transform3 t) { - f.add(shift(t*position)*src,group,above); - }); -} - -void add(picture src, triple position, bool group=true, bool above=true) -{ - add(currentpicture,src,position,group,above); -} - -// Align an arrow pointing to b from the direction dir. The arrow is -// 'length' PostScript units long. -void arrow(picture pic=currentpicture, Label L="", triple b, triple dir, - real length=arrowlength, align align=NoAlign, - pen p=currentpen, arrowbar3 arrow=Arrow3, margin3 margin=EndMargin3, - light light=nolight, light arrowheadlight=currentlight, - string name="", render render=defaultrender) -{ - Label L=L.copy(); - if(L.defaultposition) L.position(0); - L.align(L.align,dir); - L.p(p); - picture opic; - marginT3 margin=margin(b--b,p); // Extract margin.begin and margin.end - triple a=(margin.begin+length+margin.end)*unit(dir); - draw(opic,L,a--O,align,p,arrow,margin,light,arrowheadlight,name,render); - add(pic,opic,b); -} - -void arrow(picture pic=currentpicture, Label L="", triple b, pair dir, - real length=arrowlength, align align=NoAlign, - pen p=currentpen, arrowbar3 arrow=Arrow3, margin3 margin=EndMargin3, - light light=nolight, light arrowheadlight=currentlight, - string name="", render render=defaultrender, - projection P=currentprojection) -{ - arrow(pic,L,b,invert(dir,b,P),length,align,p,arrow,margin,light, - arrowheadlight,name,render); -} - -triple min3(picture pic, projection P=currentprojection) -{ - return pic.min3(P); -} - -triple max3(picture pic, projection P=currentprojection) -{ - return pic.max3(P); -} - -triple size3(picture pic, bool user=false, projection P=currentprojection) -{ - transform3 t=pic.calculateTransform3(P); - triple M=pic.max(t); - triple m=pic.min(t); - if(!user) return M-m; - t=inverse(t); - return t*M-t*m; -} - -triple point(frame f, triple dir) -{ - triple m=min3(f); - triple M=max3(f); - return m+realmult(rectify(dir),M-m); -} - -triple point(picture pic=currentpicture, triple dir, bool user=true, - projection P=currentprojection) -{ - triple min = pic.userMin(), max = pic.userMax(); - triple v=min+realmult(rectify(dir),max-min); - return user ? v : pic.calculateTransform3(P)*v; -} - -triple truepoint(picture pic=currentpicture, triple dir, bool user=true, - projection P=currentprojection) -{ - transform3 t=pic.calculateTransform3(P); - triple m=pic.min(t); - triple M=pic.max(t); - triple v=m+realmult(rectify(dir),M-m); - return user ? inverse(t)*v : v; -} - -void add(picture dest=currentpicture, object src, pair position=0, pair align=0, - bool group=true, filltype filltype=NoFill, bool above=true) -{ - if(prc()) - label(dest,src,position,align); - else if(settings.render == 0) - plain.add(dest,src,position,align,group,filltype,above); -} - -string cameralink(string label, string text="View Parameters") -{ - if(!prc() || Link == null) return ""; - return Link(label,text,"3Dgetview"); -} - -private struct viewpoint { - triple target,camera,up; - real angle; - void operator init(string s) { - s=replace(s,'\n'," "); - string[] S=split(s); - int pos(string s, string key) { - int pos=find(s,key); - if(pos < 0) return -1; - pos += length(key); - while(substr(s,pos,1) == " ") ++pos; - if(substr(s,pos,1) == "=") - return pos+1; - return -1; - } - triple C2C=X; - real ROO=1; - real ROLL=0; - angle=30; - int pos; - for(int k=0; k < S.length; ++k) { - if((pos=pos(S[k],"COO")) >= 0) - target=((real) substr(S[k],pos),(real) S[++k],(real) S[++k]); - else if((pos=pos(S[k],"C2C")) >= 0) - C2C=((real) substr(S[k],pos),(real) S[++k],(real) S[++k]); - else if((pos=pos(S[k],"ROO")) >= 0) - ROO=(real) substr(S[k],pos); - else if((pos=pos(S[k],"ROLL")) >= 0) - ROLL=(real) substr(S[k],pos); - else if((pos=pos(S[k],"AAC")) >= 0) - angle=(real) substr(S[k],pos); - } - camera=target+ROO*C2C; - triple u=unit(target-camera); - triple w=unit(Z-u.z*u); - up=rotate(ROLL,O,u)*w; - } -} - -projection perspective(string s) -{ - viewpoint v=viewpoint(s); - projection P=perspective(v.camera,v.up,v.target); - P.angle=v.angle; - P.absolute=true; - return P; -} - -private string Format(real x) -{ - // Work around movie15.sty division by zero bug; - // e.g. u=unit((1e-10,1e-10,0.9)); - if(abs(x) < 1e-9) x=0; - assert(abs(x) < 1e17,"Number too large: "+string(x)); - return format("%.18f",x,"C"); -} - -private string Format(triple v, string sep=" ") -{ - return Format(v.x)+sep+Format(v.y)+sep+Format(v.z); -} - -private string Format(real[] c) -{ - return Format((c[0],c[1],c[2])); -} - -private string format(triple v, string sep=" ") -{ - return string(v.x)+sep+string(v.y)+sep+string(v.z); -} - -private string projection(bool infinity, real viewplanesize) -{ - return "activeCamera=scene.cameras.getByIndex(0); -function asyProjection() {"+ - (infinity ? "activeCamera.projectionType=activeCamera.TYPE_ORTHOGRAPHIC;" : - "activeCamera.projectionType=activeCamera.TYPE_PERSPECTIVE;")+" -activeCamera.viewPlaneSize="+string(viewplanesize)+"; -activeCamera.binding=activeCamera.BINDING_"+(infinity ? "MAX" : "VERTICAL")+"; -} - -asyProjection(); - -handler=new CameraEventHandler(); -runtime.addEventHandler(handler); -handler.onEvent=function(event) -{ - asyProjection(); - scene.update(); -}"; -} - -private string projection(bool infinity, real viewplanesize) -{ - return "activeCamera=scene.cameras.getByIndex(0); -function asyProjection() {"+ - (infinity ? "activeCamera.projectionType=activeCamera.TYPE_ORTHOGRAPHIC;" : - "activeCamera.projectionType=activeCamera.TYPE_PERSPECTIVE;")+" -activeCamera.viewPlaneSize="+string(viewplanesize)+"; -activeCamera.binding=activeCamera.BINDING_"+(infinity ? "MAX" : "VERTICAL")+"; -} - -asyProjection(); - -handler=new CameraEventHandler(); -runtime.addEventHandler(handler); -handler.onEvent=function(event) -{ - asyProjection(); - scene.update(); -}"; -} - -private string billboard(int[] index, triple[] center) -{ - if(index.length == 0) return ""; - string s=" -var zero=new Vector3(0,0,0); -var nodes=scene.nodes; -var count=nodes.count; - -var index=new Array(); -for(i=0; i < count; i++) { - var node=nodes.getByIndex(i); - var name=node.name; - end=name.lastIndexOf(\".\")-1; - if(end > 0) { - if(name.substr(end,1) == \"\001\") { - start=name.lastIndexOf(\"-\")+1; - n=end-start; - if(n > 0) { - index[name.substr(start,n)]=i; - node.name=name.substr(0,start-1); - } - } - } -} - -var center=new Array( -"; - for(int i=0; i < center.length; ++i) - s += "Vector3("+format(center[i],",")+"), -"; - s += "); - -billboardHandler=new RenderEventHandler(); -billboardHandler.onEvent=function(event) -{ - var camera=scene.cameras.getByIndex(0); - var position=camera.position; - var direction=position.subtract(camera.targetPosition); - var up=camera.up.subtract(position); - - function f(i,k) { - j=index[i]; - if(j >= 0) { - var node=nodes.getByIndex(j); - var name=node.name; - var R=Matrix4x4(); - R.setView(zero,direction,up); - var c=center[k]; - var T=node.transform; - T.setIdentity(); - T.translateInPlace(c.scale(-1)); - T.multiplyInPlace(R); - T.translateInPlace(c); - } - } -"; - for(int i=0; i < index.length; ++i) - s += "f("+string(i)+","+string(index[i])+"); -"; - s += " - runtime.refresh(); -} - -runtime.addEventHandler(billboardHandler); - -runtime.refresh(); -"; - return s; -} - -string lightscript(light light) { - string script="for(var i=scene.lights.count-1; i >= 0; i--) - scene.lights.removeByIndex(i);"+'\n\n'; - for(int i=0; i < light.position.length; ++i) { - string Li="L"+string(i); - real[] diffuse=light.diffuse[i]; - script += Li+"=scene.createLight();"+'\n'+ - Li+".direction.set("+format(-light.position[i],",")+");"+'\n'+ - Li+".color.set("+format((diffuse[0],diffuse[1],diffuse[2]),",")+");"+'\n'; - } - // Work around initialization bug in Adobe Reader 8.0: - return script +" -scene.lightScheme=scene.LIGHT_MODE_HEADLAMP; -scene.lightScheme=scene.LIGHT_MODE_FILE; -"; -} - -void writeJavaScript(string name, string preamble, string script) -{ - file out=output(name); - write(out,preamble); - if(script != "") { - write(out,endl); - file in=input(script); - while(true) { - string line=in; - if(eof(in)) break; - write(out,line,endl); - } - } - close(out); - if(settings.verbose > 1) write("Wrote "+name); - if(!settings.inlinetex) - file3.push(name); -} - -pair viewportmargin(pair lambda) -{ - return maxbound(0.5*(viewportsize-lambda),viewportmargin); -} - -string embed3D(string label="", string text=label, string prefix, - frame f, string format="", - real width=0, real height=0, - string options="", string script="", - light light=currentlight, projection P=currentprojection) -{ - if(!prc(format) || Embed == null) return ""; - - if(width == 0) width=settings.paperwidth; - if(height == 0) height=settings.paperheight; - - if(script == "") script=defaultembed3Dscript; - - // Adobe Reader doesn't appear to support user-specified viewport lights. - string lightscript=light.on() && light.viewport ? "" : lightscript(light); - - real viewplanesize; - if(P.infinity) { - triple lambda=max3(f)-min3(f); - pair margin=viewportmargin((lambda.x,lambda.y)); - viewplanesize=(max(lambda.x+2*margin.x,lambda.y+2*margin.y))/P.zoom; - } else - if(!P.absolute) P.angle=2*aTan(Tan(0.5*P.angle)); - - int[] index; - triple[] center; - shipout3(prefix,f,index,center); - - string name=prefix+".js"; - writeJavaScript(name,lightscript+projection(P.infinity,viewplanesize)+ - billboard(index,center),script); - - if(!settings.inlinetex) - file3.push(prefix+".prc"); - - triple target=P.target; - if(P.viewportshift != 0) { - triple lambda=max3(f)-min3(f); - target -= (P.viewportshift.x*lambda.x/P.zoom, - P.viewportshift.y*lambda.y/P.zoom,0); - } - - triple v=P.vector(); - triple u=unit(v); - triple w=Z-u.z*u; - real roll; - if(abs(w) > sqrtEpsilon) { - w=unit(w); - triple up=unit(perp(P.up,u)); - roll=degrees(acos1(dot(up,w)))*sgn(dot(cross(up,w),u)); - } else roll=0; - - string options3=light.viewport ? "3Dlights=Headlamp" : "3Dlights=File"; - if(defaultembed3Doptions != "") options3 += ","+defaultembed3Doptions; - if((settings.render < 0 || !settings.embed) && settings.auto3D) - options3 += ",poster"; - options3 += ",text={"+text+"},label="+label+ - ",toolbar="+(settings.toolbar ? "true" : "false")+ - ",3Daac="+Format(P.angle)+ - ",3Dc2c="+Format(u)+ - ",3Dcoo="+Format(target)+ - ",3Droll="+Format(roll)+ - ",3Droo="+Format(abs(v))+ - ",3Dbg="+Format(light.background()); - if(options != "") options3 += ","+options; - if(settings.inlinetex) - prefix=jobname(prefix); - options3 += ",3Djscript="+prefix+".js"; - - return Embed(prefix+".prc",options3,width,height); -} - -struct scene -{ - frame f; - transform3 t; - projection P; - bool adjusted; - real width,height; - pair viewportmargin; - transform3 T=identity4; - picture pic2; - - void operator init(frame f, real width, real height, - projection P=currentprojection) { - this.f=f; - this.t=identity4; - this.P=P; - this.width=width; - this.height=height; - } - - void operator init(picture pic, real xsize=pic.xsize, real ysize=pic.ysize, - bool keepAspect=pic.keepAspect, bool is3D=true, - projection P=currentprojection) { - real xsize3=pic.xsize3, ysize3=pic.ysize3, zsize3=pic.zsize3; - bool warn=true; - - if(xsize3 == 0 && ysize3 == 0 && zsize3 == 0) { - xsize3=ysize3=zsize3=max(xsize,ysize); - warn=false; - } - - if(P.absolute) - this.P=P.copy(); - else if(P.showtarget) - draw(pic,P.target,nullpen); - - t=pic.scaling(xsize3,ysize3,zsize3,keepAspect,warn); - adjusted=false; - triple m=pic.min(t); - triple M=pic.max(t); - - if(!P.absolute) { - this.P=t*P; - if(this.P.center && settings.render != 0) { - triple target=0.5*(m+M); - this.P.target=target; - this.P.calculate(); - } - if(this.P.autoadjust || this.P.infinity) - adjusted=adjusted | this.P.adjust(m,M); - } - - bool scale=xsize != 0 || ysize != 0; - bool scaleAdjust=scale && this.P.autoadjust; - bool noAdjust=(this.P.absolute || !scaleAdjust); - - if(pic.bounds3.exact && noAdjust) - this.P.bboxonly=false; - - f=pic.fit3(t,pic.bounds3.exact ? pic2 : null,this.P); - - if(!pic.bounds3.exact) { - if(noAdjust) - this.P.bboxonly=false; - - transform3 s=pic.scale3(f,xsize3,ysize3,zsize3,keepAspect); - t=s*t; - this.P=s*this.P; - f=pic.fit3(t,pic2,this.P); - } - - if(is3D || scale) { - pic2.bounds.exact=true; - transform s=pic2.scaling(xsize,ysize,keepAspect); - - pair m2=pic2.min(s); - pair M2=pic2.max(s); - pair lambda=M2-m2; - viewportmargin=viewportmargin(lambda); - width=ceil(lambda.x+2*viewportmargin.x); - height=ceil(lambda.y+2*viewportmargin.y); - - if(!this.P.absolute) { - if(scaleAdjust) { - pair v=(s.xx,s.yy); - transform3 T=this.P.t; - pair x=project(X,T); - pair y=project(Y,T); - pair z=project(Z,T); - real f(pair a, pair b) { - return b == 0 ? (0.5*(a.x+a.y)) : - (b.x^2*a.x+b.y^2*a.y)/(b.x^2+b.y^2); - } - transform3 s=keepAspect ? scale3(min(f(v,x),f(v,y),f(v,z))) : - xscale3(f(v,x))*yscale3(f(v,y))*zscale3(f(v,z)); - s=shift(this.P.target)*s*shift(-this.P.target); - t=s*t; - this.P=s*this.P; - this.P.bboxonly=false; - if(!is3D) pic2.erase(); - f=pic.fit3(t,is3D ? null : pic2,this.P); - } - - if(this.P.autoadjust || this.P.infinity) - adjusted=adjusted | this.P.adjust(min3(f),max3(f)); - } - } - } - - // Choose the angle to be just large enough to view the entire image. - real angle(projection P) { - T=identity4; - real h=-0.5*P.target.z; - pair r,R; - real diff=realMax; - pair s; - int i; - do { - r=minratio(f); - R=maxratio(f); - pair lasts=s; - if(P.autoadjust) { - s=r+R; - if(s != 0) { - transform3 t=shift(h*s.x,h*s.y,0); - f=t*f; - T=t*T; - adjusted=true; - } - } - diff=abs(s-lasts); - ++i; - } while (diff > angleprecision && i < maxangleiterations); - real aspect=width > 0 ? height/width : 1; - real rx=-r.x*aspect; - real Rx=R.x*aspect; - real ry=-r.y; - real Ry=R.y; - if(!P.autoadjust) { - if(rx > Rx) Rx=rx; - else rx=Rx; - if(ry > Ry) Ry=ry; - else ry=Ry; - } - return (1+angleprecision)*max(aTan(rx)+aTan(Rx),aTan(ry)+aTan(Ry)); - } -} - -object embed(string label="", string text=label, string prefix=defaultfilename, - scene S, string format="", bool view=true, string options="", - string script="", light light=currentlight) -{ - object F; - transform3 modelview; - projection P=S.P; - transform3 tinv=inverse(S.t); - - projection Q; - modelview=P.T.modelview; - if(P.absolute) { - Q=modelview*P; - } else { - triple target=P.target; - S.f=modelview*S.f; - P=modelview*P; - Q=P.copy(); - light=modelview*light; - - if(P.infinity) { - triple m=min3(S.f); - triple M=max3(S.f); - triple lambda=M-m; - S.viewportmargin=viewportmargin((lambda.x,lambda.y)); - S.width=ceil(lambda.x+2*S.viewportmargin.x); - S.height=ceil(lambda.y+2*S.viewportmargin.y); - S.f=shift((-0.5(m.x+M.x),-0.5*(m.y+M.y),0))*S.f; // Eye will be at (0,0,0) - } else { - if(P.angle == 0) { - P.angle=S.angle(P); - modelview=S.T*modelview; - if(S.viewportmargin.y != 0) - P.angle=2*aTan(Tan(0.5*P.angle)-S.viewportmargin.y/P.target.z); - } - Q.angle=P.angle; - if(settings.verbose > 0) { - transform3 inv=inverse(modelview); - if(S.adjusted) - write("adjusting camera to ",tinv*inv*P.camera); - target=inv*P.target; - } - P=S.T*P; - } - if(settings.verbose > 0) { - if((P.center && settings.render != 0) || (!P.infinity && P.autoadjust)) - write("adjusting target to ",tinv*target); - } - } - - if(prefix == "") prefix=outprefix(); - bool prc=prc(format); - bool preview=settings.render > 0; - if(prc) { - // The movie15.sty package cannot handle spaces or dots in filenames. - string dir=stripfile(prefix); - prefix=dir+replace(stripdirectory(prefix), - new string[][]{{" ","_"},{".","_"}}); - if(settings.embed || nativeformat() == "pdf") - prefix += "+"+(string) file3.length; - } else - preview=false; - if(preview || (!prc && settings.render != 0)) { - frame f=S.f; - triple m,M; - real zcenter; - real r; - if(P.absolute) { - f=modelview*f; - m=min3(f); - M=max3(f); - r=0.5*abs(M-m); - zcenter=0.5*(M.z+m.z); - } else { - m=min3(f); - M=max3(f); - zcenter=P.target.z; - r=P.distance(m,M); - } - M=(M.x,M.y,zcenter+r); - m=(m.x,m.y,zcenter-r); - - if(P.infinity) { - triple margin=(S.viewportmargin.x,S.viewportmargin.y,0); - M += margin; - m -= margin; - } else if(M.z >= 0) abort("camera too close"); - - shipout3(prefix,f,preview ? nativeformat() : format, - S.width-defaultrender.margin,S.height-defaultrender.margin, - P.infinity ? 0 : 2aTan(Tan(0.5*P.angle)*P.zoom), - P.zoom,m,M,P.viewportshift, - tinv*inverse(modelview)*shift(0,0,zcenter),light.background(), - P.absolute ? (modelview*light).position : light.position, - light.diffuse,light.ambient,light.specular, - light.viewport,view && !preview); - if(!preview) return F; - } - - string image; - if((preview || (prc && settings.render == 0)) && settings.embed) { - image=prefix; - if(settings.inlinetex) image += "_0"; - if(!preview && !shipped && !S.pic2.empty2()) { - transform T=S.pic2.scaling(S.width,S.height); - shipout(image,S.pic2.fit(T),newframe,nativeformat(),false,false,null); - } - - image += "."+nativeformat(); - if(!settings.inlinetex) file3.push(image); - image=graphic(image,"hiresbb"); - } - if(prc) { - if(!P.infinity && P.viewportshift != 0) - warning("offaxis", - "PRC does not support off-axis projections; use pan instead of -shift"); - F.L=embed3D(label,text=image,prefix,S.f,format, - S.width-2,S.height-2,options,script,light,Q); - } - return F; -} - -object embed(string label="", string text=label, - string prefix=defaultfilename, - picture pic, string format="", - real xsize=pic.xsize, real ysize=pic.ysize, - bool keepAspect=pic.keepAspect, bool view=true, string options="", - string script="", light light=currentlight, - projection P=currentprojection) -{ - bool is3D=is3D(format); - scene S=scene(pic,xsize,ysize,keepAspect,is3D,P); - if(is3D) - return embed(label,text,prefix,S,format,view,options,script,light); - else { - object F; - transform T=S.pic2.scaling(xsize,ysize,keepAspect); - F.f=pic.fit(scale(S.t[0][0])*T); - add(F.f,S.pic2.fit()); - return F; - } -} - -object embed(string label="", string text=label, - string prefix=defaultfilename, - frame f, string format="", real width=0, real height=0, - bool view=true, string options="", string script="", - light light=currentlight, projection P=currentprojection) -{ - if(is3D(format)) - return embed(label,text,prefix,scene(f,width,height,P),format,view,options, - script,light); - else { - object F; - F.f=f; - return F; - } -} - -embed3=new object(string prefix, frame f, string format, string options, - string script, light light, projection P) { - return embed(prefix=prefix,f,format,options,script,light,P); -}; - -frame embedder(object embedder(string prefix, string format), - string prefix, string format, bool view, light light) -{ - frame f; - bool prc=prc(format); - if(!prc && settings.render != 0 && !view) { - static int previewcount=0; - bool keep=prefix != ""; - prefix=outprefix(prefix)+"+"+(string) previewcount; - ++previewcount; - format=nativeformat(); - if(!keep) file3.push(prefix+"."+format); - } - object F=embedder(prefix,format); - if(prc) - label(f,F.L); - else { - if(settings.render == 0) { - add(f,F.f); - if(light.background != nullpen) - box(f,light.background,Fill,above=false); - } else if(!view) - label(f,graphic(prefix,"hiresbb")); - } - return f; -} - -currentpicture.fitter=new frame(string prefix, picture pic, string format, - real xsize, real ysize, bool keepAspect, - bool view, string options, string script, - light light, projection P) { - frame f; - bool empty3=pic.empty3(); - if(!empty3) f=embedder(new object(string prefix, string format) { - return embed(prefix=prefix,pic,format,xsize,ysize,keepAspect,view, - options,script,light,P); - },prefix,format,view,light); - if(is3D(format) || empty3) add(f,pic.fit2(xsize,ysize,keepAspect)); - return f; -}; - -frame embedder(string prefix, frame f, string format, real width, real height, - bool view, string options, string script, light light, - projection P) -{ - return embedder(new object(string prefix, string format) { - return embed(prefix=prefix,f,format,width,height,view,options,script, - light,P); - },prefix,format,view,light); -} - -projection[][] ThreeViewsUS={{TopView}, - {FrontView,RightView}}; - -projection[][] SixViewsUS={{null,TopView}, - {LeftView,FrontView,RightView,BackView}, - {null,BottomView}}; - -projection[][] ThreeViewsFR={{RightView,FrontView}, - {null,TopView}}; - -projection[][] SixViewsFR={{null,BottomView}, - {RightView,FrontView,LeftView,BackView}, - {null,TopView}}; - -projection[][] ThreeViews={{FrontView,TopView,RightView}}; - -projection[][] SixViews={{FrontView,TopView,RightView}, - {BackView,BottomView,LeftView}}; - -void addViews(picture dest, picture src, projection[][] views=SixViewsUS, - bool group=true, filltype filltype=NoFill) -{ - frame[][] F=array(views.length,new frame[]); - pair[][] M=array(views.length,new pair[]); - pair[][] m=array(views.length,new pair[]); - - for(int i=0; i < views.length; ++i) { - projection[] viewsi=views[i]; - frame[] Fi=F[i]; - pair[] Mi=M[i]; - pair[] mi=m[i]; - for(projection P : viewsi) { - if(P != null) { - frame f=src.fit(P); - mi.push(min(f)); - Mi.push(max(f)); - Fi.push(f); - } else { - pair Infinity=(infinity,infinity); - mi.push(Infinity); - Mi.push(-Infinity); - Fi.push(newframe); - } - } - } - - real[] my=new real[views.length]; - real[] My=new real[views.length]; - - int Nj=0; - for(int i=0; i < views.length; ++i) { - my[i]=minbound(m[i]).y; - My[i]=maxbound(M[i]).y; - Nj=max(Nj,views[i].length); - } - - real[] mx=array(Nj,infinity); - real[] Mx=array(Nj,-infinity); - for(int i=0; i < views.length; ++i) { - pair[] mi=m[i]; - pair[] Mi=M[i]; - for(int j=0; j < views[i].length; ++j) { - mx[j]=min(mx[j],mi[j].x); - Mx[j]=max(Mx[j],Mi[j].x); - } - } - - if(group) begingroup(dest); - - real y; - for(int i=0; i < views.length; ++i) { - real x; - pair[] mi=m[i]; - for(int j=0; j < views[i].length; ++j) { - if(size(F[i][j]) != 0) - add(dest,shift(x-mx[j],y+my[i])*F[i][j],filltype); - x += (Mx[j]-mx[j]); - } - y -= (My[i]-my[i]); - } - - if(group) endgroup(dest); -} - -void addViews(picture src, projection[][] views=SixViewsUS, bool group=true, - filltype filltype=NoFill) -{ - addViews(currentpicture,src,views,group,filltype); -} - -void addStereoViews(picture dest, picture src, bool group=true, - filltype filltype=NoFill, real eyetoview=defaulteyetoview, - bool leftright=true, projection P=currentprojection) -{ - triple v=P.vector(); - triple h=0.5*abs(v)*eyetoview*unit(cross(P.up,v)); - projection leftEye=P.copy(); - leftEye.camera -= h; - leftEye.calculate(); - projection rightEye=P.copy(); - rightEye.camera += h; - rightEye.calculate(); - addViews(dest,src,leftright ? - new projection[][] {{leftEye,rightEye}} : - new projection[][] {{rightEye,leftEye}},group,filltype); -} - -void addStereoViews(picture src, bool group=true, - filltype filltype=NoFill, - real eyetoview=defaulteyetoview, bool leftright=true, - projection P=currentprojection) -{ - addStereoViews(currentpicture,src,group,filltype,eyetoview,leftright,P); -} - -// Fit an array of 3D pictures simultaneously using the sizing of picture all. -frame[] fit3(string prefix="", picture[] pictures, picture all, - string format="", bool view=true, string options="", - string script="", light light=currentlight, - projection P=currentprojection) -{ - frame[] out; - scene S=scene(all,P); - triple m=all.min(S.t); - triple M=all.max(S.t); - out=new frame[pictures.length]; - int i=0; - bool reverse=settings.reverse; - settings.animating=true; - - for(picture pic : pictures) { - picture pic2; - frame f=pic.fit3(S.t,pic2,S.P); - if(settings.interrupt) break; - add(f,pic2.fit2()); - draw(f,m,nullpen); - draw(f,M,nullpen); - out[i]=f; - ++i; - } - - while(!settings.interrupt) { - for(int i=settings.reverse ? pictures.length-1 : 0; - i >= 0 && i < pictures.length && !settings.interrupt; - settings.reverse ? --i : ++i) { - frame f=embedder(prefix,out[i],format,S.width,S.height,view,options, - script,light,S.P); - if(!settings.loop) out[i]=f; - } - if(!settings.loop) break; - } - - settings.animating=false; - settings.interrupt=false; - settings.reverse=reverse; - - return out; -} - -// Fit an array of pictures simultaneously using the size of the first picture. -fit=new frame[](string prefix="", picture[] pictures, string format="", - bool view=true, string options="", string script="", - projection P=currentprojection) { - if(pictures.length == 0) - return new frame[]; - - picture all; - size(all,pictures[0]); - for(picture pic : pictures) - add(all,pic); - - return all.empty3() ? fit2(pictures,all) : - fit3(prefix,pictures,all,format,view,options,script,P); -}; - -// Add frame src to picture dest about position. -void add(picture dest=currentpicture, frame src, triple position) -{ - if(is3D(src)) { - dest.add(new void(frame f, transform3 t, picture, projection) { - add(f,shift(t*position)*src); - },true); - } else { - dest.add(new void(frame, transform3 t, picture pic, projection P) { - if(pic != null) { - pic.add(new void(frame f, transform T) { - add(f,T*shift(project(t*position,P))*src); - },true); - } - },true); - } - dest.addBox(position,position,min3(src),max3(src)); -} - -exitfcn currentexitfunction=atexit(); - -void exitfunction() -{ - if(currentexitfunction != null) currentexitfunction(); - if(!settings.keep) - for(int i=0; i < file3.length; ++i) - delete(file3[i]); - file3=new string[]; -} - -atexit(exitfunction); |