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.asy1705
1 files changed, 448 insertions, 1257 deletions
diff --git a/Master/texmf-dist/asymptote/contour.asy b/Master/texmf-dist/asymptote/contour.asy
index 27a79439e2d..fbb4cd1c71a 100644
--- a/Master/texmf-dist/asymptote/contour.asy
+++ b/Master/texmf-dist/asymptote/contour.asy
@@ -1,922 +1,307 @@
-/*
- 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.
-*/
-
+// Contour routines written by Radoslav Marinov and John Bowman.
+
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;
+real eps=10000*realEpsilon;
+
+// 1
+// 6 +-------------------+ 5
+// | \ / |
+// | \ / |
+// | \ / |
+// | \ / |
+// 2 | X | 0
+// | / \ |
+// | / \ |
+// | / \ |
+// | / \ |
+// 7 +-------------------+ 4 or 8
+// 3
-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 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.
}
-private bool same(gridpoint gp1, gridpoint gp2)
+// Case 1: line passes through two vertices of a triangle
+private segment case1(pair p0, pair p1, int edge)
{
- return abs(gp2.z-gp1.z+(gp2.i-gp1.i,gp2.j-gp1.j)) < eps;
+ // 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;
}
-
-// 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;
+// 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;
}
+ 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;
- }
- }
+// 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;
}
+// 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;
-/*
- 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.
+ real eps=eps*max(abs(v0),abs(v1),abs(v2));
- 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
+ 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 {
- 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;
+ 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);
} 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;
+ 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
+ } 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
+ }
+ }
}
-
-// 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;
+// 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);}
- 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);
+ 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;
+ }
+ }
}
- 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();
+// 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));
}
- 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();
+ resultcnt[i]=gd;
}
}
- 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);
+ return result;
}
-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.
+// 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
+// midpoint: optional array containing values of f at cell midpoints
// 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)
+// join: interpolation operator (e.g. operator -- or operator ..)
+guide[][] contour(pair[][] z, real[][] f,
+ real[][] midpoint=new real[][], real[] c,
+ interpolate join=operator --)
{
- 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");
+ 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");
c=sort(c);
- boxdata[][] bd=new boxdata[nx][ny];
-
- segment[][] result=new segment[c.length][];
+ bool midpoints=midpoint.length > 0;
+ segment segments[][][]=new segment[nx][ny][];
+
+ // go over region a rectangle at a time
for(int i=0; i < nx; ++i) {
- boxdata[] bdi=bd[i];
+ pair[] zi=z[i];
+ pair[] zp=z[i+1];
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) {
- boxdata bdij=bdi[j]=boxdata(fi[j],fp[j],fi[j+1],fp[j+1]);
-
+ 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
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;
-
+ 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
+ // 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(f00);
- check(f10);
- check(f01);
- check(f11);
+ check(vertdat0);
+ check(vertdat1);
+ check(vertdat2);
+ check(vertdat3);
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);
+ // 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));
return 0;
}
-
+
void process(int l, int u) {
if(l >= u) return;
int i=quotient(l+u,2);
@@ -928,171 +313,150 @@ segment[][] contouredges(real[][] f, real[] c, int subsample=1)
process(i+1,u);
}
}
-
+
process(0,c.length);
}
}
-
- // Find contours and follow them
+
+ // set up return value
+ pair[][][] points=new pair[c.length][][];
+
for(int i=0; i < nx; ++i) {
- boxdata[] bdi=bd[i];
+ segment[][] segmentsi=segments[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;
- }
+ 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;
+ }
+ }
+ }
+ return -1;
}
-
- // 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;
+
+ 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;
+ }
+ }
+ }
}
+ search();
}
- });
+ }
+ }
}
- 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]);
+ // Follow contour in cell
+ int edge=forward(i,j,first=false);
- 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;
-}
+ // Follow contour forward outside of cell
+ follow(forward,edge);
-// 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);
-}
+ // Follow contour backward outside of cell
+ follow(backward,C.edge);
-// 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);
+ points[C.c].push(g);
+ }
+ }
}
- 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);
+ collect(points,c); // Required to join remaining case1 cycles.
+
+ return connect(points,c,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 ..)
-// 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)
+// 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 --)
{
int nx=f.length-1;
if(nx == 0)
@@ -1109,43 +473,37 @@ guide[][] contour(real[][] f, pair a, pair b, real[] c,
zi[j]=(xi,interp(a.y,b.y,j/ny));
}
}
- return contour(z,f,c,join,subsample);
+ return contour(z,f,midpoint,c,join);
}
// 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)
+// 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 ..)
guide[][] contour(real f(real, real), pair a, pair b,
real[] c, int nx=ngraph, int ny=nx,
- interpolate join=operator--, int subsample=1)
+ interpolate join=operator --)
{
- // evaluate function at points and subsample
+ // evaluate function at points and midpoints
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,a,b,c,join,subsample);
+ return contour(dat,midpoint,a,b,c,join);
}
-
-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)
{
@@ -1172,239 +530,6 @@ 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);
@@ -1489,3 +614,69 @@ 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);
+}