diff options
Diffstat (limited to 'Master/texmf/asymptote')
32 files changed, 1883 insertions, 1182 deletions
diff --git a/Master/texmf/asymptote/GUI/xasyVersion.py b/Master/texmf/asymptote/GUI/xasyVersion.py index a1ab6b203c5..5d5a51724f2 100755 --- a/Master/texmf/asymptote/GUI/xasyVersion.py +++ b/Master/texmf/asymptote/GUI/xasyVersion.py @@ -1,2 +1,2 @@ #!/usr/bin/env python -xasyVersion = "1.82" +xasyVersion = "1.83" diff --git a/Master/texmf/asymptote/animation.asy b/Master/texmf/asymptote/animation.asy index 6fcd2ad70e2..27f7171a639 100644 --- a/Master/texmf/asymptote/animation.asy +++ b/Master/texmf/asymptote/animation.asy @@ -8,16 +8,17 @@ // animation delay is in milliseconds real animationdelay=50; -typedef frame fit(string prefix="",picture); +typedef frame enclosure(frame); -frame NoBox(string prefix="", picture pic) { - return pic.fit(prefix); +frame NoBox(frame f) { + return f; } -fit BBox(string prefix="", real xmargin=0, real ymargin=xmargin, - pen p=currentpen, filltype filltype=NoFill) { - return new frame(string prefix, picture pic) { - return bbox(prefix,pic,xmargin,ymargin,p,filltype); +enclosure BBox(real xmargin=0, real ymargin=xmargin, + pen p=currentpen, filltype filltype=NoFill) { + return new frame(frame f) { + box(f,xmargin,ymargin,p,filltype,above=false); + return f; }; } @@ -59,11 +60,11 @@ struct animation { shipped=false; } - void add(picture pic=currentpicture, fit fit=NoBox) { + void add(picture pic=currentpicture, enclosure enclosure=NoBox) { if(global) { ++index; pictures.push(pic.copy()); - } else this.shipout(fit(prefix,pic)); + } else this.shipout(enclosure(pic.fit())); } void purge(bool keep=settings.keep) { @@ -86,30 +87,33 @@ struct animation { return rc; } + void glmovie(string prefix=prefix, projection P=currentprojection) { + if(!view() || settings.render == 0) return; + fit(prefix,pictures,view=true,P); + } + // Export all frames with the same scaling. - void export(string prefix=prefix, fit fit=NoBox, - bool multipage=false, bool view=false) { + void export(string prefix=prefix, enclosure enclosure=NoBox, + bool multipage=false, bool view=false, + projection P=currentprojection) { if(pictures.length == 0) return; if(!global) multipage=false; - rescale(pictures); - frame multi; bool inlinetex=settings.inlinetex; if(multipage) settings.inlinetex=false; - for(int i=0; i < pictures.length; ++i) { + frame multi; + frame[] fits=fit(prefix,pictures,view=false,P); + for(int i=0; i < fits.length; ++i) { + string s=name(prefix,i); if(multipage) { - add(multi,fit(prefix,pictures[i])); + add(multi,enclosure(fits[i])); newpage(multi); + files.push(s+"."+nativeformat()); } else { - if(pictures[i].empty3() || settings.render <= 0) { - real render=settings.render; - settings.render=0; - this.shipout(name(prefix,i),fit(prefix,pictures[i])); - settings.render=render; - } else { // Render 3D frames - files.push(name(prefix,i)+"."+nativeformat()); - fit(prefix,pictures[i]); - } + if(pictures[i].empty3() || settings.render <= 0) + this.shipout(s,enclosure(fits[i])); + else // 3D frames + files.push(s+"."+nativeformat()); } } if(multipage) { @@ -129,8 +133,8 @@ struct animation { return s; } - string pdf(fit fit=NoBox, real delay=animationdelay, string options="", - bool keep=settings.keep, bool multipage=true) { + string pdf(enclosure enclosure=NoBox, real delay=animationdelay, + string options="", bool keep=settings.keep, bool multipage=true) { if(settings.inlinetex) multipage=true; if(!global) multipage=false; if(settings.tex != "pdflatex") @@ -141,16 +145,17 @@ struct animation { string pdfname=filename+".pdf"; if(global) - export(filename,fit,multipage=multipage); + export(filename,enclosure,multipage=multipage); shipped=false; - if(!keep && !settings.inlinetex) { + if(!keep) { exitfcn currentexitfunction=atexit(); void exitfunction() { if(currentexitfunction != null) currentexitfunction(); - this.purge(); - if(multipage) + if(multipage || !settings.inlinetex) + this.purge(); + if(multipage && !settings.inlinetex) delete(pdfname); } atexit(exitfunction); @@ -162,15 +167,15 @@ struct animation { return load(index,delay,options,multipage); } - int movie(fit fit=NoBox, int loops=0, real delay=animationdelay, + int movie(enclosure enclosure=NoBox, int loops=0, real delay=animationdelay, string format=settings.outformat == "" ? "gif" : settings.outformat, string options="", bool keep=settings.keep) { if(global) { if(format == "pdf") { - export(fit,multipage=true,view=true); + export(enclosure,multipage=true,view=true); return 0; } - export(fit); + export(enclosure); } return merge(loops,delay,format,options,keep); } diff --git a/Master/texmf/asymptote/bezulate.asy b/Master/texmf/asymptote/bezulate.asy index 2b852623faf..d935d358989 100644 --- a/Master/texmf/asymptote/bezulate.asy +++ b/Master/texmf/asymptote/bezulate.asy @@ -238,7 +238,7 @@ path[] bezulate(path[] p) // avoid infinite recursion ++refinements; if(refinements > mantissaBits) { - write("warning: too many subdivisions"); + warning("subdivisions","too many subdivisions",position=true); } else { p=subdivide(p); i=-1; diff --git a/Master/texmf/asymptote/geometry.asy b/Master/texmf/asymptote/geometry.asy index 12e35339dcc..0db3172d365 100644 --- a/Master/texmf/asymptote/geometry.asy +++ b/Master/texmf/asymptote/geometry.asy @@ -228,24 +228,6 @@ pair[] intersectionpoints(pair A, pair B, real[] equation) // *=======================================================* // *......................COORDINATES......................* -// Copyright (c) 2007, Philippe Ivaldi. -// Version: $Id: coordinates.asy,v 0.0 2007/02/03 16:06:23 Philippe Ivaldi Exp$ -// Last modified: Wed Aug 15 15:53:01 CEST 2007 - -// This program is free software; you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation; either version 3 of the License, or -// any later version. - -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. - -// You should have received a copy of the GNU General Public License -// along with this program; if not, write to the Free Software -// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA -// 02110-1301, USA. real EPS=sqrt(realEpsilon); @@ -492,8 +474,9 @@ bool samecoordsys(bool warn=true ... point[] M) t=M[i].coordsys; } if(warn && !ret) - write("Warning, the coordinate system of two objects are not the same. -The operation will be done relatively to the default coordinate system."); + warning("coodinatesystem", + "the coordinate system of two objects are not the same. +The operation will be done relative to the default coordinate system."); return ret; } @@ -2418,8 +2401,10 @@ bool samecoordsys(bool warn=true ... bqe[] bqes) t=bqes[i].coordsys; } if(warn && !ret) - write("Warning, the coordinate system of two bivariate quadratic equations are not the same. -The operation will be done relatively to the default coordinate system."); + warning("coodinatesystem", + "the coordinate system of two bivariate quadratic equations are not +the same. The operation will be done relatively to the default coordinate +system."); return ret; } @@ -2710,7 +2695,9 @@ int circlenodesnumberfactor=100;/*<asyxml></code><documentation>Factor for the n /*<asyxml><function type="int" signature="circlenodesnumber(real)"><code></asyxml>*/ int circlenodesnumber(real r) {/*<asyxml></code><documentation>Return the number of nodes for drawing a circle of radius 'r'.</documentation></function></asyxml>*/ - if (circlenodesnumberfactor < 100) write("Warning: variable 'circlenodesnumberfactor' maybe too small."); + if (circlenodesnumberfactor < 100) + warning("circlenodesnumberfactor", + "variable 'circlenodesnumberfactor' may be too small."); int oi=ceil(circlenodesnumberfactor*abs(r)^0.1); oi=45*floor(oi/45); return oi == 0 ? 4 : conicnodesfactor*oi; @@ -2729,7 +2716,9 @@ int ellipsenodesnumberfactor=250;/*<asyxml></code><documentation>Factor for the /*<asyxml><function type="int" signature="ellipsenodesnumber(real,real)"><code></asyxml>*/ int ellipsenodesnumber(real a, real b) {/*<asyxml></code><documentation>Return the number of nodes to draw a ellipse of axis 'a' and 'b'.</documentation></function></asyxml>*/ - if (ellipsenodesnumberfactor < 250) write("Warning: variable 'ellipsenodesnumberfactor' maybe too small."); + if (ellipsenodesnumberfactor < 250) + write("ellipsenodesnumberfactor", + "variable 'ellipsenodesnumberfactor' maybe too small."); int tmp=circlenodesnumberfactor; circlenodesnumberfactor=ellipsenodesnumberfactor; int oi=circlenodesnumber(max(abs(a),abs(b))/min(abs(a),abs(b))); @@ -5481,7 +5470,7 @@ circle excircle(point A, point B, point C) } private int[] numarray={1,2,3}; -numarray.cyclic(true); +numarray.cyclic=true; /*<asyxml><struct signature="triangle"><code></asyxml>*/ struct triangle {/*<asyxml></code><documentation></documentation></asyxml>*/ @@ -6440,13 +6429,19 @@ point operator *(inversion i, point P) return inverse(i.k,i.C,P); } +void lineinversion() +{ + warning("lineinversion","the inversion of the line is not a circle. +The returned circle has an infinite radius, circle.l has been set."); +} + + /*<asyxml><function type="circle" signature="inverse(real,point,line)"><code></asyxml>*/ circle inverse(real k, point A, line l) {/*<asyxml></code><documentation>Return the inverse circle of 'l' with respect to point 'A' and inversion radius 'k'.</documentation></function></asyxml>*/ if(A @ l) { - write("Warning: the inversion of the line is not a circle."); - write("The returned circle has an infinite radius, cirlce.l have been set."); + lineinversion(); circle C=circle(A, infinity); C.l=l; return C; @@ -6467,8 +6462,7 @@ circle inverse(real k, point A, circle c) respect to point A and inversion radius 'k'.</documentation></function></asyxml>*/ if(degenerate(c)) return inverse(k,A,c.l); if(A @ c) { - write("Warning: the inversion of the circle is not a circle."); - write("The returned circle has an infinite radius, cirlce.l have been set."); + lineinversion(); point M=rotate(180,c.C)*A, Mp=rotate(90,c.C)*A; circle oc=circle(A,infinity); oc.l=line(inverse(k,A,M),inverse(k,A,Mp)); @@ -7123,7 +7117,7 @@ path square(pair z1, pair z2) // relative to the path z--z+dir. void perpendicular(picture pic=currentpicture, pair z, pair align, pair dir=E, real size=0, pen p=currentpen, - margin margin=NoMargin, filltype filltype=NoFill) + margin margin=NoMargin, filltype filltype=NoFill) { perpendicularmark(pic,(point) z,align,dir,size,p,margin,filltype); } @@ -7133,7 +7127,7 @@ void perpendicular(picture pic=currentpicture, pair z, pair align, // relative to the path z--z+dir(g,0) void perpendicular(picture pic=currentpicture, pair z, pair align, path g, real size=0, pen p=currentpen, margin margin=NoMargin, - filltype filltype=NoFill) + filltype filltype=NoFill) { perpendicularmark(pic,(point) z,align,dir(g,0),size,p,margin,filltype); } diff --git a/Master/texmf/asymptote/graph.asy b/Master/texmf/asymptote/graph.asy index 2983137b612..55974c63f7d 100644 --- a/Master/texmf/asymptote/graph.asy +++ b/Master/texmf/asymptote/graph.asy @@ -1,6 +1,6 @@ private import math; -import splinetype; +import graph_splinetype; import graph_settings; scaleT Linear; @@ -1893,7 +1893,7 @@ interpolate Hermite(splinetype splinetype) }; } -interpolate Hermite=Hermite(defaultspline); +interpolate Hermite=Hermite(Spline); guide graph(picture pic=currentpicture, real f(real), real a, real b, int n=ngraph, real T(real)=identity, interpolate join=operator --) diff --git a/Master/texmf/asymptote/graph3.asy b/Master/texmf/asymptote/graph3.asy index 8fe6c9ed5c6..c2630e10f11 100644 --- a/Master/texmf/asymptote/graph3.asy +++ b/Master/texmf/asymptote/graph3.asy @@ -1509,6 +1509,27 @@ path3[] segment(triple[] v, bool[] cond, interpolate3 join=operator --) segment.length); } +bool uperiodic(triple[][] a) { + int n=a.length; + if(n == 0) return false; + int m=a[0].length; + triple[] a0=a[0]; + triple[] a1=a[n-1]; + real epsilon=sqrtEpsilon*norm(a); + for(int j=0; j < m; ++j) + if(abs(a0[j]-a1[j]) > epsilon) return false; + return true; +} +bool vperiodic(triple[][] a) { + int n=a.length; + if(n == 0) return false; + int m=a[0].length-1; + real epsilon=sqrtEpsilon*norm(a); + for(int i=0; i < n; ++i) + if(abs(a[i][0]-a[i][m]) > epsilon) return false; + return true; +} + // return the surface described by a matrix f surface surface(triple[][] f, bool[][] cond={}) { @@ -1533,6 +1554,7 @@ surface surface(triple[][] f, bool[][] cond={}) } surface s=surface(count); + s.index=new int[nx][ny]; int k=-1; for(int i=0; i < nx; ++i) { bool[] condi,condp; @@ -1542,122 +1564,18 @@ surface surface(triple[][] f, bool[][] cond={}) } triple[] fi=f[i]; triple[] fp=f[i+1]; + int[] indexi=s.index[i]; for(int j=0; j < ny; ++j) { if(all || (condi[j] && condi[j+1] && condp[j] && condp[j+1])) s.s[++k]=patch(new triple[] {fi[j],fp[j],fp[j+1],fi[j+1]}); - } - } - return s; -} - -private surface bispline(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]; - for(int j=0; j < m; ++j) - if(condi[j]) ++count; - } - } - - surface g=surface(count); - int k=-1; - for(int i=0; i < n; ++i) { - bool[] condi=all ? null : cond[i]; - real xi=x[i]; - real[] zi=z[i]; - real[] zp=z[i+1]; - real[] ri=r[i]; - real[] rp=r[i+1]; - real[] pi=p[i]; - real[] pp=p[i+1]; - real[] qi=q[i]; - real[] qp=q[i+1]; - real xp=x[i+1]; - real hx=(xp-xi)/3; - for(int j=0; j < m; ++j) { - real yj=y[j]; - real yp=y[j+1]; - if(all || condi[j]) { - triple[][] P=array(4,array(4,O)); - real hy=(yp-yj)/3; - real hxy=hx*hy; - // first x and y directions - for(int k=0; k < 4; ++k) { - P[k][0] += xi*X; - P[0][k] += yj*Y; - P[k][1] += (xp+2*xi)/3*X; - P[1][k] += (yp+2*yj)/3*Y; - P[k][2] += (2*xp+xi)/3*X; - P[2][k] += (2*yp+yj)/3*Y; - P[k][3] += xp*X; - P[3][k] += yp*Y; - } - // now z, first the value - P[0][0] += zi[j]*Z; - P[0][3] += zp[j]*Z; - P[3][0] += zi[j+1]*Z; - P[3][3] += zp[j+1]*Z; - // 2nd, first derivative - P[0][1] += (P[0][0].z+hx*pi[j])*Z; - P[3][1] += (P[3][0].z+hx*pi[j+1])*Z; - P[0][2] += (P[0][3].z-hx*pp[j])*Z; - P[3][2] += (P[3][3].z-hx*pp[j+1])*Z; - P[1][0] += (P[0][0].z+hy*qi[j])*Z; - P[1][3] += (P[0][3].z+hy*qp[j])*Z; - P[2][0] += (P[3][0].z-hy*qi[j+1])*Z; - P[2][3] += (P[3][3].z-hy*qp[j+1])*Z; - // 3nd, second derivative - P[1][1] += (P[1][0].z+P[0][1].z-P[0][0].z+hxy*ri[j])*Z; - P[2][1] += (P[2][0].z+P[3][1].z-P[3][0].z-hxy*ri[j+1])*Z; - P[1][2] += (P[0][2].z+P[1][3].z-P[0][3].z-hxy*rp[j])*Z; - P[2][2] += (P[3][2].z+P[2][3].z-P[3][3].z+hxy*rp[j+1])*Z; - g.s[++k]=patch(P); + indexi[j]=k; } - } } - return g; -} -// return the surface described by a real matrix f, interpolated with -// xsplinetype and ysplinetype. -surface surface(real[][] f, real[] x, real[] y, - splinetype xsplinetype=null, splinetype ysplinetype=xsplinetype, - bool[][] cond={}) -{ - real epsilon=sqrtEpsilon*max(abs(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 bispline(f,p,q,r,x,y,cond); + if(uperiodic(f)) s.ucyclic(true); + if(vperiodic(f)) s.vcyclic(true); + + return s; } // return the surface described by a real matrix f, interpolated with @@ -1738,74 +1656,8 @@ surface surface(triple f(pair z), pair a, pair b, int nu=nmesh, int nv=nu, splinetype[] usplinetype, splinetype[] vsplinetype=Spline, bool cond(pair z)=null) { - real[] upt=uniform(a.x,b.x,nu); - real[] vpt=uniform(a.y,b.y,nv); - real[] ipt=sequence(nu+1); - real[] jpt=sequence(nv+1); - real[][] fx=new real[nu+1][nv+1]; - real[][] fy=new real[nu+1][nv+1]; - real[][] fz=new real[nu+1][nv+1]; - - bool[][] active; - bool all=cond == null; - if(!all) active=new bool[nu+1][nv+1]; - - real norm; - for(int i=0; i <= nu; ++i) { - real upti=upt[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=(upti,vpt[j]); - triple f=(all || (activei[j]=cond(z))) ? f(z) : O; - norm=max(norm,abs(f)); - fxi[j]=f.x; - fyi[j]=f.y; - fzi[j]=f.z; - } - } - - real epsilon=sqrtEpsilon*norm; - - if(usplinetype.length == 0) { - bool uperiodic(real[][] a) { - for(int j=0; j < nv; ++j) - if(abs(a[0][j]-a[nu][j]) > epsilon) return false; - return true; - } - 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) { - bool vperiodic(real[][] a) { - for(int i=0; i < nu; ++i) - if(abs(a[i][0]-a[i][nv]) > epsilon) return false; - return true; - } - 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"); - - surface sx=surface(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active); - surface sy=surface(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active); - surface sz=surface(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active); - - surface s=surface(sx.s.length); - for(int k=0; k < sx.s.length; ++k) { - triple[][] Q=new triple[4][4]; - for(int i=0; i < 4 ; ++i) { - for(int j=0; j < 4 ; ++j) { - Q[i][j]=(sx.s[k].P[i][j].z,sy.s[k].P[i][j].z,sz.s[k].P[i][j].z); - } - s.s[k]=patch(Q); - } - } - return s; + return surface(f,uniform(a.x,b.x,nu),uniform(a.y,b.y,nv), + usplinetype,vsplinetype,cond); } // return the surface described by a real function f over box(a,b), diff --git a/Master/texmf/asymptote/graph_splinetype.asy b/Master/texmf/asymptote/graph_splinetype.asy new file mode 100644 index 00000000000..fbcb1b00c75 --- /dev/null +++ b/Master/texmf/asymptote/graph_splinetype.asy @@ -0,0 +1,251 @@ +typedef real[] splinetype(real[], real[]); + +restricted real[] Spline(real[] x, real[] y); +restricted splinetype[] Spline; + +string morepoints="interpolation requires at least 2 points"; +string differentlengths="arrays have different lengths"; +void checklengths(int x, int y, string text=differentlengths) +{ + if(x != y) + abort(text+": "+string(x)+" != "+string(y)); +} + +// Linear interpolation +real[] linear(real[] x, real[] y) +{ + int n=x.length; + checklengths(n,y.length); + real[] d=new real[n]; + for(int i=0; i < n-1; ++i) + d[i]=(y[i+1]-y[i])/(x[i+1]-x[i]); + d[n-1]=d[n-2]; + return d; +} + +// Standard cubic spline interpolation with not-a-knot condition: +// s'''(x_2^-)=s'''(x_2^+) et s'''(x_(n_2)^-)=s'''(x_(n-2)^+) +// if n=2, linear interpolation is returned +// if n=3, an interpolation polynomial of degree <= 2 is returned: +// p(x_1)=y_1, p(x_2)=y_2, p(x_3)=y_3 +real[] notaknot(real[] x, real[] y) +{ + int n=x.length; + checklengths(n,y.length); + real[] d; + if(n > 3) { + real[] a=new real[n]; + real[] b=new real[n]; + real[] c=new real[n]; + real[] g=new real[n]; + b[0]=x[2]-x[1]; + c[0]=x[2]-x[0]; + a[0]=0; + g[0]=((x[1]-x[0])^2*(y[2]-y[1])/b[0]+b[0]*(2*b[0]+3*(x[1]-x[0]))* + (y[1]-y[0])/(x[1]-x[0]))/c[0]; + for(int i=1; i < n-1; ++i) { + a[i]=x[i+1]-x[i]; + c[i]=x[i]-x[i-1]; + b[i]=2*(a[i]+c[i]); + g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]); + } + c[n-1]=0; + b[n-1]=x[n-2]-x[n-3]; + a[n-1]=x[n-1]-x[n-3]; + g[n-1]=((x[n-1]-x[n-2])^2*(y[n-2]-y[n-3])/b[n-1]+ + b[n-1]*(2*b[n-1]+3(x[n-1]-x[n-2]))* + (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/a[n-1]; + d=tridiagonal(a,b,c,g); + } else if(n == 2) { + real val=(y[1]-y[0])/(x[1]-x[0]); + d=new real[] {val,val}; + } else if(n == 3) { + real a=(y[1]-y[0])/(x[1]-x[0]); + real b=(y[2]-y[1])/(x[2]-x[1]); + real c=(b-a)/(x[2]-x[0]); + d=new real[] {a+c*(x[0]-x[1]),a+c*(x[1]-x[0]),a+c*(2*x[2]-x[0]-x[1])}; + } else abort(morepoints); + return d; +} + +// Standard cubic spline interpolation with periodic condition +// s'(a)=s'(b), s''(a)=s''(b), assuming that f(a)=f(b) +// if n=2, linear interpolation is returned +real[] periodic(real[] x, real[] y) +{ + int n=x.length; + checklengths(n,y.length); + if(abs(y[n-1]-y[0]) > sqrtEpsilon*norm(y)) + abort("function values are not periodic"); + real[] d; + if(n > 2) { + real[] a=new real[n-1]; + real[] b=new real[n-1]; + real[] c=new real[n-1]; + real[] g=new real[n-1]; + c[0]=x[n-1]-x[n-2]; + a[0]=x[1]-x[0]; + b[0]=2*(a[0]+c[0]); + g[0]=3*c[0]*(y[1]-y[0])/a[0]+3*a[0]*(y[n-1]-y[n-2])/c[0]; + for(int i=1; i < n-1; ++i) { + a[i]=x[i+1]-x[i]; + c[i]=x[i]-x[i-1]; + b[i]=2*(a[i]+c[i]); + g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]); + } + d=tridiagonal(a,b,c,g); + d.push(d[0]); + } else if(n == 2) { + d=new real[] {0,0}; + } else abort(morepoints); + return d; +} + +// Standard cubic spline interpolation with the natural condition +// s''(a)=s''(b)=0. +// if n=2, linear interpolation is returned +// Don't use the natural type unless the underlying function +// has zero second end points derivatives. +real[] natural(real[] x, real[] y) +{ + int n=x.length; + checklengths(n,y.length); + real[] d; + if(n > 2) { + real[] a=new real[n]; + real[] b=new real[n]; + real[] c=new real[n]; + real[] g=new real[n]; + b[0]=2*(x[1]-x[0]); + c[0]=x[1]-x[0]; + a[0]=0; + g[0]=3*(y[1]-y[0]); + for(int i=1; i < n-1; ++i) { + a[i]=x[i+1]-x[i]; + c[i]=x[i]-x[i-1]; + b[i]=2*(a[i]+c[i]); + g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]); + } + c[n-1]=0; + a[n-1]=x[n-1]-x[n-2]; + b[n-1]=2*a[n-1]; + g[n-1]=3*(y[n-1]-y[n-2]); + d=tridiagonal(a,b,c,g); + } else if(n == 2) { + real val=(y[1]-y[0])/(x[1]-x[0]); + d=new real[] {val,val}; + } else abort(morepoints); + return d; +} + +// Standard cubic spline interpolation with clamped conditions f'(a), f'(b) +splinetype clamped(real slopea, real slopeb) +{ + return new real[] (real[] x, real[] y) { + int n=x.length; + checklengths(n,y.length); + real[] d; + if(n > 2) { + real[] a=new real[n]; + real[] b=new real[n]; + real[] c=new real[n]; + real[] g=new real[n]; + b[0]=x[1]-x[0]; + g[0]=b[0]*slopea; + c[0]=0; + a[0]=0; + for(int i=1; i < n-1; ++i) { + a[i]=x[i+1]-x[i]; + c[i]=x[i]-x[i-1]; + b[i]=2*(a[i]+c[i]); + g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]); + } + c[n-1]=0; + a[n-1]=0; + b[n-1]=x[n-1]-x[n-2]; + g[n-1]=b[n-1]*slopeb; + d=tridiagonal(a,b,c,g); + } else if(n == 2) { + d=new real[] {slopea,slopeb}; + } else abort(morepoints); + return d; + }; +} + +// Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) +// Modified MATLAB code +// [1] Fritsch, F. N. and R. E. Carlson, +// "Monotone Piecewise Cubic Interpolation," +// SIAM J. Numerical Analysis, Vol. 17, 1980, pp.238-246. +// [2] Kahaner, David, Cleve Moler, Stephen Nash, +// Numerical Methods and Software, Prentice Hall, 1988. +real[] monotonic(real[] x, real[] y) +{ + int n=x.length; + checklengths(n,y.length); + real[] d=new real[n]; + if(n > 2) { + real[] h=new real[n-1]; + real[] del=new real[n-1]; + for(int i=0; i < n-1; ++i) { + h[i]=x[i+1]-x[i]; + del[i]=(y[i+1]-y[i])/h[i]; + } + int j=0; + int k[]=new int[]; + for(int i=0; i < n-2; ++i) + if((sgn(del[i])*sgn(del[i+1])) > 0) {k[j]=i; j=j+1;} + + real[] hs=new real[j]; + for(int i=0; i < j; ++i) hs[i]=h[k[i]]+h[k[i]+1]; + real w1[]=new real[j]; + real w2[]=new real[j]; + real dmax[]=new real[j]; + real dmin[]=new real[j]; + for(int i=0; i < j; ++i) { + w1[i]=(h[k[i]]+hs[i])/(3*hs[i]); + w2[i]=(h[k[i]+1]+hs[i])/(3*hs[i]); + dmax[i]=max(abs(del[k[i]]),abs(del[k[i]+1])); + dmin[i]=min(abs(del[k[i]]),abs(del[k[i]+1])); + } + for(int i=0; i < n; ++i) d[i]=0; + for(int i=0; i < j; ++i) + d[k[i]+1]=dmin[i]/(w1[i]*(del[k[i]]/dmax[i])+w2[i]*(del[k[i]+1]/dmax[i])); + d[0]=((2*h[0]+h[1])*del[0]-h[0]*del[1])/(h[0]+h[1]); + if(sgn(d[0]) != sgn(del[0])) {d[0]=0;} + else if((sgn(del[0]) != sgn(del[1])) && (abs(d[0]) > abs(3*del[0]))) + d[0]=3*del[0]; + + d[n-1]=((2*h[n-2]+h[n-3])*del[n-2]-h[n-2]*del[n-2])/(h[n-2]+h[n-3]); + if(sgn(d[n-1]) != sgn(del[n-2])) {d[n-1]=0;} + else if((sgn(del[n-2]) != sgn(del[n-3])) && + (abs(d[n-1]) > abs(3*del[n-2]))) + d[n-1]=3*del[n-2]; + } else if(n == 2) { + d[0]=d[1]=(y[1]-y[0])/(x[1]-x[0]); + } else abort(morepoints); + return d; +} + +// Return standard cubic spline interpolation as a guide +guide hermite(real[] x, real[] y, splinetype splinetype=null) +{ + int n=x.length; + if(n == 0) return nullpath; + + guide g=(x[0],y[0]); + if(n == 1) return g; + if(n == 2) return g--(x[1],y[1]); + + if(splinetype == null) + splinetype=(x[0] == x[x.length-1] && y[0] == y[y.length-1]) ? + periodic : notaknot; + + real[] dy=splinetype(x,y); + for(int i=1; i < n; ++i) { + pair z=(x[i],y[i]); + real dx=x[i]-x[i-1]; + g=g..controls((x[i-1],y[i-1])+dx*(1,dy[i-1])/3) and (z-dx*(1,dy[i])/3)..z; + } + return g; +} diff --git a/Master/texmf/asymptote/interpolate.asy b/Master/texmf/asymptote/interpolate.asy index a6b820bc7d5..a600c94c125 100644 --- a/Master/texmf/asymptote/interpolate.asy +++ b/Master/texmf/asymptote/interpolate.asy @@ -29,7 +29,7 @@ // Here s is a real function that is constant on (-infinity,a] and [b,infinity). private import math; -import splinetype; +import graph_splinetype; typedef real fhorner(real); diff --git a/Master/texmf/asymptote/labelpath.asy b/Master/texmf/asymptote/labelpath.asy index fe597a1608d..39c1908efd3 100644 --- a/Master/texmf/asymptote/labelpath.asy +++ b/Master/texmf/asymptote/labelpath.asy @@ -12,7 +12,7 @@ void labelpath(frame f, Label L, path g, string justify=Centered, _labelpath(f,L.s,L.size,g,justify,(L.T.x,L.T.y+0.5linewidth(p)),p); return; } - write("warning: labelpath requires -tex latex"); + warning("labelpathlatex","labelpath requires -tex latex"); } void labelpath(picture pic=currentpicture, Label L, path g, diff --git a/Master/texmf/asymptote/labelpath3.asy b/Master/texmf/asymptote/labelpath3.asy index bf2e6d4c084..8483a76a047 100644 --- a/Master/texmf/asymptote/labelpath3.asy +++ b/Master/texmf/asymptote/labelpath3.asy @@ -41,11 +41,11 @@ triple[] nextframe(path3 p, real reltimestart, triple[] start, real return bf[subdiv]; } -surface labelpath(string txt, path3 p, real angle=90, triple optional=O) +surface labelpath(string s, path3 p, real angle=90, triple optional=O) { real Cos=Cos(angle); real Sin=Sin(angle); - path[] text=texpath(Label(txt,(0,0),Align,basealign)); + path[] text=texpath(Label(s,(0,0),Align,basealign)); text=scale(1/(max(text).x-min(text).x))*text; path[][] decompose=containmentTree(text); diff --git a/Master/texmf/asymptote/obj.asy b/Master/texmf/asymptote/obj.asy index 02f839b594d..158d417c136 100644 --- a/Master/texmf/asymptote/obj.asy +++ b/Master/texmf/asymptote/obj.asy @@ -21,7 +21,7 @@ struct obj { pen[] meshpen; path3[][] read(string datafile, bool verbose=false) { - file in=word(line(input(datafile))); + file in=input(datafile).word().line(); triple[] vert; path3[][] g; g[0]=new path3[] ; @@ -92,8 +92,8 @@ struct obj { material surfacepen, pen meshpen=nullpen) { material[] surfacepen={surfacepen}; pen[] meshpen={meshpen}; - surfacepen.cyclic(true); - meshpen.cyclic(true); + surfacepen.cyclic=true; + meshpen.cyclic=true; operator init(read(datafile,verbose),surfacepen,meshpen); } } diff --git a/Master/texmf/asymptote/plain.asy b/Master/texmf/asymptote/plain.asy index 8424c767ffa..f9e18b883ea 100644 --- a/Master/texmf/asymptote/plain.asy +++ b/Master/texmf/asymptote/plain.asy @@ -20,8 +20,8 @@ include plain_constants; access version; if(version.VERSION != VERSION) { - write(stdout,"warning: using possibly incompatible version "+ - version.VERSION+" of plain.asy"+'\n'); + warning("version","using possibly incompatible version "+ + version.VERSION+" of plain.asy"+'\n'); } include plain_pens; @@ -68,7 +68,7 @@ saveFunction[] saveFunctions={}; // When save is called, this will be redefined to do the corresponding restore. void restore() { - write("warning: restore called with no matching save"); + warning("nomatchingsave","restore called with no matching save"); } void addSaveFunction(saveFunction s) @@ -113,7 +113,8 @@ restoreThunk save() void restoredefaults() { - write("warning: restoredefaults called with no matching savedefaults"); + warning("nomatchingsavedefaults", + "restoredefaults called with no matching savedefaults"); } restoreThunk buildRestoreDefaults() diff --git a/Master/texmf/asymptote/plain_Label.asy b/Master/texmf/asymptote/plain_Label.asy index 3c11a2cedbd..b75d0de3967 100644 --- a/Master/texmf/asymptote/plain_Label.asy +++ b/Master/texmf/asymptote/plain_Label.asy @@ -335,7 +335,13 @@ struct Label { pic.add(new void (frame f, transform t) { out(f,t,point(g,position), alignrelative ? -Align*dir(t*g,position)*I : Align); - },true); + },!alignrelative); + + frame f; + pair align=alignrelative ? -Align*dir(g,position)*I : Align; + label(f,(0,0),align); + pair position=point(g,position); + pic.addBox(position,position,min(f),max(f)); } void write(file file=stdout, suffix suffix=endl) { diff --git a/Master/texmf/asymptote/plain_arrows.asy b/Master/texmf/asymptote/plain_arrows.asy index 0ef8c2b2d5c..18618c2b5b8 100644 --- a/Master/texmf/asymptote/plain_arrows.asy +++ b/Master/texmf/asymptote/plain_arrows.asy @@ -543,27 +543,36 @@ void arrow(picture pic=currentpicture, Label L="", pair b, pair dir, draw(b,pic,L,a--(0,0),align,p,arrow,margin); } -// Force an array of 2D pictures to be as least as large as picture all. -void rescale2(picture[] pictures, picture all) -{ - if(!all.empty2()) { - transform t=inverse(all.calculateTransform()*pictures[0].T); - pair m=t*min(all); - pair M=t*max(all); - for(int i=0; i < pictures.length; ++i) { - draw(pictures[i],m,nullpen); - draw(pictures[i],M,nullpen); - } - } -} - -// Force an array of pictures to have a uniform scaling. -void rescale(picture[] pictures) -{ - if(pictures.length == 0) return; +// Fit an array of pictures simultaneously using the sizing of picture all. +frame[] fit2(picture[] pictures, picture all) +{ + frame[] out; + if(!all.empty2()) { + transform t=all.calculateTransform(); + pair m=all.min(t); + pair M=all.max(t); + for(picture pic : pictures) { + frame f=pic.fit(t); + draw(f,m,nullpen); + draw(f,M,nullpen); + out.push(f); + } + } + return out; +} + +// Fit an array of pictures simultaneously using the size of the first picture. +frame[] fit(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); - rescale2(pictures,all); + + return fit2(pictures,all); } diff --git a/Master/texmf/asymptote/plain_boxes.asy b/Master/texmf/asymptote/plain_boxes.asy index 7fb82c37c1f..ff82a0b0b32 100644 --- a/Master/texmf/asymptote/plain_boxes.asy +++ b/Master/texmf/asymptote/plain_boxes.asy @@ -126,11 +126,11 @@ pair point(object F, pair dir, transform t=identity()) return z+point(F.g,T[0]); } -frame bbox(string prefix="", picture pic=currentpicture, +frame bbox(picture pic=currentpicture, real xmargin=0, real ymargin=xmargin, pen p=currentpen, filltype filltype=NoFill) { - frame f=pic.fit(prefix,max(pic.xsize-2*xmargin,0),max(pic.ysize-2*ymargin,0)); + frame f=pic.fit(max(pic.xsize-2*xmargin,0),max(pic.ysize-2*ymargin,0)); box(f,xmargin,ymargin,p,filltype,above=false); return f; } diff --git a/Master/texmf/asymptote/plain_constants.asy b/Master/texmf/asymptote/plain_constants.asy index c1d8c086f71..32138756891 100644 --- a/Master/texmf/asymptote/plain_constants.asy +++ b/Master/texmf/asymptote/plain_constants.asy @@ -67,6 +67,13 @@ void tab(file file) {write(file,'\t');} void comma(file file) {write(file,',');} typedef void suffix(file); +// Used by interactive write to warn that the outputted type is the resolution +// of an overloaded name. +void overloadedMessage(file file) { + write(file,' <overloaded>'); + endl(file); +} + void write(suffix suffix=endl) {suffix(stdout);} void write(file file, suffix suffix=none) {suffix(file);} diff --git a/Master/texmf/asymptote/plain_paths.asy b/Master/texmf/asymptote/plain_paths.asy index d457fa5bc28..bf3f2790c1d 100644 --- a/Master/texmf/asymptote/plain_paths.asy +++ b/Master/texmf/asymptote/plain_paths.asy @@ -178,7 +178,7 @@ slice cut(path p, path knife, int n) slice s; real[][] T=intersections(p,knife); if(T.length == 0) {s.before=p; s.after=nullpath; return s;} - T.cyclic(true); + T.cyclic=true; real t=T[n][0]; s.before=subpath(p,0,t); s.after=subpath(p,t,length(p)); diff --git a/Master/texmf/asymptote/plain_pens.asy b/Master/texmf/asymptote/plain_pens.asy index fd2dfadde68..a8a0d0c7461 100644 --- a/Master/texmf/asymptote/plain_pens.asy +++ b/Master/texmf/asymptote/plain_pens.asy @@ -139,11 +139,11 @@ pen[] colorPen={red,blue,green,magenta,cyan,orange,purple,brown, deepblue,deepgreen,chartreuse,fuchsia,lightred, lightblue,black,pink,yellow,gray}; -colorPen.cyclic(true); +colorPen.cyclic=true; pen[] monoPen={solid,dashed,dotted,longdashed,dashdotted, longdashdotted}; -monoPen.cyclic(true); +monoPen.cyclic=true; pen Pen(int n) { diff --git a/Master/texmf/asymptote/plain_picture.asy b/Master/texmf/asymptote/plain_picture.asy index b3836d7ef51..432a4a6ae6d 100644 --- a/Master/texmf/asymptote/plain_picture.asy +++ b/Master/texmf/asymptote/plain_picture.asy @@ -401,10 +401,18 @@ struct transformation { transform3 compute() { return projection*modelview; } + transformation copy() { + transformation T=new transformation; + T.modelview=copy(modelview); + T.projection=copy(projection); + T.infinity=infinity; + T.oblique=oblique; + return T; + } } struct projection { - transform3 t; // projection*modelview (cached) + transform3 t; // projection*modelview (cached) bool infinity; bool oblique; bool absolute=false; @@ -420,21 +428,16 @@ struct projection { bool autoadjust=true; // Adjust camera to lie outside bounding volume? bool center=false; // Center target within bounding volume? int ninterpolate; // Used for projecting nurbs to 2D Bezier curves. + transformation T; void calculate() { - transformation T=projector(camera,up,target); + T=projector(camera,up,target); t=T.compute(); infinity=T.infinity; oblique=T.oblique; ninterpolate=infinity ? 1 : 16; } - transformation transformation() { - return projector(camera,up,target); - } - - transform3 modelview() {return transformation().modelview;} - triple vector() { return camera-target; } @@ -473,6 +476,7 @@ struct projection { P.center=center; P.projector=projector; P.ninterpolate=ninterpolate; + P.T=T.copy(); return P; } @@ -500,6 +504,71 @@ struct projection { projection currentprojection; +struct light { + real[][] diffuse; + real[][] ambient; + real[][] specular; + pen background=nullpen; // Background color of the 3D canvas. + real specularfactor; + bool viewport; // Are the lights specified (and fixed) in the viewport frame? + triple[] position; // Only directional lights are currently implemented. + + transform3 T=identity(4); // Transform to apply to normal vectors. + + bool on() {return position.length > 0;} + + void operator init(pen[] diffuse, + pen[] ambient=array(diffuse.length,black), + pen[] specular=diffuse, pen background=nullpen, + real specularfactor=1, + bool viewport=false, triple[] position) { + int n=diffuse.length; + assert(ambient.length == n && specular.length == n && position.length == n); + + this.diffuse=new real[n][]; + this.ambient=new real[n][]; + this.specular=new real[n][]; + this.background=background; + this.position=new triple[n]; + for(int i=0; i < position.length; ++i) { + this.diffuse[i]=rgba(diffuse[i]); + this.ambient[i]=rgba(ambient[i]); + this.specular[i]=rgba(specular[i]); + this.position[i]=unit(position[i]); + } + this.specularfactor=specularfactor; + this.viewport=viewport; + } + + void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse, + pen background=nullpen, real specularfactor=1, + bool viewport=false...triple[] position) { + int n=position.length; + operator init(array(n,diffuse),array(n,ambient),array(n,specular), + background,specularfactor,viewport,position); + } + + void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse, + pen background=nullpen, bool viewport=false, + real x, real y, real z) { + operator init(diffuse,ambient,specular,background,viewport,(x,y,z)); + } + + void operator init(explicit light light) { + diffuse=copy(light.diffuse); + ambient=copy(light.ambient); + specular=copy(light.specular); + background=light.background; + specularfactor=light.specularfactor; + viewport=light.viewport; + position=copy(light.position); + } + + real[] background() {return rgba(background == nullpen ? white : background);} +} + +light currentlight; + triple min3(pen p) { return linewidth(p)*(-0.5,-0.5,-0.5); @@ -990,7 +1059,7 @@ struct picture { if(status == simplex.problem.OPTIMAL) { return scaling.build(p.a(),p.b()).a; } else if(status == simplex.problem.UNBOUNDED) { - if(warn) write("warning: "+dir+" scaling in picture unbounded"); + if(warn) warning("unbounded",dir+" scaling in picture unbounded"); return 0; } else { if(!warn) return 1; @@ -1001,8 +1070,8 @@ struct picture { abort("unbounded picture"); } if(userzero) return 1; - write("warning: cannot fit picture to "+dir+"size "+(string) size - +"...enlarging..."); + warning("cannotfit","cannot fit picture to "+dir+"size "+(string) size + +"...enlarging..."); return calculateScaling(dir,coords,sqrt(2)*size,warn); } } @@ -1184,11 +1253,11 @@ struct picture { pair d=size(f); static real epsilon=100*realEpsilon; if(d.x > xsize*(1+epsilon)) - write("warning: frame exceeds xlimit: "+(string) d.x+" > "+ - (string) xsize); + warning("xlimit","frame exceeds xlimit: "+(string) d.x+" > "+ + (string) xsize); if(d.y > ysize*(1+epsilon)) - write("warning: frame exceeds ylimit: "+(string) d.y+" > "+ - (string) ysize); + warning("ylimit","frame exceeds ylimit: "+(string) d.y+" > "+ + (string) ysize); return f; } @@ -1290,14 +1359,15 @@ struct picture { } static frame fitter(string,picture,string,real,real,bool,bool,string,string, - projection); + light,projection); frame fit(string prefix="", string format="", real xsize=this.xsize, real ysize=this.ysize, bool keepAspect=this.keepAspect, bool view=false, - string options="", string script="", + string options="", string script="", light light=currentlight, projection P=currentprojection) { return fitter == null ? fit2(xsize,ysize,keepAspect) : - fitter(prefix,this,format,xsize,ysize,keepAspect,view,options,script,P); + fitter(prefix,this,format,xsize,ysize,keepAspect,view,options,script, + light,P); } // In case only an approximate picture size estimate is available, return the diff --git a/Master/texmf/asymptote/plain_shipout.asy b/Master/texmf/asymptote/plain_shipout.asy index 22aeef7121f..2e86822cad0 100644 --- a/Master/texmf/asymptote/plain_shipout.asy +++ b/Master/texmf/asymptote/plain_shipout.asy @@ -27,7 +27,7 @@ typedef frame orientation(frame); orientation orientation=Portrait; // Forward references to functions defined in module three. -object embed3(string, frame, string, string, string, projection); +object embed3(string, frame, string, string, string, light, projection); string Embed(string name, string options="", real width=0, real height=0); string Link(string label, string text, string options=""); @@ -59,10 +59,10 @@ include plain_xasy; void shipout(string prefix=defaultfilename, frame f, string format="", bool wait=false, bool view=true, string options="", string script="", - projection P=currentprojection) + light light=currentlight, projection P=currentprojection) { if(is3D(f)) { - f=enclose(prefix,embed3(prefix,f,format,options,script,P)); + f=enclose(prefix,embed3(prefix,f,format,options,script,light,P)); if(settings.render != 0 && !prc(format)) { shipped=true; return; @@ -88,14 +88,14 @@ void shipout(string prefix=defaultfilename, picture pic=currentpicture, orientation orientation=orientation, string format="", bool wait=false, bool view=true, string options="", string script="", - projection P=currentprojection) + light light=currentlight, projection P=currentprojection) { if(!uptodate()) { bool inlinetex=settings.inlinetex; bool prc=prc(format); if(prc && !pic.empty3()) settings.inlinetex=settings.inlineimage; - frame f=pic.fit(prefix,format,view=view,options,script,P); + frame f=pic.fit(prefix,format,view=view,options,script,light,P); if(!pic.empty2() || settings.render == 0 || prc) shipout(prefix,orientation(f),format,wait,view); settings.inlinetex=inlinetex; diff --git a/Master/texmf/asymptote/slide.asy b/Master/texmf/asymptote/slide.asy index b154c1d4746..b23a352fc5d 100644 --- a/Master/texmf/asymptote/slide.asy +++ b/Master/texmf/asymptote/slide.asy @@ -217,11 +217,10 @@ bool checkposition() return true; } -void step() -{ +void erasestep(int erasenode) { if(!stepping || !allowstepping) return; if(!checkposition()) return; - lastnode.push(currentpicture.nodes.length-1); + lastnode.push(erasenode); nextpage(steppagenumberpen); for(int i=0; i < firstnode.length; ++i) { for(int j=firstnode[i]; j <= lastnode[i]; ++j) { @@ -233,6 +232,12 @@ void step() tex(bulletcolor(newbulletcolor)); } +void step() +{ + // Step without erasing anything. + erasestep(currentpicture.nodes.length-1); +} + void incrementposition(pair z) { currentposition += z; @@ -279,11 +284,13 @@ void remark(bool center=false, string s, pair align=0, pen p=itempen, pair M=tinv*min(f); if(abs(offset.x+M.x) > 1) - write("warning: slide too wide on page "+(string) page+':\n'+(string) s); + warning("slidetoowide","slide too wide on page "+(string) page+':\n'+ + (string) s); if(abs(offset.y+M.y) > 1) { void toohigh() { - write("warning: slide too high on page "+(string) page+':\n'+(string) s); + warning("slidetoohigh","slide too high on page "+(string) page+':\n'+ + (string) s); } if(M.y-m.y < 2) { newslide(); offset=(offset.x,currentposition.y); @@ -326,7 +333,7 @@ void skip(real n=1) } void display(frame[] f, real margin=0, pair align=S, pen p=itempen, - pen figuremattpen=figuremattpen) + pen figuremattpen=figuremattpen, bool final=true) { if(f.length == 0) return; real[] width=new real[f.length]; @@ -336,7 +343,7 @@ void display(frame[] f, real margin=0, pair align=S, pen p=itempen, sum += width[i]; } if(sum > pagewidth) - write("warning: slide too wide on page "+(string) page); + warning("toowide","slide too wide on page "+(string) page); else margin=(pagewidth-sum)/(f.length+1); real pos; frame F; @@ -347,19 +354,22 @@ void display(frame[] f, real margin=0, pair align=S, pen p=itempen, pos += w; } add(F,(0,currentposition.y),align); - real a=0.5(unit(align).y-1); - incrementposition((0,(tinv*(a*(max(F)-min(F))-itemskip*I*lineskip(p)*pt)).y)); + if (final) { + real a=0.5(unit(align).y-1); + incrementposition( + (0, (tinv*(a*(max(F)-min(F))-itemskip*I*lineskip(p)*pt)).y)); + } } void display(frame f, real margin=0, pair align=S, pen p=itempen, - pen figuremattpen=figuremattpen) + pen figuremattpen=figuremattpen, bool final=true) { - display(new frame[] {f},margin,align,p,figuremattpen); + display(new frame[] {f},margin,align,p,figuremattpen, final); } void display(string[] s, real margin=0, string[] captions=new string[], string caption="", pair align=S, pen p=itempen, - pen figuremattpen=figuremattpen) + pen figuremattpen=figuremattpen, bool final=true) { frame[] f=new frame[s.length]; frame F; @@ -374,32 +384,59 @@ void display(string[] s, real margin=0, string[] captions=new string[], if(captions[i] != "") label(f[i],captions[i],point(f[i],S).x+I*y,S); } - display(f,margin,align,p,figuremattpen); + display(f,margin,align,p,figuremattpen, final); if(caption != "") center(caption,p); } void display(string s, string caption="", pair align=S, pen p=itempen, - pen figuremattpen=figuremattpen) + pen figuremattpen=figuremattpen, bool final=true) { - display(new string[] {s},caption,align,p,figuremattpen); + display(new string[] {s},caption,align,p,figuremattpen, final); } void figure(string[] s, string options="", real margin=0, string[] captions=new string[], string caption="", - pair align=S, pen p=itempen, pen figuremattpen=figuremattpen) + pair align=S, pen p=itempen, pen figuremattpen=figuremattpen, + bool final=true) { string[] S; for(int i=0; i < s.length; ++i) { S[i]=graphic(s[i],options); } - display(S,margin,captions,caption,align,itempen,figuremattpen); + display(S,margin,captions,caption,align,itempen,figuremattpen, final); } void figure(string s, string options="", string caption="", pair align=S, - pen p=itempen, pen figuremattpen=figuremattpen) + pen p=itempen, pen figuremattpen=figuremattpen, bool final=true) +{ + figure(new string[] {s},options,caption,align,p,figuremattpen, final); +} + +void multifigure(string[] slist, string options="", string caption="", + pair align=S, pen p=itempen, pen figuremattpen=figuremattpen) +{ + if (stepping) { + int lastnode = currentpicture.nodes.length-1; + for (int i=0; i<slist.length-1; ++i) { + figure(slist[i], options, caption, align, p, figuremattpen, final=false); + erasestep(lastnode); + } + } + figure(slist[slist.length-1], options, caption, align, p, figuremattpen, + final=true); + + if(!firststep) step(); + firststep=false; +} + +void indexedfigure(string prefix, int n, string options="", string caption="", + pair align=S, pen p=itempen, pen figuremattpen=figuremattpen) { - figure(new string[] {s},options,caption,align,p,figuremattpen); + string[] s; + for (int i=0; i<n; ++i) + s.push(prefix+"+"+string(i)); + multifigure(s, options, caption, align, p, figuremattpen); } string[] codefile; diff --git a/Master/texmf/asymptote/solids.asy b/Master/texmf/asymptote/solids.asy index 0b2ffa2a508..db98af3088c 100644 --- a/Master/texmf/asymptote/solids.asy +++ b/Master/texmf/asymptote/solids.asy @@ -212,7 +212,8 @@ struct revolution { // must be recomputed if camera is adjusted path3[] silhouette(int m=64, projection P=currentprojection) { if(is3D()) - write("warning: silhouette routine is intended only for 2d projections"); + warning("2Dsilhouette", + "silhouette routine is intended only for 2d projections"); path3 G,H; int N=size(g); int M=(m == 0) ? N : m; diff --git a/Master/texmf/asymptote/syzygy.asy b/Master/texmf/asymptote/syzygy.asy index dce87805c9d..93889b252e5 100644 --- a/Master/texmf/asymptote/syzygy.asy +++ b/Master/texmf/asymptote/syzygy.asy @@ -669,16 +669,6 @@ Braid apply(Relation r, Braid b, int step, int place) { } // Tableau {{{1 -frame[] fit(picture[] pics) { - frame[] f; - for (int i=0; i<pics.length; ++i) { - frame ff=pics[i].fit(); - //label(ff, (string)i, (10,10)); - //f.push(pics[i].fit()); - f.push(ff); - } - return f; -} // Draw a number of frames in a nice circular arrangement. picture tableau(frame[] cards, bool number=false) { diff --git a/Master/texmf/asymptote/three.asy b/Master/texmf/asymptote/three.asy index 39075193073..1f9392ee094 100644 --- a/Master/texmf/asymptote/three.asy +++ b/Master/texmf/asymptote/three.asy @@ -13,13 +13,11 @@ if(prc0()) { real defaultshininess=0.25; real defaultgranularity=0; -real linegranularity=0.01; -real tubegranularity=0.003; +real linegranularity=0.005; +int linesectors=8; // Number of angular sectors. real dotgranularity=0.0001; -real viewportfactor=1.002; // Factor used to expand orthographic viewport. -real angleprecision=1e-3; // Precision for centering perspective projections. -real anglefactor=max(1.01,1+angleprecision); -// Factor used to expand perspective viewport. +real angleprecision=1e-5; // Precision for centering perspective projections. +real rendermargin=0.02; string defaultembed3Doptions; string defaultembed3Dscript; @@ -137,7 +135,9 @@ triple project(triple u, triple v) triple perp(triple v) { triple u=cross(v,Y); - return (abs(u) > sqrtEpsilon) ? unit(u) : unit(cross(v,Z)); + 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 @@ -757,7 +757,7 @@ real[] theta(triple[] v, real[] alpha, real[] beta, { real[] a,b,c,f,l,psi; int n=alpha.length; - bool cyclic=v.cyclicflag; + bool cyclic=v.cyclic; for(int i=0; i < n; ++i) l[i]=1/length(v[i+1]-v[i]); int i0,in; @@ -766,8 +766,8 @@ real[] theta(triple[] v, real[] alpha, real[] beta, 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); + l.cyclic=true; + psi.cyclic=true; } else { psi[n-1]=0; if(dir0 == O) { @@ -869,9 +869,9 @@ void aim(flatguide3 g, int N) beta[k]=g.Tension[K].in; } if(cyclic) { - v.cyclic(true); - alpha.cyclic(true); - beta.cyclic(true); + v.cyclic=true; + alpha.cyclic=true; + beta.cyclic=true; } else v[n]=g.nodes[(start+n) % N]; int final=(end-1) % N; @@ -883,8 +883,8 @@ void aim(flatguide3 g, int N) real[] theta=theta(v,alpha,beta,d0,d1,g.out[start].gamma,g.in[final].gamma, reference); - v.cyclic(true); - theta.cyclic(true); + 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); @@ -1029,6 +1029,17 @@ 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, projection P=currentprojection) +{ + return path3(p,new triple(pair z) {return invert(z,P);}); +} + // Construct a path from a path3 by applying P to each control point. path path(path3 p, pair P(triple)=xypart) { @@ -1801,6 +1812,16 @@ 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); @@ -1978,6 +1999,38 @@ void addPath(picture pic, path3 g, pen 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) @@ -1988,8 +2041,8 @@ void draw(picture pic=currentpicture, Label L="", path3 g, if(is3D()) { draw(f,G,p,light,null); if(pic != null && size(G) > 0) { - pic.addPoint(min(G,P.t)); - pic.addPoint(max(G,P.t)); + pic.addPoint(min(G,P)); + pic.addPoint(max(G,P)); } } if(pic != null) @@ -2004,7 +2057,7 @@ void draw(picture pic=currentpicture, Label L="", path3 g, addPath(pic,g,q); } -include three_arrows; +include three_tube; draw=new void(frame f, path3 g, material p=currentpen, light light=nolight, projection P=currentprojection) { @@ -2015,7 +2068,7 @@ draw=new void(frame f, path3 g, material p=currentpen, if(settings.thick) { real width=linewidth(q); if(width > 0) { - tube T=tube(g,width,p.granularity); + tube T=tube(g,width,linesectors); int L=length(g); if(L >= 0) { if(!cyclic(g)) { @@ -2050,7 +2103,7 @@ draw=new void(frame f, path3 g, material p=currentpen, else { real[] dash=(real[]) split(type," "); if(sum(dash) > 0) { - dash.cyclic(true); + dash.cyclic=true; real offset=offset(q); real L=arclength(g); int i=0; @@ -2081,6 +2134,8 @@ void draw(picture pic=currentpicture, explicit path3[] g, for(int i=0; i < g.length; ++i) draw(pic,g[i],p,margin,light); } +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, @@ -2321,7 +2376,7 @@ pair viewportmargin(pair lambda) string embed3D(string label="", string text=label, string prefix, frame f, string format="", - real width=0, real height=0, real angle=30, + real width=0, real height=0, string options="", string script="", light light=currentlight, projection P=currentprojection) { @@ -2341,7 +2396,7 @@ string embed3D(string label="", string text=label, string prefix, pair margin=viewportmargin((lambda.x,lambda.y)); viewplanesize=(max(lambda.x+2*margin.x,lambda.y+2*margin.y))/(cm*P.zoom); } else - if(!P.absolute) angle=2*aTan(Tan(0.5*angle)); + if(!P.absolute) P.angle=2*aTan(Tan(0.5*P.angle)); string name=prefix+".js"; writeJavaScript(name,lightscript+projection(P.infinity,viewplanesize),script); @@ -2375,7 +2430,7 @@ string embed3D(string label="", string text=label, string prefix, options3 += ",poster"; options3 += ",text={"+text+"},label="+label+ ",toolbar="+(settings.toolbar ? "true" : "false")+ - ",3Daac="+Format(P.absolute ? P.angle : angle)+ + ",3Daac="+Format(P.angle)+ ",3Dc2c="+Format(u)+ ",3Dcoo="+Format(target/cm)+ ",3Droll="+Format(roll)+ @@ -2387,302 +2442,350 @@ string embed3D(string label="", string text=label, string prefix, return Embed(stripdirectory(prefix),options3,width,height); } -object embed(string label="", string text=label, - string prefix=defaultfilename, - frame f, string format="", - real width=0, real height=0, real angle=30, - string options="", string script="", - light light=currentlight, projection P=currentprojection) -{ - object F; - - if(is3D(format)) - F.L=embed3D(label,text,prefix,f,format,width,height,angle,options,script, - light,P); - else - F.f=f; - 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="", real angle=0, - light light=currentlight, projection P=currentprojection) +struct scene { - object F; - real xsize3=pic.xsize3, ysize3=pic.ysize3, zsize3=pic.zsize3; - bool warn=true; - transform3 modelview; - - if(xsize3 == 0 && ysize3 == 0 && zsize3 == 0) { - xsize3=ysize3=zsize3=max(xsize,ysize); - warn=false; + frame f; + transform3 t; + projection P; + bool adjusted; + real width,height; + transform3 T=identity4; + picture pic2; + + void operator init(frame f, projection P=currentprojection) { + this.f=f; + this.t=identity4; + this.P=P; } + + 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; + } - projection P=P.copy(); + if(P.absolute) + this.P=P.copy(); + else if(P.showtarget) + draw(pic,P.target,nullpen); - if(!P.absolute && 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); - transform3 t=pic.scaling(xsize3,ysize3,zsize3,keepAspect,warn); - bool adjusted=false; - transform3 tinv=inverse(t); - triple m=pic.min(t); - triple M=pic.max(t); + if(!P.absolute) { + this.P=t*P; + if(this.P.center) { + bool recalculate=false; + triple target=0.5*(m+M); + this.P.target=target; + recalculate=true; + if(recalculate) this.P.calculate(); + } + if(this.P.autoadjust || this.P.infinity) + adjusted=adjusted | this.P.adjust(m,M); + } - if(!P.absolute) { - P=t*P; - if(P.center) { - bool recalculate=false; - triple target=0.5*(m+M); - P.target=target; - recalculate=true; - if(recalculate) P.calculate(); + f=pic.fit3(t,pic.bounds3.exact ? pic2 : null,this.P); + + if(!pic.bounds3.exact) { + 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(P.autoadjust || P.infinity) - adjusted=adjusted | P.adjust(m,M); - } - picture pic2; - - frame f=pic.fit3(t,pic.bounds3.exact ? pic2 : null,P); + bool scale=xsize != 0 || ysize != 0; + + 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; + pair viewportmargin=viewportmargin(lambda); + width=ceil(lambda.x+2*viewportmargin.x); + height=ceil(lambda.y+2*viewportmargin.y); + + if(!this.P.absolute) { + if(scale && this.P.autoadjust) { + 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); + } + pic2.erase(); + 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; + f=pic.fit3(t,is3D ? null : pic2,this.P); + } - if(!pic.bounds3.exact) { - transform3 s=pic.scale3(f,xsize3,ysize3,zsize3,keepAspect); - t=s*t; - tinv=inverse(t); - P=s*P; - f=pic.fit3(t,pic2,P); + if(this.P.autoadjust || this.P.infinity) + adjusted=adjusted | this.P.adjust(min3(f),max3(f)); + } + } } - bool is3D=is3D(format); - bool scale=xsize != 0 || ysize != 0; - - 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; - pair viewportmargin=viewportmargin(lambda); - real width=ceil(lambda.x+2*viewportmargin.x); - real height=ceil(lambda.y+2*viewportmargin.y); - - projection Q; - if(!P.absolute) { - if(scale && P.autoadjust) { - pair v=(s.xx,s.yy); - transform3 T=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); + // Choose the angle to be just large enough to view the entire image. + real angle(projection P) { + T=identity4; + int maxiterations=100; + 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; } - pic2.erase(); - 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(P.target)*s*shift(-P.target); - t=s*t; - P=s*P; - f=pic.fit3(t,is3D ? null : pic2,P); } + diff=abs(s-lasts); + ++i; + } while (diff > angleprecision && i < maxiterations); + 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)); + } +} - if(P.autoadjust || P.infinity) - adjusted=adjusted | P.adjust(min3(f),max3(f)); +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); - triple target=P.target; - modelview=P.modelview(); - f=modelview*f; - P=modelview*P; - Q=P.copy(); - light=modelview*light; - - if(P.infinity) { - triple m=min3(f); - triple M=max3(f); - triple lambda=M-m; - viewportmargin=viewportmargin((lambda.x,lambda.y)); - width=lambda.x+2*viewportmargin.x; - height=lambda.y+2*viewportmargin.y; - f=shift((-0.5(m.x+M.x),-0.5*(m.y+M.y),0))*f; // Eye will be at (0,0,0). - } else { - transform3 T=identity4; - // Choose the angle to be just large enough to view the entire image: - if(angle == 0) angle=P.angle; - int maxiterations=100; - if(is3D && angle == 0) { - 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 < maxiterations); - 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; - } - - angle=anglefactor*max(aTan(rx)+aTan(Rx),aTan(ry)+aTan(Ry)); - if(viewportmargin.y != 0) - angle=2*aTan(Tan(0.5*angle)-viewportmargin.y/P.target.z); - - modelview=T*modelview; - } - if(settings.verbose > 0) { - transform3 inv=inverse(modelview); - if(adjusted) - write("adjusting camera to ",tinv*inv*P.camera); - target=inv*P.target; - } - P=T*P; + 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; + + pair viewportmargin=viewportmargin; + if(P.infinity) { + triple m=min3(S.f); + triple M=max3(S.f); + triple lambda=M-m; + viewportmargin=viewportmargin((lambda.x,lambda.y)); + S.width=lambda.x+2*viewportmargin.x; + S.height=lambda.y+2*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) { + Q.angle=P.angle=S.angle(P); + modelview=S.T*modelview; + if(viewportmargin.y != 0) + P.angle=2*aTan(Tan(0.5*P.angle)-viewportmargin.y/P.target.z); } if(settings.verbose > 0) { - if(P.center || (!P.infinity && P.autoadjust)) - write("adjusting target to ",tinv*target); + 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 || (!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. - prefix=replace(prefix,new string[][]{{" ","_"},{".","_"}}); - if(settings.embed || nativeformat() == "pdf") - prefix += "+"+(string) file3.length; - } else - preview=false; - if(preview || (!prc && settings.render != 0)) { - frame f=f; - triple m,M; - real zcenter; - if(P.absolute) { - modelview=P.modelview(); - f=modelview*f; - Q=modelview*P; - m=min3(f); - M=max3(f); - real r=0.5*abs(M-m); - zcenter=0.5*(M.z+m.z); - M=(M.x,M.y,zcenter+r); - m=(m.x,m.y,zcenter-r); - angle=P.angle; - } else { - m=min3(f); - M=max3(f); - zcenter=P.target.z; - real d=P.distance(m,M); - M=(M.x,M.y,zcenter+d); - m=(m.x,m.y,zcenter-d); - } - - if(P.infinity) { - triple margin=(viewportfactor-1.0)*(abs(M.x-m.x),abs(M.y-m.y),0) - +(viewportmargin.x,viewportmargin.y,0); - M += margin; - m -= margin; - } else if(M.z >= 0) abort("camera too close"); - - shipout3(prefix,f,preview ? nativeformat() : format, - width,height,P.infinity ? 0 : 2aTan(Tan(0.5*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; + 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. + prefix=replace(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=(viewportmargin.x,viewportmargin.y,0); + M += margin; + m -= margin; + } else if(M.z >= 0) abort("camera too close"); + + shipout3(prefix,f,preview ? nativeformat() : format, + S.width-rendermargin,S.height-rendermargin, + 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 && settings.embed) { + image=prefix; + if(settings.inlinetex) image += "_0"; + 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; +} - string image; - if(preview && settings.embed) { - image=prefix; - if(settings.inlinetex) image += "_0"; - image += "."+nativeformat(); - if(!settings.inlinetex) file3.push(image); - image=graphic(image,"hiresbb"); - } - if(prc) { - if(!P.infinity && P.viewportshift != 0) - write("warning: PRC does not support off-axis projections; use pan instead of shift"); - F.L=embed3D(label,text=image,prefix,f,format, - width,height,angle,options,script,light,Q); - } - +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(T)); + return F; } +} - if(!is3D) { - transform T=pic2.scaling(xsize,ysize,keepAspect); - F.f=pic.fit(scale(t[0][0])*T); - add(F.f,pic2.fit(T)); +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,P),format,view,options,script,light); + else { + object F; + F.f=f; + return F; } - - return F; } embed3=new object(string prefix, frame f, string format, string options, - string script, projection P) { - return embed(prefix=prefix,f,format,options,script,P); + 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, projection P) { + 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)); - if(!empty3) { - 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=embed(prefix=prefix,pic,format,xsize,ysize,keepAspect,view, - options,script,currentlight,P); - if(prc) - label(f,F.L); - else { - if(settings.render == 0) { - add(f,F.f); - if(currentlight.background != nullpen) - box(f,currentlight.background,Fill,above=false); - } else if(!view) - label(f,graphic(prefix,"hiresbb")); - } - } return f; }; +frame embedder(string prefix, frame f, string format, bool view, + string options, string script, light light, projection P) +{ + return embedder(new object(string prefix, string format) { + return embed(prefix=prefix,f,format,view,options,script,light,P); + },prefix,format,view,light); +} + void addViews(picture dest, picture src, bool group=true, filltype filltype=NoFill) { @@ -2733,29 +2836,53 @@ void addAllViews(picture src, bool group=true, filltype filltype=NoFill) addAllViews(currentpicture,src,group,filltype); } -// Force an array of 3D pictures to be as least as large as picture all. -void rescale3(picture[] pictures, picture all, projection P=currentprojection) +// 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) { - if(!all.empty3()) { - transform3 t=inverse(all.calculateTransform3(P)*pictures[0].T3); - triple m=t*min3(all); - triple M=t*max3(all); + 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 loop=settings.loop; + for(picture pic : pictures) { + picture pic2; + frame f=pic.fit3(S.t,pic2,S.P); + if(loop && !settings.loop) break; + add(f,pic2.fit2()); + draw(f,m,nullpen); + draw(f,M,nullpen); + out[i]=loop ? f : embedder(prefix,f,format,view,options,script,light,S.P); + ++i; + } + + while(settings.loop) for(int i=0; i < pictures.length; ++i) { - draw(pictures[i],m,nullpen); - draw(pictures[i],M,nullpen); + if(!settings.loop) break; + embedder(prefix,out[i],format,view,options,script,light,S.P); } - } + + return out; } -// Force an array of pictures to have a uniform scaling using currenprojection. -rescale=new void(picture[] pictures) { - if(pictures.length == 0) return; +// 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); - rescale2(pictures,all); - rescale3(pictures,all); + + return all.empty3() ? fit2(pictures,all) : + fit3(prefix,pictures,all,format,view,options,script,P); }; exitfcn currentexitfunction=atexit(); diff --git a/Master/texmf/asymptote/three_arrows.asy b/Master/texmf/asymptote/three_arrows.asy index 841c4f731f4..161a716bca6 100644 --- a/Master/texmf/asymptote/three_arrows.asy +++ b/Master/texmf/asymptote/three_arrows.asy @@ -48,91 +48,20 @@ void bend(surface s, path3 g, real L) } } -void render(path3 s, real granularity=tubegranularity, void f(path3, real)) -{ - static int maxdepth=ceil(-log(realEpsilon)/log(2))+1; - 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); - real R=infinity; - int nintervals=3; - real fuzz=sqrtEpsilon*max(abs(z0),abs(z1)); - if(S > 0) { - --depth; - for(int i=0; i <= nintervals; ++i) { - R=min(R,radius(z0,c0,c1,z1,i/nintervals)); - if(S > max(granularity*R,fuzz)) { - 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 tube -{ - surface s; - path3 center; - - void operator init(path3 g, real width, real granularity=tubegranularity) { - real r=0.5*width; - - transform3 t=scale3(r); - - int n=length(g); - for(int i=0; i < n; ++i) { - path3 gi=subpath(g,i,i+1); - real S=straightness(g,i); - if(S < sqrtEpsilon*r) { - triple v=point(g,i); - triple u=point(g,i+1)-v; - s.append(shift(v)*align(unit(u))*scale(r,r,abs(u))*unitcylinder); - center=center&gi; - } else { - render(gi,granularity,new void(path3 q, real) { - real L=arclength(q); - surface segment=scale(r,r,L)*unitcylinder; - bend(segment,q,L); - s.s.append(segment.s); - triple a=L*Z/3; - center=center&bend(O,q,L)..controls bend(a,q,L) and bend(2a,q,L).. - bend(3a,q,L); - }); - } - - if((cyclic(g) || i > 0) && abs(dir(g,i,1)-dir(g,i,-1)) > sqrtEpsilon) - s.append(shift(point(g,i))*t*align(dir(g,i,-1))*unithemisphere); - } - } -} - // Refine a noncyclic path3 g so that it approaches its endpoint in // geometrically spaced steps. -path3 approach(path3 g, int n, real radix=3) +path3 approach(path3 p, int n, real radix=3) { guide3 G; - real L=length(g); + real L=length(p); real tlast=0; real r=1/radix; for(int i=1; i < n; ++i) { real t=L*(1-r^i); - G=G&subpath(g,tlast,t); + G=G&subpath(p,tlast,t); tlast=t; } - return G&subpath(g,tlast,L); + return G&subpath(p,tlast,L); } struct arrowhead3 @@ -176,7 +105,7 @@ struct arrowhead3 } if(draw) for(path3 g : H) - s.append(tube(g,width).s); + s.append(tube(g,width,linesectors).s); return shift(v)*s; } @@ -231,8 +160,7 @@ DefaultHead3.head=new surface(path3 g, position position=EndPoint, for(int i=0; i < 4; ++i) { for(int j=0; j < 4; ++j) { real k=1-p.P[i][j].z/remainL; - p.P[i][j]=(k*p.P[i][j].x,k*p.P[i][j].y,p.P[i][j].z); - p.P[i][j]=bend(p.P[i][j],q,l); + p.P[i][j]=bend((k*p.P[i][j].x,k*p.P[i][j].y,p.P[i][j].z),q,l); } } } diff --git a/Master/texmf/asymptote/three_light.asy b/Master/texmf/asymptote/three_light.asy index 8e5a6281795..36ba9a34c19 100644 --- a/Master/texmf/asymptote/three_light.asy +++ b/Master/texmf/asymptote/three_light.asy @@ -75,88 +75,27 @@ material emissive(material m, real granularity=m.granularity) granularity); } -struct light { - real[][] diffuse; - real[][] ambient; - real[][] specular; - pen background; // Background color of the 3D canvas. - real specularfactor; - bool viewport; // Are the lights specified (and fixed) in the viewport frame? - triple[] position; // Only directional lights are currently implemented. - - transform3 T=identity(4); // Transform to apply to normal vectors. - - bool on() {return position.length > 0;} - - void operator init(pen[] diffuse, - pen[] ambient=array(diffuse.length,black), - pen[] specular=diffuse, pen background=nullpen, - real specularfactor=1, - bool viewport=true, triple[] position) { - int n=diffuse.length; - assert(ambient.length == n && specular.length == n && position.length == n); - - this.diffuse=new real[n][]; - this.ambient=new real[n][]; - this.specular=new real[n][]; - this.background=background; - this.position=new triple[n]; - for(int i=0; i < position.length; ++i) { - this.diffuse[i]=rgba(diffuse[i]); - this.ambient[i]=rgba(ambient[i]); - this.specular[i]=rgba(specular[i]); - this.position[i]=unit(position[i]); - } - this.specularfactor=specularfactor; - this.viewport=viewport; - } - - void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse, - pen background=nullpen, real specularfactor=1, - bool viewport=true...triple[] position) { - int n=position.length; - operator init(array(n,diffuse),array(n,ambient),array(n,specular), - background,specularfactor,viewport,position); - } - - void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse, - pen background=nullpen, bool viewport=true, - real x, real y, real z) { - operator init(diffuse,ambient,specular,background,viewport,(x,y,z)); +pen color(triple normal, material m, light light, transform3 T=light.T) { + triple[] position=light.position; + if(invisible((pen) m)) return invisible; + if(position.length == 0) return m.diffuse(); + normal=unit(T*normal); + if(settings.twosided) normal *= sgn(normal.z); + real s=m.shininess*128; + real[] Diffuse=rgba(m.diffuse()); + real[] Ambient=rgba(m.ambient()); + real[] Specular=rgba(m.specular()); + real[] p=rgba(m.emissive()); + for(int i=0; i < position.length; ++i) { + triple L=light.viewport ? position[i] : T*position[i]; + real Ldotn=max(dot(normal,L),0); + p += light.ambient[i]*Ambient+Ldotn*light.diffuse[i]*Diffuse; + // Apply specularfactor to partially compensate non-pixel-based rendering. + if(Ldotn > 0) // Phong-Blinn model of specular reflection + p += dot(normal,unit(L+Z))^s*light.specularfactor* + light.specular[i]*Specular; } - - void operator init(explicit light light) { - diffuse=copy(light.diffuse); - ambient=copy(light.ambient); - specular=copy(light.specular); - background=light.background; - specularfactor=light.specularfactor; - viewport=light.viewport; - position=copy(light.position); - } - - pen color(triple normal, material m, transform3 T=T) { - if(invisible((pen) m)) return invisible; - if(position.length == 0) return m.diffuse(); - normal=unit(T*normal); - if(settings.twosided) normal *= sgn(normal.z); - real s=m.shininess*128; - real[] Diffuse=rgba(m.diffuse()); - real[] Ambient=rgba(m.ambient()); - real[] Specular=rgba(m.specular()); - real[] p=rgba(m.emissive()); - for(int i=0; i < position.length; ++i) { - triple L=viewport ? position[i] : T*position[i]; - real Ldotn=max(dot(normal,L),0); - p += ambient[i]*Ambient+Ldotn*diffuse[i]*Diffuse; - // Apply specularfactor to partially compensate non-pixel-based rendering. - if(Ldotn > 0) // Phong-Blinn model of specular reflection - p += dot(normal,unit(L+Z))^s*specularfactor*specular[i]*Specular; - } - return rgb(p[0],p[1],p[2])+opacity(opacity(m.diffuse())); - } - - real[] background() {return rgba(background == nullpen ? white : background);} + return rgb(p[0],p[1],p[2])+opacity(opacity(m.diffuse())); } light operator * (transform3 t, light light) @@ -168,14 +107,16 @@ light operator * (transform3 t, light light) light operator cast(triple v) {return light(v);} -light currentlight=light(ambient=gray(0.1),specularfactor=3, - (0.25,-0.25,1)); +light Viewport=light(ambient=gray(0.1),specularfactor=3,viewport=true, + (0.25,-0.25,1)); light White=light(new pen[] {rgb(0.38,0.38,0.45),rgb(0.6,0.6,0.67), - rgb(0.5,0.5,0.57)},specularfactor=3,viewport=false, + rgb(0.5,0.5,0.57)},specularfactor=3, new triple[] {(-2,-1.5,-0.5),(2,1.1,-2.5),(-0.5,0,2)}); -light Headlamp=light(gray(0.9),ambient=gray(0.1),specularfactor=3, - (0.5,0.5,1/sqrt(2)),specular=gray(0.7)); +light Headlamp=light(gray(0.8),ambient=gray(0.1),specular=gray(0.7), + specularfactor=3,viewport=true,dir(42,48)); + +currentlight=Headlamp; light nolight; diff --git a/Master/texmf/asymptote/three_surface.asy b/Master/texmf/asymptote/three_surface.asy index f2ad86dfe7e..e545d87cf7b 100644 --- a/Master/texmf/asymptote/three_surface.asy +++ b/Master/texmf/asymptote/three_surface.asy @@ -1,4 +1,5 @@ import bezulate; +private import interpolate; int nslice=12; real camerafactor=1.2; @@ -15,24 +16,24 @@ struct patch { path3 external() { return - P[0][0]..controls P[0][1] and P[0][2].. - P[0][3]..controls P[1][3] and P[2][3].. - P[3][3]..controls P[3][2] and P[3][1].. - P[3][0]..controls P[2][0] and P[1][0]..cycle; + P[0][0]..controls P[1][0] and P[2][0].. + P[3][0]..controls P[3][1] and P[3][2].. + P[3][3]..controls P[2][3] and P[1][3].. + P[0][3]..controls P[0][2] and P[0][1]..cycle; } triple[] internal() { - return new triple[] {P[1][1],P[1][2],P[2][2],P[2][1]}; + return new triple[] {P[1][1],P[2][1],P[2][2],P[1][2]}; } triple cornermean() { - return 0.25*(P[0][0]+P[0][3]+P[3][3]+P[3][0]); + return 0.25*(P[0][0]+P[0][3]+P[3][0]+P[3][3]); } - triple[] corners() {return new triple[] {P[0][0],P[0][3],P[3][3],P[3][0]};} + triple[] corners() {return new triple[] {P[0][0],P[3][0],P[3][3],P[0][3]};} real[] map(real f(triple)) { - return new real[] {f(P[0][0]),f(P[0][3]),f(P[3][3]),f(P[3][0])}; + return new real[] {f(P[0][0]),f(P[3][0]),f(P[3][3]),f(P[0][3])}; } triple Bu(int j, real u) {return bezier(P[0][j],P[1][j],P[2][j],P[3][j],u);} @@ -43,8 +44,10 @@ struct patch { triple BuPPP(int j) {return bezierPPP(P[0][j],P[1][j],P[2][j],P[3][j]);} path3 uequals(real u) { - return straight ? Bu(0,u)--Bu(3,u) : - Bu(0,u)..controls Bu(1,u) and Bu(2,u)..Bu(3,u); + triple z0=Bu(0,u); + triple z1=Bu(3,u); + return path3(new triple[] {z0,Bu(2,u)},new triple[] {z0,z1}, + new triple[] {Bu(1,u),z1},new bool[] {straight,false},false); } triple Bv(int i, real v) {return bezier(P[i][0],P[i][1],P[i][2],P[i][3],v);} @@ -55,8 +58,10 @@ struct patch { triple BvPPP(int i) {return bezierPPP(P[i][0],P[i][1],P[i][2],P[i][3]);} path3 vequals(real v) { - return straight ? Bv(0,v)--Bv(3,v) : - Bv(0,v)..controls Bv(1,v) and Bv(2,v)..Bv(3,v); + triple z0=Bv(0,v); + triple z1=Bv(3,v); + return path3(new triple[] {z0,Bv(2,v)},new triple[] {z0,z1}, + new triple[] {Bv(1,v),z1},new bool[] {straight,false},false); } triple point(real u, real v) { @@ -65,76 +70,76 @@ struct patch { // compute normal vectors for degenerate cases private triple normal0(real u, real v, real epsilon) { - triple n=0.5*(cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u), - bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v))+ - cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u), - bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v))); + triple n=0.5*(cross(bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v), + bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u))+ + cross(bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v), + bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u))); return abs(n) > epsilon ? n : - 0.25*cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u), - bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v))+ - 1/6*(cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u), - bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v))+ - cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u), - bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v)))+ - 1/12*(cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u), - bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v))+ - cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u), - bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v)))+ - 1/36*cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u), - bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v)); + 0.25*cross(bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v), + bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u))+ + 1/6*(cross(bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v), + bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u))+ + cross(bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v), + bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u)))+ + 1/12*(cross(bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v), + bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u))+ + cross(bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v), + bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u)))+ + 1/36*cross(bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v), + bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u)); } static real fuzz=1000*realEpsilon; triple normal(real u, real v) { - triple n=cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u), - bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v)); + triple n=cross(bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v), + bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u)); real epsilon=fuzz*change2(P); return (abs(n) > epsilon) ? n : normal0(u,v,epsilon); } triple normal00() { - triple n=9*cross(P[0][1]-P[0][0],P[1][0]-P[0][0]); + triple n=9*cross(P[1][0]-P[0][0],P[0][1]-P[0][0]); real epsilon=fuzz*change2(P); return abs(n) > epsilon ? n : normal0(0,0,epsilon); } - triple normal01() { - triple n=9*cross(P[0][3]-P[0][2],P[1][3]-P[0][3]); + triple normal10() { + triple n=9*cross(P[3][0]-P[2][0],P[3][1]-P[3][0]); real epsilon=fuzz*change2(P); - return abs(n) > epsilon ? n : normal0(0,1,epsilon); + return abs(n) > epsilon ? n : normal0(1,0,epsilon); } triple normal11() { - triple n=9*cross(P[3][3]-P[3][2],P[3][3]-P[2][3]); + triple n=9*cross(P[3][3]-P[2][3],P[3][3]-P[3][2]); real epsilon=fuzz*change2(P); return abs(n) > epsilon ? n : normal0(1,1,epsilon); } - triple normal10() { - triple n=9*cross(P[3][1]-P[3][0],P[3][0]-P[2][0]); + triple normal01() { + triple n=9*cross(P[1][3]-P[0][3],P[0][3]-P[0][2]); real epsilon=fuzz*change2(P); - return abs(n) > epsilon ? n : normal0(1,0,epsilon); + return abs(n) > epsilon ? n : normal0(0,1,epsilon); } pen[] colors(material m, light light=currentlight) { bool nocolors=colors.length == 0; if(normals.length > 0) - return new pen[] {light.color(normals[0],nocolors ? m : colors[0]), - light.color(normals[1],nocolors ? m : colors[1]), - light.color(normals[2],nocolors ? m : colors[2]), - light.color(normals[3],nocolors ? m : colors[3])}; + return new pen[] {color(normals[0],nocolors ? m : colors[0],light), + color(normals[1],nocolors ? m : colors[1],light), + color(normals[2],nocolors ? m : colors[2],light), + color(normals[3],nocolors ? m : colors[3],light)}; if(planar) { triple normal=normal(0.5,0.5); - return new pen[] {light.color(normal,nocolors ? m : colors[0]), - light.color(normal,nocolors ? m : colors[1]), - light.color(normal,nocolors ? m : colors[2]), - light.color(normal,nocolors ? m : colors[3])}; + return new pen[] {color(normal,nocolors ? m : colors[0],light), + color(normal,nocolors ? m : colors[1],light), + color(normal,nocolors ? m : colors[2],light), + color(normal,nocolors ? m : colors[3],light)}; } - return new pen[] {light.color(normal00(),nocolors ? m : colors[0]), - light.color(normal01(),nocolors ? m : colors[1]), - light.color(normal11(),nocolors ? m : colors[2]), - light.color(normal10(),nocolors ? m : colors[3])}; + return new pen[] {color(normal00(),nocolors ? m : colors[0],light), + color(normal10(),nocolors ? m : colors[1],light), + color(normal11(),nocolors ? m : colors[2],light), + color(normal01(),nocolors ? m : colors[3],light)}; } triple min3,max3; @@ -162,11 +167,19 @@ struct patch { } pair min(projection P, pair bound=project(this.P[0][0],P.t)) { - return minbezier(this.P,P.t,bound); + triple[][] Q=P.T.modelview*this.P; + if(P.infinity) + return xypart(minbezier(Q,(bound.x,bound.y,0))); + real d=P.T.projection[3][2]; + return maxratio(Q,d*bound)/d; // d is negative } pair max(projection P, pair bound=project(this.P[0][0],P.t)) { - return maxbezier(this.P,P.t,bound); + triple[][] Q=P.T.modelview*this.P; + if(P.infinity) + return xypart(maxbezier(Q,(bound.x,bound.y,0))); + real d=P.T.projection[3][2]; + return minratio(Q,d*bound)/d; // d is negative } void operator init(triple[][] P, triple[] normals=new triple[], @@ -242,11 +255,11 @@ struct patch { } else straight=false; P=new triple[][] { - {point(external,0),postcontrol(external,0),precontrol(external,1), - point(external,1)}, - {precontrol(external,0),internal[0],internal[1],postcontrol(external,1)}, - {postcontrol(external,3),internal[3],internal[2],precontrol(external,2)}, - {point(external,3),precontrol(external,3),postcontrol(external,2), + {point(external,0),precontrol(external,0),postcontrol(external,3), + point(external,3)}, + {postcontrol(external,0),internal[0],internal[3],precontrol(external,3)}, + {precontrol(external,1),internal[1],internal[2],postcontrol(external,2)}, + {point(external,1),postcontrol(external,1),precontrol(external,2), point(external,2)} }; } @@ -280,10 +293,10 @@ struct patch { delta[j]=(external[(j+1)% 4]-external[j])/3; P=new triple[][] { - {external[0],external[0]+delta[0],external[1]-delta[0],external[1]}, - {external[0]-delta[3],internal[0],internal[1],external[1]+delta[1]}, - {external[3]+delta[3],internal[3],internal[2],external[2]-delta[1]}, - {external[3],external[3]-delta[2],external[2]+delta[2],external[2]} + {external[0],external[0]-delta[3],external[3]+delta[3],external[3]}, + {external[0]+delta[0],internal[0],internal[3],external[3]-delta[2]}, + {external[1]-delta[0],internal[1],internal[2],external[2]+delta[2]}, + {external[1],external[1]+delta[1],external[2]-delta[1],external[2]} }; } } @@ -322,8 +335,243 @@ patch reverse(patch s) return S; } +// Return the Coons patch control points corresponding to path p. +pair[][] coons(path p) +{ + int L=length(p); + if(L == 1) + p=p--cycle--cycle--cycle; + else if(L == 2) + p=p--cycle--cycle; + else if(L == 3) + p=p--cycle; + + pair[] internal=new pair[4]; + for(int j=0; j < 4; ++j) { + internal[j]=nineth*(-4*point(p,j) + +6*(precontrol(p,j)+postcontrol(p,j)) + -2*(point(p,j-1)+point(p,j+1)) + +3*(precontrol(p,j-1)+postcontrol(p,j+1)) + -point(p,j+2)); + } + + return new pair[][] { + {point(p,0),precontrol(p,0),postcontrol(p,3),point(p,3)}, + {postcontrol(p,0),internal[0],internal[3],precontrol(p,3)}, + {precontrol(p,1),internal[1],internal[2],postcontrol(p,2)}, + {point(p,1),postcontrol(p,1),precontrol(p,2),point(p,2)} + }; +} + +// Decompose a possibly nonconvex cyclic path into an array of paths that +// yield nondegenerate Coons patches. +path[] regularize(path p, bool checkboundary=true) +{ + path[] s; + + if(!cyclic(p)) + abort("cyclic path expected"); + + int L=length(p); + + if(L > 4) { + for(path g : bezulate(p)) + s.append(regularize(g,checkboundary)); + return s; + } + + bool straight=piecewisestraight(p); + if(L <= 3 && straight) { + return new path[] {p}; + } + + // Split p along the angle bisector at t. + bool split(path p, real t) { + pair dir=dir(p,t); + if(dir != 0) { + path g=subpath(p,t,t+length(p)); + int L=length(g); + pair z=point(g,0); + real[] T=intersections(g,z,z+I*dir); + for(int i=0; i < T.length; ++i) { + real cut=T[i]; + if(cut > sqrtEpsilon && cut < L-sqrtEpsilon) { + pair w=point(g,cut); + if(!inside(p,0.5*(z+w),zerowinding)) continue; + pair delta=sqrtEpsilon*(w-z); + if(intersections(g,z-delta--w+delta).length != 2) continue; + s.append(regularize(subpath(g,0,cut)--cycle,checkboundary)); + s.append(regularize(subpath(g,cut,L)--cycle,checkboundary)); + return true; + } + } + } + return false; + } + + // Ensure that all interior angles are less than 180 degrees. + real fuzz=1e-4; + int sign=sgn(windingnumber(p,inside(p,zerowinding))); + for(int i=0; i < L; ++i) { + if(sign*(conj(dir(p,i,-1))*dir(p,i,1)).y < -fuzz) { + if(split(p,i)) return s; + } + } + + if(straight) + return new path[] {p}; + + pair[][] P=coons(p); + + // Check for degeneracy. + pair[][] U=new pair[3][4]; + pair[][] V=new pair[4][3]; + + for(int i=0; i < 3; ++i) { + for(int j=0; j < 4; ++j) + U[i][j]=P[i+1][j]-P[i][j]; + } + + for(int i=0; i < 4; ++i) { + for(int j=0; j < 3; ++j) + V[i][j]=P[i][j+1]-P[i][j]; + } + + int[] choose2={1,2,1}; + int[] choose3={1,3,3,1}; + + real T[][]=new real[6][6]; + for(int p=0; p < 6; ++p) { + int kstart=max(p-2,0); + int kstop=min(p,3); + real[] Tp=T[p]; + for(int q=0; q < 6; ++q) { + real Tpq; + int jstop=min(q,3); + int jstart=max(q-2,0); + for(int k=kstart; k <= kstop; ++k) { + int choose3k=choose3[k]; + for(int j=jstart; j <= jstop; ++j) { + int i=p-k; + int l=q-j; + Tpq += (conj(U[i][j])*V[k][l]).y* + choose2[i]*choose3k*choose3[j]*choose2[l]; + } + } + Tp[q]=Tpq; + } + } + + bool3 aligned=default; + bool degenerate=false; + + for(int p=0; p < 6; ++p) { + for(int q=0; q < 6; ++q) { + if(aligned == default) { + if(T[p][q] > sqrtEpsilon) aligned=true; + if(T[p][q] < -sqrtEpsilon) aligned=false; + } else { + if((T[p][q] > sqrtEpsilon && aligned == false) || + (T[p][q] < -sqrtEpsilon && aligned == true)) degenerate=true; + } + } + } + + if(!degenerate) { + if(aligned == (sign >= 0)) + return new path[] {p}; + return s; + } + + if(checkboundary) { + // Polynomial coefficients of (B_i'' B_j + B_i' B_j')/3. + static real[][][] fpv0={ + {{5, -20, 30, -20, 5}, + {-3, 24, -54, 48, -15}, + {0, -6, 27, -36, 15}, + {0, 0, -3, 8, -5}}, + {{-7, 36, -66, 52, -15}, + {3, -36, 108, -120, 45}, + {0, 6, -45, 84, -45}, + {0, 0, 3, -16, 15}}, + {{2, -18, 45, -44, 15}, + {0, 12, -63, 96, -45}, + {0, 0, 18, -60, 45}, + {0, 0, 0, 8, -15}}, + {{0, 2, -9, 12, -5}, + {0, 0, 9, -24, 15}, + {0, 0, 0, 12, -15}, + {0, 0, 0, 0, 5}} + }; + + // Compute one-ninth of the derivative of the Jacobian along the boundary. + real[][] c=array(4,array(5,0.0)); + for(int i=0; i < 4; ++i) { + real[][] fpv0i=fpv0[i]; + for(int j=0; j < 4; ++j) { + real[] w=fpv0i[j]; + c[0] += w*(conj(P[i][0])*(P[j][1]-P[j][0])).y; // v=0 + c[1] += w*(conj(P[3][j]-P[2][j])*P[3][i]).y; // u=1 + c[2] += w*(conj(P[i][3])*(P[j][3]-P[j][2])).y; // v=1 + c[3] += w*(conj(P[0][j]-P[1][j])*P[0][i]).y; // u=0 + } + } + + pair BuP(int j, real u) { + return bezierP(P[0][j],P[1][j],P[2][j],P[3][j],u); + } + pair BvP(int i, real v) { + return bezierP(P[i][0],P[i][1],P[i][2],P[i][3],v); + } + real normal(real u, real v) { + return (conj(bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v))* + bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u)).y; + } + + // Use Rolle's theorem to check for degeneracy on the boundary. + real M=0; + real cut; + for(int i=0; i < 4; ++i) { + if(!straight(p,i)) { + real[] ci=c[i]; + pair[] R=quarticroots(ci[4],ci[3],ci[2],ci[1],ci[0]); + for(pair r : R) { + if(fabs(r.y) < sqrtEpsilon) { + real t=r.x; + if(0 <= t && t <= 1) { + real[] U={t,1,t,0}; + real[] V={0,t,1,t}; + real[] T={t,t,1-t,1-t}; + real N=sign*normal(U[i],V[i]); + if(N < M) { + M=N; cut=i+T[i]; + } + } + } + } + } + } + + // Split at the worst boundary degeneracy. + if(M < 0 && split(p,cut)) return s; + } + + // Split arbitrarily to resolve any remaining (internal) degeneracy. + checkboundary=false; + for(int i=0; i < L; ++i) + if(!straight(p,i) && split(p,i+0.5)) return s; + + while(true) + for(int i=0; i < L; ++i) + if(!straight(p,i) && split(p,i+unitrand())) return s; + + return s; +} + struct surface { patch[] s; + int index[][]; + bool vcyclic; bool empty() { return s.length == 0; @@ -341,6 +589,8 @@ struct surface { this.s=new patch[s.s.length]; for(int i=0; i < s.s.length; ++i) this.s[i]=patch(s.s[i]); + this.index=copy(s.index); + this.vcyclic=s.vcyclic; } void operator init(triple[][][] P, triple[][] normals=new triple[][], @@ -352,10 +602,8 @@ struct surface { } void colors(pen[][] palette) { - for(int i=0; i < s.length; ++i) { - pen[] palettei=palette[i]; - s[i].colors=new pen[] {palettei[0],palettei[1],palettei[2],palettei[3]}; - } + for(int i=0; i < s.length; ++i) + s[i].colors=copy(palette[i]); } triple[][] corners() { @@ -376,232 +624,66 @@ struct surface { return sequence(new triple(int i) {return s[i].cornermean();},s.length); } - // A constructor for a possibly nonconvex cyclic path in a given plane. - void operator init (path p, triple plane(pair)=XYplane, - bool checkboundary=true) { - if(!cyclic(p)) - abort("cyclic path expected"); - - int L=length(p); + triple point(real u, real v) { + int U=floor(u); + int V=floor(v); + int index=index.length == 0 ? U+V : index[U][V]; + return s[index].point(u-U,v-V); + } - if(L > 4) { - for(path g : bezulate(p)) - s.append(surface(g,plane,checkboundary).s); - return; - } - - pair[][] P(path p) { - if(L == 1) - p=p--cycle--cycle--cycle; - else if(L == 2) - p=p--cycle--cycle; - else if(L == 3) - p=p--cycle; - - pair[] internal=new pair[4]; - for(int j=0; j < 4; ++j) { - internal[j]=nineth*(-4*point(p,j) - +6*(precontrol(p,j)+postcontrol(p,j)) - -2*(point(p,j-1)+point(p,j+1)) - +3*(precontrol(p,j-1)+postcontrol(p,j+1)) - -point(p,j+2)); - } - - return new pair[][] { - {point(p,0),postcontrol(p,0),precontrol(p,1),point(p,1)}, - {precontrol(p,0),internal[0],internal[1],postcontrol(p,1)}, - {postcontrol(p,3),internal[3],internal[2],precontrol(p,2)}, - {point(p,3),precontrol(p,3),postcontrol(p,2),point(p,2)} - }; - } + triple normal(real u, real v) { + int U=floor(u); + int V=floor(v); + int index=index.length == 0 ? U+V : index[U][V]; + return s[index].normal(u-U,v-V); + } + + void ucyclic(bool f) + { + index.cyclic=f; + } + + void vcyclic(bool f) + { + for(int[] i : index) + i.cyclic=f; + vcyclic=f; + } + + bool ucyclic() + { + return index.cyclic; + } + + bool vcyclic() + { + return vcyclic; + } + path3 uequals(real u) { + if(index.length == 0) return nullpath3; + int U=floor(u); + int[] index=index[U]; + path3 g; + for(int i : index) + g=g&s[i].uequals(u-U); + return vcyclic() ? g&cycle : g; + } + + path3 vequals(real v) { + if(index.length == 0) return nullpath3; + int V=floor(v); + path3 g; + for(int[] i : index) + g=g&s[i[V]].vequals(v-V); + return ucyclic() ? g&cycle : g; + } + + // A constructor for a possibly nonconvex cyclic path in a given plane. + void operator init(path p, triple plane(pair)=XYplane) { bool straight=piecewisestraight(p); - if(L <= 3 && straight) { - s=new patch[] {patch(P(p),plane,straight)}; - return; - } - - // Split p along the angle bisector at t. - bool split(path p, real t) { - pair dir=dir(p,t); - if(dir != 0) { - path g=subpath(p,t,t+length(p)); - int L=length(g); - pair z=point(g,0); - real[] T=intersections(g,z,z+I*dir); - for(int i=0; i < T.length; ++i) { - real cut=T[i]; - if(cut > sqrtEpsilon && cut < L-sqrtEpsilon) { - pair w=point(g,cut); - if(!inside(p,0.5*(z+w),zerowinding)) continue; - pair delta=sqrtEpsilon*(w-z); - if(intersections(g,z-delta--w+delta).length != 2) continue; - s=surface(subpath(g,0,cut)--cycle,plane,checkboundary).s; - s.append(surface(subpath(g,cut,L)--cycle,plane,checkboundary).s); - return true; - } - } - } - return false; - } - - // Ensure that all interior angles are less than 180 degrees. - real fuzz=1e-4; - int sign=sgn(windingnumber(p,inside(p,zerowinding))); - for(int i=0; i < L; ++i) { - if(sign*(conj(dir(p,i,-1))*dir(p,i,1)).y < -fuzz) { - if(split(p,i)) return; - } - } - - pair[][] P=P(p); - - if(straight) { - s=new patch[] {patch(P,plane,straight)}; - return; - } - - // Check for degeneracy. - pair[][] U=new pair[3][4]; - pair[][] V=new pair[4][3]; - - for(int i=0; i < 3; ++i) { - for(int j=0; j < 4; ++j) - U[i][j]=P[i+1][j]-P[i][j]; - } - - for(int i=0; i < 4; ++i) { - for(int j=0; j < 3; ++j) - V[i][j]=P[i][j+1]-P[i][j]; - } - - int[] choose2={1,2,1}; - int[] choose3={1,3,3,1}; - - real T[][]=new real[6][6]; - for(int p=0; p < 6; ++p) { - int kstart=max(p-2,0); - int kstop=min(p,3); - real[] Tp=T[p]; - for(int q=0; q < 6; ++q) { - real Tpq; - int jstop=min(q,3); - int jstart=max(q-2,0); - for(int k=kstart; k <= kstop; ++k) { - int choose3k=choose3[k]; - for(int j=jstart; j <= jstop; ++j) { - int i=p-k; - int l=q-j; - Tpq += (conj(U[i][j])*V[k][l]).y* - choose2[i]*choose3k*choose3[j]*choose2[l]; - } - } - Tp[q]=Tpq; - } - } - - bool3 aligned=default; - bool degenerate=false; - - for(int p=0; p < 6; ++p) { - for(int q=0; q < 6; ++q) { - if(aligned == default) { - if(T[p][q] < -sqrtEpsilon) aligned=true; - if(T[p][q] > sqrtEpsilon) aligned=false; - } else { - if((T[p][q] < -sqrtEpsilon && aligned == false) || - (T[p][q] > sqrtEpsilon && aligned == true)) degenerate=true; - } - } - } - - if(!degenerate) { - if(aligned == (sign >= 0)) - s=new patch[] {patch(P,plane)}; - return; - } - - if(checkboundary) { - // Polynomial coefficients of (B_i'' B_j + B_i' B_j')/3. - static real[][][] fpv0={ - {{5, -20, 30, -20, 5}, - {-3, 24, -54, 48, -15}, - {0, -6, 27, -36, 15}, - {0, 0, -3, 8, -5}}, - {{-7, 36, -66, 52, -15}, - {3, -36, 108, -120, 45}, - {0, 6, -45, 84, -45}, - {0, 0, 3, -16, 15}}, - {{2, -18, 45, -44, 15}, - {0, 12, -63, 96, -45}, - {0, 0, 18, -60, 45}, - {0, 0, 0, 8, -15}}, - {{0, 2, -9, 12, -5}, - {0, 0, 9, -24, 15}, - {0, 0, 0, 12, -15}, - {0, 0, 0, 0, 5}} - }; - - // Compute one-ninth of the derivative of the Jacobian along the boundary. - real[][] c=array(4,array(5,0.0)); - for(int i=0; i < 4; ++i) { - real[][] fpv0i=fpv0[i]; - for(int j=0; j < 4; ++j) { - real[] w=fpv0i[j]; - c[0] += w*(conj(P[0][j]-P[1][j])*P[0][i]).y; // u=0 - c[1] += w*(conj(P[i][3])*(P[j][3]-P[j][2])).y; // v=1 - c[2] += w*(conj(P[3][j]-P[2][j])*P[3][i]).y; // u=1 - c[3] += w*(conj(P[i][0])*(P[j][1]-P[j][0])).y; // v=0 - } - } - - pair BuP(int j, real u) { - return bezierP(P[0][j],P[1][j],P[2][j],P[3][j],u); - } - pair BvP(int i, real v) { - return bezierP(P[i][0],P[i][1],P[i][2],P[i][3],v); - } - real normal(real u, real v) { - return (conj(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u))* - bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v)).y; - } - - // Use Rolle's theorem to check for degeneracy on the boundary. - real M=0; - real cut; - for(int i=0; i < 4; ++i) { - if(!straight(p,i)) { - real[] ci=c[i]; - pair[] R=quarticroots(ci[4],ci[3],ci[2],ci[1],ci[0]); - for(pair r : R) { - if(fabs(r.y) < sqrtEpsilon) { - real t=r.x; - if(0 <= t && t <= 1) { - real[] U={0,t,1,t}; - real[] V={t,1,t,0}; - real[] T={t,t,1-t,1-t}; - real N=sign*normal(U[i],V[i]); - if(N < M) { - M=N; cut=i+T[i]; - } - } - } - } - } - } - - // Split at the worst boundary degeneracy. - if(M < 0 && split(p,cut)) return; - } - - // Split arbitrarily to resolve any remaining (internal) degeneracy. - checkboundary=false; - for(int i=0; i < L; ++i) - if(!straight(p,i) && split(p,i+0.5)) return; - - while(true) - for(int i=0; i < L; ++i) - if(!straight(p,i) && split(p,i+unitrand())) return; + for(path g : regularize(p)) + s.push(patch(coons(g),plane,straight)); } void operator init(explicit path[] g, triple plane(pair)=XYplane) { @@ -715,12 +797,13 @@ struct surface { // An optional surface pen color(int i, real j) may be specified // to override the color at vertex(i,j). void operator init(triple c, path3 g, triple axis, int n=nslice, - real angle1=0, real angle2= 360, + real angle1=0, real angle2=360, pen color(int i, real j)=null) { axis=unit(axis); real w=(angle2-angle1)/n; int L=length(g); s=new patch[L*n]; + index=new int[n][L]; int m=-1; transform3[] T=new transform3[n+1]; transform3 t=rotate(w,c,c+axis); @@ -728,6 +811,9 @@ struct surface { for(int k=1; k <= n; ++k) T[k]=T[k-1]*t; + typedef pen colorfcn(int i, real j); + bool defaultcolors=(colorfcn) color == null; + for(int i=0; i < L; ++i) { path3 h=subpath(g,i,i+1); path3 r=reverse(h); @@ -754,13 +840,16 @@ struct surface { for(int k=0; k < n; ++k, j += w) { transform3 Tp=T[k+1]; triple dirp=dir(j+w); - path3 G=Tk*h{dirj}..{dirp}Tp*r{-dirp}..{-dirj}cycle; + path3 G=reverse(Tk*h{dirj}..{dirp}Tp*r{-dirp}..{-dirj}cycle); Tk=Tp; dirj=dirp; - s[++m]=color == null ? patch(G) : - patch(G,new pen[] {color(i,j),color(i+1,j),color(i+1,j+w), - color(i,j+w)}); + s[++m]=defaultcolors ? patch(G) : + patch(G,new pen[] {color(i,j),color(i,j+w),color(i+1,j+w), + color(i+1,j)}); + index[k][i]=m; } + ucyclic((angle2-angle1) % 360 == 0); + vcyclic(cyclic(g)); } } @@ -784,6 +873,9 @@ surface operator * (transform3 t, surface s) S.s=new patch[s.s.length]; for(int i=0; i < s.s.length; ++i) S.s[i]=t*s.s[i]; + S.index=copy(s.index); + S.vcyclic=(bool) s.vcyclic; + return S; } @@ -894,10 +986,11 @@ triple[][] subpatchend(triple[][] P, real u, real v) {c4[2],c5[2],c6[2],c7[2]}}; } -patch subpatch(patch s, real ua, real va, real ub, real vb) +patch subpatch(patch s, pair a, pair b) { - assert(ua >= 0 && va >= 0 && ub <= 1 && vb <= 1 && ua < ub && va < vb); - return patch(subpatchbegin(subpatchend(s.P,ub,vb),ua/ub,va/vb), + assert(a.x >= 0 && a.y >= 0 && b.x <= 1 && b.y <= 1 && + a.x < b.x && a.y < b.y); + return patch(subpatchbegin(subpatchend(s.P,b.x,b.y),a.x/b.x,a.y/b.y), s.straight,s.planar); } @@ -958,8 +1051,17 @@ void draw3D(frame f, patch s, material m, light light=currentlight) bool lighton=light.on(); if(!lighton && !invisible((pen) m)) m=emissive(m); + real PRCshininess; + if(prc()) { + // Empirical translation table from Phong-Blinn to PRC shininess model: + static real[] x={0.015,0.025,0.05,0.07,0.1,0.14,0.23,0.5,0.65,0.75,0.85, + 0.875,0.9,1}; + static real[] y={0.05,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.55,0.6,0.7,0.8,0.9,1}; + static realfunction s=fspline(x,y,monotonic); + PRCshininess=s(m.shininess); + } real granularity=m.granularity >= 0 ? m.granularity : defaultgranularity; - draw(f,s.P,s.straight,m.p,m.opacity,m.shininess,granularity, + draw(f,s.P,s.straight,m.p,m.opacity,m.shininess,PRCshininess,granularity, s.planar ? s.normal(0.5,0.5) : O,lighton,s.colors); } @@ -972,7 +1074,7 @@ void tensorshade(transform t=identity(), frame f, patch s, } restricted pen[] nullpens={nullpen}; -nullpens.cyclic(true); +nullpens.cyclic=true; void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1, material[] surfacepen, pen[] meshpen=nullpens, @@ -1011,7 +1113,7 @@ void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1, depth=sort(depth); - light.T=shiftless(P.modelview()); + light.T=shiftless(P.T.modelview); // Draw from farthest to nearest while(depth.length > 0) { @@ -1033,8 +1135,8 @@ void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1, { material[] surfacepen={surfacepen}; pen[] meshpen={meshpen}; - surfacepen.cyclic(true); - meshpen.cyclic(true); + surfacepen.cyclic=true; + meshpen.cyclic=true; draw(t,f,s,nu,nv,surfacepen,meshpen,light,meshlight,P); } @@ -1082,8 +1184,8 @@ void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1, { material[] surfacepen={surfacepen}; pen[] meshpen={meshpen}; - surfacepen.cyclic(true); - meshpen.cyclic(true); + surfacepen.cyclic=true; + meshpen.cyclic=true; draw(pic,s,nu,nv,surfacepen,meshpen,light,meshlight); } @@ -1092,18 +1194,26 @@ void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1, light light=currentlight, light meshlight=light) { pen[] meshpen={meshpen}; - meshpen.cyclic(true); + meshpen.cyclic=true; draw(pic,s,nu,nv,surfacepen,meshpen,light,meshlight); } -surface extrude(path p, triple axis=Z) +surface extrude(path3 p, path3 q) { static patch[] allocate; - path3 G=path3(p); - path3 G2=shift(axis)*G; return surface(...sequence(new patch(int i) { - return patch(subpath(G,i,i+1)--subpath(G2,i+1,i)--cycle); - },length(G))); + return patch(subpath(p,i,i+1)--subpath(q,i+1,i)--cycle); + },length(p))); +} + +surface extrude(path3 p, triple axis=Z) +{ + return extrude(p,shift(axis)*p); +} + +surface extrude(path p, triple axis=Z) +{ + return extrude(path3(p),axis); } surface extrude(explicit path[] p, triple axis=Z) @@ -1179,11 +1289,11 @@ void label(frame f, Label L, triple position, align align=NoAlign, } else { if(L.filltype == NoFill) fill(f,path(L,project(position,P.t),P), - light.color(L.T3*Z,L.p,shiftless(P.modelview()))); + color(L.T3*Z,L.p,light,shiftless(P.T.modelview))); else { frame d; fill(d,path(L,project(position,P.t),P), - light.color(L.T3*Z,L.p,shiftless(P.modelview()))); + color(L.T3*Z,L.p,light,shiftless(P.T.modelview))); add(f,d,L.filltype); } } @@ -1215,11 +1325,11 @@ void label(picture pic=currentpicture, Label L, triple position, if(pic != null) { if(L.filltype == NoFill) fill(project(v,P.t),pic,path(L,P), - light.color(L.T3*Z,L.p,shiftless(P.modelview()))); + color(L.T3*Z,L.p,light,shiftless(P.T.modelview))); else { picture d; fill(project(v,P.t),d,path(L,P), - light.color(L.T3*Z,L.p,shiftless(P.modelview()))); + color(L.T3*Z,L.p,light,shiftless(P.T.modelview))); add(pic,d,L.filltype); } } @@ -1271,30 +1381,74 @@ surface extrude(Label L, triple axis=Z) restricted surface nullsurface; +surface labelsurface(Label L, surface s, real uoffset, real voffset, + real height=0, bool bottom=false, bool top=true) +{ + int nu=s.index.length; + if(nu == 0) return nullsurface; + int nv=s.index[0].length; + if(nv == 0) return nullsurface; + + path[] g=texpath(L); + pair m=min(g); + pair M=max(g); + pair lambda=inverse(L.T*scale(nu-epsilon,nv-epsilon))*(M-m); + lambda=(abs(lambda.x),abs(lambda.y)); + path[] G=bezulate(g); + + path3 transpath(path p, real height) { + return path3(unstraighten(p),new triple(pair z) { + real u=uoffset+(z.x-m.x)/lambda.x; + real v=voffset+(z.y-m.y)/lambda.y; + if(((u < 0 || u >= nu) && !s.ucyclic()) || + ((v < 0 || v >= nv) && !s.vcyclic())) + abort("cannot fit string to surface"); + return s.point(u,v)+height*unit(s.normal(u,v)); + }); + } + + surface s; + for(path p : G) { + for(path g : regularize(p)) { + path3 b; + bool extrude=height > 0; + if(bottom || extrude) + b=transpath(g,0); + if(bottom) s.s.push(patch(b)); + if(top || extrude) { + path3 h=transpath(g,height); + if(top) s.s.push(patch(h)); + if(extrude) s.append(extrude(b,h)); + } + } + } + return s; +} + private real a=4/3*(sqrt(2)-1); -private transform3 t=rotate(90,O,Z); -private transform3 t2=t*t; -private transform3 t3=t2*t; +private transform3 t1=rotate(90,O,Z); +private transform3 t2=t1*t1; +private transform3 t3=t2*t1; private transform3 i=xscale3(-1)*zscale3(-1); -restricted patch octant1=patch(X{Z}..{-X}Z..Z{Y}..{-Z}Y{X}..{-Y}cycle, - new triple[] {(1,a,a),(a,a^2,1),(a^2,a,1), - (a,1,a)}); +restricted patch octant1=patch(X{Y}..{-X}Y{Z}..{-Y}Z..Z{X}..{-Z}cycle, + new triple[] {(1,a,a),(a,1,a),(a^2,a,1), + (a,a^2,1)}); -restricted surface unithemisphere=surface(octant1,t*octant1,t2*octant1, +restricted surface unithemisphere=surface(octant1,t1*octant1,t2*octant1, t3*octant1); -restricted surface unitsphere=surface(octant1,t*octant1,t2*octant1,t3*octant1, - i*octant1,i*t*octant1,i*t2*octant1, +restricted surface unitsphere=surface(octant1,t1*octant1,t2*octant1,t3*octant1, + i*octant1,i*t1*octant1,i*t2*octant1, i*t3*octant1); restricted patch unitfrustum(real t1, real t2) { real s1=interp(t1,t2,1/3); real s2=interp(t1,t2,2/3); - return patch(interp(Z,X,t2)--interp(Z,X,t1){Y}..{-X}interp(Z,Y,t1)-- - interp(Z,Y,t2){X}..{-Y}cycle, - new triple[] {(s2,s2*a,1-s2),(s1,s1*a,1-s1),(s1*a,s1,1-s1), - (s2*a,s2,1-s2)}); + return patch(interp(Z,X,t2){Y}..{-X}interp(Z,Y,t2)--interp(Z,Y,t1){X}..{-Y} + interp(Z,X,t1)--cycle, + new triple[] {(s2,s2*a,1-s2),(s2*a,s2,1-s2),(s1*a,s1,1-s1), + (s1,s1*a,1-s1)}); } // Return a unitcone constructed from n frusta (the final one being degenerate) @@ -1306,7 +1460,7 @@ surface unitcone(int n=6) for(int i=0; i < n; ++i) { patch s=unitfrustum(i < n-1 ? r^(i+1) : 0,r^i); unitcone.s[i]=s; - unitcone.s[n+i]=t*s; + unitcone.s[n+i]=t1*s; unitcone.s[2n+i]=t2*s; unitcone.s[3n+i]=t3*s; } @@ -1316,9 +1470,9 @@ surface unitcone(int n=6) restricted surface unitcone=unitcone(); restricted surface unitsolidcone=surface(patch(unitcircle3)...unitcone.s); -private patch unitcylinder1=patch(X--X+Z{Y}..{-X}Y+Z--Y{X}..{-Y}cycle); +private patch unitcylinder1=patch(X{Y}..{-X}Y--Y+Z{X}..{-Y}X+Z--cycle); -restricted surface unitcylinder=surface(unitcylinder1,t*unitcylinder1, +restricted surface unitcylinder=surface(unitcylinder1,t1*unitcylinder1, t2*unitcylinder1,t3*unitcylinder1); private patch unitplane=patch(new triple[] {O,X,X+Y,Y}); diff --git a/Master/texmf/asymptote/three_tube.asy b/Master/texmf/asymptote/three_tube.asy new file mode 100644 index 00000000000..c79ae3f619f --- /dev/null +++ b/Master/texmf/asymptote/three_tube.asy @@ -0,0 +1,373 @@ +void render(path3 s, real granularity=linegranularity, void f(path3, real)) +{ + 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 +{ + triple p,r,t,s; + void operator init(triple p, triple r, triple t) + { + this.p=p; + this.r=r; + this.t=t; + s=cross(t,r); + } +} + +// Rotation minimizing frame +// http://www.cs.hku.hk/research/techreps/document/TR-2007-07.pdf +rmf[] rmf(path3 g, real[] t) +{ + rmf[] R=new rmf[t.length]; + triple d=dir(g,0); + R[0]=rmf(point(g,0),perp(d),d); + for(int i=1; i < t.length; ++i) { + rmf Ri=R[i-1]; + real t=t[i]; + triple p=point(g,t); + triple v1=p-Ri.p; + if(v1 != O) { + triple r=Ri.r; + triple u1=unit(v1); + triple ti=Ri.t; + triple tp=ti-2*dot(u1,ti)*u1; + ti=dir(g,t); + triple rp=r-2*dot(u1,r)*u1; + triple u2=unit(ti-tp); + rp=rp-2*dot(u2,rp)*u2; + R[i]=rmf(p,unit(rp),unit(ti)); + } else + R[i]=R[i-1]; + } + return R; +} + +surface bispline(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]; + for(int j=0; j < m; ++j) + if(condi[j]) ++count; + } + } + + surface s=surface(count); + s.index=new int[n][m]; + int k=-1; + for(int i=0; i < n; ++i) { + bool[] condi=all ? null : cond[i]; + real xi=x[i]; + real[] zi=z[i]; + real[] zp=z[i+1]; + real[] ri=r[i]; + real[] rp=r[i+1]; + real[] pi=p[i]; + real[] pp=p[i+1]; + real[] qi=q[i]; + real[] qp=q[i+1]; + real xp=x[i+1]; + real hx=(xp-xi)/3; + int[] indexi=s.index[i]; + for(int j=0; j < m; ++j) { + real yj=y[j]; + real yp=y[j+1]; + if(all || condi[j]) { + triple[][] P=array(4,array(4,O)); + real hy=(yp-yj)/3; + real hxy=hx*hy; + // x and y directions + for(int k=0; k < 4; ++k) { + P[0][k] += xi*X; + P[k][0] += yj*Y; + P[1][k] += (xp+2*xi)/3*X; + P[k][1] += (yp+2*yj)/3*Y; + P[2][k] += (2*xp+xi)/3*X; + P[k][2] += (2*yp+yj)/3*Y; + P[3][k] += xp*X; + P[k][3] += yp*Y; + } + // z: value + P[0][0] += zi[j]*Z; + P[3][0] += zp[j]*Z; + P[0][3] += zi[j+1]*Z; + P[3][3] += zp[j+1]*Z; + // z: first derivative + P[1][0] += (P[0][0].z+hx*pi[j])*Z; + P[1][3] += (P[0][3].z+hx*pi[j+1])*Z; + P[2][0] += (P[3][0].z-hx*pp[j])*Z; + P[2][3] += (P[3][3].z-hx*pp[j+1])*Z; + P[0][1] += (P[0][0].z+hy*qi[j])*Z; + P[3][1] += (P[3][0].z+hy*qp[j])*Z; + P[0][2] += (P[0][3].z-hy*qi[j+1])*Z; + P[3][2] += (P[3][3].z-hy*qp[j+1])*Z; + // z: second derivative + P[1][1] += (P[0][1].z+P[1][0].z-P[0][0].z+hxy*ri[j])*Z; + P[1][2] += (P[0][2].z+P[1][3].z-P[0][3].z-hxy*ri[j+1])*Z; + P[2][1] += (P[2][0].z+P[3][1].z-P[3][0].z-hxy*rp[j])*Z; + P[2][2] += (P[2][3].z+P[3][2].z-P[3][3].z+hxy*rp[j+1])*Z; + s.s[++k]=patch(P); + indexi[j]=k; + } + } + } + + return s; +} + +// return the surface described by a real matrix f, interpolated with +// xsplinetype and ysplinetype. +surface surface(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]); + surface s=bispline(f,p,q,r,x,y,cond); + if(xsplinetype == periodic) s.ucyclic(true); + if(ysplinetype == periodic) s.vcyclic(true); + return s; +} + +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]; + real epsilon=sqrtEpsilon*norm(a); + for(int j=0; j < m; ++j) + if(abs(a0[j]-a1[j]) > epsilon) return false; + return true; +} +bool vperiodic(real[][] a) { + int n=a.length; + if(n == 0) return false; + int m=a[0].length-1; + real epsilon=sqrtEpsilon*norm(a); + for(int i=0; i < n; ++i) + if(abs(a[i][0]-a[i][m]) > epsilon) return false; + return true; +} + +// 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) +{ + 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]); + triple f=(all || (activei[j]=cond(z))) ? f(z) : O; + 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"); + + surface sx=surface(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active); + surface sy=surface(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active); + surface sz=surface(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active); + + surface s=surface(sx.s.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.s.length; ++k) { + triple[][] Q=new triple[4][]; + for(int i=0; i < 4 ; ++i) + Q[i]=sequence(new triple(int j) { + return (sx.s[k].P[i][j].z,sy.s[k].P[i][j].z,sz.s[k].P[i][j].z); + },4); + s.s[k]=patch(Q); + } + + 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); + + 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)); +} + +struct tube +{ + surface s; + path3 center; + + void operator init(path3 p, real width, int sectors=4, + real granularity=linegranularity) { + sectors += sectors % 2; // Must be even. + int h=quotient(sectors,2); + real r=0.5*width; + + void generate(path3 p) { + int n=length(p); + if(piecewisestraight(p)) { + for(int i=0; i < n; ++i) { + triple v=point(p,i); + triple u=point(p,i+1)-v; + s.append(shift(v)*align(unit(u))*scale(r,r,abs(u))*unitcylinder); + } + center=center&p; + } else { + real[] T; + path3 G; + for(int i=0; i < n; ++i) + render(subpath(p,i,i+1),granularity, + new void(path3 g, real s) { + G=G&g; + T.push(i+s); + }); + T.push(n); + T.cyclic=cyclic(p); + rmf[] rmf=rmf(p,T); + triple f(pair t) { + rmf R=rmf[round(t.x)]; + return point(G,t.x)+r*(R.r*cos(t.y)-R.s*sin(t.y)); + } + + real[] v=uniform(0,2pi,sectors); + static splinetype[] Monotonic={monotonic,monotonic,monotonic}; + static splinetype[] Periodic={periodic,periodic,periodic}; + if(T.length > 0) { + surface S=surface(f,sequence(T.length),v,Monotonic,Periodic); + 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]; + pre[0]=point[0]=0.5*(S.s[index[0]].P[0][0]+S.s[index[h]].P[0][0]); + for(int i=0; i < n; ++i) { + index=S.index[i]; + triple [][] P=S.s[index[0]].P; + triple [][] Q=S.s[index[h]].P; + post[i]=0.5*(P[1][0]+Q[1][0]); + pre[i+1]=0.5*(P[2][0]+Q[2][0]); + point[i+1]=0.5*(P[3][0]+Q[3][0]); + } + index=S.index[n-1]; + post[n]=0.5*(S.s[index[0]].P[3][0]+S.s[index[h]].P[3][0]); + center=center&path3(pre,point,post,array(n+1,false),T.cyclic); + } + } + } + + transform3 t=scale3(r); + bool cyclic=cyclic(p); + int begin=0; + int n=length(p); + for(int i=cyclic ? 0 : 1; i < n; ++i) + if(abs(dir(p,i,1)-dir(p,i,-1)) > sqrtEpsilon) { + generate(subpath(p,begin,i)); + s.append(shift(point(p,i))*t*align(dir(p,i,-1))*unithemisphere); + begin=i; + } + generate(subpath(p,begin,n)); + } +} diff --git a/Master/texmf/asymptote/tree.asy b/Master/texmf/asymptote/tree.asy index 0f3559112c3..1e603ec4b6b 100644 --- a/Master/texmf/asymptote/tree.asy +++ b/Master/texmf/asymptote/tree.asy @@ -7,8 +7,8 @@ struct tree { - tree left = null; - tree right = null; + tree left; + tree right; int key = 0; int value = 0; } diff --git a/Master/texmf/asymptote/trembling.asy b/Master/texmf/asymptote/trembling.asy index fc50ce9b7f5..b01d572d6e2 100644 --- a/Master/texmf/asymptote/trembling.asy +++ b/Master/texmf/asymptote/trembling.asy @@ -84,7 +84,7 @@ real atime(pair m, path g, real fuzz=trembleFuzz()) // if(t.length <= 0) { // degenerated case // r=initr; // T *= shift(shx,0); - // write("warning: atime needs numerical adjustment."); + // warning("atime","atime needs numerical adjustment.",position=true); // } // } while(t.length <= 0); if(t.length > 0) ot=t[1]; diff --git a/Master/texmf/asymptote/tube.asy b/Master/texmf/asymptote/tube.asy index 1f8ef950251..844a39c92f0 100644 --- a/Master/texmf/asymptote/tube.asy +++ b/Master/texmf/asymptote/tube.asy @@ -58,20 +58,7 @@ real[] sample(path3 g, real r, real relstep=0) return t; } -struct Rmf -{ - triple p,r,t; // s=cross(t,r); - real reltime; - void operator init(triple p, triple r, triple t, real reltime) - { - this.p=p; - this.r=r; - this.t=t; - this.reltime=reltime; - } -} - -real degrees(Rmf a, Rmf b) +real degrees(rmf a, rmf b) { real d=degrees(acos1(dot(a.r,b.r))); real dt=dot(cross(a.r,b.r),a.t); @@ -79,37 +66,6 @@ real degrees(Rmf a, Rmf b) return d%360; } -private Rmf[] rmf(path3 g, Rmf U0=Rmf(O,O,O,0), real[] t) -{ - if(U0.t == O) { - triple d=dir(g,0); - U0=Rmf(point(g,0),perp(d),d,0); - } - int l=length(g); - Rmf[] R={U0}; - triple rp,v1,v2,tp,ti,p; - real c; - - for(int i=0; i < t.length-1; ++i) { - p=point(g,t[i+1]); - v1=p-R[i].p; - c=dot(v1,v1); - if(c != 0) { - rp=R[i].r-2*dot(v1,R[i].r)*v1/c; - ti=R[i].t; - tp=ti-2*dot(v1,ti)*v1/c; - ti=dir(g,t[i+1]); - v2=ti-tp; - rp=rp-2*dot(v2,rp)*v2/dot(v2,v2); - R.push(Rmf(p,unit(rp),unit(ti),t[i+1]/l)); - } else { - write("Warning: path3 has duplicated point in Rmf."); - R.push(R[R.length-1]); - } - } - return R; -} - restricted int coloredNodes=1; restricted int coloredSegments=2; @@ -158,8 +114,7 @@ coloredpath operator cast(guide p) return coloredpath(p); } -private surface surface(Rmf[] R, coloredpath cp,transform T(real)= - new transform(real t) {return identity();}, +private surface surface(rmf[] R, real[] t, coloredpath cp, transform T(real), bool cyclic) { path g=cp.p; @@ -169,22 +124,22 @@ private surface surface(Rmf[] R, coloredpath cp,transform T(real)= planar[i]=straight(g,i); surface s; - path3 sec=path3(T(R[0].reltime)*g); + path3 sec=path3(T(t[0]/l)*g); real adjust=0; if(cyclic) adjust=-degrees(R[0],R[R.length-1])/(R.length-1); - path3 sec1=shift(R[0].p)*transform3(R[0].r,cross(R[0].t,R[0].r),R[0].t)*sec, + path3 sec1=shift(R[0].p)*transform3(R[0].r,R[0].s,R[0].t)*sec, sec2; for(int i=1; i < R.length; ++i) { - sec=path3(T(R[i].reltime)*g); + sec=path3(T(t[i]/l)*g); sec2=shift(R[i].p)*transform3(R[i].r,cross(R[i].t,R[i].r),R[i].t)* rotate(i*adjust,Z)*sec; for(int j=0; j < l; ++j) { surface st=surface(subpath(sec1,j,j+1)--subpath(sec2,j+1,j)--cycle, planar=planar[j]); if(cp.usepens) { - pen[] tp1=cp.pens(R[i-1].reltime), tp2=cp.pens(R[i].reltime); - tp1.cyclic(true); tp2.cyclic(true); + pen[] tp1=cp.pens(t[i-1]/l), tp2=cp.pens(t[i]/l); + tp1.cyclic=true; tp2.cyclic=true; if(cp.colortype == coloredSegments) { st.colors(new pen[][] {{tp1[j],tp1[j],tp2[j],tp2[j]}}); } else { @@ -206,6 +161,6 @@ surface tube(path3 g, coloredpath section, real[] t=sample(g,max(M.x-m.x,M.y-m.y)/max(realEpsilon,abs(corner)), min(abs(relstep),1)); bool cyclic=cyclic(g); - t.cyclic(cyclic); - return surface(rmf(g,t),section,T,cyclic); + t.cyclic=cyclic; + return surface(rmf(g,t),t,section,T,cyclic); } diff --git a/Master/texmf/asymptote/version.asy b/Master/texmf/asymptote/version.asy index 713dc924f1c..d6c18f1b6de 100644 --- a/Master/texmf/asymptote/version.asy +++ b/Master/texmf/asymptote/version.asy @@ -1 +1 @@ -string VERSION="1.82"; +string VERSION="1.83"; |