diff options
author | Karl Berry <karl@freefriends.org> | 2010-05-16 18:38:18 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2010-05-16 18:38:18 +0000 |
commit | 9b042a93e8f489c694af2e2e121e5aedd7671d8a (patch) | |
tree | 4b4e6b79f46eb1faf6b5ca4b5318ae907daf38bc /Master/texmf/asymptote/contour.asy | |
parent | 165b3785f8aa3ff267097b479870e9a39ba5c14b (diff) |
asymptote 1.94
git-svn-id: svn://tug.org/texlive/trunk@18288 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf/asymptote/contour.asy')
-rw-r--r-- | Master/texmf/asymptote/contour.asy | 1705 |
1 files changed, 1257 insertions, 448 deletions
diff --git a/Master/texmf/asymptote/contour.asy b/Master/texmf/asymptote/contour.asy index fbb4cd1c71a..27a79439e2d 100644 --- a/Master/texmf/asymptote/contour.asy +++ b/Master/texmf/asymptote/contour.asy @@ -1,307 +1,922 @@ -// Contour routines written by Radoslav Marinov and John Bowman. - +/* + Contour routines written by Radoslav Marinov, John Bowman, and Chris Savage. + + [2009/10/15: C Savage] generate oriented contours + [2009/10/19: C Savage] use boxes instead of triangles +*/ + +/* + Contours lines/guides are oriented throughout. By convention, + for a single contour, higher values are to the left and/or lower + values are to the right along the direction of the lines/guide. +*/ + import graph_settings; -real eps=10000*realEpsilon; - -// 1 -// 6 +-------------------+ 5 -// | \ / | -// | \ / | -// | \ / | -// | \ / | -// 2 | X | 0 -// | / \ | -// | / \ | -// | / \ | -// | / \ | -// 7 +-------------------+ 4 or 8 -// 3 +private real eps=sqrtEpsilon; -private struct segment -{ - bool active; - pair a,b; // Endpoints; a is always an edge point if one exists. - int c; // Contour value. - int edge; // -1: interior, 0 to 3: edge, - // 4-8: single-vertex edge, 9: double-vertex edge. +/* + GRID CONTOURS + + Contours on a grid of points are determined as follows: + for each grid square, the function is approximated as the unique + paraboloid passing through the function values at the four + corners. The intersection of a paraboloid with the f(x,y)=c + plane is a line or hyperbola. + + Grid data structures: + + boxcontour: + Describes a particular contour segment in a grid square. + + boxdata: + Describes contours in a grid square (holds boxcontours). + + segment: + Describes a contour line. Usually a closed (interior) contour, + a line that terminates on the border, or a border segment used + to enclose a region. + + Segment: + Describes a contour line. + + Main grid routines: + + setcontour: + Determines the contours in a grid square. + + contouredges: + Determines the contour segments over a grid of function values. + + connect: + Converts contours into guides + +*/ + +private typedef int boxtype; +private boxtype exterior=-1; +private boxtype edge = 0; +private boxtype interior=+1; + +private typedef int contourshape; +private contourshape line =1; +private contourshape hyperbola=2; + +// Describe position by grid square and position in square +private struct gridpoint { + int i,j; + pair z; + void operator init(int i, int j, pair z) { + this.i=i; + this.j=j; + this.z=z; + } + void operator init(gridpoint gp) { + this.i=gp.i; + this.j=gp.j; + this.z=gp.z; + } } -// Case 1: line passes through two vertices of a triangle -private segment case1(pair p0, pair p1, int edge) +private bool same(gridpoint gp1, gridpoint gp2) { - // Will cause a duplicate guide; luckily case1 is rare - segment rtrn; - rtrn.active=true; - rtrn.a=p0; - rtrn.b=p1; - rtrn.edge=edge; - return rtrn; + return abs(gp2.z-gp1.z+(gp2.i-gp1.i,gp2.j-gp1.j)) < eps; } -// Case 2: line passes through a vertex and a side of a triangle -// (the first vertex passed and the side between the other two) -private segment case2(pair p0, pair p1, pair p2, - real v0, real v1, real v2, int edge) -{ - segment rtrn; - pair val=interp(p1,p2,abs(v1/(v2-v1))); - rtrn.active=true; - if(edge < 4) { - rtrn.a=val; - rtrn.b=p0; - } else { - rtrn.a=p0; - rtrn.b=val; + +// Describe contour in unit square(scaling to be done later). +private struct boxcontour { + bool active; + contourshape type; // Shape of contour segment(line or hyperbola) + pair a,b; // Start/end point of contour segment. + // Higher values to left along a--b. + real x0,y0,m; // For hyperbola: (x-x0)*(y-y0)=m + int signx,signy; // Sign of x-x0&y-y0 for hyperbola piece; + // identifies which direction it opens + int i,j; // Indices of lower left corner in position or + // data array. + int index; // Contour index + + void operator init(contourshape type, pair a, pair b, + real x0, real y0, real m, int signx, int signy, + int i, int j, int index) { + this.active=true; + this.type=type; + this.a=a; + this.b=b; + + this.x0=x0; + this.y0=y0; + this.m=m; + this.signx=signx; + this.signy=signy; + + this.i=i; + this.j=j; + this.index=index; + } + // Generate list of points along the line/hyperbola segment + // representing the contour in the box + gridpoint[] points(int subsample=1, bool first=true, bool last=true) { + gridpoint[] gp; + if(first) + gp.push(gridpoint(i,j,a)); + if(subsample > 0) { + // Linear case + if(type == line) { + for(int k=1; k <= subsample; ++k) { + pair z=interp(a,b,k/(subsample+1)); + gp.push(gridpoint(i,j,z)); + } + } else if(type == hyperbola) { + // Special hyperbolic case of m=0 + // The contours here are infinite lines at x=x0 and y=y0, + // but handedness always connects a semi-infinite + // horizontal segment with a semi-infinite vertical segment + // connected at (x0,y0). + // If (x0,y0) is outside the unit box, there is only one + // line segment to include; otherwise, there are both + // a horizontal and a vertical line segment to include. + if(m == 0) { + // Single line + if(a.x == b.x || a.y == b.y) { + for(int k=1; k <= subsample; ++k) { + pair z=interp(a,b,k/(subsample+1)); + gp.push(gridpoint(i,j,z)); + } + // Two lines(may get one extra point here) + } else { + int nsub=quotient(subsample,2); + pair mid=(x0,y0); + for(int k=1; k <= nsub; ++k) { + pair z=interp(a,mid,k/(nsub+1)); + gp.push(gridpoint(i,j,z)); + } + gp.push(gridpoint(i,j,mid)); + for(int k=1; k <= nsub; ++k) { + pair z=interp(mid,b,k/(nsub+1)); + gp.push(gridpoint(i,j,z)); + } + } + // General hyperbolic case (m != 0). + // Parametric equations(m > 0): + // x(t)=x0 +/- sqrt(m)*exp(t) + // y(t)=y0 +/- sqrt(m)*exp(-t) + // Parametric equations (m < 0): + // x(t)=x0 +/- sqrt(-m)*exp(t) + // y(t)=y0 -/+ sqrt(-m)*exp(-t) + // Points will be taken equally spaced in parameter t. + } else { + real sqrtm=sqrt(abs(m)); + real ta=log(signx*(a.x-x0)/sqrtm); + real tb=log(signx*(b.x-x0)/sqrtm); + real[] t=uniform(ta,tb,subsample+1); + for(int k=1; k <= subsample; ++k) { + pair z=(x0+signx*sqrtm*exp(t[k]), + y0+signy*sqrtm*exp(-t[k])); + gp.push(gridpoint(i,j,z)); + } + } + } + } + if(last) + gp.push(gridpoint(i,j,b)); + + return gp; } - rtrn.edge=edge; - return rtrn; } -// Case 3: line passes through two sides of a triangle -// (through the sides formed by the first & second, and second & third -// vertices) -private segment case3(pair p0, pair p1, pair p2, - real v0, real v1, real v2, int edge=-1) -{ - segment rtrn; - rtrn.active=true; - rtrn.a=interp(p1,p0,abs(v1/(v0-v1))); - rtrn.b=interp(p1,p2,abs(v1/(v2-v1))); - rtrn.edge=edge; - return rtrn; +// Hold data for a single grid square +private struct boxdata { + boxtype type; // Does box contain a contour segment (edge of + // contour region) or is it entirely interior/ + // exterior to contour region ? + real min,max; // Smallest/largest corner value + real max2; // Second-largest corner value + boxcontour[] data; // Stores actual contour segment data + + int count() {return data.length;} + void operator init(real f00, real f10, real f01, real f11) { + real[] X={f00,f10,f01,f11}; + min=min(X); + max=max(X); + X.delete(find(X == max)); + max2=max(X); + } + void settype(real c) { + // Interior case(f >= c) + if(min > c-eps) { + type=interior; + // Exterior case(f < c) + } else if(max < c-eps) { + type=exterior; + // Special case: only one corner at f=c, f < c elsewhere + //(no segment in this case) + } else if((max < c+eps) && (max2 < c-eps)) { + type=exterior; + // Edge of contour passes through box + } else { + type=edge; + } + } } -// Check if a line passes through a triangle, and draw the required line. -private segment checktriangle(pair p0, pair p1, pair p2, - real v0, real v1, real v2, int edge=-1) -{ - // default null return - static segment dflt; - real eps=eps*max(abs(v0),abs(v1),abs(v2)); +/* + Determine contours within a unit square box. - if(v0 < -eps) { - if(v1 < -eps) { - if(v2 < -eps) return dflt; // nothing to do - else if(v2 <= eps) return dflt; // nothing to do - else return case3(p0,p2,p1,v0,v2,v1); - } else if(v1 <= eps) { - if(v2 < -eps) return dflt; // nothing to do - else if(v2 <= eps) return case1(p1,p2,5+edge); - else return case2(p1,p0,p2,v1,v0,v2,5+edge); - } else { - if(v2 < -eps) return case3(p0,p1,p2,v0,v1,v2,edge); - else if(v2 <= eps) - return case2(p2,p0,p1,v2,v0,v1,edge); - else return case3(p1,p0,p2,v1,v0,v2,edge); - } - } else if(v0 <= eps) { - if(v1 < -eps) { - if(v2 < -eps) return dflt; // nothing to do - else if(v2 <= eps) return case1(p0,p2,4+edge); - else return case2(p0,p1,p2,v0,v1,v2,4+edge); - } else if(v1 <= eps) { - if(v2 < -eps) return case1(p0,p1,9); - else if(v2 <= eps) return dflt; // use finer partitioning. - else return case1(p0,p1,9); + Here, we approximate the function on the unit square to be a quadric + surface passing through the specified values at the four corners: + f(x,y)=(1-x)(1-y) f00+x(1-y) f10+(1-x)y f01+xy f11 + =a0+ax x+ay y+axy xy + where f00, f10, f01 and f11 are the function values at the four + corners of the unit square 0 < x < 1&0 < y < 1 and: + a0 =f00 + ax =f10-f00 + ay =f01-f00 + axy=f00+f11-f10-f01 + This can also be expressed in paraboloid form as: + f(x,y)=alpha [(x+y-cp)^2-(x-y-cn)^2]+d + where: + alpha=axy/4 + cp =-(ax+ay)/a11 + cn =-(ax-ay)/a11 + d =(a0 axy-ax ay)/axy + In the procedure below, we take f00 - > f00-c etc. for a contour + level c and we search for f=0. + + For this surface, there are two possible contour shapes: + linear: (y-y0)/(x-x0)=m + hyperbolic: (x-x0)*(y-y0)=m + The linear case has a single line. The hyperbolic case may have + zero, one or two segments within the box (there are two sides of + a hyperbola, each of which may or may not pass through the unit + square). A hyperbola with m=0 is a special case that is handled + separately below. + + If c0 is the desired contour level, we effectively find the + contours at c0-epsilon for arbitrarily small epsilon. Flat + regions equal to c0 are considered to be interior to the + contour curves. Regions that lie at the contour level are + considered to be interior to the contour curves. As a result, + contours are only constructed if they are immediately adjacent + to some region interior to the square that falls below the + contour value; in other words, if an edge falls on the contour + value, but a point within the square arbitrarily close to the + edge falls above the contour value, that edge (or applicable + portion) is not included. This requirement gives the following: + *) ensures contours on an edge are unique (do not appear in + an adjacent square with the same orientation) + *) no three line vertices (four line vertices are possible, but + are not usually an issue) + *) all segments can be joined into closed curves or curves that + terminate on the boundary (no unclosed curves terminate in + the interior region of the grid) + + Note the logic below skips cases that have been filtered out + by the boxdata.settype() routine. +*/ +private void setcontour(real f00, real f10, real f01, real f11, real epsf, + boxdata bd, int i, int j, int index) { + // SPECIAL CASE: two diagonal corners at the contour level with + // the other two below does not yield any contours within the + // unit box, but may have been previously misidentified as an + // edge containing region. + if(((f00*f11 == 0) && (f10*f01 > 0)) || ((f01*f10 == 0) && (f00*f11 > 0))) { + bd.type=exterior; + return; + } + + // NOTE: From this point on, we can assume at least one contour + // segment exists in the square. This allows several cases to + // be ignored or simplified below, particularly edge cases. + + // Form used to approximate function on unit square + real F(real x, real y) { + return interp(interp(f00,f10,x),interp(f01,f11,x),y); + } + + // Write contour as a0+ax*x+ay*y +axy*x*y=0 + real a0 =f00; + real ax =f10-f00; + if(abs(ax) < epsf) ax=0; + real ay =f01-f00; + if(abs(ay) < epsf) ay=0; + real axy=f00+f11-f01 -f10; + if(abs(axy) < epsf) axy=0; + + // Linear contour(s) + if(axy == 0) { + pair a,b; + // Horizontal + if(ax == 0) { + if(ay == 0) return; // Contour is at most an isolated point; ignore. + real y0=-a0/ay; + if(abs(y0-1) < eps) y0=1; + if((f00 > 0) || (f01 < 0)) { + a=(1,y0); + b=(0,y0); + } else { + a=(0,y0); + b=(1,y0); + } + // Vertical + } else if(ay == 0) { + real x0=-a0/ax; + if(abs(x0-1) < eps) x0=1; + if((f00 > 0) || (f10 < 0)) { + a=(x0,0); + b=(x0,1); + } else { + a=(x0,1); + b=(x0,0); + } + // Angled line } else { - if(v2 < -eps) return case2(p0,p1,p2,v0,v1,v2,4+edge); - else if(v2 <= eps) return case1(p0,p2,4+edge); - else return dflt; // nothing to do - } - } else { - if(v1 < -eps) { - if(v2 < -eps) return case3(p1,p0,p2,v1,v0,v2,edge); - else if(v2 <= eps) - return case2(p2,p0,p1,v2,v0,v1,edge); - else return case3(p0,p1,p2,v0,v1,v2,edge); - } else if(v1 <= eps) { - if(v2 < -eps) return case2(p1,p0,p2,v1,v0,v2,5+edge); - else if(v2 <= eps) return case1(p1,p2,5+edge); - else return dflt; // nothing to do + real x0=-a0/ax; + if(abs(x0-1) < eps) x0=1; + real y0=-a0/ay; + if(abs(y0-1) < eps) y0=1; + int count=0; + real[] farr={f00,f10,f11,f01}; + farr.cyclic=true; + pair[] corners={(0,0),(1,0),(1,1),(0,1)}; + pair[] sidedir={(1,0),(0,1),(-1,0),(0,-1)}; + + int count=0; + for(int i=0; i < farr.length; ++i) { + // Corner + if(farr[i] == 0) { + ++count; + if(farr[i-1] > 0) { + a=corners[i]; + } else { + b=corners[i]; + } + // Side + } else if(farr[i]*farr[i+1] < 0) { + ++count; + if(farr[i] > 0) { + a=corners[i]-(farr[i]/(farr[i+1]-farr[i]))*sidedir[i]; + } else { + b=corners[i]-(farr[i]/(farr[i+1]-farr[i]))*sidedir[i]; + } + } + } + // Check(if logic is correct above, this will not happen) + if(count != 2) { + abort("Unexpected error in setcontour routine: odd number of" + +" crossings (linear case)"); + } + } + boxcontour bc=boxcontour(line,a,b,0,0,0,1,1,i,j,index); + bd.data.push(bc); + return; + } + + // Hyperbolic contour(s) + // Described in form: (x-x0)*(y-y0)=m + real x0=-ay/axy; + if(abs(x0-1) < eps) x0=1; + real y0=-ax/axy; + if(abs(y0-1) < eps) y0=1; + real m =ay*ax-a0*axy; + m=(abs(m) < eps) ? 0 : m/axy^2; + + // Special case here: straight segments (possibly crossing) + if(m == 0) { + pair a,b; + int signx,signy; + // Assuming at least one corner is below contour level here + if(x0 == 0) { + signx=+1; + if(y0 == 0) { + a=(1,0); + b=(0,1); + signy=+1; + } else if(y0 == 1) { + a=(0,0); + b=(1,1); + signy=-1; + } else if(y0 < 0 || y0 > 1) { + a=(0,0); + b=(0,1); + signy=y0 > 0 ? -1 : +1; + } else { + if(f10 > 0) { + a=(1,y0); + b=(0,1); + signy=+1; + } else { + a=(0,0); + b=(1,y0); + signy=-1; + } + } + boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index); + bd.data.push(bc); + return; + } else if(x0 == 1) { + signx=-1; + if(y0 == 0) { + a=(1,1); + b=(0,0); + signy=+1; + } else if(y0 == 1) { + a=(0,1); + b=(1,0); + signy=-1; + } else if(y0 < 0 || y0 > 1) { + a=(1,1); + b=(1,0); + signy=y0 > 0 ? -1 : +1; + } else { + if(f01 > 0) { + a=(0,y0); + b=(1,0); + signy=-1; + } else { + a=(1,1); + b=(0,y0); + signy=+1; + } + } + boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index); + bd.data.push(bc); + return; + } else if(y0 == 0) { + signy=+1; + if(x0 < 0 || x0 > 1) { + a=(1,0); + b=(0,0); + signx=x0 > 0 ? -1 : +1; + } else { + if(f11 > 0) { + a=(x0,1); + b=(0,0); + signx=-1; + } else { + a=(1,0); + b=(x0,1); + signx=+1; + } + } + boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index); + bd.data.push(bc); + return; + } else if(y0 == 1) { + signy=-1; + if(x0 < 0 || x0 > 1) { + a=(0,1); + b=(1,1); + signx=x0 > 0 ? -1 : +1; + } else { + if(f00 > 0) { + a=(x0,0); + b=(1,1); + signx=+1; + } else { + a=(0,1); + b=(x0,0); + signx=-1; + } + } + boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index); + bd.data.push(bc); + return; + } else if(x0 < 0 || x0 > 1) { + signx=x0 > 0 ? -1 : +1; + if(f00 > 0) { + a=(1,y0); + b=(0,y0); + signy=+1; + } else { + a=(0,y0); + b=(1,y0); + signy=-1; + } + boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index); + bd.data.push(bc); + return; + } else if(y0 < 0 || y0 > 1) { + signy=y0 > 0 ? -1 : +1; + if(f00 > 0) { + a=(x0,0); + b=(x0,1); + signx=+1; + } else { + a=(x0,1); + b=(x0,0); + signx=-1; + } + boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index); + bd.data.push(bc); + return; } else { - if(v2 < -eps) return case3(p0,p2,p1,v0,v2,v1); - else if(v2 <= eps) return dflt; // nothing to do - else return dflt; // nothing to do - } - } -} + if(f10 > 0) { + a=(0,y0); + b=(x0,0); + boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,-1,-1,i,j,index); + bd.data.push(bc); + a=(1,y0); + b=(x0,1); + bc=boxcontour(hyperbola,a,b,x0,y0,m,+1,+1,i,j,index); + bd.data.push(bc); + return; + } else { + a=(x0,0); + b=(1,y0); + boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,+1,-1,i,j,index); + bd.data.push(bc); + a=(x0,1); + b=(0,y0); + bc=boxcontour(hyperbola,a,b,x0,y0,m,-1,+1,i,j,index); + bd.data.push(bc); + return; + } + } + } + + // General hyperbola case + int signc=(F(x0,y0) > 0) ? +1 : -1; + + pair[] points; + + real xB=(y0 == 0) ? infinity : x0-m/y0; + if(abs(xB) < eps) xB=0; + if(xB >= 0 && xB <= 1-eps) points.push((xB,0)); -// Collect connecting path segments. -private void collect(pair[][][] points, real[] c) -{ - // use to reverse an array, omitting the first point - int[] reverseF(int n) {return sequence(new int(int x){return n-1-x;},n-1);} - // use to reverse an array, omitting the last point - int[] reverseL(int n) {return sequence(new int(int x){return n-2-x;},n-1);} + real xT=(y0 == 1) ? infinity : x0+m/(1-y0); + if(abs(xT-1) < eps) xT=1; + if(xT >= eps && xT <= 1) points.push((xT,1)); + + real yL=(x0 == 0) ? infinity : y0-m/x0; + if(abs(yL-1) < eps) yL=1; - for(int cnt=0; cnt < c.length; ++cnt) { - pair[][] gdscnt=points[cnt]; - for(int i=0; i < gdscnt.length; ++i) { - pair[] gig=gdscnt[i]; - int Li=gig.length; - for(int j=i+1; j < gdscnt.length; ++j) { - pair[] gjg=gdscnt[j]; - int Lj=gjg.length; - if(abs(gig[0]-gjg[0]) < eps) { - gdscnt[j]=gjg[reverseF(Lj)]; - gdscnt[j].append(gig); - gdscnt.delete(i); - --i; - break; - } else if(abs(gig[0]-gjg[Lj-1]) < eps) { - gig.delete(0); - gdscnt[j].append(gig); - gdscnt.delete(i); - --i; - break; - } else if(abs(gig[Li-1]-gjg[0]) < eps) { - gjg.delete(0); - gig.append(gjg); - gdscnt[j]=gig; - gdscnt.delete(i); - --i; - break; - } else if(abs(gig[Li-1]-gjg[Lj-1]) < eps) { - gig.append(gjg[reverseL(Lj)]); - gdscnt[j]=gig; - gdscnt.delete(i); - --i; - break; - } + if(yL > eps && yL <= 1) points.push((0,yL)); + + real yR=(x0 == 1) ? infinity : y0+m/(1-x0); + if(abs(yR) < eps) yR=0; + if(yR >= 0 && yR <= 1-eps) points.push((1,yR)); + + // Check (if logic is correct above, this will not happen) + if(!(points.length == 2 || points.length == 4)) { + abort("Unexpected error in setcontour routine: odd number of" + +" crossings (hyperbolic case)"); + } + + // Lower left side + if((x0 > 0) && (y0 > 0) && (f00*signc < 0)) { + pair[] pts0; + for(int i=0; i < points.length; ++i) { + if((points[i].x < x0) && (points[i].y < y0)) { + pts0.push(points[i]); + } + } + if(pts0.length == 2) { + pair a0,b0; + if((f00 > 0) ^(pts0[0].x < pts0[1].x)) { + a0=pts0[0]; + b0=pts0[1]; + } else { + a0=pts0[1]; + b0=pts0[0]; + } + boxcontour bc=boxcontour(hyperbola,a0,b0,x0,y0,m,-1,-1,i,j,index); + bd.data.push(bc); + } + } + + // Lower right side + if((x0 < 1) && (y0 > 0) && (f10*signc < 0)) { + pair[] pts0; + for(int i=0; i < points.length; ++i) { + if((points[i].x > x0) && (points[i].y < y0)) { + pts0.push(points[i]); + } + } + if(pts0.length == 2) { + pair a0,b0; + if((f10 > 0) ^(pts0[0].x < pts0[1].x)) { + a0=pts0[0]; + b0=pts0[1]; + } else { + a0=pts0[1]; + b0=pts0[0]; + } + boxcontour bc=boxcontour(hyperbola,a0,b0,x0,y0,m,+1,-1,i,j,index); + bd.data.push(bc); + } + } + + // Upper right side + if((x0 < 1) && (y0 < 1) && (f11*signc < 0)) { + pair[] pts0; + for(int i=0; i < points.length; ++i) { + if((points[i].x > x0) && (points[i].y > y0)) { + pts0.push(points[i]); + } + } + if(pts0.length == 2) { + pair a0,b0; + if((f11 > 0) ^(pts0[0].x > pts0[1].x)) { + a0=pts0[0]; + b0=pts0[1]; + } else { + a0=pts0[1]; + b0=pts0[0]; } + boxcontour bc=boxcontour(hyperbola,a0,b0,x0,y0,m,+1,+1,i,j,index); + bd.data.push(bc); } } + + // Upper left side + if((x0 > 0) && (y0 < 1) && (f01*signc < 0)) { + pair[] pts0; + for(int i=0; i < points.length; ++i) { + if((points[i].x < x0) && (points[i].y > y0)) { + pts0.push(points[i]); + } + } + if(pts0.length == 2) { + pair a0,b0; + if((f01 > 0) ^(pts0[0].x > pts0[1].x)) { + a0=pts0[0]; + b0=pts0[1]; + } else { + a0=pts0[1]; + b0=pts0[0]; + } + boxcontour bc=boxcontour(hyperbola,a0,b0,x0,y0,m,-1,+1,i,j,index); + bd.data.push(bc); + } + } + return; } -// Join path segments. -private guide[][] connect(pair[][][] points, real[] c, interpolate join) -{ - // set up return value - guide[][] result=new guide[c.length][]; - for(int cnt=0; cnt < c.length; ++cnt) { - pair[][] pointscnt=points[cnt]; - guide[] resultcnt=result[cnt]=new guide[pointscnt.length]; - for(int i=0; i < pointscnt.length; ++i) { - pair[] pts=pointscnt[i]; - guide gd; - if(pts.length > 0) { - if(pts.length > 1 && abs(pts[0]-pts[pts.length-1]) < eps) { - guide[] g=sequence(new guide(int i) { - return pts[i]; - },pts.length-1); - g.push(cycle); - gd=join(...g); - } else - gd=join(...sequence(new guide(int i) { - return pts[i]; - },pts.length)); + +// Checks if end of first contour segment matches the beginning of +// the second. +private bool connected(boxcontour bc1, boxcontour bc2) { + return abs(bc2.a-bc1.b+(bc2.i-bc1.i,bc2.j-bc1.j)) < eps; +} + +// Returns index of first active element in bca that with beginning +// that connects to the end of bc, or -1 if no such element. +private int connectedindex(boxcontour bc, boxcontour[] bca, + bool activeonly=true) { + for(int i=0; i < bca.length; ++i) { + if(!bca[i].active) continue; + if(connected(bc,bca[i])) { + return i; + } + } + return -1; +} + +// Returns index of first active element in bca with end that connects +// to the start of bc, or -1 if no such element. +private int connectedindex(boxcontour[] bca, boxcontour bc, + bool activeonly=true) { + for(int i=0; i < bca.length; ++i) { + if(!bca[i].active) continue; + if(connected(bca[i],bc)) { + return i; + } + } + return -1; +} + + +// Processes indices for grid regions touching the +// end/start (forward=true/false) of the contour segment +private void searchindex(boxcontour bc, bool forward, void f(int i, int j)) { + pair z=forward ? bc.b : bc.a; + + int i=bc.i; + int j=bc.j; + + if(z == (0,0)) f(i-1,j-1); + if(z.y == 0) f(i,j-1); + if(z == (1,0)) f(i+1,j-1); + if(z.x == 1) f(i+1,j); + if(z == (1,1)) f(i+1,j+1); + if(z.y == 1) f(i,j+1); + if(z == (0,1)) f(i-1,j+1); + if(z.x == 0) f(i-1,j); +} + +// Contour segment +private struct segment { + gridpoint[] data; + void operator init() { + } + void operator init(boxcontour bc, int subsample=1) { + bc.active=false; + this.data.append(bc.points(subsample,first=true,last=true)); + } + void operator init(int i, int j, pair z) { + gridpoint gp=gridpoint(i,j,z); + data.push(gp); + } + void operator init(gridpoint[] gp) { + this.data.append(gp); + } + gridpoint start() { + if(data.length == 0) { + return gridpoint(-1,-1,(-infinity,-infinity)); + } + gridpoint gp=data[0]; + return gridpoint(gp.i,gp.j,gp.z); + } + gridpoint end() { + if(data.length == 0) { + return gridpoint(-1,-1,(-infinity,-infinity)); + } + gridpoint gp=data[data.length-1]; + return gridpoint(gp.i,gp.j,gp.z); + } + bool closed() { + return same(this.start(),this.end()); + } + void append(boxcontour bc, int subsample=1) { + bc.active=false; + data.append(bc.points(subsample,first=false,last=true)); + } + void prepend(boxcontour bc, int subsample=1) { + bc.active=false; + data.insert(0 ... bc.points(subsample,first=true,last=false)); + } + void append(int i, int j, pair z) { + gridpoint gp=gridpoint(i,j,z); + data.push(gp); + } + void prepend(int i, int j, pair z) { + gridpoint gp=gridpoint(i,j,z); + data.insert(0,gp); + } + segment copy() { + segment seg=new segment; + seg.data=new gridpoint[data.length]; + for(int i=0; i < data.length; ++i) { + seg.data[i]=gridpoint(data[i].i,data[i].j,data[i].z); + } + return seg; + } + segment reversecopy() { + segment seg=new segment; + seg.data=new gridpoint[data.length]; + for(int i=0; i < data.length; ++i) { + seg.data[data.length-i-1]=gridpoint(data[i].i,data[i].j,data[i].z); + } + return seg; + } +} + +// Container to hold edge and border segments that form one continuous line +private struct Segment { + segment[] edges; + segment[] borders; + void operator init() { + } + void operator init(segment seg) { + edges.push(seg); + } + void operator init(gridpoint[] gp) { + segment seg=segment(gp); + edges.push(seg); + } + gridpoint start() { + if(edges.length == 0) { + if(borders.length > 0) { + return borders[0].start(); } - resultcnt[i]=gd; + return gridpoint(-1,-1,(-infinity,-infinity)); } + return edges[0].start(); + } + gridpoint end() { + if(edges.length == 0 && borders.length == 0) { + return gridpoint(-1,-1,(-infinity,-infinity)); + } + if(edges.length > borders.length) { + return edges[edges.length-1].end(); + } else { + return borders[borders.length-1].end(); + } + } + bool closed() { + return same(this.start(),this.end()); + } + void addedge(segment seg) { + edges.push(seg); + } + void addedge(gridpoint[] gp) { + segment seg=segment(gp); + edges.push(seg); + } + void addborder(segment seg) { + borders.push(seg); + } + void addborder(gridpoint[] gp) { + segment seg=segment(gp); + borders.push(seg); + } + void append(Segment S) { + edges.append(S.edges); + borders.append(S.borders); } - return result; } +private Segment[] Segment(segment[] s) +{ + return sequence(new Segment(int i) {return Segment(s[i]);},s.length); +} -// Return contour guides for a 2D data array. -// z: two-dimensional array of nonoverlapping mesh points +private Segment[][] Segment(segment[][] s) +{ + Segment[][] S=new Segment[s.length][]; + for(int i=0; i < s.length; ++i) + S[i]=Segment(s[i]); + return S; +} + +// Return contour points for a 2D data array. // f: two-dimensional array of corresponding f(z) data values -// midpoint: optional array containing values of f at cell midpoints // c: array of contour values -// join: interpolation operator (e.g. operator -- or operator ..) -guide[][] contour(pair[][] z, real[][] f, - real[][] midpoint=new real[][], real[] c, - interpolate join=operator --) +// subsample: number of points to use in each box in addition to endpoints +segment[][] contouredges(real[][] f, real[] c, int subsample=1) { - int nx=z.length-1; - if(nx == 0) - abort("array z must have length >= 2"); - int ny=z[0].length-1; - if(ny == 0) - abort("array z[0] must have length >= 2"); + int nx=f.length-1; + if(nx <= 0) + abort("array f must have length >= 2"); + int ny=f[0].length-1; + if(ny <= 0) + abort("array f[0] must have length >= 2"); c=sort(c); - bool midpoints=midpoint.length > 0; + boxdata[][] bd=new boxdata[nx][ny]; + + segment[][] result=new segment[c.length][]; - segment segments[][][]=new segment[nx][ny][]; - - // go over region a rectangle at a time for(int i=0; i < nx; ++i) { - pair[] zi=z[i]; - pair[] zp=z[i+1]; + boxdata[] bdi=bd[i]; real[] fi=f[i]; real[] fp=f[i+1]; - real[] midpointi; - if(midpoints) midpointi=midpoint[i]; - segment[][] segmentsi=segments[i]; + for(int j=0; j < ny; ++j) { - segment[] segmentsij=segmentsi[j]; - - // define points - pair bleft=zi[j]; - pair bright=zp[j]; - pair tleft=zi[j+1]; - pair tright=zp[j+1]; - pair middle=0.25*(bleft+bright+tleft+tright); - - real f00=fi[j]; - real f01=fi[j+1]; - real f10=fp[j]; - real f11=fp[j+1]; - real fmm=midpoints ? midpoint[i][j] : 0.25*(f00+f01+f10+f11); - - // optimization: we make sure we don't work with empty rectangles + boxdata bdij=bdi[j]=boxdata(fi[j],fp[j],fi[j+1],fp[j+1]); + int checkcell(int cnt) { real C=c[cnt]; - real vertdat0=f00-C; // bottom-left vertex - real vertdat1=f10-C; // bottom-right vertex - real vertdat2=f01-C; // top-left vertex - real vertdat3=f11-C; // top-right vertex + + real f00=fi[j]; + real f10=fp[j]; + real f01=fi[j+1]; + real f11=fp[j+1]; + + real epsf=eps*max(abs(f00),abs(f10),abs(f01),abs(f11),abs(C)); + + f00=f00-C; + f10=f10-C; + f01=f01-C; + f11=f11-C; + + if(abs(f00) < epsf) f00=0; + if(abs(f10) < epsf) f10=0; + if(abs(f01) < epsf) f01=0; + if(abs(f11) < epsf) f11=0; + - // optimization: we make sure we don't work with empty rectangles int countm=0; int countz=0; int countp=0; - + void check(real vertdat) { - if(vertdat < -eps) ++countm; + if(vertdat < -eps)++countm; else { - if(vertdat <= eps) ++countz; - else ++countp; + if(vertdat <= eps)++countz; + else++countp; } } - check(vertdat0); - check(vertdat1); - check(vertdat2); - check(vertdat3); + check(f00); + check(f10); + check(f01); + check(f11); if(countm == 4) return 1; // nothing to do if(countp == 4) return -1; // nothing to do if((countm == 3 || countp == 3) && countz == 1) return 0; - // go through the triangles - - void addseg(segment seg) { - if(seg.active) { - seg.c=cnt; - segmentsij.push(seg); - } - } - real vertdat4=fmm-C; - addseg(checktriangle(bright,tright,middle, - vertdat1,vertdat3,vertdat4,0)); - addseg(checktriangle(tright,tleft,middle, - vertdat3,vertdat2,vertdat4,1)); - addseg(checktriangle(tleft,bleft,middle, - vertdat2,vertdat0,vertdat4,2)); - addseg(checktriangle(bleft,bright,middle, - vertdat0,vertdat1,vertdat4,3)); + // Calculate individual box contours + bdij.settype(C); + if(bdij.type == edge) + setcontour(f00,f10,f01,f11,epsf,bdij,i,j,cnt); return 0; } - + void process(int l, int u) { if(l >= u) return; int i=quotient(l+u,2); @@ -313,150 +928,171 @@ guide[][] contour(pair[][] z, real[][] f, process(i+1,u); } } - + process(0,c.length); } } - - // set up return value - pair[][][] points=new pair[c.length][][]; - + + // Find contours and follow them for(int i=0; i < nx; ++i) { - segment[][] segmentsi=segments[i]; + boxdata[] bdi=bd[i]; for(int j=0; j < ny; ++j) { - segment[] segmentsij=segmentsi[j]; - for(int k=0; k < segmentsij.length; ++k) { - segment C=segmentsij[k]; - - if(!C.active) continue; - - pair[] g=new pair[] {C.a,C.b}; - segmentsij[k].active=false; - - int forward(int I, int J, bool first=true) { - if(I >= 0 && I < nx && J >= 0 && J < ny) { - segment[] segmentsIJ=segments[I][J]; - for(int l=0; l < segmentsIJ.length; ++l) { - segment D=segmentsIJ[l]; - if(!D.active) continue; - if(abs(D.a-g[g.length-1]) < eps) { - g.push(D.b); - segmentsIJ[l].active=false; - if(D.edge >= 0 && !first) return D.edge; - first=false; - l=-1; - } else if(abs(D.b-g[g.length-1]) < eps) { - g.push(D.a); - segmentsIJ[l].active=false; - if(D.edge >= 0 && !first) return D.edge; - first=false; - l=-1; - } - } - } - return -1; - } - - int backward(int I, int J, bool first=true) { - if(I >= 0 && I < nx && J >= 0 && J < ny) { - segment[] segmentsIJ=segments[I][J]; - for(int l=0; l < segmentsIJ.length; ++l) { - segment D=segmentsIJ[l]; - if(!D.active) continue; - if(abs(D.a-g[0]) < eps) { - g.insert(0,D.b); - segmentsIJ[l].active=false; - if(D.edge >= 0 && !first) return D.edge; - first=false; - l=-1; - } else if(abs(D.b-g[0]) < eps) { - g.insert(0,D.a); - segmentsIJ[l].active=false; - if(D.edge >= 0 && !first) return D.edge; - first=false; - l=-1; + boxdata bd0=bdi[j]; + if(bd0.count() == 0) continue; + for(int k=0; k < bd0.count(); ++k) { + boxcontour bc0=bd0.data[k]; + + if(!bc0.active) continue; + + // Note: boxcontour set inactive when added to segment + segment seg=segment(bc0,subsample); + + // Forward direction + bool foundnext=true; + while(foundnext) { + foundnext=false; + searchindex(bc0,true,new void(int i, int j) { + if((i >= 0) && (i < nx) && (j >= 0) && (j < ny)) { + boxcontour[] data=bd[i][j].data; + int k0=connectedindex(bc0,data); + if(k0 >= 0) { + bc0=data[k0]; + seg.append(bc0,subsample); + foundnext=true; + } } - } - } - return -1; + }); } - - void follow(int f(int, int, bool first=true), int edge) { - int I=i; - int J=j; - while(true) { - static int ix[]={1,0,-1,0}; - static int iy[]={0,1,0,-1}; - if(edge >= 0 && edge < 4) { - I += ix[edge]; - J += iy[edge]; - edge=f(I,J); - } else { - if(edge == -1) break; - if(edge < 9) { - int edge0=(edge-5) % 4; - int edge1=(edge-4) % 4; - int ix0=ix[edge0]; - int iy0=iy[edge0]; - I += ix0; - J += iy0; - // Search all 3 corner cells - if((edge=f(I,J)) == -1) { - I += ix[edge1]; - J += iy[edge1]; - if((edge=f(I,J)) == -1) { - I -= ix0; - J -= iy0; - edge=f(I,J); - } - } - } else { - // Double-vertex edge: search all 8 surrounding cells - void search() { - for(int i=-1; i <= 1; ++i) { - for(int j=-1; j <= 1; ++j) { - if((edge=f(I+i,J+j,false)) >= 0) { - I += i; - J += j; - return; - } - } - } + + // Backward direction + bc0=bd0.data[k]; + bool foundprev=true; + while(foundprev) { + foundprev=false; + searchindex(bc0,false,new void(int i, int j) { + if((i >= 0) && (i < nx) && (j >= 0) && (j < ny)) { + boxcontour[] data=bd[i][j].data; + int k0=connectedindex(data,bc0); + if(k0 >= 0) { + bc0=data[k0]; + seg.prepend(bc0,subsample); + foundprev=true; } - search(); } - } - } + }); } - // Follow contour in cell - int edge=forward(i,j,first=false); - - // Follow contour forward outside of cell - follow(forward,edge); - - // Follow contour backward outside of cell - follow(backward,C.edge); - - points[C.c].push(g); + result[bc0.index].push(seg); } } } + + // Note: every segment here _should_ be cyclic or terminate on the + // boundary + return result; +} + +// Connect contours into guides. +// Same initial/final points indicates a closed path. +// Borders are always joined using--. +private guide connect(Segment S, pair[][] z, interpolate join) +{ + pair loc(gridpoint gp) { + pair offset=z[gp.i][gp.j]; + pair size=z[gp.i+1][gp.j+1]-z[gp.i][gp.j]; + return offset+(size.x*gp.z.x,size.y*gp.z.y); + } + pair[] loc(gridpoint[] gp) { + pair[] result=new pair[gp.length]; + for(int i; i < gp.length; ++i) { + result[i]=loc(gp[i]); + } + return result; + } + + bool closed=S.closed(); + + pair[][] edges=new pair[S.edges.length][]; + for(int i; i < S.edges.length; ++i) { + edges[i]=loc(S.edges[i].data); + } + pair[][] borders=new pair[S.borders.length][]; + for(int i; i < S.borders.length; ++i) { + borders[i]=loc(S.borders[i].data); + } + + if(edges.length == 0 && borders.length == 1) { + guide g=operator--(...borders[0]); + if(closed) g=g--cycle; + return g; + } + + if(edges.length == 1 && borders.length == 0) { + pair[] pts=edges[0]; + if(closed) pts.delete(pts.length-1); + guide g=join(...pts); + if(closed) g=join(g,cycle); + return g; + } + + guide[] ge=new guide[edges.length]; + for(int i=0; i < ge.length; ++i) + ge[i]=join(...edges[i]); + + guide[] gb=new guide[borders.length]; + for(int i=0; i < gb.length; ++i) + gb[i]=operator--(...borders[i]); + + guide g=ge[0]; + if(0 < gb.length) g=g&gb[0]; + for(int i=1; i < ge.length; ++i) { + g=g&ge[i]; + if(i < gb.length) g=g&gb[i]; + } + if(closed) g=g&cycle; + return g; +} - collect(points,c); // Required to join remaining case1 cycles. +// Connect contours into guides. +private guide[] connect(Segment[] S, pair[][] z, interpolate join) +{ + return sequence(new guide(int i) {return connect(S[i],z,join);},S.length); +} - return connect(points,c,join); +// Connect contours into guides. +private guide[][] connect(Segment[][] S, pair[][] z, interpolate join) +{ + guide[][] result=new guide[S.length][]; + for(int i=0; i < S.length; ++i) { + result[i]=connect(S[i],z,join); + } + return result; +} + +// Return contour guides for a 2D data array. +// z: two-dimensional array of nonoverlapping mesh points +// f: two-dimensional array of corresponding f(z) data values +// c: array of contour values +// join: interpolation operator (e.g. operator--or operator ..) +// subsample: number of interior points to include in each grid square +// (in addition to points on edge) +guide[][] contour(pair[][] z, real[][] f, real[] c, + interpolate join=operator--, int subsample=1) +{ + segment[][] seg=contouredges(f,c,subsample); + Segment[][] Seg=Segment(seg); + return connect(Seg,z,join); } // Return contour guides for a 2D data array on a uniform lattice // f: two-dimensional array of real data values -// midpoint: optional array containing data values at cell midpoints // a,b: diagonally opposite vertices of rectangular domain // c: array of contour values -// join: interpolation operator (e.g. operator -- or operator ..) -guide[][] contour(real[][] f, real[][] midpoint=new real[][], - pair a, pair b, real[] c, - interpolate join=operator --) +// join: interpolation operator (e.g. operator--or operator ..) +// subsample: number of interior points to include in each grid square +// (in addition to points on edge) +guide[][] contour(real[][] f, pair a, pair b, real[] c, + interpolate join=operator--, int subsample=1) { int nx=f.length-1; if(nx == 0) @@ -473,37 +1109,43 @@ guide[][] contour(real[][] f, real[][] midpoint=new real[][], zi[j]=(xi,interp(a.y,b.y,j/ny)); } } - return contour(z,f,midpoint,c,join); + return contour(z,f,c,join,subsample); } // return contour guides for a real-valued function -// f: real-valued function of two real variables -// a,b: diagonally opposite vertices of rectangular domain -// c: array of contour values -// nx,ny: number of subdivisions in x and y directions (determines accuracy) -// join: interpolation operator (e.g. operator -- or operator ..) +// f: real-valued function of two real variables +// a,b: diagonally opposite vertices of rectangular domain +// c: array of contour values +// nx,ny: number of subdivisions in x and y directions(determines accuracy) +// join: interpolation operator (e.g. operator--or operator ..) +// subsample: number of interior points to include in each grid square +// (in addition to points on edge) guide[][] contour(real f(real, real), pair a, pair b, real[] c, int nx=ngraph, int ny=nx, - interpolate join=operator --) + interpolate join=operator--, int subsample=1) { - // evaluate function at points and midpoints + // evaluate function at points and subsample real[][] dat=new real[nx+1][ny+1]; - real[][] midpoint=new real[nx+1][ny+1]; for(int i=0; i <= nx; ++i) { real x=interp(a.x,b.x,i/nx); - real x2=interp(a.x,b.x,(i+0.5)/nx); real[] dati=dat[i]; - real[] midpointi=midpoint[i]; for(int j=0; j <= ny; ++j) { dati[j]=f(x,interp(a.y,b.y,j/ny)); - midpointi[j]=f(x2,interp(a.y,b.y,(j+0.5)/ny)); } } - return contour(dat,midpoint,a,b,c,join); + return contour(dat,a,b,c,join,subsample); } - + +guide[][] contour(real f(pair), pair a, pair b, + real[] c, int nx=ngraph, int ny=nx, + interpolate join=operator--, int subsample=1) +{ + return contour(new real(real x, real y) {return f((x,y));}, + a,b,c,nx,ny,join,subsample); +} + void draw(picture pic=currentpicture, Label[] L=new Label[], guide[][] g, pen[] p) { @@ -530,6 +1172,239 @@ void draw(picture pic=currentpicture, Label[] L=new Label[], draw(pic,L,g,sequence(new pen(int) {return p;},g.length)); } +// Draw the contour +void draw(picture pic=currentpicture, Label L, + guide[] g, pen p=currentpen) +{ + draw(pic,g,p); + for(int i=0; i < g.length; ++i) { + if(L.s != "" && size(g[i]) > 1) { + label(pic,L,g[i],p); + } + } +} + +/* CONTOURS FOR IRREGULARLY SPACED POINTS */ +// +// +---------+ +// |\ /| +// | \ / | +// | \ / | +// | \ / | +// | X | +// | / \ | +// | / \ | +// | / \ | +// |/ \| +// +---------+ +// + +// Is triangle p0--p1--p2--cycle counterclockwise ? +private bool isCCW(pair p0, pair p1, pair p2) {return side(p0,p1,p2) < 0;} + +private struct segment +{ + bool active; + bool reversed; // True if lower values are to the left along line a--b. + pair a,b; // Endpoints; a is always an edge point if one exists. + int c; // Contour value. +} + +// Case 1: line passes through two vertices of a triangle +private segment case1(pair p0, pair p1, pair p2, + real v0, real v1, real v2) +{ + // Will cause a duplicate guide; luckily case1 is rare + segment rtrn; + rtrn.active=true; + rtrn.a=p0; + rtrn.b=p1; + rtrn.reversed=(isCCW(p0,p1,p2) ^(v2 > 0)); + return rtrn; +} + +// Cases 2 and 3: line passes through a vertex and a side of a triangle +//(the first vertex passed and the side between the other two) +private segment case2(pair p0, pair p1, pair p2, + real v0, real v1, real v2) +{ + segment rtrn; + rtrn.active=true; + pair val=interp(p1,p2,abs(v1/(v2-v1))); + rtrn.a=val; + rtrn.b=p0; + rtrn.reversed=!(isCCW(p0,p1,p2) ^(v2 > 0)); + return rtrn; +} + +private segment case3(pair p0, pair p1, pair p2, + real v0, real v1, real v2) +{ + segment rtrn; + rtrn.active=true; + pair val=interp(p1,p2,abs(v1/(v2-v1))); + rtrn.a=p0; + rtrn.b=val; + rtrn.reversed=(isCCW(p0,p1,p2) ^(v2 > 0)); + return rtrn; +} + +// Case 4: line passes through two sides of a triangle +//(through the sides formed by the first&second, and second&third vertices) +private segment case4(pair p0, pair p1, pair p2, + real v0, real v1, real v2) +{ + segment rtrn; + rtrn.active=true; + rtrn.a=interp(p1,p0,abs(v1/(v0-v1))); + rtrn.b=interp(p1,p2,abs(v1/(v2-v1))); + rtrn.reversed=(isCCW(p0,p1,p2) ^(v2 > 0)); + return rtrn; +} + +// Check if a line passes through a triangle, and draw the required line. +private segment checktriangle(pair p0, pair p1, pair p2, + real v0, real v1, real v2) +{ + // default null return + static segment dflt; + + real eps=eps*max(abs(v0),abs(v1),abs(v2),1); + + if(v0 < -eps) { + if(v1 < -eps) { + if(v2 < -eps) return dflt; // nothing to do + else if(v2 <= eps) return dflt; // nothing to do + else return case4(p0,p2,p1,v0,v2,v1); + } else if(v1 <= eps) { + if(v2 < -eps) return dflt; // nothing to do + else if(v2 <= eps) return case1(p1,p2,p0,v1,v2,v0); + else return case3(p1,p0,p2,v1,v0,v2); + } else { + if(v2 < -eps) return case4(p0,p1,p2,v0,v1,v2); + else if(v2 <= eps) + return case2(p2,p0,p1,v2,v0,v1); + else return case4(p1,p0,p2,v1,v0,v2); + } + } else if(v0 <= eps) { + if(v1 < -eps) { + if(v2 < -eps) return dflt; // nothing to do + else if(v2 <= eps) return case1(p0,p2,p1,v0,v2,v1); + else return case2(p0,p1,p2,v0,v1,v2); + } else if(v1 <= eps) { + if(v2 < -eps) return case1(p0,p1,p2,v0,v1,v2); + else if(v2 <= eps) return dflt; // use finer partitioning. + else return case1(p0,p1,p2,v0,v1,v2); + } else { + if(v2 < -eps) return case2(p0,p1,p2,v0,v1,v2); + else if(v2 <= eps) return case1(p0,p2,p1,v0,v2,v1); + else return dflt; // nothing to do + } + } else { + if(v1 < -eps) { + if(v2 < -eps) return case4(p1,p0,p2,v1,v0,v2); + else if(v2 <= eps) + return case2(p2,p0,p1,v2,v0,v1); + else return case4(p0,p1,p2,v0,v1,v2); + } else if(v1 <= eps) { + if(v2 < -eps) return case3(p1,p0,p2,v1,v0,v2); + else if(v2 <= eps) return case1(p1,p2,p0,v1,v2,v0); + else return dflt; // nothing to do + } else { + if(v2 < -eps) return case4(p0,p2,p1,v0,v2,v1); + else if(v2 <= eps) return dflt; // nothing to do + else return dflt; // nothing to do + } + } +} + +// Collect connecting path segments. +private void collect(pair[][][] points, real[] c) +{ + for(int cnt=0; cnt < c.length; ++cnt) { + pair[][] gdscnt=points[cnt]; + for(int i=0; i < gdscnt.length; ++i) { + pair[] gig=gdscnt[i]; + int Li=gig.length; + for(int j=i+1; j < gdscnt.length; ++j) { + pair[] gjg=gdscnt[j]; + int Lj=gjg.length; + if(abs(gig[0]-gjg[Lj-1]) < eps) { + gig.delete(0); + gdscnt[j].append(gig); + gdscnt.delete(i); + --i; + break; + } else if(abs(gig[Li-1]-gjg[0]) < eps) { + gjg.delete(0); + gig.append(gjg); + gdscnt[j]=gig; + gdscnt.delete(i); + --i; + break; + } + } + } + } +} + +// Join path segments. +private guide[][] connect(pair[][][] points, real[] c, interpolate join) +{ + // set up return value + guide[][] result=new guide[c.length][]; + for(int cnt=0; cnt < c.length; ++cnt) { + pair[][] pointscnt=points[cnt]; + guide[] resultcnt=result[cnt]=new guide[pointscnt.length]; + for(int i=0; i < pointscnt.length; ++i) { + pair[] pts=pointscnt[i]; + guide gd; + if(pts.length > 0) { + if(pts.length > 1 && abs(pts[0]-pts[pts.length-1]) < eps) { + guide[] g=sequence(new guide(int i) { + return pts[i]; + },pts.length-1); + g.push(cycle); + gd=join(...g); + } else + gd=join(...sequence(new guide(int i) { + return pts[i]; + },pts.length)); + } + resultcnt[i]=gd; + } + } + return result; +} + +guide[][] contour(pair[] z, real[] f, real[] c, interpolate join=operator--) +{ + if(z.length != f.length) + abort("z and f arrays have different lengths"); + + int[][] trn=triangulate(z); + + // array to store guides found so far + pair[][][] points=new pair[c.length][][]; + + for(int cnt=0; cnt < c.length; ++cnt) { + pair[][] pointscnt=points[cnt]; + real C=c[cnt]; + for(int i=0; i < trn.length; ++i) { + int[] trni=trn[i]; + int i0=trni[0], i1=trni[1], i2=trni[2]; + segment seg=checktriangle(z[i0],z[i1],z[i2],f[i0]-C,f[i1]-C,f[i2]-C); + if(seg.active) + pointscnt.push(seg.reversed ? new pair[] {seg.b,seg.a} : + new pair[] {seg.a,seg.b}); + } + } + + collect(points,c); + + return connect(points,c,join); +} + // Extend palette by the colors below and above at each end. pen[] extend(pen[] palette, pen below, pen above) { pen[] p=copy(palette); @@ -614,69 +1489,3 @@ void fill(picture pic=currentpicture, guide[][] g, pen[][] palette) } } } - -// routines for irregularly spaced points: - -// check existing guides and adds new segment to them if possible, -// or otherwise store segment as a new guide -private void addseg(pair[][] gds, segment seg) -{ - if(!seg.active) return; - // search for a path to extend - for(int i=0; i < gds.length; ++i) { - pair[] gd=gds[i]; - if(abs(gd[0]-seg.b) < eps) { - gd.insert(0,seg.a); - return; - } else if(abs(gd[gd.length-1]-seg.b) < eps) { - gd.push(seg.a); - return; - } else if(abs(gd[0]-seg.a) < eps) { - gd.insert(0,seg.b); - return; - } else if(abs(gd[gd.length-1]-seg.a) < eps) { - gd.push(seg.b); - return; - } - } - - // in case nothing is found - pair[] segm; - segm=new pair[] {seg.a,seg.b}; - gds.push(segm); - - return; -} - -guide[][] contour(real f(pair), pair a, pair b, - real[] c, int nx=ngraph, int ny=nx, - interpolate join=operator --) -{ - return contour(new real(real x, real y) {return f((x,y));},a,b,c,nx,ny,join); -} - -guide[][] contour(pair[] z, real[] f, real[] c, interpolate join=operator --) -{ - if(z.length != f.length) - abort("z and f arrays have different lengths"); - - int[][] trn=triangulate(z); - - // array to store guides found so far - pair[][][] points=new pair[c.length][][]; - - for(int cnt=0; cnt < c.length; ++cnt) { - pair[][] pointscnt=points[cnt]; - real C=c[cnt]; - for(int i=0; i < trn.length; ++i) { - int[] trni=trn[i]; - int i0=trni[0], i1=trni[1], i2=trni[2]; - addseg(pointscnt,checktriangle(z[i0],z[i1],z[i2], - f[i0]-C,f[i1]-C,f[i2]-C)); - } - } - - collect(points,c); - - return connect(points,c,join); -} |