func plc_labels (z, y, x, ireg, levs=, region=, condraw=, type=, ngradmax=, maxlabs=, labcol=, opaque=, width=, height=, spv=,concol=) /* DOCUMENT plc_labels, z, y, x, ireg, levs=levels, type=, ngradmax=, maxlabs=, labcol=, opaque=, width=, height=; draw labels (and optionally draw the contours themselves) of the point centered values Z on the 2D mesh Y,X,IREG. by using the contour vertices. Labels will not be drawn if gradients are too strong. contour vertices are currently returned by svn's fillc routine zone indices are calculated with digit2 this functionality might later be done more efficiently by direct calls to yorick routines (ie mesh_loc) named arguments levs= an array of contour levels to label (not optional) region= the region number to find labels in (need ireg array) condraw= if condraw == 1, then draw the contours themselves type= type of line if condraw defined ngradmax= a measure of the non-dimensional max gradient below which labels can be drawn (larger implies more labels), default=20 maxlabs= approximate max number labels to put on contour line segment (default 1) labcol= color to draw the label with (default black) concol= color to draw the contour with (default black) opaque= if defined, put white down behind the label height= size of characters KEYWORDS: region, SEE ALSO: fillc, plc */ { local dz, nx, ny, xmax, xmin, ymax, ymin, xc, yc, zc; require, "digit2.i"; dz = dimsof(z); nx = dz(2); ny = dz(3); if (is_void(maxlabs)) maxlabs = 1; if (is_void(type)) type = "solid"; if (is_void(width)) width = 1; if (is_void(ngradmax)) ngradmax = 20; if (is_void(labcol)) labcol="black"; if (is_void(concol)) concol="black"; if (is_void(levs)) { // error,"the levs keyword has not been set"; levs = find_nice_contours(z,nlevs=10); // print, "levs set to ", levs; } dz = dimsof(z); if (is_void(x)) { x = indgen(dz(2))(,-:1:dz(3)); print, " plc_labels: x array not set, using indices ", dz(2); } if (is_void(y)) { y = indgen(dz(3))(-:1:dz(2),); print, " plc_labels: y array not set, using indices ", dz(2); } if (sum(dimsof(z) == dimsof(y) == dimsof(x))) { error,"plc_labels: dimsof x, y, and z arrays are inconsistent"; } //maxseg= 1000; //maxpoly= 1000; //maxnode= 16317; xmax = max(x); xmin = min(x); ymax = max(y); ymin = min(y); small = (xmax-xmin)*1.e-7; if (!is_void(ireg) && !is_void(region)) { z2 = z; // make a copy we can mess with miss_value = 2*max(z2); // assign a missing value so no contours out of the region of interest; list = where (ireg != region); if (numberof(list) > 0) { z2(list) = miss_value; } fillout = fillc(z2, levs, x , y , maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1, miss_value=miss_value, nested=1); list = where(*fillout.zc != numberof(levs)); // indices where there are no special values // move the returned values into arrays for the region of interest.; ib = 1; nc = []; zc = []; yc = []; xc = []; for (i=1;i<=numberof(list);i++) { ie = ib + (*fillout.nc)(i) - 1; grow, nc, (*fillout.nc)(i); grow, zc, (*fillout.zc)(i); grow, yc, (*fillout.yc)(ib:ie); grow, xc, (*fillout.xc)(ib:ie); ib = ie+1; } numnodes_contour = nc; // read,prompt="debug",dummy; } else { // print, "using missing value: ", spv; fillout = fillc(z, levs, x, y, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode, contours=1,nested=1, miss_value=spv); zc = *fillout.zc; yc = *fillout.yc; xc = *fillout.xc; numnodes_contour = *fillout.nc; } if (numberof(numnodes_contour) == 0) { //print,"early return"; return; } // find the zones for each node in the contours nd = digit2(yc, xc, y, x); list = where(nd <= 0); if (numberof(list) > 0) nd(list) = 1; // anything outside the zones is just point in first zone id = nd%(nx-1); jd = (nd - id)/(nx-1) + 1; n = numberof(numnodes_contour); // a measure of the gradient, normalized by the dimensions of the plot; grada = (xmax-xmin)*(abs(z(dif,)/(x(dif,)+small)))(pcen,) + (ymax-ymin)*(abs(z(,dif)/(y(,dif)+small)))(,pcen) ; // print, "mingrada, max ", min(grada), max(grada); // read,prompt="debug",dummy; mgrad = max(grada); grada(1:3,) = grada(-2:0,) = 1000.*mgrad; grada(,1:3) = grada(,-2:0) = 1000.*mgrad; /* // a measure of the curvature (2nd derivative) normalized by the dimensions of the plot; curva = array(0.,dimsof(z)); xx = x(zcen,) yy = y(,zcen) zx = ((z(dif,)/x(dif,))(dif,))/xx(dif,); zy = ((z(,dif)/y(,dif))(,dif))/yy(,dif); curva(2:-1,) = abs(zx) + curva(2:-1,); curva(,2:-1) = abs(zy) + curva(,2:-1); curva(1:3,) = curva(-2:0,) = 1000.*max(curva) curva(,1:3) = curva(,-2:0) = 1000.*max(curva) // read,prompt="debug curva",dummy; */ labx = []; laby = []; lablev = []; nlabs = 0; kb = 1; for (k=1; k<=n; k++) { ke = kb + numnodes_contour(k) - 1; izone = id(kb:ke); jzone = jd(kb:ke); ndm = izone + nx*(jzone-1); xa = xc(kb:ke); ya = yc(kb:ke); cz = zc(k)+1; // print, " plotting contour, level number, val, nnodes ", k, cz, levs(cz), numnodes_contour(k); //if (!is_void(condraw)) plg, yc(kb:ke), xc(kb:ke), marks=0, if (condraw == 1) plg, yc(kb:ke), xc(kb:ke), marks=0, type=type,width=width,color=concol; // a measure of the wiggliness of the contour; wig = array(0.,dimsof(xa)); if (numnodes_contour(k) > 2) { dr = sqrt(xa(dif)^2 + ya(dif)^2)(cum); drd = dr(dif)+small; drz = dr(zcen); drzd = drz(dif)+small; // xx = xa(zcen); // yy = ya(zcen); y2 = ((ya(dif)/drd)(dif))/drzd; x2 = ((xa(dif)/drd)(dif))/drzd; wig(2:-1) = sqrt(y2*y2 + x2*x2) + wig(2:-1,); // print, "y2", y2; // print, "x2", x2; // print, "wig bef smooth", wig; wig = wig(zcen)(pcen); // print, "wig aft smooth", wig; wig = wig*10/(avg(wig)+1.e-20); wig(1:3,) = wig(-2:0,) = 1000.*max(wig); // print, "renormalized wig", wig; // read, prompt="wig debug",dummy; } // print, " at wig1, number nodes ", numberof(ndm); // a measure of the number of contours in the vicinity of this point; czp1 = min(numberof(levs),cz+1); czm1 = max(1, cz-1); // print, " czp1, czm1, levs, normal ", czp1, czm1, levs(czp1), levs(czm1), // (czp1-czm1)/(levs(czp1)-levs(czm1)); // print, "min, max of grada ", min(grada(ndm)), max(grada(ndm)); if (czp1 != czm1) { // ngrad will be a measure of the number of contours which would ; // span the x/y domain, if the gradient were constant at ; // a value of grada(ndm1) everywhere; // so a number greater than about 10 implies there are too many contours around to be good; // ie the gradient is too steep ngrad = grada(ndm)*(czp1-czm1)/(levs(czp1)-levs(czm1)); list = where( ngrad < ngradmax); } else { list = array(1,numberof(ndm)); } if (numberof(list) == 0) { // print, "ndm", ndm; // print, "grada", grada(ndm); // print, "czp1, czm1, levs ", czp1, czm1, levs(czp1), levs(czm1); // print, " the gradient is too steep for this contour, skipping ", ngrad; // read, prompt="steep debug", dummy; kb = ke + 1; continue; } ndm = ndm(list); xa = xa(list); ya = ya(list); wig = wig(list); // print, " at wig, number nodes ", numberof(ndm); // dont even consider any contours which are very small; extent = (max(xa)-min(xa))/(xmax-xmin) + (max(ya)-min(ya))/(ymax-ymin) ; // print, " extent is ", extent; numplt = 0; if (extent > 0.1) { while (numberof(ndm) > 0) { if (numplt >= maxlabs) { // print, " numplt >= maxlabs "; break; } // now check the potential label sites against previous labels; // and eliminate any that are too close; if (nlabs > 0) { dist = abs(xa-labx(-,))/(xmax-xmin) + abs(ya-laby(-,))/(ymax-ymin); mdist = dist(,min); cmax = array(.10,numberof(xa)); // read, prompt="debug", dummy; // list = where((dist - cmax)(,min) > 0.); list = where((mdist > cmax)); // print, " reducing numberof nodes closeness from ", numberof(xa), "to ", numberof(list); ndm = ndm(list); xa = xa(list); ya = ya(list); wig = wig(list); mdist = mdist(list) } if (numberof(ndm) > 0) { // now pick out the minimum gradient for the node to actually plot; // print, " reduced to ", numberof(ndm), " points"; if (czp1 != czm1) { ngrad = grada(ndm)*(czp1-czm1)/(levs(czp1)-levs(czm1)); } else { ngrad = array(1.,numberof(ndm)); } // print, " ngrad ", ngrad; // print, " wig ", wig; // a normalized version of the curvature with a mean value of order 1; // ncurva = curva(ndm)*10/avg(curva(ndm)); // print, " ncurva ", ncurva; if (nlabs == 0) { //works penalty = ngrad+ncurva; penalty = ngrad+wig; ndist = array(0.,dimsof(wig)); } else { ndist = 1./mdist; // print, " penalty for being near another label ", ndist; //works penalty = ngrad+ncurva+ndist; penalty = ngrad+wig+ndist; } // print, " total penalty ", penalty; list =penalty(mnx); // print, " min at ", list, ngrad(list), wig(list), ndist(list), penalty(list); xl = xa(list); yl = ya(list); // wig = wig(list); // print, "plotting a label", levs(cz); plt, pr1(levs(cz)), xl, yl, tosys=1, justify="CH",color=labcol, opaque=opaque, height=height; // read,prompt=" just plotted", dummy grow, labx, xl; grow, laby, yl; nlabs = nlabs + 1; numplt = numplt + 1; } else { // print, "ndm is zero"; } } // read, prompt="debug labels...", dummy; kb = ke + 1; } else { // print, " contour too short, skipping"; kb = ke + 1; } } }