summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/asymptote/contour.asy
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/asymptote/contour.asy')
-rw-r--r--Master/texmf-dist/asymptote/contour.asy1491
1 files changed, 1491 insertions, 0 deletions
diff --git a/Master/texmf-dist/asymptote/contour.asy b/Master/texmf-dist/asymptote/contour.asy
new file mode 100644
index 00000000000..27a79439e2d
--- /dev/null
+++ b/Master/texmf-dist/asymptote/contour.asy
@@ -0,0 +1,1491 @@
+/*
+ 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;
+
+private real eps=sqrtEpsilon;
+
+/*
+ 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;
+ }
+}
+
+private bool same(gridpoint gp1, gridpoint gp2)
+{
+ return abs(gp2.z-gp1.z+(gp2.i-gp1.i,gp2.j-gp1.j)) < eps;
+}
+
+
+// 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;
+ }
+}
+
+// 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;
+ }
+ }
+}
+
+
+/*
+ Determine contours within a unit square box.
+
+ 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 {
+ 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(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));
+
+ 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;
+
+ 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;
+}
+
+
+// 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();
+ }
+ 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);
+ }
+}
+
+private Segment[] Segment(segment[] s)
+{
+ return sequence(new Segment(int i) {return Segment(s[i]);},s.length);
+}
+
+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
+// c: array of contour values
+// subsample: number of points to use in each box in addition to endpoints
+segment[][] contouredges(real[][] f, real[] c, int subsample=1)
+{
+ 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);
+ boxdata[][] bd=new boxdata[nx][ny];
+
+ segment[][] result=new segment[c.length][];
+
+ for(int i=0; i < nx; ++i) {
+ boxdata[] bdi=bd[i];
+ real[] fi=f[i];
+ real[] fp=f[i+1];
+
+ for(int j=0; j < ny; ++j) {
+ boxdata bdij=bdi[j]=boxdata(fi[j],fp[j],fi[j+1],fp[j+1]);
+
+ int checkcell(int cnt) {
+ real C=c[cnt];
+
+ 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;
+
+
+ int countm=0;
+ int countz=0;
+ int countp=0;
+
+ void check(real vertdat) {
+ if(vertdat < -eps)++countm;
+ else {
+ if(vertdat <= eps)++countz;
+ else++countp;
+ }
+ }
+
+ 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;
+
+ // 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);
+ int sign=checkcell(i);
+ if(sign == -1) process(i+1,u);
+ else if(sign == 1) process(l,i);
+ else {
+ process(l,i);
+ process(i+1,u);
+ }
+ }
+
+ process(0,c.length);
+ }
+ }
+
+ // Find contours and follow them
+ for(int i=0; i < nx; ++i) {
+ boxdata[] bdi=bd[i];
+ for(int j=0; j < ny; ++j) {
+ 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;
+ }
+ }
+ });
+ }
+
+ // 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;
+ }
+ }
+ });
+ }
+
+ 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;
+}
+
+// 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);
+}
+
+// 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
+// a,b: diagonally opposite vertices of rectangular domain
+// 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(real[][] f, pair a, pair b, real[] c,
+ interpolate join=operator--, int subsample=1)
+{
+ 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");
+
+ pair[][] z=new pair[nx+1][ny+1];
+ for(int i=0; i <= nx; ++i) {
+ pair[] zi=z[i];
+ real xi=interp(a.x,b.x,i/nx);
+ for(int j=0; j <= ny; ++j) {
+ zi[j]=(xi,interp(a.y,b.y,j/ny));
+ }
+ }
+ 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 ..)
+// 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--, int subsample=1)
+{
+ // evaluate function at points and subsample
+ real[][] dat=new real[nx+1][ny+1];
+
+ for(int i=0; i <= nx; ++i) {
+ real x=interp(a.x,b.x,i/nx);
+ real[] dati=dat[i];
+ for(int j=0; j <= ny; ++j) {
+ dati[j]=f(x,interp(a.y,b.y,j/ny));
+ }
+ }
+
+ 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)
+{
+ begingroup(pic);
+ for(int cnt=0; cnt < g.length; ++cnt) {
+ guide[] gcnt=g[cnt];
+ pen pcnt=p[cnt];
+ for(int i=0; i < gcnt.length; ++i)
+ draw(pic,gcnt[i],pcnt);
+ if(L.length > 0) {
+ Label Lcnt=L[cnt];
+ for(int i=0; i < gcnt.length; ++i) {
+ if(Lcnt.s != "" && size(gcnt[i]) > 1)
+ label(pic,Lcnt,gcnt[i],pcnt);
+ }
+ }
+ }
+ endgroup(pic);
+}
+
+void draw(picture pic=currentpicture, Label[] L=new Label[],
+ guide[][] g, pen p=currentpen)
+{
+ 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);
+ p.insert(0,below);
+ p.push(above);
+ return p;
+}
+
+// Compute the interior palette for a sequence of cyclic contours
+// corresponding to palette.
+pen[][] interior(picture pic=currentpicture, guide[][] g, pen[] palette)
+{
+ if(palette.length != g.length+1)
+ abort("Palette array must have length one more than guide array");
+ pen[][] fillpalette=new pen[g.length][];
+ for(int i=0; i < g.length; ++i) {
+ guide[] gi=g[i];
+ guide[] gp;
+ if(i+1 < g.length) gp=g[i+1];
+ guide[] gm;
+ if(i > 0) gm=g[i-1];
+
+ pen[] fillpalettei=new pen[gi.length];
+ for(int j=0; j < gi.length; ++j) {
+ path P=gi[j];
+ if(cyclic(P)) {
+ int index=i+1;
+ bool nextinside;
+ for(int k=0; k < gp.length; ++k) {
+ path next=gp[k];
+ if(cyclic(next)) {
+ if(inside(P,point(next,0)))
+ nextinside=true;
+ else if(inside(next,point(P,0)))
+ index=i;
+ }
+ }
+ if(!nextinside) {
+ // Check to see if previous contour is inside
+ for(int k=0; k < gm.length; ++k) {
+ path prev=gm[k];
+ if(cyclic(prev)) {
+ if(inside(P,point(prev,0)))
+ index=i;
+ }
+ }
+ }
+ fillpalettei[j]=palette[index];
+ }
+ fillpalette[i]=fillpalettei;
+ }
+ }
+ return fillpalette;
+}
+
+// Fill the interior of cyclic contours with palette
+void fill(picture pic=currentpicture, guide[][] g, pen[][] palette)
+{
+ for(int i=0; i < g.length; ++i) {
+ guide[] gi=g[i];
+ guide[] gp;
+ if(i+1 < g.length) gp=g[i+1];
+ guide[] gm;
+ if(i > 0) gm=g[i-1];
+
+ for(int j=0; j < gi.length; ++j) {
+ path P=gi[j];
+ path[] S=P;
+ if(cyclic(P)) {
+ for(int k=0; k < gp.length; ++k) {
+ path next=gp[k];
+ if(cyclic(next) && inside(P,point(next,0)))
+ S=S^^next;
+ }
+ for(int k=0; k < gm.length; ++k) {
+ path next=gm[k];
+ if(cyclic(next) && inside(P,point(next,0)))
+ S=S^^next;
+ }
+ fill(pic,S,palette[i][j]+evenodd);
+ }
+ }
+ }
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-24 01:26:15 +0000