The R Project SVN – R – Diff – /trunk/src/main/plot3d.c

 	return 0;
     else
 	return -n;
+}
+/*
+ * Generate a list of segments for a single level
+ */
+static SEGP* contourLines(double *x, int nx, double *y, int ny,
+			 double *z, double zc, double atom)
+{
+    double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];
+    int i, j, k, l, m, nacode;
+    SEGP seglist;
+    SEGP *segmentDB;
+    /* Initialize the segment data base */
+    /* Note we must be careful about resetting */
+    /* the top of the stack, otherwise we run out of */
+    /* memory after a sequence of displaylist replays */
+    /*
+     * This reset is done out in GEcontourLines
+     */
+    segmentDB = (SEGP*)R_alloc(nx*ny, sizeof(SEGP));
+    for (i = 0; i < nx; i++)
+	for (j = 0; j < ny; j++)
+	    segmentDB[i + j * nx] = NULL;
+    for (i = 0; i < nx - 1; i++) {
+	xl = x[i];
+	xh = x[i + 1];
+	for (j = 0; j < ny - 1; j++) {
+	    yl = y[j];
+	    yh = y[j + 1];
+	    k = i + j * nx;
+	    zll = z[k];
+	    zhl = z[k + 1];
+	    zlh = z[k + nx];
+	    zhh = z[k + nx + 1];
+	    /* If the value at a corner is exactly equal to a contour level,
+	     * change that value by a tiny amount */
+	    if (zll == zc) zll += atom;
+	    if (zhl == zc) zhl += atom;
+	    if (zlh == zc) zlh += atom;
+	    if (zhh == zc) zhh += atom;
+#ifdef DEBUG_contour
+	    /* Haven't seen this happening (MM): */
+	    if (zll == zc) REprintf(" [%d,%d] ll: %g\n",i,j, zll);
+	    if (zhl == zc) REprintf(" [%d,%d] hl: %g\n",i,j, zhl);
+	    if (zlh == zc) REprintf(" [%d,%d] lh: %g\n",i,j, zlh);
+	    if (zhh == zc) REprintf(" [%d,%d] hh: %g\n",i,j, zhh);
+#endif
+	    /* Check for intersections with sides */
+	    nacode = 0;
+	    if (R_FINITE(zll)) nacode += 1;
+	    if (R_FINITE(zhl)) nacode += 2;
+	    if (R_FINITE(zlh)) nacode += 4;
+	    if (R_FINITE(zhh)) nacode += 8;
+	    k = 0;
+	    switch (nacode) {
+	    case 15:
+		if (ctr_intersect(zll, zhl, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yl; k++;
+		}
+		if (ctr_intersect(zll, zlh, zc, &f)) {
+		    yy[k] = yl + f * (yh - yl);
+		    xx[k] = xl; k++;
+		}
+		if (ctr_intersect(zhl, zhh, zc, &f)) {
+		    yy[k] = yl + f * (yh - yl);
+		    xx[k] = xh; k++;
+		}
+		if (ctr_intersect(zlh, zhh, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yh; k++;
+		}
+		break;
+	    case 14:
+		if (ctr_intersect(zhl, zhh, zc, &f)) {
+		    yy[k] = yl + f * (yh - yl);
+		    xx[k] = xh; k++;
+		}
+		if (ctr_intersect(zlh, zhh, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yh; k++;
+		}
+		if (ctr_intersect(zlh, zhl, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yh + f * (yl - yh);
+		    k++;
+		}
+		break;
+	    case 13:
+		if (ctr_intersect(zll, zlh, zc, &f)) {
+		    yy[k] = yl + f * (yh - yl);
+		    xx[k] = xl; k++;
+		}
+		if (ctr_intersect(zlh, zhh, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yh; k++;
+		}
+		if (ctr_intersect(zll, zhh, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yl + f * (yh - yl);
+		    k++;
+		}
+		break;
+	    case 11:
+		if (ctr_intersect(zhl, zhh, zc, &f)) {
+		    yy[k] = yl + f * (yh - yl);
+		    xx[k] = xh; k++;
+		}
+		if (ctr_intersect(zll, zhl, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yl; k++;
+		}
+		if (ctr_intersect(zll, zhh, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yl + f * (yh - yl);
+		    k++;
+		}
+		break;
+	    case 7:
+		if (ctr_intersect(zll, zlh, zc, &f)) {
+		    yy[k] = yl + f * (yh - yl);
+		    xx[k] = xl; k++;
+		}
+		if (ctr_intersect(zll, zhl, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yl; k++;
+		}
+		if (ctr_intersect(zlh, zhl, zc, &f)) {
+		    xx[k] = xl + f * (xh - xl);
+		    yy[k] = yh + f * (yl - yh);
+		    k++;
+		}
+		break;
+	    }
+	    /* We now have k(=2,4) endpoints */
+	    /* Decide which to join */
+	    seglist = NULL;
+	    if (k > 0) {
+		if (k == 2) {
+		    seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
+		}
+		else if (k == 4) {
+		    for (k = 3; k >= 1; k--) {
+			m = k;
+			xl = xx[k];
+			for (l = 0; l < k; l++) {
+			    if (xx[l] > xl) {
+				xl = xx[l];
+				m = l;
+			    }
+			}
+			if (m != k) {
+			    xl = xx[k];
+			    yl = yy[k];
+			    xx[k] = xx[m];
+			    yy[k] = yy[m];
+			    xx[m] = xl;
+			    yy[m] = yl;
+			}
+		    }
+		    seglist = ctr_newseg(xx[0], yy[0], xx[1], yy[1], seglist);
+		    seglist = ctr_newseg(xx[2], yy[2], xx[3], yy[3], seglist);
+		}
+		else error("k != 2 or 4");
+	    }
+	    segmentDB[i + j * nx] = seglist;
+	}
+    }
+    return segmentDB;
+}
+/* maximal number of line segments of one contour segment:
+ * for preventing infinite loops -- shouldn't be needed --> warning */
+#define MAX_ns 25000
+#define CONTOUR_LIST_STEP 100
+#define CONTOUR_LIST_LEVEL 0
+#define CONTOUR_LIST_X 1
+#define CONTOUR_LIST_Y 2
+static SEXP growList(SEXP oldlist) {
+    int i, len;
+    SEXP templist;
+    len = LENGTH(oldlist);
+    templist = PROTECT(allocVector(VECSXP, len + CONTOUR_LIST_STEP));
+    for (i=0; i<len; i++)
+	SET_VECTOR_ELT(templist, i, VECTOR_ELT(oldlist, i));
+    UNPROTECT(1);
+    return templist;
+}
+/*
+ * Store the list of segments for a single level in the SEXP
+ * list that will be returned to the user
+ */
+int addContourLines(double *x, int nx, double *y, int ny,
+		     double *z, double zc, double atom,
+		     SEGP* segmentDB, int nlines, SEXP container)
+{
+    double xend, yend;
+    int i, ii, j, jj, ns, ns2, dir, nc;
+    SEGP seglist, seg, s, start, end;
+    SEXP ctr, level, xsxp, ysxp;
+    /* Begin following contours. */
+    /* 1. Grab a segment */
+    /* 2. Follow its tail */
+    /* 3. Follow its head */
+    /* 4. Save the contour */
+    for (i = 0; i < nx - 1; i++)
+	for (j = 0; j < ny - 1; j++) {
+	    while ((seglist = segmentDB[i + j * nx])) {
+		ii = i; jj = j;
+		start = end = seglist;
+		segmentDB[i + j * nx] = seglist->next;
+		xend = seglist->x1;
+		yend = seglist->y1;
+		while ((dir = ctr_segdir(xend, yend, x, y,
+					 &ii, &jj, nx, ny))) {
+		    segmentDB[ii + jj * nx]
+			= ctr_segupdate(xend, yend, dir, TRUE,/* = tail */
+					segmentDB[ii + jj * nx], &seg);
+		    if (!seg) break;
+		    end->next = seg;
+		    end = seg;
+		    xend = end->x1;
+		    yend = end->y1;
+		}
+		end->next = NULL; /* <<< new for 1.2.3 */
+		ii = i; jj = j;
+		xend = seglist->x0;
+		yend = seglist->y0;
+		while ((dir = ctr_segdir(xend, yend, x, y,
+					 &ii, &jj, nx, ny))) {
+		    segmentDB[ii + jj * nx]
+			= ctr_segupdate(xend, yend, dir, FALSE,/* ie. head */
+					segmentDB[ii+jj*nx], &seg);
+		    if (!seg) break;
+		    seg->next = start;
+		    start = seg;
+		    xend = start->x0;
+		    yend = start->y0;
+		}
+		/* ns := #{segments of polyline} -- need to allocate */
+		s = start;
+		ns = 0;
+		/* MAX_ns: prevent inf.loop (shouldn't be needed) */
+		while (s && ns < MAX_ns) {
+		    ns++;
+		    s = s->next;
+		}
+		if(ns == MAX_ns)
+		    warning("contour(): circular/long seglist -- bug.report()!");
+		/* countour midpoint : use for labelling sometime (not yet!) */
+		if (ns > 3) ns2 = ns/2; else ns2 = -1;
+		/*
+		 * "write" the contour locations into the list of contours
+		 */
+		ctr = PROTECT(allocVector(VECSXP, 3));
+		level = PROTECT(allocVector(REALSXP, 1));
+		xsxp = PROTECT(allocVector(REALSXP, ns + 1));
+		ysxp = PROTECT(allocVector(REALSXP, ns + 1));
+		REAL(level)[0] = zc;
+		SET_VECTOR_ELT(ctr, CONTOUR_LIST_LEVEL, level);
+		s = start;
+		REAL(xsxp)[0] = s->x0;
+		REAL(ysxp)[0] = s->y0;
+		ns = 1;
+		while (s->next && ns < MAX_ns) {
+		    s = s->next;
+		    REAL(xsxp)[ns] = s->x0;
+		    REAL(ysxp)[ns++] = s->y0;
+		}
+		REAL(xsxp)[ns] = s->x1;
+		REAL(ysxp)[ns] = s->y1;
+		SET_VECTOR_ELT(ctr, CONTOUR_LIST_X, xsxp);
+		SET_VECTOR_ELT(ctr, CONTOUR_LIST_Y, ysxp);
+		/*
+		 * We're about to add another line to the list ...
+		 */
+		nlines += 1;
+		nc = LENGTH(VECTOR_ELT(container, 0));
+		if (nlines == nc)
+		    /* Where does this get UNPROTECTed? */
+		    SET_VECTOR_ELT(container, 0,
+				   growList(VECTOR_ELT(container, 0)));
+		SET_VECTOR_ELT(VECTOR_ELT(container, 0), nlines - 1, ctr);
+		UNPROTECT(4);
+	    }
+	}
+    return nlines;
+}
+/*
+ * Given nx x values, ny y values, nx*ny z values,
+ * and nl cut-values in z ...
+ * ... produce a list of contour lines:
+ *   list of sub-lists
+ *     sub-list = x vector, y vector, and cut-value.
+ */
+SEXP GEcontourLines(double *x, int nx, double *y, int ny,
+		    double *z, double *levels, int nl,
+		    GEDevDesc *dd)
+{
+    char *vmax;
+    int i, nlines, len;
+    double atom, zmin, zmax;
+    SEGP* segmentDB;
+    SEXP container, mainlist, templist;
+    /*
+     * "tie-breaker" values
+     */
+    zmin = DBL_MAX;
+    zmax = DBL_MIN;
+    for (i = 0; i < nx * ny; i++)
+	if (R_FINITE(z[i])) {
+	    if (zmax < z[i]) zmax =  z[i];
+	    if (zmin > z[i]) zmin =  z[i];
+	}
+    if (zmin >= zmax) {
+	if (zmin == zmax)
+	    warning("all z values are equal");
+	else
+	    warning("all z values are NA");
+	return R_NilValue;
+    }
+    /* change to 1e-3, reconsidered because of PR#897
+     * but 1e-7, and even  2*DBL_EPSILON do not prevent inf.loop in contour().
+     * maybe something like   16 * DBL_EPSILON * (..).
+     * see also MAX_ns above */
+    atom = 1e-3 * (zmax - zmin);
+    /*
+     * Create a "container" which is a list with only 1 element.
+     * The element is the list of lines that will be built up.
+     * I create the container because this allows me to PROTECT
+     * the container once here and then UNPROTECT it at the end of
+     * this function and, as long as I always work with
+     * VECTOR_ELT(container, 0) and SET_VECTOR_ELT(container, 0)
+     * in functions called from here, I don't need to worry about
+     * protectin the list that I am building up.
+     * Why bother?  Because the list I am building can potentially
+     * grow and it's awkward to get the PROTECTs/UNPROTECTs right
+     * when you're in a loop and growing a list.
+     */
+    container = PROTECT(allocVector(VECSXP, 1));
+    /*
+     * Create "large" list (will trim excess at the end if necesary)
+     */
+    SET_VECTOR_ELT(container, 0, allocVector(VECSXP, CONTOUR_LIST_STEP));
+    nlines = 0;
+    /*
+     * Add lines for each contour level
+     */
+    for (i = 0; i < nl; i++) {
+	/*
+	 * The vmaxget/set is to manage the memory that gets
+	 * R_alloc'ed in the creation of the segmentDB structure
+	 */
+	vmax = vmaxget();
+	/*
+	 * Generate a segment database
+	 */
+	segmentDB = contourLines(x, nx, y, ny, z, levels[i], atom);
+	/*
+	 * Add lines to the list based on the segment database
+	 */
+	nlines = addContourLines(x, nx, y, ny, z, levels[i],
+				 atom, segmentDB, nlines,
+				 container);
+	vmaxset(vmax);
+    }
+    /*
+     * Trim the list of lines to the appropriate length.
+     */
+    len = LENGTH(VECTOR_ELT(container, 0));
+    if (nlines < len) {
+	mainlist = VECTOR_ELT(container, 0);
+	templist = PROTECT(allocVector(VECSXP, nlines));
+	for (i=0; i<nlines; i++)
+	    SET_VECTOR_ELT(templist, i, VECTOR_ELT(mainlist, i));
+	mainlist = templist;
+	UNPROTECT(1);  /* UNPROTECT templist */
+    } else
+	mainlist = VECTOR_ELT(container, 0);
+    UNPROTECT(1);  /* UNPROTECT container */
+    return mainlist;
+}
+SEXP GEdrawContourLines();
 static void contour(SEXP x, int nx, SEXP y, int ny, SEXP z,
 		    double zc,
 		    SEXP labels, int cnum,
 		    Rboolean drawLabels, int method,
 		    Rboolean vectorFonts, int typeface, int fontindex,
 		    double atom, DevDesc *dd)
+{
 /* draw a contour for one given contour level `zc' */
-/* maximal number of line segments of one contour segment:
- * for preventing infinite loops -- shouldn't be needed --> warning */
-#define MAX_ns 25000
     char *vmax;
     double f, xl, xh, yl, yh, zll, zhl, zlh, zhh, xx[4], yy[4];
     double xend, yend;
     int i, ii, j, jj, k, l, m, nacode, ns, ns2, dir;
     Rboolean vectorFonts = FALSE;
     int typeface = 0;
     int fontindex = 0;
     double labcex;
     DevDesc *dd = CurrentDevice();
+    SEXP result = R_NilValue;
     GCheckState(dd);
     if (length(args) < 4)
 	errorcall(call, "too few arguments");
     ltysave = Rf_gpptr(dd)->lty;
     colsave = Rf_gpptr(dd)->col;
     lwdsave = Rf_gpptr(dd)->lwd;
     cexsave = Rf_gpptr(dd)->cex;
     labelList = PROTECT(R_NilValue);
+    /* draw contour for levels[i] */
     GMode(1, dd);
-    for (i = 0; i < nc; i++) {/* draw contour for levels[i] */
+    for (i = 0; i < nc; i++) {
 	vmax = vmaxget();
 	Rf_gpptr(dd)->lty = INTEGER(lty)[i % nlty];
 	if (Rf_gpptr(dd)->lty == NA_INTEGER)
 	    Rf_gpptr(dd)->lty = ltysave;
 	Rf_gpptr(dd)->col = INTEGER(col)[i % ncol];
 	contour(x, nx, y, ny, z, REAL(c)[i], labels, i,
 		drawLabels, method-1,
 		vectorFonts, typeface, fontindex, atom, dd);
 	vmaxset(vmax);
+    }
     GMode(0, dd);
     vmaxset(vmax0);
     Rf_gpptr(dd)->lty = ltysave;
     Rf_gpptr(dd)->col = colsave;
     Rf_gpptr(dd)->lwd = lwdsave;
     Rf_gpptr(dd)->cex = cexsave;
     UNPROTECT(5);
     /* NOTE: only record operation if no "error"  */
     /* NOTE: on replay, call == R_NilValue */
     if (GRecording(call))
 	recordGraphicOperation(op, oargs, dd);
-    return R_NilValue;
+    return result;
+}
 	/*  F i l l e d   C o n t o u r   P l o t s  */

The R Project SVN R

(root)/trunk/src/main/plot3d.c @ 28362 – Rev 20377 → 21061